.. literalinclude:: scripts/simple_KalmanFilter2.py
+.. include:: snippets/Header2Algo10.rst
+
+.. literalinclude:: scripts/simple_KalmanFilter2.res
+
+.. include:: snippets/Header2Algo11.rst
+
+.. _simple_KalmanFilter2_state:
+.. image:: scripts/simple_KalmanFilter2_state.png
+ :align: center
+ :width: 90%
+
+.. _simple_KalmanFilter2_variance:
+.. image:: scripts/simple_KalmanFilter2_variance.png
+ :align: center
+ :width: 90%
+
.. ------------------------------------ ..
.. include:: snippets/Header2Algo06.rst
.. ------------------------------------ ..
.. include:: snippets/Header2Algo07.rst
+- [Welch06]_
- [WikipediaKF]_
],
},
)
+case.setObserver(
+ Info=" Analyzed state at current observation:",
+ Template='ValuePrinter',
+ Variable='Analysis',
+ )
#
XaStep, VaStep = 0., 1.
for i in range(1,len(Yobs)):
print("")
print(" Final a posteriori variance:",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='Noisy measurements')
+plt.plot(Estimates,'r-',label='Estimated state')
+plt.axhline(Xtrue,color='b',label='Truth value')
+plt.legend()
+plt.title('Estimate of the state', fontweight='bold')
+plt.xlabel('Observation step')
+plt.ylabel('Voltage')
+plt.savefig("simple_KalmanFilter2_state.png")
+#
+plt.figure()
+iobs = range(1,len(Observations))
+plt.plot(iobs,Variances[iobs],label='A posteriori error variance')
+plt.title('Estimate of the a posteriori error variance', fontweight='bold')
+plt.xlabel('Observation step')
+plt.ylabel('$(Voltage)^2$')
+plt.setp(plt.gca(),'ylim',[0,.01])
+plt.savefig("simple_KalmanFilter2_variance.png")
--- /dev/null
+Estimation of a constant variable by filtering
+----------------------------------------------
+ Noisy measurements acquired on 50 time steps
+
+ Analyzed state at current observation: [0.]
+ Analyzed state at current observation: [-0.41804504]
+ Analyzed state at current observation: [-0.3114053]
+ Analyzed state at current observation: [-0.31191336]
+ Analyzed state at current observation: [-0.32761493]
+ Analyzed state at current observation: [-0.33597167]
+ Analyzed state at current observation: [-0.35629573]
+ Analyzed state at current observation: [-0.36840289]
+ Analyzed state at current observation: [-0.37392713]
+ Analyzed state at current observation: [-0.36331937]
+ Analyzed state at current observation: [-0.35750362]
+ Analyzed state at current observation: [-0.37963052]
+ Analyzed state at current observation: [-0.37117993]
+ Analyzed state at current observation: [-0.36732985]
+ Analyzed state at current observation: [-0.37148382]
+ Analyzed state at current observation: [-0.36798059]
+ Analyzed state at current observation: [-0.37371077]
+ Analyzed state at current observation: [-0.3661228]
+ Analyzed state at current observation: [-0.36777529]
+ Analyzed state at current observation: [-0.37681677]
+ Analyzed state at current observation: [-0.37007654]
+ Analyzed state at current observation: [-0.37974517]
+ Analyzed state at current observation: [-0.37964703]
+ Analyzed state at current observation: [-0.37514278]
+ Analyzed state at current observation: [-0.38143128]
+ Analyzed state at current observation: [-0.38790654]
+ Analyzed state at current observation: [-0.38880008]
+ Analyzed state at current observation: [-0.38393577]
+ Analyzed state at current observation: [-0.3831028]
+ Analyzed state at current observation: [-0.37680097]
+ Analyzed state at current observation: [-0.37891813]
+ Analyzed state at current observation: [-0.38147782]
+ Analyzed state at current observation: [-0.37981569]
+ Analyzed state at current observation: [-0.38274266]
+ Analyzed state at current observation: [-0.38418507]
+ Analyzed state at current observation: [-0.38923054]
+ Analyzed state at current observation: [-0.38400006]
+ Analyzed state at current observation: [-0.38562502]
+ Analyzed state at current observation: [-0.3840503]
+ Analyzed state at current observation: [-0.38775222]
+ Analyzed state at current observation: [-0.37700787]
+ Analyzed state at current observation: [-0.37328191]
+ Analyzed state at current observation: [-0.38024181]
+ Analyzed state at current observation: [-0.3815806]
+ Analyzed state at current observation: [-0.38392063]
+ Analyzed state at current observation: [-0.38539266]
+ Analyzed state at current observation: [-0.37856929]
+ Analyzed state at current observation: [-0.37744505]
+ Analyzed state at current observation: [-0.37154554]
+ Analyzed state at current observation: [-0.37405773]
+ Analyzed state at current observation: [-0.37725236]
+
+ Final a posteriori variance: [[0.00033921]]
+
iteratively. The keyword "*nextStep*", included in the execution order, allows
to not store the background in duplicate of the previous analysis.
+In a completely similar way to the reanalysis, the estimate is made in
+displaying intermediate results during iterative filtering. Thanks to these
+intermediate information, one can also obtain the graphs illustrating the
+estimation of the state and the associated *a posteriori* error covariance.
.. literalinclude:: scripts/simple_KalmanFilter2.py
+.. include:: snippets/Header2Algo10.rst
+
+.. literalinclude:: scripts/simple_KalmanFilter2.res
+
+.. include:: snippets/Header2Algo11.rst
+
+.. _simple_KalmanFilter2_state:
+.. image:: scripts/simple_KalmanFilter2_state.png
+ :align: center
+ :width: 90%
+
+.. _simple_KalmanFilter2_variance:
+.. image:: scripts/simple_KalmanFilter2_variance.png
+ :align: center
+ :width: 90%
+
.. ------------------------------------ ..
.. include:: snippets/Header2Algo06.rst
.. ------------------------------------ ..
.. include:: snippets/Header2Algo07.rst
+- [Welch06]_
- [WikipediaKF]_
],
},
)
+case.setObserver(
+ Info=" État analysé à l'observation courante :",
+ Template='ValuePrinter',
+ Variable='Analysis',
+ )
#
XaStep, VaStep = 0., 1.
for i in range(1,len(Yobs)):
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")
--- /dev/null
+Estimation par filtrage d'une variable constante
+------------------------------------------------
+ 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]]
+
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, 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.
if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon !
_Innovation = _Innovation - Cm * Un
#
- _A = R + numpy.dot(Ht, Pn_predicted * Ha)
- _u = numpy.linalg.solve( _A , _Innovation )
- Xn = Xn_predicted + Pn_predicted * Ha * _u
- Kn = Pn_predicted * Ha * (R + numpy.dot(Ht, Pn_predicted * Ha)).I
+ Kn = Pn_predicted * Ha * numpy.linalg.inv(R + numpy.dot(Ht, Pn_predicted * Ha))
+ Xn = Xn_predicted + Kn * _Innovation
Pn = Pn_predicted - Kn * Ht * Pn_predicted
- Xa, _HXa = Xn, _HX # Pointeurs
+ Xa = Xn # Pointeurs
#
self.StoredVariables["CurrentIterationNumber"].store( len(self.StoredVariables["Analysis"]) )
# ---> avec analysis
self.StoredVariables["Analysis"].store( Xa )
if self._toStore("SimulatedObservationAtCurrentAnalysis"):
- self.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( _HXa )
+ self.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( Hm((Xn, None)) )
if self._toStore("InnovationAtCurrentAnalysis"):
self.StoredVariables["InnovationAtCurrentAnalysis"].store( _Innovation )
# ---> avec current state
if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon !
_Innovation = _Innovation - Cm * Un
#
- _A = R + numpy.dot(Ht, Pn_predicted * Ha)
- _u = numpy.linalg.solve( _A , _Innovation )
- Xn = Xn_predicted + Pn_predicted * Ha * _u
- Kn = Pn_predicted * Ha * (R + numpy.dot(Ht, Pn_predicted * Ha)).I
+ Kn = Pn_predicted * Ha * numpy.linalg.inv(R + numpy.dot(Ht, Pn_predicted * Ha))
+ Xn = Xn_predicted + Kn * _Innovation
Pn = Pn_predicted - Kn * Ht * Pn_predicted
- Xa, _HXa = Xn, _HX # Pointeurs
+ Xa = Xn # Pointeurs
#
self.StoredVariables["CurrentIterationNumber"].store( len(self.StoredVariables["Analysis"]) )
# ---> avec analysis
self.StoredVariables["Analysis"].store( Xa )
if self._toStore("SimulatedObservationAtCurrentAnalysis"):
- self.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( _HXa )
+ self.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( Ht * Xn )
if self._toStore("InnovationAtCurrentAnalysis"):
self.StoredVariables["InnovationAtCurrentAnalysis"].store( _Innovation )
# ---> avec current state
