13.7. Algorithme de calcul « EnsembleKalmanFilter »¶

13.7.1. Description¶

Cet algorithme réalise une estimation de l’état d’un système dynamique par un filtre de Kalman d’ensemble (EnKF), permettant d’éviter de devoir calculer les opérateurs tangent ou adjoint pour les opérateurs d’observation ou d’évolution, comme dans les filtres de Kalman simple ou étendu.

Il s’applique aux cas d’opérateurs d’observation et d’évolution incrémentale (processus) non-linéaires et présente d’excellentes qualités de robustesse et de performances. Il peut être interprété comme une réduction d’ordre du filtre de Kalman classique, avec une remarquable qualité d’assimilation de ce filtrage pour les problèmes de grande taille. Il peut être rapproché d’un Algorithme de calcul « UnscentedKalmanFilter » dont les qualités sont similaires pour les systèmes non-linéaires.

On remarque qu’il n’y a pas d’analyse effectuée au pas de temps initial (numéroté 0 dans l’indexage temporel) car il n’y a pas de prévision à cet instant (l’ébauche est stockée comme pseudo-analyse au pas initial). Si les observations sont fournies en série par l’utilisateur, la première n’est donc pas utilisée. 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.

Dans le cas d’opérateurs linéaires ou « faiblement » non-linéaire, on peut aisément utiliser l”Algorithme de calcul « ExtendedKalmanFilter » ou même l”Algorithme de calcul « KalmanFilter », qui sont souvent largement moins coûteux en évaluations sur de petits systèmes. On peut vérifier la linéarité des opérateurs à l’aide de l”Algorithme de vérification « LinearityTest ».

Il existe de nombreuses variantes déterministes ou stochastiques de cet algorithme, permettant en particulier d’effectuer de la réduction de taille des problèmes algébriques à différents niveaux (en utilisant des méthodes de rang réduit, de la réduction de dimension, des changements d’espace de calcul, conduisant à des schémas de type Ensemble Square Root Kalman Filters (EnSRKF) ou Reduced-Rank Square Root Filters (RRSQRT), à des transformations déterministes…). On ne rentre pas ici dans le détail complexe des classifications et des équivalences algorithmiques, qui sont disponibles dans la littérature. On propose ici les formulations stables et robustes suivantes :

« EnKF » (Ensemble Kalman Filter, voir [Evensen94]), algorithme stochastique original, permettant de traiter de manière consistante un opérateur d’évolution non-linéaire,
« ETKF » (Ensemble-Transform Kalman Filter), algorithme déterministe d’EnKF, permettant de traiter un opérateur d’évolution non-linéaire avec beaucoup moins de membres (on recommande d’utiliser un nombre de membres de l’ordre de 10 ou même parfois moins),
« ETKF-N » (Ensemble-Transform Kalman Filter of finite size N), algorithme d’ETKF dit de « taille finie N », évitant de recourir à une inflation souvent nécessaire avec les autres algorithmes,
« MLEF » (Maximum Likelihood Kalman Filter, voir [Zupanski05]), algorithme déterministe d’EnKF, permettant en plus de traiter de manière consistante un opérateur d’observation non-linéaire,
« IEnKF » (Iterative EnKF), algorithme déterministe d’EnKF, améliorant le traitement des non-linéarités des opérateurs,
« E3DVAR » (EnKF 3DVAR, ou 3D-Var-Ben), algorithme couplant assimilation d’ensemble et variationnelle, qui utilise en parallèle une assimilation variationnelle 3DVAR pour l’estimation d’un unique meilleur état et un algorithme d’ensemble EnKF pour améliorer l’estimation des covariances d’erreurs a posteriori,
« EnKS » (Ensemble Kalman Smoother), algorithme de lissage avec un décalage temporel fixe de taille L.

Sans pouvoir prétendre à l’universalité, on recommande d’utiliser la formulation « EnKF » comme référence, la formulation « ETKF-N » ou « IEnKF » pour une performance robuste, et les autres algorithmes (dans l’ordre) comme des moyens pour obtenir une assimilation de données plus économique et de qualité (éventuellement) similaire.

13.7.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 locale du minimum, permettant en théorie d’atteindre un état localement optimal (par opposition à un état « globalement optimal »).

Les méthodes proposées par cet algorithme ne requièrent pas de dérivation de la fonction objectif ou de l’un des opérateurs, permettant d’éviter ce coût 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 présentent un parallélisme interne, et peuvent donc profiter de ressources informatiques de répartition de calculs. L’interaction potentielle, entre le parallélisme interne des méthodes, 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 statiques, fixés par des propriétés algorithmiques particulières. En pratique, il peut y avoir plusieurs critères de convergence actifs simultanément.

La propriété algorithmique la plus courante est celle des calculs directs, qui évaluent la solution à convergence sans itération contrôlable. Il n’y a aucun seuil de convergence à régler dans ce cas.

13.7.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 :

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 « State ».

