Projet Idopt

previous up next contents
Précédent : Grands domaines d'application Remonter : Projet IDOPT, Identification et optimisation Suivant : Actions industrielles

Résultats nouveaux

Techniques de résolution de problèmes inverses et de contrôle optimal



Participants : Anestis Antoniadis , Emmanuel Blanchard , Jacques Blum , Aïcha Bounaïm , Isabelle Charpentier , Didier Girard , Alain Le Breton , François-Xavier Le Dimet , Gérard Grégoire , Zouhir Hamrouni , Dinh Tuan Pham , Junqinq Yang


Mots-clés : analyse de sensibilité, contrôle optimal, décomposition de domaine, estimation non paramétrique, filtrage, ondelettes, problème inverse, régression, régularisation, spline, validation croisée


Résumé : Des méthodes ou outils sont développés dans le projet pour la résolution des problèmes inverses ou des problèmes de contrôle optimal. Ils se divisent en deux grandes classes : ceux qui sont déterministes et en particulier utilisent l'adjoint d'un système, et ceux qui sont stochastiques, qui vont du filtrage de Kalman à la validation croisée en passant par les méthodes d'ondelettes ou de splines. Ces méthodes seront utilisées dans les §4.2 et 4.3, associées à des applications ciblées.


Dérivation automatique d'un code adjoint; application au modèle atmosphérique Meso-NH

Les techniques d'assimilation de données fondées sur le contrôle optimal utilisent l'adjoint du modèle pour estimer le gradient de l'écart entre les solutions du modèle et les observations. Ce gradient permet de mettre en oeuvre des algorithmes d'optimisation menant à l'estimation de l'état optimal par rapport aux observations. Une difficulté essentielle provient de la dérivation de cet adjoint qui est un travail difficile et long.

Meso-NH est le modèle atmosphérique développé conjointement par le Centre National de Recherches Météorologiques et le Laboratoire d'Aérologie de Toulouse. Ce modèle, susceptible de simuler des phénomènes allant de quelques mètres (micro-échelle) à 100 km (meso-échelle), offre la possibilité de réaliser des simulations 1D vertical, 2D ou 3D en espace. À ce jour, ce code permet de simuler des événements météorologiques très fins, sans être pour autant un code de prévision.

Afin de le compléter, nous avons construit le modèle adjoint de la partie adiabatique du code (dynamique du modèle). Mené dans le cadre du projet IDOPT et avec le soutien de l'Inria (action incitative Mode Inverse Opérationnel), le développement de l'adjoint a été réalisé par utilisation du logiciel de différentiation automatique Odyssée (Inria, projet SAFIR). Cependant la taille mémoire et le temps calcul requis pour l'exécution du code adjoint, sont trop importants pour envisager une utilisation en mode opérationnel de ce code adjoint.

Induites par un recalcul local de la trajectoire du modèle direct, ces deux limitations ont été levées par la mise en oeuvre quasi automatique d'une sauvegarde de la trajectoire du modèle direct sur fichier. Les fichiers dérivés ont ensuite été optimisés manuellement; cette première étape franchie, de nombreuses applications sont envisageables : analyse de sensibilité, assimilation de données, couplages de modèle...

Ce travail présenté devant les communautés de météorologie [33] et de différentiation automatique [34], est plus largement décrit dans [46].

Etudes de sensibilité et méthodes du second ordre

Un modèle physique pour lequel une partie de l'information est manquante (coefficients, conditions initiales ou aux limites à identifier) peut être fermé par un principe variationnel minimisant, par exemple, un écart aux observations. Le problème est alors posé en utilisant la méthodologie du contrôle optimal et la solution est obtenue par résolution du système d'optimalité. De nombreuses applications en physique requièrent des études de sensibilité, qui ne sont rien d'autre que l'estimation du gradient d'une fonction réponse par rapport aux paramètres. On est donc amené, pour le calcul de sensibilité, à dériver le système d'optimalité et donc à tenir compte des propriétés au second ordre. De façon générale, on a introduit un adjoint du second ordre qui permet de résoudre ce problème. De plus ce système permet d'avoir accès à des propriétés du second ordre : spectre du hessien, produit hessien-vecteur avec des applications pour des algorithmes numériques performants de type Newton, etc. [12] [40] [41] [42]. Ces travaux font l'objet d'une collaboration avec Florida State university pour la météorologie.

Méthodes de décomposition en sous-domaines

Dans sa thèse, A. BOUNAiM poursuit un travail sur les Méthodes de Décomposition en sous-Domaines (MDD) appliquées à un problème de contrôle optimal.

L'idée motivant cette thèse est de minimiser simultanément la fonctionnelle coût et les raccordements à l'interface entre les sous-domaines en utilisant les méthodes de lagrangien augmenté. Dans le cas de la décomposition en deux sous-domaines, diverses considérations des contraintes à l'interface ont été prises en compte et des tests sur les lagrangiens augmentés associés à chaque cas ont été faits ainsi que la comparaison de leurs convergences [28].

