WikiPrépaLivrets

X ENS Modélisation PSI 2024Sujet et corrigé

3,0(2 votes)
  • Modélisation cinématique des mécanismes
  • Optimisation numérique et méthode d'Euler explicite
  • Moteur à courant continu
  • Systèmes linéaires, matrices d'état et filtre de Kalman
  • Fonctions de transfert et diagrammes de Bode
  • Stabilité et correction des systèmes asservis

Téléchargements

  • Rapport du jury : non disponible

Présentation du sujet

Modélisation et amélioration des performances de vol d'un drone bio-inspiré du colibri, projet Colibri
Afficher ou masquer la section

Le sujet modélise un drone miniature à battement d'ailes inspiré du colibri. La première partie porte sur la conception mécanique de la transmission moteur-ailes, l'optimisation du rapport de réduction pour limiter l'échauffement, et l'estimation en temps réel de l'état de charge de la batterie par un filtre de Kalman implémenté sur microcontrôleur. La seconde partie étudie et corrige la stabilité du vol en lacet et en tangage par asservissement.

  1. 1Partie 1, I - Étude de la transmissionModélise le mécanisme cinématique transformant la rotation continue du moteur en battement alternatif des ailes et optimise deux paramètres géométriques par minimisation numérique d'une fonction coût.
  2. 2Partie 1, II - Étude de l'association moteur/réducteurDétermine, à partir des équations du moteur à courant continu, le rapport de transmission qui minimise la puissance perdue par effet Joule en vol stationnaire.
  3. 3Partie 1, III - Étude de la batterieConstruit un modèle d'état de la batterie puis un filtre de Kalman permettant d'estimer en temps réel son état de charge, avec une implémentation en Python sur microcontrôleur.
  4. 4Partie 2, I - Stabilisation du vol en lacetÉtablit la fonction de transfert du mouvement de rotation en lacet, détermine le coefficient de frottement fluide à partir d'un essai en oscillation libre, puis règle un correcteur proportionnel intégral pour respecter les critères de précision et de stabilité.
  5. 5Partie 2, II - Stabilisation du vol en tangageÉtablit la fonction de transfert du mouvement de tangage à partir des équations de la dynamique, étudie sa stabilité selon les paramètres de conception, puis évalue une correction proportionnelle dérivée.

L'épreuve en chiffres

Moyenne 8,9 / 20 · écart-type 3,6 · 1 007 présents · où vous situez-vous ?
Afficher ou masquer la section
Moyenne
8,9/ 20
Écart-type
3,6
Présents
1 007
Durée
5 h
moyenne 8,905101520
Deux tiers des copies environ (moyenne ± écart-type)

Votre note sur 20 à ce sujet, en conditions de concours.

Source : document officiel du concours, épreuve du 19 avril 2024. Notes publiées par le concours (après harmonisation le cas échéant). Courbe : estimation par une loi normale.

Ces sujets peuvent vous intéresser

Lecture du sujet en ligne

L'énoncé complet, avec les formules et les figures, sans ouvrir le PDF.
Afficher ou masquer la section

ECOLES NORMALES SUPERIEURES ECOLE POLYTECHNIQUE

CONCOURS D'ADMISSION 2024

VENDREDI 19 AVRIL 2024 08h00-13h00

FILIERE PSI - Epreuve n ^∘8

MODELISATION (XUSR)

Le projet COLIBRI

Figure 1 : Colibri en vol stationnaire
Un des moteurs de l'innovation technologique réside dans le biomimétisme. La nature ayant mis des millénaires à adapter les organismes vivants afin d'optimiser leurs performances, il paraît tout naturel de s'en inspirer pour créer des objets innovants. À titre d'exemples, on peut citer un système d'éclairage public bioluminescent inspiré des micro-organismes marins, ou la construction d'une éolienne compacte dont les pales sont inspirées des ailes des libellules capables de capter des vents très faibles de l'ordre de 6 km/h.
Ce sujet s'inspire d'un projet de recherche mené à l'Université Libre de Bruxelles. Le but est de concevoir un drone miniature et léger, qui pourrait se déplacer dans des endroits isolés et accidentés de façon stable, en s'inspirant de la dynamique de vol du colibri. Cet oiseau, l'un des plus petits au monde puisque sa taille varie entre 6 et 30 cm , possède une masse d'une petite dizaine de grammes. Malgré cela, il est capable de battre des ailes jusqu'à une fréquence de 200 battements par seconde. Il possède un vol stationnaire très stable et peut se déplacer en moyenne autour de 56 km/h. Mais il est surtout connu pour être le seul oiseau à pouvoir voler en arrière grâce à un mouvement particulier de ces ailes.
De nombreux laboratoires de recherche développent des robots qui se déplacent par battements d'ailes. La figure 2 classe quatre projets différents (AVFL Robotic Hummingbird, Nano Hummingbird, KU-Beetle et Colibri) par rapport aux espèces vivantes de colibris, pour 3 paramètres importants de la conception (longueur des ailes, masse, fréquence de battements).
Figure 2 : Répartition des projets
Question 1. Comparer les différentes technologies de robots, par rapport aux espèces de colibris vivants. Quel paramètre doit encore être optimisé ?
Le sujet se compose de deux parties principales. Dans la première, vous allez être amenés à comprendre les choix qui ont été faits sur la conception du drone nommé Colibri. On cherchera à maintenir en équilibre le drone. Les ailes auront donc un mouvement dans le plan horizontal ( x⃗_0, y_0^(→−) ). La seconde partie cherchera à analyser les performances du drone en mouvement de lacet et de tangage afin de les améliorer.

