DGL lösen: Kalman-Filter im kontinuierlichen Zeitbereich

WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

EliteTUM

Forum-Fortgeschrittener


	Beiträge: 70

	Anmeldedatum: 21.04.11

	Wohnort: München

	Version: ---

Verfasst am: 14.03.2012, 17:44 Titel: DGL lösen: Kalman-Filter im kontinuierlichen Zeitbereich

Hi Leute,

ich sitz gerade vor einer Fragestellung bei der ich mir nicht sicher bin.

Ich habe sowohl im zeitdiskreten als auch zeitkontinuierlichen Bereich mich in den Kalman-Filter (genauer Extendend Kalman-Filter, EKF) eingelesen. Nun ist das Rechnen des EKF für diskrete Modelle ja rel. einfach, da man im Predikt-Schritt einfach die Ergebnisse von vorher nimmt. Die Gleichungen findet man z.B. hier Wiki.

Für den zeitkontinuierlichen Fall, also dt -> 0 ist das ganze schon etwas schwieriger, weil hier ja nicht mehr zwischen einem Predict und einem update-Step unterschieden wird. Die ganzen Gleichungen schrumpfen quasi zusammen und beziehen sich jeweils auf die Werte zum atkuellen Zeitpunkt t. Zu finden unter sind die Gleichungen hier: Wiki.
Laut Särkkä (Quelle, siehe z.B. S. 133) handelt es sich dabei um 3 voneinander abhängige DGLs.

Meine Frage nun, wie löst man diese korrekt?

1. Was ich nach viiiiel viiiel Recherche mitbekommen habe, wäre der wohl korrekte Ansatz das Lösen der allgemeinen Riccati-Matrix-DGL, in Matlab glaube ich mit dem Vefehl care(...) machbar. Korrekt?

2. Auf der Homepage von Prof. Dan Simon (Link) findet man zwei M-Files mit Beispielen für zeitkontinuierliche Kalman-Filter, also Kalman-Bucy-Filter. Die beiden Files wären ExtendedBody.m und Mismodel.m, beide unter dem Link zum Download zu finden. Hier berechnet er das aber ganz anders, viel vereinfachter.

Er geht wie folgt vor:

Code:

% 1. Berechne die Kalman-Verstärkung K mit P vom letzten und F vom aktuellen Zeitschritt
K = P * H' * inv(R);

% 2. Berechne P_Punkt
Pdot = F * P + P * F' - P * H' * inv(R) * H * P;

% 3. Integriere ausgehend vom letzten P und dem gerade berechneten P_Punkt das neue P auf
P = P + Pdot * dt;

% 4. Berechne das aktuelle xHatDot aus den letzten Werten für xHat und der aktuellen Stellgröße u
% siehe File

% 5. Korrigiere xHatdot mit der aktuellen Messung nach der Formel
xhatdot = xhatdot + K * (z - H * xhat);

% 6. Integriere die Schätzung für den Zustand xHat aus der letzten Schätzung xHat und der Änderung xhatdot auf
xhat = xhat + xhatdot * dt;

Funktion ohne Link?

Das ganze entspricht doch Formal nicht wirklich dem Lösen der Riccati-Gleichung, oder? Ich glaube in einer der beiden Dateien kommentiert er es so, als wäre das Berechnen von P_dot wie er es macht das Lösen der Riccati-Gleichung, aber ich hab das so meine Zweifel. Hat wer Rat? Ist das vllt eine mehr oder weniger gültige Näherung für sehr sehr kleine dt?

Abschließende Frage: Wie würdet ihr diese drei gekoppelten DGLs in MatLab für das aktuelle t lösen? Bin nicht sicher, ob/wie ich dafür den Befehl care(...) einsetzen kann.

Bin für jede Hilfe dankbar!

Viele Grüße!
_________________

- EliteTUM
_____________________________________

EliteTUM

Themenstarter

Forum-Fortgeschrittener


	Beiträge: 70

	Anmeldedatum: 21.04.11

	Wohnort: München

	Version: ---

Verfasst am: 15.03.2012, 15:36 Titel:

Hallo zusammen, so ich hab versucht mit folgendem hilfreichen Beitrag von Mathworks (Link) eine Lösung zu finden. Dummerweise hätte ich dafür alle meine nichtlinearen Systeme, dich ich als Funktionen geschrieben habe, umschreiben müssen, damit Sie damit funktionieren.