De plus, ces algorithmes ont été comparés d'une part à une MDD avec recouvrement en utilisant la méthode alternée de Schwarz. D'autre part, des conditions de transmission aux interfaces obtenues sont proches de celles des travaux effectués par J.-D. Benamou dans le projet IDENT de l'INRIA Rocquencourt.

Des résultats numériques pour un problème-test ont été obtenus sur la machine parallèle SP1 de l'IMAG en utilisant la librairie de passage de message MPI et ont fait l'objet de deux communications [29,30]. D'autres tests ont été effectués sur le CRAY T3E du CEA de Grenoble.

Problèmes inverses, ondelettes et méthodes non-paramétriques

Dans un ensemble de travaux maintenant publiés ou en cours de publication, [10],[23], nous avons développé des méthodes d'identification de signaux noyés dans du bruit aléatoire par des méthodes de sélection de modèles et de régression non paramétrique fondées sur les notions d'analyse multirésolution et de décomposition en ondelettes. Ces travaux ont introduit de nouveaux types d'estimateur pour des classes de fonctions ``hétérogènes'' (espaces de Besov) et ont également permis de comprendre plus profondément le mécanisme d'autres méthodes d'identification plus classiques.

La plupart des résultats théoriques obtenus sont de nature asymptotique dans le sens où le nombre d'observations dont on dispose tend vers l'infini. Comme tout résultat asymptotique, il y a certains doutes sur le bienfondé des propriétés asymptotiques lorsque l'on ne dispose que d'un nombre fini et limité d'observations. La limite du nombre d'observations nécessaire pour obtenir des résultats proches des résultats asymptotiques est une question importante que nous avons abordée. Nous avons également abordé le problème d'échantillonnage dans le cadre des analyses multirésolutions permettant la réalisation d'algorithmes d'identification aussi rapides que la FFT. En effet, alors qu'en traitement du signal il est usuel de supposer que l'échantillonnage des signaux observés est à pas régulier, ceci est rarement le cas en régression où le plus souvent l'échantillonnage est aléatoire. Nous avons donc développé dans  [9] une méthode d'estimation adaptée á l'échantillonnage aléatoire.

Enfin nous avons établi les liens qui pouvaient exister entre algorithmes de régularisation pour des problèmes mal conditionnés [10], les méthodes de sélection de modèles [23] et les techniques de lissage à l'aide d'ondelettes.

Nous avons également utilisé la transformation continue en ondelettes pour détecter des discontinuités éventuellement présentes dans des signaux observés avec du bruit [22]. Ces méthodes sont actuellement poursuivies pour estimer l'intensité de certains processus ponctuels intervenant lors de la modélisation de données de survie en médecine.

Par ailleurs, pour étudier plus généralement le problème d'une fonction discontinue dont les observations sont bruitées, nous avons utilisé la méthode de régression linéaire locale. La motivation pour le recours à cette méthode non-paramétrique repose sur sa propriété d'absence d'effets de bord. Les résultats obtenus ont fait l'objet d'un séminaire [36] et d'une communication [38].

En résumé l'ensemble de ces travaux nous ont permis :

Techniques de validation croisée dans les problèmes inverses

D. Girard a poursuivi l'étude des versions randomisées de la validation croisée et d'autres critères de choix de modèles ou de paramètres de régularisation, dans plusieurs directions, notamment :