Exemple : {"EstimationOf":"State"}

HybridCostDecrementTolerance

Valeur réelle. Cette clé indique une valeur limite, conduisant à arrêter le processus itératif d’optimisation dans la partie variationnelle du couplage, lorsque la fonction coût décroît moins que cette tolérance au dernier pas. Le 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 : {"HybridCostDecrementTolerance":1.e-7}

HybridCovarianceEquilibrium

Valeur réelle. Cette clé indique, en optimisation hybride variationnelle, le facteur d’équilibre entre la covariance statique a priori et la covariance d’ensemble. Ce facteur est compris entre 0 et 1, et sa valeur par défaut est 0.5.

Exemple : {"HybridCovarianceEquilibrium":0.5}

HybridMaximumNumberOfIterations

Valeur entière. Cette clé indique le nombre maximum d’itérations internes possibles en optimisation hybride, pour la partie variationnelle. Le défaut est 15000, qui est très similaire à une absence de limite sur les 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 : {"HybridMaximumNumberOfIterations":100}

InflationFactor

Valeur réelle. Cette clé indique le facteur d’inflation dans les méthodes d’ensemble, à appliquer sur la covariance ou les anomalies selon le choix du type d’inflation. Sa valeur doit être positive si l’inflation est additive, ou supérieure à 1 si l’inflation est multiplicative. La valeur par défaut est 1, qui conduit à une absence d’inflation multiplicative. L’absence d’inflation additive est obtenue en indiquant une valeur de 0.

Exemple : {"InflationFactor":1.}

InflationType

Nom prédéfini. Cette clé permet d’indiquer la méthode d’inflation dans les méthodes d’ensemble, pour celles qui nécessitent une telle technique. L’inflation peut être appliquée de diverses manières, selon les options suivantes : multiplicative ou additive du facteur d’inflation spécifié, appliquée sur l’ébauche ou sur l’analyse, appliquée sur les covariances ou sur les anomalies. Une inflation multiplicative sur les anomalies, qui sont obtenues en retranchant la moyenne d’ensemble, est effectuée en multipliant ces anomalies par le facteur d’inflation, puis en reconstituant les membres de l’ensemble par ajout de la moyenne préalablement calculée. Un seul type d’inflation est appliqué à la fois, et la valeur par défaut est « MultiplicativeOnAnalysisAnomalies ». Les noms possibles sont dans la liste suivante : [ « MultiplicativeOnAnalysisAnomalies », « MultiplicativeOnBackgroundAnomalies », ].

Exemple : {"InflationType":"MultiplicativeOnAnalysisAnomalies"}

NumberOfMembers

Valeur entière. Cette clé indique le nombre de membres utilisés pour réaliser la méthode d’ensemble. Le défaut est de 100, et il est recommandé de l’adapter aux besoins pour des problèmes réels.

Exemple : {"NumberOfMembers":100}

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}

SmootherLagL

Valeur entière. Cette clé indique le nombre d’intervalles de temps de lissage dans le passé pour l’EnKS. C’est bien un nombre d’intervalles, et non pas une durée fixe. La valeur par défaut est 0, qui conduit à une absence de lissage.

Exemple : {"SmootherLagL":0}

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 », « ForecastCovariance », « ForecastState », « IndexOfOptimum », « InnovationAtCurrentAnalysis », « InnovationAtCurrentState », « SimulatedObservationAtCurrentAnalysis », « SimulatedObservationAtCurrentOptimum », « SimulatedObservationAtCurrentState », ].

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 la formulation « EnKF » d’origine, et les choix possibles sont « EnKF » (Ensemble Kalman Filter), « ETKF » (Ensemble-Transform Kalman Filter), « ETKF-N » (Ensemble-Transform Kalman Filter), « MLEF » (Maximum Likelihood Kalman Filter), « IEnKF » (Iterative_EnKF), « E3DVAR » (EnKF 3DVAR), « EnKS » (Ensemble Kalman Smoother).

Il est conseillé d’essayer les variantes « ETKF-N » ou « IEnKF » pour une performance robuste, et de réduire le nombre de membres à une dizaine ou moins pour toutes les variantes autres que la formulation « EnKF » originale.

Exemple : {"Variant":"EnKF"}

13.7.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, de la variable « ADD » du post-processing en interface graphique, ou du 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 l”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]

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")[:]

ForecastCovariance

Liste de matrices. Chaque élément est une matrice de covariance d’erreur sur l’état prévu par le modèle au cours du déroulement itératif temporel de l’algorithme utilisé.

Exemple : pf = ADD.get("ForecastCovariance")[-1]

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]

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]

SimulatedObservationAtCurrentAnalysis

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. Cette quantité est identique au vecteur d’observation simulé à l’état courant dans le cas d’une assimilation mono-état.

Exemple : hxs = ADD.get("SimulatedObservationAtCurrentAnalysis")[-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]

13.7.5. Voir aussi¶

Références vers d’autres sections :

Références bibliographiques :