Keywords: estimation, vitesse angulaire, équations d'Euler, traitement du signal estimation, angular rate, Euler's equations, signal processing

Cette thèse étudie l'estimation de vitesse de rotation d'un corps rigide à partir de mesures de directions (par exemple champ magnétique, direction du soleil) embarquées. L'objectif est de remplacer les gyromètres qui sont chers comparés aux autres capteurs inertiels et sujets à des saturations et à des dysfonctionnements.

Dans une première partie de la thèse, on traite les cas spécifiques d'une rotation à axe fixe ou légèrement variable.

Dans une seconde partie, on traite le cas d'une rotation quelconque par un observateur asymptotique non-linéaire. On construit l'observateur à partir de mesures de deux vecteurs de référence non colinéaires, ou bien d'un seul vecteur. La connaissance des coordonnées inertielles des vecteurs de référence n'est pas nécéssaire. On étend ensuite l'observateur pour estimer en plus le couple et les paramètres d'inertie. Les équations d'Euler jouent un rôle central dans les travaux présentés ici.

Il apparaît que, du moins pour les illustrations considérées, les gyromètres peuvent être remplacés par un algorithme d'estimation basé sur des capteurs de direction qui sont bien moins chers et plus robustes.

Chapter 1 Introduction and presentation of the problem

This thesis addresses the general question of estimating the angular rate (also referred to as angular velocity) of a rigid body thanks to direction sensors attached to it. This broad question has applications in various fields of engineering and applied science. Let us start the discussion by a brief description of specific examples. i) In aerospace, the deployment phase of spinning satellites starts by a detumbling maneuver during which the angular rate is controlled in an active way until it reaches a zero value. The control strategy employs an estimation of this variable, in closed-loop. The satellite uses its own on-board sensors to evaluate this variable. ii) High velocity spinning objects are very common in ballistics. For example, the XM25 air-burst rifle (smart-weapon) fires smart shells which estimate their rotation to determine the traveled distance (so that explosion of the projectile can be activated at any user-defined distance). The number of spiral rotations after fire gives a direct reckoning of the range1 . The sensors are embedded inside the shell. iii) Finally, the problem of angular rate estimation can also be found in the emerging field of smart devices for sports. It is important for athletes in many sports to train their skills to spin a ball. This estimation can be achieved from on-board sensors, without requiring any external monitoring systems.

Brief mathematical description of the problem

To understand the scientific question at stake, we now sketch a mathematical formulation for it. The attitude of a rigid body is characterized by a rotation matrix R(t) defined as the mapping from an inertial frame of reference to a body frame, which is, by definition, attached to the rigid body. This rotation matrix satisfies the differential equation

Ṙ = R    0 -ω 3 ω 2 ω 3 0 -ω 1 -ω 2 ω 1 0    R [ω × ] (1.1)
where ω (ω 1 ω 2 ω 3 ) is the angular rate (or angular velocity) of the rigid body expressed in the body frame. The motion of ω itself is governed by the famed Euler's equations

[LL82] ω = J -1 (Jω × ω + τ ) (1.2)
where J is the matrix of inertia of the rigid body and τ is the external torque applied to the rigid body.

In the applications where it is desired to know R and ω with some degree of accuracy, the rigid body is usually equipped with inertial sensors attached to it 2 such as rate gyros and direction sensors.

Rate gyros measure the quantity ω up to noise and biases. Direction sensors work as follows. Consider a reference vector a expressed in the inertial frame. Then, the expression of a in the body frame at time t is

a(t) = R(t) a (1.3)
The variable a(t) is called a vector observation or vector measurement. This is the variable produced by a direction sensor such as magnetometers or Sun sensors, among several possibilities (see Table 1.1). The reference vector a itself may be known or unknown depending on the application. For example, the expression of the Earth magnetic field in an Earth-centered reference frame may be calculated with a high degree of accuracy from the International Geomagnetic Reference Field (IGRF) model which is built from actual experimental measurements, updated every five years.

vector sensor unit magnetic field magnetometer [G,T] light direction Sun sensors -(unit vector) gravity field accelerometer [m.s -2 ] Earth direction Earth sensor -(unit vector) Table 1.1: Examples of vector measurements and corresponding sensors It is well known that these sensors are sufficient to estimate R and ω. In fact, many attitude estimation algorithms have been proposed over the last years to estimate the attitude from biased gyro measurements together with vector observations. The sensor fusion can be done by Kalman Filter-like algorithms [CBIO06, SRK + 08] or by non-linear observers [MHP08, MS10, VSO08, TRB11, SMS11, GFJS12, TMHL12].

An important fact is that the gyros have numerous drawbacks. They are very expensive compared to direction sensors, prone to saturation during high rate rotations [START_REF] Bar-Itzhack | Classification of algorithms for angular velocity estimation[END_REF] and subject to failure. A famous example is the Hubble Space Telescope, which was put into "safe hold" (i.e. sleep mode) on November 13 rd , 1999 as four of its six gyros malfunctioned [START_REF] Oshman | Spacecraft angluar velocity estimation using sequential observations of a single directional vector[END_REF]. Eventually, the six gyros were replaced in 1999 and again in 2009 during the highly expensive missions STS-103 and STS-125 [START_REF] Lallo | Experience with the Hubble Space Telescope: 20 years of an archetype[END_REF].

Fortunately, they are not strictly necessary in the applications considered here. A recent trend is to replace (or consolidate) the gyro measurements with other sensors. In particular, vector measurements have been investigated for estimating ω for the last 15 years. An example of this "gyro-less" trend is the Solar Anomalous and Magnetospheric Particle Explorer (SAMPEX) spacecraft, which was launched in July 1992 and whose inertial unit only comprises Sun sensors and magnetometers [START_REF] Tsai | SAMPEX spin stabilized mode[END_REF].

This thesis follows this trend and proposes new results in this area. Formally, we address the following problem.

Problem 1.1: Main problem of the thesis

From one or several measurements of the type (1.3), find an estimate of the angular rate ω appearing in (1.1), assuming it satisfies (1.2).

State of the art

Several algorithms have been proposed over the last few years to address the Problem mentioned above. The survey [START_REF] Bar-Itzhack | Classification of algorithms for angular velocity estimation[END_REF] classifies them into two categories: derivative methods, which employ a time differentiation of the measurements, and estimation methods, which are essentially asymptotic estimators for an augmented state. In this thesis, we further categorize them as indirect methods, which require algebraic manipulation of the measurements to express the attitude, and direct methods, which do not. Due to the overlapping of categories, we finally identify four types of methods which are detailed below.

Indirect methods

Many angular rate estimation methods consider that the attitude matrix R is known3 . Implicitly, this assumes that the vector observations have already been processed to produce R. Presuming that at least two non-collinear reference vectors a, b are known, together with their expression (measurements) in the body frame a, b, the rotation matrix can be found as the solution of the classic Wahba problem [START_REF] Wahba | Problem 65-1: a least squares estimate of spacecraft attitude[END_REF]. This is a minimization problem on R. Many solutions to the Wahba problem have been proposed such as QUEST or TRIAD [START_REF] Shuster | Approximate algorithms for fast optimal attitude computation[END_REF][START_REF] Shuster | Kalman filtering of spacecraft attitude and the QUEST model[END_REF]. Regardless of the solution method, we now continue with the classification exposed in [START_REF] Bar-Itzhack | Classification of algorithms for angular velocity estimation[END_REF].

1. The indirect derivative approach consists in solving (1.1) to obtain4

[ω × ] = R Ṙ (1.4)
which identifies ω, unambiguously. In practice, R, whose value is obtained from the vector measurements, is corrupted by noise. The difficulty of this approach hence lies in the time-differentiation, which introduces high-frequency noise in the estimation of ω. As a result, the value ω obtained from (1.4) needs to be smoothed to produce the estimation ω. The smoothing can be achieved by a passive low-pass filter or by an active filter which combines ω with the (state) Euler's equations to produce an estimation ω that is optimal with respect to some noise model. A schematic view of the indirect derivative approach is represented in Figure 1.1.

Very generally, the active filter can be an Extended Kalman Filter (EKF) or one of its numerous derived forms [ABIH98, ABIDH99, HBI99, TM12].

2. The indirect estimation approach consists in computing an asymptotic observer for the augmented state S (R, ω) with state dynamics given by (1.1)-(1.2) and output y = R. This assumes that the matrix of inertia J and the torque τ of the Euler's equations are known. A schematic view of the indirect estimation approach is represented in Figure 1.2.

Solving Wahba a, b

Asymptotic observer The asymptotic observer may be a non-linear observer [Sal91, BPG05, KGK05, TS07, JG11], an EKF [Sun05, [START_REF] Bar-Itzhack | Rigid body rate inference from attitude variation[END_REF][START_REF] Thienel | Comparison of angular velocity estimation methods for spinning spacecraft[END_REF] or one of its many derivatives [CKN97, [START_REF] Oshman | Sequential attitude and attitude-rate estimation using integrated-rate parameters[END_REF][START_REF] Dellus | Estimation of satellite angular velocity using sequential measurements of a single inertial vector[END_REF].

Ṙ = R [ω × ] J ω = Jω × ω + τ R ω a,b
The main drawback of the indirect methods is the preliminary step of solving the Wahba problem. It requires the knowledge of the inertial expressions of the vectors a, b (e.g. the Earth magnetic field or the direction of the Sun in an inertial frame), which are not always available. It relies on possibly heavy algebraic manipulations with non-negligible computational burden. Moreover, the resulting noise model on R is inaccurate. As a result, the "optimality" of the Kalman Filter algorithms cited above may be debated.

Direct methods

Other methods from the literature reconstruct the angular rate from vector measurements directly, namely by bypassing the relatively heavy first step of solving the Wahba problem. These methods do not use the inertial coordinates a, b of the reference vectors. They are generally based on the assumption that the reference vectors are constant, so that the measurements satisfy the differential equation ȧ = a × ω (1.5) Again, we classify these methods using two categories: derivative approaches and estimation approaches.

1. The direct derivative approach consists in solving an algebraic equation stemming from (1.5) to find an estimate of ω. For example, when two non-collinear unit vector measurements are available a, b, then ω satisfies ȧ × a + ḃ × b = P (a, b)ω

where P (a, b) is a symmetric definite positive matrix defined as P (a, b) I -aa + I -bb so that ω = P (a, b) -1 ȧ × a + ḃ × b (1.6) [START_REF] Oshman | Spacecraft angluar velocity estimation using sequential observations of a single directional vector[END_REF] presents another direct derivative approach from a single vector measurements using equation (1.5) and its higher time derivative

ä = ȧ × ω + a × (Jω × ω + τ )
The difficulty of this approach lies in performing the time differentiation of the measurements. This induces high-frequency noise in the value ω obtained by applying (1.6), so that, again, post-filtering is necessary [START_REF] Tortora | Spacecraft angular rate estimation from magnetometer data only using an analytic predictor[END_REF]. A schematic view of the direct derivative approach is represented in Figure 1. 2. The direct estimation approach consists in computing an asymptotic observer for the augmented state S (a, b, ω)

with state dynamics given by (1.5)-(1.2) and output y = (a, b). This assumes that the matrix of inertia J and the torque τ of the Euler's equations are known. A schematic view of the indirect estimation approach is represented in Figure 1.4. To the best of our knowledge, the only direct estimation estimator existing in the literature is an EKF [START_REF] Bar-Itzhack | Classification of algorithms for angular velocity estimation[END_REF]. EKF [START_REF] Bar-Itzhack | Rigid body rate inference from attitude variation[END_REF] direct Algebraic methods [START_REF] Oshman | Spacecraft angluar velocity estimation using sequential observations of a single directional vector[END_REF] EKF post-filtering [START_REF] Tortora | Spacecraft angular rate estimation from magnetometer data only using an analytic predictor[END_REF] EKF [START_REF] Bar-Itzhack | Classification of algorithms for angular velocity estimation[END_REF] Part II: non-linear observer Table 1.2: Classification of the literature and contributions of the thesis The direct methods do not require the knowledge of a, b, which is useful for practical implementation. Yet, they usually consider that a, b are constant. This assumption is acceptable for many applications, but not for all. For the example of an orbiting satellite, the Earth magnetic field significantly varies with the latitude. This is not a problem for short-term estimation, but, obviously becomes problematic for long-term purposes. Kalman Filters (and their extensions) present numerous advantages. They have a general formulation, which is well known and understood. They have been tested extensively and have earned the trust of the scientific community. Yet, proofs of convergence are rare and, when they exist, they are only local. Moreover, they are considered as relatively heavy from a computational perspective 5 . Non-linear observers are way simpler to implement and allow for precise convergence analysis,

The classification we propose is summarized in Table 1.2. Our contributions are indirect derivative methods or direct estimation methods.

Contribution

In this thesis, we propose simple, easy to implement methods to solve the Problem of estimating the angular rate from vector measurements, along with detailed analysis of the error they produce.

We assume that the reference vectors are constant. As R(t) is a rotation matrix, the measurements (1.3) have constant norm, thus lie on a sphere. The difficulty of the estimation problem at stake is mostly defined by the nature of the rotation. The rotation can be single-axis, defining a periodic motion and periodic measurement signals. It can also be slightly varying, generating almost periodic signals, or it can be very complex for an arbitrary rotation. Typical measurements are represented in Figure 1.5 for each of these three cases. In this thesis, we propose three types of solution addressing each of the cases discussed above. This thesis consists in two distinct parts proposing solutions tailored to these cases.

In the first part of the thesis, we provide simple indirect derivative methods in cases where a priori information on the rotation motion is known. In Chapter 3, we consider a single-axis rotation parametrized by a single angle. We give a straightforward estimator for it. The relative simplicity of the underlying equations allows us to give analytical expression of the impact of sensor discrepancies on the estimator. In Chapter 4, we consider the slightly more complex case where the axis of rotation is not constant (e.g. tilting rotation). Some a priori information on its motion is known. We employ a Euler angles parametrization of the attitude. Two cases are treated: i) For a slow drift motion single-axis tilting complex rotation Figure 1.5: Measurement a(t) evolving on a sphere for single-axis (left), slightly varying (middle) and complex (right) rotation motions.

of the axis of rotation, an extension of the single-axis estimator is employed, which relies on the separate analysis of the phase and the modulus of a complex quantity ii) For a small amplitude tilting of the axis of rotation, we propose a two-phase estimator obtained by a careful decomposition of the signals built on the natural frequency separation found in the spectrogram of the measurement signals.

In the second part of the thesis, we use a direct estimation approach and present a non-linear observer of ω, which does not require knowledge of the inertial reference vectors a, b. This is the main contribution of this thesis. In Chapter 5, we make the usual assumptions that the torque τ and the matrix of inertia J are known. We build the observer from two non-collinear vector measurements. We provide a convergence proof and an estimation of the basin of attraction. To our knowledge, no non-linear observer for angular rate exists in the literature. Compared to a Kalman Filter-like estimator, the observer we propose in Chapter 5 does not require solving a differential (or algebraic) Riccati equation, which reduces the computational burden for implementation. It also allows convergence analysis and estimation of the basin of attraction which, to our knowledge, is also new for an angular rate estimator. In Chapter 6, we give a proof of convergence in the case where only one vector is available for measurements. The proof relies on a persistent excitation assumption on the measurement a(t), which is eventually shown to hold in almost all cases. In Chapter 7, we relax the assumption that τ is known and extend the observer to further estimate a torque (assumed to be polynomial in time). Lastly, in Chapter 8, we relax the assumption that J is known and we extend the observer to further estimate the ratios of inertia. The method proposed in this part allows one to estimate the angular rate of a rigid body by using on-board measurements of one or two constant vectors. Interestingly, these vectors need not be known in any inertial frame of reference.

To summarize our approach, it should be noted that the Euler's equations play a central role in all the works developed in this thesis. In the methods of Part 1, they are not directly exploited but grant applicability of frequency analysis. In the methods of Part 2, the Euler's equations directly appear in the augmented state observer governing equations.

As a conclusion of our investigations, we wish to highlight that, at least for the illustrative cases considered here, it appears that rate gyros could be replaced with an estimation algorithm employing direction sensors which are much cheaper, more rugged and more resilient sensors. Further, we believe that this work suggests that, in the classic problem of attitude reconstruction, as performed in e.g. [CBIO06, SRK + 08], or using non-linear observers in [MHP08, MS10, VSO08, TRB11, GFJS12, TMHL12], the rate gyros could be replaced with more in-depth analysis of the vector measurements. This question is left opened.

Note.

For convenience, most of the technical results from the literature used in this thesis are given a proof in Appendix.

The works presented in this thesis are mostly of theoretical nature. However, some experimental results are presented in Appendix.

The works presented in this thesis have been the subject of the following publications.

• Journal papers

Introduction et présentation du problème

Brève description mathématique du problème

Pour comprendre la question scientifique en jeu, nous en ébauchons ici une formulation mathématique. L'orientation du corps rigide est caractérisée par une matrice de rotation R(t) de passage entre un repère inertiel et un repère engin solidaire du corps rigide. Cette matrice de rotation satisfait l'équation différentielle Il est important de savoir que les gyromètres présentent de nombreux désavantages. Ils sont beaucoup plus coûteux que les capteurs de directions, sujets à saturation durant les rotations à grande vitesse [START_REF] Bar-Itzhack | Classification of algorithms for angular velocity estimation[END_REF] et à des dysfonctionnements. La mise en veille du Télescope Spatial Hubble le 13 novembre 1999 suite à une panne de quatre de ses six gyromètres [START_REF] Oshman | Spacecraft angluar velocity estimation using sequential observations of a single directional vector[END_REF] en est un célèbre exemple. Les gyromètres ont été remplacé en 1999 et à nouveau en 2009 lors des très coûteuses missions STS-103 et STS-125 [START_REF] Lallo | Experience with the Hubble Space Telescope: 20 years of an archetype[END_REF].

Ṙ = R    0 -ω 3 ω 2 ω 3 0 -ω 1 -ω 2 ω 1 0    R [ω × ] (1.7) où ω (ω 1 ω 2 ω 3 )
Heureusement, les gyromètres ne sont pas strictement nécessaire dans les applications considérées ici. Une tendance récente consiste à remplacer (ou consolider) les mesures de gyromètres par d'autre capteurs. En particulier, l'estimation de ω à partir de capteurs de direction est à l'étude depuis environ 15 ans. Un exemple de cette tendance "sans gyromètre" est le satellite SAMPEX (Solar Anomalous and Magnetospheric Particle Explorer) lancé en juillet 1992, dont la centrale inertielle ne comprend que des capteurs solaires et des magnétomètres [START_REF] Tsai | SAMPEX spin stabilized mode[END_REF].

Cette thèse suit cette tendance et propose de nouveaux résultats dans ce domaine. Formellement, nous considérons le problème suivant.

Problème 1.1: Problème principal de la thèse

A partir de mesures d'une ou plusieurs directions du type (1.9), trouver une estimation de la vitesse angulaire ω définie par (1.7), sous l'hypothèse qu'elle satisfait (1.8).

Etat de l'art

Plusieurs algorithmes ont été proposés ces dernières années pour répondre au Problème mentionné ci-dessus. L'étude [START_REF] Bar-Itzhack | Classification of algorithms for angular velocity estimation[END_REF] les classifie dans deux catégories : les méthodes dérivatives basées sur une dérivation temporelle des mesures et les méthodes d'estimation qui sont des observateurs asymptotiques pour un état augmenté. Dans cette thèse, nous poursuivons la classification en considérant les méthodes indirectes qui requièrent une manipulation algébrique préalable des mesures pour exprimer l'attitude, par opposition aux méthodes directes. En raison du chevauchement des deux dichotomies, nous identifions finalement quatre types de méthodes détaillés ci-dessous.

Méthodes indirectes

De nombreuses méthodes d'estimation de ω considèrent que la matrice d'attitude R est connue6 . Implicitement, ceci suppose que les mesures de direction ont déjà été traitées pour produire R. En supposant connus au moins deux vecteurs de référence non-colinéaires a, b et leur expression (mesures) dans le repère engin a, b, la matrice de rotation peut être déterminée comme solution du célèbre problème de Wahba [START_REF] Wahba | Problem 65-1: a least squares estimate of spacecraft attitude[END_REF]. C'est un problème de minimisation portant sur R. De nombreuses solution au problème de Wahba ont été proposées comme QUEST ou TRIAD [START_REF] Shuster | Approximate algorithms for fast optimal attitude computation[END_REF][START_REF] Shuster | Kalman filtering of spacecraft attitude and the QUEST model[END_REF]. Indépendamment de la méthode employée, nous poursuivons maintenant avec la classification établie dans [START_REF] Bar-Itzhack | Classification of algorithms for angular velocity estimation[END_REF].

Résolution Wahba

a, b

Observateur asymptotique Le principal désavantage des méthodes indirectes est l'étape préliminaire de résolution du problème de Wahba. Elle requière la connaissance des coordonnées inertielles des vecteurs de référence a, b (par exemple le champ magnétique terrestre ou la direction du soleil dans un repère inertiel), qui ne sont pas toujours disponibles. Elle repose sur des manipulations algébriques des mesures nécessitant un temps de calcul potentiellement long. En outre, le modèle de bruit qui en résulte sur R n'est pas précis. En conséquence, l'"optimalité" du filtre de Kalman discutée plus haut peut être remise en question.

Ṙ = R [ω × ] J ω = Jω × ω + τ R ω a,b

Méthodes directes

D'autres méthodes de la littérature reconstruisent la vitesse angulaire directement à partir de mesures de direction, c'est-à-dire en court-circuitant l'étape préalable relativement lourde de résolution du problème de Ainsi, ω est déterminée de manière univoque par la formule

ω = P (a, b) -1 ȧ × a + ḃ × b (1.12)
[OD03] présente une autre approche dérivative directe à partir d'une seule mesure de direction en utilisant l'équation (1.11) et sa dérivée

ä = ȧ × ω + a × (Jω × ω + τ )
La difficulté de cette approche réside dans la différentiation temporelle des mesures. Celle-ci induit un bruit à haute fréquence dans la quantité ω obtenue depuis (1.12) de sorte que, à nouveau, un post-filtrage est nécessaire [START_REF] Tortora | Spacecraft angular rate estimation from magnetometer data only using an analytic predictor[END_REF]. Une vue schématique de l'approche dérivative directe est représentée sur la Figure 1.8. Observateur asymptotique

ω = P (a, b) -1 ȧ × a + ḃ × b Filtre ω a, b
ȧ = a × ω ḃ = b × ω J ω = Jω × ω + τ ω a, b

Contribution

Dans cette thèse, nous proposons des méthodes simples et faciles à implémenter pour résoudre le Problème d'estimation de vitesse angulaire par mesures de direction, ainsi qu'une analyse détaillée de l'erreur d'estimation associée.

Nous faisons l'hypothèse que les vecteurs de référence sont constants. Comme R(t) est une matrice de rotation, les mesures (1.9) ont une norme constante, donc évoluent sur une sphère. La difficulté du problème d'estimation en jeu est principalement définie par la nature de la rotation. Elle peut-être à axe fixe, à axe légèrement variable ou bien quelconque. Des exemples typiques de mesures sont représentées sur la Figure 1.10 pour chacun de ces trois cas. Dans cette thèse, nous proposons trois types de solutions permettant de traiter ces cas. Elle est constituée de deux parties.

Dans la première partie de cette thèse, nous présentons des méthodes dérivatives indirectes simples dans des cas où on connait une information a priori sur le mouvement de rotation. Dans le Chapitre 3, nous considérons une rotation à axe fixe, paramétrée par un angle. Nous construisons un estimateur statique de cet angle. La simplicité des équations sous-jacentes nous permet de donner une expression analytique de l'impact des défauts de capteurs sur l'estimateur. Dans le Chapitre 4, nous considérons le cas d'un axe de rotation légèrement perturbé. Nous employons une paramétrisation de l'attitude par les angles d'Euler. Deux sous-cas du mouvement de l'axe sont envisagés : i) pour un mouvement de dérive lente, nous employons une extension de l'estimateur à axe fixe qui repose sur les études séparées de la phase et du module d'un nombre complexe bien choisi ii) pour une oscillation à faible amplitude, nous proposons un estimateur de phases qui repose sur une séparation naturelle de deux fréquences caractéristiques du mouvement contenues dans le spectrogramme des mesures.

Dans la seconde partie de la thèse, nous utilisons une approche d'estimation directe et présentons un observateur non-linéaire de ω ne nécessitant pas la connaissance des coordonnées inertielles des vecteurs de référence a, b. C'est la contribution principale de axe fixe axe perturbé rotation quelconque Figure 1.10: Mesure a(t) évoluant sur une sphère pour un mouvement de rotation à axe fixe (gauche), à axe légèrement perturbé (milieu) ou bien quelconque (droite). cette thèse. Dans le Chapitre 5, nous faisons l'hypothèse usuelle que la matrice d'inertie J et le couple τ sont connus. Nous construisons l'observateur à partir de deux vecteurs de référence non-colinéaires. Nous établissons une preuve de convergence et une estimation du bassin d'attraction. A notre connaissance, aucun observateur de ce type n'existe dans la littérature. Comparé à un estimateur de type filtre de Kalman, l'observateur que nous proposons ne nécessite pas de résoudre une équation différentielle (ou algébrique) de Riccati, ce qui réduit le temps de calcul en vue d'une implémentation. Il permet en outre une analyse de convergence et une estimation du bassin d'attraction ce qui, à notre connaissance, est aussi nouveau pour un estimateur de vitesse angulaire. Dans le Chapitre 6, nous donnons une preuve de convergence dans le cas où un seul vecteur de référence est disponible. La preuve repose sur une hypothèse d'excitation persistante sur la mesure a(t). Nous démontrons que cette hypothèse est valide pour presque toutes les conditions initiales. Dans le Chapitre 7, nous écartons l'hypothèse de couple connu et étendons l'observateur pour estimer en outre un couple (supposé polynomial en temps). Enfin, dans le Chapitre 8, nous écartons l'hypothèse que J est connu et étendons l'observateur pour estimer en outre les ratios d'inertie. La méthode proposée dans cette partie permet d'estimer la vitesse angulaire d'un corps rigide en rotation en utilisant des mesures embarquées d'un ou deux vecteurs de référence. De manière intéressante, la connaissance des coordonnées inertielles de ces vecteurs n'est pas nécessaire.