Filtrage de Kalman

  Le filtrage est l'outil de base dans l'approche séquentielle pour le problème de l'assimilation de données dans les modèles numériques, et plus particulièrement en météorologie et en océanographie. Cette approche est de type stochastique : l'état initial cherché est supposé aléatoire et par suite la dynamique du système est elle-même un processus stochastique (un terme de bruit peut être également introduit pour tenir compte de l'imperfection du modèle). Les données apparaissent comme les valeurs d'un processus lié au processus d'état contaminées par un bruit d'observation. Il s'agit alors de déterminer une bonne approximation de l'espérance conditionnelle de l'état du système au vu des données. Le caractère non linéaire des équations de dynamique dans les applications qui nous intéressent conduit à l'utilisation d'un filtre sous-optimal dit de Kalman-Bucy étendu (KBE) dans lequel on linéarise les équations au voisinage de l'estimation courante de l'état. Cependant, vu la très grande dimension de l'état du système, le filtre KBE usuel conduit a priori à des calculs prohibitifs. De plus, la grande taille de la matrice de covariance des erreurs de dynamique du modèle pose le problème de sa spécification d'une manière adéquate.

L'objet du travail amorcé dans le cadre de ce projet consiste à mener une étude approfondie des possibilités de l'approche par filtrage et à terme d'appliquer la méthode sur des données réelles. À cette fin, nous avons proposé un filtre de type Kalman étendu basé sur l'utilisation d'une matrice de covariance des erreurs singulière et de rang faible. Ce filtre opère selon le principe de ne pas faire de corrections dans les directions d'atténuation naturelle des erreurs. Les corrections sont effectuées uniquement dans des directions appartenant à un sous-espace vectoriel. Celui-ci est construit au départ par la méthode des fonctions orthogonales empiriques (EOF), mais évolue par la suite selon le modèle. Le filtre est baptisé SEEK (Singular Evolutive Extended Kalman) [6] et a été expérimenté avec succès par Gourdeau, Verron et Pham pour l'assimilation des données altimétriques dans un cadre réaliste d'un modèle aux équations primitives pour l'océan Pacifique tropical [17], [15]. De plus, une version améliorée et baptisée filtre SEIK (Singular Evolutive Interpolated Kalman) a également été proposée [43], [15]. Une étude de faisabilité basée sur un modèle quasi-géostrophique de taille réduite a montré sa bonne performance. Nous explorons également actuellement diverses autres pistes de réduction de coût de calcul pour le filtre, sans trop dégrader sa performance.

Identification et optimisation en physique



Participants : Anne Bagnérés , Emmanuel Blanchard , Jacques Blum , Stéphane Bouchereau , Isabelle Charpentier , Stéphane Despreaux , Ioana Paren , Patrick Witomski


Mots-clés : calcul parallèle, capillarité, contrôle optimal, cristallogénèse, électromagnétisme, frontière libre, fusion nucléaire, micromagnétisme, optimisation de formes, plasma, problème inverse


Résumé : Les domaines applicatifs sont la physique des plasmas (pour la fusion nucléaire), le micromagnétisme (pour les têtes de lecture magnétiques), la cristallogénèse (procédé Bridgman de fabrication de cristaux) et la capillarité. Dans chacun de ces thèmes, une optimisation est faite, que ce soit pour identifier une quantité physique non connue, pour minimiser une énergie ou pour optimiser une forme ou une frontière libre.


Physique des plasmas

Identification de sources non linéaires dans un plasma de fusion

Ce problème est motivé par l'interprétation des mesures expérimentales dans le plasma (gaz ionisé) d'un Tokamak (dispositif expérimental visant à confiner le plasma dans un champ magnétique).

L'équilibre axisymérique du plasma est régi par une EDP elliptique non-linéaire qui s'écrit :

\begin{eqnarray*}A \Psi = r f_1(\Psi) + 1/r\ f_2(\Psi) \ \ \ \ \ \ \ \ \ (1)\\ ... ... } A = - \mbox{div} ( 1/r\ \mbox{grad} .)\ \ \ \ \ \ \ \ \ \ \ \\ end{eqnarray*}

dans le plan de section méridienne du tore de coordonnées $(r,\ z)$.L'inconnue $\Psi(r,z)$ est le flux du champ magnétique; l'opérateur $A$est elliptique linéaire; le second membre de l'équation (1) représente la densité de courant du plasma. Le problème est d'identifier les fonctions $f_1(\Psi)$ et $f_2(\Psi)$qui ne peuvent être mesurées directement dans le plasma. Pour ce faire, on dispose d'informations surabondantes:

* la mesure expérimentale du flux $\Psi$et de sa dérivée normale $\partial \Psi/\partial n$ sur le bord du domaine (conditions de Dirichlet et Neumann),

* la connaissance des intégrales sur un certain nombre de cordes verticales de la composante verticale du champ magnétique, à savoir $1/r\ \partial \Psi/\partial r$.

De nombreux problèmes ouverts demeurent, comme le problème mathématique de l'identifiabilité de $f_1$ et $f_2$ à partir des mesures, et leur étude fait l'objet de la thèse d'E. Blanchard. Nous nous sommes tout d'abord intéressés au cas cylindrique où l'équation (1) devient $- \Delta \Psi = f (\Psi)$, l'identifiabilité de $f$ à partir de la condition de Neumann est un problème ouvert dans le cas général. Des études mathématiques et numériques ont été effectuées pour améliorer la compréhension de ce problème. Le problème est formulé par la minimisation de l'écart quadratique entre les mesures expérimentales et les grandeurs calculées. Une régularisation de Tikhonov est utilisée pour rendre le problème stable. La source $f$ est identifiée par décomposition dans une base de $B$-splines cubiques. Un algorithme de choix automatique du paramètre de régularisation est développé, à l'aide des techniques de validation croisée, qui s'avèrent cependant assez coûteuses dans la pratique. À l'aide d'un changement de base par rapport à une norme dérivée de la norme $H^2$, l'identification peut être réalisée de façon satisfaisante dans une base d'ondelettes à support compact de I. Daubechies, sans terme de régularisation dans la fonctionnelle à minimiser [18]. Dans le cadre d'un nouveau contrat avec le CEA, nous testons ce changement de norme de type $H^2$ pour l'identification de la densité de courant dans l'équation de Grad-Shafranov, ainsi que l'introduction de bases d'ondelettes pour que cela puisse être mis en oeuvre dans le code opérationnel sur TORE SUPRA.

Modélisation d'un jet magnétohydrodynamique en configuration axisymétrique. Applications au vent stellaire.

On considère le problème de modélisation d'un jet de plasma, applicable à la fois aux objets stellaires et aux objets extragalactiques. L'objet de la thèse d'I. Paun, co-encadrée par J. Blum et B. Michaux, est d'étudier la configuration magnétique, supposée axisymétrique, qui correspond à l'éjection du plasma à partir d'un disque de matière.

