13.10. Algorithme de calcul « KalmanFilter »¶

13.10.1. Description¶

Cet algorithme réalise une estimation de l’état d’un système dynamique par un filtre de Kalman. Sous forme discrète, c’est un estimateur itératif (ou récursif) de l’état courant à l’aide de l’état précédent et des observations actuelles. Le temps (ou pseudo-temps) entre deux pas est celui qui sépare les observations successives. Chaque pas d’itération est composé de deux étapes successives dites classiquement de « prédiction » puis de « correction ». L’étape de prédiction utilise un opérateur d’évolution incrémentale pour établir une estimation de l’état courant à partir de l’état estimé au pas précédent. L’étape de correction (ou de mise à jour) utilise les observations courantes pour améliorer l’estimation en corrigeant l’état prédit.

Il est théoriquement réservé aux cas d’opérateurs d’observation et d’évolution incrémentale (processus) linéaires, même s’il fonctionne parfois dans les cas « faiblement » non-linéaire. On peut vérifier la linéarité de l’opérateur d’observation à l’aide d’un Algorithme de vérification « LinearityTest ».

Conceptuellement, on peut représenter le schéma temporel d’action des opérateurs d’évolution et d’observation dans cet algorithme de la manière suivante, avec x l’état, P la covariance d’erreur d’état, t le temps itératif discret :

Fig. 13.2 Schéma temporel des étapes en assimilation de données par filtre de Kalman¶

Dans ce schéma, l’analyse (x,P) est obtenue à travers la « correction » par l’observation de la « prévision » de l’état précédent. 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.

Ce filtre peut aussi être utilisé pour estimer (conjointement ou uniquement) des paramètres et non pas l’état, auquel cas ni le temps ni l’évolution n’ont plus de signification. Les pas d’itération sont alors liés à l’insertion d’une nouvelle observation dans l’estimation récursive. On consultera la section Approfondir l’assimilation de données pour la dynamique pour les concepts de mise en oeuvre.

En cas de non-linéarité des opérateurs, même peu marquée, on lui préférera un Algorithme de calcul « ExtendedKalmanFilter », ou un Algorithme de calcul « EnsembleKalmanFilter » et un Algorithme de calcul « UnscentedKalmanFilter » qui sont plus stables, supportent des bornes sur l’état, etc. On peut vérifier la linéarité des opérateurs à l’aide d’un Algorithme de vérification « LinearityTest ».

13.10.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 requièrent la dérivation de la fonction objectif ou de l’un des opérateurs. Il 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 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 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 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.10.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"}

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

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

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

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]

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.10.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.10.5.1. Premier exemple¶

Le filtre de Kalman peut être utilisé pour une réanalyse des observations d’un modèle dynamique donné. C’est parce que l’ensemble de l’historique complet de l’observation est déjà connu au début des fenêtres temporelles qu’on parle de réanalyse, même si l’analyse itérative conserve comme inconnues les observations futures à un pas de temps donné.

Cet exemple décrit l’estimation itérative d’une quantité physique constante (une tension électrique) selon l’exemple de [Welch06] (pages 11 et suivantes, aussi disponible dans le SciPy Cookbook). Ce modèle permet d’illustrer l’excellent comportement de cet algorithme vis-à-vis du bruit de mesure lorsque le modèle d’évolution est simple. Le problème physique est l’estimation d’une tension électrique, observée sur 50 pas de temps, avec du bruit, ce qui implique ensuite 50 étapes d’analyses par le filtre. L’état idéalisé (valeur dite « vraie », inconnu dans un cas réel) est désigné par Xtrue dans l’exemple. Les observations $\mathbf{y}^o$ (désignée par Yobs dans l’exemple) sont à renseigner, en utilisant le mot-clé « VectorSerie », comme une série chronologique de mesures. On privilégie les observations au détriment de l’ébauche, par l’indication d’une importante variance d’erreur d’ébauche par rapport à la variance d’erreur d’observation. La première observation n’est pas utilisée car l’ébauche $\mathbf{x}^b$ sert de première estimation de l’état.

L’estimation s’effectue en affichant des résultats intermédiaires lors du filtrage itératif. Grâce à ces informations intermédiaires, on peut aussi obtenir les graphiques illustrant l’estimation de l’état et de la covariance d’erreur a posteriori associée.

# -*- 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 par filtrage 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 = 1.)
#
case.setObservationOperator(Matrix             = [1.])
case.setObservation        (VectorSerie        = Yobs)
case.setObservationError   (ScalarSparseMatrix = 0.1**2)
#
case.setEvolutionModel     (Matrix             = [1.])
case.setEvolutionError     (ScalarSparseMatrix = 1e-5)
#
case.setAlgorithmParameters(
    Algorithm="KalmanFilter",
    Parameters={
        "StoreSupplementaryCalculations":[
            "Analysis",
            "APosterioriCovariance",
            ],
        },
    )
case.setObserver(
    Info="  État analysé à l'observation courante :",
    Template='ValuePrinter',
    Variable='Analysis',
    )
