13.1. Algorithme de calcul « 3DVAR »¶

13.1.1. Description¶

Cet algorithme réalise une estimation d’état par minimisation variationnelle de la fonctionnelle $J$ d’écart classique en assimilation de données statique:

$J(\mathbf{x})=(\mathbf{x}-\mathbf{x}^b)^T.\mathbf{B}^{-1}.(\mathbf{x}-\mathbf{x}^b)+(\mathbf{y}^o-H(\mathbf{x}))^T.\mathbf{R}^{-1}.(\mathbf{y}^o-H(\mathbf{x}))$

qui est usuellement désignée comme la fonctionnelle « 3D-Var » (voir par exemple [Talagrand97]). Les dénominations « 3D-Var », « 3D-VAR » et « 3DVAR » sont équivalentes.

Il existe diverses variantes de cet algorithme. On propose ici les formulations stables et robustes suivantes :

« 3DVAR » (3D Variational analysis, voir [Lorenc86], [LeDimet86], [Talagrand97]), algorithme classique d’origine, extrêmement robuste, opérant dans l’espace du modèle,
« 3DVAR-VAN » (3D Variational Analysis with No inversion of B, voir [Lorenc88]), algorithme similaire, opérant dans l’espace du modèle, permettant d’éviter l’inversion de la matrice de covariance B (sauf dans le cas où il y des bornes),
« 3DVAR-Incr » (Incremental 3DVAR, voir [Courtier94]), algorithme plus économique que les précédents, impliquant une approximation des opérateurs non-linéaires,
« 3DVAR-PSAS » (Physical-space Statistical Analysis Scheme for 3DVAR, voir [Courtier97], [Cohn98]), algorithme parfois plus économique car opérant dans l’espace des observations, impliquant une approximation des opérateurs non-linéaires, ne permettant pas la prise en compte de bornes.

On recommande fortement d’utiliser le « 3DVAR » d’origine. Les algorithmes « 3DVAR » et « 3DVAR-Incr » (et pas les autres) permettent explicitement la modification du point initial de leur minimisation, même si ce n’est pas recommandé.

Cet algorithme d’optimisation mono-objectif est naturellement écrit pour une estimation unique, sans notion dynamique ou itérative (il n’y a donc pas besoin dans ce cas d’opérateur d’évolution incrémentale, ni de covariance d’erreurs d’évolution). Dans le cadre traditionnel de l’assimilation de données temporelle ou itérative que traite ADAO, il peut aussi être utilisé sur une succession d’observations, plaçant alors l’estimation dans un cadre récursif similaire à un Algorithme de calcul « KalmanFilter ». Une estimation standard « 3D-VAR » est effectuée à chaque pas d’observation sur l’état prévu par le modèle d’évolution incrémentale, sachant que la covariance d’erreur d’état reste la covariance d’ébauche initialement fournie par l’utilisateur. Pour être explicite, contrairement aux filtres de type Kalman, la covariance d’erreurs sur les états n’est pas remise à jour.

Une extension du 3DVAR, couplant en parallèle une méthode 3DVAR, pour l’estimation d’un unique meilleur état, avec un filtre de Kalman d’ensemble pour l’estimation des covariances d’erreurs, permet d’améliorer l’estimation de ces covariances d’erreurs a posteriori. On atteint cette extension en utilisant le variant « E3DVAR » de l’algorithme de filtrage Algorithme de calcul « EnsembleKalmanFilter ».

On remarque que les statistiques d’erreurs d’observation et d’évolution sont par hypothèse représentées ici sous forme gaussienne. Alors, dans le cas particulier où l’opérateur d’observation $H$ est linéaire, cet algorithme est strictement équivalent à l’interpolation optimale (OI ou « Optimal Interpolation »). De plus, il réalise à la fois l’estimation de variance minimale (MV ou « Minimum Variance estimator ») et l’estimation de maximum a posteriori (MAP ou « Maximum A Posteriori estimator »), qui coïncident dans ce cas précis.

13.1.2. Quelques propriétés notables des méthodes implémentées¶

Pour compléter la description on synthétise ici quelques propriétés notables, des méthodes de l’algorithme ou de leurs implémentations. Ces propriétés peuvent avoir une influence sur la manière de l’utiliser ou sur ses performances de calcul. Pour de plus amples renseignements, on se reportera aux références plus complètes indiquées à la fin du descriptif de cet algorithme.

Les méthodes d’optimisation proposées par cet algorithme effectuent une recherche non locale du minimum, sans pour autant néanmoins assurer une recherche globale. C’est le cas lorsque les méthodes d’optimisation de recherche locale présentent de plus des capacités permettant d’éviter de rester bloquées par le premier minimum local trouvé. Ces capacités sont parfois heuristiques.

Les méthodes proposées par cet algorithme requièrent la dérivation de la fonction objectif ou de l’un des opérateurs. Cela nécessite que l’un au moins des opérateurs d’observation ou d’évolution soit différentiable voire les deux, et cela implique un temps de calcul supplémentaire dans le cas où les dérivées sont calculées numériquement par de multiples évaluations.

Les méthodes proposées par cet algorithme ne présentent pas de parallélisme interne, mais utilisent la dérivation numérique d’opérateur(s) qui est, elle, parallélisable. L’interaction potentielle, entre le parallélisme de la dérivation numérique, et le parallélisme éventuellement présent dans les opérateurs d’observation ou d’évolution intégrant les codes de l’utilisateur, doit donc être soigneusement réglée.

Les méthodes proposées par cet algorithme atteignent leur convergence sur un ou plusieurs critères de résidu ou de nombre. En pratique, il peut y avoir plusieurs critères de convergence actifs simultanément.

Le résidu peut être une mesure standard basée sur un écart (« écart calculs-mesures » par exemple), ou une valeur remarquable lié à l’algorithme (« nullité d’un gradient » par exemple).

Le nombre est fréquemment un élément remarquable lié à l’algorithme, comme un nombre d’itérations ou un nombre d’évaluations, mais cela peut aussi être par exemple un nombre de générations pour un algorithme évolutionnaire.

Il convient de régler soigneusement les seuils de convergence, pour limiter le coût calcul global de l’algorithme, ou pour assurer une adaptation de la convergence au cas physique traité.

13.1.3. Commandes requises et optionnelles¶

Les commandes générales requises, disponibles en édition dans l’interface graphique ou textuelle, sont les suivantes :

Background: Vecteur. La variable désigne le vecteur d’ébauche ou d’initialisation, usuellement noté $\mathbf{x}^b$ . Sa valeur est définie comme un objet de type « Vector » ou « VectorSerie ». Sa disponibilité en sortie est conditionnée par le booléen « Stored » associé en entrée.

BackgroundError: Matrice. La variable désigne la matrice de covariance des erreurs d’ébauche, usuellement notée $\mathbf{B}$ . Sa valeur est définie comme un objet de type « Matrix », de type « ScalarSparseMatrix », ou de type « DiagonalSparseMatrix », comme décrit en détail dans la section Conditions requises pour décrire des matrices de covariance. Sa disponibilité en sortie est conditionnée par le booléen « Stored » associé en entrée.

EvolutionError: Matrice. La variable désigne la matrice de covariance des erreurs a priori d’évolution, usuellement notée $\mathbf{Q}$ . Sa valeur est définie comme un objet de type « Matrix », de type « ScalarSparseMatrix », ou de type « DiagonalSparseMatrix », comme décrit en détail dans la section Conditions requises pour décrire des matrices de covariance. Sa disponibilité en sortie est conditionnée par le booléen « Stored » associé en entrée.

EvolutionModel: Opérateur. La variable désigne l’opérateur d’évolution du modèle, usuellement noté $M$ , qui décrit un pas élémentaire d’évolution dynamique ou itérative. Sa valeur est définie comme un objet de type « Function » ou de type « Matrix ». Dans le cas du type « Function », différentes formes fonctionnelles peuvent être utilisées, comme décrit dans la section Conditions requises pour les fonctions décrivant un opérateur. Si un contrôle $U$ est inclus dans le modèle d’observation, l’opérateur doit être appliqué à une paire $(X,U)$ .