Partie 1 : Modélisation du mécanisme

Le mouvement des ailes du Colibri est imposé par un motoréducteur à courant continu, associé à un système de transmission optimisé pour obtenir un mouvement le plus symétrique possible entre l'aller et le retour des ailes sur un cycle de battement et ainsi éviter les moments parasites. Dans une première partie, on étudiera le mécanisme de transmission. On montrera ensuite que le rapport de réduction du réducteur influe sur l'échauffement du moteur. On cherchera donc à déterminer un rapport optimal. Enfin, l'étude du mécanisme se terminera par une modélisation en temps réel de la batterie afin de connaître le plus précisément possible le temps de vol restant du drone Colibri.

I. Étude de la transmission

Au fur et à mesure de l'évolution du projet de recherche, le système de transmission de puissance du moteur aux ailes a beaucoup évolué. On cherche notamment à optimiser certains paramètres dimensionnels de la transmission afin d'obtenir le mouvement le plus sinusoïdal possible des ailes sur un cycle de battement.
Figure 3.a : Système réel
Figure 3.b : Schéma cinématique de la partie droite
Figure 3 : Système de transmission
La dernière conception prend la forme de la figure 3 . Un premier étage d'amplification à la sortie du motoréducteur est constitué du bâti et des solides S_1, S_2 et S_3. Il permet de modifier la rotation continue en sortie du motoréducteur modélisée par l'angle θ en rotation alternative du solide S_3 modélisée par l'angle φ. Un second étage est constitué d'un système poulie/courroie permettant la rotation des ailes modélisée par l'angle ψ.
On note :
  • (x⃗_0, y⃗_0, z⃗_0) la base associée au solide S_0, corps du Colibri, telle que (OC^–, y⃗_0^–) = γ et OC = L
  • (x⃗_1, y⃗_1, z⃗_1) la base associée au solide S_1 telle que (x⃗_0, x⃗_1) = θ et OA = L_1
  • (x⃗_2, y⃗_2, z⃗_2) la base associée au solide S_2 telle que (x⃗_1, x⃗_2) = δ et AB = L_2
  • (x⃗_3, y⃗_3, z⃗_3) la base associée au solide S_3 telle que (x⃗_0, x⃗_3) = φ et BC = μ. On note D_3 = 17 mm le diamètre de la roue 3 .
  • (x⃗_4, y⃗_4, z⃗_4) la base associée au solide S_4 telle que (x⃗_0, x⃗_4) = ψ. On note D_4 = 5, 2 mm le diamètre de la roue 4.
L'encombrement étant un critère primordial pour ce projet, on se laisse libre les deux paramètres L_2 et γ. Le premier paramètre influence la trajectoire des points des ailes lors
d'un battement de ces dernières, tandis que le second influence l'amplitude de battement. On pose OB^(→−) = h(t)y⃗_0.
On désire à ce stade une amplitude du premier étage de 20^∘ et une fréquence de battement de 20 Hz des ailes; on prendra L = 14, 75 mm et L_1 = 2, 25 mm. Pour des questions d'encombrement, on choisit les deux paramètres à optimiser dans des intervalles tels que : L_2 ∈ [13, 17] et γ ∈ [22^∘, 26^∘] (soit Lcosγ ∈ [13, 25; 13, 68] et Lsinγ ∈ [5, 52; 6, 47] ).
Question 2. Réaliser sur votre copie deux schémas du premier étage d'amplification correspondant aux angles θ = θ_1 = − π/2 et θ = θ_2 = π/2.
Question 3. Dans le cas général, et à partir de deux fermetures géométriques, montrer que :
L_1^2 cos^2 θ + (L_1 sinθ − Lcosγ + Lsinγtanφ)^2 = L_2^2
Question 4. Donner dans chacun des cas de la question 2 :
  • le signe de φ et l'expression de h .
  • l'expression de tanφ en fonction de L, L_1, L_2, et γ.
Question 5. En déduire l'expression de tanφ en fonction de L, L_1, L_2, γ et θ.
Question 6. Déterminer la relation liant ψ, θ et les paramètres constants du problème.
La détermination des paramètres L_2, et γ est réalisée en minimisant la fonction coût suivante :
J(L_2, γ) = ∫_0^(2π)[(φ − φ_m) + Asinθ]^2 dθ
avec A = 20^∘ et φ_m = (φ_1 + φ_2)/2 où φ_1 = φ(θ_1) et φ_2 = φ(θ_2).
Question 7. Justifier la forme de la fonction coût.
Pour la suite, on suppose que le module numpy a été importé (avec la commande import numpy as np), et on utilisera un schéma d'intégration basé sur la méthode des trapèzes. Les variables globales suivantes seront supposées déjà définies :
  • thetac : tableau numpy de flottants de 200 valeurs régulièrement réparties entre 0^∘ et 360^∘, représentant l'angle θ;
  • paramL2 : tableau numpy de flottants de 50 valeurs entre 13 et 17 inclus représentant la longueur L_2 à optimiser.
  • paramGamma : tableau numpy de flottants de 50 valeurs entre 22^∘ et 26^∘ inclus représentant l'angle γ à optimiser.