La modélisation mathématique est dérivée des équations de la magnétohydrodynamique. Elle correspond au couplage entre une équation aux dérivées partielles 2D de type Grad-Shafranov généralisée pour le flux poloïdal, décrivant la collimation de l'écoulement, et l'équation de Bernoulli généralisée, décrivant la conservation de l'énergie de l'écoulement.

La nature mathématique de ce couplage d'équations est complexe, du fait du changement de type de l'équation de Grad-Shafranov (elliptique-hyperbolique) et de la dégénérescence de l'opérateur du second ordre au voisinage de la surface d'Alfven et des surfaces magnétosoniques rapide et lente.

Par ailleurs, l'étude des conditions aux limites rendant le problème bien posé est encore un problème ouvert à l'heure actuelle.

Le travail est mené en collaboration avec l'équipe de Guy Pelletier, du laboratoire d'Astrophysique de l'Observatoire de Grenoble, qui travaille depuis plusieurs années sur la modélisation physique des jets MHD, et avec K. Morawetz du Courant Institute.

Problèmes d'optimisation en micromagnétisme

Les matériaux ferromagnétiques présentent une aimantation que l'on peut imaginer comme un champ de vecteurs de direction variable et de norme constante. Dans les matériaux que l'on étudie, elle a tendance à s'organiser en domaines, régions où le vecteur est constant, délimités par des parois, zones de transition d'épaisseur très fine. À l'équilibre, la configuration minimise une énergie, l'énergie du micromagnétisme, composée de quatre termes : l'échange, l'anisotropie, l'énergie démagnétisante et celle de Zeeman. Ce minimum est éventuellement local, ce qui se traduit par différentes configurations de paroi possibles.

Lorsqu'on applique un champ externe, il y a deux types d'évolution possibles : soit l'aimantation dans les domaines change de sens, soit les parois se déplacent. Ces phénomènes sont modélisés par l'équation de Landau-Lifshitz-Gilbert (LLG), équation parabolique non linéaire.

On s'intéresse, d'une part, au calcul numérique de la configuration de paroi et de son déplacement sous l'effet d'un champ appliqué dans un matériau à anisotropie perpendiculaire, d'autre part à la structure en domaine dans des matériaux de géométrie quelconque.

Pour ce premier travail, la structure modélisée est une plaque. L'aimantation est supposée périodique dans les deux directions du plan de la plaque. On a obtenu, à l'aide d'un code de résolution de l'équation de LLG en 3D implémenté sur Paragon, une configuration d'aimantation à l'équilibre composée de deux parois rectilignes, l'une contenant une ligne de Bloch $2\pi$,l'autre sans structure. Ce travail a fait l'objet d'un rapport de recherche [45].

La résolution du second problème a nécessité de refaire le calcul du champ démagnétisant de façon à prendre en compte des structures non périodiques. Ce travail a fait l'objet d'un stage de DEA [47]. On utilise la formulation du champ en produit de convolution que l'on calcule dans le domaine de Fourier. Les effets dûs à la périodicité sont neutralisés en ajoutant une zone où le vecteur aimantation est mis à (0,0,0) et par troncature du tenseur démagnétisant.

Des calculs de ce type ont déjà été publiés pour des matériaux de géométrie cubique ou parallélépipédique. On a généralisé la technique au calcul dans des matériaux de géométrie quelconque. Le matériau est plongé dans une boîte d'air parallélépipédique : le domaine fictif. L'algorithme présente la qualité première de la technique initiale : son efficacité - la complexité est en O(N log N) pour N noeuds de maillage au lieu de N$^2$ lorsqu'on fait le produit de convolution directement - avec, en plus, de la souplesse pour le choix de la géométrie.

On a testé ce code de calcul du champ démagnétisant sur un cube, une sphère et un pôle de tête d'enregistrement magnétique. Les résultats ont été comparés avec des calculs faits avec le logiciel Flux3D et publiés dans [44].

Un code résolvant l'équation de LLG complète 3D a été implémenté sur un calculateur parallèle, le Cray T3E.

Contrôle d'interfaces

Optimisation en cristallogenèse

Nous nous intéressons à une méthode de cristallogenèse : le procédé Bridgman. Il s'agit de contrôler l'interface de cristallisation à partir de paramètres agissant sur la thermique du procédé : la vitesse de tirage et l'orientation des réflecteurs pour les échanges de chaleur par rayonnement. À l'instant initial, le matériau doit être fondu de telle façon que le germe soit juste à la température de fusion. À l'instant final le matériau doit être tout entier en dessous de la température de solidification. En cours d'évolution la forme de l'interface doit rester stable. Le procédé est régi par une EDP non linéaire (Stefan avec changement de phase et conditions de flux non linéaire sur la frontière). Ce problème non linéaire est traité avec les méthodes développées par E. Magenes, C. Verdi et A. Visintin. La recherche des paramètres optimaux est un problème d'optimisation non convexe avec contraintes non linéaires.

Durant sa thèse, Mathias Marguliès a écrit un modèle du procédé et a étudié la résolution du problème d'optimisation par pénalisation et également par les multiplicateurs de Lagrange. Il utilise le calcul de l'adjoint adapté à la méthode de résolution choisie pour le problème direct. Il a mis au point un logiciel aidant au contrôle du procédé. Sa thèse a été soutenue en octobre 96. Les résultats font apparaître la prédominance des paramètres de contrôle du rayonnement et la nécessité d'affiner le modèle en prenant en compte la thermique dans le creuset.