Il importe de préciser que les équations d'Euler jouent un rôle central dans les travaux développés dans cette thèse. Pour les méthodes de la Partie 1, elles ne sont pas exploitées directement mais garantissent l'applicabilité de l'analyse de fréquences. Dans les méthodes de la Partie 2, elles apparaissent directement dans les équations différentielles régissant la dynamique de l'état augmenté.

En conclusion de nos investigations, nous souhaitons souligner que, du moins pour les illustrations considérés ici, les gyromètres pourraient être remplacés par un algorithme d'estimation utilisant des capteurs de direction qui se trouvent être moins chers et plus résistants. En outre, nous pensons que, dans le problème classique d'estimation d'attitude tel que réalisé dans [CBIO06, SRK + 08], ou bien par des observateurs non-linéaires dans [MHP08, MS10, VSO08, TRB11, GFJS12, TMHL12], les gyromètres pourraient être remplacés par une analyse approfondie des capteurs de direction. Cette question reste ouverte.

Note Par commodité, des preuves de la plupart des résultats techniques provenant de la littérature utilisés dans cette thèse sont fournies en Appendice.

Les travaux présentés dans cette thèse sont principalement de nature théorique. Toutefois, des résultats expérimentaux sont présentés en Appendice.

Les travaux présentés dans cette thèse ont fait l'objet des publications suivantes. Chapter 2

Mathematical formulation and notations

Formulation mathématique et notations. Dans ce chapitre, nous introduisons en détail les notions mathématiques et les notations utiles à la compréhension du problème presenté dans l'introduction.

We now present a mathematical formulation of the problem discussed in the introduction.

Inertial frame, body frame

The ambient space is the usual three dimensional Euclidean real space R 3 equipped with an inertial frame R i = e 1 , e 2 , e 3 (2.1)

Any vector x ∈ R 3 is defined by its coordinates in R i x =    x 1 x 2 x 3   
x 1 e 1 + x 2 e 2 + x 3 e 3 where, by definition,

e 1 =    1 0 0    , e 2 =    0 1 0    , e 3 =    0 0 1   
We consider an orthonormal body frame

R b (b 1 , b 2 , b 3 ) (2.2)
attached to a rigid body. For any vector x ∈ R 3 , we denote

x    x 1 x 2 x 3    (2.3)
Chapter 2. Mathematical formulation and notations the coordinates of x in R b , so that

x = x 1 b 1 + x 2 b 2 + x 3 b 3
and we will loosely write that "x is x expressed in R b ". The motion of the rigid body has six degrees of freedom, three for the translation of its mass center, and three for its attitude. In this thesis, we are only interested in the attitude, which is unambiguously characterized by the rotation matrix R from R i to R b defined as

Re 1 = b 1 , Re 2 = b 2 , Re 3 = b 3
For any x ∈ R 3 , the coordinates of x in R i and R b are related by the following formula we have

x = R x
The rotation matrix R(t), which we assume to be a smooth function of time t, evolves in the Lie-group

SO 3 M ∈ R 3×3 , M M = I, det M = 1
The tangent space to SO 3 at point M is

T M (SO 3 ) = H ∈ R 3×3 , M H + H M = 0 = M H, H = -H
It follows that, at all times, the derivative Ṙ(t), which belongs to T R(t) (SO 3 ), can be written

Ṙ(t) = R(t)H(t) (2.4)
where H(t) is a skew-symmetric 3 × 3 matrix, i.e. H(t) can be written

H(t) =    0 -ω 3 (t) ω 2 (t) ω 3 (t) 0 -ω 1 (t) -ω 2 (t) ω 1 (t) 0    , with ω 1 (t), ω 2 (t), ω 3 (t) ∈ R
We denote