Question 8. Écrire les instructions en Python permettant de générer un tableau numpy nommé TabJ donnant les valeurs de la fonction J pour tout couple ( L_2, γ ) pris dans les listes paramL2 et paramGamma. Les constantes L_1 et L seront supposées déjà définies respectivement par les variables L1 et L.
Le tracé du résultat donne la figure suivante :
Figure 4 : Simulation numérique de la fonction coût J(L_2, γ)
Le minimum de cette fonction coût permet de choisir L_2 = 15 mm et γ = 23, 75^∘. On peut alors tracer l'angle φ − φ_m en fonction de θ pour un cycle de battement sur la figure 5.
Pour optimiser les problèmes de dynamique, il est nécessaire d'avoir une accélération angulaire la plus symétrique possible (entre l'aller et le retour de l'aile) sur un cycle de battement.
Figure 5 : Évolution de l'angle φ − φ_m durant un cycle de battement
Pour la fin de cette partie, on suppose connus :
  • phic : tableau numpy de flottants de même taille que thetac contenant les valeurs de φ − φ_m pour chaque valeur de θ.
  • D_3 = 17 mm et D_4 = 5, 2 mm seront considérées comme variables globales.
Question 9. Écrire l'(es) instruction(s) Python nécessaire(s) pour construire la liste psic des valeurs de ψ − ψ_m représentant l'angle des ailes recentré ( ψ_m.correspond à l'angle des ailes pour un angle φ_m de la roue 3 ).
Le cycle de battement est considéré stationnaire à f = 20 Hz.
Question 10. On appelle a = ψ¨ l'accélération angulaire des ailes. Montrer que :
a = 4π^2 f^2(d^2(ψ − ψ_m))/(dθ^2)
On cherche à déterminer une approximation numérique, notée TabA, de l'accélération angulaire a = ψ¨.
Question 11. En utilisant le schéma d'Euler explicite, déterminer la relation de récurrence à deux termes liant les éléments de TabA avec les éléments de psic.
Question 12. En déduire les instructions Python permettant de calculer l'accélération angulaire durant un cycle de battement à une fréquence de 20 Hz .
On obtient alors la courbe présentée en annexe 1. De plus, la précédente conception générait une accélération angulaire dont la forme est donnée aussi en annexe 1.
Question 13. Comparer les graphes d'accélération angulaire du solide lié aux ailes du drone Colibri, et conclure sur l'intérêt de cette nouvelle cinématique.

II. Étude de l'association moteur/réducteur

Les équations du moteur à courant continu sont données ci-dessous :
{C_m = C_(em) − C_p; C_(em) = k.i; C_p = k.i_0; u_m = kω_m + R.i
où:
  • C_(em) est le couple électromoteur du moteur à courant continu ;
  • C_p est le couple de perte constant dans le moteur ( i_0 = 0, 1 A est le courant associé, donné par le constructeur) ;
  • C_m est le couple en sortie du moteur ;
  • i est le courant consommé par le moteur ;
  • ω_m et u_m respectivement la vitesse de rotation du moteur et sa tension d'alimentation ;
  • k = 7, 87 ⋅ 10^(− 4) N ⋅ m/A et R = 1, 28Ω, respectivement la constante de couple et la résistance de l'induit (rotor).
La principale limitation dans la durée de vol de ce genre de drone est l'échauffement prématuré du moteur. On cherche à déterminer le rapport de transmission r_(transm) = ω_m/θ˙ optimal permettant de limiter cet échauffement tout en assurant une fréquence de rotation correcte des ailes.
On se place dans le cadre d'un essai de vol stationnaire où le drone doit supporter son poids comme charge. La masse du drone sera prise à 24 g . Une campagne d'essais a permis de déterminer la fréquence de battement des ailes à 20 Hz pour une puissance mécanique nécessaire en sortie du moteur de P^∘ = 4 W.
Question 14. Donner l'expression de la puissance électrique consommée P_e et de la puissance mécanique fournie P_m en fonction de u_m, ω_m et des constantes du moteur.
Question 15. En déduire l'expression de la puissance perdue par effet Joule en fonction de u_m, ω_m et des constantes du moteur.
On cherche à déterminer graphiquement l'ensemble des points de fonctionnement possibles du moteur dans un plan ( u_m, ω_m ) .
Question 16. Déterminer l'expression de u_m en fonction de ω_m pour un fonctionnement à vide du moteur.
Question 17. Déterminer l'expression de u_m en fonction de ω_m pour un fonctionnement à P_m maximale.
On se fixe une tension maximale U_(max) aux bornes du moteur imposée par la batterie à 6 V .
Question 18. Représenter sur votre copie la surface, dans le plan ( ω_m, u_m ), qui représente l'ensemble des points de fonctionnement possibles du moteur.
Dans notre application, la tension u_m et la vitesse de rotation du moteur ω_m ne sont pas indépendantes mais doivent respecter la condition de l'essai définie par P_m(u_m, ω_m) = P^∗.
Question 19. Donner l'équation de u_m = f(ω_m) correspondant à la condition P_m = P^∗.
La courbe représentative de cette fonction non linéaire coupe la représentation graphique de la condition issue de la question 17 pour la vitesse ω = ω_(mini) et la courbe de la tension maximale pour ω = ω_(maxi). On cherche à déterminer ces 2 valeurs en utilisant l'algorithme de dichotomie.
Question 20. Écrire une fonction dichotomie qui prend en argument une fonction monotone f et qui renvoie la solution de l'équation f(x) = 0 à 10^(− 3) près entre les valeurs 0 et 7000rad/s.
Question 21. Écrire les fonctions fmini et fmaxi arguments de la fonction dichotomie qui vont permettre de déterminer respectivement ω_(mini) et ω_(maxi). Toutes les constantes du problème k, R, i0, Petoile, Umax représentant respectivement
k, R, i_0, P^∗, U_(max) sont supposées être déjà initialisées en mémoire comme variables globales.
Question 22. Écrire les lignes de code permettant de déterminer ω_(mini) et ω_(maxi).
L'utilisation des fonctions précédentes permet de déterminer les deux vitesses extrémales admissibles ω_(mini) et ω_(maxi) pour l'application du drone. On trouve :
ω_(min1) = 2875rad/s et ω_(max1) = 6014rad/s
La recherche du meilleur rapport de réduction revient à déterminer la vitesse qui minimise la puissance perdue sachant que P_m(u_m, ω_m) = P^∗.
Question 23. Montrer que la puissance perdue peut alors se réécrire avec cette contrainte sous la forme :
P_j = R((P^∗)/(k ⋅ ω_m) + i_0)^2
Question 24. En déduire la vitesse optimale et le rapport de transmission r_(transm) ( r_(transm) > 1 ) associé. On prendra π ≈ 3, 1 pour les calculs.
Question 25. Suite à cette étude, le laboratoire de recherche qui conçoit le drone Colibri a réalisé un prototype avec un réducteur à roues dentées de rapport de transmission de r_(transm) = 39. Quelle(s) justification(s) pouvez-vous avancer pour le choix de ce rapport?

III. Étude de la batterie

Les chercheurs à l'origine du projet ont fait le choix d'une batterie Lithium-Polymère associée au motoréducteur à courant continu étudié dans la partie précédente. La gestion de la vitesse de battement d'ailes est assurée par une carte électronique positionnée sur le «nez» du robot. On donne sur la figure 6 le positionnement de cette technologie de batterie.
Figure 6 : Technologies de batterie dans le plan Puissance/Énergie spécifique
Question 26. En quoi cette technologie est-elle intéressante pour cette application?
Malgré toutes les optimisations possibles dans la structure du mini drone, le point déterminant pour le temps de vol du Colibri est l'autonomie de la batterie. Cela passe par une connaissance assez fine de l'état de la batterie en temps réel. La quantité de charge d'une batterie est définie principalement par la variable SOC («State Of Charge», pour état de charge) exprimée en pourcentage :
sOC = (Q_m − ∫idt)/(Q_m) × 100
où Q_m est la quantité de charge maximale et i l'intensité, imposée par la charge, que doit délivrer la batterie. Sous cette forme, l'état de charge est difficilement évaluable en temps réel, car Q_m évolue en fonction des charges/décharges de la batterie.
On se base alors sur le modèle de batterie de la figure 7 :
Figure 7 : Modèle de la batterie
  • C_s : capacité de la double couche de la batterie ;
  • R_t : résistance de transfert ;
  • R_1 : résistance interne de la batterie ;
  • V_0 : tension à vide (noté souvent «OCV» en anglais pour «Open Circuit Voltage »);
  • C_0 : capacité de la batterie ;
  • V_b : tension de sortie de la batterie ;
  • i : courant délivré par la batterie .
Avec le modèle de la figure 7, des essais ont montré une très bonne corrélation entre la fonction V_0 = f(SOC) et un modèle linéaire autour d'un point de fonctionnement. On posera donc pour la suite :
V_0 = a_0 ⋅ SOC
avec a_0 ∈ ℝ.

A) Modèle dynamique de la batterie

On pose X le vecteur dit d'état, représentant l'état du système dynamique batterie, tel que X^⊤ = (SOC, V_s).
Question 27. Montrer, à partir des lois de l'électrocinétique, que X vérifie l'équation:
X˙ = A ⋅ X + B ⋅ i
où A ∈ M_(2 × 2)(ℝ) et B ∈ M_(2 × 1)(ℝ). Vous donnerez les expressions de A et B.
Question 28. À partir de la loi des mailles, donner une relation liant V_b, V_s, SOC et le courant i.
Comme le courant i est imposé par la charge mécanique en sortie du moteur, il est supposé connu quelque soit le temps. On pose alors V = V_b + R_i i.
Question 29. Montrer que V = C ⋅ X, où C ∈ M_(1 × 2)(ℝ). Vous donnerez l'expression de C La connaissance de l'état de charge de la batterie n'est pas accessible directement. Mais la mesure de V_b et la connaissance de i (donc la mesure de V ) permet indirectement une mesure de la variation de X à tout instant. On est donc amené à résoudre le système d'équation suivant :
{x˙ = A ⋅ X + B ⋅ i; V = C ⋅ X
La première équation traduit la dynamique d'évolution du vecteur d'état en fonction de la commande i supposée connue (le courant est imposé par la charge mécanique). La seconde équation traduit la connaissance que l'on a du système à tout instant, à travers la mesure de V_b.
La fréquence d'échantillonnage de la mesure de V_b est notée f_e, elle est supposée élevée.
Question 30. Discrétiser l'équation (1) du système précédent par la méthode d'Euler explicite et donner les matrices A_e et B_e telles que :
X_(k + 1) = A_e ⋅ X_k + B_e ⋅ i_k
où X_k représente l'état du système au temps t_k.
Finalement, on cherche à résoudre le système discrétisé suivant :
{X_(k + 1) = A_e ⋅ X_k + B_e ⋅ i_k, (3a); V_(k + 1) = C ⋅ X_(k + 1), (4a)

B) Filtre de Kalman

Il est classique dans ce genre d'application d'utiliser le filtre de Kalman pour résoudre ce problème : avoir une estimation de l'état de charge de la batterie à tout instant. La force de cet algorithme est de pouvoir prendre en compte les incertitudes de mesures et de modèles.
Dans toute la suite de cette partie, on appellera E(u) l'espérance mathématique de la variable aléatoire u.
On modifie donc le système d'équation (3a)-(4a) pour le réécrire sous la forme :
{X_(k + 1) = A_e ⋅ X_k + B_e ⋅ i_k + Φ_k; V_(k + 1) = C ⋅ X_(k + 1) + Ψ_(k + 1)
  • Φ_k ∈ M_(2 × 1)(ℝ) est un vecteur de deux variables aléatoires représentant les incertitudes de modèle sur les variables SOC et V_s. Ces variables sont censées être gaussiennes, stationnaires (indépendantes du temps), centrées et indépendantes. On peut donc écrire : E(Φ_k) = 0∀k ∈ ℕ et E(Φ_k ⋅ Φ_k^⊤) = Q une matrice diagonale dont les termes représentent la variance des variables SOC et V_s. Q est appelée matrice de covariance.
  • Ψ_k est une variable aléatoire gaussienne, stationnaire, centrée et de variance σ^2. On peut donc écrire : E(Ψ_k) = 0∀k ∈ ℕ et E(Ψ_k^2) = σ^2. Elle représente les erreurs de mesure de la tension de sortie de la batterie V_b, et est indépendante de Φ_k.
L'algorithme de Kalman consiste à construire un estimateur de l'état X_k = (SOC_k, V_(sk))^⊤ du système à l'instant t_k. Cet estimateur noté X^_k = (SOिC _k V^_(sk))^⊤ va dépendre de l'estimateur de l'état et de la commande à l'instant précédent, et de la mesure à l'instant t_k : X^_(k + 1) = F(X^_k, i_k, V_(k + 1)). On choisit un estimateur linéaire de telle sorte que l'on peut écrire:
X^_(k + 1) = A_f ⋅ X^_k + B_f ⋅ ⋅ _k + K_(k + 1) ⋅ V_(k + 1)
La suite de cette partie consiste à déterminer les matrices A_f, B_f et K_(k + 1).
On définit l'erreur d'estimateur par le vecteur :
ε_(k + 1) = X^_(k + 1) − X_(k + 1)
Un estimateur parfait imposerait ε_k = 0∀k ∈ ℕ. Mais cela est impossible vu le caractère aléatoire de la mesure; de plus l'état X_k n'est jamais accessible, on ne connaît que la valeur de V_k. On cherche donc un estimateur optimal, au sens où on veut annuler la moyenne de l'erreur et minimiser sa variance. Dans la suite on pourra noter l_d la matrice identité.
Question 31. Exprimer ε_(k + 1) en fonction de ε_k, X^_k, i_k, Φ_k et Ψ_(k + 1). En déduire, en justifiant, l'expression de E(ε_(k + 1)) en fonction de E(ε_k), X^_k et i_k.
Question 32. Choisir l'expression de A_f et B_f permettant d'obtenir une équation récursive de la moyenne de l'erreur d'estimation sous la forme: E(ε_(k + 1)) = M ⋅ A_e ⋅ E(ε_k). Vous donnerez l'expression de la matrice M.
Question 33. Reprendre alors l'expression de l'estimateur et montrer que X^_(k + 1) vérifie:
X^_(k + 1) = A_e ⋅ X^_k + B_e ⋅ i_k + K_(k + 1)[V_(k + 1) − C(A_e ⋅ X^_k + B_e i_k)]
K_(k + 1) est appelé le gain du filtre de Kalman, et apparaît comme une pondération entre la fidélité du modèle et la fidélité de la mesure. La suite va consister à déterminer K_(k + 1) en travaillant sur la dispersion de l'erreur d'estimation.
on note P_k = E(ε_k ⋅ ε_k^⊤).
Question 34. Montrer que P_k est une matrice symétrique.
Question 35. Calculer P_(k + 1), et montrer que son expression peut se mettre sous la forme :
P_(k + 1) =, (A_e ⋅ P_k ⋅ A_e^⊤ + Q) − K_(k + 1) ⋅ C ⋅ (A_e ⋅ P_k ⋅ A_e^⊤ + Q) − (A_e ⋅ P_k ⋅ A_e^⊤ + Q) ⋅ C^⊤ ⋅ K_(k + 1)^⊤; + K_(k + 1)[C(A_e ⋅ P_k ⋅ A_e^⊤ + Q)C^⊤ + σ^2]K_(k + 1)^⊤
Par la suite, et afin d'alléger les notations, on pose P_(k + 1)^((k)) = A_e ⋅ P_k ⋅ A_e^⊤ + Q. Cette matrice s'interprète comme la prédiction au temps t_(k + 1) de la matrice de covariance de l'erreur d'estimation.
On cherche l'expression du gain de Kalman K_(k + 1) qui minimise la somme des variances de l'erreur d'estimation sur les variables SOC et V_s du vecteur d'état. Cette condition revient à vérifier :
(∂tr(P_(k + 1)))/(∂K_(k + 1)) = 0
Cette notation de dérivée s'interprète comme les dérivées partielles du scalaire tr(P_(k + 1)) par rapport aux termes de la matrice K_(k + 1), rangés dans une matrice symétrique de même taille que P_(k + 1). Il est possible de démontrer les propriétés suivantes :
∀A ∈ M_(n × n)(ℝ) symétrique, B ∈ M_(n × n)(ℝ) et C ∈ M_(n × n)(ℝ) symétrique :
(∂tr(AB))/(∂B) = A et (∂tr(ACA))/(∂A) = 2 ⋅ A ⋅ C
Question 36. Montrer que si A est symétrique: (∂tr(AB))/(∂A) = B^⊤
Question 37. En utilisant les propriétés précédentes, montrer que K_(k + 1) s'écrit sous la forme :
K_(k + 1) = P_(k + 1)^((k)) ⋅ C^⊤(C ⋅ P_(k + 1)^((k)) ⋅ C^⊤ + σ^2)^(− 1)
Question 38. Montrer alors que :
P_(k + 1) = (I_d − K_(k + 1) ⋅ C) ⋅ P_(k + 1)^((k))
L'algorithme de Kalman est donc un algorithme récursif qui peut se résumer de la façon suivante :
  1. Phase d'initialisation:
Choix de X^_0, P_0
2. À chaque pas de temps :
a) Prédiction
X^_(k + 1)^((k)) = A_e ⋅ X^_k + B_e ⋅ i_k; P_(k + 1)^((k)) = A_e ⋅ P_k ⋅ A_e^⊤ + Q
b) Mise à jour, correction à partir des mesures :
K_(k + 1) = P_(k + 1)^((k)) ⋅ C^⊤ ⋅ (C ⋅ P_(k + 1)^((k)) ⋅ C^⊤ + σ^2)^(− 1); P_(k + 1) = (I_d − K_(k + 1) ⋅ C) ⋅ P_(k + 1)^((k)); X^_(k + 1) = X^_(k + 1)^((k)) + K_(k + 1)[V_(k + 1) − C ⋅ X^_(k + 1)^((k))]

C) Implémentation