La difficulté majeure du point de vue numérique a été le problème de minimisation. Si la fonctionnelle est différentiable, elle est non convexe et les contraintes portent sur les paramètres de contrôle et sur la variable d' état. Dans la lignée, J. Monnier et P. Witomski ont proposé un sujet de DEA pour essayer d'améliorer la minimisation pour ce type de fonctionnelles possédant de nombreux minima locaux. Dans son rapport [51] S. Putot a proposé un algorithme par pénalisation ``adaptative'' avec régularisation des profils de température de contrôle admissibles. Cette étude sera approfondie.

Problème de la détermination d'une forme optimale en cristallogenèse

C'est le sujet de thèse de S. Despreaux qui a débuté en octobre 1995. Les méthodes classiques de cristallogenèse à partir d'un bain liquide ne produisent en général que des monocristaux de forme cylindrique. La réalisation de géométries complexes nécessite la mise en oeuvre a posteriori d'étapes d'usinages délicats. Une alternative est la technique de préformage du cristal. Elle est basée sur l'utilisation des forces capillaires. Le matériau est d'abord fondu dans un creuset placé à l'intérieur d'un four. On réalise ensuite un ménisque liquide qui s'appuie sur une filière (frontière extérieure) et sur le cristal (frontière intérieure) qui se forme à partir d'un germe préorienté. Le problème est de savoir si l'on peut déterminer une forme de filière qui assure avec le cristal un contact liquide-solide à hauteur constante.

Cela nécessite évidemment de savoir résoudre le problème direct de la détermination de la forme du ménisque à frontières données. Il faut résoudre l'équation de la capillarité (elliptique non linéaire) avec des conditions aux limites de type Dirichlet sur la filière et Neumann non linéaire sur le cristal.

La première partie du travail a porté sur l'aspect optimisation de formes. Une étude des méthodes existantes, a permis de conclure que la technique basée sur les transformations de domaines était trés bien adaptée à ce cas. Après avoir mené à bien l'étude du cas discret, S. Despreaux a implanté cette méthode en modélisant les frontières à l'aide de splines cubiques, ce qui a permis d'appliquer différents algorithmes de descentes (Quasi-Newton, plus profonde descente, recherche du pas optimal par la méthode de Wolfe et par interpolation...). Des résultats encourageants ont ainsi pu être obtenus. La deuxième partie porte sur l'étude du cas continu. Pour résoudre l'équation de la capillarité en 2-D nous avons été amenés à d'abord étudier le cas axisymétrique. En utilisant un principe du maximum, il a ensuite été possible de montrer l'existence et l'unicité d'une solution moyennant une condition portant sur la courbure et l'épaisseur du domaine. On retrouve des conditions que nous avions imposé dans les essais numériques. Il s'agit d'une amélioration des résultats connus imposant des conditions de convexité sur le domaine. Nous travaillons actuellement sur l'existence d'un domaine optimal, dans la classe des domaines admissibles vérifiant ces conditions. La thèse de S. Despreaux devrait être soutenue dans le milieu de l'année 98.

Problèmes autour du phénomène de capillarité

La thèse que S. Bouchereau [7] a soutenue le 21 avril 1997 porte sur la modélisation et la simulation numérique de l'évolution d'un fluide sous l'effet conjugué des forces de tension superficielle et d'un champ électrique. Une goutte de liquide est posée sur un support conducteur. En l'absence de champ électrique la forme de l'interface résulte de l'équilibre entre les forces de tension superficielle et de gravitation. Lorsque l'on met la goutte sous une tension V que l'on fait croître, on observe une déformation de la surface libre. On cherche à calculer cette déformation et à mettre en évidence le blocage de la goutte que les physiciens observent expérimentalement pour certaines valeurs du potentiel.

Il s'agit d'un problème de contrôle optimal gouverné par une EDP elliptique dans un domaine non borné. La variable d'optimisation est la forme du domaine. S. Bouchereau a tout d'abord envisagé le cas axisymétrique, ce qui est raisonnable pour des valeurs du potentiel entre 0 et 600 volts. Il a écrit un code de calcul par éléments finis pour la résolution du problème direct et une méthode de type cheminement avec pour paramètre V. Les résultats peuvent être comparés pour des petites valeurs du potentiel avec ceux obtenus par une approximation condensateur plan habituellement faite par les physiciens. Cependant on constate les limites de la méthode numérique lorsque l'angle de contact entre la goutte et le substrat devient très faible. Le calcul du potentiel est un problème de transmission dans un ouvert à coin qui devient singulier. Pour aller au delà des 600 volts il était donc nécessaire de passer en 3-D et de considérer les singularités portées par la frontière. C'est la principale contribution de S. Bouchereau à ce problème. Il a étudié les singularités du problème direct et a mis au point une méthode de calcul par éléments frontières et intégrales singulières. La base des éléments finis est enrichie par des fonctions singulières sur le bord du domaine. Pour l'optimisation de formes la technique de transport par rapport à un domaine fixe est utilisée. La dérivation de la fonctionnelle discrète à minimiser est délicate en raison de la dépendance des singularités par rapport au domaine. La mise en oeuvre numérique représente un travail important et la prise en compte des intégrales singulières est délicate et coûteuse en calcul. Des résultats numériques probants ont été obtenus et le blocage de la goutte a pu être mis en évidence. Par contre les phénomènes d'ondulation de la ligne de mouillage n'ont pas pu être obtenus et restent un problème ouvert.