Observation: Liste de vecteurs. La variable désigne le vecteur d’observation utilisé en assimilation de données ou en optimisation, et usuellement noté $\mathbf{y}^o$ . Sa valeur est définie comme un objet de type « Vector » si c’est une unique observation (temporelle ou pas) ou « VectorSerie » si c’est une succession d’observations. Sa disponibilité en sortie est conditionnée par le booléen « Stored » associé en entrée.

ObservationError: Matrice. La variable désigne la matrice de covariance des erreurs a priori d’ébauche, usuellement notée $\mathbf{R}$ . Cette matrice est définie comme un objet de type « Matrix », de type « ScalarSparseMatrix », ou de type « DiagonalSparseMatrix », comme décrit en détail dans la section Conditions requises pour décrire des matrices de covariance. Sa disponibilité en sortie est conditionnée par le booléen « Stored » associé en entrée.

ObservationOperator: Opérateur. La variable désigne l’opérateur d’observation, usuellement noté $H$ , qui transforme les paramètres d’entrée $\mathbf{x}$ en résultats $\mathbf{y}$ qui sont à comparer aux observations $\mathbf{y}^o$ . Sa valeur est définie comme un objet de type « Function » ou de type « Matrix ». Dans le cas du type « Function », différentes formes fonctionnelles peuvent être utilisées, comme décrit dans la section Conditions requises pour les fonctions décrivant un opérateur. Si un contrôle $U$ est inclus dans le modèle d’observation, l’opérateur doit être appliqué à une paire $(X,U)$ .

Les commandes optionnelles générales, disponibles en édition dans l’interface graphique ou textuelle, sont indiquées dans la Liste des commandes et mots-clés pour un cas d’assimilation de données ou d’optimisation. De plus, les paramètres de la commande « AlgorithmParameters » permettent d’indiquer les options particulières, décrites ci-après, de l’algorithme. On se reportera à la Description des options d’un algorithme par « AlgorithmParameters » pour le bon usage de cette commande.

Les options sont les suivantes :

Bounds

Liste de paires de valeurs réelles. Cette clé permet de définir des paires de bornes supérieure et inférieure pour chaque variable d’état optimisée. Les bornes doivent être données par une liste de liste de paires de bornes inférieure/supérieure pour chaque variable, avec une valeur None chaque fois qu’il n’y a pas de borne. Les bornes peuvent toujours être spécifiées, mais seuls les optimiseurs sous contraintes les prennent en compte. Si la liste est vide, cela équivaut à une absence de bornes.

Exemple : {"Bounds":[[2.,5.],[1.e-2,10.],[-30.,None],[None,None]]}

CostDecrementTolerance

Valeur réelle. Cette clé indique une valeur limite, conduisant à arrêter le processus itératif d’optimisation lorsque la fonction coût décroît moins que cette tolérance au dernier pas. La valeur par défaut est de 1.e-7, et il est recommandé de l’adapter aux besoins pour des problèmes réels. On peut se reporter à la partie décrivant les manières de Contrôler la convergence pour des cas de calculs et algorithmes itératifs pour des recommandations plus détaillées.

Exemple : {"CostDecrementTolerance":1.e-7}

EstimationOf

Nom prédéfini. Cette clé permet de choisir le type d’estimation à réaliser. Cela peut être soit une estimation de l’état, avec la valeur « State », ou une estimation de paramètres, avec la valeur « Parameters ». Le choix par défaut est « Parameters ».

Exemple : {"EstimationOf":"Parameters"}

GradientNormTolerance

Valeur réelle. Cette clé indique une valeur limite, conduisant à arrêter le processus itératif d’optimisation lorsque la norme du gradient est en dessous de cette limite. C’est utilisé uniquement par les optimiseurs sans contraintes. Le défaut est 1.e-5 et il n’est pas recommandé de le changer.

Exemple : {"GradientNormTolerance":1.e-5}

InitializationPoint

Vecteur. La variable désigne un vecteur à utiliser comme l’état initial autour duquel démarre un algorithme itératif. Par défaut, cet état initial n’a pas besoin d’être fourni et il est égal à l’ébauche $\mathbf{x}^b$ . Sa valeur doit permettre de construire un vecteur de taille identique à l’ébauche. Dans le cas où il est fourni, il ne remplace l’ébauche que pour l’initialisation.

Exemple : {"InitializationPoint":[1, 2, 3, 4, 5]}

MaximumNumberOfIterations

Valeur entière. Cette clé indique le nombre maximum d’itérations internes possibles en optimisation itérative. Le défaut est 15000, qui est très similaire à une absence de limite sur ce nombre d’itérations. Il est ainsi recommandé d’adapter ce paramètre aux besoins pour des problèmes réels. Pour certains optimiseurs, le nombre de pas effectif d’arrêt peut être légèrement différent de la limite à cause d’exigences de contrôle interne de l’algorithme. On peut se reporter à la partie décrivant les manières de Contrôler la convergence pour des cas de calculs et algorithmes itératifs pour des recommandations plus détaillées.

Exemple : {"MaximumNumberOfIterations":100}

Minimizer: Nom prédéfini. Cette clé permet de changer le minimiseur pour l’optimiseur. Le choix par défaut est « LBFGSB », et les choix possibles sont « LBFGSB » (minimisation non linéaire sous contraintes, voir [Byrd95], [Morales11], [Zhu97]), « TNC » (minimisation non linéaire sous contraintes), « CG » (minimisation non linéaire sans contraintes), « BFGS » (minimisation non linéaire sans contraintes), Il est fortement conseillé de conserver la valeur par défaut.

NumberOfSamplesForQuantiles

Valeur entière. Cette clé indique le nombre de simulations effectuées pour estimer les quantiles. Cette option n’est utile que si le calcul supplémentaire « SimulationQuantiles » a été choisi. Le défaut est 100, ce qui suffit souvent pour une estimation correcte de quantiles courants à 5%, 10%, 90% ou 95%.

Exemple : {"NumberOfSamplesForQuantiles":100}

ProjectedGradientTolerance

Valeur réelle. Cette clé indique une valeur limite, conduisant à arrêter le processus itératif d’optimisation lorsque toutes les composantes du gradient projeté sont en-dessous de cette limite. C’est utilisé uniquement par les optimiseurs sous contraintes. Le défaut est -1, qui désigne le défaut interne de chaque optimiseur (usuellement 1.e-5), et il n’est pas recommandé de le changer.

Exemple : {"ProjectedGradientTolerance":-1}

Quantiles

Liste de valeurs réelles. Cette liste indique les valeurs de quantile, entre 0 et 1, à estimer par simulation autour de l’état optimal. L’échantillonnage utilise des tirages aléatoires gaussiens multivariés, dirigés par la matrice de covariance a posteriori. Cette option n’est utile que si le calcul supplémentaire « SimulationQuantiles » a été choisi. La valeur par défaut est une liste vide.

Exemple : {"Quantiles":[0.1,0.9]}

SetSeed

Valeur entière. Cette clé permet de donner un nombre entier pour fixer la graine du générateur aléatoire utilisé dans l’algorithme. Par défaut, la graine est laissée non initialisée, et elle utilise ainsi l’initialisation par défaut de l’ordinateur, qui varie donc à chaque étude. Pour assurer la reproductibilité de résultats impliquant des tirages aléatoires, il est fortement conseillé d’initialiser la graine. Une valeur simple est par exemple 123456789. Il est conseillé de mettre un entier à plus de 6 ou 7 chiffres pour bien initialiser le générateur aléatoire.

Exemple : {"SetSeed":123456789}

SimulationForQuantiles

Nom prédéfini. Cette clé indique le type de simulation, linéaire (avec l’opérateur d’observation tangent appliqué sur des incréments de perturbations autour de l’état optimal) ou non-linéaire (avec l’opérateur d’observation standard appliqué aux états perturbés), que l’on veut faire pour chaque perturbation. Cela change essentiellement le temps de chaque simulation élémentaire, usuellement plus long en non-linéaire qu’en linéaire. Cette option n’est utile que si le calcul supplémentaire « SimulationQuantiles » a été choisi. La valeur par défaut est « Linear », et les choix possibles sont « Linear » et « NonLinear ».

