|
DYNAMIQUE D'UN SATELLITE ET QUATERNION D'ATTITUDE (suite) |
|
|
|
Mises à jour : septembre 2004, sept 2011, avril 2013 Cours de mathématiques appliquées à la mécanique I Dérivation et rotation instantanée II Calcul de la matrice de passage III Equation de comportement d'un quaternion IV Equations du mouvement d'un satellite autour de son centre d'inertie V Initialisation du quaternion d'attitude |
|
|
|
|
|
|
|
|
Voir un article spécialement
consacré au produit des quaternions suivant la base de représentation
BIBLIOGRAPHIE : Vous trouverez une présentation et surtout
les exemples d'utilisation des quaternions dans : Techniques et technologies
des véhicules spatiaux ( CNES )
Dans
le cours précédent, nous avons associé à toute rotation géométrique (
déplacement ) un quaternion Q, d'un espace S muni d'une algèbre.
Les
développements qui suivent se proposent de relier la rotation instantanée à la
dérivée d'un quaternion.
I DERIVATION ET ROTATION INSTANTANEE W.
Nous
désignerons par I, J, K les quaternions de base de l'espace de référence absolu
et ceux i, j, k d'une base R mobile ( liée au satellite, par exemple , mis ce
n'est pas une nécessité ici ). Nous rappelons que la rotation
instantanée ne concerne que le mouvement de R par rapport à Ra.
Le
quaternion Q représente une rotation axiale, par exemple en mécanique
classique, la rotation géométrique permettant
de passer à l'instant t de la base I, J, K absolue à la base relative i, j, k.
Si la base relative est mobile, la rotation est évolutive, tant par son angle
que par son axe de rotation.
Ce
quaternion, variable à priori, est donc fonction du temps t.
1°) LA DERIVATION TEMPORELLE D'UN QUATERNION Q(t) :
Cette
dérivation est définie comme pour les vecteurs, par exemple par la relation :
En
pratique la dérivée d'un quaternion, tout comme le quaternion lui-même n'a pas
de lien avec la base, puisqu'elle caractérise l'évolution d'une transformation
géométrique intrinsèque. C'est une notion
intrinsèque. On veillera à ne pas confondre avec la dérivée d'un vecteur
qui, elle est associée à une base de référence, tout comme la vitesse.
|
Le lecteur vérifiera que |
|
est un quaternion pur en utilisant les règles opératoires sur les quaternions |
2°) INTRODUCTION DE LA ROTATION INSTANTANEE W .
Nous
allons montrer sur un exemple simple, puis de manière générale que le
quaternion
purement
imaginaire ou pur représente le vecteur rotation instantanée du repère R
( ou du solide si la base i, j, k est attachée à ce solide), exprimé dans la base
d'écriture du quaternion
a) Sur un exemple:
Soit le quaternion représentant la rotation géométrique
d'angle q autour de l'axe fixe K
et son dérivé, défini comme plus haut
Le lecteur démontrera aisément que le vecteur rotation instantanée du mouvement
absolu, exprimé dans la base I J K est :
b) Généralisation :
En
utilisant les quaternions Q(t) et son conjugué, représentant les rotation RQ
et son inverse, rotations variables de R3, donc fonction du temps t.
[
i(t), j(t), k(t) ] est la base mobile déduite de la base
canonique de R3 [ I, J, K ] inertielle , par la
rotation géométrique associée à Q(t). Nous notons WS (représentation en axes relatifs )
et WI (représentation en axes inertiels )
les quaternions représentant la rotation géométrique W dans chaque
base.
Ainsi
en exprimant i comme image de I, puis en le dérivant , naturellement dans la
base inertielle, il vient:
|
mais |
|
est
pur, donc posant : |
|
Nous
obtenons ainsi la dérivée absolue de i par :
Nous
retrouvons la définition classique du vecteur rotation de la mécanique et
découvrons que ce vecteur est donné en axes absolus par:
En
utilisant ( voir cours précédent) les propriétés du quaternion inverse, il
vient une relation importante :
où
le vecteur rotation ( au sens mécanique de vitesse angulaire ) est exprimé, en
tant que quaternion de S, dans la base de référence absolue I J K des
quaternions de S et ensuite comme vecteur sur la base absolue X, Y, Z de R3.
Rappelons
que la rotation instantanée est une notion intrinsèque ne dépendant que du
mouvement et pas de la base d'expression du vecteur. Seule l'expression de ce
vecteur dépend de la base de décomposition.
c) REMARQUE CAPITALE ET NON EVIDENTE :
NB : Vous ferez bien la différence entre
le quaternion et le vecteur associé. Rappelons aussi que la rotation est une
donnée intrinsèque qui ne dépend que du mouvement d'une base par rapport à
l'autre. Ce sont les expressions des quaternions ou des vecteurs qui changent
avec la base.
Comprenons bien que dans l'espace affine R3, le vecteur
rotation est une notion intrinsèque. Donc si le quaternion Q transforme la base
des quaternions I J K ( ou des vecteurs X Y Z ) en la base des quaternions i j
k ( ou des vecteurs x y z ), alors le vecteur Wconsidéré comme lié à I J K, devient
W'et n'a donc pas changé de
coordonnées
Ainsi pour retrouver le vecteur W intrinsèque, il faut effectuer la rotation inverse de
quaternion conjugué de Q.
Donc
exprimé dans la base mobile, le vecteur rotation instantanée, en tant que
quaternion, a un quaternion WS qui se
déduit de WI par la rotation inverse de Q donc
par la rotation de quaternion conjugué de Q.
où
le vecteur rotation est exprimé, en tant que quaternion, dans la base de
référence relative i, j, k des quaternions de S et ensuite comme vecteur sur la
base relative x, y, z de R3.
Donc
comme plus haut, mais dans la base relative on a :
En
utilisant ( voir cours précédent) les propriétés du quaternion inverse, il
vient une relation importante symétrique de celle obtenue plus haut:
CONCLUSION: Le même vecteur rotation s'exprime dans chaque
base comme ci-dessous :
|
|
Rotation absolue exprimée en axes inertiels I, J, K ou X, Y, Z |
|
|
Rotation absolue exprimée en axes mobiles i, j, k ou x, y, z liés au satellite |
II CALCUL DE LA MATRICE DE
PASSAGE :
Appelons
P( XYZ --> x y z ) la matrice de passage de la base de référence absolue à la base mobile
relative, liée au satellite. P est aussi appelée matrice de l'opérateur
ou matrice de la rotation qui transforme XYZ en xyz.
ATTENTION :
Nous savons que les anciennes composantes s'expriment à l'aide de P et
des nouvelles composantes. Soit
Soit
Q un quaternion unité élément de S, donc associé à une rotation RQ
|
où
traduit en quaternions, appliqué à V |
|
Pour
calculer la matrice de changement de base on écrit:
faisons
le calcul pour la première relation, à partir du transformé vectoriel de X qui
donne le vecteur x
donc
tous calculs terminés, dans la base absolue :
|
Les
éléments de P sont notés Pij et aij ceux de P-1
avec |
|
NB : Vous pourrez vérifier que l'écriture de P est la même que ce
soit avec Q exprimé dans I J K ou Q exprimé dans i j k.
II bis
CALCUL DU QUATERNION CONNAISSANT P :
Donnons
une procédure de calcul qui pourrait être utile dans un codage informatique:
Remarque initiale : La norme du quaternion étant égale à 1, au moins un des q²i >= 1/4
|
On pose et on calcule : |
|
||
|
D'après la remarque, il existe au moins un Si positif ou nul. On commencera donc par rechercher le plus grand des qi ou encore le Si le plus grand. On ne gardera que le signe +, le signe étant indifférent. |
|||
|
Si |
|
|
|
|
Si |
|
|
|
|
Si |
|
|
|
|
Si |
|
|
|
III EQUATION DE COMPORTEMENT DU
QUATERNION EN AXES SATELLITE:
A
l'expérience, les équations en axes satellite sont les plus importantes dans
les applications pratiques, puisqu'en mécanique, il apparaît clairement que les
axes satellites sont évidemment mieux placés par rapport au corps que ceux du
référentiel fixe..
L'évolution
temporelle du quaternion est donnée par l'équation :
donnant
ainsi une relation fondamentale matricielle d'évolution du quaternion
exprimé en axes satellite avec une matrice Ms qui se calcule avec les
composantes de WS :
avec
une matrice Ms qui vaut :
IV EQUATIONS DU MOUVEMENT D'UN
SATELLITE AUTOUR DE SON CENTRE D'INERTIE:
La
résolution du problème de l'attitude d'un corps dans l'espace consiste donc à
résoudre en parallèle 7 équations totalement couplées :
·
Une
équation vectorielle équivalente à 3 équations différentielles scalaires,
résultant du théorème du moment cinétique
appliqué au centre d'inertie du solide et fournissant p, q, r, donc
la rotation instantanée absolue.
·
Une équation vérifiée par le quaternion Q,
équivalente à 4 équations différentielles scalaires, traduisant les équations
d'évolution du quaternion d'attitude, nécessitant la connaissance de p, q, r
calculés en parallèle. Son intégration fournit donc le quaternion donnant
l'attitude du satellite.
|
|
Théorème
du moment cinétique, donnant p, q, r grâce aux conditions initiales sur W. CONDITIONS: repère x, y, z lié au satellite et moment des forces en G, connu sur les axes satellite |
|
|
|
|
|
Equation
d'évolution du quaternion d'attitude Q et conditions initiales permettant le
calcul de Q. La connaissance de Q permet le calcul de la matrice de passage de la base fixe à celle du satellite. |
|
|
|
|
|
Equation de redondance, permettant la vérification du bon fonctionnement de l'algorithme de calcul |
La
dernière relation est une équation redondante à vérifier et surtout destinée à
rappeler qu'il faut veiller à normer le quaternion d'attitude tout au long du
calcul.
NB
: Suivant les bons conseils de spécialistes de l'utilisation des quaternions, il est recommandé de normer le quaternion Q tout au long
des calculs d'intégration lorsque nécessaire.
De
plus, il vaut mieux travailler en axes mobiles qu'en axes fixes, comme souligné
plus haut ( simple respect des symétries du solide ou de ses axes principaux )
·
L'attitude
du satellite en référentiel absolu ( c'est à dire les angles de repérage, Euler
ou autres ), s'en déduit alors par l'intermédiaire de la matrice de passage P, elle même fonction de Q.
·
Si l'on s'intéresse à un mouvement rapporté à un référentiel non
absolu, il suffit d'effectuer les changements de repère, concernant les
positions et les vitesses ( attention à la composition des vecteurs rotation ).
NB : On peut aussi évaluer le vecteur rotation relatif du
satellite par rapport à ce repère R, et calculer alors la matrice MS,
qui permet de déterminer le quaternion d'attitude Q du satellite par rapport à
R.
Par exemple, si le repère relatif R a une rotation instantanée
connue dans le repère absolu par ses composantes (WX, WY, WZ), alors il faudra
remplacer p, q, r par p*, q*, r* dans le calcul de Ms, puisque c'est l'attitude
dans le repère relatif R qui est alors recherchée:
QUATERNION ASSOCIE A UNE ROTATION
AXIALE UNIFORME :
On
peut être amené à suivre l'évolution d'un repère tournant à vitesse
constante autour d'un axe fixe D.
On appelle w la vitesse angulaire D l'axe et p q r les composantes du vecteur
rotation instantanée W constant.
V INITIALISATION DU QUATERNION D'ATTITUDE :
1°) GENERALITES SUR L'INITIALISATION : Si, comme il est préconisé, le
quaternion d'attitude Q est inclus dans le système différentiel du mouvement,
alors il doit être initialisé par Q(t=to)=Qo
L'initialisation
doit donc prendre en compte une position initiale So du repère de
référence lié au satellite.
Soit en
se donnant Qo qui transforme le repère de référence Ro en So , le
problème est réglé
Soit en
se donnant les éléments de la rotation Ro ----> So avec l'angle
et l'axe, ce qui revient au cas précédent
Soit en
se donnant la matrice de passage (Pij) de Ro ----> So CAS 1
Soit en
se donnant les angles d'Euler CAS 2 au temps to
Soit en
se donnant les angles de Cardan CAS 3 au temps to
être
faite à l'aide des angles de repérage et donc de la matrice P de passage qui
s'exprime en fonction de ces angles. Nous supposons donc connue la matrice P
par ses coefficients Pij.
Le
lecteur vérifiera, en observant attentivement l'écriture de P et notamment la
forme des coefficients hors diagonale, que l'on peut calculer le quaternion
d'attitude avec la matrice P, par les relations suivantes:
Une
série de tests est nécessaire sur les éléments diagonaux pour calculer le
quaternion à une indétermination près de signe. Dans tous les cas on supposera
cependant que qo est dans l'intervalle 0, 1.
En
pratique dès que le signe d'une composante de Q est choisi, il détermine
automatiquement le signe des autres, mais la matrice P de passage n'en est pas
affectée. Il n'y a donc pas lieu de s'inquiéter, et si qo est différent de 0
nous prendrons toujours cette composante > 0 sinon ce sera q1 etc....
Exemple de programmes de calcul d'initialisation du quaternion Q = [q0
q1 q2 q3 ]
REMARQUE INITIALE :
Il est capital que la matrice de passage P=(Pij) soit
orthogonale, si des imprécisions de calcul laissent présager un doute, il faut
alors rechercher la matrice orthogonale P* la plus proche de P, au sens de la
norme de Frobenius ( || A ||= S(aij)² )
L'algorithme de recherche par itération est le suivant :
a) PROGRAMME D'INSPIRATION PERSONNELLE : Voici un exemple de programme sous Matlab. VOIR OU FAIRE UNE APPLICATION NUMERIQUE
NB
: Le dernier indice 0 des variables à 2 indices indique la valeur initiale, par
exemple qt30 est la valeur initiale de la troisième composante q3 du quaternion
d'attitude.
NB
: Pmat est la matrice de passage de la base absolue à la relative.
NB
: le symbole ^ est celui de la puissance.
Un
problème se pose, celui du calcul de la racine carrée de nombres calculés
numériquement et pouvant donner à cause des arrondis, des valeurs voisines de
0, donc au signe douteux. Des tests doivent être prévus pour pallier cette
difficulté.
Nous
prendrons toujours q0 >0 ou nul
trac=trace(Pmat);
if (1+trac)>0,
qt00=(1+trac)^0.5/2;
else
qt00=0;
end
if
(1-trac+2*Pmat(1,1))>0,
qt10=(1-trac+2*Pmat(1,1))^0.5/2;
else
qt10=0;
end
if
(1-trac+2*Pmat(2,2))>0,
qt20=(1-trac+2*Pmat(2,2))^0.5/2;
else
qt20=0;
end
if (1-trac+2*Pmat(3,3))>0,
qt30=(1-trac+2*Pmat(3,3))^0.5/2;
else
qt30=0;
end
if abs(qt00)>1e-8,
qt10=qt10*sign((Pmat(3,2)-Pmat(2,3))/4/qt00);
qt20=qt20*sign((Pmat(1,3)-Pmat(3,1))/4/qt00);
qt30=qt30*sign((Pmat(2,1)-Pmat(1,2))/4/qt00);
else if
abs(qt10)>1e-8,
qt20=qt20*sign((Pmat(1,2)+Pmat(2,1))/4/qt10);
qt30=qt30*sign((Pmat(1,3)+Pmat(3,1))/4/qt10);
qt00=qt00*sign((Pmat(3,2)-Pmat(2,3))/4/qt10);
else if
abs(qt20)>1e-8,
qt10=qt10*sign((Pmat(1,2)+Pmat(2,1))/4/qt20);
qt30=qt30*sign((Pmat(2,3)+Pmat(3,2))/4/qt20);
qt00=qt00*sign((Pmat(1,3)-Pmat(3,1))/4/qt20);
else if
abs(qt30)>1e-8,
qt20=qt20*sign((Pmat(3,2)+Pmat(2,3))/4/qt30);
qt10=qt10*sign((Pmat(1,3)+Pmat(3,1))/4/qt30);
qt00=qt00*sign((Pmat(2,1)-Pmat(1,2))/4/qt30);
end
Q=[qt00 qt10 qt20
qt30];
b) PROGRAMME RECUPERE SUR LE NET ( SOURCE SÛRE) :
REMARQUE :
La logique précédente ci-dessus m'est personnelle et pourrait être moins bonne
que celle qui est fournie par le CNES dans le site http://logiciels.cnes.fr/MARMOTTES/marmottes-mathematique.pdf , notamment quand on calcule une
racine d'une quantité proche de 0. Allez donc voir ce site, pour plus de
détails ou récupérez la version PDF chez moi ou sur le site Marmottes et
lisez les pages 43et 44. ( NB : on se méfiera du choix du signe opposé au nôtre des 3
dernières composantes du quaternion, dans la définition du quaterniondit
d'orientation par rapport à celui de rotation( le nôtre ), dans le site CNES )
:
Voici la logique réécrite avec nos notations : L'idée de base étant de calculer la composante la plus grande
du quaternion, vu qu'on est sûr d'en avoir au moins une plus grande que 0.25,
car la somme des carrés des 4 est 1.
Calcul de So=P11+P22+P33
Si So > - 0.1
Alors qo=0.5*(1+So)^0.5; q1=(P32-P23)/4/qo; q2=(P31-P13)/4/qo;
q3=(P21-P12)/4/qo;
Sinon calculer S1=P11-P22-P33
Si S1 > - 0.1
Alors q1=0.5*(1+S1)^0.5; q2=(P12+P21)/4/q1; q3=(P31+P13)/4/q1;
q0=(P32-P23)/4/q1;
Sinon calculer S2= - P11+P22-P33
Si S2 > - 0.1
Alors q2=0.5*(1+S2)^0.5; q1=(P12+P21)/4/q2;
q3=(P31+P13)/4/q2;q0=(P13-P31)/4/q2;
Sinon calculer S3= -P11-P22+P33
Si S3 > - 0.1
Alors q3=0.5*(1+S3)^0.5; q0=(P21-P12)/4/q3;
q1=(P31+P13)/4/q3;q2=(P32+P23)/4/q3;
Fin si
Fin si
Fin si
Elle
est plus simple et plus facile à exprimer, donc à garder.
CAS 2 : INITIALISATION AVEC LES ANGLES D'EULER OU
CALCUL DES ANGLES D'EULER:
MON POINT DE VUE :
[1]
P11 = cosj cosy - sinj cosq siny
[2]
P21 = cosj siny + sinj cosq cosy
[3]
P31 = sinj sinq
[4]
P12 = -sinj
cosy - cosj cosq siny
[5]
P22 = -sinj
siny + cosj cosq cosy
[6] P32 = cosj sinq
[7] P13 = sinq siny
[8] P23 = -sinq cosy
[9]
P33 = cosq
On
peut alors calculer le quaternion à l'aide de a) ou b), comme ci-dessus.
Si
l'on considère qu'un angle est parfaitement défini par son sinus et son
cosinus, on peut considérer qu'il y a 6 inconnues, siny, cosy, sinj, cosj, sinq, cosq, et donc que
6 équations sont nécessaires sur les 9.
Comme
dans chaque paquet de 3 équations, les 2 premières entraînent la troisième, il
y a 3 équations surabondantes.
Issu
du site CNES http://logiciels.cnes.fr/MARMOTTES/marmottes-mathematique.pdf ou lu dans la version PDF page 44.
Ci-dessous
la représentation clasique des " angles d'Euler faisant passer du repère
XYZ à xyz, grâce aux trois rotations successives y d'axe Z, q d'axe u, j d'axe z, respectivement de quaternions associés Q1,
Q2, Q3:
Le
lecteur à titre d'exercice d'entraînement pourra établir les formules de calcul
des composantes du quaternion produit, en fonction des angles d'Euler.
Le
quaternion produit des 3 rotations est :
où
a est l'angle de la rotation
résultante et D l'axe de cette rotation, éléments
qi se calculent assez aisément en fonction des angles d'Euler, grâce au
résultat ci-dessous.
La
donnée des angles d'Euler permet donc une initialisation très rapide du
quaternion.
NB : Comparez les 2 points de vue. VOIR OU FAIRE UNE APPLICATION NUMERIQUE
CALCUL DES ANGLES D'EULER :
Le
problème se pose en pratique en sens contraire, vu que dans les applications
c'est l'attitude qui est intéressante. La fin d'un calcul avec les quaternions
donne la matrice instantanée P, avec laquelle il faut
reconstituer l'attitude, c'est à dire par exemple les angles d'évolution du
solide en mouvement, quelquefois simplement en suivant l'axe satellite en
repère inertiel.
a) Par exemple, il faut retrouver les angles d'Euler. Il s'agit donc de résoudre le système d'équations (1) à (9). ATTENTION : q est ici un des angles d'Euler celui de la nutation, y celui de précession et j de rotation propre. Voir angles d'Euler
b)
On peut aussi résoudre 3 des 4 équations à 3 inconnues q, y, j où naturellement une des équations
est surabondante (puisque le quaternion à une norme unité <---> somme des
carrés égale à 1).
c)
On peut aussi calculer ainsi y et j:
et q par :
ce qui achève le calcul.
CAS 3 : INITIALISATION AVEC LES
ANGLES D'EULER OU CALCUL DES ANGLES DE CARDAN:
La
succession de repères est :
XYZ --- Y --> abZ -- q ---> xbg --- F --> xyz
Nous
avons ainsi défini les angles conventionnels, également appelés ANGLES DE
CARDAN :
- Lacet y mesuré autour de Z
- Tangage q mesuré autour de b ( voisin de Y lorsque les
angles sont petits)
- Roulis f mesuré autour de x (
voisin de X lorsque les angles sont petits)
La
matrice P de passage de XYZ à xyz s'explicite classiquement :
a) Expression du quaternion en fonction des angles de Cardan :
Un
calcul long, dont j'assume la responsabilité ( car je ne l'ai pas retrouvé
ailleurs pour l'instant en 2011) donne un résultat d'une étonnante symétrie,
pour l'expression du quaternion à l'aide des angles de Cardan.
NB
: En 2013 une nouvelle recherche sur la toile me permet de conforter mes
calculs et en plus le cours associé est très formateur :
http://www.morris.umn.edu/academic/math/Ma4901/gravelle.pdf
Les
quaternions associés aux rotations élémentaires et à la rotation résultante
sont :
Tous
calculs effectués ( ce serait un excellent entraînement pour le lecteur), on
obtient les relations remarquables suivantes:
qui permettent une initialisation sans problème du
quaternion Q.
b) Calcul des angles
de Cardan en fonction du quaternion:
Le lecteur utilisera la
matrice de passage soit en termes d'angles de Cardan soit de quaternion pour
établir les relations suivantes:
REMARQUE D'UN LECTEUR COMPETENT (M Alban Leroyer à Nantes, en sept 2004)
"Concernant
la composition des rotations en utilisant les quaternions, notamment sur les
dernières modifications que vous avez apportées concernant les initialisations
des quaternions pour les angles d'Euler et de Cardan, je trouve qu'il est plus
simple de travailler à partir de l'expression des quaternions dans les bases
courantes (équivalent des matrices de passage d'une base à une autre). Dans ce
cas, il ne faut pas oublier d'inverser le sens du produit des quaternions pour
trouver le quaternion correspondant à la composition des rotations
successives."
Raisonnons
sur un exemple vu plus haut, où toutes les quantités ( quaternions et vecteurs
) sont exprimées en base canonique de départ Ro:
Voilà ce que dit ce
correspondant en clair : au lieu d'exprimer les composantes des axes successifs
de rotation dans la base canonique de référence de Ro, soit :
On exprime les composantes
de l'axe de rotation n° i dans la base relative R(i -1) précédente lors de la
rotation R(i -1) ---> Ri, soit :
On peut alors parler des
quaternions relatifs associés aux expressions relatives :
Le lecteur vérifiera, au
besoin en s'exerçant aux mathématiques des quaternions et des matrices de
changement de base, ou par le calcul direct, que le quaternion Q cherché
peut aussi s'écrire :
Le résultat peut bien
évidemment être généralisé au cas de n rotations successives.
Merci donc
pour cette remarque fort utile. On aura aussi bien noté que l'ordre des
quaternions dans le produit a été inversé.
Guiziou Robert et Chiavassa Rolland /1994,
révision novembre 2002, révision sept 2004, révision sept2011 et avec l'aide de
M Alban Leroyer