Question 39. On cherche à implanter l'algorithme de Kalman sur un microcontrôleur embarqué. Le problème de ce genre de carte est le peu de mémoire et de puissance de calcul. Faire une étude succincte de complexité à chaque pas de temps de l'algorithme de Kalman appliqué à l'estimation de l'état de charge de la batterie, et conclure sur la faisabilité d'implantation dans un microcontrôleur embarqué.
On envisage d'implémenter l'algorithme précédent sur un microcontrôleur embarquant un langage MicroPython. Les fonctionnalités de ce langage sont identiques à Python, mais les matrices seront traitées comme des listes de listes, la bibliothèque numpy n'étant pas installée nativement.
Figure 8 : Exemple d'illustration de carte intégrant MicroPython
Un module Bluetooth sera ajouté pour communiquer les résultats de l'algorithme en temps réel. Cette partie du projet n'est pas étudiée ici.
Dans la suite, on considère que l'on peut connaître à tout instant, grâce à deux capteurs, le courant i délivré par la batterie, ainsi que la tension à ses bornes V_b avec les fonctions MesureCourant() et MesureTension(). Ces fonctions ne prennent pas d'argument et renvoient chacune un flottant. De plus, on importe le module time dont la méthode time() renvoie la valeur du temps courant sous la forme d'un flottant. Enfin, quelque soit le résultat trouvé à la question 30 , on pourra considérer que les matrices A_e et B_e s'écrivent sous la forme :
A_e = (1, 0; 0, a_(11)) B_e = ((b_0)/(b_1))
Les valeurs a_(11), b_0 et b_1 dépendent de l'incrément de temps entre deux mesures et des constantes du problème. On suppose pour la suite que l'on a défini 3 fonctions nommées respectivement a11, b0, b1 qui prennent en argument un incrément de temps et qui renvoient la valeur du coefficient correspondant.
On présente ci-dessous le code implanté dans la carte :
1## Importation du module et valeurs numériques##
2 import time
3 C0,Cs=11252.4,27.31 # (F)
4 Ri,Rt=0.117,0.164 #(Ohm)
5 a,b=1.076,4.75 #(V)
6 ## Main ##
7 Initialisation()
8tps=time.time()
9 while True :
10 ik=MesureCourant()
11 vk=MesureTension()
12 dt=time.time()-tps
13 X,P=Prediction(X,P)
14 X,P=MiseAJour(X,P)
15 tps=time.time()
16 #Cette ligne permet de transmettre le résultat de l'algorithme mais n'est pas traité ici.
17 time.sleep(0.05) #permet d'attendre la nouvelle mesure.
Question 40. Écrire la fonction Initialisation() qui initialise toutes les variables du problème sous la forme de variables globales. On initialisera les variables X^_0^⊤ = (1, 0) et P_0 = (0, 0; 0, 0), et on s'impose : σ = 0, 01 et Q = (0, 03^2, 0; 0, 0, 03^2).
Dans la suite des questions, on pourra avantageusement utiliser la forme particulière de la matrice A_e.
Question 41. Écrire la fonction Prediction( X, P ) qui renvoie une estimation de X^_(k + 1) et P_(k + 1) à partir de la phase de prédiction de l'algorithme.
Question 42. Écrire la fonction MiseAJour( X, P ) qui renvoie les valeurs de X^_(k + 1) et met à jour la valeur de la matrice de covariance de l'erreur P_(k + 1).
Question 43. On rappelle que l'objectif de cette partie était de pouvoir évaluer l'état de charge de la batterie au cours du temps. Comment obtenir cette information à partir du problème ci-dessus ?