Exemple : {"SimulationForQuantiles":"Linear"}

StateBoundsForQuantiles

Liste de paires de valeurs réelles. Cette clé permet de définir des paires de bornes supérieure et inférieure pour chaque variable d’état utilisée dans la simulation des quantiles. Les bornes doivent être données par une liste de liste de paires de bornes inférieure/supérieure pour chaque variable, avec une valeur None chaque fois qu’il n’y a pas de borne.

En l’absence de définition de ces bornes pour la simulation des quantiles et si des bornes d’optimisation sont définies, ce sont ces dernières qui sont utilisées pour la simulation des quantiles. Si ces bornes pour la simulation des quantiles sont définies, elles sont utilisées quelles que soient les bornes d’optimisation définies. Si cette variable est définie à None, alors aucune borne n’est utilisée pour les états utilisés dans la simulation des quantiles quelles que soient les bornes d’optimisation définies.

Exemple : {"StateBoundsForQuantiles":[[2.,5.],[1.e-2,10.],[-30.,None],[None,None]]}

StoreSupplementaryCalculations

Liste de noms. Cette liste indique les noms des variables supplémentaires, qui peuvent être disponibles au cours du déroulement ou à la fin de l’algorithme, si elles sont initialement demandées par l’utilisateur. Leur disponibilité implique, potentiellement, des calculs ou du stockage coûteux. La valeur par défaut est donc une liste vide, aucune de ces variables n’étant calculée et stockée par défaut (sauf les variables inconditionnelles). Les noms possibles pour les variables supplémentaires sont dans la liste suivante (la description détaillée de chaque variable nommée est donnée dans la suite de cette documentation par algorithme spécifique, dans la sous-partie « Informations et variables disponibles à la fin de l’algorithme ») : [ « Analysis », « APosterioriCorrelations », « APosterioriCovariance », « APosterioriStandardDeviations », « APosterioriVariances », « BMA », « CostFunctionJ », « CostFunctionJAtCurrentOptimum », « CostFunctionJb », « CostFunctionJbAtCurrentOptimum », « CostFunctionJo », « CostFunctionJoAtCurrentOptimum », « CurrentIterationNumber », « CurrentOptimum », « CurrentState », « CurrentStepNumber », « EnsembleOfSimulations », « EnsembleOfStates », « ForecastState », « IndexOfOptimum », « Innovation », « InnovationAtCurrentAnalysis », « InnovationAtCurrentState », « JacobianMatrixAtBackground », « JacobianMatrixAtOptimum », « KalmanGainAtOptimum », « MahalanobisConsistency », « OMA », « OMB », « SampledStateForQuantiles », « SigmaObs2 », « SimulatedObservationAtBackground », « SimulatedObservationAtCurrentOptimum », « SimulatedObservationAtCurrentState », « SimulatedObservationAtOptimum », « SimulationQuantiles », ].

Exemple : {"StoreSupplementaryCalculations":["CurrentState", "Residu"]}

Variant

Nom prédéfini. Cette clé permet de choisir l’une des variantes possibles pour l’algorithme principal. La variante par défaut est le « 3DVAR » d’origine, et les choix possibles sont « 3DVAR » (3D Variational analysis classique), « 3DVAR-VAN » (3D Variational Analysis with No inversion of B), « 3DVAR-Incr » (Incremental 3DVAR), « 3DVAR-PSAS » (Physical-space Statistical Analysis Scheme for 3DVAR), Il est fortement recommandé de conserver la valeur par défaut.

Exemple : {"Variant":"3DVAR"}

13.1.4. Informations et variables disponibles à la fin de l’algorithme¶

En sortie, après exécution de l’algorithme, on dispose d’informations et de variables issues du calcul. La description des Variables et informations disponibles en sortie indique la manière de les obtenir, par la méthode nommée get, depuis la variable « ADD » du post-processing en interface graphique, ou depuis le cas en interface textuelle. Les variables d’entrée, mises à disposition de l’utilisateur en sortie pour faciliter l’écriture des procédures de post-processing, sont décrites dans un Inventaire des informations potentiellement disponibles en sortie.

Sorties permanentes (non conditionnelles)

Les sorties non conditionnelles de l’algorithme sont les suivantes :

Analysis

Liste de vecteurs. Chaque élément de cette variable est un état optimal $\mathbf{x}^*$ en optimisation, une interpolation ou une analyse $\mathbf{x}^a$ en assimilation de données.

Exemple : xa = ADD.get("Analysis")[-1]

CostFunctionJ

Liste de valeurs. Chaque élément est une valeur de fonctionnelle d’écart $J$ choisie.

Exemple : J = ADD.get("CostFunctionJ")[:]

CostFunctionJb

Liste de valeurs. Chaque élément est une valeur de fonctionnelle d’écart $J^b$ , c’est-à-dire de la partie écart à l’ébauche. Si cette partie n’existe pas dans la fonctionnelle, sa valeur est nulle.

Exemple : Jb = ADD.get("CostFunctionJb")[:]

CostFunctionJo

Liste de valeurs. Chaque élément est une valeur de fonctionnelle d’écart $J^o$ , c’est-à-dire de la partie écart à l’observation.

Exemple : Jo = ADD.get("CostFunctionJo")[:]

Ensemble des sorties à la demande (conditionnelles ou non)

L’ensemble des sorties (conditionnelles ou non) de l’algorithme, classées par ordre alphabétique, est le suivant :

Analysis

Liste de vecteurs. Chaque élément de cette variable est un état optimal $\mathbf{x}^*$ en optimisation, une interpolation ou une analyse $\mathbf{x}^a$ en assimilation de données.

Exemple : xa = ADD.get("Analysis")[-1]

APosterioriCorrelations

Liste de matrices. Chaque élément est une matrice de corrélations des erreurs a posteriori de l’état optimal, issue de la matrice $\mathbf{A}$ des covariances. Pour en disposer, il faut avoir en même temps demandé le calcul de ces covariances d’erreurs a posteriori.

Exemple : apc = ADD.get("APosterioriCorrelations")[-1]

APosterioriCovariance

Liste de matrices. Chaque élément est une matrice $\mathbf{A}$ de covariances des erreurs a posteriori de l’état optimal.

Exemple : apc = ADD.get("APosterioriCovariance")[-1]

APosterioriStandardDeviations

Liste de matrices. Chaque élément est une matrice diagonale d’écarts-types des erreurs a posteriori de l’état optimal, issue de la matrice $\mathbf{A}$ des covariances. Pour en disposer, il faut avoir en même temps demandé le calcul de ces covariances d’erreurs a posteriori.

Exemple : aps = ADD.get("APosterioriStandardDeviations")[-1]

APosterioriVariances

Liste de matrices. Chaque élément est une matrice diagonale de variances des erreurs a posteriori de l’état optimal, issue de la matrice $\mathbf{A}$ des covariances. Pour en disposer, il faut avoir en même temps demandé le calcul de ces covariances d’erreurs a posteriori.

Exemple : apv = ADD.get("APosterioriVariances")[-1]

BMA

Liste de vecteurs. Chaque élément est un vecteur d’écart entre l’ébauche et l’état optimal.

Exemple : bma = ADD.get("BMA")[-1]

CostFunctionJ

Liste de valeurs. Chaque élément est une valeur de fonctionnelle d’écart $J$ choisie.

Exemple : J = ADD.get("CostFunctionJ")[:]

CostFunctionJAtCurrentOptimum

Liste de valeurs. Chaque élément est une valeur de fonctionnelle d’écart $J$ . A chaque pas, la valeur correspond à l’état optimal trouvé depuis le début.

Exemple : JACO = ADD.get("CostFunctionJAtCurrentOptimum")[:]

CostFunctionJb

Liste de valeurs. Chaque élément est une valeur de fonctionnelle d’écart $J^b$ , c’est-à-dire de la partie écart à l’ébauche. Si cette partie n’existe pas dans la fonctionnelle, sa valeur est nulle.