ω(t) =    ω 1 (t) ω 2 (t) ω 3 (t)    (2.
d dt R 0 R = R 0 Ṙ = R 0 R [ω × ]
Thus, the corresponding angular rate is again ω. This means that ω does not depend on the choice of the inertial frame.

Euler's equations for a rotating rigid body

The motion of ω itself is governed by the famed Euler's equations

1 [LL82] ω = J -1 (Jω × ω + τ ) (2.7)
where J is the matrix of inertia of the rigid body, expressed in R b , and τ is the external torque applied to the rigid body, also expressed in R b . For brevity we denote E(ω) J -1 (Jω × ω) , χ J -1 τ so that the Euler's equations (2.7) are simply written as

ω = E(ω) + χ (2.8)
If, without loss of generality, the axes of R b are aligned with the principal axes of inertia of the rigid body, then the matrix J is diagonal

J =    J 1 0 0 0 J 2 0 0 0 J 3   
(2.9)

The parameters J 1 , J 2 , J 3 are the principal moments of inertia of the rigid body. Without loss of generality, we will always assume

J 1 ≥ J 2 ≥ J 3 > 0 (2.10)
With these notations, the expression of E(ω) can be simplified into

E(ω) =    J 2 -J 3 J 1 ω 2 ω 3 J 3 -J 1 J 2 ω 3 ω 1 J 1 -J 2 J 3 ω 1 ω 2   
The coefficients

d 1 J 2 -J 3 J 1 ≥ 0, d 2 J 3 -J 1 J 2 ≤ 0, d 3 J 1 -J 2 J 3 ≥ 0 (2.11)
are the ratios of inertia of the rigid body. They represent its inertial asymmetry. We call discordance of the rigid body the quantity d max defined as

d max max(d 1 , |d 2 |, d 3 ) (2.12)
The discordance is zero if and only if J 1 = J 2 = J 3 , i.e. the rigid body is an inertial sphere. The torque τ (and equivalently χ) gathers the contributions from control inputs and disturbances2 . In the absence of torques, the Euler's equations reduce to

d dt    ω 1 ω 2 ω 3    = E(ω) =    d 1 ω 2 ω 3 d 2 ω 3 ω 1 d 3 ω 1 ω 2   
(2.13) and are called free Euler's equations. The rigid body is said to be in free-rotation.

Angular momentum

The angular momentum of the rigid body is a vector

M ∈ R 3 defined by M RJω (2.14) Its expression in R b is R M = Jω.
If the rigid body is in free-rotation, we have

Ṁ = ṘJω + RJ ω = R [ω × ] Jω + RJω × ω = R (ω × Jω + Jω × ω) = 0
so that M is constant. As a result, the trajectory of ω lies on the surface of the ellipsoid defined by

|Jω| 2 = |M| 2 ⇔ ω J 2 ω = |M| 2 (2.15)
The shape of the ellipsoid only depends on the norm |M|.

Measurements

In this thesis, we consider that the rigid body is equipped with sensors measuring threedimensional vectors expressed in R b . Such measurements are called vector observations or vector measurements. In details, consider a vector a ∈ R 3 . Then, the corresponding measurement is, up to sensor discrepancies, the expression of a in R b , namely at time t

a(t) = R(t) a (2.16)
In the specific case where a is constant, then we have

ȧ = Ṙ a = -[ω × ] R a = -ω × a = a × ω
so that a satisfies the linear time-varying (LTV) differential equation ȧ(t) = a(t) × ω(t) (2.17)

Notations relative to rigid body dynamics

For vector or matrix quantities, the unit applies to each coordinate.

notation meaning unit R i = (e 1 , e 2 , e 3 ) inertial frame - R b = (b 1 , b 2 , b 3 ) body frame, orthonormal base of R 3 - x, x coordinates of a generic vector of R 3 in R i , R b - R rotation matrix from R i to R b , i.e. Re i = b i - ω angular rate expressed in R b [rad/
u × v cross product of vectors u, v ∈ R 3 [u × ] cross product matrix of vector u ∈ R 3 , i.e. ∀v ∈ R 3 , [u × ] v = u × v S 2
unit vectors of R 3 R n×n set of real matrices of size n × n X transpose of matrix or vector X SO 3 set of rotations matrices of size 3 × 3 r u (ψ) rotation matrix of axis u ∈ S 2 and angle ψ ∈ R I identity matrix, the size is understood from the context C field of complex numbers z, z real and imaginary parts of a complex number z z * conjugate of a complex number z In this chapter, we consider the case of single-axis (or planar) rotations. Such rotations correspond to a matrix R(t) of the form

arg -π z argument determination in [-π, π) of a complex number z z 1 • z 2 dot product of complex numbers z 1 , z 2 , i.e. z 1 • z 2 = z 1 z 2 + z 1 z 2 D(z,
R(t) = r u (ψ(t)) cos ψ(t)I + sin ψ(t) [u × ] + (1 -cos ψ(t))uu
where u ∈ S 2 is a constant unit vector and ψ(t) ∈ R is the angle of the rotation (see Appendix A.1 for more details). The corresponding angular rate is

ω(t) = ψ(t)u
The measurement can be written

a(t) = R(t) a = r u (-ψ(t))a
If a(t) is constant over time it does not provide any information on the rotation. This can happen only in the two (non-exclusive) following cases:

• u and a are collinear, in which case a(t) = a for all t • ψ is constant over time, in which case a(t) = r u (-ψ)a for all t Leaving out these cases, in the rest of the chapter we assume that a and u are linearly independent and that ψ varies.

Mathematical formulation, problem statement

Decoupling the information on u and ψ. The coordinates of a(t) are a parametric representation of a curve contained in a plane orthogonal to u as represented in Figure 3 Thus, the direction of u is unambiguously defined. Without loss of generality, we consider that u is the vertical inertial vector, namely u = (0 0 1) and we have

a(t) = cos ψ(t) a -sin ψ(t)    0 0 1    × a + (1 -cos ψ(t))a 3    0 0 1    so that a 1 (t) = cos ψ(t) a 1 + sin ψ(t) a 2 a 2 (t) = -sin ψ(t) a 1 + cos ψ(t) a 2 a 3 (t) = a 3
Thus, all information on the angle is contained in the first two coordinates of the measurements. For convenience, we gather them into a single complex-valued variable

y(t) a 1 (t) -ia 2 (t) = e iψ (a 1 -ia 2 ) ∈ C (3.1)
where i is the usual unit imaginary number. The natural estimate ψ of ψ

ψ = arg -π y
where arg -π is the argument determination in [-π, π), determines ψ up-to the argument of a 1 -ia 2 , which is constant.

Sensor discrepancies, sampling, noises Formula (3.1) does not account for sensor discrepancies such as biases and non-linearities 1 . In general, the measurement equation writes y(t) = f (ψ(t)) where f is a 2π-periodic complex function. This recasts the problem of estimating ψ(t) as the phase estimation of a 2-dimensional vector. The signal is sampled at frequency ν s = 1 ∆t and corrupted by noises. As a result, the measurements can be written

y[k] = f (ψ[k]) + n[k] ∈ C, 1 ≤ k ≤ N (3.2)
where 1 see Appendix D.1 for an experimental example on Sun sensors

• f is an (unknown) 2π-periodic function valued in C (or equivalently in R 2 ) • ψ[k] ψ(k∆t) • n[k] is a measurement noise Problem 3.1
Find an estimate of ψ from measurements of the form (3.2).

A naive estimate ψ[k] of ψ[k] is simply ψ[k] = arg -π y[k]
Yet, to account for sensor bias, it might be useful to estimate the argument with respect to an origin z 0 which is not necessarily 0. Also, it could be desirable to estimate the cumulative angle (e.g. to count turns), and not only its value in [-π, π). A simple solution for this is to introduce the following estimator

ψz 0 [k] = k-1 j=1 arg -π y[j + 1] -z 0 y[j] -z 0 (3.3)
for a given choice of an origin z 0 . Choosing z 0 will be addressed in Section 3.3.3.

Remark.

In the following, we assume that the first sample ψ[1] = 0, so that the estimation results are to be understood up to a constant value.

This chapter proposes an analysis of the estimation error produced by (3.3).

e z 0 [k] = ψz 0 [k] -ψ[k] (3.4)
In Section 3.2, we first consider an idealized case without noise. Starting with a continuous phase description, and introducing sampling, we lay the basis of the analysis. In Section 3.3, we take the noises into account and quantify its impact on the estimation error. Simulation results are provided in Section 3.5. The results presented here can serve as a basis for an estimation of the angular rate ω = ψu.

Estimation principle for a periodic signal

In this section, we consider a case without any noise and study general properties of periodic signals. We start by introducing some notations and by defining the phase variation of f around an origin in Section 3.2.1. In Section 3.2.2 we see how this phase variation gives an estimate of ψ as it varies continuously. Finally, we relate this study with the definition of ψ[k] in Section 3.2.3.

Preliminaries and notations

We introduce handy notations and recall some basic complex analysis results (as exposed in details in [START_REF] Rudin | Real and complex analysis[END_REF] e.g.). Here, the ambient space is the complex plane C counterclockwise orientated, its origin is denoted by O.

We assume that f is continuously differentiable. For any ψ 1 , ψ 2 ∈ R, we denote by

C [ψ 1 ,ψ 2 ] the curve described by f (ζ) for ζ ∈ [ψ 1 , ψ 2 ].
The closed curve C [0,2π] is simply denoted by C. To simplify, we assume that C is a Jordan curve, i.e. is non-self-intersecting. Thus, by the classic Jordan curve theorem [START_REF] Hales | The Jordan curve theorem, formally and informally[END_REF], C separates C in two regions (connected components). We denote by I the (interior) bounded region.

For any ψ ∈ R and any z 0 ∈ C\C, the phase variation P of f (ζ) with respect to z 0 as 0 ≤ ζ ≤ ψ is defined as

P (ψ, z 0 ) = C [0,ψ] dz z -z 0 ψ 0 f (ζ) f (ζ) -z 0 dζ (3.5)
where designates the imaginary part. This definition is illustrated in Fig. 3.2. By defini-

P (ψ, z 1 ) f (ψ) z 0 ∈ I z 1 / ∈ I f (0) P (ψ, z 0 ) Figure 3.2: Phase variations around origins z 0 or z 1 when f (ζ) describes C [0,ψ] .
tion, P (2π, z 0 ) is equal to 2πI z 0 , where I z 0 ∈ Z is the algebraic number of counterclockwise rotation of C about z 0 . Namely, I z 0 = 0 if and only if z 0 ∈ I, and then it equals to 1 or -1. By assumption, we consider that I z 0 = 1 for all z 0 ∈ I throughout the chapter, i.e. that C is positively oriented.

Remark 3.1. For an oriented segment

[z 1 , z 2 ] and z 0 / ∈ [z 1 , z 2 ], we have [z 1 ,z 2 ] dz z -z 0 = arg -π z 2 -z 0 z 1 -z 0
We will abundantly use this result in the rest of this chapter.

Estimate of a continuous phase

We take z 0 ∈ I. We will see that P (ψ, z 0 ) gives an estimate of ψ under some assumptions on f . For now, let us denote the error variable

e z 0 (ψ) = P (ψ, z 0 ) -ψ (3.6)
and derive bound for its magnitude. Let {c n } n∈Z be the coefficients of the Fourier expansion of f . We define

g(ζ) exp(-iζ) (f (ζ) -z 0 ) which satisfies g(ζ) = c 1 + n =0,1 c n e i(n-1)ζ + (c 0 -z 0 )e -iζ
Using this new function, we can directly prove the following result.

Proposition 3.1. For any z 0 ∈ I,

e z 0 (2nπ) = 0, ∀n ∈ Z (3.7)
and e z 0 (ψ) is bounded by

|e z 0 | ∞ = max 0≤ψ≤2π ψ 0 g (ζ) g(ζ) dζ (3.8)
Proof. For any ψ ∈ R, we have

P (ψ, z 0 ) = ψ 0 i + g (ζ) g(ζ) dζ = ψ + ψ 0 g (ζ) g(ζ) dζ
Since Z 0 ∈ I and we have P (2π, z 0 ) = 2π which shows that

2π 0 g (ζ) g(ζ) dζ = 0
The result follows immediately.

Equation (3.7) means that the estimate matches ψ perfectly, at least once every cycle. It is illustrated on experimental data in Appendix D.1. The error bound (3.8) states that the error is bounded by the phase variation of the curve g with respect to the origin O. It is of theoretical but of no practical value. To derive a more concrete bound, we can make further assumptions on g.

Let us assume that there exists r < |c 1 | such that the set {g(ζ), ζ ∈ [0; 2π]} lies in the closed circle D(c 1 , r) of center c 1 and radius r. Then, as is illustrated in Fig. 3.3, for any ψ the angle g(0)Og(ψ) is included in the angle AOB, whose value is 2 arcsin r |c 1 | . Thus,

|e z 0 (ψ)| is bounded by |e z 0 | ∞ ≤ 2 arcsin r |c 1 | (3.9)
This inequality allows us to directly derive the following result, obtained with

r = |g -c 1 | ∞ . O B A c 1 r g(0) g(ψ)
Figure 3.3: Bound on the estimation error when O lies outside a disk containing the curve of g.

Proposition 3.2. If f is such that its Fourier expansion satisfies

|c 1 | > n =0,1 |c n | + |c 0 -z 0 | (3.10)
then one has

|e z 0 | ∞ ≤ 2 arcsin n =0,1 |c n | + |c 0 -z 0 | |c 1 | (3.11)
Example 3.1. This property stems from the "almost circular" nature of the parametric curve defined by the measurements. It is an implicit consequence of assumption (3.10).

To stress this, consider a positively oriented ellipse with eccentricity e between 0 and 1

f (ψ) = A cos ψ + i √ 1 -e 2 A sin ψ
In this case, f has only two Fourier coefficients

c 1 = 1 + √ 1 -e 2 2 , c -1 = 1 - √ 1 -e 2 2 As 0 < c -1 < c 1 , condition (3.10) is met for z 0 = O and we have |e O | ∞ ≤ 2 arcsin 1 - √ 1 -e 2 1 + √ 1 -e 2 = e 2 2 + O(e 4 )
Even for a flat ellipse with e = 0.99, this estimation error does not exceed π 2 radians.

Sampling

We now take into account the fact that measurements are only available in discrete time.

To study this fact, we consider a noise-less version of (3.2). To estimate ψ[k], we apply the preceding method, to the polygonal line joining the vertices

f (ψ[1]) . . . f (ψ[k]
), which are measured. The case of noisy measurements is considered later in Section 3.3. Here, we will see that the integral term (3.5) is exactly ψ[k] as defined in (3.3). To formalize this, let us first introduce a few notations specific to this section.

For any 1

≤ k ≤ N -1, we simply denote by γk = C [ψ[k],ψ[k+1]] and γ k the oriented segment [f (ψ[k + 1]), f (ψ[k])]
, so that the concatenation γ k • γk is a closed curve. We also denote by Ω k the unbounded region defined by this closed curve. These notations are illustrated in Fig. 3.4. The next result states that, if z 0 is chosen in the intersection of all the Ω k , then the integral (3.5) can be calculated from the sampled values, without loss of accuracy.

z 0 ∈ Ω k γ k f (ψ[k]) f (ψ[k + 1]) z 1 / ∈ Ω k γk
Proposition 3.3. Assuming that z 0 satisfies

z 0 ∈ N -1 k=1 Ω k ∩ I (3.12) then for all 2 ≤ k ≤ N , P (ψ[k], z 0 ) = k-1 j=1 arg -π f (ψ[j + 1]) -z 0 f (ψ[j]) -z 0 (3.13) Proof. Let k be any 1, . . . , N -1. The condition z 0 ∈ Ω k implies that 0 = γ k dz z -z 0 + γk dz z -z 0
Using the formula recalled in Remark 3.1, we have

γk dz z -z 0 = - γ k dz z -z 0 = arg -π f (ψ[k + 1]) -z 0 f (ψ[k]) -z 0
Similarly, for 2 ≤ k ≤ N , the calculations above are readily generalized to

P (ψ[k], z 0 ) = k-1 j=1 γj dz z -z 0 = k-1 j=1 arg -π f (ψ[j + 1]) -z 0 f (ψ[j]) -z 0 Proposition 3.3 guarantees that, in the absence of noise i.e. n[k] = 0 for 1 ≤ k ≤ N , the estimate error is exactly e z 0 [k] = e z 0 (ψ[k]
) as defined in (3.4). Hence, the bounds (3.8)-(3.9)-(3.11) on e z 0 listed in Section 3.2.2 also apply on e z 0 [k].

Impact of measurement noise

Critical value of the noise magnitude

We now consider the measurement equation (3.2) in its full version, i.e. with noise. For any 1

≤ k ≤ N -1, we denote by Q k the (non-necessarily convex) quadrilateral with ordered vertices f (ψ[k]), y[k], y[k + 1], f (ψ[k + 1]).
The next result provides a bound on the estimate error.

Proposition 3.4. Assuming that z 0 satisfies condition (3.12) and

z 0 / ∈ Q k , ∀k ∈ {1, . . . , N -1} (3.14) then for all 2 ≤ k ≤ N , e z 0 [k] = e z 0 (ψ[k]) + arg -π y[k] -z 0 f (ψ[k]) -z 0 -arg -π y[1] -z 0 f (ψ[1]) -z 0 (3.15) z 0 / ∈ Q k γ k f (ψ[k]) f (ψ[k + 1]) Q k γk y[k + 1] y[k] γ n- k+1 γ n k Figure 3.5: Q k 's boundary is the closed curve γ n k • γk • γ n- k+1 • γ k . Proof. For any 1 ≤ k ≤ N -1, denote • γ n k the oriented segment [f (ψ[k]), y[k]], and γ n- k the opposed oriented segment [y[k], f (ψ[k])] • γk the oriented segment [y[k], y[k + 1]]
These notations are illustrated in Fig. 3.5. As z 0 / ∈ Q k , we have

γ n k dz z -z 0 + γk dz z -z 0 + γ n- k+1 dz z -z 0 + γ k dz z -z 0 = 0
Let us now consider k ∈ {2, . . . , N } and let us sum the terms of this equation for 1 ≤ j ≤ k -1. This yields

γ n 1 dz z -z 0 + γ1 •...•γ k-1 dz z -z 0 + γ n- k dz z -z 0 + γ 1 •...•γ k-1 dz z -z 0 = 0
Integrating dz z-z 0 on the polygonal lines and taking the imaginary part of the obtained result, we get

ψz 0 [k] = k-1 j=1 arg -π f (ψ[j + 1]) -z 0 f (ψ[j]) -z 0 + arg -π y[k] -z 0 f (ψ[k]) -z 0 + arg -π y[1] -z 0 f (ψ[1]) -z 0
Under the assumption (3.12), we have

k-1 j=1 arg -π f (ψ[j + 1]) -z 0 f (ψ[j]) -z 0 = e z 0 (ψ[k]) + ψ[k]
The result follows immediately.

We now wish to define a set of allowed origins in the sense of (3.12) and (3.14), independent of the values of ψ[k]. To do so, we introduce ρ the maximum noise amplitude and ∆ the maximum ψ variation between two consecutive samples such that for all k,

|y[k] -f (ψ[k])| ≤ ρ and |ψ[k + 1] -ψ[k]| ≤ ∆. It is assumed that ∆ < π.
For any ψ, ζ, let E(ψ, ζ, ρ) be the set defined as

E(ψ, ζ, ρ) = z∈[f (ψ),f (ζ)]

D(z, ρ)

We can now define the set of forbidden origins

F ρ,∆ = 2π ψ=0 ψ+∆ ζ=ψ-∆ E(ψ, ζ, ρ) (3.16)
which contains all the Q k and Ω c k , and its complement, the set of allowed origins

A ρ,∆ = I\F ρ,∆ (3.17)
depending only on ρ and ∆. For z 0 ∈ I, we define

δ(z 0 ) min ζ |f (ζ) -z 0 | If δ max max z 0 ∈I δ(z 0 ) is less than ρ, then A ρ,∆ is
δ max = A √ 1 -e 2 .
Hence, regardless of sampling, the minimum signal-to-noise ratio (SNR) allowed is:

SNR min = 10 log A √ 1 -e 2 A 2 = -10 log(1 -e 2 ) [dB]
For eccentricities e = 0.2, 0.7, or 0.99, we have respective minimum SNR of 0.4, 7, and 39 [dB].

Study of the set of allowed origins

For any z ∈ C, we denote by z * its complex conjugate, z its real part and z its imaginary part. For any z 1 , z 2 ∈ C one has

z * 1 z 2 = z 1 • z 2 , z * 1 z 2 = det[z 1 , z 2 ]
where • designates the usual dot product in the plane. We will extensively use these notations in this section.

Proposition 3.5. If I is a convex set, the set of forbidden origins (3.16) can simply be expressed as

F ρ,∆ = 2π ψ=0 E(ψ, ψ + ∆, ρ) (3.18) v(ψ) u(ψ) ρ ρ f (ψ + ∆) f (ψ) f (ψ + ∆) + ρv(ψ) f (ψ) + ρv(ψ) f (ψ + ∆) -ρv(ψ) f (ψ) -ρv(ψ) E(ψ, ψ + ∆, ρ) D - ψ D ψ D - ψ Figure 3
.6: Notations for the proof of Proposition 3.5.

Proof. To shorten the derivation of this result, we simply denote

E ψ = E(ψ, ψ + ∆, ρ)
By definition of F ρ,∆ , we have

2π ψ=0 E ψ ⊂ F ρ,∆
We now prove the reciprocal inclusion. For any ψ, any ζ ∈ [ψ -∆, ψ + ∆] and any

w ∈ E(ψ, ζ, ρ), one will prove that w ∈ 2π ψ=0 E ψ . Case 1: if ζ is in [ψ, ψ + ∆]. As ∆ < π, f (ψ) = f (ψ + ∆).
We define the orthogonal unit vectors u(ψ) and v(ψ) as

u(ψ) = f (ψ + ∆) -f (ψ) |f (ψ + ∆) -f (ψ)| , v(ψ) = e i π 2 u(ψ)
We denote by D ψ the line going through f (ψ) and f (ψ + ∆). We also denote the parallel lines

D ± ψ D ψ ± ρv(ψ).
Those definitions are illustrated in Fig. 3.6. Respective equations of these lines are

v(ψ) * (z -f (ψ)) = 0 v(ψ) * (z -f (ψ)) = ρ v(ψ) * (z -f (ψ)) = -ρ
The set P ψ lying between D + ψ and D - ψ is characterized by

-ρ ≤ v(ψ) * (z -f (ψ)) ≤ ρ
If w ∈ P ψ , one can easily show that either w ∈ E ψ , or δ(w) < ρ. In both cases, we have

w ∈ 2π ψ=0 E ψ
We now prove that there exists ξ ∈ [ψ, ζ] so that w ∈ P ξ , which, according to the preceding remark, concludes the proof under Case 1. As I is convex, C is positively oriented and v(ψ) is the positive normal vector to the oriented segment [f (ψ), f (ψ + ∆)], we have

v(ψ) * (f (ζ) -f (ψ)) ≤ 0 Thus, for any z ∈ D(f (ζ), ρ) we have v * (ψ)(z -f (ψ)) ≤ v * (ψ)(z -f (ζ) ≤ |v(ψ)||z -f (ζ)| ≤ ρ
We also have for any

z ∈ D(f (ψ), ρ), v * (ψ)(z -f (ψ)) ≤ ρ. Hence, the set D(f (ζ), ρ) ∪ D(f (ψ), ρ)
lies in a convex closed half-plane. As any element of its convex hull, w also satisfies

v * (ψ)(w -f (ψ)) ≤ ρ
Applying the same reasoning to ψ, ζ and ζ + ∆ shows that w satisfies

v * (ζ)(w -f (ζ)) ≥ -ρ
If w belongs to P ψ ∪ P ζ , we have nothing to prove. Otherwise, both following inequalities

hold v * (ψ)(w -f (ψ)) < -ρ, v * (ζ)(w -f (ζ)) > ρ
The intermediate value theorem implies that the continuous real function

ξ → v * (ψ)(w -f (ξ))
reaches the value 0 for some ψ ∈ [ψ, ζ], and w ∈ P ψ .

Case 2: if ζ ∈ [ψ -∆, ψ[, we have w ∈ E(ζ, ψ, ρ) with ψ ∈ [ζ, ζ + ∆]
. Thus, we can apply the result of Case 1. This establishes the desired inclusion and concludes the proof.

Thanks to Proposition 3.5, we can give a more intuitive expression of the set of allowed origins. Indeed, denote

S 2π ψ=0 {z ∈ C, v * (ψ)(z -f (ψ)) > ρ}. One can prove that I ∩ S = I\F ρ,∆ . The only difficulty is to establish that I\S ⊂ I ∩ F ρ,∆ . Let z be in I\S. For any ψ ∈ [0, 2π], z is on the same side of D ψ as f (ψ) -ρv(ψ). The line going through f (ψ) and z intersects C on a point f (ζ) with ψ < ζ < ψ + ∆, thus z ∈ E(ψ, ζ, ρ) ⊂ F ψ,ρ .
We have proven the following result. Proposition 3.6. If I is a convex set, we have

A ρ,∆ = I ∩ 2π ψ=0 {z ∈ C, v * (ψ)(z -f (ψ)) > ρ} (3.19)
In particular, A ρ,∆ is also a convex set.

Formula (3.19) helps us understand the set of allowed origins. In particular, one can show that it inherits the axial symmetry properties of f . If f has two axis of symmetry, thus a center of symmetry C at their intersection, and if the convex set A ρ,∆ is not empty, then it contains C. Thus, a simple criterion on (ρ, ∆) for A ρ,∆ to be not empty is that

C satisfies v(ψ) * (C -f (ψ)) > ρ, ∀ψ
For example, if f describes a circle of radius A, we easily deduce that A ρ,∆ is not empty if and only if ρ < A cos ∆ 2 More generally, for the ellipse of Example 3.1, the criterion becomes

ρ < A √ 1 -e 2 cos ∆ 2
Let us now assume that A ρ,∆ is not empty and that I is strictly convex. Thus, the phase of v(ψ) is strictly increasing with ψ and (v(ψ), v (ψ)) is a direct frame of reference of the plane. Namely, v(ψ) * v (ψ) > 0. A ρ,∆ is also strictly convex and its boundary can be derived from its tangent lines. Therefore, it is included in (but not necessarily equal to)

2π ψ=0 lim δψ→0,δψ =0 D + ψ ∩ D + ψ+δψ
The next result gives an explicit expression of this boundary.

Proposition 3.7. If A ρ,∆ is not empty and I is strictly convex, then A ρ,∆ is an interior region defined by the (not necessarily Jordan) closed curve

h : [0, 2π] ψ → f (ψ) + ρv(ψ) + i v(ψ) * f (ψ) |v (ψ)| v(ψ) (3.20)
Proof. Let ψ be in [0, 2π]. For δψ = 0 sufficiently small, v(ψ) and v(ψ+δψ) are independent, namely v(ψ) * v(ψ + δψ) = 0. Thus D + ψ and D + ψ+δψ have a unique intersection z(ψ, ψ+δψ),

solution of v(ψ) * (z -f (ψ)) = ρ = v(ψ + δψ) * (z -f (ψ + δψ))
which gives

z(ψ, ψ + δψ) = f (ψ) + ρv(ψ) + i v(ψ) * (f (ψ + δψ) -f (ψ)) v(ψ) * v(ψ + δψ) v(ψ) + i v(ψ) * (v(ψ + δψ) -v(ψ)) v(ψ) * v(ψ + δψ) ρv(ψ) We have v(ψ) * v(ψ + δψ) = v(ψ) * v (ψ)δψ + o(δψ). Differentiating v(ψ) * v(ψ) = 1 yields v(ψ) * v (ψ) = 0, i.e. v(ψ) and v (ψ) are orthogonal. Hence v(ψ) * v (ψ) = |v (ψ)| = 0
and we have

z(ψ, ψ + δψ) = h(ψ) + o(δψ) Thus, lim δψ→0,δψ =0 D + ψ ∩ D + ψ+δψ = h(ψ), which concludes the proof. Example 3.3. If f represents a circle f (ψ) = Ae iψ and if ρ < A cos ∆ ψ 2 , then formula (3.20) yields h(ψ) = (cos ∆ 2 -ρ)e i(ψ+ ∆ 2 )
Thus, A ρ,∆ is the open disk of center O and radius cos ∆ 2 -ρ.

Practical determination of an origin

We focus on the case of an ellipse. The best choice for an origin would be the ellipse center C. A practical difficulty is that the parameters of the ellipse can not be directly computed by the sampled noisy measurements. One simple alternative is to chose the mean of the measurements to average the noise contribution to zero. Though this solution is easy to implement, it might reveal troublesome as the average of the samples can be very different from C (thus not in A ρ,∆ ). The culprit can be that the f (ψ[k]) values may be not homogeneously distributed along the ellipse. One way to circumvent this problem is to compute an origin that only depends on the geometric figure defined by the samples, disregarding potential aggregates of values.

. , y[k],

say z 1 , . . . , z n k so that the polygon line Γk of ordered vertices z 1 , . . . , z n k is convex and so that its shape reflects the one of Γ k . For example, one could chose Γk as the convex hull of Γ k . An other way could be to compute the largest convex polygon included in the set y[1], . . . , y [k]. Still another way could be to consider not the Chebyshev center but the gravity center of the polygonal Γk (centroid). We illustrate the first solution and we call Chebyshev center of Γ k the Chebyshev center of its convex hull.

On Fig. 3.7, we simulated a non-homogeneous

ψ[k] distribution around a circle, with ρ = 1. ∆ is computed as the largest |ψ[k + 1] -ψ[k]|.

Summary and main result

We now summarize the estimation results presented in this chapter.

Theorem 3.1: Estimation of the phase of a single axis rotation from samples. Solution to Problem 3.1

Consider measurements of the form

y[k] = f (ψ[k]) + n[k] ∈ C, 1 ≤ k ≤ N where • f is a 2π-periodic function valued in C parametrizing a Jordan curve C • ψ[k] = ψ(k∆t) • n[k] is a measurement noise Assume that f is such that its Fourier expansion {c n } n∈Z satisfies |c 1 | > n =0,1 |c n | + |c 0 -z 0 |
Assume that the interior region I defined by the boundary C is strictly convex. Assume that the noise n is (uniformly) bounded by ρ and that |ψ

[k + 1] -ψ[k]| is (uniformly) bounded by ∆ < π. Then consider A ρ,∆ = I ∩ 2π ψ=0 {z ∈ C, v * (ψ)(z -f (ψ)) > ρ} where v(ψ) = e +iπ/2 f (ψ+∆)-f (ψ) |f (ψ+∆)-f (ψ)| . If A ρ,∆ is not empty, then the following sequence ψz 0 [k] = k-1 j=1 arg -π y[j + 1] -z 0 y[j] -z 0 (3.21)
where z 0 ∈ A ρ,∆ provides an estimate of ψ[k] with an error that is bounded by

|e z 0 | ∞ ≤ 2 arcsin n =0,1 |c n | + |c 0 -z 0 | |c 1 | + 2 arcsin ρ δ(z 0 ) where δ(z 0 ) = min ζ |f (ζ) -z 0 |.
In practice, a recommendation it to select z 0 as one of the following: i) the Chebyshev center of measurements, ii) the polygon centroid.

Corollary (Particular case of an ellipse). Further assume that f defines an ellipse

f (ψ) = C + A cos ψ + i √ 1 -e 2 A sin ψ. Then A ρ,∆ is not empty if and only if ρ < A √ 1 -e 2 cos ∆
|e C | ∞ ≤ 2 arcsin 1 - √ 1 -e 2 1 + √ 1 -e 2 + 2 arcsin ρ A √ 1 -e 2
In this chapter we have studied a simple estimation procedure to reconstruct the phase of a 2-dimensional vector data representing periodic dynamics. The sensor discrepancies such as biases and non-linearities were represented by a complex valued function f . The underlying idea stems from the classic notion of index of a curve and relates to the mathematics of complex and curve analysis. The application of the method to real data naturally raises the problem of noise and sampling. These issues have been treated by careful analysis, and sufficient conditions have been proposed to guarantee that the estimation is not jeopardized. Explicit bounds have been found for the error estimate. A critical notion of the proposed approach is the definition of suitable or allowed origins. Necessary conditions on noise level and on sampling rate have been derived to guarantee their existence. These conditions can be made explicit for the particular case of an ellipse, which is a natural approximation for numerous real situations.

Sensor non-linearities are illustrated through experiments for coarse Sun sensors in Appendix D.1. To reduce these non-linearities, off-line calibration is necessary. These non-linearities are easy to apprehend for single-axis rotation, where only one angle is involved. Their impact for full 3D rotation is much more complicated. In the rest of the

ν s = 100 [Hz] ν s = 50 [Hz] ν s = 10 [Hz] SNR = 30 [dB] σ = 5.7 • σ = 6.3 • σ = 6.5 • SNR = 13 [dB] σ = 14.2 • σ = 13.5 • σ = 14.4 • SNR = 5 [dB] σ = 24.5 • σ = 23.8 • σ = 22.9 •
thesis, we assume that the sensor have been calibrated, so that only additive measurement noise remains. To illustrate these results, we now perform simulations for the rotation of a satellite during a rest-to-rest reorientation maneuver. For t ∈ [0, 3], ψ(t) = 1, for t ∈ [3, 6], ψ(t) = -1. We consider several sampling rates ranging from 10 [Hz] to 100 [Hz] and SNR ranging from 30 [dB] (noise over signal amplitude 22%) to 5 [dB] (noise over signal amplitude 78%). Typical results are pictured in Fig. 3.8. Numerical results are presented in Table 3.1. As expected the estimate error grows with the noise amplitude. The sampling rate however does not seem to affect the error standard deviation.

Simulation illustrations

It appears that the angle ψ is well estimated and that the spin motion of the satellite is relatively accurately estimated, enabling monitoring strategies. Post filtering of the estimate can also be used on-line or off-line (to avoid filtering lag) to smooth out visible artifacts. To obtain the presented results, the Chebyshev center of the samples was computed only once (using linear programming techniques with an interior point method).

Alternatively, the centroid point could be used, which reveals significantly lighter in terms of computational effort. The estimation algorithm itself consists in simple arithmetic operations. No pre-filtering of data was employed.

Chapter 4 Rotation with small variations of the main axis

Rotation à axe légèrement variable. Dans ce chapitre nous considérons le cas d'un axe de rotation légèrement perturbé. Nous employons une paramétrisation de l'attitude par les angles d'Euler. Deux sous-cas du mouvement de l'axe sont envisagés : i) pour un mouvement de dérive lente, nous employons une extension de l'estimateur à axe fixe qui repose sur les études séparées de la phase et du module d'un nombre complexe bien choisi ii) pour une oscillation à faible amplitude, nous proposons un estimateur de phases qui repose sur une séparation naturelle de deux fréquences caractéristiques du mouvement contenues dans le spectrogramme des mesures. We now consider the more complex problem of small displacement of the rotation axis. This is still not a general case since it does not account for every possible rotation. Yet, the problem is of practical importance and can not be fully reduced to a simple phase estimation as in Chapter 3. It requires a more complex analysis, which we perform in this chapter.

In Section 4.1, we introduce a parametrization of the rotation matrix with Euler angles and formulate the problem. We consider two kinds of such motions, requiring specific analysis tools. In Section 4.2, we consider the motion of a high-spin rigid body whose rotation axis is slowly drifting. A typical example is a ballistic launch. An extension of estimator (3.3) gives a solution to the formulated problem. In Section 4.3, we consider the motion of a rigid body whose rotation axis is tilting with a small amplitude. This is typically the case of a long-axis mode celestial body such as the Toutatis 4179 asteroid [START_REF] Mueller | The diagnosis of complex rotation in the lightcurve of 4179 Toutatis and potential applications to other asteroids and bare cometary nuclei[END_REF].

Problem statement

We use the 'ZXZ' Euler angles parametrization of the rotation matrix R [LL82], as defined in Appendix A.2. For convenience, their pictorial representation is recalled below.

The angular velocity and the Euler angles are related by (see e.g. [START_REF] Landau | Mechanics. MIR Moscou[END_REF] (35,1)) As in Chapter 3, we use the complex valued variable y a 1 -ia 2 which is a projection of the measurement a onto a horizontal plane. Using the Euler angles expression of the rotation matrix yields

ω 1 = φ sin ψ sin θ + θ cos ψ ω 2 = φ cos ψ sin θ -θ sin ψ ω 3 = φ cos θ + ψ      (4.1)
y = a 1 (cos ϕ cos ψ -sin ϕ sin ψ cos θ + i(cos ϕ sin ψ + sin ϕ cos ψ cos θ)) + a 2 (sin ϕ cos ψ + cos ϕ sin ψ cos θ -i(-sin ϕ sin ψ + cos ϕ cos ψ cos θ)) + a 3 (sin ψ sin θ -i cos ψ sin θ) = a 1 (cos(ϕ + ψ) + i sin(ϕ + ψ) + (1 -cos θ)(sin ϕ sin ψ -i sin ϕ cos ψ)) + a 2 (sin(ϕ + ψ) -i cos(ϕ + ψ) -(1 -cos θ)(cos ϕ sin ψ -i cos ϕ cos ψ)) -ia 3 sin θ (cos ψ + i sin ψ)
= a 1 e i(ϕ+ψ) -ia 1 (1 -cos θ) sin ϕe iψ -ia 2 e i(ϕ+ψ) + ia 2 (1 -cos θ) cos ϕe iψ -ia 3 sin θe iψ = (a 1 + ia 2 )e i(ϕ+ψ) + 1 -cos θ 2 e iψ (a 1 e -iϕa 1 e iϕ + ia 2 e -iϕ + ia 2 e iϕ ) -ia 3 sin θe iψ

In the end, the measurement can be written

y = a 1 -ia 2 2 (1 + cos θ)e i(ϕ+ψ) -ia 3 sin θe iψ + a 1 + ia 2 2 (1 -cos θ)e i(ψ-ϕ) (4.2)
Note that for θ = ϕ = 0 this formula reduces to

y = (a 1 -ia 2 )e iψ
which is exactly (3.1). In this Chapter only we assume that the inertial coordinates of the reference vector are known.

Assumption 4.1. a 1 , a 2 , a 3 are known.

Mathematically, we can now formulate an estimation problem as follows.

Problem 4.1

From measurements of the form (4.2) where ϕ, ψ, θ are the respective precession, spin and nutation angle of a rigid body, and under Assumption 4.1, find estimates φ, ψ, θ, ω of the Euler angles ϕ, θ, ψ and of the angular rate ω. Remark 4.1. One can make several choices to define the inertial frame of reference, depending on the analysis one wishes to perform. For example, one can exhibit the nonobservability of the rotation angle around the reference vector a as follows. If the inertial frame is chosen so that the a = e 3 , then a 1 = a 2 = 0, a 3 = 1. The measurement equation (4.2) would reduce to y = -i sin θe iψ . Under this form, it is clear that angles ψ and θ are observable and that ϕ is not, as it does not impact the measurement. Thus, Problem 4.1 does not have a satisfying solution in the general case. Typically, troubles will arise if the rigid body is animated with precession around the reference vector.

Remark 4.2 (on Assumption 4.1). As explained in Remark 2.1, ω does not depend on the choice of the inertial frame of reference. Thus, we do not lose generality by considering that the inertial frame and body frame are aligned at initial time t 0 , i.e

e 1 = b 1 (t 0 ), e 2 = b 2 (t 0 ), e 3 = b 3 (t 0 )
It follows that, up to measurement noise,

a = a(t 0 )
As a result, Assumption 4.1 is not restrictive for the problem of estimating ω. Moreover, we can consider that the initial values of the Euler angles are known

ϕ(t 0 ) = 0, ψ(t 0 ) = 0, θ(t 0 ) = 0

Example: a high-spin projectile in ballistic mode

In this section, we consider the example of a ballistic launch where the rotation axis is slowly moving in a plane while the rigid body has a high-spin angular rate. A typical example is a high-spin projectile traveling in ballistic mode after a shot. Typical measurements a(t) on the unit sphere and their projection in a horizontal plane are represented in Figure 4.2.

Estimation principle

For the motion under consideration, we formulate the following assumption Assumption 4.2. The precession ϕ is negligible compared to the spin ψ Under this assumption, we rewrite an approximation version of (4.2) with ϕ = 0, which yields

y (a 1 -ia 2 cos θ + ia 3 sin θ) e iψ = a 2 1 + (a 3 sin θ -a 2 cos θ) 2 e i(ψ+ζ(θ)) (4.3) where ζ(θ) arg -π (a 1 -ia 2 cos θ + ia 3 sin θ) is a known function of θ. Equation (4.1) is approximated into ω 1 θ cos ψ ω 2 -θ sin ψ ω 3 ψ      (4.4)
The quantity ψ + ζ(θ) is given from y by the estimator (3.3). The modulus of y is

|y| = a 2 1 + (a 3 sin θ -a 2 cos θ) 2
which allows one to find θ and, in turn, ψ. An estimate of ω directly follows from equation (4.4). A schematic view of this simple algorithm is represented in

Simulation results

This simple algorithm was tested on data for a projectile produced with the simulator FQW98 developed by V. Fleck at French-German Research Institute of Saint-Louis1 . This high fidelity simulator accounts for [START_REF] Fleck | Introduction à la Balitistique Extérieure[END_REF] gravity, wind, aerodynamic forces (drag, lift, Magnus effect) and torques, and the Coriolis effect. The parameters of this simulation are as described in [START_REF] Changey | Modélisation et estimation par filtrage non linéaire, de l'attitude d'un projectile à partir de magnétomètres[END_REF]. The inertial frame is radial to the Earth's surface and defined by the launch site latitude and launch azimuth. the ones represented in Figure 4.2. The angular velocity values and estimation errors are represented in Figure 4.6 and 4.7. The plots for ω 2 are not represented, as they are strongly similar to the ones for ω 1 . The relative error amplitude for ω 1 starts around 100%, remains quite large for the first 5 seconds and shrinks down to about 10%. The culprit is that ψ shows fast dynamics during this time line, so that the approximation ψ = 0 in (4.4) is rough. The relative error for ω 3 is much lower, around 2.10 -4 .

Tilting of the rotation axis

We now consider a second kind of motion where the rotation axis is tilting. This is typically the case of a long-axis mode (LAM) celestial bodies such as 4179 Toutatis [START_REF] Mueller | The diagnosis of complex rotation in the lightcurve of 4179 Toutatis and potential applications to other asteroids and bare cometary nuclei[END_REF]. It can also be observed on smaller objects such as space debris [START_REF] Bonnal | Active debris removal: Recent progress and current trends[END_REF]. Such motions are characterized by a small nutation. Assumption 4.2 is replaced with the following one.

Assumption 4.3. The nutation θ is small and satisfies

a 2 1 + a 2 2 θ |a 3 |
In the following, we invoke arguments of frequency separation in signal processing to explain how the sought-after information on the angles can be obtained from the measurements, under careful assumptions.

Frequency analysis of the measurement equation

We rewrite the measurement equation (4.2) as

y(t) = y 1 (t)e iφ 1 (t) + y 2 (t)e iφ 2 (t) + y 3 (t)e iφ 3 (t) (4.5) with φ 1 ϕ + ψ y 1 (a 1 -ia 2 ) 1+cos θ 2 φ 2 ψ y 2 ia 3 sin θ φ 3 ψ -ϕ y 3 (a 1 + ia 2 ) 1-cos θ 2
Under this form, the measurement is the sum of three terms with phase φ i and modulated amplitude y i . Interestingly, if the φ i , y i and their derivatives satisfy some desirable properties, one can recover the instantaneous frequencies dφ i dt by mean of windowed Fourier transforms (see [START_REF] Mallat | A wavelet tour of signal processing[END_REF] Chapter 4.4) as we will explain below. This is a solution for this source separation problem.

According to Assumption (4.3) we have

y 1 = (a 1 -ia 2 ) 1 + O(θ 2 )) y 2 = ia 3 θ 1 + O(θ 2 ) y 3 = (a 1 + ia 2 )θ 2 1 4 + O(θ 2 )
which implies |y 3 | |y 1 |, |y 2 |. Thus, we consider an approximation of the measurements (4.5) of the form y = y 1 (t)e iφ 1 (t) + y 2 (t)e iφ 2 (t) (4.6)

Though we apparently lose generality, one should bear in mind that the following analysis could, with additional tedious but relatively easy calculations, extend to the general case (4.5) (practical implementation remains simple).

Let g be a (window) function even, positive, with finite support in [-1 2 , 1 2 ] and unit L 2 norm

1/2 -1/2 g 2 (t)dt = 1
Its Fourier transform G defined by

G(ν) = ∞ -∞ g(t)e -iνt dt satisfies G(0) ≥ G(ν), ∀ν
and decays rapidly to 0 when ν grows. This decay is characterized by (among other) the bandwidth ∆ν defined as

|G( ∆ν 2 )| = |G(0)| √ 2
In general, we have

G(ν) G(0), ∀|ν| ≥ 2∆ν
Take a scaling factor s > 0 and define the windowed Fourier transform Sy (a.k.a. spectrogram [START_REF] Mallat | A wavelet tour of signal processing[END_REF]) of the signal y at time u and frequency ν using g as window scaled by s as

Sy(u, ν) = 1 s s/2 -s/2 y(u + t)g( t s )e -i(u+t)ν dt = 1 s s/2 -s/2 y 1 (u + t)g( t s )e i[φ 1 (u+t)-(u+t)ν] dt + 1 s s/2 -s/2 y 2 (u + t)g( t s )e i[φ 2 (u+t)-(u+t)ν] dt S 1 (u, ν) + S 2 (u, ν)
The spectrogram satisfies the following property (adapted from [START_REF] Mallat | A wavelet tour of signal processing[END_REF], Theorem 4.5).

Proposition 4.1. For j = 1, 2 we have

S j (u, ν) = y j (u)e i(φ j (u)-νu) G s(ν -φj (u)) + ε j (u, ν) with |ε j (u, ν)| ≤ s 2 √ 3 | ẏj (u)| + s 2 8 √ 5 |ÿ j | u + |y j (u)| s 2 8 √ 5 | φj | u (4.7)
where

| • | u designates the supremum over [u -s/2, u + s/2]
Before proving this result, let us sketch an interpretation. Assuming the error terms ε j are small, the spectrogram reduces to two main terms Sy(u, ν) y 1 (u)e i(φ 1 (u)-νu) G s(ν -φ1 (u)) + y 2 (u)e i(φ 2 (u)-νu) G s(ν -φ2 (u)) Thus, if the instantaneous frequencies φ1 and φ2 are separated enough with respect to the bandwidth of G, namely if they satisfy

s| φ1 -φ2 | ≥ 2∆ν (4.8)
then |S(u, •)| has two main lobes for

ν j φj (u), j = 1, 2
More precisely, we have

Sy(u, φj (u)) y j (u)e i[φ j (u)-u φj (u)] G(0)
Hence, one can recover the instantaneous frequencies φj (u) by detecting the peaks of the spectrogram [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. Further, the corresponding amplitude of the lobe gives us y j (u). As will appear, these informations are instrumental in estimating the three Euler angles.

Proof. (of Proposition 4.1) For a fixed value (u, ν) and t ∈ [us 2 , u + s 2 ], denote

y j (u + t) = y j (u) + t ẏj (u) + t 2 2 α j (t) φ j (u + t) = φ j (u) + t φj (u) + t 2 2 β j (t) with |α j (t)| ≤ |ÿ j | u , |β j (t)| ≤ | φj | u
For brevity, we omit the bounds ± s 2 of the integrals below. For j = 1, 2 we have

S j (u, ν) = y j (u) e i[φ j (u)-νu] s g( t s
)e -it(ν-φj (u)) dt

+ y j (u) e i[φ j (u)-νu] s g( t s )e -it(ν-φj (u)) (e i t 2 2 β j (t) -1)dt + ẏj (u) 1 s tg( t s )e i[φ j (u+t)-ν(u+t)] dt + 1 2s t 2 α j (t)g( t s )e i[φ j (u+t)-ν(u+t)] dt = y j (u)e i[φ j (u)-νu] G s(ν -φj (u)) + ε j (u, ν)
where

|ε j (u, ν)| ≤ |y j (u)| g( t s ) e i t 2 2 β j (t) -1 dt s + | ẏj (u)| |t|g( t s ) dt s + |ÿ j | u 2 t 2 g( t s ) dt s
To conclude, we use the following two technical points.

1. For any x ∈ R, we have

|e ix -1| 2 = (cos x -1) 2 + sin 2 x = 2(1 -cos x) = 2 n≥1 (-1) n-1 x 2n (2n)! ≤ x 2
Thus, for all t

|e i t 2 2 β j (t) -1| ≤ t 2 2 |β j (t)| ≤ t 2 2 | φj | u 2.
For any n ∈ N, we have

s/2 -s/2 |t n |g( t s ) dt s = s n 1/2 -1/2 |t| n g(t)dt ≤ s n 1/2 -1/2 t 2n dt 1/2 -1/2 g(t) 2 dt = s n 2 n √ 2n + 1
Using the latter point for n = 1, 2 concludes the proof.

If the separating condition (4.8) is met, the decay of G suffices to isolate the contributions of S 1 (u, φ1 (u)) and S 2 (u, φ2 (u)) in the frequency domain. To safely consider that the spectrogram peaks detection will give satisfactory result, one should verify that none of the 6 terms of (4.7) for j = 1, 2 disturbs the location of any of the peaks ν = φ1 (u), φ2 (u). Considering (4.7) and the fact that the peak at frequency ν = φj (u) has amplitude |y j (u)|G(0) allows to derive the following set of conditions that guarantees, together with separating condition (4.8), reliability of the peaks detection

s| ẏj | 2C 1 |y j | (4.9) s 2 |ÿ j | u C 2 |y j (u)|, ∀u (4.10) s 2 | φj | C 2 (4.11)
where we have denoted • compute estimates

C 1 2 √ 3G(0) , C 2 8 √ 5G (
dϕ dt [k] = ν 1 -ν 2 , dψ dt [k] = ν 2 • estimate ϕ[k], ψ[k]
by cumulative numerical integration of the estimates of their respective derivative defined above

• estimate θ[k] from cos θ[k] = 2A 1 a 2 1 + a 2 2 G(0) -1, sin θ[k] = A 2 |a 3 |G(0) • estimate ω[k] from (4.1) Remark 4.3.
• The algorithm provides average values (over a time window of size s) of φ, ψ and θ.

Thus, it does not allow to observe high-frequency variations (compared to 1 s ) of these quantities.

• The algorithm needs initial values of ϕ and ψ to estimate the full rotation. This is not restrictive, as explained in Remark 4.2.

spectrogram + peaks detection

a 1 , a 2 , a 3 y ω 3 = φ cos θ + ψ φ, ψ ϕ 0 , ψ 0 ω 1 = φ sin ψ sin θ + θ cos ψ ω 2 = φ cos ψ sin θ -θ sin ψ φ, ψ θ ω3 d dt ω1 , ω2

• The computation of Sy(t[k], •) requires the values of y over

[t[k] -s/2, t[k] + s/2].
Therefore, it introduces a short lag.

• The estimation of θ[k] can be implemented only if |a 3 | is known. It is the only place of this section where we invoke Assumption 4.1. This is not restrictive, as explained in Remark 4.2.

In the next section, we determine sufficient conditions such that (4.8)-(4.11) hold in the case of a free-rotation.

Frequency separation for long-axis mode

We consider that the rigid body satisfies the free Euler's equation (2.13) (with J 1 > J 2 > J 3 ) and that the long-axis mode 2 criterion holds

d 3 |ω 1 (0)| < d 1 |ω 3 (0)|
Without loss of generality, we assume that the body frame is chosen so that ω 1 (0) ≥ 0, ω 3 (0) > 0. We establish sufficient conditions bearing on the inertia parameters and the initial condition guaranteeing practicability of the previously described technique of estimation.

As explained in details in Appendix C, there exists t 0 ∈ R such that for all t

ω 1 (t) = A 1 cn(w(t -t 0 )), ω 2 (t) = A 2 sn(w(t -t 0 )), ω 3 (t) = A 3 dn(w(t -t 0 ))
2 See Appendix C for an explanation on "long-axis mode" and "short-axis mode"

where A 1 , A 2 , A 3 , w are defined as follows

m 1 d 1 ω 1 (0) 2 + 1 |d 2 | ω 2 (0) 2 1 d 3 ω 3 (0) 2 + 1 |d 2 | ω 2 (0) 2 ∈ [0, 1) A 1 ω 1 (0) 2 + d 1 |d 2 | ω 2 (0) 2 ≥ 0 A 2 |d 2 | d 1 ω 1 (0) 2 + ω 2 (0) 2 ≥ 0 A 3 ω 3 (0) 2 + d 3 |d 2 | ω 2 (0) 2 > 0 w -d 1 |d 2 |ω 3 (0) 2 + d 1 d 3 ω 2 (0) 2 < 0
and cn, sn, dn are the Jacobi elliptic functions 3 of parameter m ∈ [0, 1). The angular momentum of the rigid body M is constant. We denote by M its norm

M |M|

Conveniently, we choose the inertial frame so that M is aligned with the vertical vector e 3 , namely e 3 = M M Then, we have Jω = RM = M Re 3 , which yields in term of Euler angles

J 1 ω 1 = M sin ψ sin θ (4.12)
J 2 ω 2 = M cos ψ sin θ (4.13)

J 3 ω 3 = M cos θ (4.14)
Nutation angle: meeting conditions (4.9)-(4.10)

For brevity, we now assume (without loss of generality) that t 0 = 0. In our problem, the nutation angle is characterized by its cosine value. From (4.14) we have

cos θ(t) = J 3 M ω 3 (t) = J 3 A 3 M dn(wt)
Remark 4.4. dn(•) oscillates between 1 and √ 1 -m > 0. Thus, θ(t), which is the angle between the angular momentum M and the third body frame vector b 3 , oscillates between a minimum value θ 0 arccos

J 3 A 3 M (4.15)
and a maximum value

θ max = arccos J 3 A 3 M √ 1 -m < π 2
As J 3 is the smallest principal moment of inertia, b 3 corresponds to the longest dimension of the rigid body. This explains why this kind of rotation-motion is called long-axis mode.

We now differentiate cos θ(t) = cos θ 0 dn(wt) with respect to time, using the differential system satisfied by the Jacobi elliptic functions recalled in Proposition B.1.

d dt cos θ(t) = -mw cos θ 0 sn(wt)cn(wt) d 2 dt 2 cos θ(t) = mw 2 cos θ 0 dn(wt) sn 2 (wt) -cn 2 (wt)
Hence, we obtain the following inequalities

cos θ ≥ cos θ 0 √ 1 -m, |(cos θ) | ≤ cos θ 0 m|w| 2 , |(cos θ) | ≤ cos θ 0 mw 2
Thus, conditions (4.9)-(4.10) are satisfied for

y 1 = a 1 -ia 2 2 (1 + cos θ) if swm 2 √ 1 -mC 1 s 2 w 2 m √ 1 -mC 2 (4.16)
We now differentiate sin θ. We have

sin θ(t) = √ 1 -cos 2 θ = 1 -cos 2 θ 0 dn 2 (wt) = 1 -cos 2 θ 0 (1 -msn 2 (wt)) = sin θ 0 1 + m cos 2 θ 0 sin θ 0 sn 2 (wt)
Denote η m cos 2 θ 0 sin θ 0 (4.17)

The differentiation yields (after calculations)

sin θ ≥ sin θ 0 , |(sin θ) | ≤ ηw sin θ 0 2 , |(sin θ) | ≤ (1 + η)ηw 2 sin θ 0
Thus, conditions (4.9)-(4.10) are satisfied for

y 2 = ia 3 sin θ if swη 4C 1 s 2 w 2 (1 + η)η C 2 (4.18)
Precession rate and spin rate: meeting conditions (4.8) and (4.11)

The spin angle is recovered via the value of its tangent from (4.12),(4.13)

tan ψ(t) = J 1 J 2 ω 1 (t) ω 2 (t) = J 1 A 1 J 2 A 2 cn (wt) sn (wt) = J 1 √ d 1 J 2 |d 2 | cn (wt) sn (wt)
which allows to compute ψ as

(1 + tan 2 ψ(t)) ψ(t) = d dt    J 1 √ d 1 J 2 |d 2 | cn(wt) sn(wt)   
Equations (4.1) give the derivatives of the precession angle

φ(t) = ω 1 (t) sin ψ(t) + ω 2 (t) cos ψ(t) sin θ(t) = J 1 ω 2 1 (t) + J 2 ω 2 2 (t) J 2 1 ω 2 1 (t) + J 2 2 ω 2 2 (t)
M ϕ and ψ are strictly increasing and, after some heavy calculations, one can show that

φ ≥ M J 1 ψ ≤ |d 2 | 3 2 d 3 d 1 J 2 J 1 M 2 J 2 1 φ + ψ ≤ 1 2 |d 2 | 3 2 d 3 d 1 J 2 J 1 M 2 J 2 1
Thus, condition (4.11) is guaranteed for φ 1 = ϕ + ψ and φ 2 = ψ if

|d 2 | 3 2 d 3 d 1 J 2 J 1 M 2 J 2 1 (4.19)
The frequency separation condition (4.8) holds for

φ 1 -φ 2 = ϕ if s M J 1 ≥ 2∆ν (4.20)

Simplification for the case of nearly axis-symmetric rigid body

This case, which is relevant in many application, corresponds to

J 2 = J 1 (1 -) > J 3 , 0 < 1
In this case, the third ratio of inertia can be written

d 3 = J 1 -J 2 J 3 = J 3 1
After some first order calculations, one can regroup and simplify conditions (4.16) to (4.20) into the following result.

Proposition 4.2. If 1 -J 2 J 1 satisfies 1 then the following conditions guarantee that conditions (4.8) to (4.11) hold

s M J 1 ≥ 2∆ν, s M J 1 2C 1 , s M J 1 2 1 - J 3 J 1 J 3 J 1 C 2
Remark 4.5. These assumptions rule out rotation motions with strong precession around a mentioned in Remark 4.1.

Simulation results

In this section we illustrate the sensitivity of the algorithm with respect to the time-window size s. We simulate a rigid body as a homogeneous ellipsoid with semi-axes 0.

a = √ 3 (1 1 1) , ω(0) = (101 0 630) [ • /s], ϕ(0) = 0, ψ(0) = π 2 , θ(0) = π 8
For the windowed Fourier transforms we use a Hann window [START_REF] Mallat | A wavelet tour of signal processing[END_REF]

g(t) = 2 2 3 cos 2 πt
We simulate the Euler equations without torque and with a damping torque.

No torque. The free Euler equations (2.13) and measurements equation (2.17) are solved using a Runge Kutta 4 numerical scheme with fixed time step 0.01 [s]. The corresponding Euler angles are reported in Figure 4.9. θ(t) is periodic. ψ and ϕ grow in an almost time-affine trend, up to a bounded term. This last term (hardly visible in Figure 4.9) is periodic for ψ and aperiodic for ϕ. Therefore, the rigid body never goes back to its initial attitude (see [START_REF] Landau | Mechanics. MIR Moscou[END_REF] page 182).

The corresponding measurements are represented in Figure 4.10 evolving on the unit sphere and projected onto a horizontal plane. The angular rate history is represented in Damping torque. We now add a damping torque so that the Euler's equation can be written ω = E(ω) -0.02ω

The angular rate history is represented in Figure 4.13 together with the corresponding

Perspectives

In this chapter, we have exposed two cases of 3D rotation of a rigid body in free motion which can be estimated from vector measurements. Some assumptions guaranteeing the non-ambiguity of the measurements have been established by a careful study of the solutions of the free rotation dynamics and have been formulated in terms of inertia parameters, initial angular rate, and maximal nutation.

The problem under study in this contribution is relatively general. Specific studies could be conducted in a similar manner for other governing dynamics: the aerodynamics of atmospheric flight could be considered explicitly. Interestingly, one could draw some

t [s] ν [rad/s] ν [rad/s]

Part II Angular rate non-linear observers Chapter 5 A non-linear observer using two vector measurements

Observateur non-linéaire à partir de mesures de deux directions. Dans ce chapitre nous utilisons une approche d'estimation directe et présentons un observateur non-linéaire de ω ne nécessitant pas la connaissance des coordonnées inertielles des vecteurs de référence a, b. C'est la contribution principale de cette thèse. Nous faisons l'hypothèse usuelle que la matrice d'inertie J et le couple τ sont connus. Nous construisons l'observateur à partir de deux vecteurs de référence non-colinéaires. Nous établissons une preuve de convergence et une estimation du bassin d'attraction.

In this chapter, we consider that the matrix of inertia J and the torque τ are known and we derive angular rate non-linear observer from two reference vector measurements. Extensions of the observer presented in Chapters 7, 8 will permit to relax these assumptions and to estimate the torque and the ratios of inertia.

The chapter is organized as follows. In Section 5.1, we formulate the problem statement. In Section 5.2, we define a non-linear observer with extended state and output injection. To prove its convergence, the error equation is identified as a linear time-varying (LTV) system disturbed by a linear-quadratic term. The dominant part of the LTV dynamics can be shown, by a scaling resulting from a high gain design, to generate an arbitrarily fast exponentially convergent dynamics. In turn, this property reveals instrumental to conclude on the exponential uniform convergence of the error dynamics. Illustrative simulation results are given in Section 5.3.

Problem statement

In this chapter we consider the following assumptions. is constant for all times. Without loss of generality, we assume a b ≥ 0 (if not, one can simply consider -a instead of a). This Assumption also implies that the measurements a, b satisfy the differential system (2.17), namely ȧ = a × ω, ḃ = b × ω

The problem we address in this paper is the following.

Problem 5.1

Under Assumptions 5.1-5.2-5.3, find an estimate ω of ω from the measurements a, b.

Observer definition and analysis of convergence

Observer dynamics

To solve Problem 5.1, the main idea developed here is to consider the reconstruction of the augmented 9-dimensional state S by its estimate Ŝ

S =    a b ω    , Ŝ =    â b ω   
According to (2.8)-(2.17), the state is governed by

Ṡ =    a × ω b × ω E(ω) + χ    (5.1)
and the following observer is proposed

Ṡ =    a × ω + αk(a -â) b × ω + αk(b -b) E(ω) + χ + k 2 a × â + k 2 b × b    (5.2) where α ∈ (0, 2 √ 1 -p) and k > 0 are constant (tuning) parameters. Denote S S -Ŝ =    a - â b - b ω - ω   
(5.3) the error state. We have

Ṡ =    -αkI 0 [a × ] 0 -αkI [b × ] k 2 [a × ] k 2 [b × ] 0    S +    0 0 E(ω) -E(ω)   
(5.4) In Section 5.2.4 we will exhibit, for each value α ∈ (0, 2 √ 1 -p), a threshold value k * such that for k > k * , S converges locally uniformly exponentially to zero.

Preliminary change of variables and properties

The subsequent study of the dynamics (5.4) employs a preliminary change of coordinates. Denote

X    a - â b - b ω-ω k   
(5.5) yielding a reformulation of (5.4)

Ẋ = kA(t)X +    0 0 E(ω)-E(ω) k   
(5.6)

with

A(t)      -αI 0 a(t) × 0 -αI b(t) × a(t) × b(t) × 0     
(5.7) which we will analyze as an ideal linear time-varying (LTV) system Ẋ = kA(t)X (5.8) disturbed by the input term

ξ    0 0 E(ω)-E(ω) k   
(5.9)

The main idea for the analysis is that for sufficiently large values of k, the rate of convergence of (5.8) will ensure stability of system (5.6). As a preliminary step, we start by upper-bounding A(t) and the disturbance (5.9).

Proposition 5.1 (Bound on the unforced LTV system). A(t) defined in (5.7) is upperbounded by

A max max √ 2 + 2α 2 , √ 3 + α 2 Proof. Let Y ∈ R 9 such that |Y | = 1. For convenience, denote Y =    Y 1 Y 2 Y 3    with Y i ∈ R 3 .
One has

|A(t)Y | 2 =| -αY 1 + a × Y 3 | 2 + | -αY 2 + b × Y 3 | 2 + |a × Y 1 + b × Y 2 | 2 ≤(1 + α 2 ) |Y 1 | 2 + |a × Y 3 | 2 + |Y 2 | 2 + |b × Y 3 | 2 + 2 |a × Y 1 | 2 + |b × Y 2 | 2 ≤ max 2 + 2α 2 , 3 + α 2 |Y | 2 = A 2 max Hence, A(t) = max |Y |=1 |A(t)Y | ≤ A max .
Proposition 5.2 (Bound on the disturbance). For any X, ξ is bounded by

|ξ| ≤ d max √ 2ω max |X| + k|X| 2
Proof. We have

|ξ| = 1 k |E(ω) -E(ω)|
with, due to the quadratic nature of E(•),

E(ω) -E(ω) = J -1 (J ω × ω + Jω × ω -J ω × ω) =    d 1 (ω 2 ω3 + ω2 ω 3 ) d 2 (ω 3 ω1 + ω3 ω 1 ) d 3 (ω 1 ω2 + ω1 ω 2 )    -    d 1 ω2 ω3 d 2 ω3 ω1 d 3 ω1 ω2    δ 1 -δ 2
The following inequality is straightforward

|δ 2 | ≤ d 2 max |ω| 2 Moreover, by Cauchy-Schwarz inequality (ω 2 ω3 + ω2 ω 3 ) 2 ≤ (ω 2 2 + ω 2 3 )(ω 2 2 + ω2 3 ) ≤ (ω 2 2 + ω 2 3 )|ω| 2
Using similar inequalities for all the coordinates of δ 1 yields

|δ 1 | 2 ≤ 2d 2 max |ω| 2 |ω| 2 ≤ 2d 2 max ω 2 max |ω| 2 (5.10) Hence, |ξ| ≤ |δ 1 | + |δ 2 | k ≤ d max √ 2ω max ω k + k ω k 2 ≤ d max √ 2ω max |X| + k|X| 2
Remark 5.1 (on the quantity d max ). The quantity d max , which is less than 1 is not critical. Thus, for brevity we will only consider the conservative inequality

|ξ| ≤ √ 2ω max |X| + k|X| 2 (5.11)
d max is however essential to the analysis presented in Chapter 6. More details are given in Remark 6.2.

Analysis of the underlying LTV dynamics

We will now use a result on the exponential stability of LTV systems. The claim of [HI11] Theorem 2.1, which is instrumental in the proof of the next result, is as follows: consider a LTV system Ẋ = M (t)X such that

• M (•) is l-Lipschitz, with l > 0

• there exists K ≥ 1, c ≥ 0 such that for any t and any s ≥ 0, e M (t)s ≤ Ke -cs ("frozen time" asymptotic stability)

Then, for any t 0 , X 0 , the solution of Ẋ = M (t)X with initial condition X(t 0 ) = X 0 satisfies, for any t ≥ t 0 ,

|X(t)| ≤ Ke ( √ Kl ln K-c)(t-t 0 ) |X 0 |
Using this result, whose proof is given for convenience in Appendix B.3, we will show that the convergence of (5.8) can be tailored by choosing k to arbitrarily increase the rate of convergence, while keeping the overshoot constant.

Theorem 5.1: uniform exponential stability of system (5.8)

Let α ∈ (0, 2 √ 1 -p) be fixed. There exists a continuous function γ(k) satisfying

lim k→+∞ γ(k) = +∞
such that the solution of (5.8) satisfies

|X(t)| ≤ Ke -γ(k)(t-t 0 ) |X(t 0 )| with K 1 + α 2 √ 1-p 1 -α 2 √ 1-p
(5.12)

for any initial condition t 0 , X(t 0 ) and any t ≥ t 0 .

Proof. Consider any fixed value of t. We start by studying the frozen-time matrix A(t). Denote

µ 8(1 -p 2 ) > 0
Introduce the following (time-varying) matrices

P 1 =      a 0 b-pa √ 2(1-p 2 ) 0 b a-pb √ 2(1-p 2 ) 0 0 0      ∈ R 9×3 P 2 = 1 µ    2(pa -b) 0 2(a -pb) 0 α a × b - √ 8 -α 2 a × b    ∈ R 9×2 P 3 = 1 µ     2a × b 0 2a × b 0 α(b -a) 4(1 + p) -α 2 (a -b)     ∈ R 9×2 P 4 = 1 µ     2a × b 0 -2a × b 0 α(a + b) -4(1 -p) -α 2 (a + b)     ∈ R 9×2 and P = (P 1 |P 2 |P 3 |P 4 ) ∈ R 9×9
We have

P -1 A(t)P =      M 1 0 0 0 0 M 2 0 0 0 0 M 3 0 0 0 0 M 4      with M 1 = -αI, M j = 1 2   -α -α j -α 2 α j -α 2 -α   for j = 2, 3, 4 with α 2 2 √ 2 > α 3 2 √ 1 + p ≥ α 4 2 √ 1 -p > α
For all s ≥ 0 e A(t)s ≤ P P -1 e -α 2 s

Moreover P P -1 = λ max (P P ) λ min (P P )

where λ max (•), λ min (•) respectively designate the maximum and minimum eigenvalues. Besides,

P P =      I 0 0 0 0 Q 2 0 0 0 0 Q 3 0 0 0 0 Q 4     
with, for j = 2, 3, 4

Q j =   1 + α 2 α 2 j α α j 1 -α α j α α j 1 -α α j 1 -α 2 α 2 j   yielding the eigenvalues eig(P P ) = 1, 1 ± α 2 √ 2 , 1 ± α 2 √ 1 + p , 1 ± α 2 √ 1 -p ,
Thus, for all s ≥ 0 e A(t)s ≤ Ke -α

2 s with K = λ max (P P ) λ min (P P ) = 1 + α 2 √ 1-p 1 -α 2 √ 1-p
Let the tuning gain k > 0 be fixed. The scaled matrix kA(•) satisfies

e kA(t)s ≤ Ke -kα 2 s , ∀t, ∀s ≥ 0 Moreover, for any Y =    Y 1 Y 2 Y 3    ∈ R 9 and any t, s ∈ R, one has (kA(s) -kA(t))Y = k s t Ȧ(u)du Y = k    s t a(u) × ω(u)du Y 3 s t b(u) × ω(u)du Y 3 s t a(u) × ω(u)du Y 1 + s t b(u) × ω(u)du Y 2    Hence |(kA(s) -kA(t))Y | 2 ≤ 2ω 2 max k 2 |s -t| 2 |Y | 2 Thus, kA(•) is kL-Lipschitz with L √ 2ω max (5.13)
We now apply [HI11], Theorem 2.1 recalled earlier. For any t 0 and any X 0 , the solution of (5.8) with initial condition X(t 0 ) = X 0 satisfies for all t ≥ t 0

|X(t)| ≤ Ke ( √ kLK ln K-kα 2 )(t-t 0 ) |X 0 |
which concludes the proof with

γ(k) kα 2 - √ kLK ln K (5.14)
Remark 5.2. Additionally, one can note that

• γ(k) > 0 ⇔ k > 4LK ln K α 2
in which case Theorem 5.2.3 ensures exponential stability of system (5.8).

• γ(•) is strictly increasing for

k > 4K ln KL α 2 .

Convergence of the observer

Define r as

r(k) 1 √ A max K 3 1 - K 2 √ 2ω max γ(k) γ(k) k 3 2
(5.15) and k * as

k * = √ ln K + √ ln K + 2αK 2 α 2 √ 2Kω max > 0 (5.16)
The following holds Proposition 5.

r(k) > 0 if and only if k > k *

Proof. A simple rewriting of r(k) > 0 yields, successively,

r(k) > 0 ⇔ γ(k) > K 2 √ 2ω max = K 2 L ⇔ α 2 k - √ LK ln K √ k -K 2 L > 0 ⇔ √ k > √ LK ln K + √ LK ln K + 2αK 2 L α = √ k *
which concludes the proof.

We can now state the main result of the chapter.

Theorem 5.2: solution to Problem 5.1

For any α ∈ (0, 2 √ 1 -p), there exists k * defined by (5.16)-(5.12) such that for k > k * , the observer (5.2) defines an error dynamics (5.4) for which the equilibrium 0 is locally uniformly exponentially stable. The basin of attraction of this equilibrium contains the ellipsoid

S(t 0 ), |a(t 0 ) -â(t 0 )| 2 + |b(t 0 ) -b(t 0 )| 2 + |ω(t 0 ) -ω(t 0 )| 2 k 2 < r(k) 2
(5.17)

where r(k) is defined by (5.15).

Proof. Let k > k * . Consider the candidate Lyapunov function V 1 (t, X) kX +∞ t φ(s, t) φ(s, t)ds X (5.18)
where φ is the transition matrix of system (5.8).

Lemma 5.1. For all t, X we have

V 1 (t, X) ≥ 1 2Amax |X| 2 Proof. We have V 1 (t, X) = k +∞ t |φ(s, t)X| 2 ds
Consider s ≥ t. We have

d ds |φ(s, t)X| 2 = (φ(s, t)X) k(A(s) + A(s) )φ(s, t)X ≤ 2kA max |φ(s, t)X| 2 Thus d ds |φ(s, t)X| 2 ≥ -2kA max |φ(s, t)X| 2 from which we deduce |φ(s, t)X| 2 ≥ e -2kAmax(s-t) |φ(t, t)X| 2 = e -2kAmax(s-t) |X| 2
The result follows immediately.

Moreover, Theorem 5.2.3 implies that for all s ≥ t

|φ(s, t)X| ≤ Ke -γ(k)(s-t) |X| which gives V 1 (t, X) ≤ kK 2 +∞ t e -2γ(k)(s-t) ds|X| 2 = kK 2 2γ(k) |X| 2
and additionally

+∞ t φ(s, t) φ(s, t)ds ≤ K 2 2γ(k) We denote c 1 1 2A max , c 2 kK 2 2γ(k) (5.19) so that c 1 |X| 2 ≤ V 1 (t, X) ≤ c 2 |X| 2 By construction, V 1 satisfies ∂V 1 ∂t (t, X) + ∇V 1 (t, X) kA(t)X = -k|X| 2
Hence, the derivative of V 1 along the trajectories of (5.6) is

d dt V 1 (t, X) = -k|X| 2 + ∇V 1 (t, X) ξ Using |∇V 1 (t, X)| = 2k +∞ t φ(s, t) φ(s, t)dsX ≤ 2k +∞ t φ(s, t) φ(s, t)ds |X| ≤ kK 2 γ(k) |X| = 2c 2 |X|
together with inequality (5.11) yields

∇V 1 (t, X) ξ ≤ kK 2 γ(k) √ 2ω max |X| 2 + k|X| 3 Hence d dt V 1 (t, X) ≤ -k|X| 2 1 - K 2 √ 2ω max γ(k) - kK 2 γ(k) |X| -W (X)
As k > k * , we have

1 - K 2 √ 2ω max γ(k) > 0
We proceed as in [START_REF] Khalil | Nonlinear Systems[END_REF] Corollary 3.4. If the initial condition of (5.6) satisfies |X(t 0 )| < r(k) or equivalently

|X(t 0 )| < γ(k) kK 2 1 - K 2 √ 2ω max γ(k) × c 1 c 2 then W 3 (X(t 0 )) > 0 and, while W (X(t)) > 0, X(•) remains bounded by |X(t)| 2 ≤ V 1 (t) c 1 ≤ V 1 (t 0 ) c 1 ≤ c 2 c 1 |X(t 0 )| 2 which shows that W (X) ≥ k 1 - K 2 √ 2ω max γ(k) - kK 2 γ(k) c 2 c 1 |X(t 0 )| |X| 2 c 3 |X| 2 so that d dt V 1 (t, X) ≤ -c 3 |X| 2 ≤ - c 3 c 2 V 1 (t, X)
and in the end

|X(t)| 2 ≤ 1 c 1 e - c 3 c 2 (t-t 0 ) V 1 (t 0 ) ≤ c 2 c 1 e - c 3 c 2 (t-t 0 ) |X(t 0 )| 2
This shows that the equilibrium 0 of system (5.6) is locally uniformly exponentially stable. From (5.5), one directly deduce that the basin of attraction contains the ellipsoid (5.17).

Remark 5.3. The limitations imposed on â(t 0 ) and b(t 0 ) in (5.17) are not truly restrictive, as the actual values a(t 0 ), b(t 0 ) are assumed known, so the observer may be initialized with â(t 0 ) = a(t 0 ), b(t 0 ) = b(t 0 ). What matters is that the error on the unknown quantity ω(t 0 ) can be large in practice. Interestingly, when k goes to infinity r(k) tends to the limit

1 √ A max K 3 α 2 3 2 > 0
and arbitrarily large initial error ω(t 0 ) -ω(t 0 ) is thus allowed from (5.17).

Remark 5.4. The threshold k * depends linearly on ω max , which gives helpful hint in the tuning of observer (5.2).

Remark 5.5. As c 3 c 2 goes to infinity with k, the convergence rate is arbitrary fast for k large enough.

Simulation results

In this section, we illustrate the sensitivity of the observer with respect to three parameters • p which quantifies the linear independence of (a, b)

• ω max the maximal rotation rate of the rigid body

• the tuning gain k

Simulations were run for a model of a CubeSat [START_REF]The CubeSat program[END_REF]. The rotating rigid body under consideration is a rectangular parallelepiped of dimensions 20 × 20 × 10 [cm 3 ] and mass 2 [kg] assumed to be homogeneously distributed. No torque is applied on this system, which is thus in free-rotation. In this simulation, the two reference unit vectors are the Sun direction a and normalized magnetic field b. The satellite is equipped with • 6 Sun sensors providing at all times a measure of the Sun direction y a in a Sun sensor frame R s , as explained in Appendix D.1.

b(t) a(t)

• 3 magnetometers able to measure the normalized magnetic field y b in a magnetometer frame R m Typical normalized sensor outputs evolving on the unit sphere are given in Figures 5.1. Because the initial angular velocity vector is not aligned with any of the principal axes of inertia, the rotation motion is not periodic. Additive random noise have been added on each channel.

It shall be noted that, in practical applications, the sensor frames R s and R m need not coincide and can also differ from the body frame R b (defined along the principal axes of inertia) through a constant rotation R m,b , respectively R s,b . With these notations, we have

a = R m,b y a , b = R s,b y b
which is a simple change of coordinates of the measurements.

For sake of accuracy in the implementation, reference dynamics (5.1) and state observer (5.2) were simulated using Runge-Kutta 4 method with sample period 0.1 [s] for various values of p (with α = √ 1 -p) and various trajectories ω(•).

Figure 5.2 shows the convergence of the observer with parameter values corresponding to the measurements shown in Figure 5.1. Note that the vector measurement noise is filtered by the observer, thanks to the relatively low value of the gain k. the influence of p. When p gets close to 1, the rate of convergence is decreased. This was to be expected as the measured vector become collinear. To the limit, when p = 1, all the matrices A(t) become singular and the proof of convergence can not be applied anymore. In Figure 5.4 we report the behavior of the observer for increasing values of ω max . The faster the rotation, the slower the convergence. A faster convergence can be achieved by increasing the gain k. This increases the sensitivity to noise, as represented in Figure 5.5.

This observer has been tested on experimental data, as presented in Appendix D.2.

ω max = [ • /s] ω max = [ • /s] ω max = [ • /s] |ω - ω| [ • /s] t [s]

Chapter 6 A non-linear observer using single vector measurements

Observateur non-linéaire à partir de mesures d'une seule direction. Dans ce chapitre, nous donnons une preuve de convergence dans le cas où un seul vecteur de référence est disponible. La preuve repose sur une hypothèse d'excitation persistante sur la mesure a(t). Nous démontrons que cette hypothèse est valide, en rotation libre, pour presque toutes les conditions initiales. This chapter proposes an observer using single vector measurements. This addresses the case where only one direction sensor is available, and also the case p = 1 (collinearity of the two measured directions) in the Chapter 5 which led to a failure of the previous convergence proof. Interestingly, the singular case is treated with an observer having the same structure. However, its study is completely different.

The chapter is organized as follows. In Section 6.1, we introduce assumptions and the problem statement. In Section 6.2, we define the proposed non-linear observer. To prove its convergence, the error equation is identified, again, as a linear time-varying (LTV) system disturbed by a linear-quadratic term. Under a persistent excitation (PE) assumption, the LTV dynamics is shown to generate an exponentially convergent dynamics. This property, together with assumptions on the inertia parameters of the rigid body, reveal instrumental to conclude on the exponential uniform convergence of the error dynamics. Importantly, the PE assumption is proven to be automatically satisfied in the particular case of free-rotation. In details, in Section 6.3, we establish that for almost all initial conditions, the PE assumption holds. This result stems from a detailed analysis of the various types of solutions to the free-rotation dynamics. Illustrative simulation results are given in Section 6.4.

Problem statement

In this chapter we make the following assumptions. Assumption 6.1. a is constant Assumption 6.2. ω is bounded: |ω(t)| ≤ ω max at all times Assumption 6.3. J and τ (or equivalently χ = J -1 τ ) are known Assumption 6.4 (persistent excitation). There exist constant parameters T > 0 and 0 < µ < 1 such that a(•) satisfies

1 T t+T t a(s) × a(s) × ds ≥ µI, ∀t (6.1)
The problem we address in this chapter is the following.

Problem 6.1

Under Assumptions 6.1-6.2-6.3-6.4, find an estimate ω of ω from the measurements a.

Remark 6.1 (on the persistent excitation). (6.1) is equivalent to

1 T t+T t x a(s) 2 ds ≤ 1 -µ, ∀t, ∀|x| = 1 (6.2)
which is only possible if a(•) varies uniformly on every interval [t, t + T ]. Without the PE assumption, Problem 6.1 may not have a solution. For example, the initial conditions

a(t 0 ) =    1 0 0    , ω(t 0 ) =    w 0 0   
yield a(t) = a(t 0 ) for all t, regardless of the value of w. Hence, the system is clearly not observable. Such a case is discarded by the PE assumption. Note that this assumption bears on the trajectory, hence on the initial condition a(t 0 ), ω(t 0 ) and on the torque τ only.

Observer definition and analysis of convergence

Observer dynamics and properties

The observer is analog to (5.2), except that only one measurement is used. As a result, we consider a 6-dimensional augmented state.

ȧ ω = a × ω + k(a -â) E(ω) + χ + k 2 a × â (6.3)
where k > 0 is a constant (tuning) parameter. Denote S a -â ω -ω (6.4) the error state. We have

Ṡ = -kI [a × ] k 2 [a × ] 0 S + 0 E(ω) -E(ω) (6.5)
The study of the dynamics (6.5) employs a preliminary change of coordinates similar to (5.5). Denote

X a - â ω-ω k (6.6) yielding Ẋ = kA(t)X + 0 E(ω)-E(ω) k (6.7) with A(t)   -I a(t) × a(t) × 0   =   -I a(t) × -a(t) × 0   (6.8)
which we will analyze as an ideal LTV system Ẋ = kA(t)X (6.9) disturbed by the input term

ξ 0 E(ω)-E(ω) k (6.10)
We start by upper-bounding the disturbance (6.10). The inequality

|ξ| ≤ d max ( √ 2ω max |X| + k|X| 2 ) (6.11)
of Proposition 5.2 still holds.

Remark 6.2 (on the quantity d max ). As J 1 , J 2 , J 3 are the main moments of inertia of the rigid body, we have [START_REF] Landau | Mechanics. MIR Moscou[END_REF] ( §32,9)

J i ≤ J j + J k
for all permutations i, j, k and hence 0 ≤ d max ≤ 1. Moreover, d max = 0 if and only if J 1 = J 2 = J 3 . d max appears as a measurement of how far the rigid body is from an ideal symmetric body. For this reason, we call it discordance of the rigid body. Examples:

• For a homogeneous parallelepiped of size l × l × L, with L ≥ l, we have

d max = L 2 -l 2 L 2 + l 2
• For a homogeneous straight cylinder of radius r and height h we have

d max = |h 2 -3r 2 | h 2 + 3r 2

Analysis of the underlying LTV dynamics

The shape of A(t) given in (6.8) is well known in the field of adaptive control. We now proceed according to the usual procedure in this area. Along the trajectories of (6.9) we have

d dt |X| 2 = -2k|X 1 | 2 = -X C CX (6.12) with C ( √ 2kI 0)
As will be seen in the proof of the following Theorem 6.2.2, the PE assumption will imply, in turn, that the pair (kA(•), C) is uniformly completely observable (UCO), which guarantees uniform exponential stability of the LTV system (see below).

Theorem 6.1: uniform exponential stability of system (6.9) There exists 0 < c < 1, depending only on T, µ, k and ω max , such that the solution of (6.9) satisfies, for all integer N ≥ 0,

|X(t)| 2 ≤ c N |X(t 0 )| 2 , ∀t ∈ [t 0 + N T, t 0 + (N + 1)T ]
for any initial condition t 0 , X(t 0 ).

Proof. Equation (6.12) gives the result for N = 0. Further, for all t

|X(t + T )| 2 = |X(t)| 2 -X(t) 2 W (t, t + T )X(t)
where

W (t, t + T ) t+T t φ(s, t) C Cφ(s, t)ds
is the observability Gramian of the pair (kA(•), C) and φ is the transition matrix associated with (6.9). Computing W is no easy task. However, the output injection UCO equivalence result presented in [START_REF] Ioannou | Robust Adaptive Control[END_REF] allows us to consider a much simpler system. Denote

K(t) √ k √ 2 I -a(t) × and M (t) kA(t) + K(t)C = 0 k a(t) × 0 0
The observability Gramian W of the pair (M (•), C) is easily computed as

W (t, t + T ) = 2k t+T t I A(s, t) A(s, t) A(s, t) A(s, t) ds where A(s, t) k s t a(u) × du
Such a Gramian is often considered in adaptive control and has been extensively studied e.g. in [START_REF] Khalil | Nonlinear Systems[END_REF], Lemma 13.4. Here, we have

• t+T t k a(s) × k a(s) × ds ≥ T k 2 µI, ∀t • k a(•) × is bounded by k • d dt k a(•)
× is bounded by kω max from which we deduce that there exists 0 < β 1 < 1 depending on T, µ, k, ω max such that

W (t, t + T ) ≥ β 1 I, ∀t
There also exists β 2 > 0 depending on k, T such that W (t, t + T ) ≤ β 2 I. From [IS95], Lemma 4.8.1 (output injection UCO equivalence), W (t, t + T ) is also lower-bounded. More precisely, we have

W (t, t + T ) ≥ β 1 2(1 + β 2 T k) I (1 -c)I
with 0 < c < 1. Assume the result is true for an integer N ≥ 0. Then for any t ∈ [t 0 + N T, t 0 + (N + 1)T ] we have

|X(t + T )| 2 = |X(t)| 2 -X(t) W (t, t + T )X(t) ≤ c|X(t)| 2 ≤ c N +1 |X(t 0 )| 2
which concludes the proof by induction.

Convergence of the observer

Consider the quantity

d * 1 -c 2 √ 2T ω max (6.13)
where c is defined in Theorem 6.2.2. The following Theorem, which is the main result of this chapter, shows that if d max < d * , the observer (6.3) provides a solution to Problem 6.1.

Theorem 6.2: solution to Problem 6.1

We suppose that Assumptions 6.2-6.4 are satisfied and that

d max < d *
where d * is defined in (6.13). The observer (6.3) defines an error dynamics (6.5) for which the equilibrium 0 is locally uniformly exponentially stable. The basin of attraction of this equilibrium contains the ellipsoid

S(t 0 ), |a(t 0 ) -â(t 0 )| 2 + |ω(t 0 ) -ω(t 0 )| 2 k 2 < r 2 (6.14) with r 2 (1 -c) 3 8 √ 3d 2 max T 3 k 3 1 - 2 √ 2d max T ω max 1 -c 2 (6.15)
Proof. Consider the candidate Lyapunov function

V (t, X) X +∞ t φ(s, t) φ(s, t)ds X
where φ is the transition matrix of system (6.9). Let (t, X) be fixed. One easily shows that kA(t) ≤ k √ 3, for all t. Thus (see Lemma 5.1))

V (t, X) ≥ 1 2k √ 3 |X| 2 c 1 |X| 2 W 1 (X)
Moreover, Theorem 6.2.2 implies that

V (t, X) = +∞ N =0 t+(N +1)T t+N T X φ(s, t) φ(s, t)X ≤ T +∞ N =0 c N |X| 2 = T 1 -c |X| 2 c 2 |X| 2 W 2 (X) By construction, V satisfies ∂V ∂t (t, X) + ∇V (t, X) kA(t)X = -|X| 2
Hence, the derivative of V along the trajectories of the disturbed dynamics (6.7) is

d dt V (t, X) = -|X| 2 + ∇V (t, X) ξ Using |∇V (t, X)| = 2 +∞ t φ(s, t) φ(s, t)dsX ≤ 2T 1 -c |X|
together with inequality (6.11) yields

∇V (t, X) ξ ≤ 2d max T 1 -c √ 2ω max |X| 2 + k|X| 3 Hence d dt V (t, X) ≤ -|X| 2 1 - 2 √ 2d max T ω max 1 -c - 2d max T k 1 -c |X| -W 3 (X)
By assumption

d max < d * , which implies 1 - 2 √ 2d max T ω max 1 -c > 0
We proceed as in [START_REF] Khalil | Nonlinear Systems[END_REF] Corollary 3.4. If the initial condition of (6.7) satisfies |X(t 0 )| < r or equivalently

|X(t 0 )| < 1 -c 2d max kT 1 - 2 √ 2d max T ω max 1 -c × c 1 c 2 then W 3 (X(t 0 )) > 0 and, while W 3 (X(t)) > 0, X(•) remains bounded by |X(t)| 2 ≤ V (t) c 1 ≤ V (t 0 ) c 1 ≤ c 2 c 1 |X(t 0 )| 2 which shows that W 3 (X) ≥ 1 - 2 √ 2d max T ω max 1 -c - 2d max kT 1 -c c 2 c 1 |X(t 0 )| |X| 2 c 3 |X| 2 Consequently |X(t)| 2 ≤ c 2 c 1 e - c 3 c 2 (t-t 0 ) |X(t 0 )| 2
which shows that the equilibrium 0 of system (6.7) is locally uniformly exponentially stable. From (6.6), one directly deduces that the basin of attraction contains the ellipsoid (6.14).

Remark 6.3. Again, the limitations imposed on â(t 0 ) in (6.14) are not truly restrictive because, as the actual value a(t 0 ) is assumed known, the observer may be initialized with â(t 0 ) = a(t 0 ). What matters is that the error on the unknown quantity ω(t 0 ) can be large in practice.

It follows that for all t, a(t) can be written a(t) = R(t) a = cos wtR 0 a -sin wtR 0 (u × a) + (1 -cos wt)u aR 0 u (6.16) Remark 6.4. The direction u of the rotation can be simply computed from M. We have

M = RJω = J i Rω = wJ i RR 0 u = wJ i u which implies that u = M |M|
The impact of the single-axis nature of the rotation on the PE assumption is as explained in the next two subsections.

Type A solution with M aligned with a or ω = 0 Consider that a is aligned with M = R(t 0 )Jω(t 0 ). In this case, u = ±a (see Remark 6.4). Thus, (6.16) yields a(t) = R 0 a constant over time. This also holds when ω(t 0 ) = 0. For any T > 0 we have, for the unit vector

x = R 0 a 1 T 0 (a(s) x) 2 ds = 1
Thus, condition (6.2) is not satisfied.

Type A solution with M not aligned with a

Conversely, consider that a is not aligned with M. Define v, z such that (u, v, z) is a direct orthonormal basis of R 3 . The decomposition of the unit vector a in this basis is given as

a = a 1 u + a 2 v + a 3 z, a 2 1 + a 2 2 + a 2 3 = 1, with a 2 1 < 1 We have a(t) = R 0 (a 1 u + (a 2 cos wt + a 3 sin wt)v + (a 3 cos wt -a 2 sin wt)z) Note T = 2π
w . Consider any t and any unit vector x. Note x 1 , x 2 , x 3 such that

x = R 0 (x 1 u + x 2 v + x 3 z)
We have

1 T t+T t a(s) T x 2 ds = 1 T t+T t (a 1 x 1 + (a 2 cos wt + a 3 sin wt)x 2 + (a 3 cos wt -a 2 sin wt)x 3 ) 2 ds = a 2 1 x 2 1 + a 2 2 + a 2 3 2 (x 2 2 + x 2 3 ) ≤ (1 -µ) with µ min 1 -a 2 1 , 1 + a 2 1 2 ∈ (0, 1)
Thus, condition (6.2) is satisfied.

Type A' solutions

As shown in Appendix C, the Type A' solutions are characterized by J 1 > J 2 > J 3 and

d 3 |ω 1 (t 0 )| = d 1 |ω(t 0 )| = 0
which defines a zero-measure set. For this reason, they are called singular solutions. In this case, ω(t) converges to an eigenvector of J when t goes to infinity. The rotation R(t) is thus asymptotically arbitrarily close to a single-axis rotation around M = R(t 0 )Jω(t 0 ). The arguments already employed for the Type A solutions show that condition (6.2) is satisfied unless R(t 0 )Jω(t 0 ) and a are aligned.

Study of Type B and Type C solutions

In this section we show that the Type B and Type C solutions satisfy the PE assumption. Both proofs relies on the following technical result.

Proposition 6.2 (preliminary result). If condition (6.1) is not satisfied, then for all T > 0 and all ε > 0 small enough, there exists t such that for all y ∈ R 3 , and all s ∈ [t, t + T ],

• R(s)y remains between two planes orthogonal to a and distant by ε|y| • R(s) y remains between two parallel planes distant by ε|y|.

Proof. Consider T > 0 and µ such that 0 < µ < min 1 4T ω max , T ω max 4 < 1 (6.17)

Assume that (6.2) is not satisfied. Then there exists t ∈ R and x ∈ R 3 such that |x| = 1 and 1

T t+T t a(s) T x 2 ds ≥ 1 -µ (6.18)
As will appear, one can use the bounded variations of a(•) stemming from its governing dynamics to establish a lower-bound on the integrand. Denote

h(s) a(s) x 2 , ∀s
We will now show by contradiction that

h(s) ≥ 1 -2 T ω max µ, ∀s ∈ [t, t + T ]
Assume that there exists s 0 such that

h(s 0 ) < 1 -2 T ω max µ
We have, for all s,

| ḣ(s)| = 2 ȧ(s) xa(s) x = 2(a(s) × ω) xa(s) x ≤ 2ω max
Assume s 0 ≤ t + T 2 and denote

s 1 s 0 + T µ ω max ≤ t + T We have, for any s ∈ [s 0 , s 1 ] ⊂ [t, t + T ] h(s) ≤ h(s 0 ) + 2ω max (s -s 0 ) < 1 -2 T ω max µ + 2ω max (s -s 0 ) Hence 1 T t+T t a(s) T x 2 ds < 1 - 1 T T µ ω max + 1 T s 1 s 0 1 -2 T ω max µ + 2ω max (s -s 0 ) ds = 1 -2µ + µ = 1 -µ
which contradicts (6.18). The case

s 0 > t + T 2 is analog with s ∈ [s 0 - T µ ω max , s 0 ] ⊂ [t, t + T ]
Finally, we have, for all s,

0 < 1 -2 T ω max µ ≤ a(s) x 2 ≤ 1
which shows that the continuous function s → a(s) x is of constant sign, strictly positive without loss of generality. Thus, we have

0 < 1 -2 T ω max µ ≤ a(s) x ≤ 1 and in turn |a(s) -x| 2 = 2 -2a(s) x ≤ 4 T ω max µ δ √ µ (6.19)
Denote by R 1 a rotation matrix so that

a = R 1 x
and, for all s, u(s), ξ(s) such that

R(s) r u(s) (ξ(s))R 1
Note that R(s)x = r u(s) (ξ(s))a. The next Lemma formulates that the rotation R(s) is uniformly close to r a (ξ(s))R 1 .

Lemma 6.1. We have, for all s ∈ [t, t + T ] and all y ∈ R 3

|R(s)y -r a (ξ(s))R 1 y| 2 ≤ 30δ √ µ|y| 2 (6.20)
where δ is defined by (6.19).

Proof. Let s ∈ [t, t + T ].
For clarity we may omit the s dependency of u and ξ.

Denote ∆ R(s) -r a (ξ)R 1 = sin ξ ([u × ] -[a × ]) + (1 -cos ξ) uu -aa R 1 If a = u(s), ∆ = 0 ≤ 30δ √ µ.
Otherwise, for A = x we have, from (6.19)

|∆A| 2 = |R(s)x -r a (ξ)R 1 x| 2 = |R(s)x -a| 2 = x -R (s)a 2 = |x -a(s)| 2 ≤ δ √ µ Denote v, z so that (u, v, z) is an orthonormal basis of R 3 . Write a = a 1 u + a 2 v + a 3 z, a 2 1 + a 2 2 + a 2 3 = 1 We have δ √ µ ≥ |R(s)x -a| 2 = |(r u (ξ) -I)a| = |(a 2 (cos ξ -1) -a 3 sin ξ)v + (a 2 sin ξ + a 3 (cos ξ -1))z| 2 = 4(a 2 2 + a 2 3 ) sin 2 ξ 2 Now, denote B = R 1 u×a |u×a| .
We have

B A = 1 |u × a| (u × a) R 1 x = 1 |u × a| (u × a) a = 0
so that A and B are orthonormal. We have

|∆B| 2 = sin 2 ξ a 2 2 + a 2 3 |u × (u × a) -a × (u × a)| 2 = sin 2 ξ a 2 2 + a 2 3 (1 -a u) 2 |u + a| 2 = sin 2 ξ a 2 2 + a 2 3 4(1 -a 2 1 ) 2 ≤ 16(a 2 2 + a 2 3 ) sin 2 ξ 2 ≤ 4δ √ µ For C = A × B we have |∆C| 2 = |R(s)(A × B) -r a (ξ)R 1 (A × B)| 2 = |R(s)A × R(s)B -r a (ξ)R 1 A × r a (ξ)R 1 B| 2 = |R(s)A × ∆B + ∆A × r a (ξ)R 1 B| 2 ≤ 2(δ √ µ + 4δ √ µ) = 10δ √ µ
and the vectors A, B, C are orthonormal. Finally, for any vector

y = y 1 A + y 2 B + y 3 C, y 2 1 + y 2 2 + y 2 3 = 1 we have |∆y| 2 = |y 1 ∆A + y 2 ∆B + y 3 ∆C| 2 ≤ 3 y 2 1 |∆A| 2 + y 2 2 |∆B| 2 + y 2 3 |∆C| 2 ≤ 3(y 2 1 + y 2 2 + y 2 3 )10δ √ µ = 30δ √ µ|y| 2
which concludes the proof of Lemma 6.1.

Denote ε = 2 30δ √ µ and consider any y in R 3 and any s in [t, t + T ]. On the one hand, r a (ξ(s))R 1 y lies on a circle orthogonal to a. On the other hand,

|R(s)y -r a (ξ(s))R 1 y| ≤ ε 2 |y|
This yields the first item of Proposition 6.2 as µ > 0 is arbitrary small. Rewriting the result of Lemma 6.1 as

R 1 r a (-ξ(s))y -R(s) y 2 ≤ 30δ √ µ|y| 2
for any s ∈ [t, t + T ] and any y yields the second item and concludes the proof.

Type C solutions

These solutions are characterized by J 1 > J 2 > J 3 and

d 3 |ω 1 (t 0 )| = d 1 |ω 3 (t 0 )|
In this case the trajectory of ω(•) is closed and thus periodic of a certain period T 0 > 0, and not contained in a plane. Assume that condition (6.2) is not satisfied. We apply the second item of Proposition 6.2 with T = T 0 . For any ε small enough, there exists t such that for all s

∈ [t, t + T ] Jω(s) = R (s)M
remains between two parallel planes and distant by ε |M|. As ω(•) is T 0 -periodic, this is true for all s ∈ R. When ε goes to 0, we conclude that the trajectory of ω(•) remains in a plane, which is a contradiction. Thus, condition (6.2) is satisfied, unconditionally on R(t 0 ).

Type B solutions

We now consider the case where ω(t 0 ) is not an eigenvector of J and two moments of inertia are equal. Since the trajectory is a circle, it is contained in a plane and we can not apply directly the same technique as for Type C solutions. Without loss of generality, we study the case J 1 = J 2 > J 3 (the case J 1 > J 2 = J 3 is analog). We thus consider a trajectory ω such that ω(t 0 ) satisfies (ω 1 (t 0 ), ω 2 (t 0 )) = (0, 0), ω 3 (t 0 ) = 0

Following the extensive analysis exposed in [START_REF] Landau | Mechanics. MIR Moscou[END_REF], we conveniently chose the inertial frame so that e 3 is aligned with M, namely

e 3 = M |M|
For this choice of e 3 and in the case where J 1 = J 2 , equations (2.6)-(2.13) simplify considerably and one can show that the rotation matrix satisfies for all t

R(t) =    (. . . ) (. . . ) η cos ξ 1 (t -t 1 ) (. . . ) (. . . ) η sin ξ 1 (t -t 1 ) η cos ξ 2 (t -t 2 ) η sin ξ 2 (t -t 2 ) √ 1 -η 2    (6.21)
where (. . . ) designates terms that are irrelevant in the following analysis, t 1 , t 2 are constant and

η J 2 1 ω 1 (t 0 ) 2 + J 2 1 ω 2 (t 0 ) 2 J 2 1 ω 1 (t 0 ) 2 + J 2 2 ω 2 (t 0 ) 2 + J 2 3 ω 3 (t 0 ) 2 ∈ (0, 1) ξ 1 ω 1 (t 0 ) 2 + ω 2 (t 0 ) 2 + J 2 3 J 2 1 ω 3 (t 0 ) 2 > 0 ξ 2 J 3 J 1 -1 ω 3 (t 0 ) = 0
We now show that condition (6.2) is satisfied by contradiction. Assuming that it is not, one can apply the first item of Proposition 6.2 with

T = max 2π ξ 1 , 2π |ξ 2 |
For ε small enough, there exists t such that for all s ∈ [t, t + T ] R(s)e 3 remains between two planes orthogonal to a and distant by ε. Moreover, expression (6.21) yields for all s

R(s)e 3 =    η cos ξ 1 (s -t 1 ) η sin ξ 1 (s -t 1 ) √ 1 -η 2    Simple geometric considerations show that 1 -(a e 3 ) 2 ≤ ε 2η
which yields a = ±e 3 when ε goes to 0. Hence, for ε small enough, and all s ∈ [t,

t + T ] R(s)e 1 =    (. . . ) (. . . ) η cos ξ 2 (s -t 2 )   
remains between two planes orthogonal to a = ±e 3 . Taking ε < 2η yields a contradiction. The trajectories Re 1 and Re 3 are represented in Figure 6.1 for better visual understanding of the proof.

Conclusion on the PE assumption

We have shown the following result.

Theorem 6.3: PE property

Consider the vector a(t) = R(t) a

where R(t) is a rotation matrix defined as the solution of the free-rotation dynamics (2.6),(2.13). Assumption 6.4 is satisfied for almost all initial conditions (R(t 0 ), ω(t 0 )). It fails only in the cases listed below 1. ω(t 0 ) is an eigenvector of J and R(t 0 )Jω(t 0 ) is aligned with a, or 2. the eigenvalues of J are of the form J 1 > J 2 > J 3 , the coordinates of ω(t 0 ) in the trihedron of orthonormal eigendirections of J satisfy

d 3 |ω 1 (t 0 )| = d 1 |ω 3 (t 0 )|
and R(t 0 )Jω(t 0 ) is aligned with a.

It follows that, except for the initial conditions listed in items 1,2, the conclusion of Theorem 6.2.3 holds without requiring Assumption 2, which is automatically satisfied. Therefore, in almost all cases, observer (6.3) asymptotically reconstructs the desired angular rate ω.

Simulation results

In this section we illustrate the convergence of the observer and sketch the sensitivity with respect to the tuning gain k.

Simulations were run for a model of a CubeSat [START_REF]The CubeSat program[END_REF]. The rotating rigid body under consideration is a rectangular parallelepiped of dimensions about 20 × 20 × 10 [cm 3 ] and mass 2 [kg] assumed to be slightly non-homogeneously distributed. The resulting moments of inertia are

J 1 = 87 [kg.cm 2 ], J 2 = 83 [kg.cm 2 ], J 3 = 37 [kg.cm 2 ]
No torque is applied on this system, which is thus in free-rotation. Referring to Section 6.3, we will consider Type A and Type C trajectories.

In this simulation the reference unit vector is the normalized magnetic field a. The satellite is equipped with 3 magnetometers able to measure the normalized magnetic field y a in a magnetometer frame R m .

It shall be noted that, in practical applications, the sensor frame R m can differ from the body frame R b (defined along the principal axes of inertia) through a constant rotation R m,b . With these notations, we have which is a simple change of coordinates of the measurements.

a = R m,b y a
For sake of accuracy in the implementation, reference dynamics and state observer (6.3) were simulated using Runge-Kutta 4 method with sample period 0.01 [s]. The generated trajectories correspond to ω max 100 [ • /s].

Noise-free simulations

To emphasize the role of the tuning gain k, we first assume that the sensors are perfect i.e. without noises. Typical measurements for a general Type C trajectory are represented in Figure 6.2. Figure 6.3 shows the convergence of the observer for various values of k. Interestingly, large values of k produce undesirable effects. This is a structural difference with the two reference vectors based observer previously introduced in Chapter 5. The reason is that the convergence is guaranteed by a PE argument and not by a uniformly negative bound on eigenvalues.

In Figure 6.4 we represented the observer error for a case where the PE assumption is not satisfied, namely for a constant ω with M and a = (1 0 0) aligned. This is a singular case, as discussed earlier. Interestingly, the coordinates ω 2 -ω2 and ω 3 -ω3 converge to zero, while ω 1 -ω1 converges to a constant value. This can easily be proved by using LaSalle's invariance principle. Indeed, in this case, ω and a are constant and the observer system (6.3) is time-invariant.

Measurement noise

We now study the impact of measurement noises on the observer performance. The simulation parameters remain the same but we add Gaussian measurement noise with standard deviation σ = 0.03 [Hz -1 2 ]. Typical measurements are represented in Figure 6.5. Because of the noise, a sphere visualization would not be easy to apprehend. Thus, we represent the three measurement axes separately. The observer yields a residual error, about 5% in Figure 6.6 for k = 1. Note that the measurement noise is filtered, thanks to a relatively low value of the gain k. For large values of k, the observer does not converge anymore (not represented).

ω 1 -ω1 ω 2 -ω2 ω 3 -ω3 ω - ω [ • /s] t [s]

Chapter 7

Angular rate and torque observer from two vector measurements Estimateur de vitesse angulaire et de couple. Dans ce chapitre, nous écartons l'hypothèse de couple connu et étendons l'observateur pour estimer en outre un couple supposé polynomial en temps.

In this chapter, we assume that the the matrix of inertia J is known. The torque τ is unknown but slowly varying. We extend the two vector rate velocity observer (5.2) to estimate it.

The chapter is organized as follows. In Section 7.1, we introduce the problem statement. In Section 7.2, we define a non-linear observer with extended state and output injection and prove its convergence. Illustrative simulation results are given in Section 7.3. The scenarios under consideration include impulsive torques such as those generated by reaction control systems (RCS) to provide attitude control of a satellite, sublimation torques on an asteroid, and eddy current damping torques on space debris. In these scenarios, the assumption of piecewise constant or slowly varying torque is instrumental in advancing the performance of the reconstruction.

Problem statement

The assumptions of this chapter are as follows.

Assumption 7.1. The reference vectors a, b are constant. Assumption 7.2. J is known. Assumption 7.3. τ is slowly varying (so is χ). It is not known but it generates a rotation such that ω is bounded: |ω(t)| ≤ ω max at all times t. Problem 7.1 Under Assumptions 7.1-7.2-7.3, find estimates (ω, χ) of (ω, χ) from the vector measurements a, b.

Observer design and analysis of convergence

Observer definition

We extend the proposed observer (5.2) under the form

ȧ = a × ω + αk(a -â) ḃ = b × ω + αk(b -b) ω = E(ω) + χ + k 2 a × â + b × b ˙ = E(ω) + γ 1 √ k(ω -ˆ ) + χ χ = γ 2 k(ω -ˆ )                        (7.1)
where k, α, γ 1 , γ 2 > 0 are constant tuning parameters. The inputs of this observer are the two vector measurements a and b; the outputs are ŵ and χ, which are estimates of the angular rate ω and the external acceleration χ = J -1 τ .

Error dynamics

To study the convergence of this observer, we introduce the scaled errors

X =    a - â b - b ω-ω k       X 1 X 2 X 3    ∈ R 9 , Y = ω-ˆ k χ- χ k √ k Y 1 Y 2 ∈ R 6 (7.2)
which are governed by

Ẋ1 = -αkX 1 + ka × X 3 Ẋ2 = -αkX 2 + kb × X 3 Ẋ3 = k (a × X 1 + b × X 2 ) + E(ω) -E(ω) k + √ kY 2 Ẏ1 = -γ 1 √ kY 1 + √ kY 2 + γ 1 √ kX 3 + E(ω) -E(ω) k Ẏ2 = -γ 2 √ kY 1 + γ 2 √ kX 3                            (7.3)
In the next section, we show that, for sufficiently large k, system (7.3) has local uniform exponential stability, thus providing a solution to Problem 7.1. Indeed, when X 3 and Y 2 asymptotically converge to 0, we have

ω(t) -ω(t) → 0, χ -χ(t) → 0, as t → +∞
Following [START_REF] Khalil | Nonlinear Systems[END_REF]Th. 3.13], we establish the exponential stability of the linearization about the origin X = 0, Y = 0 and conclude on the non-linear dynamics. Linearization yields Ẋ = kA 1 (t)X + ξ

Ẏ = √ kA 2 Y + ζ    (7.4) with A 1 (t)      -I 0 a(t) × 0 -I b(t) × a(t) × b(t) × 0      , ξ    0 0 J ac E(ω)X 3 + √ kY 2    A 2 -γ 1 I I -γ 2 I 0 , ζ γ 1 √ kX 3 + J ac E(ω)X 3 γ 2 √ kX 3
Thus, (7.4) appears as the interconnection of two systems

Ẋ = kA 1 (t)X (7.5) Ẏ = √ kA 2 Y (7.6)
disturbed by the respective input terms ξ, ζ, as represented in Fig. 7.1.

Ẋ = kA 1 (t)X + ξ X ζ( √ k, X) ξ( √ k, X, Y ) ξ Ẏ = √ kA 2 Y + ζ ζ Y

Convergence proof

Asymptotic behavior of the X-system

We assume that 0 < α < 2 √ 1 -p and we build on the proof of Theorem 5.2.4. We have exhibited in (5.18)-(5.19) a Lyapunov function V 1 (t, X) and constants 0 < c 1 ≤ c 2 depending only on k, α and the value of the dot product p = a b, such that for all (t, X)

c 1 |X| 2 ≤ V 1 (t, X) ≤ c 2 |X| 2 |∇V 1 (t, X)| ≤ 2c 2 |X| d dt V 1 (t, X) = -k|X| 2
where the total derivative is taken along the trajectories of system (7.5).

Study of the Y -system

For any choice of γ 1 , γ 2 > 0, A 2 has eigenvalues λ 1 , λ 2 with strictly negative real parts. We assume that γ 2 1 = 4γ 2 so that λ 1 = λ 2 . Denote by v 1 , v 2 two associated eigenvectors,

P = [v 1 v 2 ]
√ kA 2 t Y | ≤ P P -1 e - √ kλt |Y | Define V 2 as V 2 (Y ) √ k Y ∞ 0 e √ kA 2 t e √ kA 2 t dt Y = √ k ∞ 0 |e √ kA 2 t Y | 2 dt
We have, for all Y (see Lemma 5.1),

V 2 (Y ) ≥ 1 2 A 2 |Y | 2 c 1 |Y | 2
and also

V 2 (Y ) ≤ P 2 P -1 2 2λ c 2 |Y | 2 (7.7)
Finally, by upper-bounding the gradient of the quadratic function V 2 , we obtain

|∇V 2 (Y )| ≤ 2c 2 |Y |
Moreover, along the trajectories of (7.6)

d dt V 2 (Y ) = - √ k|Y | 2
Remark 7.1. c 1 , c 2 depend only on the choice of the tuning parameters γ 1 , γ 2 .

Connecting the systems

We now investigate the convergence of the interconnection. Consider the candidate Lyapunov function

V (t, (X, Y )) V 1 (t, X) + V 2 (Y ) Denote C 1 min(c 1 , c 1 ), C 2 max(c 2 , c 2 ), Z X Y We have C 1 |Z| 2 ≤ V (t, Z) ≤ C 2 |Z| 2
and the derivative of V along the trajectories of (7.2) satisfies

d dt V (t, Z) = -k|X| 2 + ∇V 1 (t, X) ξ - √ k|Y | 2 + ∇V 2 (Y ) ζ Equation (5.10) yields J ac E(ω) ≤ √ 2ω max
Hence, the perturbation (coupling) terms are bounded by

|ξ| ≤ √ 2ω max |X| + √ k|Y | |ζ| ≤ √ 2ω max + (γ 1 + γ 2 ) √ k |X| It follows that d dt V (t, Z) ≤ -k|X| 2 + 2c 2 √ 2ω max |X| 2 + 2 √ kc 2 |X||Y | - √ k|Y | 2 + 2c 2 √ 2ω max + (γ 1 + γ 2 ) √ k |X||Y | = √ kZ - √ kI (c 2 + c 2 (γ 1 + γ 2 ))I (c 2 I + c 2 (γ 1 + γ 2 ))I -I Z + 2c 2 √ 2ω max |X| 2 + 2c 2 √ 2ω max |X||Y | (7.8)
Interestingly, for sufficiently large k, the symmetric matrix

- √ k c 2 + c 2 (γ 1 + γ 2 ) c 2 + c 2 (γ 1 + γ 2 )
-1 For any choice of γ 1 , γ 2 > 0, the observer (7.1) defines an error (7.2) which, for sufficiently large k > 0, converges locally uniformly exponentially to 0.

is
Remark 7.2. As C 2 is bounded, the convergence is arbitrarily fast when k grows to infinity. For applications, choosing large values for k increases the sensitivity to noise, so a natural recommendation is to consider only reasonable values for k.

Observer extensions

To account for more general torque models (e.g. piecewise affine time-varying signals, as considered below) such as χ = χ 1 , χ1 = 0 the observer can be further extended into the form ȧ = a × ω + αk(a -â)

ḃ = b × ω + αk(b -b) ω = E(ω) + χ + k 2 a × â + b × b ˙ = E(ω) + γ 1 √ k(ω -ˆ ) + χ χ = χ1 + γ 2 k(ω -ˆ ) χ1 = γ 3 k 3 2 (ω -ˆ )                                (7.9)
where k, γ 1 , γ 2 , γ 3 > 0 are constant tuning parameters. The convergence analysis is identical. The generalization to signals having a zero n-th order derivative χ = χ 1 , χ1 = χ 2 , . . . , χn = 0 is straightforward.

Simulation results on aerospace applications

In this section, we illustrate the observer in various examples from the field of aerospace engineering. For the sake of accuracy of implementation, the reference dynamics (2.7),(2.17) and state observer (7.1) were simulated using the Runge-Kutta 4 method. The sampling (measurement) period is specified in each example.

Estimating torques during a satellite reorientation maneuver (Heaviside step torques)

The first problem we address is the estimation of torques during a satellite reorientation maneuver. We consider a satellite which is a parallelepiped of size 90 gain k is seen in Fig. 7.4, which shows the observer performance for increasing values of k and fixed values γ 1 = 1, γ 2 = 0.2. As expected, the convergence is better for larger values of the gain. Thus, the tuning of k will result from a trade-off between performance and noise reduction. Finally, Fig. 7.5 stresses the role of the discussed observer extension. Torque signals that are not piecewise constant but slowly varying are well estimated using the extension, by reducing asymptotic bias.

×

Estimating sublimation torques on an asteroid (impulsive and noisy torques)

In this second example, we wish to estimate a torque signal with no particular signature except that it is non-zero for short periods only. This property is representative of torques produced by sublimation on asteroids, that is, the sudden transformation into water vapor of ice on the surface of an asteroid upon exposure to sunlight. Very generally, the monitoring of rotation states of asteroids and comets has been the subject of numerous studies. This subject is of importance as the spin of an asteroid can affect its orbit in quite a subtle way due to the Yarkovsky effect [START_REF] Rubincam | Radiative spin-up and spin down of small asteroids[END_REF], causing uncertainty in the prediction of future trajectories. Usually, the time-series analysis of emitted light, or "light curves" [START_REF] Lorenz | Spinning Flight[END_REF], is employed to determine the spin. The frequent cases of pure single-axis mode rotation is (to some extent) a resolved problem; however, numerous recent studies have reported that several comets and asteroids have more complex rotations (e.g., 1P/Halley, and 4179 Toutatis [START_REF] Mueller | The diagnosis of complex rotation in the lightcurve of 4179 Toutatis and potential applications to other asteroids and bare cometary nuclei[END_REF][START_REF] Samarasinha | Long-term evolution of rotational states and nongravitational effects for Halley-like cometary nuclei[END_REF]). In such cases, the light curves are difficult to analyze. A suggested alternative technique is as follows: assuming that a set of suitable sensors could be attached to the object (as achieved by the Rosetta lander [UBD + 15]) and that their measurements could be processed, at least for some time interval, it would be possible to reconstruct the torques produced by cometary activity, in particular the sublimation-induced torques, using the approach proposed here. Figure 7.6 shows the estimated torque for a simulated Halley-like comet considered here as an ellipsoid with semi-axes 16 × 8 × 7.999 [km 3 ] and mass 2.10 14 [kg] homogeneously distributed, and an initial rotation period of 12 hours. The magnitude of the torque was chosen in accordance with the calculations described in [START_REF] Samarasinha | Long-term evolution of rotational states and nongravitational effects for Halley-like cometary nuclei[END_REF]. For realism, some measurement noise was added, which explains the shape of the dashed curve.

Estimating eddy-current braking in space debris (linear damping torques)

This third example considers a slowly varying torque, proportional to the angular velocity of the body. This is representative of torques produced by magnetic braking inside space debris. Space debris orbiting in the magnetosphere are subject to magnetic braking stemming from eddy currents created in their spinning bodies. As developed in [PPBLV11, PHB + 12], the torques driving the spin motion of an empty, thin-walled body (e.g., the H10 stage of the Ariane 4 rocket), are of the form

τ    T 1 0 0 0 T 2 0 0 0 T 3    ω
where T 1 , T 2 , T 3 are positive parameters proportional to the magnitude of the local magnetic field. These parameters also depend on the rotational axis of the long-term rotation of the object (axial spin or flat spin). Along typical orbits, the torques defined above yield a quasi-linear asymptotic decay, with decay times ranging from 20 days (axial spin) to 250 days (flat spin). [START_REF] Praly | Study of the eddy current damping of the spin dynamics of spatial debris from the Ariane launcher[END_REF] conducted a finite element method simulation to estimate the induction phenomenon, in which the model of the spinning object was simplified to a symmetrical cylinder with semi-spherical ends. Certainly, some error was generated by this approximation, and it would be interesting to use experimental data to improve the model. The error we take into consideration appears in the torque generation law, where we assume that the torque is actually linearly dependent on the angular velocity through an unknown matrix T τ = T w

We can modify the above observer to estimate the whole matrix T , or equivalently P J -1 T . To address this problem, a further extension of the proposed observer is needed, with the idea of interconnecting the X system (7.5) with another exponentially stable system. In detail, while the Euler equations become ω = E(ω) + P ω we propose the following observer

ȧ = a × ω + αk(a -â) ḃ = b × ω + αk(b -b) ω = E(ω) + P ω + k 2 a × â + b × b ˙ = E(ω) + P ω + γ 1 (ω -ˆ ) Ṗ = γ 2 (ω -ˆ )ω                    (7.10) with k, γ 1 , γ 2 > 0. Consider the (scaled) observer error X    a - â b - b ω-ω k       X 1 X 2 X 3    ∈ R 9 , Y 1 k ω -ˆ P -P Y 1 Y 2 ∈ R 3 × R 3×3
Note that Y 2 is a 3 × 3 matrix. The linearized error system around X = 0, Y = 0 is

Ẋ = k      -αI 0 a(t) × 0 -αI b(t) × a(t) × b(t) × 0      X +    0 0 J ac E(ω)X 3 + P X 3 + Y 2 ω    Ẏ1 = -γ 1 Y 1 + Y 2 ω + J ac E(ω)X 3 + P X 3 + γ 1 X 3 Ẏ2 = -γ 2 Y 1 ω + γ 2 X 3 ω
In this form, the error system appears as an interconnection of the exponentially stable system (7.5) with the Y system

Ẏ1 = -γ 1 Y 1 + Y 2 ω Ẏ2 = -γ 2 Y 2 ω (7.11)
The convergence analysis is similar: establish the exponential stability of the Y system, and derive conditions such that the interconnection maintains stability. Briefly, a candidate Lyapunov function for system (7.11) would be

V (Y ) = γ 2 2 |Y 1 | 2 + γ 1 2 Tr Y 2 Y 2
(where Tr (.) is the trace operator) which is positive definite and satisfies

V (Y ) = -γ 1 γ 2 |Y 1 | 2
along the trajectories of (7.11). Uniform exponential stability requires a PE-like assumption on ω guaranteeing that ω(t) persistently reaches all the directions in R 3 . In the absence of measurement noise, the relative error for P decreases from 10 6 to 10 -10 within a few hours. With a noise level of 10 -7 , the relative error reaches the asymptotic value of 10 -1 , meaning that the coefficient of the P matrix (or T ) is reconstructed with a 10% error.

Chapter 8

Angular rate and ratios of inertia observer from two vector measurements Estimateur de vitesse angulaire et de ratios d'inertie. Dans ce chapitre, nous écartons l'hypothèse que J est connu et étendons l'observateur pour estimer en outre les ratios d'inertie.

In this chapter, we assume that the torque τ is known and that the matrix of inertia J is unknown. We extend the two vector angular rate observer (5.2) to further estimate the ratios of inertia (d 1 , d 2 , d 3 ).

The chapter is organized as follows. In Section 8.1, we introduce the assumptions and the problem statement. In Section 8.2, we define a non-linear observer with extended state and output injection and prove its convergence. Illustrative simulation results are given in Section 8.3.

Problem statement

In this Chapter we denote

D(ω)    ω 2 ω 3 0 0 0 ω 3 ω 1 0 0 0 ω 1 ω 2   
so that the Euler's equation can conveniently be written

ω = E(ω) + χ = D(ω)d + χ
The assumptions of this chapter are as follows. In the next section, we show that, for sufficiently large k, system (8.5) demonstrates local uniform exponential stability, thus providing a solution to Problem 8.1. Indeed, when X 3 and Y 2 go to 0, we have

ω(t) -ω(t) → 0, d -d(t) → 0, as t → +∞
Following [Kha96, Th. 3.13], we establish the exponential stability of the linearization about the origin X = 0, Y = 0 and conclude on the non-linear dynamics. Linearization yields Ẋ = kA 1 (t)X + ξ

Ẏ = A 2 Y + ζ    (8.6) with A 1 (t)      -I 0 a(t) × 0 -I b(t) × a(t) × b(t) × 0      , ξ    0 0 J ac E(ω)X 3 + D(ω)Y 2    A 2 -γ 1 I D(ω) -γ 2 D(ω) 0 , ζ γ 1 X 3 + J ac E(ω)X 3 γ 2 D(ω)X 3
Thus, (8.6) appears as the interconnection of two systems

Ẋ = kA 1 (t)X (8.7) Ẏ = A 2 Y (8.8)
disturbed by the respective input term ξ, ζ, as represented in Fig. 8.1.

Ẋ = kA 1 (t)X + ξ X ζ(X) ξ(X, Y ) ξ Ẏ = A 2 Y + ζ ζ Y

Convergence proof

Asymptotic behavior of the X-system

It is exactly the same as in Section 7.2.3.

Study of the Y -system

Following a similar analysis as in Section 6.2.2, one can show that system Y is uniformly exponentially stable thanks to the PE assumption 8.4. Thus, there exists a function V 2 (t, Y ) and constants 0 < c 1 ≤ c 2 such that, for all Y (see [START_REF] Khalil | Nonlinear Systems[END_REF] Theorem 3.12)

c 1 |Y | 2 ≤ V 2 (t, Y ) ≤ c 2 |Y | 2 ∇V 2 (t, Y ) ≤ 2c 2 |Y | d dt V 2 (t, Y ) = -|Y | 2
where the total derivative is taken along the trajectories of (8.8).

Connecting the systems

We now investigate the convergence of the interconnection. Consider the candidate Lyapunov function

V (t, (X, Y )) V 1 (t, X) + V 2 (t, Y ) Denote C 1 min(c 1 , c 1 ), C 2 max(c 2 , c 2 ), Z X Y We have C 1 |Z| 2 ≤ V (t, Z) ≤ C 2 |Z| 2
and the derivative of V along the trajectories of (8.4) satisfies

d dt V (t, Z) = -k|X| 2 + ∇V 1 (t, X) ξ -|Y | 2 + ∇V 2 (Y ) ζ
Using inequality (8.2) together with

J ac E(ω) ≤ √ 2ω max
yields the following bound on the coupling terms

|ξ| ≤ √ 2ω max |X| + ω 2 max |Y | |ζ| ≤ √ 2ω max + γ 1 + γ 2 ω 2 max |X| It follows that d dt V (t, Z) ≤ -k|X| 2 + 2c 2 √ 2ω max |X| 2 + 2c 2 ω 2 max |X||Y | -|Y | 2 + 2c 2 √ 2ω max + γ 1 + γ 2 ω 2 max |X||Y | = Z (-k + 2c 2 √ 2ω max )I ηI ηI -I Z P Z (8.9) with η c 2 ω 2 max + c 2 √ 2ω max + γ 1 + γ 2 ω 2 max
Interestingly, for sufficiently large k, the symmetric matrix P is definite negative. Under this condition, d dt V is definite negative and system (8.4) is uniformly exponentially stable.

Chapter 9 Perspectives

Perspectives. Dans ce chapitre nous présentons des pistes de recherche pour élargir les résultats exposés dans le reste du manuscrit. Ces pistes ont été abordées lors de la soutenance, notamment en réponse aux questions du jury.

Adaptive gain tuning

Observer (5.2) allows to estimate the angular rate ω for a sufficiently large value of k. The convergence rate increases with k, but so does the sensitivity to noise. One way to solve this trade-off would be to chose k adaptively, building on the results presented in [START_REF] Ioannou | Robust Adaptive Control[END_REF].

Internal torques

The mass of a satellite evolves in time with its fuel consumption. This induces a change of the matrix of inertia, that we have not taken into account in this manuscript. A solution to tackle this problem is to model this effect as internal torques in addition to the external torque τ .

Observer on S 2

The measurement vectors a, b lie on the unit sphere S 2 . One criticism that could be done to observer (5.2) is that â, b are not constrained to remain on S 2 (although they asymptotically converge to it). Thus, the observer does not account for the intrinsic geometry of the state. We could rewrite it to take this fact into account, constructively.

Nonlinear small-gain design

The proofs of local convergence of the observers presented in Chapters 7 and 8 rely on a linearization argument. This does not allow to estimate the basin of attraction of the 0 equilibrium point. A more detailed analysis, such as the one performed in Chapters 5 and 6 could be done to estimate this basin of attraction.

Chapter 9. Perspectives Figure 9.1: The smart-phone is to be inserted in a wooden ball via a hole.

Varying reference

In real-life applications, it may happen that the reference vectors a, b vary in time. For example the Earth magnetic field, which mostly depends on the latitude, is varying along the orbit of a satellite. Under this assumption, the measurements equation (2.17) becomes

ȧ(t) = a(t) × ω(t) + R(t) ȧ(t)
Assuming that ȧ is known, we can still estimate ω if we add R in the augmented state. This can be done by merging the angular rate observer (5.2) with an attitude observer such as presented e.g. in [START_REF] Trumpf | Analysis of non-linear attitude observers for time-varying reference measurements[END_REF].

Experiments

We are conducting experiments to validate the convergence of observer (5.2) using the sensors available on the smart-phone presented in Appendix D.2. In the following we give a brief report of one of those experiments.

The smart-phone is inserted inside a wooden ball comprising a hole. Figure 9.1 pictures the smart-phone partially inserted. The hole is then plugged with foam to ensure stability of the smart-phone. The rigid body under consideration is the system {ball + smart-phone}. This device gives the user a lot of flexibility on the choice of rotation motion, as it can freely roll on the ground, according to its initial conditions.

In details, we implement observer (5.2) with a the data produced by the accelerometer, b the data produced by the magnetometer. The tuning parameters are set to α = 1, k = 5. In the observer implementation, the torque governing the dynamics is arbitrarily set to τ = 0, which assumes that the rigid body is in free-rotation. To compute the matrix of inertia, we roughly assimilate the system to a homogeneous ball of radius 15 [cm] and mass density 0.69. This yields

J =    0.0878 0 0 0 0.0878 0 0 0 0.0878   
In this expression we have neglected the contribution of the smart-phone to the matrix of inertia.

The ball is then given some initial velocity and left rolling on the ground until it is stopped. Once the data are recovered from the smart-phone, they are processed in pseudo real-time (i.e. from initial time to end-time) by the observer. The estimation results can be compared against the measurements of the angular rate given by the rate gyros of the smart-phone. In accordance with the theory, the gyros data are not used in the observer, they only serve as reference for final comparisons.

The estimation results are reported in Figure 9.2. It appears that the estimates are unbiased (though a bit noisy).

u 1 = 0 1 -1 0 x 2 -x 1 |x 2 -x 1 | , g 1 = x 1 u 1 . . . u n = 0 1 -1 0 x 1 -x n |x 1 -x n | , g n = x n u n Remark B.2.
In general the Chebyshev center of P is not unique. For a rectangle e.g., there are infinitely many Chebyshev centers.

B.2 Jacobi elliptic functions

In this appendix we properly define the Jacobi elliptic functions and state some useful properties. Consider m ∈ [0, 1) and φ ∈ R, or m = 1 and φ ∈ (-π 2 , π 2 ). Define t(φ) as

t(φ) = φ 0 dζ 1 -m sin 2 ζ t(•)
is a strictly increasing one-to-one mapping from R (or (-π 2 , π 2 ) for m = 1) to R. Denote (abusively) by φ(t) its reciprocal function. The Jacobi elliptic functions of parameter m are defined for all t ∈ R as cn(t) cos φ(t), sn(t) sin φ(t), dn(t) 1 -m sin 2 φ(t)

We have represented in

Remark.

A direct corollary of this proposition is that for all t, (cn(t), sn(t), dn(t)) remains in the intersection of two cylinders defined in R 3 by

x 2 1 + x 2 2 = 1, mx 2 2 + x 2 3 = 1 The next proposition formulates a reciprocal result. Proposition B.2. Consider x 1 > 0, x 2 ∈ R and x 3 ≥ 0 satisfying x 2 1 + x 2 2 = 1, mx 2 2 + x 2 3 = 1 There exists t 1 ∈ R such that cn(t 1 ) = x 1 , sn(t 1 ) = x 2 , dn(t 1 ) = x 3 If m < 1, it is also true for x 1 = 0.

Proof. The equality x 2 1 + x 2 2 = 1 implies that there exists φ 0 ∈ R such that x 1 = cos φ 0 , x 2 = sin φ 0

If x 1 > 0, one can choose φ 0 ∈ (-π 2 , π 2 ). φ(•) is a one-to-one mapping from R to R (or to (-π 2 , π 2 ) for m = 1). Hence, there exists t 1 ∈ R such that φ(t 1 ) = φ 0 , which yields x 1 = cos φ(t 1 ) = cn(t 1 ), x 2 = sin φ(t 1 ) = sn(t 1 ) Finally x 3 ≥ 0 yields mx 2 2 + x 2 3 = 1 ⇒ x 3 = 1 -mx 2 2 = 1 -m sn(t 1 ) 2 = dn(t 1 ) which concludes the proof.

1 0 dx √ 1 -x 2 √ 1 -mx 2 = t(φ) + T 2
where we have used the fact that sin 2 is π-periodic and even, and the change of variables x = sin ζ. Thus, φ(•) satisfies, for all t ∈ R φ(t + T 2 ) = φ(t) + π, φ(t + T ) = φ(t) + 2π

The conclusion follows immediately.

Remark. In the limit case m = 1, the Jacobi elliptic functions are not periodic. φ(t) converges to ± π 2 when t goes to ±∞ and we have the respective limits lim t→±∞ cn(t) = 0, lim t→±∞ sn(t) = ±1, lim t→±∞ dn(t) = 0 Further properties on Jacobi elliptic functions, such as numerical approximations, can be found in [AS64] §16.

These functions allow one to compute an analytical solution to the free Euler's equations as explained in Appendix C.

B.3 A result on LTV systems

In this appendix we give a proof of [START_REF] Hill | Exponential stability of time-varying linear systems[END_REF] Theorem 2.1. The claim is as follows: consider a LTV system Ẋ = M (t)X such that • M (•) is l-Lipschitz, with l > 0 • there exists K ≥ 1, c ≥ 0 such that for any t and any s ≥ 0, e M (t)s ≤ Ke -cs Then, for any t 0 , X 0 , the solution of Ẋ = M (t)X with initial condition X(t 0 ) = X 0 satisfies, for any t ≥ t 0 ,

|X(t)| ≤ Ke ( √ Kl ln 2-c)(t-t 0 ) |X 0 |
Proof. Without loss of generality we can take t 0 = 0. The result is obvious for X = 0. Consider a non-zero solution (i.e. X 0 = 0) and define u(t) and v(t) as We start by the following Lemma. in which case, since ω(0) is not an eigenvector of J, ω 1 (0) and ω 3 (0) differ from 0.

u
m 1 d 1 ω 1 (0) 2 + 1 |d 2 | ω 2 (0) 2 1 d 3 ω 3 (0) 2 + 1 |d 2 | ω 2 (0) 2 (C.1) A 1 sign(ω 1 (0)) ω 1 (0) 2 + d 1 |d 2 | ω 2 (0) 2 = 0 A 2 |d 2 | d 1 |A 1 | = |d 2 | d 1 ω 1 (0) 2 + ω 2 (
The normalized functions defined by

Ω i (t) 1 A i ω i t w , i = 1, 2, 3
satisfy the following differential system (for brevity we omit the t dependency of the Ω i and the t w dependency of the ω i )

Ω1 = 1 wA 1 ω1 = d 1 wA 1 ω 2 ω 3 = - ω 2 ω 3 A 2 A 3 = -Ω 2 Ω 3 Ω2 = 1 wA 2 ω2 = d 2 wA 2 ω 3 ω 1 = ω 3 ω 1 A 3 A 1 = Ω 3 Ω 1 Ω3 = 1 wA 3 ω3 = d 3 wA 3 ω 1 ω 2 = -m ω 1 ω 2 A 1 A 2 = -mΩ 1 Ω 2
which is exactly system (B.1). Moreover, their initial conditions satisfy

Ω 1 (0) 2 + Ω 2 (0) 2 = ω 1 (0) 2 A 2 1 + ω 2 (0) 2 A 2 2 = 1 A 2 1 ω 1 (0) 2 + d 1 |d 2 | ω 2 (0) 2 = 1 mΩ 2 (0) 2 + Ω 3 (0) 2 = m ω 2 (0) 2 A 2 2 + ω 3 (0) 2 A 2 3 = m A 2 1 d 1 |d 2 | ω 2 (0) 2 + d 1 d 3 ω 3 (0) 2 = 1
Finally, we have Ω 1 (0), Ω 3 (0) ≥ 0 and if m = 1, as stated in Remark C.1, ω 1 (0) = 0, which implies that Ω 1 (0) > 0. Proposition B.2 guarantees that there exists t 1 ∈ R such that Ω 1 (0) = cn(-t 1 ), Ω 2 (0) = sn(-t 1 ), Ω 3 (0) = dn(-t 1 ) Thus, (Ω 1 , Ω 2 , Ω 3 ) and (cn(• -t 1 ), sn(• -t 1 ), dn(• -t 1 )) satisfy the same differential system (B.1) and the same initial conditions for t = 0. By Cauchy-Lipschitz uniqueness theorem, they are equal and we have for all t ω 1 (t) = A 1 cn(w(t -t 1 )) ω 2 (t) = A 2 sn(w(t -t 1 )) ω 3 (t) = A 3 dn(w(t -t 1 ))

Remark. The periodicity of the Jacobi elliptic functions 1 implies that ω 1 and ω 2 are T 0 -periodic, ω 3 is T 0 2 -periodic, where T 0 is defined by

T 0 4 |w| 1 0 1 √ 1 -x 2 √ 1 -mx 2 (C.3)
where m and w are respectively defined by (C.1) and (C.2).

This case is called long-axis mode 2 (LAM) [START_REF] Mueller | The diagnosis of complex rotation in the lightcurve of 4179 Toutatis and potential applications to other asteroids and bare cometary nuclei[END_REF]. We just have to exchange the roles of indices 1 and 3 to find ω 1 (t) = A 1 dn(w (t -t 1 )) ω 2 (t) = A 2 sn(w (t -t 1 )) ω 3 (t) = A 3 cn(w (t -t 1 ))

where A 1 , A 3 are described in the previous case and

m 1 d 3 ω 3 (0) 2 + 1 |d 2 | ω 2 (0) 2 1 d 1 ω 1 (0) 2 + 1 |d 2 | ω 2 (0) 2 ∈ (0, 1) A 2 |d 2 | d 3 ω 3 (0) 2 + ω 2 (0) 2 w -sign(ω 1 (0)ω 3 (0)) d 3 |d 2 |ω 1 (0) 2 + d 1 d 3 ω 2 (0) 2 = 0
with the convention sign(0) = 1. This case is called short-axis mode (SAM) [START_REF] Mueller | The diagnosis of complex rotation in the lightcurve of 4179 Toutatis and potential applications to other asteroids and bare cometary nuclei[END_REF].

Remark. Again, ω 1 , ω 2 are T 0 -periodic, ω 3 is T 0 2 -periodic, where T 0 is defined by

T 0 4 |w | 1 0 1 √ 1 -x 2 √ 1 -m x 2

C.3.3 Conclusion

To sum up, in the case where We call these solutions Type A' as they converge to an eigenvector of J, which is a Type A solution. The set defined by (C.5) has zero Lebesgue measure and such solutions are never observed in practice. For this reason, they are called pathological solutions.

C.4 Graphical representations

We have represented in Figures C.1 and C.2 the four types of solution. For convenience of visualization, we chose solutions having the same angular momentum norm |M|, so that they evolve on the surface of the same ellipsoid, as defined in (2.15). function of this angle, namely proportional to a n, as long as the sensor is exposed to the light. If the photocell lies in the shadow of the rigid body, the output signal is zero. Hence, one Sun sensor yields (after normalization) "half" a scalar information max(a n, 0) Consider now 6 photocells distributed onto the rigid body whose outward-pointing normals are the 6 body frame directions ±b 1 , ±b 2 , ±b 3 , as represented in The simple cosine description model is known to be incomplete [START_REF] Larson | Space mission analysis and design[END_REF][START_REF] Winetraub | Attitude determinationadvanced sun sensors for pico-satellites[END_REF] as it does not account for sensor non-linearities. We verified this fact on the coarse Sun sensor TSL250 from Texas Instrument. This is a simple light to voltage sensor consisting in a photodiode and a transimpedance amplifier, represented in Figure D.2.

We conducted the following experiment. Four sensors are mounted on a wheel, disposed in a dark room with a light source, as represented in Figure D.3.

The wheel is mounted on a gear connected to a motorized rail (not represented), which allows to accurately control its phase position ψ or velocity ψ. If one uses the raw measurements from a complete rotation of the wheel to feed equation (D.1), then the curve presented in Figure D.4 is obtained. There is a large difference with the theoretical circle that should be obtained. The experimental setup allows us to establish the response of every cell. To do so, we observed the following protocol. This thesis addresses the general question of estimating the angular rate of a rigid body from on-board direction sensors (e.g. magnetometers, Sun sensors). The objective is to replace rate gyros which are very expensive compared to direction sensors, prone to saturation during high rate rotations and subject to failure. In a first part of the thesis, we address the specific cases of single-axis and slightly perturbed axis rotations. In a second part, we address the general case by an asymptotic non-linear observer. We build the observer from two non-collinear vector measurements or from a single vector measurements. The knowledge of the inertial coordinates of the reference vectors is not necessary. We then extend the observer to further estimate unknown torques and inertia parameters. The Euler's equations play a central role in all the works developed in this thesis. It appears that, at least for the illustrative cases considered, rate gyros could be replaced with an estimation algorithm employing direction sensors which are much cheaper, more rugged and more resilient sensors.

Keywords: estimation, angular rate, Euler's equations, signal processing

Figure 1 . 1 :

 11 Figure 1.1: Schematic representation of an indirect derivative approach.

Figure 1 . 2 :

 12 Figure 1.2: Schematic representation of an indirect estimation approach

  Figure 1.3: Schematic representation of a direct derivative approach.

  Figure 1.4: Schematic representation of an direct estimation approach.

  post-filtering[START_REF] Thienel | Comparison of angular velocity estimation methods for spinning spacecraft[END_REF] Part I: phase and modulus analysis Part I: frequency analysis non-linear observer[START_REF] Jorgensen | Observer based sliding mode attitude control: Theoretical and experimental results[END_REF]

Figure 1

 1 Figure 1.6: Schéma d'une approche dérivative indirecte. Le filtre actif peut être un filtre de Kalman étendu (EKF) ou l'une de ses nombreuses variantes [ABIH98, ABIDH99, HBI99, TM12].

Figure 1

 1 Figure 1.7: Schéma d'une approche d'estimation indirecte

  Wahba. Ces méthodes ne nécessitent pas la connaissance des coordonnées inertielles a, b des vecteurs de référence. Ceux-ci sont sont généralement considérés comme constant, de sorte que les mesures satisfont l'équation différentielle ȧ = a × ω (1.11) A nouveau, nous classifions ces méthodes en deux catégories : approche dérivative et approche d'estimation. 1. L'approche dérivative directe consiste à résoudre une équation algébrique provenant de (1.11) pour estimer ω. Par exemple, lorsqu'on dispose de mesures de deux vecteurs de référence non-colinéaires unitaires a, b, ω satisfait ȧ × a + ḃ × b = P (a, b)ω où P (a, b) est la matrice symétrique définie positive suivante P (a, b) I -aa + I -bb

  Figure 1.8: Schéma d'une approche dérivative directe.

Figure 1 . 9 :

 19 Figure 1.9: Schéma d'une approche d'estimation directe.

Figure 3

 3 Figure 3.1: a(t) describes a curve in a plane orthogonal to u.

Figure 3 . 4 :

 34 Figure 3.4: Curves γ k and γk .

  For example, we can try to compute a point as far as possible from the polygonal line drawn by the values y[k], namely its Chebyshev center. It is defined as the center of the maximal-radius disk inscribed into the polygonal line defined by the y[k] and presented in details in Appendix B.1. Denote by Γ k the polygon of ordered vertices y[1], . . . , y[k], y[1]. If Γ k is convex, one can compute the Chebyshev center of Γ k with convex optimization methods (linear programming). If not, the task is more difficult. We propose to pick a subset of y[1], . .

Figure 3

 3 Figure 3.7: The Chebyshev center position (bubble) does not depend on the density of points distribution, but rather on the geometric shape of the samples. Therefore, it remains in the allowed origins set, while the average (cross) is dragged out by a non-homogeneous distribution of the ψ[k] [simulation data].

2

  and for z 0 in A ρ,∆ , (3.21) gives an estimate of ψ[k]. For z 0 = C, the error is bounded by

Figure 3

 3 Figure 3.8: ψ[k] matches ψ[k] from the start to the end of the maneuver [simulation results].

Figure 4

 4 Figure 4.1: Representation of the 'ZXZ' Euler angles.

Figure 4

 4 Figure 4.2: Left: a(t) draws circles slowly moving up the unit sphere. Right: projection onto a horizontal plane.

Figure= θ cos ψ ω 2 = -θ sin ψ ω1 , ω2 Figure 4 . 3 :

 ω243 Figure 4.3: Estimation of ψ, θ, ω under Assumption 4.2.

Figure 4

 4 Figure 4.4: Trajectory of the (center of gravity of the) projectile [simulation data].

Figure 4 Figure 4 Figure 4

 444 Figure 4.5: The small precession ϕ approximation is relevant [simulation data].

  0) When g is e.g. a normalized Hann window [Mal98], we have ∆ν 9.05 [rad], C 1 2.83, C 2 14.61 If the conditions (4.9)-(4.10)-(4.11) hold, then the study of the local maxima (their arguments and value) of Sy will give a convenient and reliable solution to the stated problem. In detail, we use the following algorithm, represented in Figure 4.8 for convenience. Algorithm. Solution to Problem 4.1 under Assumption 4.3 and conditions (4.8)-(4.9)-(4.10)-(4.11) Inputs: • sampled data y[k], sampling time t[k] • window function g • window size s Steps for each t[k]: • calculate the spectrogram Sy(t[k], •) (e.g. by mean of a fast Fourier Transform) • find the two main lobes of |Sy(t[k], •)| corresponding to frequencies ν 1 , ν 2 with |ν 1 | > |ν 2 | and amplitudes A 1 , A 2 (e.g. by an exhaustive search of local maxima)

Figure 4 . 8 :

 48 Figure 4.8: Estimation of ϕ, θ, ψ, ω under Assumption 4.3 and conditions (4.8)-(4.11).

Figure 4

 4 Figure 4.9: Euler angles [simulation data].

  5 [m], 0.75 [m], 1 [m] and mass 200 [kg]. We use the following parameters

Figure 4

 4 Figure 4.11: Angular rate (left) and relative estimation error (right) [simulation].

Figure 4 .

 4 Figure 4.11 together with the corresponding relative estimation error. The estimation error remains below 10%. It worsens when the scaling factor s decreases. To explain this fact, we compute the corresponding spectrograms |Sy| in Figure 4.12 . In both figures, one can clearly see the two main lobes of the spectrogram and the peaks corresponding to frequencies ψ and ϕ + ψ. For s = 3.5 [s], the peaks corresponding to ψ are corrupted by side oscillations of the main lobe. The estimation becomes less accurate.

Figure 4 .

 4 Figure 4.12: Spectrogram for s = 6 [s] (left) and s = 4 [s] (right).

  Figure 4.13: Angular rate (left) and relative estimation error (right) [simulation].

Figure 4 .

 4 Figure 4.14: Spectrogram for s = 5 [s] (left) and s = 3.5 [s] (right).

  Assumption 5.1. a, b are constant and linearly independent Assumption 5.2. J and τ (or equivalently χ = J -1 τ ) are known Assumption 5.3. ω is bounded : |ω(t)| ≤ ω max at all times Assumption 5.1 implies that p a b = a b

Figure 5 . 1 :

 51 Figure 5.1: Normalized measurements a(t) (Sun measurements) and b(t) (magnetic measurements) during rotation motion [simulation data].

FigureFigure 5

 5 Figure 5.2 shows the convergence of the observer with parameter values corresponding to the measurements shown in Figure 5.1. Note that the vector measurement noise is filtered by the observer, thanks to the relatively low value of the gain k.Figure 5.3 shows

Figure 5

 5 Figure 5.3: k = 0.25, ω max = 6 [ • /s]. The rate of convergence degrades when p increases [simulation results].

Figure 5

 5 Figure 5.4: k = 0.25, p = 0.2. Impact of ω max on the convergence rate [simulation results].

Figure 5

 5 Figure 5.5: ω max = 62 [ • /s], p = 0.2. When k increases, the convergence is faster but the measurement noise filtering degrades [simulation results].

Figure 6 . 1 :

 61 Figure 6.1: Re 3 and Re 1 evolving on the unit sphere.

Figure 6 . 2 :

 62 Figure 6.2: Typical measurements in the ideal noise-free case [simulation data].

Figure 6 . 3 :

 63 Figure 6.3: Convergence of the observer for increasing values of k [simulation results].

Figure 6 . 4 :

 64 Figure 6.4: Without the PE assumption asymptotic convergence of the observer is lost, a bias remains [simulation results].

Figure 6 . 5 :

 65 Figure 6.5: Vector measurements with additive noise [simulation data].

Figure 6

 6 Figure 6.6: k = 1. Observer performance under noisy measurement [simulation results].

Figure 7 . 1 :

 71 Figure 7.1: Interconnection of systems (7.5)-(7.6).

Figure 7

 7 Figure 7.2: k = 4, γ 1 = 1, γ 2 = 0.2, ω max = 255 [ • /s]. Convergence of the observer: angular velocity error [simulation results].

Figure 7 Figure 7

 77 Figure 7.3: k = 4, γ 1 = 1, γ 2 = 0.2, ω max = 255 [ • /s]. Convergence of the observer: torque (solid) and estimated torque (dashed) [simulation results].

Figure 7

 7 Figure 7.5: k = 4, γ 1 = 1.5, γ 2 = 1. Left: Observer (7.1) tracks slowly varying torques with phase shift. Right: Extended observer (7.9) (γ 3 = 0.15) reduces the residual error [simulation results].

Figure 7 . 6 :

 76 Figure 7.6: Estimation of sublimation torques for a Halley-like comet [simulation results].

Figure 7

 7 Figure 7.7: k = 3, γ 1 = 1, γ 2 = 0.2. Estimation of P without noise (left) and with 10 -7 noise level (right) [simulation results].

Figure 8 . 1 :

 81 Figure 8.1: Interconnection of systems (8.7)-(8.8).

ω

  Figure 8.2: k = 5, γ = 1, γ 2 = 0.8. The faster the trajectory, the faster the convergence of the ratios of inertia estimation [simulation results].

Figure 8

 8 Figure 8.3: Above a certain value, the parameter k has little impact on the convergence [simulation results].

  Figure B.1: Left: function φ(•). Right: elliptic functions.

Proposition B. 3 .

 3 For m < 1, denote 2 cn and sn are T -periodic. dn is T 2 -periodic.Proof. We havet (φ + π) = sin 2 ζ = t(φ) + 2

0 2

 0 r(t)| u(s)ds, ∀t ≥ 0 Proof. Consider ρ ≥ 0 fixed. We haveẊ = M (ρ)X + (M -M (ρ)) XThus, for any t ≥ 0X( t ) = e M (ρ) t X 0 + t 0 e M (ρ)( t -s) (M (s) -M (ρ)) X(s)ds |X( t )| ≤ e M (ρ) t |X 0 | + t 0 e M (ρ)( t -s) M (s) -M (ρ) |X(s)|ds ≤ Ke -c t |X 0 | + t 0 Ke -c( t -s) l|s -ρ| |X(s)|dswhich can be written, usingKl = 2 u(t) ≤ 1 + t |s -ρ| u( s)ds = 1 + t 0|s -ρ| u(s)dsSetting ρ = r(t) concludes the proof of the Lemma.

  (t)|u(s)ds + 1 -v(t), ∀t ≥ 0 (B.2)We now show a second lemmaLemma B.2. v satisfies v(t) ≥ 1 + t 0 |s -r(t)| v(s)ds, ∀t ≥ 0 C.3.1 Case √ d 3 |ω 1 (0)| ≤ √ d 1 |ω 3 (0)|Introduce the following parameters

  0) 2 > 0 A 3 sign(ω 3 (0)) d 3 d 1 m |A 1 | = sign(ω 3 (0)) ω 3 (0) 2 + d 3 |d 2 | ω 2 (0) 2 = 0 w -A 2 A 3 d 1 A 1 = -sign(ω 1 (0)ω 3 (0)) d 1 |d 2 |ω 3 (0) 2 + d 1 d 3 ω 2 (0) 2 = 0 (C.2)with the convention sign(0) = 1. Remark C.1. Inequality √ d 3 |ω 1 (0)| ≤ √ d 1 |ω 3 (0)| implies that m ∈ (0, 1] and m = 1 if and only if d 3 |ω 1 (0)| = d 1 |ω 3 (0)|