#
case.execute()
Xa = case.get("Analysis")
Pa = case.get("APosterioriCovariance")
#
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_KalmanFilter1_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_KalmanFilter1_variance.png")

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

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

  État analysé à l'observation courante : [0.]
  État analysé à l'observation courante : [-0.41804504]
  État analysé à l'observation courante : [-0.3114053]
  État analysé à l'observation courante : [-0.31191336]
  État analysé à l'observation courante : [-0.32761493]
  État analysé à l'observation courante : [-0.33597167]
  État analysé à l'observation courante : [-0.35629573]
  État analysé à l'observation courante : [-0.36840289]
  État analysé à l'observation courante : [-0.37392713]
  État analysé à l'observation courante : [-0.36331937]
  État analysé à l'observation courante : [-0.35750362]
  État analysé à l'observation courante : [-0.37963052]
  État analysé à l'observation courante : [-0.37117993]
  État analysé à l'observation courante : [-0.36732985]
  État analysé à l'observation courante : [-0.37148382]
  État analysé à l'observation courante : [-0.36798059]
  État analysé à l'observation courante : [-0.37371077]
  État analysé à l'observation courante : [-0.3661228]
  État analysé à l'observation courante : [-0.36777529]
  État analysé à l'observation courante : [-0.37681677]
  État analysé à l'observation courante : [-0.37007654]
  État analysé à l'observation courante : [-0.37974517]
  État analysé à l'observation courante : [-0.37964703]
  État analysé à l'observation courante : [-0.37514278]
  État analysé à l'observation courante : [-0.38143128]
  État analysé à l'observation courante : [-0.38790654]
  État analysé à l'observation courante : [-0.38880008]
  État analysé à l'observation courante : [-0.38393577]
  État analysé à l'observation courante : [-0.3831028]
  État analysé à l'observation courante : [-0.37680097]
  État analysé à l'observation courante : [-0.37891813]
  État analysé à l'observation courante : [-0.38147782]
  État analysé à l'observation courante : [-0.37981569]
  État analysé à l'observation courante : [-0.38274266]
  État analysé à l'observation courante : [-0.38418507]
  État analysé à l'observation courante : [-0.38923054]
  État analysé à l'observation courante : [-0.38400006]
  État analysé à l'observation courante : [-0.38562502]
  État analysé à l'observation courante : [-0.3840503]
  État analysé à l'observation courante : [-0.38775222]
  État analysé à l'observation courante : [-0.37700787]
  État analysé à l'observation courante : [-0.37328191]
  État analysé à l'observation courante : [-0.38024181]
  État analysé à l'observation courante : [-0.3815806]
  État analysé à l'observation courante : [-0.38392063]
  État analysé à l'observation courante : [-0.38539266]
  État analysé à l'observation courante : [-0.37856929]
  État analysé à l'observation courante : [-0.37744505]
  État analysé à l'observation courante : [-0.37154554]
  État analysé à l'observation courante : [-0.37405773]
  État analysé à l'observation courante : [-0.37725236]

  Variance a posteriori finale : [[0.00033921]]

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

_images/simple_KalmanFilter1_variance.png

13.10.5.2. Deuxième exemple¶

Le filtre de Kalman peut aussi être utilisé pour une analyse courante 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.

L’exemple suivant porte sur le même système dynamique simple issu de la référence [Welch06]. La différence essentielle consiste à effectuer l’exécution d’une étape de Kalman à l’arrivée de chaque observation fournie itérativement. Le mot-clé « nextStep », inclut dans l’ordre d’exécution, permet de ne pas stocker l’ébauche en double de l’analyse précédente.

De manière entièrement similaire à la réanalyse (sachant que l’on peut afficher des résultats intermédiaires qui sont ici omis par simplicité), l’estimation donne les mêmes résultats lors du filtrage itératif. Grâce à ces informations intermédiaires, on peut aussi obtenir les graphiques illustrant l’estimation de l’état et de la covariance d’erreur a posteriori associée.

# -*- 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 par filtrage 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 = 1.)
#
case.setObservationOperator(Matrix             = [1.])
case.setObservationError   (ScalarSparseMatrix = 0.1**2)
#
case.setEvolutionModel     (Matrix             = [1.])
case.setEvolutionError     (ScalarSparseMatrix = 1e-5)
#
case.setAlgorithmParameters(
    Algorithm="KalmanFilter",
    Parameters={
        "StoreSupplementaryCalculations":[
            "Analysis",
            "APosterioriCovariance",
            ],
        },
    )
#
# Boucle pour obtenir une analyse à l'arrivée de chaque observation
#
for i in range(1,len(Yobs)):
    case.setObservation(Vector = Yobs[i])
    case.execute( nextStep = True )
#
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_KalmanFilter2_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_KalmanFilter2_variance.png")

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

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

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

  Variance a posteriori finale : [[0.00033921]]

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

_images/simple_KalmanFilter2_variance.png

13.10.6. Voir aussi¶

Références vers d’autres sections :

Références bibliographiques :