Exemple : Jb = ADD.get("CostFunctionJb")[:]

CostFunctionJbAtCurrentOptimum

Liste de valeurs. Chaque élément est une valeur de fonctionnelle d’écart $J^b$ , c’est-à-dire de la partie écart à l’ébauche. A chaque pas, la valeur correspond à l’état optimal trouvé depuis le début. Si cette partie n’existe pas dans la fonctionnelle, sa valeur est nulle.

Exemple : JbACO = ADD.get("CostFunctionJbAtCurrentOptimum")[:]

CostFunctionJo

Liste de valeurs. Chaque élément est une valeur de fonctionnelle d’écart $J^o$ , c’est-à-dire de la partie écart à l’observation.

Exemple : Jo = ADD.get("CostFunctionJo")[:]

CostFunctionJoAtCurrentOptimum

Liste de valeurs. Chaque élément est une valeur de fonctionnelle d’écart $J^o$ , c’est-à-dire de la partie écart à l’observation. A chaque pas, la valeur correspond à l’état optimal trouvé depuis le début.

Exemple : JoACO = ADD.get("CostFunctionJoAtCurrentOptimum")[:]

CurrentIterationNumber

Liste d’entiers. Chaque élément est l’index d’itération courant au cours du déroulement itératif de l’algorithme utilisé. Il y a une valeur d’index d’itération par pas d’assimilation correspondant à un état observé.

Exemple : cin = ADD.get("CurrentIterationNumber")[-1]

CurrentOptimum

Liste de vecteurs. Chaque élément est le vecteur d’état optimal au pas de temps courant au cours du déroulement itératif de l’algorithme d’optimisation utilisé. Ce n’est pas nécessairement le dernier état.

Exemple : xo = ADD.get("CurrentOptimum")[:]

CurrentState

Liste de vecteurs. Chaque élément est un vecteur d’état courant utilisé au cours du déroulement itératif de l’algorithme utilisé.

Exemple : xs = ADD.get("CurrentState")[:]

CurrentStepNumber

Liste d’entiers. Chaque élément est l’index du pas courant au cours du déroulement itératif, piloté par la série des observations, de l’algorithme utilisé. Cela correspond au pas d’observation utilisé. Remarque : ce n’est pas l’index d’itération courant d’algorithme même si cela coïncide pour des algorithmes non itératifs.

Exemple : csn = ADD.get("CurrentStepNumber")[-1]

EnsembleOfSimulations

Liste de vecteurs ou matrice. Chaque élément est une collection ordonnée de vecteurs d’état physique ou d’état simulé éventuellement observé $\mathbf{y}$ . Ce sont des sorties d’opérateur $H$ , c’est-à-dire des états d’observation simulés (nommés « snapshots » en terminologie de bases réduites). A chaque index de pas, il y a 1 état par colonne si cette liste est sous forme matricielle, ou 1 état par élément si c’est effectivement une liste. Important : la numérotation du support ou des points, sur lequel ou auxquels sont fournis une valeur d’état dans chaque vecteur, est implicitement celle de l’ordre naturel de numérotation du vecteur d’état, de 0 à la « taille moins 1 » de ce vecteur.

Exemple : {"EnsembleOfSimulations":[y1, y2, y3...]}

EnsembleOfStates

Liste de vecteurs ou matrice. Chaque élément est une collection ordonnée de vecteurs d’état physique ou d’état paramétrique $\mathbf{x}$ . Ce sont des entrées d’opérateur $H$ , c’est-à-dire des états courants avant observation. A chaque index de pas, il y a 1 état par colonne si cette liste est sous forme matricielle, ou 1 état par élément si c’est effectivement une liste. Important : la numérotation du support ou des points, sur lequel ou auxquels sont fournis une valeur d’état dans chaque vecteur, est implicitement celle de l’ordre naturel de numérotation du vecteur d’état, de 0 à la « taille moins 1 » de ce vecteur.

Exemple : {"EnsembleOfStates":[x1, x2, x3...]}

ForecastState

Liste de vecteurs. Chaque élément est un vecteur d’état (ou un ensemble de vecteurs d’états selon l’algorithme) prévu(s) par le modèle au cours du déroulement itératif temporel de l’algorithme utilisé.

Exemple : xf = ADD.get("ForecastState")[:]

IndexOfOptimum

Liste d’entiers. Chaque élément est l’index d’itération de l’optimum obtenu au cours du déroulement itératif de l’algorithme d’optimisation utilisé. Ce n’est pas nécessairement le numéro de la dernière itération.

Exemple : ioo = ADD.get("IndexOfOptimum")[-1]

Innovation

Liste de vecteurs. Chaque élément est un vecteur d’innovation, qui est en statique l’écart de l’optimum à l’ébauche, et en dynamique l’incrément d’évolution.

Exemple : d = ADD.get("Innovation")[-1]

InnovationAtCurrentAnalysis

Liste de vecteurs. Chaque élément est un vecteur d’innovation à l’état analysé courant. Cette quantité est identique au vecteur d’innovation à l’état analysé dans le cas d’une assimilation mono-état.

Exemple : da = ADD.get("InnovationAtCurrentAnalysis")[-1]

InnovationAtCurrentState

Liste de vecteurs. Chaque élément est un vecteur d’innovation à l’état courant avant analyse.

Exemple : ds = ADD.get("InnovationAtCurrentState")[-1]

JacobianMatrixAtBackground

Liste de matrices. Chaque élément est une matrice jacobienne de dérivées partielles de la sortie de l’opérateur d’observation par rapport aux paramètres d’entrée, une colonne de dérivées par paramètre. Le calcul est effectué à l’état initial.

Exemple: gradh = ADD.get("JacobianMatrixAtBackground")[-1]

JacobianMatrixAtOptimum

Liste de matrices. Chaque élément est une matrice jacobienne de dérivées partielles de la sortie de l’opérateur d’observation par rapport aux paramètres d’entrée, une colonne de dérivées par paramètre. Le calcul est effectué à l’état optimal.

Exemple: gradh = ADD.get("JacobianMatrixAtOptimum")[-1]

KalmanGainAtOptimum

Liste de matrices. Chaque élément est une matrice de gain de Kalman standard, évaluée à l’aide de l’opérateur d’observation linéarisé. Le calcul est effectué à l’état optimal.

Exemple: kg = ADD.get("KalmanGainAtOptimum")[-1]

MahalanobisConsistency

Liste de valeurs. Chaque élément est une valeur de l’indicateur de qualité de Mahalanobis.

Exemple : mc = ADD.get("MahalanobisConsistency")[-1]

OMA

Liste de vecteurs. Chaque élément est un vecteur d’écart entre l’observation et l’état optimal dans l’espace des observations.

Exemple : oma = ADD.get("OMA")[-1]

OMB

Liste de vecteurs. Chaque élément est un vecteur d’écart entre l’observation et l’état d’ébauche dans l’espace des observations.

Exemple : omb = ADD.get("OMB")[-1]

SampledStateForQuantiles

Liste de séries de vecteurs. Chaque élément est une série de vecteurs d’état colonnes, généré pour estimer par simulation et/ou observation les valeurs de quantiles requis par l’utilisateur. Il y a autant d’états que le nombre d’échantillons requis pour cette estimation de quantiles.

Exemple : xq = ADD.get("SampledStateForQuantiles")[:]

SigmaObs2

Liste de valeurs. Chaque élément est une valeur de l’indicateur de qualité $(\sigma^o)^2$ de la partie observation.

Exemple : so2 = ADD.get("SigmaObs")[-1]

SimulatedObservationAtBackground

Liste de vecteurs. Chaque élément est un vecteur d’observation simulé par l’opérateur d’observation à partir de l’ébauche $\mathbf{x}^b$ . C’est la prévision à partir de l’ébauche, elle est parfois appelée « Dry ».

Exemple : hxb = ADD.get("SimulatedObservationAtBackground")[-1]

SimulatedObservationAtCurrentOptimum