1

  see Proposition B.3 2 see Remark 4.4 for an explanation C.3.2 Case √ d 3 |ω 1 (0)| > √ d 1 |ω 3 (0)|

d 3

 3 |ω 1 (0)| = d 1 |ω 3 (0)| (C.4)the solutions are periodic and draw simple closed curves. The period depends on the initial conditions. We call these solutions Type C. Their trajectory is not contained in a plane, which is the difference with the Type B. They are also called general solutions, as in practice two moments of inertia are never equal and the condition (C.4) is satisfied for almost-all initial conditions.The case whered 3 |ω 1 (0)| = d 1 |ω 3 (0)| (C.5)corresponds to the limit case m = 1 the Jacobi elliptic function converge and we have the following limitslim ∞ ω 1 (t) = lim ∞ ω 3 (t) = 0 lim ∞ ω 2 (t) = -sign(w)A 2 = -sign (ω 1 (0)ω 3 (0)) |d 2 | d 1 ω 1 (0) 2 + ω 2 (0) 2

  Figure D.2: A Sun sensor (TSL250 from Texas Instrument).

Figure D. 4 :

 4 Figure D.4: Curve obtained from the raw measurements for a complete rotation of the wheel [experimental results].

N

  , where N is a positive integer, and number of measurement N m for each step. Steps for any integer -N ≤ j ≤ N : • set the value ψ = j∆ψ • record N m measurements of the sensor output • store the mean of the N m measurements m[j]. This allows to cancel a large part of measurement noiseThe values m[j] are then interpolated to produce a continuous phase response.

  Figure D.6: Angular response of the four sensors without calibration [experimental results].