Partie 2 : Amélioration des performances de vol

Les parties précédentes ont permis de proposer un modèle pour le fonctionnement de l'oiseau. L'objectif de cette partie est d'améliorer les performances de vol du prototype du projet de recherche. Les performances attendues pour le vol sont les suivantes :
Figure 9: Prototype du projet Colibri et définition des mouvements
Performance Niveau
Précision en vitesse de lacet Erreur nulle en échelon de vitesse
Stabilité en vitesse de lacet Marge de phase de 60^∘
Précision en angle de tangage Moins de 10^∘ d'erreur de position angulaire
La définition des différents mouvements est donnée figure 9 .

I. Stabilisation du vol en lacet

On considère la stabilisation de vol en lacet (en rotation autour de ( G, z→ )). On appelle I_z l'inertie de l'oiseau autour de ( G, z→ ) et ω_z sa vitesse angulaire. En vol, lorsqu'il tourne sur lui-même, l'oiseau est soumis à un couple de frottement fluide − D_z ω_z. Il dispose de barres de stabilisation qui peuvent créer, si besoin, un couple C_(za)z⃗ avec le souffle de l'air. Ce couple n'est pas instantanément créé quand on le souhaite. L'évolution temporelle de C_(za), quand on souhaite créer un échelon de couple de stabilisation C_z, est représentée sur la figure 10.
Figure 10: Evolution du couple C_(za) quand un échelon de couple C_z est demandé au servomoteur qui pilote les barres de stabilisation. L'échelle des ordonnées a été normalisée à 1
Question 44. Déterminer la valeur numérique de la constante de temps du modèle caractéristique reliant C_(za) à C_z si on considère un modèle du premier ordre: C_(za) = 1/(1 + T_z p)C_z.
Question 45. Déterminer les expressions analytiques de la matrice [A] et du vecteur colonne [B] telles que la dynamique du système puisse s'écrire [ω˙_z; C˙_(za)] = [A][ω_z; C_(za)] + [B]C_z Vous donnerez ces expressions en fonction de D_z, I_z et T_z.
Question 46. Déterminer l'expression analytique de la fonction de transfert Ω_z(p)/C_z(p).
Afin de déterminer la valeur numérique de D_z, l'oiseau est attaché à une barre verticale, qui oscille en torsion. La raideur en torsion de la barre est K_z. On mesure l'angle de lacet φ_z(t) de l'oiseau au cours du temps (φ˙_z(t) = ω_z(t)). Cette mesure est représentée sur la figure 11.
Figure 11: Evolution de l'angle de lacet pour une oscillation en torsion libre
Question 47. Déterminer l'expression analytique de la valeur (ln(φ_(z, i)^M))/(ln(φ_(z, i + 1)^M)) où φ_(z, i)^M est la valeur d'un maximum de la courbe φ_z(t), en fonction de D_z, I_z et K_z. Cette expression permet de calculer la valeur de D_z.
Question 48. Indiquer si la fonction de transfert caractéristique du vol de l'oiseau Ω_z(p)/C_z(p), autour de l'axe de lacet, est stable ou non. Justifier votre réponse.
Au cours d'un essai en vol, on observe une dérive dans la vitesse de rotation. Pour la supprimer, on décide d'intégrer un correcteur PI dans la boucle d'asservissement de ω_z, comme indiqué sur la figure 12.
Figure 12 : Asservissement en vitesse de l'angle de lacet du drone
Question 49. Indiquer si le critère de précision en vitesse de lacet est satisfait ou non.
Les diagrammes de Bode en boucle ouverte du système corrigé sont représentés sur la figure 13.
Figure 13: Diagrammes de Bode de la boucle ouverte pour l'asservissement en vitesse de rotation sur l'angle de lacet de l'oiseau. Les valeurs choisies pour le tracé sont K_i = 1 N.m.s et T_i = 0, 15 s
Question 50. Déterminer la valeur de K_i pour satisfaire le critère de stabilité en vitesse de lacet.