Liste de vecteurs. Chaque élément est un vecteur d’observation simulé par l’opérateur d’observation à partir de l’état optimal au pas de temps courant au cours du déroulement de l’algorithme d’optimisation, c’est-à-dire dans l’espace des observations.

Exemple : hxo = ADD.get("SimulatedObservationAtCurrentOptimum")[-1]

SimulatedObservationAtCurrentState

Liste de vecteurs. Chaque élément est un vecteur d’observation simulé par l’opérateur d’observation à partir de l’état courant, c’est-à-dire dans l’espace des observations.

Exemple : hxs = ADD.get("SimulatedObservationAtCurrentState")[-1]

SimulatedObservationAtOptimum

Liste de vecteurs. Chaque élément est un vecteur d’observation obtenu par l’opérateur d’observation à partir de la simulation d’analyse ou d’état optimal $\mathbf{x}^a$ . C’est l’observation de la prévision à partir de l’analyse ou de l’état optimal, et elle est parfois appelée « Forecast ».

Exemple : hxa = ADD.get("SimulatedObservationAtOptimum")[-1]

SimulationQuantiles

Liste de séries de vecteurs. Chaque élément est une série de vecteurs colonnes d’observation, correspondant, pour un quantile particulier requis par l’utilisateur, à l’état observé qui réalise le quantile demandé. Chaque vecteur colonne d’observation est restitué dans le même ordre que les valeurs de quantiles requis par l’utilisateur.

Exemple : sQuantiles = ADD.get("SimulationQuantiles")[:]

13.1.5. Exemples d’utilisation en Python (TUI)¶

Voici un ou des exemples très simple d’usage de l’algorithme proposé et de ses paramètres, écrit en [DocR] Interface textuelle pour l’utilisateur (TUI/API). De plus, lorsque c’est possible, les informations indiquées en entrée permettent aussi de définir un cas équivalent en interface graphique [DocR] Interface graphique pour l’utilisateur (GUI/EFICAS).

13.1.5.1. Premier exemple¶

Cet exemple décrit le recalage des paramètres $\mathbf{x}$ d’un modèle d’observation $H$ quadratique. Ce modèle $H$ est représenté ici comme une fonction nommée QuadFunction pour les besoins de l’exemple. Cette fonction accepte en entrée le vecteur de coefficients $\mathbf{x}$ de la forme quadratique, et fournit en sortie le vecteur $\mathbf{y}$ d’évaluation du modèle quadratique aux points de contrôle internes, prédéfinis de manière statique dans le modèle.

Le recalage s’effectue sur la base d’un jeu initial de coefficients (état d’ébauche désigné par Xb dans l’exemple), et avec l’information $\mathbf{y}^o$ (désignée par Yobs dans l’exemple) de 5 mesures obtenues aux points de contrôle internes. On se place ici en expériences jumelles (voir Pour tester une chaîne d’assimilation de données : les expériences jumelles), et les mesures sont considérées comme parfaites. On choisit de privilégier les observations, au détriment de l’ébauche, par l’indication artificielle d’une très importante variance d’erreur d’ébauche, ici de $10^{6}$ .

L’ajustement s’effectue en affichant des résultats intermédiaires lors de l’optimisation itérative grâce à des « observer » (pour mémoire, voir Obtenir des informations sur des variables spéciales au cours d’un calcul ADAO).

# -*- coding: utf-8 -*-
#
from numpy import array, ravel
def QuadFunction( coefficients ):
    """
    Simulation quadratique aux points x : y = a x^2 + b x + c
    """
    a, b, c = list(ravel(coefficients))
    x_points = (-5, 0, 1, 3, 10)
    y_points = []
    for x in x_points:
        y_points.append( a*x*x + b*x + c )
    return array(y_points)
#
Xb   = array([1., 1., 1.])
Yobs = array([57, 2, 3, 17, 192])
#
print("Résolution du problème de calage")
print("--------------------------------")
print("")
from adao import adaoBuilder
case = adaoBuilder.New()
case.setBackground( Vector = Xb, Stored=True )
case.setBackgroundError( ScalarSparseMatrix = 1.e6 )
case.setObservation( Vector = Yobs, Stored=True )
case.setObservationError( ScalarSparseMatrix = 1. )
case.setObservationOperator( OneFunction = QuadFunction )
case.setAlgorithmParameters(
    Algorithm="3DVAR",
    Parameters={
        "MaximumNumberOfIterations": 100,
        "StoreSupplementaryCalculations": [
            "CurrentState",
            "OMA",
            ],
        },
    )
case.setObserver(
    Info="  État intermédiaire en itération courante :",
    Template="ValuePrinter",
    Variable="CurrentState",
    )
case.execute()
print("")
#
#-------------------------------------------------------------------------------
#
print("Calage de %i coefficients pour une forme quadratique 1D sur %i mesures"%(
    len(case.get("Background")),
    len(case.get("Observation")),
    ))
print("--------------------------------------------------------------------")
print("")
print("Vecteur d'observation.............:", ravel(case.get("Observation")))
print("État d'ébauche a priori...........:", ravel(case.get("Background")))
print("")
print("Coefficients théoriques attendus..:", ravel((2,-1,2)))
print("")
print("Nombre d'itérations...............:", len(case.get("CurrentState")))
print("Nombre de simulations.............:", len(case.get("CurrentState"))*4)
print("Écart Observation-Analyse (OMA)...:", ravel(case.get("OMA")[-1]))
print("Coefficients résultants du calage.:", ravel(case.get("Analysis")[-1]))
#
Xa = case.get("Analysis")[-1]
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (10, 4)
#
plt.figure()
plt.plot((-5,0,1,3,10),QuadFunction(Xb),"b--",label="Simulation à l'ébauche")
plt.plot((-5,0,1,3,10),Yobs,            "kX", label="Observation",markersize=10)
plt.plot((-5,0,1,3,10),QuadFunction(Xa),"r-", label="Simulation à l'optimum")
plt.legend()
plt.title("Calage de coefficients", fontweight="bold")
plt.xlabel("Coordonnée arbitraire")
plt.ylabel("Observations")
plt.savefig("simple_3DVAR1.png")

Le résultat de son exécution est le suivant :

Résolution du problème de calage
--------------------------------

  État intermédiaire en itération courante : [1. 1. 1.]
  État intermédiaire en itération courante : [1.99739508 1.07086406 1.01346638]
  État intermédiaire en itération courante : [1.83891966 1.04815981 1.01208385]
  État intermédiaire en itération courante : [1.8390702  1.03667176 1.01284797]
  État intermédiaire en itération courante : [1.83967236 0.99071957 1.01590445]
  État intermédiaire en itération courante : [1.84208099 0.8069108  1.02813037]
  État intermédiaire en itération courante : [ 1.93711599 -0.56383145  1.12097995]
  État intermédiaire en itération courante : [ 1.99838848 -1.00480573  1.1563713 ]
  État intermédiaire en itération courante : [ 2.0135905  -1.04815933  1.16155285]
  État intermédiaire en itération courante : [ 2.01385679 -1.03874809  1.16129657]
  État intermédiaire en itération courante : [ 2.01377856 -1.03700044  1.16157611]
  État intermédiaire en itération courante : [ 2.01338902 -1.02943736  1.16528951]
  État intermédiaire en itération courante : [ 2.01265633 -1.0170847   1.17793974]
  État intermédiaire en itération courante : [ 2.0112487  -0.99745509  1.21485091]
  État intermédiaire en itération courante : [ 2.00863696 -0.96943284  1.30917045]
  État intermédiaire en itération courante : [ 2.00453385 -0.94011716  1.51021882]
  État intermédiaire en itération courante : [ 2.00013539 -0.93313893  1.80539433]
  État intermédiaire en itération courante : [ 1.95437244 -0.76890418  2.04566856]
  État intermédiaire en itération courante : [ 1.99797362 -0.92538074  1.81674451]
  État intermédiaire en itération courante : [ 1.99760514 -0.95929669  2.01402091]
  État intermédiaire en itération courante : [ 1.99917565 -0.99152672  2.03171791]
  État intermédiaire en itération courante : [ 1.99990376 -0.99963123  2.00671578]
  État intermédiaire en itération courante : [ 1.99999841 -1.00005285  2.00039699]
  État intermédiaire en itération courante : [ 2.00000014 -1.00000307  2.00000221]
  État intermédiaire en itération courante : [ 2.         -0.99999992  1.99999987]