Figure D. 7 :

 7 Figure D.7: Output of the four cells during a rotation with ψ = 0. Raw and calibrated data [experimental data].

Figure D. 10 :

 10 Figure D.10: Sun sensor TSL250 mounted on a bicycle wheel.

Figure

  Figure D.11: The wheel is mounted on a gear connected to a trolley on a motorized rail which allows to control its phase.

Figure

  Figure D.17: Gyros measurements (solid) compared with observer based estimation (dashed) [experimental results].

  Cette thèse traite du problème d'estimation de vitesse angulaire d'un corps rigide à partir de capteurs de direction embarqués. Cette question très générale a des applications

dans des domaines variés d'ingénierie et de science appliquée. Démarrons par une brève description d'exemples spécifiques. i) En aérospatial, le déploiement de satellites en orbite commence par une phase dite de detumbling au cours de laquelle la vitesse angulaire est réduite à zéro par un contrôle actif. La stratégie de contrôle en boucle fermée nécessite une estimation de cette variable. Le satellite utilise ses propres capteurs embarqués pour réaliser l'estimation. ii) Les projectiles à fort spin sont très communs en balistique. Par exemple, l'arme intelligente XM25 envoie des grenades qui estiment leur rotation pour déterminer la distance parcourue (de sorte que l'explosion du projectile puisse être activée à une distance définie par un utilisateur). Les capteurs sont embarqués dans la grenade. iii) Enfin, le problème d'estimation de vitesse angulaire peut aussi être rencontré dans le domaine émergent des dispositifs intelligents pour le sport et le divertissement. Dans de nombreux sports, il est important pour les athlètes de développer leur capacité à faire tourner une balle. L'estimation peut être faite par des capteurs embarqués, sans avoir recours à un dispositif externe.

Table 1

 1 est la vitesse angulaire du corps rigide exprimée dans le repère engin.Les gyromètres mesurent ω à un biais et un bruit de mesure près. Les capteurs de direction fonctionnent de la manière suivante. Considérons un vecteur de référence a exprimé dans le repère inertiel, par exemple le champ magnétique terrestre. L'expression de a dans le repère engin à l'instant t est Il s'agit de la variable produite par un capteur de direction comme un magnétomètre (d'autres exemples sont présentés en Table1.3). Selon l'application, le vecteur de référence a peut lui-même être connu ou non. Par exemple, l'expression du champ magnétique terrestre dans un repère géocentrique peut être calculé avec un haut degré de précision à partir du modèle IGRF (International Geomagnetic Reference Field). Ce modéle est construit à partir de mesures expérimentales et actualisé tous les cinq ans. De nombreux algorithmes ont été proposés ces dernières années pour estimer l'attitude à partir de mesures de gyromètres biaisées et de mesures de direction. La fusion des capteurs peut être faite par un algorithme de type filtre de Kalman [CBIO06, SRK + 08] ou bien par un observateur non-linéaire [MHP08, MS10, VSO08, TRB11, SMS11, GFJS12, TMHL12].

	a(t) = R(t) a	(1.9)
	vecteur	capteur	unité
	champ magnétique	magnétomètre	[G,T]
	direction d'une source lumineuse capteurs solaires	-
	champ de gravité	accéléromètre	[m.s -2 ]
	direction de la Terre	Earth sensor	-

