13.3. Algorithme de calcul « Blue »

Description

Cet algorithme réalise une estimation de type BLUE (Best Linear Unbiased Estimator) de l’état d’un système. C’est une estimation linéaire, sans biais et optimale. De manière technique, c’est ici un estimateur d’Aitken. Il réalise la meilleure estimation linéaire de l’état à l’aide de l’état d’ébauche initial et des observations. Il est théoriquement réservé aux cas d’opérateurs d’observation linéaires, même s’il fonctionne parfois dans les cas « faiblement » non-linéaires. On peut vérifier la linéarité de l’opérateur d’observation à l’aide d’un Algorithme de vérification « LinearityTest ». Cet algorithme est toujours le plus rapide de l’ensemble des algorithmes d’assimilation d’ADAO.

Cet algorithme 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 ADAO, il peut aussi être utilisé sur une succession d’observations, plaçant alors l’estimation dans un cadre récursif en partie similaire à un Algorithme de calcul « KalmanFilter ». Une estimation standard 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.

En cas de non-linéarité, même peu marquée, on lui préférera aisément un Algorithme de calcul « ExtendedBlue » ou un Algorithme de calcul « 3DVAR ».

Remarque complémentaire : une simplification algébrique du BLUE conduit à la méthode d’interpolation dite optimale nommée « Optimal Interpolation » ou « OI ». C’est une méthode très simple et peu coûteuse, spécialement adaptée aux problèmes de très (très) grande taille, mais dont l’inconvénient est de fournir un résultat d’analyse globalement sous-optimal et bruité, voire incohérent. Le moyen d’éviter ces désavantages est d’adapter très précisément les éléments de la méthode à chaque modèle physique. Pour ces raisons, cette méthode n’est donc pas proposée.

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

Exemple : {"EstimationOf":"Parameters"}

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}

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 conseiller 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 », « CurrentOptimum », « CurrentState », « CurrentStepNumber », « ForecastState », « Innovation », « InnovationAtCurrentAnalysis », « MahalanobisConsistency », « OMA », « OMB », « SampledStateForQuantiles », « SigmaBck2 », « SigmaObs2 », « SimulatedObservationAtBackground », « SimulatedObservationAtCurrentOptimum », « SimulatedObservationAtCurrentState », « SimulatedObservationAtOptimum », « SimulationQuantiles », ].

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

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 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 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 : C = 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 : A = 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 : S = 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 : V = 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")[:]

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 : i = ADD.get("CurrentStepNumber")[-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")[:]

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 : ds = ADD.get("InnovationAtCurrentAnalysis")[-1]

MahalanobisConsistency

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

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

SigmaBck2

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

Exemple : sb2 = ADD.get("SigmaBck")[-1]

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

Exemples d’utilisation en Python (TUI)

Voici un exemple très simple d’usage de l’algorithme proposé et de ses paramètres, écrit en [DocR] Interface textuelle pour l’utilisateur (TUI/API), et dont les informations indiquées en entrée permettent de définir un cas équivalent en interface graphique.

Cet exemple décrit l’interpolation entre deux états physiques. Ces deux champs vectoriels, de discrétisation identique, sont l’observation \mathbf{y}^o et l’état d’ébauche a priori \mathbf{x}^b. Les confiances dans les erreurs sur les deux informations sont considérées comme identiques. Le modèle H observe complètement le champ disponible, c’est un opérateur de sélection matriciel.

Le champ interpolé résultant est simplement le « milieu » entre les deux champs, avec une confiance améliorée sur les erreurs.

# -*- coding: utf-8 -*-
#
from numpy import array, ravel
from adao import adaoBuilder
case = adaoBuilder.New('')
case.setBackground( Vector = array([0., 1., 2.]), Stored=True )
case.setBackgroundError( ScalarSparseMatrix = 1. )
case.setObservation( Vector = array([10., 11., 12.]), Stored=True )
case.setObservationError( ScalarSparseMatrix = 1. )
case.setObservationOperator( Matrix = array([[1., 0., 0.],
                                             [0., 1., 0.],
                                             [0., 0., 1.]]), )
case.setAlgorithmParameters(
    Algorithm='Blue',
    Parameters={
        'StoreSupplementaryCalculations': [
            'APosterioriCovariance',
            ],
        },
    )
case.execute()
#
#-------------------------------------------------------------------------------
#
print("Interpolation entre deux états vectoriels, observation et ébauche")
print("-----------------------------------------------------------------")
print("")
print("Vecteur d'observation.........:", ravel(case.get('Observation')))
print("État d'ébauche a priori.......:", ravel(case.get('Background')))
print("")
print("État théorique attendu........:", ravel([5., 6., 7.]))
print("")
print("État obtenu par interpolation.:", ravel(case.get('Analysis')[-1]))
print("Covariance a posteriori.......:\n", case.get('APosterioriCovariance')[-1])

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

Interpolation entre deux états vectoriels, observation et ébauche
-----------------------------------------------------------------

Vecteur d'observation.........: [10. 11. 12.]
État d'ébauche a priori.......: [0. 1. 2.]

État théorique attendu........: [5. 6. 7.]

État obtenu par interpolation.: [5. 6. 7.]
Covariance a posteriori.......:
 [[0.5 0.  0. ]
 [0.  0.5 0. ]
 [0.  0.  0.5]]