II. Stabilisation du vol en tangage

on considère maintenant la stabilisation du vol en tangage (en rotation autour de ( G, y⃗ )). À l'équilibre, la portance F_L créée par le battement des ailes compense parfaitement le poids mg . Néanmoins, si une perturbation d'angle θ apparaît soudainement, F_L tourne légèrement, et la composante horizontale F_L sinθ ≃ F_L θ induit un mouvement horizontal, à la vitesse u. Par conséquent, la vitesse w de l'aile par rapport à l'air (due au battement de l'aile) va augmenter à la valeur w + u quand elle descend, et va diminuer à la valeur w-u quand elle monte.
Figure 14: Modélisation du mouvement de lacet
Question 51. On suppose que la force de traînée induite par le mouvement d'une aile est proportionnelle au carré de sa vitesse absolue, quand elle bat dans l'air. Quand elle descend, elle est orientée selon − x⃗_B. Quand elle monte, elle est orientée selon x⃗_B. Montrer alors que la force de traînée F_D qui s'applique à l'oiseau, sur un cycle de battement d'aile rapide, peut être approchée par l'expression F_D = − Ku ( K est un coefficient de proportionnalité) quand une perturbation d'angle θ apparaît.
Si on prend en compte la vitesse angulaire de tangage de l'oiseau, la force de traînée qui lui est appliquée devient alors F_D = − Ku − Kz_d θ˙.
Question 52. Expliquer l'origine physique de ce nouveau terme qui apparaît dans F_D.
Question 53. Déterminer l'équation de résultante dynamique appliqué à l'oiseau, en projection sur X_B^(→−), linéarisée à l'ordre 1.
L'équation de moment dynamique appliqué à l'oiseau, en projection sur ( G, y⃗ ), s'écrit quant à elle: I_z θ¨ = − M_u u − M_θ θ˙ + C_y où I_z est l'inertie de l'oiseau autour de ( G, y⃗ ), M_u et M_θ deux constantes de proportionnalité modélisant les couples de frottement visqueux, et C_y le couple créé par les barres de stabilisation.
Question 54. Démontrer que la fonction de transfert peut s'écrire sous la forme
(θ(p))/(C_y(p)) = (mp + K)/(a_3 p^3 + a_2 p^2 + a_1 p + a_0)
Exprimer les constantes a_0, a_1, a_2 et a_3 en fonction de m, l_2, K, M_u, M_θ, z_d et g.
Les pôles de la fonction de transfert sont représentés sur la figure 15. Deux cas sont considérés : z_d < 0 et z_d > 0.

Figure 15: Position, dans le plan complexe, des pôles de la fonction de transfert du vol de l'oiseau en rotation en tangage
Question 55. Indiquer si le vol peut être stable, dans une certaine configuration de fabrication de l'oiseau.
Afin d'augmenter les marges de stabilité pour le vol de l'oiseau, on réalise un contrôle du mouvement de tangage par un correcteur proportionnel dérivé. Le résultat de la correction obtenue est représenté sur la figure 16.
Figure 16: À gauche : position, dans le plan complexe, des pôles de la fonction de transfert en boucle fermée du vol de l'oiseau (pour le tangage), avec correction proportionnelle dérivée. À droite, mesure lors d'un essai de l'évolution de l'angle de tangage
Question 56. Conclure sur la capacité de la commande de l'oiseau à vérifier le critère de précision en angle de tangage du cahier des charges.

partie 3 : Conclusion

Question 57. Rappeler les enjeux de cette recherche. Quelles études ont été menées afin de répondre partiellement aux différents défis soulevés par ce type de conception?

Fin du sujet

Annexe 1

Courbe de l'accélération angulaire sur un cycle de battement:
Figure 17 : Accélération angulaire durant un cycle de battement
Courbe de l'accélération angulaire de l'ancienne conception:
Figure 18: Accélération angulaire de la précédente conception

Questions fréquentes

4 questions
Sur quels chapitres porte le sujet de modélisation X-ENS PSI 2024 sur le drone Colibri ?
Afficher ou masquer la section

Sur quels chapitres porte le sujet de modélisation X-ENS PSI 2024 sur le drone Colibri ?

Il porte sur la modélisation cinématique et la conception mécanique, le moteur à courant continu, le filtre de Kalman en représentation d'état, et l'automatique des systèmes asservis (fonctions de transfert, stabilité, correction).

Ce sujet demande-t-il de programmer en Python ?

Oui, plusieurs questions demandent d'écrire des instructions ou des fonctions Python, notamment pour l'optimisation numérique et l'implémentation du filtre de Kalman sur microcontrôleur.

Quelles parties sont indépendantes dans ce sujet ?

Le sujet se compose de deux parties principales : la première sur la conception mécanique, motorisation et batterie du drone, la seconde sur la stabilisation de son vol, cette dernière pouvant être abordée sans avoir traité entièrement la première.

Faut-il connaître le filtre de Kalman avant de traiter ce sujet ?

Non, l'énoncé introduit et fait démontrer pas à pas les relations de l'algorithme de Kalman à partir des équations d'état de la batterie.

Pas de description pour le moment