Ich habe nun den folgenden Ansatz verfolgt:

1. Ich nehme für F und H die Jakobi-Matrizen mit den Werten vom letzten Zeitschritt. Unter Annahme kleiner Änderungen durch sehr kleine Schrittweite eine hoffentlich akzeptable Annäherung.

2. Lösung der Matrix-Riccati-DGL mit ode45.

3. Korrektur von xHatDot mit Innovation, also der Messung.

4. Aufintegration wie beim letzten post gezeigt.

Im Endeffekt ist es wohl fast das gleiche Vorgehen wie von Prof. Dr. Dan Simon, lediglich verwende ich als Integration ode45 (expliziter Runge-Kutta-4,5-Solver) während er Euler-Integration verwendet. Habe beispielhaft mal einen Vergleich für ein Halteglied 1. Ordnung (PT1-Glied), welches verrauscht wurde, gefiltert und den fehler für beide geplottet. EKFDS steht für EKF nach Prof. Dr. Dan Simon, EKF ist für das obige Vorgehen die Abkürzung.

Was ist der Unterschied zwischen meiner Lösung und einem wirklichen Extended-Kalman-Bucy-Filter: So wie ich das verstanden habe müsste man beim Kalman-Bucy-Filter die DGL für d/dt x^(t) gleichzeitig mit der Riccati-Gleichung lösen, weil in ihr die Verstärkung K(t) vorkommt, in K(t) kommt aber P(t) vor, P(t) benötigt wiederum das aktuelle H(t) und F(t) was wiederum zwei Jakobi-Matrizen abhängig von der ursprünglichen DGL d/dt x^(t) ist. Zwar liese sich das ganze machen, aber dafür hätte ich zuviel extra Arbeit einbringen müssen, in Simulink lässt sich das ganze zumindest für den linearen Fall schön zusammenklicken, Yi Cao hat auf der Mathworks-Homepage ein nettes Beispiel online gestellt. Hoffe die hier beschriebene Approximation ist zulässig. Hat da evtl. jemand Ideen dazu die das unterstreichen oder widerlegen?

Was ich interessant finde und nicht ganz beantworten kann: Am Anfang ist der Ansatz mit ode45 langsamer, manchmal sogar sehr deutlich, aber dafür siehts auf Dauer minimal genauer aus, also der Fehler kleiner. Könnte am Euler-Verfahren liegen, das am Anfang bei großen Steigungen schneller (zufällig) den Fehler ausgleicht, aber Gefahr läuft instabil zu werden.

DS_vs_Riccati.png

Beschreibung:

Vergleich: EKF mit Ansatz nach Prof. Dr. Dan Simon und Lösung der Matrix-Riccati-DGL mit ode45

Download

Dateiname:

DS_vs_Riccati.png

Dateigröße:

18.42 KB

Heruntergeladen:

1575 mal

_________________

- EliteTUM
_____________________________________

EliteTUM

Themenstarter

Forum-Fortgeschrittener


	Beiträge: 70

	Anmeldedatum: 21.04.11

	Wohnort: München

	Version: ---

Verfasst am: 16.03.2012, 10:05 Titel:

muss meine Aussage zu dem obigen Bild revidieren: Hatte einen Fehler im Code weshalb die EKF unterschiedlich berechnet wurden. Jetzt sieht die simulation bei beiden nahezu identisch aus, nur minimale Abweichungen, da muss man aber schon extremst ranzoomen. Das dürften dann Fehler aus dem Unterschied Euler <--> ode45 sein.
_________________

- EliteTUM
_____________________________________

der_leo

Forum-Newbie


	Beiträge: 6

	Anmeldedatum: 08.10.10

	Wohnort: Kassel

	Version: ---

Verfasst am: 07.06.2013, 11:12 Titel:

Hi,

ich bin eben auf daselbe Problem gestoßen und wollte nur sagen, dass mir deine ausführliche Beschreibung sehr weitergeholfen hat.

Gruß!!


Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:	Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen

Impressum | Nutzungsbedingungen | Datenschutz | FAQ |

RSS

Hosted by:

Copyright © 2007 - 2024 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.