L'évolution de ω est régie par les célèbres équations d'Euler

[START_REF] Landau | Mechanics. MIR Moscou[END_REF]

ω = J -1 (Jω × ω + τ ) (1.8)

où J est la matrice d'inertie du corps rigide et τ est le couple qui lui est appliqué. Lorsqu'il est important de connaître R et ω à un certain degré de précision, le corps rigide est généralement équipé de capteurs inertiels embarqués comme des gyromètres et des capteurs de direction.

La quantité a(t) est une mesure de direction (ou mesure de vecteur). .3: Exemples de mesures de direction et capteurs associés Il est connu que de tels capteurs suffisent à estimer R et ω.

Table 1 .

 1 4: Classification de la littérature et contributions de la thèse comprise. Ils ont été testés en profondeur et ont gagné la confiance de la communauté scientifique. Toutefois, les preuves de leur convergence sont rares et, quand elles existent, sont seulement locales. De plus, leur implémentation est considérée comme relativement gourmande en temps de calcul 8 . Les observateurs non-linéaires sont bien plus simples à implémenter et se prêtent à une analyse de convergence précise. La classification que nous proposons est résumée dans la Table 1.4. Nos contributions sont des méthodes dérivatives indirectes et des méthodes d'estimation directe.

Part I Analysis of the information contained in single vector measurements Chapter 3 Single-axis rotation Axe de rotation fixe. Dans

  ce chapitre, nous considérons une rotation à axe fixe, paramétrée par un angle. Nous construisons un estimateur statique de cet angle en le reliant à la variation d'indice de la courbe des mesures par rapport à une origine bien choisie. Il en résulte un estimateur très simple pour lequel on dispose d'une formule

	r)	closed disk of center z ∈ C (or equivalently R 2 ) and radius r > 0
	Ω c	complementary set of a subset Ω ⊂ C
	γ • γ	concatenation of curves γ and	γ
	N	set of non-negative integers
	Z	set of relative integers
	S/	quotient of set S by the equivalence relation
	| • |	Euclidean norm in any dimension, or modulus for a complex number
	M	matrix norm induced by the Euclidean norm, i.e. M = max |X|=1 |M X|
	|f | ∞	infinite norm of a function f i.e. max t |f (t)|
	J ac f	Jacobian matrix of a function f
	cn, sn, dn Jacobi elliptic functions
	O(x)	Landau's notation i.e. O(x) |x| is bounded
	V (t, X)	generally designates a real function V of time t ∈ R and state X ∈ R n
	∂V ∂t (t, X) derivative of V with respect to t ∈ R ∇V (t, X) gradient of V with respect to X ∈ R n
	•	angular degree
	rad	angular radian
		end of the proof of a Theorem or a Proposition
		end of the poof of a Lemma
	2.6 Acronyms
	EKF extended Kalman filter
	IGRF International Geomagnetic Reference Field
	LAM long-axis mode
	LTI	linear time-invariant
	LTV	linear time-varying
	PE	persistent excitation
	SAM short-axis mode
	SNR signal-to-noise ratio
	UCO uniformly completely observable