Calage de 3 coefficients pour une forme quadratique 1D sur 5 mesures
--------------------------------------------------------------------

Vecteur d'observation.............: [ 57.   2.   3.  17. 192.]
État d'ébauche a priori...........: [1. 1. 1.]

Coefficients théoriques attendus..: [ 2 -1  2]

Nombre d'itérations...............: 25
Nombre de simulations.............: 100
Écart Observation-Analyse (OMA)...: [ 5.82916385e-07  1.30744229e-07  5.87730300e-08 -6.67061357e-08
 -3.12019267e-07]
Coefficients résultants du calage.: [ 2.         -0.99999992  1.99999987]

Les graphiques illustrant le résultat de son exécution sont les suivants :

Sur la figure, les traits permettent simplement de joindre les valeurs simulées aux points de contrôle et d’en faciliter la lisibilité.

On constate, dans ce cas particulièrement simple, que l’optimum est identique aux valeurs de coefficients [2, -1, 2] utilisés pour générer les pseudo-observations de l’expérience jumelle. Naturellement, la simulation finale obtenue avec les coefficients recalés coïncide donc avec les observations aux 5 points de contrôle.

Pour illustrer plus clairement le recalage effectué, en utilisant les mêmes informations, on peut tracer la version « continue » des courbes quadratiques dont on dispose de l’état d’ébauche pour les coefficients, et de l’état optimal pour ces mêmes coefficients, en même temps que les valeurs obtenues aux points de contrôle. De plus, pour améliorer la lisibilité, on a retiré les traits qui joignent les simulations dans le graphique précédent.

13.1.5.2. Deuxième exemple¶

Le 3DVAR peut aussi être utilisé pour une analyse temporelle des observations d’un modèle dynamique donné. Dans ce cas, l’analyse est conduite de manière itérative, lors de l’arrivée de chaque observation. On utilise pour cet exemple le même système dynamique simple [Welch06] que celui qui est analysé dans les Exemples d’utilisation en Python (TUI) du Filtre de Kalman. Pour une bonne compréhension de la gestion du temps, on se reportera au Schéma temporel d’action des opérateurs pour l’assimilation de données en dynamique et aux explications décrites dans la section pour Approfondir l’assimilation de données pour la dynamique.

A chaque étape, l’analyse 3DVAR classique remet à jour uniquement l’état du système. Moyennant une modification des valeurs de covariances a priori par rapport aux hypothèses initiales du filtrage, cette réanalyse 3DVAR permet de converger vers la trajectoire vraie, comme l’illustre la figure associée, de manière ici un peu plus lente qu’avec un Filtre de Kalman.

Note

Remarque concernant les covariances a posteriori : classiquement, l’analyse itérative 3DVAR remet à jour uniquement l’état et non pas sa covariance. Comme les hypothèses d’opérateurs et de covariance a priori restent inchangées ici au cours de l’évolution, la covariance a posteriori est constante. Le tracé de cette covariance a posteriori, sur la seconde figure qui suit, permet d’insister sur cette propriété tout à fait attendue de l’analyse 3DVAR. Une hypothèse plus évoluée est proposée dans l’exemple d’après.

# -*- coding: utf-8 -*-
#
from numpy import array, random
random.seed(1234567)
Xtrue = -0.37727
Yobs = []
for i in range(51):
    Yobs.append([random.normal(Xtrue, 0.1, size=(1,)),])
#
print("Estimation variationnelle d'une trajectoire temporelle")
print("------------------------------------------------------")
print("  Observations bruitées acquises sur %i pas de temps"%(len(Yobs)-1,))
print("")
from adao import adaoBuilder
case = adaoBuilder.New()
#
case.setBackground         (Vector             = [0.])
case.setBackgroundError    (ScalarSparseMatrix = 0.1**2)
#
case.setObservationOperator(Matrix             = [1.])
case.setObservation        (VectorSerie        = Yobs)
case.setObservationError   (ScalarSparseMatrix = 0.3**2)
#
case.setEvolutionModel     (Matrix             = [1.])
case.setEvolutionError     (ScalarSparseMatrix = 1e-5)
#
case.setAlgorithmParameters(
    Algorithm="3DVAR",
    Parameters={
        "EstimationOf":"State",
        "StoreSupplementaryCalculations":[
            "Analysis",
            "APosterioriCovariance",
            ],
        },
    )
#
case.execute()
Xa = case.get("Analysis")
Pa = case.get("APosterioriCovariance")
#
print("  État analysé à l'observation finale :", Xa[-1])
print("")
print("  Variance a posteriori finale :", Pa[-1])
print("")
#
#-------------------------------------------------------------------------------
#
Observations = array([yo[0]   for yo in Yobs])
Estimates    = array([xa[0]   for xa in case.get("Analysis")])
Variances    = array([pa[0,0] for pa in case.get("APosterioriCovariance")])
#
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (10, 4)
#
plt.figure()
plt.plot(Observations,"kx",label="Mesures bruitées")
plt.plot(Estimates,"r-",label="État estimé")
plt.axhline(Xtrue,color="b",label="Valeur vraie")
plt.legend()
plt.title("Estimation de l'état", fontweight="bold")
plt.xlabel("Pas d'observation")
plt.ylabel("Tension")
plt.savefig("simple_3DVAR2_state.png")
#
plt.figure()
iobs = range(1,len(Observations))
plt.plot(iobs,Variances[iobs],label="Variance d'erreur a posteriori")
plt.title("Estimation de la variance d'erreur a posteriori", fontweight="bold")
plt.xlabel("Pas d'observation")
plt.ylabel("$(Tension)^2$")
plt.setp(plt.gca(),"ylim",[0,.01])
plt.savefig("simple_3DVAR2_variance.png")

Le résultat de son exécution est le suivant :

Estimation variationnelle d'une trajectoire temporelle
------------------------------------------------------
  Observations bruitées acquises sur 50 pas de temps

  État analysé à l'observation finale : [-0.37110687]

  Variance a posteriori finale : [[0.009]]

Les graphiques illustrant le résultat de son exécution sont les suivants :

13.1.5.3. Troisième exemple¶

A partir de l’exemple précédent, si l’on veut adapter la convergence temporelle du 3DVAR, on peut modifier par exemple les hypothèses de covariance a priori des erreurs d’ébauche au cours des itérations. Cette remise à jour est une hypothèse de l’utilisateur, et il y a de multiples possibilités qui vont dépendre de la physique du cas. On en illustre une ici.

On choisit, arbitrairement, de faire décroître la covariance a priori des erreurs d’ébauche d’un facteur constant $0.9^2=0.81$ tant qu’elle reste supérieure à une valeur limite de $0.1^2=0.01$ (qui est la valeur fixe de covariance a priori des erreurs d’ébauche de l’exemple précédent), sachant qu’elle commence à la valeur 1 (qui est la valeur fixe de covariance a priori des erreurs d’ébauche utilisée pour le premier pas du filtrage de Kalman). Cette valeur est remise à jour à chaque pas, en la réinjectant comme covariance a priori de l’état qui est utilisé comme ébauche lors du pas suivant d’analyse, dans une boucle explicite.

On constate dans ce cas que l’estimation d’état converge plus vite vers la valeur vraie, et que l’assimilation se comporte ensuite de manière similaire aux Exemples d’utilisation en Python (TUI) pour le Filtre de Kalman, ou à l’exemple précédent avec les covariances a priori manuellement adaptée. De plus, la covariance a posteriori décroît tant que l’on force la décroissance de la covariance a priori.

Note

On insiste sur le fait que les variations de covariance a priori, qui conditionnent les variations de covariance a posteriori, relèvent d’une hypothèse arbitraire de l’utilisateur et non pas d’une obligation. Cette hypothèse doit donc être adaptée en fonction du cas physique.