Méthodes numériques performantes pour l'environnement



Participants : Eric Blayo , Jacques Blum , Sylvain Carme , Mohamed Ghemires , Alain Le Breton , Laurent Debreu , François-Xavier Le Dimet , Marie-Christine Roubaud , Dinh Tuan Pham , Jacques Verron , Junqinq Yang


Mots-clés : assimilation de données, calcul parallèle, contrôle optimal, décomposition de domaine, filtrage, hydrologie, maillage adaptatif, météorologie, océanographie, prédicibilité, prévision, problème inverse


Résumé : Depuis les années 60, la modélisation numérique est utilisée en météorologie en vue de la prévision. Avec les progrès de la physique des modèles et le développement du potentiel de calcul, les avancées réalisées ont été très importantes. Dans la chaîne de la prévision numérique, le problème de l'assimilation de données (c'est-à-dire comment intégrer les données d'observation dans les modèles) est devenu un point crucial tant du point de vue de la nécessité de mettre en oeuvre des méthodes performantes que de celui du coût de calcul. Dans le domaine météorologique les méthodes variationnelles semblent devoir s'imposer pour le futur dans la plupart des centres importants : Météo-France, Centre Européeen à Reading, NCEP à Washington. Dans les autres disciplines de la géophysique, la modélisation numérique a été jusqu'alors principalement utilisée en vue de la compréhension des phénomènes. Il y a maintenant une orientation plus forte vers la prévision, comme par exemple en océanographie où un projet opérationnel de prévision (Mercator) se met en place, supporté par de nouvelles observations satellitaires (Topex-Poseidon). Le problème d'assimilation de données se pose donc maintenant dans ces disciplines avec des spécificités nouvelles. Il est clair qu'il y a là source d'importantes études méthodologiques et algorithmiques en raison de la difficulté des problèmes rencontrés (non-linéarité et très grandes dimensions). Le problème sera abordé à la fois par des méthodes variationnelles et par des méthodes de filtrage, le but final étant de comparer l'efficacité de ces deux approches.


\begin{figure} \begin{center} \includegraphics[width=6.cm]{idopt2.ps}\end{center}\end{figure}
Vue instantanée de la circulation océanique en surface dans la couronne circumpolaire antarctique (calculs effectués sur CRAY T3D)

Maillage adaptatif et méthodes de zoom pour les modèles de circulation océanique

Le problème de la résolution spatiale est un point crucial en modélisation des circulations océaniques. En effet, non seulement l'activité océanique de méso-échelle, c'est-à-dire aux échelles inférieures à 100-200 km, est importante, mais de plus ses interactions avec les mouvements à grande échelle sont fortes. Il est donc indispensable de résoudre effectivement cette turbulence de méso-échelle pour pouvoir espérer représenter correctement la circulation générale océanique.

Par ailleurs, une fine résolution spatiale peut également être nécessaire pour mieux prendre en compte des effets locaux : irrégularité d'un trait de côte, gradients bathymétriques, etc.