d'erreur explicite. La simplicité des équations sous-jacentes nous permet de donner une expression analytique de l'impact des défauts de capteurs sur l'estimateur.

Table 3 .

 3 1: Standard deviation of the error estimate for various SNR and sampling rates.

  the corresponding invertible matrix and λ = -max( (λ 1 ), (λ 2 )). We have, for all Y ∈ R 2 ,

	e	√	kA 2 t Y = P	e	√ 0 kλ 1 t	e	0 √ kλ 2 t	P -1 Y
	from which we deduce
			|e

  definite negative. Therefore, by choosing sufficiently large k, the first term in (7.8) is made dominant over the other terms that are not scaled by √ k. Under these conditions,

	d dt V is definite negative and system (7.2) is uniformly exponentially stable.
	Theorem 7.1: Solution to Problem 7.1

  since P is bounded. Thus, P is bounded and closed in R 2 , it is therefore compact. If P is defined by a set of ordered vertices x 1 , . . . , x n counter-clockwise oriented, then the parameters u i and g i are determined by the following formulas

	• P is obviously non-empty since, e.g. for any x ∈ P, X =	x 0	∈ P
	• for any X =	x r	in P , x ∈ P is bounded and, since the circle C of radius r and
	center x is included in P, we have
			r =	1 2	max x,y∈C	|x -y| ≤	1 2	max x,y∈P	|x -y| < ∞
	Remark B.1.
	angular rate [rad/s]							True 1 Estimate 1 True 2 Estimate 2 True 3 Estimate 3
						t [s]
	Figure 9.2: Observer estimates using only direction sensors versus true angular velocity
	measured by the embedded gyros [experimental results].

  (t) e ct |X( t )| K|X 0 | , Kl ln 2-c)t |X 0 | ⇔ e ct |X(t)| ≤ Ke √ ln 2t |X 0 | ⇔ u( t) ≤ v( t)Thus, it is enough to prove that for all t ≥ 0,

				√
				Kl
		v(t) e	√	ln 2t
	We have
	|X(t)| ≤ Ke ( √
		w(t) u(t) -v(t) ≤ 0
	Define r(t) as	r(t) max(t -	√ ln 2, 0)

Indeed, the barrel twist relates spin rate to range independent of muzzle velocity

Or one of its many possible parametrizations, see Appendix A

Or an equivalent equation if a parametrization of R was used

Ou l'une de ces nombreuses paramétrisations, voir Appendice A

Ou une équation équivalente si R est donnée par une autre paramétrisation

This is true only if and only if the reference frame R i is inertial

In the case of a satellite e.g., the torque could be generated by inertia wheels, magnetorquers, gravity gradient, among other possibilities

We thank S. Changey from ISL for providing the dataset

Remerciements

Nicolas fut un directeur disponible, ouvert et bienveillant. Outre ses conseils ines-

PE assumption in free-rotation

The PE Assumption 6.4 is the cornerstone of the proof of the main result. It is interesting to investigate whether it is often satisfied in practice (we have seen in Remark 6.1 that it might fail). In this section we consider that the angular rate ω satisfies the free Euler's equations (2.13). We will prove that Assumption 6.4, or equivalently condition (6.2), is satisfied for almost all initial conditions.

As proved and illustrated in Appendix C, there are four types of solutions to the free Euler's equations.

Type A ω is constant, which is observed if and only if ω(t 0 ) is an eigenvector of J.

Type A' J 1 > J 2 > J 3 singular case: ω 1 (t) and ω 3 (t) vanish, ω 2 (t) tends to a constant when t goes to infinity. This situation is observed only for a zero-measure set of initial conditions ω(t 0 ).

Type B the trajectory is periodic and draws a non-zero diameter circle. This situation is observed if and only if two moments of inertia are equal and ω(t 0 ) is not an eigenvector of J.

Type C J 1 > J 2 > J 3 regular case: the trajectory is periodic and not contained in a plane. This situation is observed for almost all initial conditions ω(t 0 ).

Study of Type A and Type A' solutions

Assume that ω(t 0 ) = 0 (or simply ω) is an eigenvector of J, namely for i = 1, 2 or 3

Proposition 6.1. For all t, R(t) can be written

Proof. R(t) and r u (wt)R 0 have the same value R(t 0 ) for t = t 0 . Moreover,

Thus both functions satisfy the rotation differential equation (2.6), which concludes the proof by Cauchy-Lipschitz uniqueness theorem.

Assumption 8.4 (persistent excitation

Observer design and analysis of convergence

Observer definition

We extend the observer (5.2) under the form

where k, α, γ 1 , γ 2 > 0 are constant tuning parameters. The inputs of this observer are the two vector measurements a and b; the outputs are ŵ and d, which are estimates of the angular velocity ω and the ratios of inertia d.

Error dynamics

To study the convergence of this observer, we introduce the scaled errors For any choice of γ 1 , γ 2 > 0 and for 0 < α < 2 √ 1 -p, the observer (8.3) defines an error (8.4) which, for sufficiently large k > 0, converges locally uniformly exponentially to 0.

Remark 8.3. Since the analysis builds on a persistent excitation argument, it is not possible to obtain arbitrary fast convergence. This is a structural difference with observer (7.1).

Further extensions

Relaxing the assumption that the torque is known, it is possible to connect the ratio of inertia observer (Y system (8.8)) with the torque and angular rate estimation observer (7.1), yielding the super extended form

which produces estimates ω, χ, d of the angular rate, external torque and ratios of inertia. As shown in Section 7.2, the rate of convergence of the torque observer error can be made arbitrary fast by sufficiently large k. This is enough to show the local convergence of the super-extended observer (8.10) (the analysis is similar to Section 8.2.3).

Simulation results

In this section we illustrate the convergence of the observer and sketch the sensitivity with respect to the tuning gain k and the maximal angular rate ω max .

Simulations were run for a model of a CubeSat [START_REF]The CubeSat program[END_REF]. The rotating rigid body under consideration is a rectangular parallelepiped of dimensions about 20 × 20 × 10 [cm 3 ] and mass 2 [kg] assumed to be slightly non-homogeneously distributed. The resulting moments of inertia are

No torque is applied on this system, which is thus in free-rotation. The dynamics are simulated using RK4 with sampling period 0.01 [s]. As it is represented in Figures 8.2 and 8.3, the estimation error d -d convergence rates highly depends on the value ω max , thus on the trajectory. The parameter k is of importance since it guarantees convergence of the estimator. Yet, arbitrary fast convergence is not achieved. This is an illustration of Remark 8.3.

Appendices Appendix A

Parametrization of a rotation matrix

In this appendix we explain in details the various parametrizations of a rotation matrix used in the document. We mention when a parametrization is redundant (i.e non-injective), in which case we indicate the appropriate quotient to make it one-to-one.

A.1 Axis and angle

For any unit vector u ∈ S 2 and angle ζ ∈ R, the rotation matrix of axis positively led by u and angle ψ is denoted by r u (ψ). Its exact expression is

where [u × ] is the matrix of the cross-product by u. Namely

This skew-symmetric matrix satisfies the following property

This parametrization is redundant as shown by the following equalities

It is easy to see that is an equivalence relation. The following mapping is one-to-one

where / designates the set quotient by relation .

A.2 Euler angles

The parametrization by a set of three Euler angles is useful because they usually bear physical meaning handy for practitioners. There are 12 different sets of Euler angles. In this document we use the 'ZXZ' nomenclature as presented in [START_REF] Landau | Mechanics. MIR Moscou[END_REF], §35. The rotation matrix can be written

Chebyshev center of a convex polygonal

Let P be a non-empty convex polygonal in the R 2 plane. In this appendix we consider the following problem.

Problem B.1

Find a circle included in P of maximum radius.

The center of such a circle is called Chebyshev center of P. In the following we will see that this is actually a linear programming problem, which has at least one solution. For each side i of P, we denote by u i the outward pointing normal unit vector. There exists g 1 , . . . , g n ∈ R such that P = 1≤i≤n

x, u i x ≤ g i Consider now a circle C of center x and radius r ≥ 0. Namely

C is included in P if and only if for all i max v∈S 2 u i (x + rv) ≤ g i From Cauchy-Schwarz inequality, the maximum value of

on the closed set P defined by the linear constraints

This highlights the linear programming nature of this problem. To show that it has at least one solution, it is enough to show, since f is a continuous function, that the constraint set P is a non-empty compact set.

Proof As W (0) = 0, we have for all t ≥ 0, W (t) ≤ 0 and w(t) ≤ t W (t) ≤ 0. This concludes the proof.

Remark B.3.

For K ≤ 2 one can show, with many more calculations, that the better following inequality holds:

Analytical solutions to the free-rotation Euler's equations

In this appendix we derive an analytical solution to the Euler's equation with no torques (2.13). The analysis depends on conditions bearing on the moments of inertia J 1 , J 2 , J 3 and the initial conditions ω(0). Some solutions are easy to compute (Types A and B). For the other ones (Types A' and C), we use the results on Jacobi elliptic functions presented in Appendix B.2 Without loss of generality we consider that

The free-rotation Euler's equations can be written

Below, we distinguish four types of solutions, which are characterized by the moments of inertia of the rigid body and the initial conditions. For brevity, and without loss of generality, we assume that the intial time is t 0 = 0.

C.1 Constant solutions: Type A

The simplest case one can imagine is when ω(0) is an eigenvector of J. In this case, the solution is constant

This happens exclusively in the following cases.

• if J 1 > J 2 > J 3 , when at least two coordinates of ω(0) are zero

The constant solutions are called Type A.

C.2 Circle solutions: Type B

We assume that exactly two moments of inertia are equal, e.g.

which yields for the ratios of inertia

Thus ω 3 is constant and ω 1 , ω 2 are easily given for all t by

Hence, the solutions draw circles in planes orthogonal to the vertical body reference vector b 3 . The case J 1 > J 2 = J 3 is analog. The solutions draw circles in planes orthogonal to the first body reference vector b 1 . These circle solutions are called Type B.

C.3 General and pathological solutions: Types C and A'.

Appendix D

Experimental results

D.1 Sun vector measurements

Below, we investigate how coarse Sun sensors (photocells) can be used to provide vector measurements and the kind of non-linearities that one has to deal with in practice. The estimator studied in Chapter 3 is employed. a designates the direction of a source of light (e.g. the Sun in the case of a satellite). Experimental results obtained on a simple testbed are provided.

D.1.1 6 Sun sensors for one vector measurement

Consider a photocell attached to the rigid body and n ∈ S 2 its outward-pointing normal. The energy deposited in the photocell is proportional to the cosine of the angle of incidence of the Sun. Thus, the output signal (voltage or current) is, roughly speaking, a cosine

D.1.3 Pictures of the experimental testbed

The experimental testbed we designed and used to obtain the presented results is pictured in the Figures D.10

D.2 Testing the angular-rate observer on a smartphone

We test the angular rate observer (5.2) on data produced by sensors available on a smart-phone.

D.2.1 Sensor framework

The Android Nexus 5 is equipped with

• tri-axis gyros and accelerometers MPU-6500

• tri-axis magnetometers AK8963

These sensors are identified in

D.2.2 Observer performance

We test the angular rate observer (5.2) on data produced by the aforementioned sensors. The reference vectors are the gravity field and the magnetic field. The direction sensors are thus the accelerometers and magnetometers. The experiment is very simple. The • the rotation motion is fully damped after only one period, so the excitation is anything but persistent.

D.2.3 Perspectives

We are currently working on a more sophisticated experiments, such as the one presented in Chapter 9.