# -*- coding: utf-8 -*-
#
from numpy import array, random
random.seed(1234567)
Xtrue = -0.37727
Yobs = []
for i in range(51):
    Yobs.append([random.normal(Xtrue, 0.1, size=(1,)),])
#
print("Estimation variationnelle d'une trajectoire temporelle")
print("------------------------------------------------------")
print("  Observations bruitées acquises sur %i pas de temps"%(len(Yobs)-1,))
print("")
from adao import adaoBuilder
case = adaoBuilder.New()
#
case.setBackground         (Vector             = [0.])
#
case.setObservationOperator(Matrix             = [1.])
case.setObservationError   (ScalarSparseMatrix = 0.3**2)
#
case.setEvolutionModel     (Matrix             = [1.])
case.setEvolutionError     (ScalarSparseMatrix = 1e-5)
#
case.setAlgorithmParameters(
    Algorithm="3DVAR",
    Parameters={
        "EstimationOf":"State",
        "StoreSupplementaryCalculations":[
            "Analysis",
            "APosterioriCovariance",
            ],
        },
    )
#
# Boucle pour obtenir une analyse à l'arrivée de chaque observation
#
Ca = 1.
for i in range(1,len(Yobs)):
    case.setObservation(Vector = Yobs[i])
    case.setBackgroundError(ScalarSparseMatrix = Ca) # Covariance remise à jour
    case.execute( nextStep = True )
    #
    # Hypothèse de décroissance de la covariance a priori de l'erreur d'ébauche
    #
    if float(Ca) <= 0.1**2:
        Ca = 0.1**2
    else:
        Ca = Ca * 0.9**2
#
Xa = case.get("Analysis")
Pa = case.get("APosterioriCovariance")
#
print("  État analysé à l'observation finale :", Xa[-1])
print("")
print("  Variance a posteriori finale :", Pa[-1])
print("")
#
#-------------------------------------------------------------------------------
#
Observations = array([yo[0]   for yo in Yobs])
Estimates    = array([xa[0]   for xa in case.get("Analysis")])
Variances    = array([pa[0,0] for pa in case.get("APosterioriCovariance")])
#
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (10, 4)
#
plt.figure()
plt.plot(Observations,"kx",label="Mesures bruitées")
plt.plot(Estimates,"r-",label="État estimé")
plt.axhline(Xtrue,color="b",label="Valeur vraie")
plt.legend()
plt.title("Estimation de l'état", fontweight="bold")
plt.xlabel("Pas d'observation")
plt.ylabel("Tension")
plt.savefig("simple_3DVAR3_state.png")
#
plt.figure()
iobs = range(1,len(Observations))
plt.plot(iobs,Variances[iobs],label="Variance d'erreur a posteriori")
plt.title("Estimation de la variance d'erreur a posteriori", fontweight="bold")
plt.xlabel("Pas d'observation")
plt.ylabel("$(Tension)^2$")
plt.savefig("simple_3DVAR3_variance.png")

Le résultat de son exécution est le suivant :

Estimation variationnelle d'une trajectoire temporelle
------------------------------------------------------
  Observations bruitées acquises sur 50 pas de temps

  État analysé à l'observation finale : [-0.37334336]

  Variance a posteriori finale : [[0.009]]

Les graphiques illustrant le résultat de son exécution sont les suivants :

13.1.5.4. Correction de l’état d’un système dynamique de Lorenz par 3DVAR¶

Cet exemple prolonge les cas simples d’usage de l’algorithme 3DVAR (voir la partie Exemples avec l’algorithme de « 3DVAR », ou un exemple similaire dans [Asch16]). Il décrit l’effet d’une assimilation de données sur la prévision temporelle d’un système dynamique. C’est un cas très simplifié qui illustre une démarche classique en météorologie (prévision du temps à court terme, en-deçà de la quinzaine de jours) ou en prévision saisonnière (prévision du temps à plus long terme, au-delà de la quinzaine de jours).

On utilise ici un modèle simple classique de système dynamique, nommé système de Lorenz, oscillateur de Lorenz, Lorenz3D ou Lorenz63 d’après le nom de son auteur, le mathématicien et météorologue Edward Lorenz ([Lorenz63], [WikipediaL63]). Il est disponible dans les modèles de tests intégrés d’ADAO sous le nom Lorenz1963.

Ce système dynamique tridimensionnel non-linéaire est un modèle extrêmement simplifié des équations de Navier-Stokes, en vue d’étudier le couplage de l’atmosphère et de l’océan, dans une configuration physique particulière de convection (convection de Rayleigh-Bénard). Il présente un comportement de chaos déterministe, ce qui signifie que le système dynamique est déterministe (du point de vue continu ou discret), mais que sa sensibilité aux conditions initiales ou aux paramètres permet d’obtenir au bout d’un temps fini une trajectoire arbitrairement différente pour de petites variations de ses conditions initiales. Ce système est connu pour être à l’origine de la notion d’effet papillon [Butterfly72] et pour l’illustration qu’il en donne.

Ce système est dépendant de ses conditions initiales et de 3 paramètres physiques $(\sigma,\rho,\beta)$ . Les équations classiques pour l’état $u=(x,y,z)$ le décrivant sont :