Enfin, on assiste actuellement au développement de l'océanographie opérationnelle avec la mise en place de plusieurs programmes de grande ampleur. C'est le cas notamment en France du projet MERCATOR, et de sa phase initiale, le projet CLIPPER (simulation haute résolution de l'océan Atlantique), dont le coût informatique est à la limite des possibilités actuelles. Un des besoins clairement identifiés dans le cadre de ce type de simulations est la possibilité de réaliser dans les modèles numériques des zooms à très haute résolution de régions présentant un intérêt particulier.

Dans ce contexte, l'intérêt de méthodes numériques permettant de réduire les coûts de calcul des modèles tout en conservant la possibilité de haute résolution locale est évident. Peu de travaux se sont à l'heure actuelle intéressés à ce problème pour les modèles océaniques, et il s'agissait dans tous les cas d'un couplage fixe (passif ou actif) entre un modèle grande échelle à basse ou moyenne résolution, centré sur la zone d'intérêt, et un modèle local à haute résolution de cette zone. Nous avons commencé cette année, dans le cadre de la thèse de L. Debreu, à travailler sur les méthodes de raffinement adaptatif pour les maillages structurés dans le contexte des modèles océaniques. De premiers tests dans des configurations académiques (modon barotrope, modèle QG multi-couches) ont mis clairement en évidence les potentialités de la méthode que nous utilisons. D'une part, nous obtenons pour un coût nettement inférieur (typiquement d'un facteur 3) des statistiques (circulation moyenne, niveaux d'énergie...) identiques à 10-20% près à celles obtenues avec une résolution uniformément fine. Par ailleurs, nos tests ont montré qu'on parvenait à une meilleure prévision locale à court terme des circulations grâce à un raffinement adaptatif plutôt qu'avec les méthodes de zoom employées classiquement ([25],[35]). Il convient maintenant de mener des validations plus poussées. Ceci sera réalisé notamment dans le cadre du modèle CLIPPER, pour tenter de répondre à deux questions : (i) les méthodes adaptatives peuvent-elles être utilisées dans le cadre d'intégrations longues des modèles ? et (ii) ces méthodes peuvent-elles fournir une alternative aux méthodes de zoom classiques pour le problème de la prédiction locale à court terme des circulations ? Par ailleurs, d'autres applications du raffinement de maillage sont envisagées. Ainsi par exemple, l'aspect adaptatif de ces méthodes peut fournir un moyen de calculer à coût réduit des statistiques approchées nécessaires dans le cadre des méthodes d'assimilation de données de type filtrage. Une collaboration sur ce thème avec l'équipe du LEGI est en cours d'élaboration.

Assimilation de données en océanographie

Le problème se pose pour l'océan dans les mêmes termes que pour l'atmosphère : on cherche à réaliser des simulations numériques en contraignant les solutions du modèle avec des observations afin de répondre à des objectifs de prévision en un sens déterministe ou probabiliste.

L'atmosphère et l'océan sont des fluides géophysiques différenciés dans leurs comportements dynamiques par des échelles caractéristiques différentes (les échelles spatiales sont plus fines dans l'océan tandis que les échelles temporelles y sont plus longues que dans l'atmosphère). En outre, les observations océaniques sont beaucoup moins denses (en temps et en espace) que les observations atmosphériques. Les nouvelles techniques satellitaires et notamment les mesures altimétriques doivent permettre d'améliorer notablement la connaissance des circulations océaniques.

Techniques déterministes de contrôle optimal

Dans la continuité des travaux menés ces dernières années sur les techniques d'assimilation variationnelle dans un modèle de circulation océanique ([13],[27]), nous avons abordé cette année le problème de la réduction d'ordre, dans le but de diminuer le coût de la méthode. Il faut rappeler en effet que le coût d'une assimilation de données par méthode adjointe est supérieur d'un ordre de grandeur à celui du modèle direct, ce qui rend l'utilisation de cette technique prohibitive dans des modèles océaniques opérationnels. Dans le cadre d'un stage de DEA [48], nous avons testé la possibilité de réduire la dimension de l'espace de contrôle, en travaillant non plus sur la condition initiale du modèle, mais sur une troncature de sa décomposition dans une base d'EOFs. On passe donc ainsi d'un espace de contrôle de dimension égale au nombre de points de grille du modèle à un espace de contrôle de dimension égale au nombre d'EOFs utilisées pour la reconstitution d'une pseudo condition initiale (typiquement 50 à 100). Les résultats obtenus sur un cas simple sont très encourageants. En effet, on diminue ainsi le coût de l'assimilation d'un facteur de l'ordre de 20, tout en conservant une bonne reconstitution de la trajectoire du modèle. En fait, les EOFs apportant une information statistique sur la structure verticale des écoulements, la trajectoire obtenue par cette méthode est même meilleure que celle obtenue par méthode classique en ce qui concerne la reconstitution des circulations océaniques profondes. Par contre, la reconstitution des écoulements de surface est légèrement dégradée, ce qui est logique dans la mesure où les données assimilées sont ici uniquement des données en surface.

Ce travail sera poursuivi cette année, et la méthode sera mise en pratique dans des cas plus complexes.

Techniques de filtrage stochastique

Une autre approche du problème consiste à utiliser les méthodes de type filtrage. Comme indiqué au §[*], un nouveau filtre (SEEK) a été développé dans le but de réduire les coûts de calcul très importants associés aux filtres classiques de type Kalman.

Dans ce cadre, S. Carme prépare une thèse de doctorat sur le développement et l'utilisation de méthodes de filtrage de rang réduit pour l'assimilation de données dans les modèles océaniques. L'objectif à court terme est de réaliser un logiciel opérationnel d'assimilation de données altimétriques satellitaires en Atlantique Nord, dans un modèle QG et utilisant le filtre mentionné plus haut. Cette réalisation est plus difficile que celle de Verron, Gourdeau et Pham [27] [25] mentionnée au §[*], car les modèles océaniques aux latitudes moyennes (par exemple en Atlantique Nord) présentent des non linéarités plus fortes. De plus on se sert ici de données réelles provenant du satellite. À plus long terme, nous envisageons d'implanter notre filtre dans des modèles aux équations primitives, à la physique plus complète. Une deuxième thèse est envisagée sur ce sujet.

Application de l'analyse de sensibilité à l'optimisation d'un système d'observation de l'océan

La prévision numérique de l'évolution de l'océan dépend des observations qui sont effectuées, et donc du système d'observation. Nous nous intéressons aux observations d'altimétrie satellitaires. Le but de ce travail est, dans un premier temps, la mise au point d'outils d'évaluation de la sensibilité de la reconstitution des champs par rapport aux observations. Traditionnellement, l'évaluation de la sensibilité se fait en utilisant l'adjoint du modèle. Mais l'approche que l'on a de l'assimilation de données nécessite de tenir compte aussi des données dans les études de sensibilité. Toute l'information (modèle et données) est contenue dans le système d'optimalité. C'est donc sur ce système que l'on doit mettre en oeuvre l'analyse de sensibilité. On a ainsi montré comment calculer la sensibilité de la condition initiale optimale (variable de contrôle), ainsi que celle d'une réponse fonctionnelle sous forme intégrale par rapport aux observations.

Nous nous plaçons ici dans le cadre d'un modèle déterministe d'océan (modèle QG de l'IMG, dont la version adjointe a été développée par B. Luong) et supposons que le modèle mathématique représente parfaitement l'océan. Des hypothèses moins restrictives pourront être envisagées par la suite si les méthodes mises au point sont satisfaisantes. Les outils mis en oeuvre sont les outils classiques de l'assimilation variationnelle (modèle adjoint et estimation du hessien), et ils demandent des moyens de calcul considérables.

Des expériences d'assimilation variationnelle de données altimétriques simulées ont d'ores et déjà été menées en 1994 et 1995, et une première analyse de sensibilité des champs reconstitués par rapport aux paramètres d'observation a été effectuée [12].

On cherche aussi une évaluation optimale des caractéristiques orbitales des satellites. Elles sont implicitement prises en compte au cours de l'assimilation par la projection de la solution du modèle sur les observations. Une étude est alors menée sur la sensibilité de la reconstitution des champs par rapport à cet opérateur de projection, le critère de sensibilité utilisé dans cette partie étant le conditionnement du hessien de la fonction-coût à minimiser.

Ce travail s'inscrit dans une démarche de prospection et nous a fourni des outils qui peuvent être exploités dans des applications, notamment sur la sensibilité des quantités diagnostiques de l'écoulement océanique par rapport aux observations. On peut aussi envisager le choix d'autres critères de sensibilité pour l'optimisation de la configuration spatio-temporelle des observations.

Couplage d'un modèle météorologique et d'un modèle chimique : étude d'algorithmes numériques performants

La prévision de l'évolution de la composition chimique de l'atmosphère nécessite :

$\bullet$ un modèle décrivant la dynamique de l'atmosphère;

$\bullet$ un modèle de diffusion des constituants chimiques;

$\bullet$ un modèle de cinétique chimique décrivant les interactions entre polluants;

$\bullet$ des observations des champs météorologiques et des champs de concentration des composants.

La prévision ne pourra se réaliser que si l'on associe toutes ces informations. Une des difficultés essentielles provient du volume important de calcul nécessaire pour mettre en oeuvre ces techniques. Notre premier travail était de tester et d'adapter un certain nombre d'algorithmes susceptibles d'être effectivement mis en oeuvre. Eu égard à l'importance des problèmes de calcul et à la nature du problème, les recours aux techniques de calcul parallèle sont indispensables.

Pour la partie convection et diffusion, nous avons utilisé une technique directe, une décomposition de domaine et nous avons présenté une méthode hybride basée sur la combinaison des deux premières. Pour la parallélisation de la cinétique chimique, nous avons utilisé un schéma maître-esclaves. Les performances obtenues sur Cray T3D montrent ces limites quand le nombre de processeurs augmente. Nous utilisons actuellement une répartition de charge pour essayer d'augmenter l'efficacité. Deux possiblités restent à traiter: décomposition de la chimie selon les vitesses de réactions et décomposition avec interface [8].

Identification de coefficients en hydrologie

Les modèles utilisés en hydrologie comportent des paramètres, tels que les coefficients de conductivité hydraulique, qui ne sont pas directement accessibles à la mesure. Les techniques d'identification de type contrôle optimal permettent une estimation de ces grandeurs à partir de ces observations. Un travail d'identification de coefficients a été poursuivi par J. Joergensen (étudiant stagiaire de l'université technique du Danemark) et Yann Le Floch (étudiant de DEA). Ce travail a été réalisé au laboratoire des Transferts en Hydrologie et Environnement de l'IMG avec la participation de G. Vachaud et R. Haverkamp dans le cadre du programme européen ECRASE. Par ailleurs les problèmes d'assimilation de données dans des modèles d'écoulement de surface ont été étudiés et mis en oeuvre par deux étudiants de l'IUP Mathématiques Appliquéees et Industrielles (S. Messina et A. Pastor au laboratoire de Génie de l'Environnement de l'université d'Oklahoma [50]). Dans ce même domaine de l'hydrologie Mlle Y. Junqing, doctorante codirigéee avec l'université de Wuhan (Chine) travaille sur l'assimilation de données pour des problèmes de sédimentation. La difficulté par rapport aux problèmes plus classiques de l'assimilation provient du fait que la géométrie du domaine change avec le temps d'où la nécessité de mettre en oeuvre des techniques d'optimisation par rapport au domaine. Ce travail est mené en collaboration avec l'Institut de Physique de l'Atmosphère de l'Académie des Sciences de Chine.


previous up next contents Précédent : Grands domaines d'application Remonter : Projet IDOPT, Identification et optimisation Suivant : Actions industrielles