$\left\{ \begin{array}{lcl} \displaystyle\frac{dx}{dt} & = & \sigma (y - x)\\~\\ \displaystyle\frac{dy}{dt} & = & \rho\, x - y - x\, z\\~\\ \displaystyle\frac{dz}{dt} & = & x\, y - \beta\, z \end{array} \right.$

avec le temps $t\in{I\!R}^+$ , et avec $\sigma=10$ , $\rho=28$ , $\beta=8/3$ des valeurs courantes de paramètres caractérisant un état chaotique. La condition initiale est quelconque.

La figure suivante illustre l’attracteur de Lorenz par la simulation directe de ce système dynamique, que l’on représente ici pour $t\in[0,40]$ et avec la condition initiale $(0,1,0)$ . L’attracteur est la structure qui correspond au comportement long terme de l’oscillateur de Lorenz, qui se présente donc ici comme deux « ailes de papillon ». La figure montre que la variable d’état $u=(x,y,z)$ de ce système dynamique évolue sur une trajectoire déterministe non périodique, qui s’enroule sur une aile du papillon pour ensuite « sauter » sur l’autre aile, et ainsi de suite de façon en apparence erratique. La figure ci-dessous présente une unique trajectoire, colorée de manière différente pour sa première moitié et sa seconde moitié pour illustrer les enroulements et sauts successifs.

Si on observe la même trajectoire en traçant séparément les 3 composantes de l’état du système (sur une durée plus courte $[0,10]$ ), avec les mêmes paramètres de simulation, on obtient la figure suivante :

En remarquant bien que les amplitudes et temporalités des deux premières composantes $x$ et $y$ ne sont effectivement pas identiques, malgré leur ressemblance, cette figure illustre la propriété de chaos par une absence de régularité spatiale et temporelle.

Pour illustrer la démarche d’analyse, on se place pour la suite en expériences jumelles (voir la partie Pour tester une chaîne d’assimilation de données : les expériences jumelles).

On choisit ici de perturber uniquement l’état initial de la simulation, en utilisant la valeur d’ébauche $u^b=[2,3,4]$ , dont les effets sont à comparer à ceux de l’état initial, dit idéal ou vrai, non perturbé, valant $u^t=[1,1,1]$ .

Le modèle est considéré comme parfait, et il est observé sur l’intervalle temporel [0,2]. Les pseudo-observations $\mathbf{y}^o$ sont au nombre de 10, construites par échantillonnage au pas de temps $\delta t=0.2$ sur cet intervalle temporel, sur la base d’une simulation depuis l’état initial non perturbé $u^t=[1,1,1]$ . A ces valeurs exactes échantillonnées est ensuite ajouté, sur chaque composante, un bruit gaussien d’amplitude $\sigma_m=0.15$ pour obtenir une pseudo-observation. L’ordre de grandeur du bruit est celui du bruit expérimental de mesure de variables réelles correspondant au modèle simplifié de Lorenz.

La figure suivante illustre sur la première variable (les autres sont similaires) ces différentes informations réparties sur l’intervalle temporel [0,2] de mesure :

Note

On insiste fortement sur le fait que les observations successives ne sont pas disponibles simultanément, à l’instant initial par exemple, mais qu’elles sont disponibles à chaque fois que la simulation atteint l’un des instants de mesure. On rappelle de plus que la courbe de simulation bleue en pointillés, obtenue à partir de l’état idéal ou vrai $u^t=[1,1,1]$ , est inconnue en dehors d’expériences jumelles. Seules les trajectoires tracées comme continues sont connues par simulation.

Sous forme numérique, enregistrées dans un fichier nommé simple_3DVAR4Observations.csv pour relecture ultérieure, les valeurs d’observations sont les suivantes :

# Observations Lorenz1963 (format: t x y z)
# 
20 +6.6255693323143952e+00 +1.3512204021575638e+01 +3.9860034533510582e+00
40 +1.5140047444332868e+01 +1.3495388075296226e+00 +4.6611660816669115e+01
60 -4.7609818075165542e+00 -7.9674154050023951e+00 +2.6781548199549242e+01
80 -8.4989430323515158e+00 -1.0087805440794597e+01 +2.5436627475898202e+01
00 -9.4032667798914851e+00 -8.4574726949541308e+00 +2.9469200645262447e+01
20 -7.0450610945351828e+00 -6.7568835708764148e+00 +2.6136357537101407e+01
40 -8.4087739643843644e+00 -9.7426531698975598e+00 +2.5181746694435667e+01
60 -9.4866303357699397e+00 -8.8142676152554458e+00 +2.9449435747857493e+01
80 -6.9109551287087747e+00 -6.3753292483483364e+00 +2.6504633755059768e+01
00 -8.0717069604712552e+00 -9.7394819137223507e+00 +2.4481237978123101e+01

L’assimilation de données est ensuite utilisée pour corriger la trajectoire initiale d’ébauche du système à chaque nouvelle observation acquise. A chaque étape $n$ , l’état courant prévu $\mathbf{u}^f_n$ est modifié pour produire un nouvel état analysé $\mathbf{u}^a_n$ , tenant compte de cette nouvelle information $\mathbf{y}^o_n$ , puis la simulation se poursuit à partir de ce nouvel état. L’assimilation des observations se produit grâce au simple script ADAO ci-dessous. Pour une meilleure lisibilité, le script présente de manière explicite la boucle temporelle sur les observations :

# -*- coding: utf-8 -*-
#
import numpy
from adao import adaoBuilder
from Models.Lorenz1963 import EDS as Lorenz1963
numpy.set_printoptions(precision=5)
#
u0Back = numpy.array([2, 3, 4])    # Ébauche (différente de l'état vrai)
sigma_m = 0.15                     # Écart-type du bruit de mesure
sigma_b = 0.1                      # Écart-type du bruit d'ébauche
H = numpy.eye(3)                   # Opérateur d'observation
ODE = Lorenz1963(dt = 0.01)        # Modèle dynamique
ODE.ObservationStep = 0.2          # Intervalle d'observation
#
observations = numpy.loadtxt("simple_3DVAR4Observations.csv")[:, 1:]
#
case = adaoBuilder.New()
case.setAlgorithmParameters(
    Algorithm="3DVAR",
    Parameters={"EstimationOf": "State"},
)
case.setBackground(Vector=u0Back)
case.setBackgroundError(ScalarSparseMatrix=sigma_b**2)
case.setObservationError(ScalarSparseMatrix=sigma_m**2)
case.setEvolutionError(ScalarSparseMatrix=1.0e-8)
case.setObservationOperator(Matrix=H)
#
print("Analyses successives corrigeant l'état prévu en utilisant l'observation :")
for i in range(observations.shape[0]):
    case.setEvolutionModel(OneFunction=ODE.StateTransition)
    case.setObservation(Vector=observations[i, :])
    case.execute(nextStep=True)
    #
    Xa = case.get("Analysis")[-1]
    print(
        "  Après l'observation no %02i : Xa[%02i] = %+9.5f %+9.5f %+9.5f"
        % (i + 1, i + 1, Xa[0], Xa[1], Xa[2])
    )

Le résultat de son exécution est le suivant :

Analyses successives corrigeant l'état prévu en utilisant l'observation :
  Après l'observation no 01 : Xa[01] = +10.81803 +20.13078 +12.79257
  Après l'observation no 02 : Xa[02] = +10.62741  -3.02604 +41.26296
  Après l'observation no 03 : Xa[03] =  -4.28903  -6.99542 +24.84772
  Après l'observation no 04 : Xa[04] =  -8.76412 -10.93891 +24.68112
  Après l'observation no 05 : Xa[05] =  -9.70093  -8.19724 +30.32881
  Après l'observation no 06 : Xa[06] =  -6.73955  -6.29483 +25.70542
  Après l'observation no 07 : Xa[07] =  -8.38183  -9.99790 +24.60690
  Après l'observation no 08 : Xa[08] =  -9.76835  -8.91467 +29.73469
  Après l'observation no 09 : Xa[09] =  -7.01017  -6.31548 +26.40657
  Après l'observation no 10 : Xa[10] =  -8.05253  -9.61682 +24.32317

La figure suivante illustre alors, sur la première variable, les états prévus, mesurés et analysés à chaque pas d’observation, ainsi que les états simulés par le modèle entre deux observations. Les ruptures temporelles de la trajectoire proviennent par nature de chaque opération d’assimilation avec une nouvelle information observée.

On constate clairement que l’état prévu $\mathbf{u}^f_n$ est de plus en plus conforme à chaque pseudo-observation $\mathbf{y}^o_n$ , et que chaque portion de trajectoire prévue est de plus en plus proche de la trajectoire idéale inconnue issue de l’état vrai $\mathbf{u}^t$ . Dans la figure suivante, la mesure classique des écarts par RMSE (Root-Mean-Square Error) décrit de manière similaire cette amélioration itérative des prévisions et des analyses intégrant l’information d’observations.

On peut aussi comparer les prévisions d’état obtenues sur l’intervalle temporel [2,10], qui se situe au-delà de la fenêtre d’observation concernée par la démarche d’assimilation de données. On dispose de trois manières d’établir cette simulation temporelle après l’instant $t=2$ :

on peut calculer la prévision issue de l’état à $t=2$ obtenu lui-même par simulation à partir de l’état initial perturbé d’ébauche $u^b=[2,3,4]$ à l’instant $t=0$ , prévision dans laquelle on n’utilise que l’information d’ébauche ;
on peut calculer la prévision issue de l’état à $t=2$ obtenu lui-même par simulation à partir de l’état initial non perturbé (état idéal ou état vrai) $u^t=[1,1,1]$ , qui est considéré comme la référence mais connue uniquement à travers les observations issues de l’échantillonnage avec bruit ;
on peut calculer la prévision issue de l’analyse par assimilation de données, qui provient de la correction d’état par les observations.

Ces trois manières de faire sont illustrées sur la figure suivante, où l’on a rappelé la fenêtre d’assimilation sur l’intervalle temporel [0,2] et où les prévisions sont établies sur l’intervalle temporel [2,10] :

Il est manifeste que la prévision issue de l’analyse $\mathbf{u}^a$ , obtenue par assimilation de données à l’instant $t=2$ , est beaucoup plus conforme à la simulation $\mathbf{u}^t$ idéale choisie pour l’expérience jumelle, que la prévision issue de l’ébauche $\mathbf{u}^b$ simulée à $t=2$ . Même si la prolongation temporelle de la simulation conduit obligatoirement à un écart avec l’état vrai qui s’accroît, à cause de la propriété chaotique du système de Lorenz, la correction par assimilation de données permet d’obtenir une prévision acceptable sur une durée beaucoup plus longue que l’absence de correction intrinsèquement présente dans la simulation issue de l’ébauche.

13.1.6. Voir aussi¶

Références vers d’autres sections :

Références bibliographiques :