\chapter{Robuste stochastische Zeitreihenmodelle in der End-of-Line-Prüfung} \label{ch:timeseries} In diesem Kapitel beheben wesentlicher Schwächen der statischen Klassifikation für die \ac{EOL}-Prüfung. Erstens systematische Analyse von Langzeit-Veränderungen der Istwerte von Prüfmerkmalen (Drift). Zweitens stochastischer Modellierungsansatz, um zeitabhängigen Unsicherheiten gezielt in die Prüfgrenzen einfließen zu lassen und die Prüfschärfe frei einstellen zu können. Dabei Erhalt der Robustheit gegen falsche Labels in den historischen Messdaten und Berücksichtigung relevanter statistischer Eigenschaften der Prüfmerkmale. Die Inhalte basieren auf der Veröffentlichung \textfixme{INES2018?} \cite{Leitner2018b}. Langsame zeitliche Veränderungen im Verhalten industrieller Herstellungsprozesse sind in der Literatur bekannt und können wie in \chref{ch:static} beobachtet die Spezifität einer automatisierten Fehlererkennung verringern \cite{Qin2012}. Eine erste mögliche Maßnahme zur Erhöhung der Spezifität ist, die gesetzten Grenzwerte (hier: Prüfgrenzen) von Anfang an so weit zu wählen, dass trotz Veränderungen in den Istwerten zumindest mittelfristig nur wenige Fehlalarme auftreten. Das Aufweiten der Prüfgrenzen hat sich bei der in \chref{ch:static} untersuchten \ac{SVDD} allerdings als problematisch erwiesen, da außerhalb der \ac{iO}-Trainingsdaten keine Messwerte zur Festlegung des neuen Schwellwertes vorhanden sind. Um die automatisch erzeugten Prüfgrenzen im Sinne einer Extrapolation über den bekannten \ac{iO}-Bereich hinaus auszudehnen bietet sich die Nutzung stochastischer Verfahren an. Diese quantifizieren die Unsicherheiten in den Istwerten durch die genutzten Wahrscheinlichkeitsverteilungen, sodass eine beliebige, definierte Einstellung der Prüfschärfe möglich ist. Der wesentliche Nachteil beim einfachen Aufweiten der Prüfgrenzen zeigt sich in Form einer zwangsweise verringerten Sensitivität der Prüfung. Ein zweiter möglicher Ansatz zum Umgang mit zeitlichen Veränderungen in den Istwerten ist die regelmäßige Wiederholung des Trainingsvorgangs, um dem aktuellen Stand der Istwerte möglichst schnell zu folgen. Wie in \secref{subsec:static_results} diskutiert besteht diese Möglichkeit aber nur, wenn regelmäßig ohne längere Unterbrechungen Messdaten aus dem Produktionsprozess generiert werden. Diese Voraussetzung ist in der \ac{EOL}-Prüfung regelmäßig nicht erfüllt. Stochastische Zeitreihenmodelle bieten das Potential, die Vorteile der beiden zuvor genannten Ansätze miteinander zu verbinden. Die Streuung der Istwerte quantifizieren Wahrscheinlichkeitsverteilungen, wodurch die Prüfschärfe einstellbar wird. Gleichzeitig berücksichtigen die Modelle explizit die zeitliche Entwicklung der betrachteten Größen und die damit verbundenen Unsicherheiten. Das Ziel sind dynamisch angepasste Prüfgrenzen, die zu jedem Zeitpunkt optimal an die aktuelle Ist-Situation angepasst sind. Dieses Kapitel beginnt mit einer einführenden Betrachtung des zeitlichen Driftverhaltens anhand von beispielhafen Prüfmerkmalen. Danach wird zunächst das gängige linear-gaußsche Zustandsraummodell aufgegriffen und im Rahmen der \ac{EOL}-Prüfung interpretiert. Darauf aufbauend erfolgt die Untersuchung robuster Zustandsschätzer, die mit nicht-normalverteilten Rauschannahmen arbeiten. Sie ermöglichen eine Zustandsschätzung und -prädiktion, die robust gegen Ausreißer ist. Eine Erweiterung auf schiefe Merkmalsverteilungen ist ebenfalls Inhalt dieses Kapitels. Aus den so gefundenen Schätz- und Prädiktionsgrößen werden Prüfgrenzen abgeleitet, die im Prüfbetrieb zum Einsatz kommen. Die Basis dafür bildet die robuste numerische Optimierung der Modellparameter, welche auf den robusten Zustandsschätzern aufbaut und ebenfalls vorgestellt wird. Das Kapitel schließt mit einer Auswertung der Anwendungsergebnisse. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Driftverhalten der Istwerte von Prüfmerkmalen} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Das Langzeitverhalten der Istwerte von zwei beispielhaften Prüfmerkmalen aus dem Anwendungsbeispiel der \ac{EOL}-Prüfung von Verbrennungsmotoren zeigt \figref{fig:tsExampleDrift}. Beide Prüfmerkmale stammen aus der gleichen Teilprüfung (Aktivierung der Nockenwellenverstellung, vgl. \figref{fig:drehzahlverlauf_pruefung}). Obwohl es sich bei beiden Merkmalen um Luftdruckmessungen am Prüfling handelt, ist das individuelle Driftverhalten deutlich unterschiedlich. Merkmal \num{213} weist außer der Drifteigenschaft zusätzlich eine starke Schiefe der Verteilung der Merkmalswerte auf: die Istwerte streuen um einen gedachten gleitenden Mittelwert deutlich stärker hin zu größeren Merkmalswerten als zu kleineren. Bei beiden Prüfmerkmalen liegt die Schwankung des Mittelwertes in der Größenordnung der Streuung der Einzelprüflinge und ist damit nicht mehr vernachlässigbar. Die potentiellen Probleme statischer Prüfgrenzen werden in \figref{fig:tsExampleDrift} nochmals deutlich: tritt eine Produktionspause wie hier zwischen Tag \num{81} und Tag \num{100} auf, sind zuvor erzeugte statische Prüfgrenzen nach der Pause unter Umständen nicht mehr repräsentativ für die neue Lage der Istwerte. Da die Prüfung nach der Unterbrechung sofort wieder starten soll, ist eine Berücksichtigung der zeitlichen Veränderungen und insbesondere von deren Unsicherheiten in den Prüfgrenzen erforderlich. %Auswahl Merkmale %Klassiker: idWindow 52 $\rightarrow$ Ansaugdruck dynamisch, Stufe NWVerstellungEin %ICIT idWindow 19 %INES idWindow 233 $\rightarrow$ Unterdruck Vak, Stufe NWVerstellungEin \begin{figure} \centering \footnotesize \setlength\figurewidth{.9\columnwidth} \setlength\figureheight{.3\columnwidth} \renewcommand\figureXLabel{\mathSymbolWithUnit{t}{d}} \renewcommand\figureYLabel{$\meas_\indexS$} \subcaptionbox{Prüfmerkmal \num{47}: rel. Ansaugdruck am Ansaugadapter\label{fig:tsExampleDriftOne}}{% Merkmal 47 = idWindow 52 \centering \footnotesize \mytikzexton \input{"graphics/P1015 02(all) 11a/001 windowGroup=001, time series, time-based.tikz"} \mytikzextoff }% % --------------- comment magic (tm) \\% \vspace{4mm}% \hspace{5mm}%special hspace for horizontal alignment of the second plot % --------------- comment magic (tm) \subcaptionbox{Prüfmerkmal \num{213}: Unterdruck der Vakuumpumpe\label{fig:tsExampleDriftTwo}}{% Merkmal 213 = idWindow 233 \centering \footnotesize \mytikzexton \input{"graphics/P1015 02(all) 11b/001 windowGroup=001, time series, time-based.tikz"} \mytikzextoff }% \caption[FIXME]{Zwei Beispiele für stark driftende Prüfmerkmale. Dargestellt sind standardisierte Istwerte \circleWithParens{myLineOne} des Motortyps M1 aus \num{2348} Prüfläufen auf Prüfstand P2 jeweils aus dem selben Zeitraum von \num{138} Tagen sowie die manuellen Prüfgrenzen \dashedLineWithParens{myManualLimits}. Verwendet wurden nur Daten von Prüflingen, die ein \ac{iO}-Prüfergebnis gemäß der manuellen Prüfgrenzen erhalten haben.} \label{fig:tsExampleDrift} \end{figure} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Schema einer \ac{EOL}-Prüfung mit Zeitreihenmodellen} \label{sec:tsSystemScheme} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Hinweis Vorverarbeitung Untergruppen wie in chapter static \textfixme{Bild für Schema erstellen und beschreiben, Fokus Aufbau Gesamtsystem aus Teilblöcken} \textfixme{Bild für Ablauf Filter, Berechnung Prüfgrenzen, Berechnung Likelihood, Optimierung erstellen und beschreiben, Fokus Berechnungen qualitativ} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Linear-gaußsches Zustandsraummodell für die \ac{EOL}-Prüfung} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Das lineare-gaußsche stochastische Zustandsraummodell aus \secref{sec:basics_recursive} dient im Folgenden als Basis für mehrere robuste Erweiterungen. Ausgangspunkt ist das Modell \begin{subequations} \label{eq:tsLGSS} \begin{IEEEeqnarray}{rCl} \stateV_{\indexS+1} &=& \stateM \stateV_{\indexS} + \inputM \inputV_{\indexS} + \procNV_{\indexS}, \IEEEyessubnumber* \label{eq:tsLGSSStateTransition}\\ \measV_{\indexS} &=& \outputM \stateV_{\indexS} + \feedM \inputV_{\indexS} + \measNV_{\indexS}, \label{eq:tsLGSSMeasurement}\\ \procNV_\indexS &\sim& \pdfNBr{\nullV, \procNCM_\indexS}, \label{eq:tsLGSSDistProcN} \\ \measNV_\indexS &\sim& \pdfNBr{\nullV, \measNCM}. \label{eq:tsLGSSDistMeasN} \end{IEEEeqnarray} \end{subequations} mit unabhängigen Zufallsvariablen $\procNV_\indexS$ und $\measNV_\indexS$. Die Parameter $\stateM \in \dimRTwo{\nStateV}{\nStateV}$, $\inputM \in \dimRTwo{\nStateV}{\nInputV}$, $\outputM \in \dimRTwo{\nMeasV}{\nStateV}$, $\feedM \in \dimRTwo{\nMeasV}{\nInputV}$, $\procNCM_\indexS \in \dimRTwo{\nStateV}{\nStateV}$ und $\measNCM \in \dimRTwo{\nMeasV}{\nMeasV}$ sind fest vorgegeben, außerdem sind $\procNCM_\indexS$ und $\measNCM$ positiv semidefinit. Die Istwerte der Prüfmerkmale zum Zeitpunkt $\indexS \in \N_+$ bilden wie in den vorangegangenen Kapiteln den Vektor $\measV_\indexS$, d.h. jedem Prüflauf ist ein $\measV_\indexS$ zugeordnet. Die Modellgleichungen \eqref{eq:tsLGSS} erlauben eine anschauliche Interpretation im Zusammenhang mit der \ac{EOL}-Prüfung. Der Zustandsvektor $\stateV_\indexS$ entspricht dem aktuellen, nicht messbaren Mittel der produzierten Prüflinge. Der Zustandsvektor (bzw. genauer: seine lineare Abbildung $\outputM \stateV_\indexS + \feedM \inputV_\indexS$) kann als Vektor der Istwerte eines hypothetischen Prüflings betrachtet werden, der zum Zeitpunkt $\indexS$ produziert worden wäre wenn keine unbekannten äußeren Einflüsse auf Produktion und Prüfung eingewirkt hätten. Zufällige Schwankungen um den Zustand welche nur den aktuellen Prüfling betreffen sind im Messrauschen $\measNV_\indexS$ zusammengefasst. Das Messrauschen umfasst dabei nicht nur Rauscheffekte im Prüfvorgang selbst (z.B. Sensorrauschen), sondern auch zufällige Abweichungen in der Fertigung (Abweichungen von Sollmaßen der Bauteile) und in der Montage (Montagetoleranzen). Im Gegensatz dazu beschreibt das Prozessrauschen $\procNV_\indexS$ längerfristige Abweichungen in der Entwicklung des Zustandsvektors, die somit einen anhaltenden Effekt auf die Messwerte mehrerer aufeinanderfolgender Prüflinge haben. Das Prozessrauschen ist daher geeignet um Drifteffekte der Produktion (z.B. Werkzeugverschleiß) und langsam veränderliche Umgebungseinflüsse im Modell aufzunehmen, wenn diese nicht als messbare Größen verfügbar sind. Zusammenfassend beschreibt das Prozessrauschen also die Unsicherheit durch langsam veränderliche, persistierende Prozesse, während das Messrauschen die momentanen Schwankungen der Einzelmessungen darstellt. Diese Interpretation entspricht dem in der Prozessüberwachung üblichen Vorgehen \cite[][295]{Venkatasubramanian2003}. Deterministische Einflüsse können soweit sie bekannt sind durch die Eingänge $\inputV_\indexS$ in das Modell einfließen. Die Normalverteilungsannahmen für $\procNV_\indexS$ \eqref{eq:tsLGSSDistProcN} und $\measNV_\indexS$ \eqref{eq:tsLGSSDistMeasN} sind Teil des stochastischen Zustandsraummodells \cite[][631\psq]{Murphy2013}. Die Modellgrößen $\stateV_\indexS$ und $\measV_\indexS$ geben dabei das Normalverhalten wieder \cite[][306]{Venkatasubramanian2003}, sie beschreiben im Rahmen der \ac{EOL}-Prüfung also die \ac{iO}-Prüfläufe. Das rekursive Filterproblem zur Schätzung und Prädiktion des Zustandsvektors $\stateV_\indexS$ für \eqref{eq:tsLGSS} löst das Kalman-Filter exakt (vgl. \secref{sec:basics_recursive}), wenn zusätzlich zu \eqref{eq:tsLGSS} die A-priori-Verteilung des geschätzten Zustandsvektors als $\pdfBr{\stateV_1} = \pdfNBr{\est{\stateV}_{\given{ 1 }{ 0 }}, \stateCM_{\given{ 1 }{ 0 }}}$ angenommen wird. %============================================================================== \subsection{Random-Walk-Modell als Spezialfall} %============================================================================== Da die Prozesse im Gesamtsystem Fertigung-Montage-Prüfung die zur letztendlichen Ausprägung der Istwerte der Prüfmerkmale führen äußerst komplex sind, ist eine Modellierung auf Basis physikalischer Gleichungen erstens aufwendig und zweitens auch nur eingeschränkt möglich. Für viele Einflussgrößen z.B. aus der Fertigung sind im Rahmen der \ac{EOL}-Prüfung die nötigen Messwerte gar nicht verfügbar, so dass schon der Systemeingang $\inputV_\indexS$ nicht vorliegt. Entsprechend steigt die Bedeutung der stochastischen Modellkomponenten $\procNV_\indexS$ und $\measNV_\indexS$ und ihrer Verteilungsannahmen gegenüber den deterministischen Parametern $\stateM$, $\inputM$, $\outputM$ und $\feedM$. Für die \ac{EOL}-Prüfung bietet sich die Nutzung eines Random-Walk-Modells \cite[][65]{BarShalom2001} für den Zustandsvektor an. Dabei weist $\stateV_\indexS$ keine eigene Dynamik mehr auf: \begin{subequations} \label{eq:tsLGSSRandomWalk} \begin{IEEEeqnarray}{rCl} \stateV_{\indexS+1} &=& \stateV_{\indexS} + \inputM \inputV_\indexS + \procNV_{\indexS}, \IEEEyessubnumber* \label{eq:tsLGSSRandomWalkState}\\ \measV_{\indexS} &=& \outputM \stateV_{\indexS} + \feedM \inputV_\indexS + \measNV_{\indexS}, \\ \procNV_\indexS &\sim& \pdfNBr{\nullV, \procNCM_\indexS}, \label{eq:tsLGSSRandomWalkProcN} \\ \measNV_\indexS &\sim& \pdfNBr{\nullV, \measNCM}. \end{IEEEeqnarray} \end{subequations} Es handelt sich also um einen Spezialfall von \eqref{eq:tsLGSS} mit $\stateM = \eye$. Die Rauschterme sind nun wesentlich für die Modellierung der Veränderungen der Merkmalswerte verantwortlich. Der Eingang $\inputV_\indexS$ kann auch ganz entfallen, wenn die dadurch entstehenden größeren Unsicherheiten durch eine größere Varianz der Rauschterme berücksichtigt werden. %============================================================================== \subsection{Unregelmäßige Zeitabstände zwischen Prüfläufen} \label{sec:tsTimeDistUnequal} %============================================================================== Wie in \figref{fig:tsExampleDrift} bereits erkennbar ist, schwanken die Zeitabstände zwischen den einzelnen Prüfvorgängen teils erheblich. Dieser Umstand ist in der Modellierung mit Zeitreihenmodellen zu berücksichtigen. % ------------------------------------- \subsubsection{Berücksichtigung im Zustandsraummodell} % ------------------------------------- Klassische regelungstechnische Anwendungen von Zustandsmodellen betrachten die Zeitabstände zwischen den einzelnen Abtastschritten als konstant \cite[][171\psq, 200\psq]{Schroder2009} \cite[][477\psq]{Lutz2010}. Für fest getaktete Algorithmen auf Echtzeitsystemen ist diese Annahme zulässig und sinnvoll. Die Parameter des zeitdiskreten Modells können damit z.B. aus einer zeitkontinuierlichen Zustandsraum-Systembeschreibung hergeleitet und fest im Regelalgorithmus parametriert werden. Sind die Zeitschritte hingegen nicht konstant, dann ist eine zeitdiskrete Beschreibung mit konstanter Systemmatrix $\stateM$ wie in \eqref{eq:tsLGSSStateTransition} im Allgemeinen nicht mehr sinnvoll, da sich das Übergangsverhalten des Zustandsvektors von einem Zeitschritt zum nächsten durch die unterschiedlich langen Zeitschritte ändert. Die Berechnung des Zustandsübergangs erfordert in diesem Fall die Integration des zeitkontinuierlichen Systemverhaltens, was auf die aufwendige Berechnung eines Matrix-Exponentials und dessen Integral in jedem Zeitschritt führt \cite[][32-35]{Astrom1997}. %\cite{vanLoan1978} \cite[][635-640]{Lutz2010} Außerdem ist neben dem deterministischen Systemverhalten auch der Einfluss des Prozessrauschens über die Länge des jeweiligen Zeitschrittes zu integrieren, was letztendlich dem Lösen einer linearen stochastischen Differentialgleichung entspricht. %Diskretisierung LTI-SDE: %Eine lineare zeitinvariante stochastische Differentialgleichung hat die Form \cite[][43-45]{Sarkka2006} %\begin{IEEEeqnarray*}{c} %\diff \stateV = \stateM_\cont \stateV(t) \diff t + \inputV(t) + \dispM_\cont \diff \wienV(t) %\end{IEEEeqnarray*} Ein einfaches Beispiel für ein zeitkontinuierliches Zustandsraummodell in Form einer solchen stochastischen Differentialgleichung ist \cite[][398]{Astrom1997} \begin{IEEEeqnarray*}{c} \frac{\diff \stateV(t)}{\diff t} = \stateM_\cont \stateV(t) + \frac{\diff \wienV_\cont(t)}{\diff t} \end{IEEEeqnarray*} oder in der dafür üblicheren Schreibweise \begin{IEEEeqnarray*}{c} \diff \stateV(t) = \stateM_\cont \stateV(t) \diff t + \diff \wienV_\cont(t). \end{IEEEeqnarray*} Die Elemente des Vektors $\diff \wienV_\cont / \diff t$ enthalten weißes Rauschen. Im Folgenden sei zur Veranschaulichung $\wienV_\cont$ eine Brownsche Bewegung (Wiener-Prozess), d.h. die Inkremente zwischen Zeitpunkt $t$ und $t^\prime$, $t < t^\prime$ sind normalverteilt mit \cite[][397\psq]{Astrom1997} \begin{IEEEeqnarray*}{c} \wienV_\cont(t) - \wienV_\cont(t^\prime) \sim \pdfNBr{\nullV, \procNCM \left(t - t^\prime \right)}. \end{IEEEeqnarray*} Die exakte Lösung der stochastischen Differentialgleichung zu diskreten Zeitpunkten $t_\indexS$ ist \cite[][402\psq]{Astrom1997} \begin{subequations} \label{eq:tsSDEDiscrete} \begin{IEEEeqnarray}{rCl} \stateV(t_{\indexS+1}) &=& \expBr{\stateM_\cont \left(t_{\indexS+1} - t_{\indexS} \right)} \stateV(t_\indexS) + \wienV(t_\indexS), \IEEEyessubnumber* \\ \wienV(t_\indexS) &=& \int_{t_\indexS}^{t_{\indexS+1}}{\expBr{\stateM_\cont \left( t_{\indexS+1} - \tau \right)} \diff \wienV_\cont(\tau)}. \end{IEEEeqnarray} \end{subequations} % Dafür können näherungsweise zeitdiskrete Lösungen mit Approximationsverfahren wie der Euler-Maruyama-Approximation als stochastischer Variante der Euler-Approximation gefunden werden \cite[][305]{Kloeden1992}. Für das zeitkontinuierliche Random-Walk-Modell gilt aufgrund der fehlenden Zustands-Eigendynamik $\stateM_\cont = \nullM$, und die Lösung \eqref{eq:tsSDEDiscrete} vereinfacht sich sich zu \begin{subequations} \label{eq:tsSDEWalkSolution} \begin{IEEEeqnarray}{rCl} \stateV(t_{\indexS+1}) &=& \stateV(t_\indexS) + \wienV(t_\indexS), \IEEEyessubnumber* \label{eq:tsSDEWalkSolutionX}\\ \wienV(t_\indexS) &=& \int_{t_\indexS}^{t_{\indexS+1}}{ \diff \wienV_\cont(\tau)} = \wienV_\cont\left(t_{\indexS+1}\right) - \wienV_\cont\left(t_\indexS\right) \\ \Rightarrow \wienV(t_\indexS) &\sim& \pdfNBr{\nullV, \procNCM \left(t_{\indexS+1} - t_\indexS \right)}. \label{eq:tsSDEWalkSolutionV} \end{IEEEeqnarray} \end{subequations} % Aus \eqref{eq:tsSDEWalkSolutionX} und \eqref{eq:tsSDEWalkSolutionV} lässt sich nun bereits das zeitdiskrete Random-Walk-Modell des Zustandsvektors in \eqref{eq:tsLGSSRandomWalkState} und \eqref{eq:tsLGSSRandomWalkProcN} für den Fall $\inputV_\indexS = \nullV$ ablesen. Es gilt daher \begin{IEEEeqnarray*}{rCl} \procNCM_\indexS &=& \procNCM \cdot \Delta t_\indexS, \\ \Delta t_\indexS &=& t_{\indexS+1} - t_\indexS. \end{IEEEeqnarray*} Zusammengefasst ermöglicht also das Random-Walk-Modell eine Berücksichtigung variabler Zeitabstände in den einzelnen Zeitschritten mit geringem Rechenaufwand. Dazu muss die auf eine Zeiteinheit bezogene Kovarianzmatrix $\procNCM$ des Prozessrauschens mit der Länge des jeweiligen Zeitschrittes skaliert werden. % ------------------------------------- \subsubsection{Verteilung der Zeitabstände im Anwendungsbeispiel} % ------------------------------------- Die Zeitstempel der vorhandenen Messdaten aus dem Anwendungsbeispiel ermöglichen eine Analyse der Verteilung der Zeitabstände $\Delta t_\indexS$ zwischen dem Ende von je zwei aufeinanderfolgenden Prüfläufen. \Figref{fig:tsTimeSteps} zeigt das Histogramm einer Langzeitauswertung für einen \ac{EOL}-Prüfstand und einen Motortyp. Das Minimum der Zeitabstände beträgt $\Delta t_{\min} = \SI{97}{s}$ und entspricht zwei ohne Pause direkt aufeinanderfolgenden Prüfläufen. Längere Abstände entstehen durch unter anderem durch nicht voll ausgelastete Produktionskapazitäten, aufeinanderfolgende Produktionslose verschiedener Motortypen, Wartungszeitfenster und generelle Produktionspausen. Der längste Abstand im Datensatz beträgt \num{19.17} Tage. \begin{figure} \centering \footnotesize \setlength\figurewidth{.9\columnwidth} \setlength\figureheight{.3\columnwidth} \renewcommand\figureXLabel{\mathSymbolWithUnit{\Delta t_\indexS}{s}} \renewcommand\figureYLabel{$\pdfBr{\Delta t_\indexS}$} \mytikzexton \input{"graphics/P1015 02(all) 12a/001 windowGroup=001, time steps.tikz"} \mytikzextoff \caption[FIXME]{Histogramm der Zeitabstände $\Delta t_\indexS$ zwischen je zwei Prüfläufen~\squareWithParens{myLineOne} (normiert als Schätzung der Wahrscheinlichkeitsdichte) und die Maximum-Likelihood-Schätzung der logarithmischen Normalverteilung~\lineWithParens{myLineTwo}. Die Daten umfassen alle \num{2548} Prüfläufe des Motortyps M1 auf Prüfstand P2 in einem Zeitraum von \num{138} Tagen.} \label{fig:tsTimeSteps} \end{figure} Um nachfolgend realistische simulierte Datensätze für Analysen an Zeitreihenmodellen und Algorithmen erzeugen zu können bietet es sich an, die Zeitabstände mit einer Wahrscheinlichkeitsverteilung nachzubilden. Für die Modellierung von Wartezeiten und Zeitabständen kommen in der Literatur unter anderem die Exponentialverteilung \cite[][100-102]{Bamberg2012} und die logarithmische Normalverteilung \cite[][2178\psq]{Nist2012} zum Einsatz. Um auch das Auftreten längerer Zeitabstände wiederzugeben wird nachfolgend die logarithmische Normalverteilung verwendet, da diese lange Ausläufer (heavy tails) hin zu langen Wartezeiten aufweist \cite[][1006\psq]{Bronstein2013}. Eine Zufallsvariable $\rva \in \R_{>0}$ mit logarithmischer Normalverteilung hat die Wahrscheinlichkeitsdichtefunktion \cite[][576\psq]{Gelman2014} \begin{IEEEeqnarray*}{c} \pdfBr{\given{\rva}{\meanSym, \stdSym^2}} = \frac{1}{\sqrt{2\pi} \stdSym \rva} \expBr{ - \frac{1}{2 \stdSym^2} \left( \ln \rva - \meanSym \right)^2 }. \end{IEEEeqnarray*} Sie entspricht der Wahrscheinlichkeitsdichte einer Zufallsvariablen $\rva \in \R_{>0}$, deren Logarithmus mit $\ln \rva \sim \pdfNBr{\meanSym, \stdSym^2}$ normalverteilt ist. Entsprechend einem Ansatz von Spinoso et al. \cite{Spinoso2014} erfolgt die Modellierung von $\Delta t_\indexS$ hier allerdings nicht direkt als logarithmische Normalverteilung. Anstatt dessen wird der $\Delta t_{\min}$ überschreitende Teil $\Delta t_\indexS - \Delta t_{\min}$ der Zeitabstände als logarithmisch normalverteilt modelliert, um die Mindestdauer der Prüfvorgänge einfließen zu lassen. Eine Maximum-Likelihood-Schätzung der Parameter ergibt \begin{IEEEeqnarray*}{c} \meanSym = \num{4.31}, \quad \stdSym = \num{2.80}. \end{IEEEeqnarray*} Die zugehörige Wahrscheinlichkeitsdichtefunktion ist in \figref{fig:tsTimeSteps} dargestellt. %Beispiel für variable Zeitabstände tStepSamples mit Gamma-Verteilung (alt) bzw. log-Normalverteilung (neu): 05 Drift / AdaptiveTracking / EstimateKalmanWindowGroupBatch Test.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Robuste Filter für stochastische Zustandsraummodelle} \label{sec:tsRobustFilter} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Aufbauend auf dem linear-gaußschen Zustandsraummodell \eqref{eq:tsLGSS} werden nun Modellvarianten und zugehörige Filteralgorithmen vorgestellt, die robust gegen Ausreißer sind. Die primäre Quelle für Ausreißer in der \ac{EOL}-Prüfung stellen dabei die \ac{niO}-Prüfläufe dar. Diese führen durch ihre extremen Merkmalswerte bei nicht robusten Verfahren wie z.B. dem Kalman-Filter zu einer starken Verfälschung der Zustandsschätzung. Die Auswirkungen betreffen über die Prädiktion der nachfolgenden Merkmalswerte zunächst die Prüfgrenzen und damit das Prüfergebnis nachfolgender Prüfläufe. Genauso beeinflusst wird aber auch die Schätzung der Modellparameter im Lernvorgang. In \secref{sec:basics_robust} des Grundlagenkapitels wurden bereits zwei Möglichkeiten zum Umgang mit Ausreißern in Schätzverfahren vorgestellt: zum einen der Ausschluss aus dem Schätzvorgang z.B. durch einen Schwellwert, und zum anderen die Einführung von datenabhängigen Gewichtungsfaktoren. Beide Verfahren werden nachfolgend verwendet, wobei der Fokus auf verschiedenen Ansätzen zur geringeren Gewichtung von Ausreißern liegt. %Die Gewichtungsfaktoren werden dabei anhand von Verteilungsannahmen hergeleitet, um weiterhin ein stochastisches Modell zur Prädiktion der Merkmalswerte und Ableitung der Prüfgrenzen zur Verfügung zu haben. Dabei kommen auch die in \secref{sec:basics_robust} vorgestellten Heavy-Tail-Verteilungen - insbesondere die Student-t-Verteilung - zum Einsatz \cite{Meinhold1989}. Die folgenden Abschnitte beschreiben robuste Verfahren zur Zustandsschätzung, wobei der Fokus zunächst auf der Robustheit gegen Ausreißer liegt und anschließend um den Aspekt der Schiefe erweitert wird. Der Vergleich der Verfahren erfolgt bereits konkret im Hinblick auf die \ac{EOL}-Prüfung, da die Eignung der robuster Verfahren von der konkreten Anwendung abhängt \cite[][62]{Pearson2002}. Als Ausgangspunkt dient das Kalman-Filter aus \secref{sec:basics_recursive}. \Algoref{alg:tsKalmanPrediction} und \algoref{alg:tsKalmanUpdate} zeigen nochmals explizit die Filtergleichungen dafür. \begin{algorithm}[htb] \caption{Prädiktionsschritt des \acl{KF}s (\acs{KF})} \label{alg:tsKalmanPrediction} \begin{algorithmic}[1] \setstretch{1.2} \State $ \est{\stateV}_{\given{\indexS}{\indexS-1}} \gets \stateM \est{\stateV}_{\given{\indexS-1}{\indexS-1}} + \inputM \inputV_{\indexS-1} $ \State $ \stateCM_{\given{\indexS}{\indexS-1}} \gets \stateM \stateCM_{\given{\indexS-1}{\indexS-1}} \trans{\stateM} + \procNCM_{\indexS-1} $ \State $ \est{\measV}_{\given{\indexS} {\indexS-1} } \gets \outputM \est{\stateV}_{\given{\indexS}{\indexS-1}} + \feedM \inputV_{\indexS} $ \State $ \measCM_{\indexS} \gets \outputM \stateCM_{\given{\indexS}{\indexS-1}} \trans{\outputM} + \measNCM $ \end{algorithmic} \end{algorithm} \begin{algorithm}[htb] \caption{Update-Schritte des \acl{KF}s (\acs{KF})} \label{alg:tsKalmanUpdate} \begin{algorithmic}[1] \setstretch{1.2} \State $ \resV_{\indexS} \gets \measV_\indexS - \est{\measV}_{\given{\indexS} {\indexS-1} } $ \State $ \kalmM_{\indexS} \gets \stateCM_{\given{\indexS}{\indexS-1}} \trans{\outputM} \inv{\measCM}_{\indexS} $ \State $ \est{\stateV}_{\given{\indexS}{\indexS}} \gets \est{\stateV}_{\given{\indexS}{\indexS-1}} + \kalmM_{\indexS} \resV_{\indexS} $ %\stateCM_{\given{\indexS}{\indexS}} &=& \left( \eye - \kalmM_{\indexS} \outputM \right) \stateCM_{\given{\indexS}{\indexS-1}} \left( \eye - \kalmM_{\indexS} \outputM \right)^T \IEEEnonumber \IEEEnosubnumber \\ % &&{}+{} \kalmM_{\indexS} \measNCM \kalmM_{\indexS}^T \label{eq:KalmanJosephFormUpdate} \State $ \stateCM_{\given{\indexS}{\indexS}} \gets \stateCM_{\given{\indexS}{\indexS-1}} - \kalmM_{\indexS} \measCM_{\indexS} \trans{\kalmM}_{\indexS} $ \end{algorithmic} \end{algorithm} %============================================================================== \subsection{\acf{GKF}} %============================================================================== Bereits für die historisch ersten Anwendungen des Kalman-Filters in der Raumfahrt kamen Techniken zum Einsatz, um mit Hilfe der vom Filter ohnehin gelieferten Schätzvarianz Ausreißer-Messungen vom Update-Schritt auszuschließen \cite[][12]{McGee1985}. Liegt dabei die Messung $\measV_\indexS$ außerhalb des auf Basis der vorangegangenen Messungen $\measV_{1:\indexS-1}$ zu erwartenden Wertebereichs, dann wird $\measV_\indexS$ vollständig ignoriert. Diese Vorgehensweise ist heute als (Validation-) Gating in der Literatur bekannt \cite[][153\psq]{BarShalom2001}. % sec. 3.5.3 Unter Annahme des linear-gaußschen Zustandsraummodells \eqref{eq:tsLGSS} gilt für Prädiktion des Zustandsvektors $\stateV_\indexS$ \begin{IEEEeqnarray*}{rCl} %\IEEEyesnumber \label{eq:tsKalmanPredictionPDFs} \pdfBr{ \given{ \stateV_\indexS }{ \measV_{ 1:\indexS-1 } }} &=& \int{ \pdfBr{ \given{\stateV_\indexS, \stateV_{\indexS-1} }{ \measV_{1:\indexS-1} }} \dif \stateV_{\indexS-1} } = \\ &=& \int{ \pdfBr{ \given{\stateV_\indexS }{ \stateV_{\indexS-1}} } \pdfBr{ \given{\stateV_{\indexS-1} }{ \measV_{1:\indexS-1}} } \dif \stateV_{\indexS-1} } = \\ &=& \int{ \pdfNBr{ \given{\stateV_\indexS }{ \stateM \stateV_{\indexS-1} + \inputM \inputV_{\indexS-1}, \procNCM_{\indexS-1} } } \pdfNBr{ \given{\stateV_{\indexS-1} }{ \est{\stateV}_{\given{\indexS-1}{\indexS-1}}, \stateCM_{\given{\indexS-1}{\indexS-1}} }} \dif \stateV_{\indexS-1} } = \\ &=& \pdfNBr{ \given{\stateV_{\indexS}}{\est{\stateV}_{\given{\indexS}{\indexS-1}}, \stateCM_{\given{\indexS}{\indexS-1}}} }. \end{IEEEeqnarray*} Analog gilt für die Messung $\measV_\indexS$ \begin{subequations} \label{eq:tsGKFPredictMeasurement} \begin{IEEEeqnarray}{rCl} \pdfBr{ \given{\measV_{\indexS}} {\measV_{1:\indexS-1}} } &=& \int{ \pdfBr{ \given{\measV_\indexS, \stateV_\indexS}{\measV_{1:\indexS-1}} } \dif \stateV_\indexS } \IEEEyessubnumber* \\ &=& \int{ \pdfBr{ \given{\measV_\indexS }{ \stateV_{\indexS}} } \pdfBr{ \given{\stateV_{\indexS} }{ \measV_{1:\indexS-1}} } \dif \stateV_\indexS}\\ &=& \int{ \pdfNBr{ \given{\measV_\indexS }{ \outputM \stateV_{\indexS} + \feedM \inputV_{\indexS}, \measNCM } } \pdfNBr{ \given{\stateV_{\indexS} }{ \est{\stateV}_{\given{\indexS}{\indexS-1}}, \stateCM_{\given{\indexS}{\indexS-1}} }} \dif \stateV_{\indexS} } = \\ &=& \pdfNBr{ \given{\measV_{\indexS}} {\est{\measV}_{\given{\indexS} {\indexS-1} } }, \measCM_{\indexS} }. \end{IEEEeqnarray} \end{subequations} % Die benötigten Prädiktionsgrößen $\est{\stateV}_{\given{\indexS}{\indexS-1}}$ und $\stateCM_{\given{\indexS}{\indexS-1}}$ für den Zustand sowie $\est{\measV}_{\given{\indexS}{\indexS-1}}$ und $\measCM_\indexS$ für die Messung liefert der Prädiktionsschritt des Kalman-Filters (\algoref{alg:tsKalmanPrediction}). Für das Residuum $\resV_\indexS = \measV_\indexS - \est{\measV}_{\given{\indexS} {\indexS-1} }$ gilt entsprechend \begin{IEEEeqnarray}{c} \pdfBr{\given{\resV_\indexS}{\measV_{1:\indexS-1}}} = \pdfNBr{ \given{\resV_{\indexS}} {\nullV, \measCM_{\indexS} }}. \label{eq:tsGKFResidualNormalPDF} \end{IEEEeqnarray} Für jede normalverteilte Zufallsvariable $\rvaV \in \dimR{\nDim}$, $\rvaV \sim \pdfNBr{\meanVSym, \covMSym}$ mit Mittelwert $\meanVSym$ und positiv definiter Kovarianzmatrix $\covMSym$ hat die quadratische Form \begin{IEEEeqnarray*}{c} \SNPE(\rvaV) = \trans{\left(\rvaV - \meanVSym \right)} \inv{\covMSym} \left(\rvaV - \meanVSym \right) \end{IEEEeqnarray*} eine Chi-Quadrat-Verteilung mit $\nDim$ Freiheitsgraden \cite[][57-59]{BarShalom2001}: \begin{IEEEeqnarray*}{c} \pdfBr{\SNPE(\rvaV)} = \pdfChiSqBr{\nDim}{\SNPE(\rvaV)}. \end{IEEEeqnarray*} Aufgrund von \eqref{eq:tsGKFResidualNormalPDF} gilt daher beim Kalman-Filter für das Residuum $\resV_\indexS$, dass \begin{IEEEeqnarray}{c} \SNPE(\resV_\indexS) = \trans{\resV}_\indexS \inv{\measCM}_\indexS \resV_\indexS \sim \pdfChiSqBr{\nMeasV}{\SNPE(\resV_\indexS)} \label{eq:tsGKFMahalanobisSquared} \end{IEEEeqnarray} eine Chi-Quadrat-Verteilung aufweist \cite[][236]{BarShalom2001}. % sec. 5.4.2 $\sqrt{\SNPE(\resV_\indexS)}$ ist auch als Mahalanobis-Distanz bekannt, und konstante Werte von $\SNPE(\resV_\indexS)$ entsprechen Hyperellipsoiden mit konstanter Wahrscheinlichkeitsdichte $\pdfBr{\given{\resV_\indexS}{\measV_{1:\indexS-1}}}$ \cite[][78-81]{Bishop2006}. %[BarShalom2001, sec. 3.7.4, p. 166] confidence region, "g-sigma" ellipsoid, probability concentration ellipsoid -> connection to Mahalanobis distance Eine Visualisierung dieser Zusammenhänge zeigt \figref{fig:tsMahalanobis} für ein Beispiel mit zweidimensionalem Merkmalsvektor $\measV_\indexS$. Alle Werte von $\resV_\indexS$ im schattierten Bereich von \fullsubfigref{fig:tsMahalanobisNormalContours} treten mit gleicher oder höherer Wahrscheinlichkeit auf wie die eingetragene beispielhafte Realisierung von $\measV_\indexS$. Die schattierten Flächen in \fullsubfigref{fig:tsMahalanobisNormalContours} und \fullsubfigref{fig:tsMahalanobisChiSquaredPDF} entsprechen daher beide dem selben Wert der Verteilungsfunktion in \fullsubfigref{fig:tsMahalanobisChiSquaredCDF}. \texttodo{Referenz [Roth2013, p. 7]} \begin{figure} \centering \footnotesize \subcaptionbox{Konturen der Wahrscheinlichkeitsdichtefunktion $\pdfBr{\given{\resV_\indexS}{\measV_{1:\indexS-1} }}$ und eine beispielhafte Realisierung \starWithParens{myMediumBlue} von $\resV_\indexS$. \label{fig:tsMahalanobisNormalContours}}[\columnwidth]{% \centering \footnotesize \setlength\figurewidth{.35\columnwidth} \setlength\figureheight{.35\columnwidth} \renewcommand\figureXLabel{$\resV_{\indexS,1}$} \renewcommand\figureYLabel{$\resV_{\indexS,2}$} \mytikzexton \input{"graphics/P9021 01/001 Normal_PDF_Contours.tikz"} \mytikzextoff }% % --------------- comment magic (tm) \\% \vspace{4mm}% % --------------- comment magic (tm) \subcaptionbox{Wahrscheinlichkeitsdichtefunktion von $\SNPE(\resV_\indexS)$ in Form einer Chi-Quadrat-Verteilung \lineWithParens{black} und deren Wert für die beispielhafte Realisierung \starWithParens{myMediumBlue} von $\resV_\indexS$. \label{fig:tsMahalanobisChiSquaredPDF}}[\columnwidth]{% \centering \footnotesize \setlength\figurewidth{.35\columnwidth} \setlength\figureheight{.25\columnwidth} \renewcommand\figureXLabel{$\SNPE(\resV_\indexS)$} \renewcommand\figureYLabel{$\pdfBr{\given{\SNPE(\resV_\indexS)}{\measV_{1:\indexS-1}}}$} \mytikzexton \input{"graphics/P9021 01/002 Chi_squared_PDF.tikz"} \mytikzextoff }% % --------------- comment magic (tm) \\% \vspace{4mm}% % --------------- comment magic (tm) \subcaptionbox{Verteilungsfunktion von $\SNPE(\resV_\indexS)$ und deren Wert für die beispielhafte Realisierung \starWithParens{myMediumBlue} von $\resV_\indexS$. \label{fig:tsMahalanobisChiSquaredCDF}}[\columnwidth]{% \centering \footnotesize \setlength\figurewidth{.35\columnwidth} \setlength\figureheight{.25\columnwidth} \renewcommand\figureXLabel{$\SNPE(\resV_\indexS)$} \renewcommand\figureYLabel{$\cdfBr{\given{\SNPE(\resV_\indexS)}{\measV_{1:\indexS-1}}}$} \mytikzexton \input{"graphics/P9021 01/003 Chi_squared_CDF.tikz"} \mytikzextoff }% \caption[FIXME]{Darstellung der Verteilung von $\resV_\indexS$ und $\SNPE(\resV_\indexS)$ beim Kalman-Filter im zweidimensionalen Fall $\nMeasV = \num{2}$ mit $\measCM_\indexS = \mmbegin{3 & 1\\1 & 2}$.} \label{fig:tsMahalanobis} \end{figure} Praktisch nutzt man beim \ac{GKF} die beschriebenen Eigenschaften der Residuen, um extreme Realisierungen von $\measV_\indexS$ anhand ihrer entsprechend großen Werte von $\SNPE(\resV_\indexS)$ zu erkennen und bei den Filter-Updates auszuschließen. Dabei berücksichtigt die Skalierung mit der inversen Kovarianzmatrix $\inv{\measCM}_\indexS$ die Unsicherheit des prädizierten Merkmalsvektors. Ein Schwellwert $\SNPE_\gate$ (Validation-Gate) legt fest, mit welcher Wahrscheinlichkeit $\Prob_\gate$ der Merkmalsvektor für den Update-Schritt akzeptiert wird: %\begin{subequations} \label{eq:tsGKFCDF} \begin{IEEEeqnarray}{c} \Prob_\gate = \ProbBr{\SNPE(\resV_\indexS) \leq \SNPE_\gate} = \cdfBr{\SNPE_\gate}, \IEEEnonumber \\ \Rightarrow \SNPE_\gate = \invCdfBr{\Prob_\gate}. \IEEEyesnumber \label{eq:tsGKFCDFInv} \end{IEEEeqnarray} %\end{subequations} Der Schwellwert $\SNPE_\gate$ wird also benutzerdefiniert mit Hilfe der inversen Verteilungsfunktion $\invCdf$ der Chi-Quadrat-Verteilung auf Basis von $\Prob_\gate$ festgelegt. Ein gängiger Wert ist z.B. $\Prob_\gate = \num{0.9973}$, welcher im eindimensionalen Fall $\nMeasV = 1$ mit $\delta_\gate = 3^2$ der dreifachen Standardabweichung der Residuen entspricht. Der Ausschluss von erkannten Ausreißern mit $\SNPE(\resV_\indexS) > \SNPE_\gate$ aus dem Update wird durch einen Wert des Kalman-Gains von $\kalmM_\indexS = \nullM$ realisiert. Innerhalb des Validation-Gate wird der Gain des normalen Kalman-Filters verwendet: \begin{IEEEeqnarray}{c} \kalmM_{\indexS} = \begin{cases} \stateCM_{\given{\indexS}{\indexS-1}} \trans{\outputM} \inv{\measCM}_{\indexS} & \text{wenn } \SNPE(\resV_\indexS) \leq \SNPE_\gate, \\ \nullM & \text{sonst}. \end{cases} \label{eq:tsGKFGainCases} \end{IEEEeqnarray} \Algoref{alg:tsGKFUpdate} zeigt den vollständigen Algorithmus des Update-Schritts beim \ac{GKF}, der direkt aus dem gewöhnlichen Kalman-Filter in \algoref{alg:tsKalmanUpdate} hervorgeht. Der Prädiktionsschritt erfolgt wie bekannt mittels \algoref{alg:tsKalmanPrediction}. \begin{algorithm}[htb] \caption{Update-Schritt des \acl{GKF}s (\acs{GKF})} \label{alg:tsGKFUpdate} \begin{algorithmic}[1] \setstretch{1.2} \State $ \resV_{\indexS} \gets \measV_\indexS - \est{\measV}_{\given{\indexS} {\indexS-1} } $ \State $ \SNPE(\resV_\indexS) \gets \trans{\resV}_\indexS \inv{\measCM}_\indexS \resV_\indexS $ \If {$\SNPE(\resV_\indexS) \leq \SNPE_\gate$} \State $ \kalmM_{\indexS} \gets \stateCM_{\given{\indexS}{\indexS-1}} \trans{\outputM} \inv{\measCM}_{\indexS} $ \Else \State $ \kalmM_{\indexS} \gets \nullM $ \EndIf \State $ \est{\stateV}_{\given{\indexS}{\indexS}} \gets \est{\stateV}_{\given{\indexS}{\indexS-1}} + \kalmM_{\indexS} \resV_{\indexS} $ \State $ \stateCM_{\given{\indexS}{\indexS}} \gets \stateCM_{\given{\indexS}{\indexS-1}} - \kalmM_{\indexS} \measCM_{\indexS} \trans{\kalmM}_{\indexS} $ \end{algorithmic} \end{algorithm} Das schaltende Verhalten von $\kalmM_\indexS$ im GKF führt zu einer Unstetigkeit der Schätzergebnisse bei $\SNPE(\resV_\indexS) = \SNPE_\gate$. Dies wird sich später als nachteilig für die Parameterschätzung und die Behandlung von sprungförmigen Änderungen in den Merkmalswerten erweisen. %============================================================================== \subsection{Zustandsraummodell mit Student-t-Messrauschen} %============================================================================== Die nachfolgenden Abschnitte stellen Zustandsraummodell auf Basis nicht-normalverteilter Rauschterme und darauf aufbauende robuste Filter vor, die das schaltende Verhalten des \ac{GKF} vermeiden. Bezüglich der Verteilungsannahmen im Zustandsraummodell \eqref{eq:tsLGSS} besteht sowohl beim Prozessrauschen $\procNV_\indexS$ als auch beim Messrauschen $\measNV_\indexS$ die Möglichkeit, von der Normalverteilungsannahme abzuweichen und so entsprechende Ausreißer zu berücksichtigen. Ausreißer im Prozessrauschterm sind in der Literatur als Innovations-Ausreißer (Innovation Outlier) bekannt, Ausreißer des Messrauschens hingegen als additive Ausreißer (Additive Outlier) \cite{Martin1987}. %Def. Kategorisierung Messung/Innovation auch in [Muirhead1986] mit Verweis auf Zeitreihen-Ausreißern in Messung/Innovation in [Fox1972, p. 350] und Quelle (Burman, 1965) darin Wie die beiden Rauschterme haben auch deren extreme Realisierungen in Form von Ausreißern eine unterschiedliche Wirkung auf das Modellverhalten. Ein Ausreißer im Prozessrauschen zum Zeitpunkt $\indexS$ modelliert über \eqref{eq:tsLGSSStateTransition} eine nachhaltige Veränderung des Zustandsvektors ab $\stateV_{\indexS+1}$. Ein Ausreißer im Messrauschen zum Zeitpunkt $\indexS$ hingegen hat nur Einfluss auf die Messung $\measV_\indexS$, der Zustandsvektor und damit auch zukünftige Messungen bleiben davon unberührt. %Wie sich beim Vergleich des Kalman-Filters mit den robusten Varianten zeigen wird, gilt dieser auf die jeweilige Messung beschränkte Einfluss von Ausreißern im Messrauschen jedoch nicht für die Schätzgrößen. Als Basis für die nachfolgenden Filteralgorithmen gelten die Modellgleichungen \begin{subequations} \label{eq:tsHeavySS} \begin{IEEEeqnarray}{rCl} \stateV_{\indexS+1} &=& \stateM \stateV_{\indexS} + \inputM \inputV_{\indexS} + \procNV_{\indexS}, \IEEEyessubnumber* \\ \measV_{\indexS} &=& \outputM \stateV_{\indexS} + \feedM \inputV_{\indexS} + \measNV_{\indexS}, \label{ts:HeavySSMeasV} \\ \procNV_\indexS &\sim& \pdfNBr{\nullV, \procNCM_\indexS}, \\ \measNV_\indexS &\sim& \pdftBr{\nullV, \measNCM, \dofMeasN}. \label{eq:tsHeavySSMeasN} \end{IEEEeqnarray} \end{subequations} Dabei wurde die Verteilung des Messrauschens $\measNV_\indexS$ durch eine Student-t-Verteilung als konkrete Form einer Heavy-Tail-Verteilung ersetzt. Die Eigenschaften der Student-t-Verteilung sind bereits aus \chref{ch:basics} bekannt. Ziel ist es, eine Robustheit der Filter gegen einzelne Ausreißer in Form von \ac{niO}-Prüfläufen zu erreichen. Auf eine Modellierung des Prozessrauschens mit einer Heavy-Tail-Verteilung wurde hingegen bewusst verzichtet. Der Grund dafür ist, dass ansonsten im Filtervorgang erst mit einigen Samples Verzögerung festgestellt werden könnte, ob ein Ausreißer dem Prozess- oder Messrauschen zuzuordnen ist (vgl. Diskussion zur \glqq{}Outlier Confusion\grqq{} in \cite{Meinhold1989}). Die Aussagen zur Berücksichtigung unregelmäßiger Zeitabstände aus \secref{sec:tsTimeDistUnequal} gelten für \eqref{eq:tsHeavySS} unverändert, da sich die Annahmen über den Zustand $\stateV_\indexS$ und das Prozessrauschen $\procNV_\indexS$ nicht geändert haben. Für $\dofMeasN \rightarrow \infty$ geht \eqref{eq:tsHeavySS} in das linear-gaußsche Zustandsraummodell \eqref{eq:tsLGSS} über. %Outlier Confusion: %Einige Quellen in [Agamennoni2012, p. 5028] \glqq{}The moment of indecision\grqq wenn Prior+Likelihood heavy tailed; auch ohne Student-t-Prozessrauschen! %[Meinhold1989] Outlier confusion entsteht wenn Prior (Prädiktion) und Likelihood (Messung) beide heavy-tailed sind (siehe Ende sec. 3) -> posterior "converges to neither the prior nor the likelihood" %[Fahrmeir1999, p. 174]: "Outlier confusion" bei gleichzeitiger Suche nach Measurement und Innovation Outliern %============================================================================== \subsection{Einordnung der Zustandsschätzer für nicht-gaußsche Verteilungsannahmen} %============================================================================== Die Aufgabe der stochastischen Zustandsschätzung in Form des Filter- bzw. Glättungsproblems löst das Kalman-Filter bzw. der \ac{RTS}-Algorithmus für das linear-gaußsche Zustandsraummodell gleichermaßen recheneffizient als auch analytisch exakt. Sowohl die Einführung von Nichtlinearitäten beim Zustandsübergang oder der Messung als auch die Annahme von nicht-normalverteilten Rauschtermen führt allerdings dazu, dass eine exakte Lösung des Filter- und Glättungsproblems im Allgemeinen nicht mehr möglich ist \cite[][644]{Bishop2006}. Beide Fälle erfordern daher den Einsatz von Näherungslösungen. Für nichtlineare Modellgleichungen bei Erhalt der Annahme normalverteilter Rauschterme und Schätzgrößen haben sich das Extended-Kalman-Filter \cite{Smith1962} und das Unscented-Kalman-Filter \cite{Julier1995, Julier2004} als approximative Schätzverfahren etabliert. %[Agamennoni2011]: \glqq{}Unfortunately, any distribution other than the Gaussian renders the filter analytically intractable and requires some form of approximation.\grqq %Beschreibung von Problemen bei Verletzung der Normalverteilungsannahmen (Sprünge, Ausreißer): [Aravkin2016, p. 3] und Quellen 83, 85, 86 Weichen hingegen die Modellannahmen für die Rauschterme (auch bei weiterhin linearen Modellgleichungen) von der Normalverteilung ab, dann ist das Kalman-Filter zwar weiterhin der \emph{lineare} Schätzer mit der geringsten Varianz der Schätzergebnisse \cite[][3]{Aravkin2016}. Die Linearität des Schätzers bezieht sich hierbei auf den Einfluss des Residuums im jeweiligen Update-Schritt. Beim Kalman-Filter geht das Residuum über die Verstärkung immer linear in das Update des Zustandsvektors ein. Durch verschiedene Näherungsverfahren lassen sich allerdings an die Modelleigenschaften angepasste \emph{nichtlineare} Schätzverfahren herleiten, die an nicht-gaußsche Verteilungsannahmen besser angepasst sind und zu diesem Zweck nichtlineare Zustands-Updates nutzen. In der Literatur existieren zahlreiche Näherungsverfahren für Filter mit nicht-normalverteilten Rauschtermen. Ein früher Ansatz ist das Gaussian-Sum-Filter \cite{Alspach1972}, bei dem die Wahrscheinlichkeitsdichtefunktionen durch gewichtete Summen von Normalverteilungen genähert werden. Die Anzahl der Komponenten nimmt dabei exponentiell über die Zeit zu, sodass eine Heuristik zur Reduktion der Anzahl der Komponenten nötig ist. Eine Näherung durch gewichtete Student-t-Verteilungen beschreiben Meinhold und Singpurwalla \cite{Meinhold1989}. Einen Ansatz zur numerischen Integration hat z.B. Kitagawa \cite{Kitagawa1987} über stückweise lineare Näherungen der Wahrscheinlichkeitsdichtefunktionen aufgestellt. \texttodo{Weitere Referenzen einfügen aus [Roth2017c, p. 1f]} Monte-Carlo-Lösungen mit Varianten des sogenannten Particle-Filters behandeln zahlreiche Veröffentlichungen seit der Vorstellung durch Gordon \cite{Gordon1993}. Anstatt durch Näherungsformeln werden die Wahrscheinlichkeitsdichten hier durch eine große Anzahl an diskreten Samples (Partikel) dargestellt. Einen Überblick über diese auch als Sequential-Monte-Carlo bekannten Verfahren geben Doucet et al. \cite{Doucet2000, Andrieu2010}. Die Implementierung einer robusten Zustands- und Parameterschätzung auf Basis eines Particle-Filters findet sich z.B. bei Schön et al. \cite{Schon2011b}. Die Particle-Filter weisen den Vorteil auf, dass durch Erhöhung des Rechenaufwands (Anzahl der Partikel) die Filterergebnisse beliebig genau genähert werden können. Die tatsächlich benötigte Anzahl an Partikeln kann dabei das in der Anwendung tolerierbare Maß an Rechenaufwand übersteigen \cite{Agamennoni2012}. Schwächen zeigen die Particle-Filter bei Modellen mit geringem Prozessrauschen im Vergleich zum Messrauschen \cite[][283]{Storvik2002}. Wie sich qualitativ bereits am langsamen Driftvorgang in \figref{fig:tsExampleDrift} erkennen lässt und genauer in \secref{sec:tsParamEst} aufgezeigt wird, ist dieser Fall für die \ac{EOL}-Prüfung besonders relevant. Außerdem liegt wie bei allen Monte-Carlo-Verfahren nur ein kleiner Teil der Partikel in den für die Bestimmung von Prüfgrenzen interessanten Randbereichen der Wahrscheinlichkeitsdichtefunktionen, sodass der Näherungsfehler hier relativ groß ist \cite[][Abschnitt 7.5.2]{Kruschke2015}. %Gerade diese Randbereiche sind aber für die Festlegung der Prüfgrenzen in \secref{sec:tsTestLimits} wichtig, da hier der Übergang vom \ac{iO}- zum \ac{niO}-Bereich stattfindet. %Agamennoni2012: "prohibitively expensive in high-dimensional problems"; "actual size of the sample required to achieve desired degree of accuracy may be prohibitively large"; "hard to assess the convergence and the reliability of the estimates even for the most carefully engineered algorithm" Deshalb werden Particle-Filter nachfolgend nicht weiter betrachtet, sondern analytische Näherungsverfahren herangezogen. %Vergleichende Auflistung von robusten Filter/Smoothern in [Agamennoni2011] und [Agamennoni2012] %Kurze Auflistung Möglichkeiten Behandlung nicht-Gaussche Modelle in [Agamennoni2011, p. 1551]: %* one-dimensional methods that are not extendable %* ad-hoc cost functions, involved, require hand tuning (meint: M-estimators) %* numerical methods: numerical integration; nonparametric: particle filter %============================================================================== \subsection{\acf{TKF}} %============================================================================== Roth et al. haben in ihren Veröffentlichungen ein robustes stochastisches Filter vorgestellt, das sowohl das Prozess- als auch das Messrauschen mit Student-t-Verteilungen modelliert \cite{Roth2013b, Roth2017c}. Es wird nachfolgend als \ac{TKF} bezeichnet. Dabei verfolgen sie den Ansatz des \ac{ADF}, bei dem die A-posteriori-Wahrscheinlichkeitsdichte $\pdfBr{ \given{ \stateV_\indexS }{ \measV_{1:\indexS } }}$ nach jedem Update-Schritt durch eine Dichte $\pdfqBr{\stateV_\indexS}$ aus einer vorgegebenen Familie von Verteilungen genähert wird. Konkret kann beispielsweise angenommen werden, dass $\pdfqBr{\stateV_\indexS}$ eine Normalverteilung aufweist, auch wenn der exakt ausgeführte Update-Schritt zunächst nicht in einer Normalverteilung von $\stateV_\indexS$ resultiert. Die Näherung erfolgt dann durch Anpassen der Parameter von $\pdfqBr{\stateV_\indexS}$ an $\pdfBr{ \given{ \stateV_\indexS }{ \measV_{1:\indexS } }}$ \cite[][652\psq]{Murphy2013}. Zweck der Näherungen ist es, die analytische Lösbarkeit der rekursiven Filtergleichungen aufrecht zu erhalten. %[Murphy2013, p. 652f] ... Miminize the Kullback-Leibler divergence % % Bishop2006, hier leider wenig hilfreich da primär auf Batch bezogen: %If q is in the exponential family, one can show that this KL minimization can be done by moment matching. => siehe [Bishop2006, p. 505f] (Achtung min KL(p||q) ) %[Bishop2006, p. 510]: ADF is a special form of Expectation Propagation! (ADF: nur ein Durchlauf durch den Datensatz (Filtern), EP: Batch); Details zu ADF -> EP in [Minka2001, Minka2001c] % %[Minka2001] eigener Abschnitt zu Assumed Density Filtering; siehe auch Diss [Minka2001c] % %Gaussian assumed density filter (ADF; "General Gaussian filtering", Moment Matching) [Särkkä2013, p. 96 ff] = Verallgemeinerung EKF, UKF %Herleitung hier für Nichtlineare Systeme %General Gaussian smoothing (= Assumed Density Smoothing?) [Särkkä2013, p. 154ff] %[Särkkä2013, p. 192] im Rahmen Parameterschätzung: Gaussian Filtering/Smoothing = Approximation with Gaussian Moment Matching %----------------------------------------------------------------------------- \subsubsection{Prädiktions- und Update-Schritt mit gemeinsamer Student-t-Verteilung} %----------------------------------------------------------------------------- Die nachfolgenden Herleitungen für das \ac{TKF} basieren auf \cite{Roth2013b} und wurden um den Systemeingang $\inputV_\indexS$ erweitert, um zu den anderen Abschnitten dieses Kapitels konsistent zu sein. Konkret werden beim \ac{TKF} die A-posteriori-Verteilungen nicht als Normalverteilungen, sondern als Student-t-Verteilungen angenommen. Es weisen also nicht nur die Rauschterme $\procNV_\indexS$ und $\measNV_\indexS$ eine Student-t-Verteilung auf, sondern auch die resultierenden Schätzungen des Zustandsvektors $\stateV_\indexS$. %Zwischen den einzelnen Filterschritten werden die Student-t-Verteilungen durch andere Student-t-Verteilungen genähert, um die Freiheitsgrade konstant zu halten. Die nachfolgenden Prädiktions- und Update-Schritte arbeiten mit den aus dem linear-gaußschen Zustandsraummodell bekannten Modellgleichungen \begin{IEEEeqnarray*}{rCl} \stateV_{\indexS+1} &=& \stateM \stateV_{\indexS} + \inputM \inputV_\indexS + \procNV_{\indexS}, \\ \measV_{\indexS} &=& \outputM \stateV_{\indexS} + \feedM \inputV_\indexS + \measNV_{\indexS}. \end{IEEEeqnarray*} Der erste wesentliche Schritt bei der Herleitung des Prädiktions-Schrittes beim TKF ist die Annahme einer gemeinsamen Student-t-Verteilung des geschätzten Zustandsvektors $\stateV_\indexS$ und des Prozessrauschens $\procNV_\indexS$ mit gleicher Anzahl der Freiheitsgrade $\dofPred$, \begin{IEEEeqnarray}{c} \pdfBr{\given{ \stateV_\indexS, \procNV_\indexS }{ \measV_{1:\indexS }}} = \pdftBr{ \givenBigThree{ \mvbegin{\stateV_\indexS \\ \procNV_\indexS} }{ \mvbegin{\est{\stateV}_{\given{ \indexS }{ \indexS }} \\ \nullV }, \mmbegin{ \stateCM_{\given{ \indexS }{ \indexS }} & \nullM \\ \nullM & \procNCM_\indexS}, \dofPred }} \label{eq:tsTKFPredictionJointPDF} \end{IEEEeqnarray} Die Eigenschaften der Student-t-Verteilung resultieren darin, dass $\stateV_\indexS$ und $\procNV_\indexS$ in \eqref{eq:tsTKFPredictionJointPDF} zwar unkorreliert, aber nicht unabhängig sind \cite{Roth2013b}. %Unter Annahme eines gemeinsamen Parameters $\dofPred$ für die Anzahl der Freiheitsgrade. Durch die Annahme einer gemeinsamen Verteilung und die Eigenschaften der Student-t-Verteilung ist es möglich, die lineare Transformation \begin{IEEEeqnarray*}{c} \stateV_{\indexS+1} = \mmbegin{\stateM & \nullM \\ \nullM & \eye} \cdot \mvbegin{\stateV_\indexS \\ \procNV_\indexS} + \inputM \inputV_\indexS \end{IEEEeqnarray*} auszuwerten, die wiederum eine Student-t-Verteilung aufweist: \begin{subequations} \label{eq:tsTKFStatePrediction} \begin{IEEEeqnarray}{rCl} \pdfBr{\given{ \stateV_{\indexS+1} }{ \measV_{1:\indexS }}} &=& \pdftBr{\given{ \stateV_{\indexS+1} }{ \est{\stateV}_{\given{ \indexS+1 }{ \indexS }}, \stateCM_{\given{ \indexS+1 }{ \indexS }} , \dofPred} }, \IEEEyessubnumber* \\ \est{\stateV}_{\given{ \indexS+1 }{ \indexS }} &=& \stateM \est{\stateV}_{\given{ \indexS }{ \indexS }} + \inputM \inputV_\indexS, \label{eq:tsTKFStateVPrediction}\\ \stateCM_{\given{ \indexS+1 }{ \indexS }} &=& \stateM \stateCM_{\given{ \indexS }{ \indexS }} \trans{\stateM} + \procNCM_\indexS . \label{eq:tsTKFStateCMPrediction} \end{IEEEeqnarray} \end{subequations} Die Berechnungsvorschriften \eqref{eq:tsTKFStateVPrediction} und \eqref{eq:tsTKFStateCMPrediction} für die Parameter der Zustandsprädiktion sind dabei exakt gleich wie beim Kalman-Filter. Um den Update-Schritt herzuleiten wird angenommen, dass die Prädiktion des Zustandsvektors und das Messrauschen eine gemeinsame Student-t-Verteilung mit gleicher Anzahl der Freiheitsgrade $\dofUpd$ aufweisen: \begin{IEEEeqnarray}{c} \pdfBr{\given{ \stateV_\indexS, \measNV_\indexS }{ \measV_{1:\indexS-1 }}} = \pdftBr{ \givenBigThree{ \mvbegin{\stateV_\indexS \\ \measNV_\indexS} }{ \mvbegin{\est{\stateV}_{\given{ \indexS }{ \indexS-1 }} \\ \nullV }, \mmbegin{ \stateCM_{\given{ \indexS }{ \indexS-1 }} & \nullM \\ \nullM & \measNCM}, \dofUpd }} \label{eq:tsTKFUpdateJointAnsatz} \end{IEEEeqnarray} Wie bereits beim Prädiktionsschritt wird nun eine lineare Transformation angewendet: \begin{IEEEeqnarray*}{c} \mvbegin{\stateV_\indexS \\ \measV_\indexS} = \mmbegin{ \eye & \nullM \\ \outputM & \eye } \cdot \mvbegin{\stateV_\indexS \\ \measNV_\indexS} + \mmbegin{\nullM \\ \feedM} \cdot \inputV_\indexS . \end{IEEEeqnarray*} Die gemeinsame Verteilung von $\stateV_\indexS$ und $\measV_\indexS$ ergibt sich aus dieser linearen Transformation zu \begin{subequations} \label{eq:tsTKFUpdateJoint} \begin{IEEEeqnarray}{rCl} \pdfBr{\given{ \stateV_\indexS, \measV_\indexS }{ \measV_{1:\indexS-1} }} %&=& \pdfBr{\givenBigThree{ \mvbegin{\stateV_\indexS \\ \measV_\indexS} }{ \measV_{1:\indexS-1} } } = \IEEEnonumber \\ &=& \pdftBr{ \givenBigThree{ \mvbegin{\stateV_\indexS \\ \measV_\indexS} }{ \mvbegin{ \est{\stateV}_{\given{ \indexS }{ \indexS-1 }} \\ \est{\measV}_{\given{ \indexS }{ \indexS-1 }} }, \mmbegin{ \stateCM_{\given{ \indexS }{ \indexS-1 }} & \stateCM_{\given{ \indexS }{ \indexS-1 }} \trans{\outputM} \\ \outputM \stateCM_{\given{ \indexS }{ \indexS-1 }} & \measCM_\indexS }, \dofUpd }}, \IEEEyessubnumber* \label{eq:tsTKFUpdateJointPDF} \\ \est{\measV}_{\given{ \indexS }{ \indexS-1 }} &=& \outputM \est{\stateV}_{\given{ \indexS }{ \indexS-1 }} + \feedM \inputV_\indexS, \label{eq:tsTKFUpdateJointMeasV} \\ \measCM_\indexS &=& \outputM \stateCM_{\given{ \indexS }{ \indexS-1 }} \trans{\outputM} + \measNCM. \label{eq:tsTKFUpdateJointMeasCM} \end{IEEEeqnarray} \end{subequations} Die Verteilungsparameter in \eqref{eq:tsTKFUpdateJointMeasV} und \eqref{eq:tsTKFUpdateJointMeasCM} stimmen wieder mit den aus dem Kalman-Filter bekannten Berechnungen überein. Der wesentliche Unterschied zum Kalman-Filter ergibt sich nun, wenn aus der gemeinsamen Verteilung \eqref{eq:tsTKFUpdateJointPDF} die bedingte Wahrscheinlichkeitsdichte von $\stateV_\indexS$ berechnet wird. Für die bedingte Dichte gilt aufgrund der gemeinsamen Student-t-Verteilung in \eqref{eq:tsTKFUpdateJointPDF}: \begin{subequations} \label{eq:tsTKFUpdate} \begin{IEEEeqnarray}{rCl} \pdfBr{\given{\stateV_\indexS }{ \measV_{1:\indexS} }} &=& \pdftBr{ \given{ \stateV_\indexS }{ \est{\stateV}_{\given{ \indexS }{ \indexS }}, \stateCM_{\given{ \indexS }{ \indexS }} , \dofUpd^\prime }}, \IEEEyessubnumber* \\ \resV_\indexS &=& \measV_\indexS - \est{\measV}_{\given{ \indexS }{ \indexS-1 }}, \\ \kalmM_\indexS &=& \stateCM_{\given{ \indexS }{ \indexS-1 }} \trans{\outputM} \inv{\measCM}_\indexS, \label{eq:tsTKFUpdateKalmM} \\ \est{\stateV}_{\given{ \indexS }{ \indexS }} &=& \est{\stateV}_{\given{ \indexS }{ \indexS-1 }} + \kalmM_\indexS \resV_\indexS, \label{eq:tsTKFUpdateStateV} \\ \SNPE(\resV_\indexS) &=& \trans{\resV}_\indexS \inv{\measCM}_\indexS \resV_\indexS, \\ \stateCM_{\given{ \indexS }{ \indexS }} &=& \frac{\dofUpd + \SNPE(\resV_\indexS)}{\dofUpd + \nMeasV} \left( \stateCM_{\given{ \indexS }{ \indexS-1 }} - \kalmM_\indexS \outputM \stateCM_{\given{ \indexS }{ \indexS-1 }} \right), \label{eq:tsTKFUpdateStateCM} \\ \dofUpd^\prime &=& \dofUpd + \nMeasV . \label{eq:tsTKFUpdateDoF} \end{IEEEeqnarray} \end{subequations} Hier tritt erstmals ein Unterschied in den Filtergleichungen zum Kalman-Filter auf. Der Kalman-Gain $\kalmM_\indexS$ \eqref{eq:tsTKFUpdateKalmM} und das Update des Zustandsvektors \eqref{eq:tsTKFUpdateStateV} sind unverändert, aber die Kovarianzmatrix \eqref{eq:tsTKFUpdateStateCM} des Zustandsvektors nach dem Updateschritt hängt nun von der quadrierten Mahalanobis-Distanz $\SNPE(\resV_\indexS)$ des Residuums ab. Außerdem hat die A-posteriori-Verteilung von $\stateV_\indexS$ eine höhere Anzahl an Freiheitsgraden $\dofUpd^\prime$ die Prädiktion. Die bisherigen Schritte zur Herleitung des \ac{TKF} waren analytisch exakt, die Annahme gleicher Freiheitsgrade aller stochastischen Komponenten ist aber einschränkend. Die Berechnungen führen außerdem bei jeder Ausführung des Updateschrittes zu einer Zunahme der Anzahl der Freiheitsgrade \eqref{eq:tsTKFUpdateDoF}. Somit ist das Verhalten dieser exakten Form des \ac{TKF} nach einigen hundert Updatesschritten praktisch nicht mehr von dem des gewöhnlichen Kalman-Filters unterscheidbar. %----------------------------------------------------------------------------- \subsubsection{Näherungsweiser Prädiktions- und Update-Schritt} %----------------------------------------------------------------------------- Beide negativen Aspekte der exakten Lösung lassen sich durch Näherungen im Sinne des \ac{ADF} beseitigen. Der initiale Zustandsvektor, das Prozessrauschen und das Messrauschen erhalten hierfür in einem ersten Schritt zunächst jeweils eigenständige Student-t-Verteilungen \begin{subequations} \begin{IEEEeqnarray}{rCl} \pdfBr{\stateV_1} &=& \pdftBr{\given{ \stateV_1 }{ \est{\stateV}_{\given{ 1 }{ 0 }}, \stateCM_{\given{ 1 }{ 0 }}, \dofState }}, \IEEEyessubnumber* \label{eq:tsTKFRealNoiseAssumptionStateV} \\ \pdfBr{\procNV_\indexS} &=& \pdftBr{\given{ \procNV_\indexS }{ \nullV, \procNCM_\indexS, \dofProcN }}, \label{eq:tsTKFRealNoiseAssumptionProcNV} \\ \pdfBr{\measNV_\indexS} &=& \pdftBr{\given{ \measNV_\indexS }{ \nullV, \measNCM, \dofMeasN }}. \label{eq:tsTKFRealNoiseAssumptionMeasNV} \end{IEEEeqnarray} \end{subequations} Für die A-posteriori-Verteilung des Zustandsvektors im Zeitschritt $\indexS$ wird eine Student-t-Verteilung mit $\dofState$ Freiheitsgraden angenommen: \begin{IEEEeqnarray*}{c} \pdfBr{\given{ \stateV_\indexS }{ \measV_{1:\indexS} }} = \pdftBr{\est{\stateV}_{\given{ \indexS }{ \indexS }}, \stateCM_{\given{ \indexS }{ \indexS }}, \dofState}. \end{IEEEeqnarray*} Da \eqref{eq:tsTKFRealNoiseAssumptionStateV} und \eqref{eq:tsTKFRealNoiseAssumptionProcNV} durch die unterschiedlichen Freiheitsgrade $\dofState$ und $\dofProcN$ nicht in einer gemeinsamen Dichtefunktion \eqref{eq:tsTKFPredictionJointPDF} kombiniert werden können, erfolgt vor dem Prädiktionsschritt eine Näherung durch \begin{subequations} \label{eq:tsTKFApproxBeforePrediction} \begin{IEEEeqnarray}{rCl} \pdfBr{\given{\stateV_\indexS }{ \measV_{1:\indexS} }} &\approx& \pdfqBr{\given{\stateV_\indexS }{ \measV_{1:\indexS} }} = \pdftBr{ \given{ \stateV_\indexS }{ \est{\stateV}_{\given{ \indexS }{ \indexS }}, \apr{\stateCM}_{\given{ \indexS }{ \indexS }} , \dofPred }} , \IEEEyessubnumber* \\ \pdfBr{\procNV_\indexS} &\approx& \pdfqBr{\procNV_\indexS} = \pdftBr{\nullV, \apr{\procNCM}_\indexS, \dofPred} . \end{IEEEeqnarray} \end{subequations} Die beiden näherungsweisen Parameter $\apr{\stateCM}_{\given{ \indexS }{ \indexS }}$ und $\apr{\procNCM}_\indexS$ werden in \eqref{eq:tsTKFPredictionJointPDF} eingesetzt und der Prädiktionsschritt gemäß \eqref{eq:tsTKFStatePrediction} ausgeführt. Die resultierende Zustandsprädiktion hat wiederum $\dofPred$ Freiheitsgrade. Für den Update-Schritt erfolgt eine Näherung für die Verteilungen der Zustandsprädiktion und das Messrauschen nach dem gleichen Prinzip mit \begin{subequations} \label{eq:tsTKFApproxBeforeUpdate} \begin{IEEEeqnarray}{rCl} \pdfBr{\given{\stateV_\indexS }{ \measV_{1:\indexS-1} }} &\approx& \pdfqBr{\given{\stateV_\indexS }{ \measV_{1:\indexS-1} }} = \pdftBr{ \given{ \stateV_\indexS }{ \est{\stateV}_{\given{ \indexS }{ \indexS-1 }}, \apr{\stateCM}_{\given{ \indexS }{ \indexS-1 }} , \dofUpd }} , \IEEEyessubnumber* \\ \pdfBr{\measNV_\indexS} &\approx& \pdfqBr{\measNV_\indexS} = \pdftBr{\nullV, \apr{\measNCM}, \dofUpd} . \end{IEEEeqnarray} \end{subequations} Die Parameter $\apr{\stateCM}_{\given{ \indexS }{ \indexS-1 }}$ und $\apr{\measNCM}$ der Näherung \eqref{eq:tsTKFApproxBeforeUpdate} können nun wiederum in den Ansatz \eqref{eq:tsTKFUpdateJointAnsatz} für den Update-Schritt eingesetzt werden. Die A-posteriori-Verteilung von $\stateV_\indexS$ hat $\dofState^\prime = \dofState + \nMeasV$ Freiheitsgrade. Dieser Wert dient nun als neues $\dofState$ für die Näherung \eqref{eq:tsTKFApproxBeforePrediction} vor dem nächsten Prädiktionsschritt, und der Zyklus beginnt von vorne. Die Anzahl der Freiheitsgrade $\dofPred$ der Näherung kann dabei konstant bleiben, sodass ein wiederholter Anstieg der Freiheitsgrade durch die Update-Schritte immer wieder auf $\dofPred$ Freiheitsgrade zurückgeführt wird. \texttodo{Bild einfügen mit Ablauf Näherung - Prädiktion - Näherung - Update - ...} %----------------------------------------------------------------------------- \subsubsection{Realisierung der Näherungen für die \ac{EOL}-Prüfung} %----------------------------------------------------------------------------- Für die Umsetzung des \ac{TKF} stellt sich die Frage, wie die Näherungen \eqref{eq:tsTKFApproxBeforePrediction} und \eqref{eq:tsTKFApproxBeforeUpdate} realisiert werden können. Im Allgemeinen ist die Frage, wie sich für die Zufallsvariable $\rva \in \dimR{\nDim}$ die Verteilung \begin{IEEEeqnarray*}{c} \pdfBr{\rva} = \pdftBr{\given{ \rva }{ \nullV, \covMSym_1, \dofSym_1 }} \end{IEEEeqnarray*} durch eine Verteilung \begin{IEEEeqnarray*}{c} \pdfqBr{\rva} = \pdftBr{\given{ \rva }{ \nullV, \covMSym_2, \dofSym_2 }} \end{IEEEeqnarray*} nähern lässt. Diese Frage wurde bereits in \secref{sec:basicsMathStoch} diskutiert. Roth schlägt vor \cite{Roth2013b}, die Freiheitsgrade $\dofSym_2$ der nähernden Verteilung fest vorzugeben und die Matrix $\covMSym_2$ mittels eines Skalierungsfaktors $\scaleDOF \in \R_{>0}$ aus der Matrix $\covMSym_1$ zu erzeugen, um die Korrelationsstruktur von $\covMSym_1$ auf $\covMSym_2$ zu übertragen. Basierend auf den Erläuterungen in \secref{sec:basicsMathStoch} kommt nachfolgend die Näherung \begin{IEEEeqnarray*}{c} \covMSym_2 = \scaleDOF_{\dofSym_2, \dofSym_1}^2 \cdot \covMSym_1 \end{IEEEeqnarray*} zum Einsatz, und der Skalierungsfaktor wird mittels \begin{IEEEeqnarray}{c} \scaleDOF_{\dofSym_2,\dofSym_1} = \arg \min_{\scaleDOF_{\dofSym_2,\dofSym_1} \in \R_{>0}} \begin{cases} \klBr{\pdf}{\pdfq} & \text{wenn } \dofSym_1 > \dofSym_2, \\ \klBr{\pdfq}{\pdf} & \text{wenn } \dofSym_1 < \dofSym_2 \end{cases} \label{eq:tsTKFRescaleStudentT} \end{IEEEeqnarray} bestimmt \cite{Leitner2018b}. Für die Anwendung des \ac{TKF} im Rahmen der \ac{EOL}-Prüfung gelten entsprechend zum Modell \eqref{eq:tsHeavySS} die Annahmen \begin{IEEEeqnarray*}{c} \dofPred = \dofProcN \rightarrow \infty, \quad \dofUpd = \dofMeasN < \infty. \end{IEEEeqnarray*} Da die Berechnung der Verteilungsparameter im Prädiktionsschritt des \ac{TKF} \eqref{eq:tsTKFStatePrediction} genau wie beim konventionellen Kalman-Filter erfolgt, werden außerdem für die praktische Umsetzung die Näherungen \eqref{eq:tsTKFApproxBeforePrediction} und \eqref{eq:tsTKFApproxBeforeUpdate} im Update-Schritt des \ac{TKF} zusammengezogen. Den resultierenden Algorithmus für den Update-Schritt des \ac{TKF} in der \ac{EOL}-Prüfung zeigt \algoref{alg:tsTKFUpdate}. Als Prädiktionsschritt kommt weiterhin \algoref{alg:tsKalmanPrediction} zum Einsatz. \begin{algorithm}[bth] \caption{Update-Schritt des \acs{TKF} für die \acs{EOL}-Prüfung} \label{alg:tsTKFUpdate} \begin{algorithmic}[1] \setstretch{1.2} \State $ \apr{\stateCM}_{\given{\indexS}{\indexS-1}} \gets \scaleDOF_{\dofUpd,\infty}^2 \cdot \stateCM_{\given{\indexS}{\indexS-1}} $ \label{alg:tsTKFUpdateEstCMOne} \State $ \apr{\measCM}_{\indexS} \gets \outputM \apr{\stateCM}_{\given{\indexS}{\indexS-1}} \outputM^T + \measNCM $ \State $ \resV_{\indexS} \gets \measV_\indexS - \est{\measV}_{\given{\indexS} {\indexS-1} } $ \State $ \apr{\SNPE}_\indexS \gets \resV_{\indexS}^T \apr{\measCM}_{\indexS}^{-1} \resV_{\indexS} $ \State $ \kalmM_{\indexS} \gets \apr{\stateCM}_{\given{\indexS}{\indexS-1}} \outputM^T \apr{\measCM}_{\indexS}^{-1} $ \State $ \apr{\stateCM}_{\given{\indexS}{\indexS}} \gets \frac{\dofUpd + \apr{\SNPE}_{\indexS}}{\dofUpd + \nMeasV} ( \apr{\stateCM}_{\given{\indexS}{\indexS-1}} - \kalmM_{\indexS} \outputM \apr{\stateCM}_{\given{\indexS}{\indexS-1}} ) $ \label{alg:tsTKFUpdateCM} \State $ \est{\stateV}_{\given{\indexS}{\indexS}} \gets \est{\stateV}_{\given{\indexS}{\indexS-1}} + \kalmM_\indexS \resV_\indexS $ \State $ \stateCM_{\given{\indexS}{\indexS}} \gets \scaleDOF_{\infty,\dofUpd+\nMeasV}^2 \cdot \apr{\stateCM}_{\given{\indexS}{\indexS}} $ \label{alg:tsTKFUpdateEstCMTwo} \end{algorithmic} \end{algorithm} Mit konstanten Freiheitsgraden $\dofMeasN$ können außerdem die beiden Skalierungsfaktoren $\scaleDOF_{\dofUpd,\infty}^2$ und $\scaleDOF_{\infty,\dofUpd+\nMeasV}^2$ vorab offline berechnet werden, wodurch der Mehraufwand an Rechenoperationen für das \ac{TKF} gegenüber dem Kalman-Filter im Online-Betrieb vernachlässigbar ist. Abschließend sei zum \ac{TKF} angemerkt, dass für $\dofUpd \rightarrow \infty$ der Update-Schritt in \algoref{alg:tsTKFUpdate} in das gewöhnliche Kalman-Filter übergeht. Dies ist eine erwünschte Eigenschaft, da die Student-t-Verteilung für $\dofSym \rightarrow \infty$ der Normalverteilung entspricht und damit als Grenzfall ein linear-gaußsches Zustandsraummodell vorliegt. %============================================================================== \subsection{\acf{RME}} %============================================================================== Bei dem im vorherigen Abschnitt vorgestellten \ac{TKF} wurden die Wahrscheinlichkeitsdichtefunktionen der Rausch- und Schätzterme bei jedem Prädiktions- und Update-Schritt durch andere Dichtefunktionen genähert, um den eigentlichen Prädiktions- bzw. Updateschritt analytisch exakt berechnen zu können. Einen anderen Ansatz verfolgen rekursive M-Schätzer, bei denen der Update-Schritt mittels eines gewichteten Least-Squares-Problem gelöst werden kann. Der Name M-Schätzer leitet sich von Maximum-Likelihood-Type-Estimator ab, sie wurden ursprünglich von Huber \cite{Huber1964} für statische Regressionsprobleme mit $\numS$ \ac{iid} Datenpunkten $\indexS$, $\indexS \in \msbegin{1, \ldots, \numS}$ eingeführt. Die nachfolgende Herleitung basiert auf einer Arbeit von Brunet \cite[][75-78]{Brunet2010}. Den Ausgangspunkt bilden sogenannte $\MEstFun$-Funktion, auf deren Basis das Schätzproblem für einen gesuchten Parametervektor $\paramV \in \R$ zunächst allgemein als \begin{IEEEeqnarray}{c} \est{\paramV} = \arg \min_{\paramV} \sum_{\indexS=1}^{\numS}{ \MEstFun(\resV_\indexS(\paramV))} \label{eq:tsMEstStart} \end{IEEEeqnarray} formuliert wird. Die $\MEstFun$-Funktion gewichtet dabei die einzelnen Residuen $\resV_\indexS$ und bestimmt so wesentlich das Verhalten des Schätzers. Im einfachsten Fall ist $\MEstFun(\res_\indexS(\paramV)) = \normbegin{\res_\indexS(\paramV)}_2^2$, was in \eqref{eq:tsMEstStart} eingesetzt auf den gewöhnlichen Least-Squares-Schätzer führt. Andere $\MEstFun$-Funktionen ermöglichen es, große Residuen geringer zu gewichten als dies beim Least-Squares-Schätzer der Fall ist und so eine Robustheit gegenüber Ausreißern zu erreichen. Als $\MEstFun$-Funktionen kommen dabei einerseits heuristisch festgelegte Funktionen in Frage. Ein Beispiel mit skalarem Residuum ist der Huber-M-Schätzer \begin{IEEEeqnarray}{c} \MEstFun(\res_\indexS(\paramV)) = \begin{cases} \frac{\res_\indexS(\paramV)^2}{2} & \text{wenn } \abs{\res_\indexS(\paramV)} < \huberParam, \\ \huberParam \cdot \abs{\res_\indexS(\paramV)} - \frac{\huberParam^2}{2} & \text{sonst}. \end{cases} \end{IEEEeqnarray} Der Einstellparameter $\huberParam \in \R_{>0}$ legt hier fest, an welcher Stelle der Übergang von einer quadratischen zu einer weniger auf Ausreißer empfindlichen linearen Bewertung des jeweiligen Residuums erfolgt. Neben einer heuristischen Festlegung kann andererseits über die Beziehung \begin{IEEEeqnarray*}{c} \pdfBr{\given{\resV_\indexS }{ \paramV }} = \expBr{ - \MEstFun(\resV_\indexS(\paramV)) } \end{IEEEeqnarray*} von einer vorgegebenen Likelihood-Funktion $\pdfBr{\given{\resV_\indexS }{ \paramV }}$ eines stochastischen Modells auf eine dazu passende $\MEstFun$-Funktion \begin{IEEEeqnarray}{c} \MEstFun(\resV_\indexS(\paramV)) = - \log \pdfBr{\given{\resV_\indexS }{ \paramV }} \label{eq:tsRMELikelihood} \end{IEEEeqnarray} geschlossen werden. Durch Einsetzen von \eqref{eq:tsRMELikelihood} in \eqref{eq:tsMEstStart} wird der Zusammenhang zum Maximum-Likelihood-Schätzer offensichtlich. Das Optimierungsproblem lautet dann \begin{IEEEeqnarray*}{c} \est{\paramV} = \arg \min_{\paramV} -\sum_{\indexS=1}^{\numS}{ \log \pdfBr{\given{\resV_\indexS }{ \paramV }} } . \end{IEEEeqnarray*} Um das Optimierungsproblem zu lösen, wird wie beim Maximum-Likelihood-Schätzer der Gradient zu minimierenden Funktion gleich Null gesetzt. Bei skalaren Residuen $\res_\indexS$ gilt \begin{IEEEeqnarray*}{rCl} &&\frac{\partial}{\partial \paramV} \sum_{\indexS=1}^{\numS}{ \MEstFun(\res_\indexS(\paramV))} = \nullV \\ &\Leftrightarrow& \sum_{\indexS=1}^{\numS} \frac{\partial \MEstFun(\res_\indexS)}{ \partial \res_\indexS } \frac{\partial \res_\indexS(\paramV)}{\partial \paramV} = \nullV. \end{IEEEeqnarray*} Unter Zuhilfenahme der sogenannten Einflussfunktion (Influence Function) $\MEstInfluenceFun(\res_\indexS)$, \begin{IEEEeqnarray*}{c} \MEstInfluenceFun(\res_\indexS) = \frac{\partial \MEstFun(\res_\indexS)}{ \partial \res_\indexS } \end{IEEEeqnarray*} erfolgt eine weitere Umformung zu \begin{subequations} \begin{IEEEeqnarray}{rCl} && \sum_{\indexS=1}^{\numS} \MEstInfluenceFun(\res_\indexS) \frac{\partial \res_\indexS(\paramV)}{\partial \paramV} = \nullV \IEEEnonumber* \\ &\Leftrightarrow& \sum_{\indexS=1}^{\numS} \frac{\MEstInfluenceFun(\res_\indexS)}{\res_\indexS(\paramV)} \res_\indexS(\paramV) \frac{\partial \res_\indexS(\paramV)}{\partial \paramV} = \nullV \\ &\Leftrightarrow& \sum_{\indexS=1}^{\numS} \MEstWeight(\res_\indexS) \res_\indexS(\paramV) \frac{\partial \res_\indexS(\paramV)}{\partial \paramV} = \nullV \\ &\Leftrightarrow& \sum_{\indexS=1}^{\numS} \MEstWeight(\res_\indexS) \frac{1}{2} \frac{\partial}{\partial \paramV} \res^2_\indexS(\paramV) = \nullV. \IEEEyesnumber \label{eq:tsRMEWeightedM} \end{IEEEeqnarray} \end{subequations} Dabei wurde eine Gewichtungsfunktion $\MEstWeight(\res_\indexS)$ als \begin{IEEEeqnarray}{c} \MEstWeight(\res_\indexS) = \frac{\MEstInfluenceFun(\res_\indexS)}{\res_\indexS(\paramV)} = \frac{1}{\res_\indexS(\paramV)} \frac{\partial \MEstFun(\res_\indexS)}{ \partial \res_\indexS } \label{eq:tsRMEWeightGeneral} \end{IEEEeqnarray} festgelegt. Es fällt nun auf, dass \eqref{eq:tsRMEWeightedM} die gleiche Struktur aufweist wie die notwendige Bedingung zur Lösung des gewichteten Least-Squares-Problems (vgl. \secref{sec:basicsApproachesToEst}): \begin{IEEEeqnarray*}{rCl} && \min_{\paramV} \frac{1}{2} \sum_{\indexS=1}^{\numS}{ \MEstWeight_\indexS \res_\indexS^2 } \\ &\Rightarrow& \sum_{\indexS=1}^{\numS}{ \MEstWeight_\indexS \frac{\partial}{\partial\paramV} \res_\indexS^2 } = \nullV. \end{IEEEeqnarray*} Ein wesentlicher Unterschied zum gewichteten Least-Squares-Problem ist aber, dass in \eqref{eq:tsRMEWeightedM} die Gewichte $\MEstWeight$ von den Residuen $\res_\indexS$ abhängig sind, die selbst wiederum eine Funktion des gesuchten Parametervektors $\paramV$ darstellen. Um \eqref{eq:tsRMEWeightedM} zu lösen, kann daher der \ac{IRLS}-Algorithmus verwendet werden \cite[][44]{Brunet2010}. Dabei erfolgt ausgehend von einer initialen Schätzung $\paramV$ abwechselnd zuerst die Berechnung der Gewichte $\MEstWeight(\res_\indexS)$ und die darauf aufbauende gewichtete Least-Squares-Schätzung, bis das Schätzergebnis konvergiert. Um nun das Filterproblem für die Modellgleichungen \eqref{eq:tsHeavySS} zu lösen, ist erstens die Herleitung der Gewichtungsfunktion $\MEstWeight(\res_\indexS)$ für Residuen mit Student-t-Verteilung erforderlich. Zweitens wird für die Anwendung ein rekursiver Update-Algorithmus benötigt. Beides wird nachfolgend für $\nMeasV = 1$ angeführt. Bezüglich der Gewichtungsfunktion gilt mit \eqref{eq:tsHeavySSMeasN} entsprechend \eqref{eq:tsRMELikelihood}, dass die $\MEstFun$-Funktion \begin{IEEEeqnarray*}{c} \MEstFun(\res_\indexS) = - \log \left( \frac{\GammaFunBr{\frac{\dofMeasN + 1}{2} }}{ \GammaFunBr{ \frac{\dofMeasN}{2} } \sqrt{\dofMeasN \pi \measNC }} \right) + \frac{\dofMeasN+1}{2} \log \left( 1 + \frac{\res_\indexS^2}{\dofMeasN \measNC} \right) . \end{IEEEeqnarray*} ist. Aus dem von $\res_\indexS$ abhängigen Anteil wird mittels \eqref{eq:tsRMEWeightGeneral} die Gewichtungsfunktion \begin{IEEEeqnarray*}{c} \MEstWeight(\res_\indexS) = \frac{1}{\res_\indexS} \frac{2}{1 + \frac{\res_\indexS^2}{\dofMeasN \measNC}} \frac{\res_\indexS}{\dofMeasN \measNC} = \frac{\dofMeasN + 1}{\dofMeasN \measNC + \res_\indexS^2} \end{IEEEeqnarray*} ermittelt. Neben dem M-Schätzer führt der \ac{EM}-Algorithmus für die Student-t-Verteilung auf die gleiche Gewichtungsfunktion \cite[][444\psq]{Gelman2014}. Eine rekursive Variante des M-Schätzers findet sich bei Rhode \cite[][48\psq]{Rhode2016b}. Dort wird der Prädiktions- und Update-Schritt gemeinsam ausgeführt und der \ac{IRLS}-Algorithmus nach einer Iteration beendet. \Algoref{alg:tsRMEUpdate} zeigt einen Vorschlag für den Update-Schritt zur Verwendung gemeinsamen Verwendung mit \algoref{alg:tsKalmanPrediction}. Dabei wurde wie bereits zuvor beim \ac{TKF} darauf geachtet, dass auch der Update-Schritt des \ac{RME} für $\dofMeasN \rightarrow \infty$ in den Update-Schritt des Kalman-Filters übergeht. \begin{algorithm}[bth] \caption{Update-Schritt des \acs{RME} für die \acs{EOL}-Prüfung} \label{alg:tsRMEUpdate} \begin{algorithmic}[1] \setstretch{1.2} \State $ \res_{\indexS} \gets \meas_\indexS - \est{\meas}_{\given{\indexS} {\indexS-1} } $ \State $ \MEstWeight_\indexS \gets \left( \dofMeasN + 1 \right) / \left( \dofMeasN \measNC + \res_\indexS^2 \right) $ \State $ \kalmM_{\indexS} \gets \MEstWeight_\indexS \stateCM_{\given{\indexS}{\indexS-1}} \trans{\outputM} \left( 1 + \MEstWeight_\indexS \outputM \stateCM_{\given{\indexS}{\indexS-1}} \trans{\outputM} \right)^{-1} $ \State $ \est{\stateV}_{\given{\indexS}{\indexS}} \gets \est{\stateV}_{\given{\indexS}{\indexS-1}} + \kalmM_{\indexS} \res_{\indexS} $ \State $ \stateCM_{\given{\indexS}{\indexS}} \gets \stateCM_{\given{\indexS}{\indexS-1}} - \kalmM_{\indexS} \outputM \stateCM_{\given{\indexS}{\indexS-1}} $ \end{algorithmic} \end{algorithm} %EM für Student's t als Mixture mit DOF fest/frei: [Murphy2013, p. 359ff] %Nurminen Thesis 2012: p. 10-13 Erklärung/Herleitung von IRLS = Gauss-Newton; p. 16 Gauss-Newton Covariance Estimation %============================================================================== \subsection{\acf{ORKF}} %============================================================================== Beim \ac{TKF} wurden die A-posteriori-Verteilung $\pdfBr{\given{\stateV_\indexS }{ \measV_{1:\indexS} }}$ und das Prozessrauschen $\pdfBr{\procNV_\indexS}$ vor dem Prädiktionsschritt sowie die prädiktive Verteilung $\pdfBr{\given{\stateV_\indexS }{ \measV_{1:\indexS-1} }}$ und das Messrauschen $\pdfBr{\measNV_\indexS}$ vor dem Update-Schritt jeweils durch eine gemeinsame Verteilung genähert. Im Gegensatz dazu nutzt der \ac{RME} für den Update-Schritt einen näherungsweisen Lösungsalgorithmus auf Basis eines M-Schätzers. Eine Alternative Herangehensweise für die Herleitung einer Näherungslösung für den Update-Schritt bietet die \acf{BV}. Entsprechende Filter wurden von Agamennoni et al. \cite{Agamennoni2011, Agamennoni2012}, Ardeshiri et al. \cite{Ardeshiri2015}, Piché et al. \cite{Piche2012} sowie Särkkä und Hartikainen \cite{Sarkka2009, Sarkka2013b} vorgestellt. Das in diesem Abschnitt vorgestellte Filter beruht auf \cite{Agamennoni2012}. Als Ausgangsbasis dient das Zustandsraummodell mit Student-t-Messrauschen \eqref{eq:tsHeavySS}. Die Filterherleitung basiert auf einer geschickten Zerlegung der Student-t-Verteilung. Dabei wird zunächst angenommen, dass bei bekanntem Zustandsvektor $\stateV_\indexS$ und bekannter Kovarianzmatrix $\measNCM_\indexS$ die Messungen normalverteilt sind: \begin{IEEEeqnarray}{c} \pdfBr{\given{\measV_\indexS }{ \stateV_\indexS, \measNCM_\indexS }} = \pdfNBr{\given{ \measV_\indexS }{ \outputM \stateV_\indexS + \feedM \inputV_\indexS, \measNCM_\indexS }}. \label{eq:tsORKFConditionallyNormalMeasurement} \end{IEEEeqnarray} Nun wird $\measNCM_\indexS$ aber nicht fest vorgegeben, sondern selbst für jeden Zeitschritt als Zufallsvariable betrachtet (vgl. $\measNCM_\indexS$ in \eqref{eq:tsORKFConditionallyNormalMeasurement} gegenüber dem festen Parameter $\measNCM$ in \eqref{eq:tsHeavySSMeasN}). Weiters erhält $\measNCM_\indexS$ eine A-priori-Verteilung \begin{IEEEeqnarray}{c} \pdfBr{\measNCM_\indexS} = \pdfInvWBr{\given{\measNCM_\indexS }{ \dofMeasN \measNCM, \dofMeasN }}. \label{eq:tsORKFWishartPrior} \end{IEEEeqnarray} Die inverse Wishart-Verteilung $\pdfInvW$ ist eine multivariate Verallgemeinerung der inversen Gamma-Verteilung. Sie kann u.a. als konjugierte A-priori-Verteilung für die Kovarianzmatrix einer Normalverteilung eingesetzt werden \cite[][102]{Bishop2006}. Der Parameter $\measNCM$ in \eqref{eq:tsORKFWishartPrior} gibt den harmonischen Mittelwert von $\measNCM_\indexS$ vor. Je kleiner die Anzahl an Freiheitsgraden $\dofMeasN$, desto stärker schwanken die Realisierungen von $\measNCM_\indexS$ um $\measNCM$. Die Relevanz von \eqref{eq:tsORKFConditionallyNormalMeasurement} und \eqref{eq:tsORKFWishartPrior} zeigt sich bei einer Marginalisierung über $\measNCM_\indexS$. Die resultierende Dichte \begin{IEEEeqnarray*}{c} \pdfBr{\given{\measV_\indexS }{ \stateV_\indexS }} = \int{ \pdfBr{\given{\measV_\indexS }{ \stateV_\indexS, \measNCM_\indexS }} \pdfBr{\measNCM_\indexS} \dif \measNCM_\indexS} = \pdftBr{\given{ \measV_\indexS }{ \outputM \stateV_\indexS + \feedM \inputV_\indexS, \measNCM, \dofMeasN }} \end{IEEEeqnarray*} entspricht einer Student-t-Verteilung und damit dem gewünschten Modell aus \eqref{ts:HeavySSMeasV} und \eqref{eq:tsHeavySSMeasN}. Um nun eine Lösung für das Filter-Update herzuleiten wird mit der Methode der \ac{BV} eine näherungsweise A-posteriori-Verteilung $\pdfqBr{\stateV_\indexS, \measNCM_\indexS}$ für jeden Zeitschritt $\indexS$ ermittelt. Als Maß für die Qualität der Näherung dient wie bei der \ac{BV} üblich die inverse Kullback-Leibler-Divergenz \begin{IEEEeqnarray}{c} \klBr{\pdfqBr{\stateV_\indexS, \measNCM_\indexS} }{ \pdfBr{ \given{ \stateV_\indexS, \measNCM_\indexS }{ \measV_{1:\indexS} }} }. \label{eq:tsORKFKL} \end{IEEEeqnarray} Die gewünschte Minimierung von \eqref{eq:tsORKFKL} durch Anpassung der Parameter von $\pdfqBr{\stateV_\indexS, \measNCM_\indexS}$ ist nur indirekt möglich, da $\pdfBr{ \given{ \stateV_\indexS, \measNCM_\indexS }{ \measV_{1:\indexS} }}$ ja gerade die gesuchte A-posteriori-Verteilung nach dem Update-Schritt ist. Hier hilft der Zusammenhang \cite[][463]{Bishop2006} \begin{subequations} \begin{IEEEeqnarray}{c} \log \pdfBr{\measV_{1:\indexS}} = \elboBr{\pdfqBr{\stateV_\indexS, \measNCM_\indexS}} + \klBr{\pdfqBr{\stateV_\indexS, \measNCM_\indexS} }{ \pdfBr{ \given{ \stateV_\indexS, \measNCM_\indexS }{ \measV_{1:\indexS} }} }, \IEEEyessubnumber* \label{eq:tsORKFLogLElboKl} \\ \elboBr{\pdfqBr{\stateV_\indexS, \measNCM_\indexS}} = \iint{ \pdfqBr{\stateV_\indexS, \measNCM_\indexS} \log \frac{ \pdfBr{ \stateV_\indexS, \measNCM_\indexS, \measV_{1:\indexS} } }{ \pdfqBr{\stateV_\indexS, \measNCM_\indexS} } \dif \stateV_\indexS \dif } \measNCM_\indexS . \label{eq:tsORKFElbo} \end{IEEEeqnarray} \end{subequations} Die Zerlegung in \eqref{eq:tsORKFLogLElboKl} ist exakt. Da $\pdfBr{\measV_{1:\indexS}}$ bezüglich der gesuchten Schätzgrößen $\stateV_\indexS$ und $\measNCM_\indexS$ konstant ist und $\klBr{\cdot}{\cdot} \geq 0$ gilt, ist die Maximierung von $\elboBr{\pdfqBr{\stateV_\indexS, \measNCM_\indexS}}$ gleichbedeutend mit der gewünschten (lokalen) Minimierung von $\klBr{\pdfqBr{\stateV_\indexS, \measNCM_\indexS} }{ \pdfBr{ \given{ \stateV_\indexS, \measNCM_\indexS }{ \measV_{1:\indexS} }} }$. In \eqref{eq:tsORKFElbo} wird jedoch im Gegensatz zu \eqref{eq:tsORKFKL} die gemeinsame Dichte $\pdfBr{ \stateV_\indexS, \measNCM_\indexS, \measV_{1:\indexS} }$ verwendet, die sich leichter auswerten lässt als $\pdfBr{ \given{ \stateV_\indexS, \measNCM_\indexS }{ \measV_{1:\indexS} }}$. Um die Optimierung von \eqref{eq:tsORKFElbo} zu ermöglichen, wird bei der \ac{BV} in der Regel zumindest eine teilweise Faktorisierung der genäherten Dichtefunktion $\pdfqBr{\stateV_\indexS, \measNCM_\indexS}$ angenommen. Dieses Vorgehen ist in der Literatur als (structured) Mean-Field-Approximation bekannt \cite[][52,60]{Beal2003}. Im Fall des \ac{ORKF} lautet die gewählte Zerlegung \begin{IEEEeqnarray*}{c} \pdfqBr{\stateV_\indexS, \measNCM_\indexS} = \pdfqBr{\stateV_\indexS} \pdfqBr{\measNCM_\indexS}. \end{IEEEeqnarray*} Der Vorteil einer solchen Zerlegung in unabhängige Faktoren ist, dass die \ac{BV} eine allgemeingültige Regel zur Bestimmung der Komponenten der genäherten Dichtefunktion liefert \cite[][464-466]{Bishop2006}. Für $\pdfqBr{\stateV_\indexS, \measNCM_\indexS}$ lautet die Lösung \begin{subequations} \label{eq:tsORKFVBR} \begin{IEEEeqnarray}{rCl} \log \pdfqBr{\stateV_\indexS} &=& \expecSubBr{\pdfqBr{\measNCM_\indexS}}{\log \pdfBr{ \stateV_\indexS, \measNCM_\indexS, \measV_{1:\indexS} }} + \text{const.}, \IEEEyessubnumber* \label{eq:tsORKFVBR1} \\ \log \pdfqBr{\measNCM_\indexS} &=& \expecSubBr{\pdfqBr{\stateV_\indexS }}{\log \pdfBr{ \stateV_\indexS, \measNCM_\indexS, \measV_{1:\indexS} }} + \text{const.}. \label{eq:tsORKFVBR2} \end{IEEEeqnarray} \end{subequations} Die beiden Gleichungen sind durch die Bildung des Erwartungswertes über die jeweils andere Dichtefunktion voneinander abhängig. Daher ist die praktische Herangehensweise zur Lösung das wechselweise iterative Berechnen von \eqref{eq:tsORKFVBR1} und \eqref{eq:tsORKFVBR2}, bis die Parameter der beiden genäherten Dichtefunktionen konvergieren. Die Eigenschaften der \ac{BV} garantieren die Konvergenz u einer (lokalen) Lösung, die \eqref{eq:tsORKFElbo} maximiert und somit \eqref{eq:tsORKFKL} minimiert. Die Herleitung der Parameter-Updates aus \eqref{eq:tsORKFVBR} erfordert zahlreiche Umformungen, es sei daher auf \cite{Agamennoni2012} verwiesen. Die genäherten A-posteriori-Verteilungen lauten \begin{IEEEeqnarray*}{rCl} \pdfqBr{\stateV_\indexS} &=& \pdfNBr{\given{ \stateV_\indexS }{ \est{\stateV}_{\given{ \indexS }{ \indexS }}, \stateCM_{\given{\indexS}{\indexS}} }}, \\ \pdfqBr{\measNCM_\indexS} &=& \pdfInvWBr{\given{\measNCM_\indexS }{ \left(\dofMeasN + 1\right) \tempORKFCM_\indexS, \left(\dofMeasN + 1\right) }} . \end{IEEEeqnarray*} % Der Parameter $\tempORKFCM_\indexS$ ist dabei das harmonische Mittel der geschätzten Kovarianz-Matrix des Messrauschens nach dem Update-Schritt $\indexS$. Den Update-Algorithmus für das \ac{ORKF} zeigt \algoref{alg:tsORKFUpdate}. \begin{algorithm}[H] \caption{Update-Schritt des \acs{ORKF} für die \acs{EOL}-Prüfung} \label{alg:tsORKFUpdate} \begin{algorithmic}[1] \setstretch{1.2} \State $ \est{\stateV}_{\given{\indexS}{\indexS}} \gets \est{\stateV}_{\given{\indexS}{\indexS-1}} $ \State $ \stateCM_{\given{\indexS}{\indexS}} \gets \stateCM_{\given{\indexS}{\indexS-1}} $ \While{not converged} \State $ \apr{\resV}_\indexS \gets \measV_\indexS - \outputM \est{\stateV}_{\given{\indexS}{\indexS}} - \feedM \inputV_\indexS $ \State $ \tempORKFCM_\indexS \gets \frac{\dofMeasN}{\dofMeasN+1} \measNCM + \frac{1}{\dofMeasN+1} ( \apr{\resV}_\indexS \apr{\resV}_\indexS^T + \outputM \stateCM_{\given{\indexS}{\indexS}} \outputM^T) $ \State $ \kalmM_\indexS \gets \stateCM_{\given{\indexS}{\indexS-1}} \outputM^T ( \outputM \stateCM_{\given{\indexS}{\indexS-1}} \outputM^T + \tempORKFCM_\indexS)^{-1} $ \State $ \est{\stateV}_{\given{\indexS}{\indexS}} \gets \est{\stateV}_{\given{\indexS}{\indexS-1}} + \kalmM_\indexS ( \measV_\indexS - \outputM \est{\stateV}_{\given{\indexS}{\indexS-1}} - \feedM \inputV_\indexS ) $ \State $ \stateCM_{\given{\indexS}{\indexS}} \gets \kalmM_\indexS \tempORKFCM_\indexS \kalmM_\indexS^T + (\eye - \kalmM_\indexS \outputM ) \stateCM_{\given{\indexS}{\indexS-1}} (\eye - \kalmM_\indexS \outputM )^T $ \EndWhile \end{algorithmic} \end{algorithm} %VB grundsätzlich erklärt im Rahmen von Filtern in [Nurminen2017, p. 18-21], [Beal2003], [Ostwald2014] %============================================================================== \subsection{Reste für analytische Filter} %============================================================================== [Meinhold1989] "a robust KF model would be one for which the posterior distribution of the state of nature would return to ist prior as the observation departs significantly from ist predicted value" Interessanter Vergleich heaviness prior und likelihood p. 480 links Vorteile Student's t measurement noise: [Piche2012, p. 1] Beispiele für andere Verteilungen als Student's t laut [Roth2017c]: Sornette2001 und Gordon2003 Optimization Viewpoint for robust / non-Gaussian filtering: siehe Aravkin / Burke / Pillonetto [Durbin2001 Kap. 9.4] [Harvey2013] [Fahrmeir1999]: State Space; Student's t; EM-type Algorithm; Smoother; Hyperparameter estimation (CV, ML, EM) %============================================================================== \subsection{Vergleich robuster Filteralgorithmen auf simulierten Datensätzen} \label{sec:tsFilterCompareSim} %============================================================================== Der Vergleich der vorgestellten Filteralgorithmen bezüglich ihrer Eignung für die Zustandsschätzung erfolgt auf simulierten Datensätzen. Der Grund für die Nutzung einer Simulation ist, dass hierfür die Modellparameter des Zustandsraummodells genau bekannt sind. Bei einem Datensatz aus der Anwendung ist hingegen eine Schätzung der Modellparameter erforderlich (siehe \secref{sec:tsParamEst}), die eine von anderen Verfahren unabhängige Bewertung der Filteralgorithmen verhindert. Darüber hinaus sind bei einem simulierten Datensatz die tatsächlichen Werte des Zustands $\stateV_\indexS$ bekannt, die in der Anwendung nicht als Referenz zur Verfügung stehen. %----------------------------------------------------------------------------- \subsubsection{Interpretation der Zeitverläufe der Schätzgrößen} %----------------------------------------------------------------------------- Einen Eindruck für das Verhalten der Filter liefert der Vergleich auf einem kurzen Zeitabschnitt, bei dem die Resultate der einzelnen Update-Schritte erkennbar sind. Dazu wurde ein Random-Walk mit Student-t-Messrauschen nach \eqref{eq:tsHeavySS} mit $\numS = \num{50}$ Samples simuliert, die Parameter zeigt \tabref{tab:tsSimStudentShort}. % EC1112 03_01 \begin{table} \caption[FIXME]{Parameter für die Simulation eines kurzen Random-Walks mit Student-t-Messrauschen} \label{tab:tsSimStudentShort} \centering \footnotesize \begin{tabular}{ll} \toprule Parameter & Wert \\ \midrule $\procNC$ & \num{0.1} \\ $\measNC$ & \num{0.1} \\ $\dofMeasN$ & \num{5} \\ $\state_1$ & \num{0} \\ $\Delta t_{\indexS}$ & \num{1} \\ $\stateG$ & \num{1} \\ $\inputG$ & \num{0} \\ $\outputG$ & \num{1} \\ $\feedG$ & \num{0} \\ \bottomrule \end{tabular} \end{table} Der simulierten Random-Walk und die Schätzergebnisse der Filteralgorithmen sind in \figref{fig:tsRobustFilterTS} dargestellt. In die Messwerte wurde dabei ein Ausreißer künstlich eingefügt, um die Reaktion der Filter darauf zu testen. Alle vier Filter starten mit einer A-priori-Verteilung mit den Parametern $\state_{\given{1}{0}} = \num{0}$ und $\stateC_{\given{1}{0}} = \num{1e-9}$. Der initiale Zustand ist den Filter damit praktisch perfekt bekannt. Das \ac{GKF} baut auf einer Normalverteilungsannahme für $\measN_\indexS$ auf, das Messrauschen in der Simulation weist aber eine Student-t-Verteilung auf. Entsprechend wurde beim \ac{GKF} gemäß der Eigenschaften der Student-t-Verteilung \cite[][578\psq]{Gelman2014} die Varianz des Messrauschens auf \begin{IEEEeqnarray*}{c} \varBr{\measN_\indexS} = \frac{\dofMeasN}{\dofMeasN-2} \measNC \end{IEEEeqnarray*} eingestellt. Alle vier Filter sind grundsätzlich in der Lage, der Entwicklung von $\state_\indexS$ zu folgen. Beim \ac{GKF} ist innerhalb von fünf Update-Schritten in etwa der stationäre Wert des Gain $\kalm_\indexS$ erreicht, die Varianz $\stateC_{\given{ \indexS }{ \indexS }}$ der Zustandsschätzung ist nicht mehr von ihrem stationären Wert unterscheidbar. Diesen Wert behält das \ac{GKF} bei, bis das Validation-Gate durch den Ausreißer bei $t = \num{15}$ aktiviert wird. Die anderen drei Filteralgorithmen weisen keine stationär konstante Kovarianz auf, da deren Update-Schritt beim \ac{TKF}, \ac{RME} und \ac{ORKF} immer abhängig von der jeweiligen Messung $\meas_\indexS$ ist. Die charakteristischen Eigenschaften der Filter treten bei der Behandlung des eingefügten Ausreißers hervor. Der geschätzte Zustand $\est{\state}_{\given{ 14 }{ 14 }}$ nach dem letzten Update vor dem Ausreißer ist bei allen Filtern praktisch gleich. Das \ac{GKF} schließt den Ausreißer bei $t = \num{15}$ aus dem Update-Schritt aus, der entsprechende Gain ist $\kalm_{15} = \num{0}$. Entsprechend bleibt die Zustandsschätzung $\est{\state}_{\given{ 15 }{ 15 }} = \est{\state}_{\given{ 14 }{ 14 }}$ unverändert und die Varianz $\stateC_{\given{ 15 }{ 15 }} > \stateC_{\given{ 14 }{ 14 }} $ steigt gemäß dem Prozessrauschen an. Alle nachfolgenden $\meas_\indexS, \indexS > \num{15}$ liegen wieder innerhalb des Gate des \ac{GKF}, sodass sich Gain und Schätzvarianz schnell wieder auf die stationären Werte zubewegen. Das Verhalten des \ac{TKF} unterscheidet sich deutlich vom \ac{GKF}. Das Update des Zustands erfolgt wie beim herkömmlichen Kalman-Filter proportional zum Residuum $\res_\indexS$, daher ist die Zustandsschätzung $\est{\state}_{\given{ 15 }{ 15 }}$ gegenüber dem wahren Zustand stark verfälscht. Gleichzeitig erhöht das \ac{TKF} die Varianz $\stateC_{\given{ 15 }{ 15 }}$ stark, wodurch der Gain im nachfolgenden Updateschritt $\indexS = 16$ höher ist als bei allen anderen Filtern. Somit wird der Einfluss des Ausreißers auf die Zustandsschätzung beim \ac{TKF} erst einen Update-Schritt nach dem Ausreißer wieder weitgehend korrigiert. Der \ac{RME}- und \ac{ORKF}-Algorithmus weisen ein zueinander sehr ähnliches Update-Verhalten bezüglich des Ausreißers auf. Der Gain $\kalm_{15}$ ist gegenüber dem \ac{TKF} deutlich verringert und verhindert somit ähnlich dem \ac{GKF} einen starken Einfluss des Ausreißers auf die Zustandsschätzung. Für extreme Ausreißer gilt im Grenzfall $\abs{\meas_\indexS} \rightarrow \infty$ beim \ac{RME} und \ac{ORKF} $\kalm_\indexS \rightarrow 0$, sodass das Verhalten in das des \ac{GKF} übergeht. \ac{RME} und \ac{ORKF} weisen im Falle von Ausreißern somit ein Verhalten zwischen dem des \ac{GKF} und des \ac{TKF} auf. % ------------------------------------------------------------------- \begin{figure} \centering \footnotesize \setlength\figurewidth{.9\columnwidth} \setlength\figureheight{.35\columnwidth} \subcaptionbox{Messungen $\meas_\indexS$ \circleWithParens{myMediumBlue} und wahrer Zustand $\state_\indexS$ \lineWithParens{myMediumBlue} sowie die Zustandsschätzungen $\est{\state}_{\given{ \indexS }{ \indexS }}$ der Filter. In die Messwerte wurde bei $t = \num{15}$ ein Ausreißer \circleWithParens{red} eingefügt.\label{fig:tsRobustFilterTSState}}[\columnwidth]{% \centering \footnotesize \renewcommand\figureXLabel{$t$} \renewcommand\figureYLabel{$\state_\indexS$, $\meas_\indexS$, $\est{\state}_{\given{\indexS}{\indexS}}$} \mytikzexton \input{"graphics/P1112 03 01/301 thesis t common state.tikz"} \mytikzextoff }% % --------------- comment magic (tm) \\% \vspace{4mm}% % --------------- comment magic (tm) \subcaptionbox{Gain in den Update-Schritten\label{fig:tsRobustFilterTSGain}}[\columnwidth]{% \centering \footnotesize \renewcommand\figureXLabel{$t$} \renewcommand\figureYLabel{$\kalm_\indexS$} \mytikzexton \input{"graphics/P1112 03 01/303 thesis t common K.tikz"} \mytikzextoff }% % --------------- comment magic (tm) \\% \vspace{4mm}% % --------------- comment magic (tm) \subcaptionbox{Standardabweichung nach den Update-Schritten\label{fig:tsRobustFilterTSCov}}[\columnwidth]{% \centering \footnotesize \renewcommand\figureXLabel{$t$} \renewcommand\figureYLabel{$\sqrt{\stateC_{\given{\indexS}{\indexS}}}$} \mytikzexton \input{"graphics/P1112 03 01/302 thesis t common cov.tikz"} \mytikzextoff }% \caption[FIXME]{Schätzergebnisse der Filteralgorithmen \ac{GKF} \lineWithParens{black}, \ac{TKF} \lineWithParens{myLineTwo}, \ac{RME} \lineWithParens{myLineThree} und \ac{ORKF} \lineWithParens{myLineFour} auf einem simulierten Random-Walk mit Student-t-Messrauschen.} \label{fig:tsRobustFilterTS} \end{figure} %----------------------------------------------------------------------------- \subsubsection{Bewertung der Schätzgrößen anhand eines langen Random-Walks} %----------------------------------------------------------------------------- Für eine quantitative Bewertung der Schätzergebnisse ist der Random-Walk zu kurz. Zusätzlich erfolgt daher nun die Auswertung anhand eines langen Random-Walks mit $\numS = \num{100000}$ Samples. Um einen Einfluss der Filterinitialisierung zu vermeiden finden die ersten \num{100} Samples keine Berücksichtigung in den Ergebnissen. Die Modellparameter zeigt \tabref{tab:tsSimStudentLong}, die Ergebnisse sind in \figref{fig:tsRobustFilterHist} in Form von Histogrammen für die einzelnen Filter dargestellt. % EC1112 04_01 \begin{table} \caption[FIXME]{Parameter für die Simulation eines langen Random-Walks mit Student-t-Messrauschen} \label{tab:tsSimStudentLong} \centering \footnotesize \begin{tabular}{ll} \toprule Parameter & Wert \\ \midrule $\procNC$ & \num{0.1} \\ $\measNC$ & \num{1.0} \\ $\dofMeasN$ & \num{5} \\ $\state_1$ & \num{0} \\ $\Delta t_{\indexS}$ & \num{0.1} \\ $\stateG$ & \num{1} \\ $\inputG$ & \num{0} \\ $\outputG$ & \num{1} \\ $\feedG$ & \num{0} \\ \bottomrule \end{tabular} \end{table} Alle vier Filter modellieren den geschätzten Zustand nach dem Update-Schritt intern als Normalverteilung $\pdfBr{\given{\state_\indexS }{ \meas_\indexS }} = \pdfNBrS{ \given{ \state_\indexS }{ \est{\state}_{\given{ \indexS }{ \indexS }}, \stateC_{\given{ \indexS }{ \indexS }} } }$. %Daher bietet es sich für einen Vergleich der Filter an, die mit einem Cholesky-Faktor $\LTM_\indexS$ von $\stateCM_{\given{ \indexS }{ \indexS }}$ skalierte Abweichung der Zustandsschätzung vom simulierten Zustand zu betrachten (vgl. \secref{sec:basicsMathStoch}): %\begin{IEEEeqnarray*}{rCl} %\stateCM_{\given{ \indexS }{ \indexS }} &=& \LTM \trans{\LTM}, \\ %\res_{\state, \indexS} &=& \inv{\LTM} \left(\stateV_\indexS - \est{\stateV}_{\given{ \indexS }{ \indexS }}\right) \sim \pdfNBr{\nullV, \eye}. %\end{IEEEeqnarray*} Daher bietet es sich für eine Untersuchung der Konsistenz der gelieferten Schätzwerte $\est{\state}_{\given{ \indexS }{ \indexS }}$ und $\stateC_{\given{ \indexS }{ \indexS }}$ an, die mit der inversen Standardabweichung $\nicefrac{1}{ \sqrt{\stateC_{\given{ \indexS }{ \indexS }} }}$ skalierte Abweichung der Zustandsschätzungen vom simulierten Zustand zu betrachten (vgl. \secref{sec:basicsMathStoch}): \begin{IEEEeqnarray}{c} \frac{ \state_\indexS - \est{\state}_{\given{ \indexS }{ \indexS }} }{\sqrt{\stateC_{\given{ \indexS }{ \indexS }} } } \sim \pdfNBr{0, 1}. \label{eq:tsIdealStateResidual} \end{IEEEeqnarray} Gleichung \eqref{eq:tsIdealStateResidual} stellt dabei den Idealzustand dar, der nur erreicht wird wenn der jeweilige Filteralgorithmus exakt arbeitet. Aufgrund der bei den einzelnen Näherungen erreichen die tatsächlichen Schätzgrößen diese Verteilung nicht genau. \Figref{fig:tsRobustFilterHist} zeigt links jeweils ein Histogramm der mit \eqref{eq:tsIdealStateResidual} standardisierten Abweichung und als Referenz die Standardnormalverteilung $\pdfNBr{0, 1}$. In der rechten Spalte ist die zugehörige Standardabweichung $\sqrt{\stateC_{\given{ \indexS }{ \indexS }}}$ wie sie vom jeweiligen Filter geliefert wird ebenfalls als Histogramm dargestellt. Als Referenz dient hier die tatsächliche Stichproben-Standardabweichung \begin{IEEEeqnarray*}{c} \est{\stdSym}_\state = \sqrt{ \frac{1}{\numS-1} \sum_{\indexS=1}^{\numS}{\left(\state_\indexS - \est{\state}_{\given{ \indexS }{ \indexS }} \right)^2 } }. \end{IEEEeqnarray*} Sie wurde aus den Realisierungen von $\state_\indexS - \est{\state}_{\given{ \indexS }{ \indexS }} $ ermittelt. \texttodo{In \figref{fig:tsRobustFilterHist} die Referenzkurven Schwarz und Blau färben} Bezüglich der standardisierten Abweichung bei der Zustandsschätzung %$\left(\state_\indexS - \est{\state}_{\given{ \indexS }{ \indexS }}\right) / \sqrt{\stateC_{\given{ \indexS }{ \indexS }} } $ zeigt das \ac{GKF} eine nahezu perfekte Übereinstimmung zwischen der theoretischen und der tatsächlichen Verteilung. \ac{TKF}, \ac{RME} und \ac{ORKF} unterschätzen jeweils die Varianz ihrer eigenen Zustandsschätzung. Beim \ac{TKF} ist dieser Effekt gleichzeitig mit einer großen Streuung von $\sqrt{\stateC_{\given{ \indexS }{ \indexS }}}$ verbunden, die aufgrund des Verhaltens in \figref{fig:tsRobustFilterTS} bereits zu erwarten war. Große Residuen $\res_\indexS$ führen hier zu einer starken Abweichung des Zustands $\est{\state}_{\given{ \indexS }{ \indexS }}$ in Verbindung mit einer gleichzeitigen Erhöhung von $\stateC_{\given{ \indexS}{ \indexS}}$. Die im Mittel etwas zu kleinen Werte von $\sqrt{\stateC_{\given{ \indexS}{ \indexS}}}$ beim \ac{ORKF} sind typisch für Schätzverfahren auf Basis der Kullback-Leibler-Divergenz $\klBr{\pdfq}{\pdf}$ \cite[][466-468]{Bishop2006}. Der \ac{RME} zeigt ein vergleichbares Verhalten. % ------------------------------------------------------------------- \begin{figure} \centering \footnotesize \setlength\figurewidth{.40\columnwidth} \setlength\figureheight{.23\columnwidth} \renewcommand\figureXLabel{$\left(\state_\indexS - \est{\state}_{\given{ \indexS }{ \indexS }}\right) / \sqrt{\stateC_{\given{ \indexS }{ \indexS }}}$ } \subcaptionbox{\ac{GKF}}[.99\columnwidth]{% \centering \footnotesize \mytikzexton %\renewcommand\figureYLabel{$\pdfBr{\state_\indexS - \est{\state}_{\given{ \indexS }{ \indexS }} / \sqrt{\stateC_{\given{ \indexS }{ \indexS }}} }$ } \renewcommand\figureYLabel{$\pdf$} \input{"graphics/P1112 04 01/401 Student Walk GKF x Histogram.tikz"} \hspace{0mm} % manual fine-tuning of horizontal alignment %\renewcommand\figureYLabel{$\pdfBr{\sqrt{\stateC_{\given{ \indexS }{ \indexS }}} }$ } \renewcommand\figureYLabel{$\pdf$} \input{"graphics/P1112 04 01/406 Student Walk P_diag_k Histogram GKF.tikz"} \mytikzextoff }% % --------------- comment magic (tm) \\% \vspace{3mm}% % --------------- comment magic (tm) \subcaptionbox{\ac{TKF}}[.99\columnwidth]{% \centering \footnotesize \mytikzexton %\renewcommand\figureYLabel{$\pdfBr{\state_\indexS - \est{\state}_{\given{ \indexS }{ \indexS }} / \sqrt{\stateC_{\given{ \indexS }{ \indexS }}} }$ } \renewcommand\figureYLabel{$\pdf$} \input{"graphics/P1112 04 01/402 Student Walk TKF x Histogram.tikz"} \hspace{3.5mm} % manual fine-tuning of horizontal alignment %\renewcommand\figureYLabel{$\pdfBr{\sqrt{\stateC_{\given{ \indexS }{ \indexS }}} }$ } \renewcommand\figureYLabel{$\pdf$} \input{"graphics/P1112 04 01/407 Student Walk P_diag_k Histogram TKF.tikz"} \mytikzextoff }% % --------------- comment magic (tm) \\% \vspace{3mm}% % --------------- comment magic (tm) \subcaptionbox{\ac{RME}}[.99\columnwidth]{% \centering \footnotesize \mytikzexton %\renewcommand\figureYLabel{$\pdfBr{\state_\indexS - \est{\state}_{\given{ \indexS }{ \indexS }} / \sqrt{\stateC_{\given{ \indexS }{ \indexS }}} }$ } \renewcommand\figureYLabel{$\pdf$} \input{"graphics/P1112 04 01/403 Student Walk RKF x Histogram.tikz"} \hspace{2mm} % manual fine-tuning of horizontal alignment %\renewcommand\figureYLabel{$\pdfBr{\sqrt{\stateC_{\given{ \indexS }{ \indexS }}} }$ } \renewcommand\figureYLabel{$\pdf$} \input{"graphics/P1112 04 01/408 Student Walk P_diag_k Histogram RKF.tikz"} \mytikzextoff }% % --------------- comment magic (tm) \\% % --------------- comment magic (tm) \subcaptionbox{\ac{ORKF}}[.99\columnwidth]{% \centering \footnotesize \mytikzexton %\renewcommand\figureYLabel{$\pdfBr{\state_\indexS - \est{\state}_{\given{ \indexS }{ \indexS }} / \sqrt{\stateC_{\given{ \indexS }{ \indexS }}} }$ } \renewcommand\figureYLabel{$\pdf$} \input{"graphics/P1112 04 01/404 Student Walk ORKF x Histogram.tikz"} \hspace{2mm} % manual fine-tuning of horizontal alignment \renewcommand\figureXLabel{$\sqrt{\stateC_{\given{ \indexS }{ \indexS }}}$ } %\renewcommand\figureYLabel{$\pdfBr{\sqrt{\stateC_{\given{ \indexS }{ \indexS }}} }$ } \renewcommand\figureYLabel{$\pdf$} \input{"graphics/P1112 04 01/409 Student Walk P_diag_k Histogram ORKF.tikz"} \mytikzextoff }% \caption[FIXME]{Histogramme der Zustandsschätzung auf einem simulierten Random-Walk mit Student-t-Messrauschen. Den standardisierten Zustandsschätzungen ist eine Normalverteilung $\pdfNBr{0,1}$ \lineWithParens{myLineOne} überlagert. Den Histogrammen von $\sqrt{\stateC_{\given{ \indexS }{ \indexS }}}$ ist der empirische Wert $\est{\stdSym}_\state$ überlagert \lineWithParens{myLineTwo}. } \label{fig:tsRobustFilterHist} \end{figure} Einen abschließenden zahlenmäßigen Vergleich der absoluten Genauigkeit der Verfahren bei der robusten Zustandsschätzung liefert der RMS-Wert \begin{IEEEeqnarray*}{c} \RMSBr{ \state_\indexS - \est{\state}_{\given{ \indexS }{ \indexS }} } = \sqrt{ \frac{1}{\numS} \sum_{\indexS=1}^{\numS}{ \left( \state_\indexS - \est{\state}_{\given{ \indexS }{ \indexS }} \right)^2 }}. \end{IEEEeqnarray*} Die Werte für den langen Random-Walk zeigt \tabref{tab:tsSimStudentLongRMS}. Das \ac{GKF} erreicht dabei den besten (geringsten) Wert, das \ac{TKF} den schlechtesten (größten). \ac{RME} und \ac{ORKF} liegen dazwischen, insgesamt sind die absoluten Unterschiede in der Genauigkeit der Zustandsschätzung zwischen den Verfahren klein. \begin{table} \caption[FIXME]{RMS der Zustandsschätzung robuster Filteralgorithmen} \label{tab:tsSimStudentLongRMS} \centering \footnotesize \begin{tabular}{ll} \toprule Filter & $\RMSBr{ \state_\indexS - \est{\state}_{\given{ \indexS }{ \indexS }} }$ \\ \midrule \ac{GKF} & \num{0.3472} \\ \ac{TKF} & \num{0.3636} \\ \ac{RME} & \num{0.3388} \\ \ac{ORKF} & \num{0.3380} \\ \bottomrule \end{tabular} \end{table} Die robuste Zustandsschätzung ist letztlich ein Hilfsmittel zur Berechnung adaptiver Prüfgrenzen, die Zustandsschätzung stellt also ein Zwischenergebnis dar. Die Vorgehensweise zur Ermittlung der Prüfgrenzen aus den Schätzergebnissen der Filter und die spezifischen Einflüsse der Filtereigenschaften darauf beschreibt der nächste Abschnitt. %Anschauliche graphische Darstellung was bei zunehmendem Residuum Messung - Prädiktion beim nicht-robusten und beim robusten Filter (verschiedene DOF) mit der Posterior der Zustandsschätzung passiert -> vgl. Diskussion [Meinhold1989, p. 481, p. 482, Fig. 1, Fig. 2] %siehe Plot [Roth2017b, p. 29 und weitere] %Kalman-Filter ist der beste lineare Schätzer bei gegebenem ($\mu$, Var) des Rauschens; Sichtweise: Outlier verändern die Varianz des Rauschens stark %[Aravkin2016, p. 5] \glqq{}In this case, the smoother is aware of the true variance of the signal; nonetheless, the reconstruction is still not satisfactory, since it cannot track the true output profile given the high measurement variance; the best linear estimate essentially averages the signal. Manipulating noise statistics is clearly not enough; to improve the estimator performance, we must change our model for the underlying distribution of the errors et.\grqq %Vollständige Filterbewertung nur bei Simulationen möglich [BarShalom2001, p. 234] %[BarShalom2001, p. 58] The chi-square distribution is often used to check state estimators for “consistency”- that is, whether their actual errors are consistent with the variances calculated by the estimator %[BarShalom2001, p. 233-234] %Filter bewerten durch Auswertung von Zustandsschätzung und Kovarianzschätzung %Chi-Squared-Verteilung der (mit der Kovarianz standardisierten) Zustands-Schätzfehler %Chi-Squared-Verteilung der (mit der Kovarianz standardisierten) Innovationen %Testen der Zustands-Schätzfehler und Innovations-Schätzfehler [BarShalom2001, p. 234] %Weitere Möglichkeiten für die Filter-Bewertung: mean/std/skewness des Fehlers, root-mean-square-error [Nurminen2015, sec. V] %Mehra1995: "It can be shown [Mehra1972], [Anderson1979] that for an optimal linear filter, the innovation sequence $\nu_k$, given by ... is a zero mean Gaussian white noise sequence with covariance $\Sigma_k$ ... . %This innovations property is useful in testing for optimality of a linear Kalman filter, to detect changes in the process model, and to build adaptive and robust Kalman filters." %Schöne Plots zu Verfahrens-Vergleichen [Nurminen2016] %Sensitivitätsuntersuchung Kalman-Filter auf Modellfehler [BarShalom2001, Ch. 5.6] %Allgemein Eigenschaften der untersuchten genäherten Filter zeigen => Check posterior (predictive) distribution, Simulation mit vielen Walks etc., vgl [Gelman2014, Ch. 6 Model Checking] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Berechnung von Prüfgrenzen mit robusten Filteralgorithmen} \label{sec:tsTestLimits} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% In \secref{sec:tsRobustFilter} wurden insgesamt vier stochastische Filteralgorithmen beschrieben, die eine robuste Zustandsschätzung für Zustandsraummodelle ermöglichen. Der Aspekt der Robustheit hat sich dabei auf die Fähigkeit bezogen, zuverlässige Zustandsschätzungen zu liefern wenn in den Messdaten einzelne unbekannten Ausreißern vorhanden sind. Die Robustheit wurde beim \ac{GKF} auf Basis des linear-gaußschen Zustandsraummodells durch ein schaltendes Verhalten des Update-Schritts erzielt. Das \ac{TKF}, \ac{RME} und \ac{ORKF} beruhen hingegen auf dem Zustandsraummodell mit Student-t-Messrauschen und analytischen Approximationen für die Herleitung der jeweiligen Update-Gleichungen. Für die Berechnung von Prüfgrenzen auf Basis der Zustandsschätzungen der Filter wird nun wie bereits bei der Untersuchung der statischen Klassifikationsverfahren in \chref{ch:static} auf das Prinzip der Auffälligkeitserkennung zurückgegriffen. Eine Messung $\measV_\indexS$ (und damit der zugehörige Prüflauf) gilt als auffällig, wenn die Messwerte $\measV_\indexS$ der Prüfmerkmale stark von der Prädiktion abweichen, die wiederum anhand der vorangegangenen Messungen $\measV_{1:\indexS-1}$ erstellt wurde. Wie groß eine Abweichung der Merkmalswerte absolut gesehen sein muss um als auffällig zu gelten definiert sich dabei erstens anhand der Unsicherheit der Prädiktion aus dem Filter-Algorithmus. Je weniger \glqq{}Wissen\grqq{} über die zu erwartenden Merkmalswerte $\measV_\indexS$ in Form einer eng umgrenzten Prädiktion $\pdfBr{ \given{ \measV_\indexS }{ \measV_{1:\indexS-1} } }$ vorhanden ist, umso unschärfer muss auch die Prüfgrenze ausfallen. Zweitens soll der Nutzer des Prüfsystems vorgeben können, ab welcher Abweichung relativ zur Unsicherheit der Prädiktion eine Realisierung von $\measV_\indexS$ nicht mehr akzeptiert wird. Die Berücksichtigung beider Einflussgrößen -- die Unsicherheit der Prädiktion und die nutzerdefinierte relative Mindestauffälligkeit -- geht auf wesentliche Kritikpunkte am Klassifikationsverfahren in \chref{ch:static} ein. %============================================================================== \subsection{Prüfgrenzen für das \ac{GKF}} %============================================================================== \texttodo{Rel. einfaches Bild mit Zeitachse nach rechts, Werte nach oben; Dichte x k k -> Dichte x k+1 k -> Dichte y k+1 k, diese in zwei Varianten normal und heavy} Tatsächlich ist genau dieses Vorgehen zur Erkennung auffälliger Messwerte bereits vom \ac{GKF} bekannt. Gleichung \eqref{eq:tsGKFPredictMeasurement} beschreibt die Prädiktion der Messung $\measV_\indexS$ auf Basis der vorangegangenen Messungen als %\begin{subequations} %\label{eq:tsGKFPredictMeasurementAgain} \begin{IEEEeqnarray}{rCl} \pdfBr{ \given{\measV_{\indexS}} {\measV_{1:\indexS-1}} } &=& \int{ \pdfBr{ \given{\measV_\indexS }{ \stateV_{\indexS}} } \: \pdfBr{ \given{\stateV_{\indexS} }{ \measV_{1:\indexS-1}} } \dif \stateV_\indexS} \IEEEnonumber* \\ &=& \int{ \pdfNBr{ \given{\measV_\indexS }{ \outputM \stateV_{\indexS} + \feedM \inputV_{\indexS}, \measNCM } } \: \pdfNBr{ \given{\stateV_{\indexS} }{ \est{\stateV}_{\given{\indexS}{\indexS-1}}, \stateCM_{\given{\indexS}{\indexS-1}} }} \dif \stateV_{\indexS} } \\ &=& \pdfNBr{ \given{\measV_{\indexS}} {\est{\measV}_{\given{\indexS} {\indexS-1} } }, \measCM_{\indexS} }. \IEEEyesnumber \label{eq:tsGKFPredictMeasurementAgainResult} \end{IEEEeqnarray} %\end{subequations} Die Parameter der Dichtefunktionen liefert der Prädiktions-Schritt des Filters. Hierbei ist zu beachten, dass \eqref{eq:tsGKFPredictMeasurementAgainResult} den \ac{iO}-Fall beschreibt, also jene Messungen die dem linear-gaußschen Zustandsraummodell \eqref{eq:tsLGSS} entsprechen. Ein Ausreißer im Sinne einer auffälligen Abweichung vom linear-gaußschen Zustandsraummodell wird beim \ac{GKF} über die inverse Verteilungsfunktion $\invCdf$ der Chi-Quadrat-Verteilung gemäß \eqref{eq:tsGKFCDFInv} erkannt. Mit diesem Zusammenhang lässt sich für das \ac{GKF} analog zum Validation-Gate ein Schwellwert $\SNPE_\test$ für $\SNPE(\resV_\indexS)$ finden, dessen Überschreitung eine Verletzung der Prüfgrenze anzeigt: \label{eq:tsGKFBoundary} \begin{IEEEeqnarray*}{c} \SNPE_\test = \invCdfBr{\Prob_\test}. \end{IEEEeqnarray*} $\Prob_\test$ definiert dabei den genannten relativen Grad der Abweichung der Messung von der Prädiktion, mit dem der Nutzer die Prüfschärfe im Verhältnis zur Unschärfe der Prädiktion festlegt. Als Resultat ergibt sich wie beim Validation-Gate ein Hyper-Ellipsoid um $\est{\measV}_{\given{\indexS} {\indexS-1} }$, das nun die Grenze zwischen einem \ac{iO}- (innen) und einem \ac{niO}-Prüfergebnis bildet: \begin{IEEEeqnarray*}{c} \text{Prüfergebnis für } \measV_\indexS = \begin{cases} \text{\ac{iO}} & \text{wenn } \SNPE(\resV_\indexS) \leq \SNPE_\test, \\ \text{\ac{niO}} & \text{sonst}. \end{cases} \end{IEEEeqnarray*} Für die praktische Anwendung wird nachfolgend immer $\Prob_\test = \Prob_\gate$ gesetzt. Im skalaren Fall $\nMeasV = 1$ ist der Zusammenhang sehr anschaulich, da die Festlegung der Prüfgrenze anstatt über den Umweg über $\SNPE(\resV_\indexS)$ direkt mittels der inversen Verteilungsfunktion von \eqref{eq:tsGKFPredictMeasurementAgainResult} erfolgen kann. Ein Beispiel zeigt \figref{fig:tsGKFBoundaryIllustration}. Der um $\res_\indexS = 0$ liegende \ac{iO}-Bereich umfasst eine Gesamtwahrscheinlichkeit $\Prob_\test$. Der \ac{iO}-Bereich ist das kleinstmögliche solche Intervall für diese Dichtefunktion bei gegebenem $\Prob_\test$. Ein solches Intervall einer Dichtefunktion ist in der Literatur auch als \ac{HDI} bekannt, da alle Punkte innerhalb des Intervalls einen gleich großen oder größeren Wert der Wahrscheinlichkeitsdichtefunktion aufweisen wie dessen Rand. Bei multivariaten Verteilungen ist die Verallgemeinerung die \ac{HDR} \cite[][32-34]{Gelman2014} \cite[][87]{Kruschke2015}. % ------------------------------------------------------------------- \begin{figure} \centering \footnotesize \setlength\figurewidth{.4\columnwidth} \setlength\figureheight{.3\columnwidth} \renewcommand\figureXLabel{$\res_\indexS / \sqrt{\measC_\indexS} $} \renewcommand\figureYLabel{$\pdfBr{\given{ \res_\indexS / \sqrt{\measC_\indexS} }{ \meas_{1:\indexS-1} }} $ } \subcaptionbox{Wahrscheinlichkeitsdichtefunktion \lineWithParens{myLineOne} des standardisierten Residuums}[\columnwidth]{% \centering \footnotesize \mytikzexton \input{"graphics/P9031 01/001 PDF.tikz"} \mytikzextoff }% % --------------- comment magic (tm) \\% \vspace{0mm}% % --------------- comment magic (tm) \subcaptionbox{Verteilungsfunktion \lineWithParens{myLineOne} des standardisierten Residuums}[\columnwidth]{% % \centering \footnotesize \setlength\figurewidth{.528\columnwidth} % width2 = width1 * 1.319 \setlength\figureheight{.3\columnwidth} \renewcommand\figureXLabel{$\res_\indexS / \sqrt{\measC_\indexS} $} \renewcommand\figureYLabel{$\cdfBr{\given{ \res_\indexS / \sqrt{\measC_\indexS} }{ \meas_{1:\indexS-1} }} $ } \mytikzexton \hspace{-9mm} \input{"graphics/P9031 01/002 CDF.tikz"} \mytikzextoff }% \caption[FIXME]{Darstellung der Verteilung des standardisierten Residuums und der daraus abgeleiteten Prüfgrenze \lineWithParens{black} mit \ac{iO}-Bereich \squareWithParens{myGray75} im eindimensionalen Fall für $\Prob_\test = \num{0.955}$.} \label{fig:tsGKFBoundaryIllustration} \end{figure} %============================================================================== \subsection{Prüfgrenzen für \ac{TKF}, \ac{RME} und \ac{ORKF}} %============================================================================== Da die Filteralgorithmen \ac{TKF}, \ac{RME} und \ac{ORKF} auf dem Random-Walk-Modell mit Student-t-Messrauschen \eqref{eq:tsHeavySS} beruhen, stellen sich die Verhältnisse bei der Berechnung von Prüfgrenzen dort anders dar als beim \ac{GKF}: \begin{IEEEeqnarray}{rCl} \pdfBr{\given{\measV_\indexS}{\measV_{1:\indexS-1}}} &=& \int \pdfBr{\given{\measV_{\indexS}}{\stateV_{\indexS}}} \: \pdfBr{\given{\stateV_{\indexS}}{\measV_{1:\indexS-1}}} \dif \stateV_{\indexS} = \IEEEnonumber \\ &=& \int \pdftBr{\given{\measV_{\indexS}}{ \outputM \stateV_\indexS + \feedM \inputV_\indexS, \measNCM, \dofMeasN }} \: \pdfNBr{\given{\stateV_{\indexS}}{ \est{\stateV}_{\given{\indexS}{\indexS-1}}, \stateCM_{\given{\indexS}{\indexS-1}} }} \dif \stateV_{\indexS}. \label{eq:tsStudentPredictionMarginalResult} \end{IEEEeqnarray} Eine einfache analytische Marginalisierung wie bei \eqref{eq:tsGKFPredictMeasurementAgainResult} ist in diesem Fall nicht mehr möglich. Es könnten nun wiederum analytische oder numerische Näherungsmethoden angewendet werden, um aus \eqref{eq:tsStudentPredictionMarginalResult} eine Dichtefunktion und daraus eine \ac{HDR} analog zum \ac{GKF} zu berechnen. Ein Beispiel für einen analytischen Ansatz findet sich bei Agamennoni \cite[][5027\psq]{Agamennoni2012}. Zu beachten ist allerdings, dass die Student-t-Verteilung des Messrauschens in \eqref{eq:tsHeavySSMeasN} gewählt wurde, um Ausreißer im stochastischen Modell der robusten Filteralgorithmen von Anfang an zu zu berücksichtigen. Das Zustandsraummodell \eqref{eq:tsHeavySS} umfasst also sowohl \ac{iO}- als auch \ac{niO}-Prüfläufe. Eine aus den Prädiktionen gemäß \eqref{eq:tsStudentPredictionMarginalResult} berechnete Prüfgrenze würde daher kein Abbild der \ac{iO}-Prüfläufe liefern, das Resultat wären deutlich zu weite Prüfgrenzen. Der Grund dafür liegt in den heavy Tails der Student-t-Verteilung. Als Ausweg besteht die Möglichkeit, die \ac{iO}-Messungen für die Prüfung als normalverteilt um den prädizierten Zustand $\stateV_{\given{ \indexS }{ \indexS-1 }}$ zu modellieren. Der Filtervorgang an sich beruht weiterhin auf der Student-t-Verteilung. Die genäherte Dichtefunktion der Messung bei gegebenem Zustand lautet dann \begin{IEEEeqnarray}{c} \pdfqSubBr{\test}{\given{\measV_{\indexS}}{\stateV_{\indexS}}} = \pdfNBr{\given{\measV_{\indexS}}{ \outputM \stateV_\indexS + \feedM \inputV_\indexS, \apr{\measNCM}_\test}}. \label{eq:tsStudentTestPredictionApprox} \end{IEEEeqnarray} Eine solche Näherung einer Student-t-Verteilung durch eine Normalverteilung kommt bereits beim Update-Schritt des \ac{TKF} zur Anwendung. Mit dem entsprechenden Skalierungsfaktor aus \eqref{eq:tsTKFRescaleStudentT} gilt \begin{IEEEeqnarray*}{c} \apr{\measNCM}_\test = \scaleDOF_{\infty,\dofMeasN}^2 \cdot \measNCM . \end{IEEEeqnarray*} % Die Marginalisierung über den Zustand liefert wie bei \eqref{eq:tsGKFPredictMeasurementAgainResult} nun eine analytisch berechenbare Dichtefunktion für die Prädiktion \begin{subequations} \label{eq:tsRobustBoundaryApprox} \begin{IEEEeqnarray}{rCl} \pdfqSubBr{\test}{\given{\measV_\indexS}{\measV_{1:\indexS-1}}} &=& \int \pdfqSubBr{\test}{\given{\measV_{\indexS}}{\stateV_{\indexS}}} \: \pdfBr{\given{\stateV_{\indexS}}{\measV_{1:k-1}}} \, \dif \stateV_{\indexS} = \pdfNBr{\given{\measV_{\indexS}}{\est{\measV}_{\given{\indexS}{\indexS-1} }, \apr{\measCM}_{\indexS,T}}}, \IEEEyessubnumber* \\ \est{\measV}_{\given{ \indexS }{ 1:\indexS-1 }} &=& \outputM \stateV_\indexS + \feedM \inputV_\indexS, \\ \apr{\measCM}_{\indexS,\test} &=& \outputM \stateCM_{\given{\indexS}{\indexS-1}} \trans{\outputM} + \apr{\measNCM}_\test . \end{IEEEeqnarray} \end{subequations} Die nötigen Verteilungsparameter für die Auswertung von \eqref{eq:tsRobustBoundaryApprox} liefern die Filter \ac{TKF}, \ac{RME} und \ac{ORKF} bereits. Es ist beim \ac{RME} und \ac{ORKF} lediglich einmalig vor Start des Filtervorgangs der Skalierungsfaktor $\scaleDOF_{\infty,\dofMeasN}$ zu bestimmen. %[Roth2013b, p. 5771]: "Quadratic forms $\trans{x} \Sigma^{-1} x$ … follow an F-distribution if p(x) = St(…), a fact which can be used to define probability regeions and gates in target tracking (similar to the use of the chi-squared distribution in the Gaussian case" %aber: Student-t gibt nicht die Verteilung der iO-Prüflinge wieder... %============================================================================== \subsection{Vergleich der Prüfgrenzen der einzelnen Filteralgorithmen} %============================================================================== Als Abschluss zur Diskussion des Verhaltens der Filter bei bekannten Systemparametern folgt nun ein Vergleich der Prüfgrenzen, die aus den Schätzergebnissen der vier Filteralgorithmen resultieren. Dafür wird wie zuvor ein skalarer Random-Walk genutzt, die Parameter zeigt \tabref{tab:tsSimNormalShort}. Um einen typischen Fall der \ac{EOL}-Prüfung mit vielen aufeinanderfolgenden \ac{iO}-Prüfläufen zu simulieren ist das Messrauschen nun jedoch normalverteilt. Für die Einstellung der Prüfschärfe wurde $\Prob_\test = \num{0.9973}$ gesetzt, die Filter nutzen $\dofMeasN = \num{5}$. Die Varianz $\measNC$ des normalverteilten Messrauschens wird für die auf Student-t basierenden Filter mittels der Skalierung $\scaleDOF_{\dofMeasN,\infty}^2 \cdot \measNCM$ auf die Modellannahmen angepasst. Eine manuell eingefügte längere Pause und ein manuell eingefügter Ausreißer zeigen das Verhalten in diesen beiden relevanten Fällen. \Figref{fig:tsRobustFilterTSBoundary} stellt die Ergebnisse der Filter und die resultierenden Prüfgrenzen dar. Die Prüfgrenzen wurden dabei durch zeitlich fein abgestufte Prädiktionsschritte zwischen den einzelnen Messungen mehrfach ausgewertet, um deren tatsächliche zeitliche Entwicklung besser darzustellen. Die Filter folgen wie bei den Untersuchungen in \secref{sec:tsFilterCompareSim} zuverlässig dem simulierten Zustand $\state_\indexS$. Für die normalen Messungen ist das Verhalten der Filter bezüglich der Zustandsschätzung und der erzeugten Prüfgrenzen praktisch identisch. Die verlängerte Prädiktionsphase zwischen $t = \num{10}$ und $t = \num{15}$ wird ebenfalls von allen Filtern ähnlich überbrückt. Es zeigt sich hier deutlich die ansteigende Unsicherheit bzgl. der Prädiktion der Messung bei $t = \num{15}$ anhand der sich aufweitenden Prüfgrenzen. Mit dem Update-Schritt bei $t = \num{15}$ schließen sich die Prüfgrenzen durch die aktualisierte Zustandsschätzung wieder. Den Ausreißer bei $t = \num{29}$ haben alle Filter erkannt. Das \ac{GKF} wird in seiner Zustandsschätzung durch den Ausreißer nicht beeinflusst, da dieser außerhalb des Validation-Gate liegt. Beim \ac{TKF} führt der Ausreißer wie in \secref{sec:tsFilterCompareSim} zu einer starken Verfälschung der Zustandsschätzung, die erst beim darauffolgenden Messwert ($t = \num{30}$) korrigiert wird. Die sprungartig ansteigende Kovarianz der Schätzung führt zu sehr weiten Prüfgrenzen, so dass beim \ac{TKF} keine ausreichende Robustheit der Prüfgrenze gegen Ausreißer vorliegt. Der \ac{RME} und das \ac{ORKF} zeigen hingegen ein robustes Verhalten, die Prüfgrenze weicht trotz des Ausreißers nur geringfügig vom Verhalten des \ac{GKF} ab. % EC1112 05_01 \begin{table} \caption[FIXME]{Parameter für die Simulation eines kurzen Random-Walks mit normalverteiltem Messrauschen} \label{tab:tsSimNormalShort} \centering \footnotesize \begin{tabular}{ll} \toprule Parameter & Wert \\ \midrule $\procNC$ & \num{0.1} \\ $\measNC$ & \num{0.1} \\ $\state_1$ & \num{0} \\ $\Delta t_{\indexS}$ & \num{1} \\ $\stateG$ & \num{1} \\ $\inputG$ & \num{0} \\ $\outputG$ & \num{1} \\ $\feedG$ & \num{0} \\ \bottomrule \end{tabular} \end{table} % ------------------------------------------------------------------- \begin{figure} \centering \footnotesize \setlength\figurewidth{.9\columnwidth} \setlength\figureheight{.25\columnwidth} \renewcommand\figureXLabel{$t$} \renewcommand\figureYLabel{$\state_\indexS, \meas_\indexS$} \subcaptionbox{\ac{GKF}}[\columnwidth]{% \centering \footnotesize \mytikzexton \input{"graphics/P1112 05 01/002 ICIT Gaussian Walk GKF.tikz"} \mytikzextoff }% % --------------- comment magic (tm) \\% \vspace{0mm}% % --------------- comment magic (tm) \subcaptionbox{\ac{TKF}}[\columnwidth]{% \centering \footnotesize \mytikzexton \input{"graphics/P1112 05 01/005 ICIT Gaussian Walk TKF.tikz"} \mytikzextoff }% % --------------- comment magic (tm) \\% \vspace{0mm}% % --------------- comment magic (tm) \subcaptionbox{\ac{RME}}[\columnwidth]{% \centering \footnotesize \mytikzexton \input{"graphics/P1112 05 01/003 ICIT Gaussian Walk RKF.tikz"} \mytikzextoff }% % --------------- comment magic (tm) \\% \vspace{0mm}% % --------------- comment magic (tm) \subcaptionbox{\ac{ORKF}}[\columnwidth]{% \centering \footnotesize \mytikzexton \input{"graphics/P1112 05 01/004 ICIT Gaussian Walk ORKF.tikz"} \mytikzextoff }% \caption[FIXME]{Simulierter Random-Walk mit Messwerten $\meas_\indexS$ \circleWithParens{myLineOne} und Zuständen $\state_\indexS$ \lineWithParens{myLineOne} sowie einem künstlich eingefügter Ausreißer \circleWithParens{myNIO}. Die Filterergebnisse sind in Form der Prädiktionen $\est{\meas}_{\given{ \indexS }{ 1:\indexS-1 }}$ \dashDotLineWithParens{black} und Prüfgrenzen \lineWithParens{black} dargestellt. Der schattierte Bereich \squareWithParens{myGray75} stellt den Bereich für \ac{iO}-Prüfergebnisse dar.} \label{fig:tsRobustFilterTSBoundary} \end{figure} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Parameterschätzung für robuste Zustandsraummodelle} \label{sec:tsParamEst} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% In der Literatur gilt bei der der Herleitung von Filteralgorithmen wie in \secref{sec:tsRobustFilter} und \secref{sec:tsTestLimits} dieser Arbeit üblicherweise die Annahme, dass die Modellparameter bekannt sind. Die Parameter umfassen bei den stochastischen Zustandsraummodellen auch die Verteilungsparameter der Rauschterme. Bei Systemen mit bekannten inneren physikalischen Zusammenhängen ist die Annahme bekannter Parameter vertretbar, wenn durch Systemkenntnis z.B. aus Konstruktionsdaten auf die Modellparameter geschlossen werden kann. Im Rahmen der \ac{EOL}-Prüfung operieren die verwendeten Zustandsraummodelle aus \secref{sec:tsRobustFilter} auf Messdaten in Form von Prüfmerkmalen, die durch ein komplexes Zusammenspiel aus Fertigungs- und Montageanlagen und Prüfständen entstehen. Daher ist neben der Zustandsschätzung und der Berechnung der Prüfgrenzen die Schätzung der Modellparameter unverzichtbar. Bei den nachfolgenden allgemeinen Betrachtungen fasst der Parametervektor \begin{IEEEeqnarray*}{c} \paramV \in \msbegin{\stateM, \inputM, \outputM, \feedM, \procNCM, \measNCM} \end{IEEEeqnarray*} alle unbekannten Parameter des Zustandsraummodells zusammen. Sollten einzelne Parameter vorab bekannt sein, können diese aus $\paramV$ ausgeschlossen und fest vorgegeben werden. Die Parameterschätzung kann für die \ac{EOL}-Prüfung offline erfolgen, siehe dazu die Struktur des Prüfsystems in \secref{sec:tsSystemScheme} und generelle Aufteilung Offline-Online in \secref{sec:EOLSplit}. %============================================================================== \subsection{Übersicht zu Möglichkeiten der Parameterschätzung} %============================================================================== Die grundlegende Problematik bei der Parameterschätzung für Zustandsraummodelle ist, dass die Realisierungen der Zustände $\stateV_\indexS$ auch in der Offline-Parameterschätzung nicht bekannt sind. Wären die $\stateV_\indexS$ bekannt, dann könnte beispielsweise die Parameterschätzung für das linear-gaußsche Zustandsraummodell als einfaches Regressionsproblem aufgefasst werden \cite[][647]{Murphy2013}. Für den Zustandsübergang gilt dann \begin{IEEEeqnarray*}{c} \mvbegin{\trans{\stateV}_2 \\ \trans{\stateV}_3 \\ \vdots \\ \trans{\stateV}_\numS } = \mmbegin{ \trans{\stateV}_1 & \trans{\inputV}_1 \\ \trans{\stateV}_2 & \trans{\inputV}_2 \\ \vdots & \vdots \\ \trans{\stateV}_{\numS-1} & \trans{\inputV}_{\numS-1}} \cdot \mvbegin{\trans{\stateM} \\ \trans{\inputM}} + \mvbegin{\trans{\procNV}_1 \\ \trans{\procNV}_2 \\ \vdots \\ \trans{\procNV}_{\numS-1}}. \end{IEEEeqnarray*} Dies entspricht bereits der Struktur eines Regressionsproblems für die zu schätzenden Parameter $\stateM$ und $\inputM$. Für die Schätzung von $\outputM$ und $\feedM$ lässt sich zwischen $\stateV_\indexS$ und $\measV_\indexS$ ein äquivalenter Zusammenhang aufstellen. Die Kovarianzmatrizen $\procNCM$ und $\measNCM$ der Rauschterme können so ebenfalls bestimmt werden. Der Ansatz scheitert allerdings an den unbekannten Zustandsvektoren $\stateV_\indexS$. %Außerdem wächst die Größe des Regressionsproblems mit der Anzahl der Messungen, was die rechnerische Lösung immer schwieriger (und langsamer) macht. Bereits Kalman hat bei der Vorstellung des Kalman-Filters \cite{Kalman1960} als rekursiven Zustandsschätzer angemerkt, dass der Filtervorgang und die Schätzung der Modellparameter nach Möglichkeit gemeinsam zu betrachten sind. Kalman hat dieses Problem jedoch offen gelassen. Seitdem wurden eine Reihe von Techniken zur Parameterschätzung für Zustandsraummodelle vorgestellt. Die nachfolgende Übersicht basiert zum Teil auf \cite[][Kapitel 12]{Sarkka2013} und \cite[][646\psq]{Murphy2013}. %----------------------------------------------------------------------------- \subsubsection{State-Augmentation} %----------------------------------------------------------------------------- Eine zunächst einfach erscheinende Möglichkeit ist, den Zustandsvektor $\stateV_\indexS$ um den Parametervektor $\paramV_\indexS$ zum Zeitpunkt $\indexS$ zu erweitern. Dieser Ansatz ist als State-Augmentation \cite[][281-285]{Jazwinski1970} oder Joint-Kalman-Filtering \cite[][64-70]{Nelson2000} bekannt und erlaubt auch die Online-Schätzung zeitveränderlicher Modellparameter. Dazu werden die zu $\paramV$ korrespondierenden Zustände im erweiterten Zustandsvektor mit einem Random-Walk-Modell in das Zustandsraummodell eingebunden. Da die Elemente von $\paramV$ die Dynamik- und Messgleichung des ursprünglichen Zustandsvektors $\stateV_\indexS$ vorgeben und nun aber selbst Teil des erweiterten Zustandsvektors sind, ist auch bei linearen Modellen für die Nutzung der State-Augmentation ein Zustandsschätzer für nichlineare Modelle wie z.B. ein Extended-Kalman-Filter erforderlich. Ähnlich zur State-Augmentation funktioniert das Dual-Kalman-Filtering. Hier arbeiten zwei Filter parallel, eines für die Schätzung der Zustände und eines für die Schätzung der Modellparameter. Die Filter nutzten die Schätzergebnisse des jeweils anderen Filters als Parameter in der jeweils eigenen Schätzung. Der Ansatz wird dadurch verkompliziert, dass die Schätzunsicherheit des jeweils anderen Filters zu berücksichtigen ist \cite[][70-82]{Nelson2000}. %In seiner allgemeinen Form lautet das Problem der gemeinsamen A-posteriori-Schätzung von Parametern und Zuständen \cite[][174\psq]{Sarkka2013} %\begin{IEEEeqnarray*}{c} %\pdfBr{\given{ \stateV_{1:\numS}, \paramV }{ \measV_{1:\numS} }} = \frac{ \pdfBr{\given{ \measV_{1:\numS} }{ \stateV_{1:\numS}, \paramV }} \: \pdfBr{\given{ \stateV_{1:\numS} }{ \paramV }} \: \pdfBr{\paramV} }{ \pdfBr{\measV_{1:\numS}} } %\end{IEEEeqnarray*} %mit %\begin{IEEEeqnarray*}{rCl} %\pdfBr{\given{ \stateV_{1:\numS} }{ \paramV }} &=& \pdfBr{\given{ \stateV_1 }{ \paramV }} \prod_{\indexS=2}^{\numS}{\pdfBr{\given{ \stateV_\indexS }{ \stateV_{\indexS-1}, \paramV }}}, \\ %\pdfBr{\given{ \measV_{1:\numS} }{ \stateV_{1:\numS}, \paramV }} &=& \prod_{\indexS=1}^{\numS}{\pdfBr{\given{ \measV_\indexS }{ \stateV_\indexS, \paramV }}}. %\end{IEEEeqnarray*} %Volle A-posteriori-Verteilung aufwendig %----------------------------------------------------------------------------- \subsubsection{Maximum-Likelihood- und Maximum-a-posteriori-Schätzung} % Achtung keine subsection %----------------------------------------------------------------------------- Für die A-posteriori-Verteilung der Modellparameter gilt \cite[][174-176]{Sarkka2013} \begin{IEEEeqnarray}{c} \pdfBr{\given{ \paramV }{ \measV_{1:\numS} }} = \frac{ \pdfBr{\given{ \measV_{1:\numS} }{ \paramV }} \: \pdfBr{\paramV} }{ \pdfBr{\measV_{1:\numS}} } = \frac{ \pdfBr{\given{ \measV_{1:\numS} }{ \paramV }} \: \pdfBr{\paramV}}{ \int \pdfBr{\given{ \measV_{1:\numS} }{ \paramV }} \: \pdfBr{\paramV} \: \dif \paramV } . \label{eq:tsParamPosterior} \end{IEEEeqnarray} Für die Ermittlung einer \ac{MAP}- bzw. \ac{ML}-Lösung ist die Auswertung der marginalen Likelihood $\pdfBr{\given{ \measV_{1:\numS} }{ \paramV }}$ von zentraler Bedeutung. Angestrebt wird eine rekursive Berechnung, um den Berechnungsaufwand pro Messung unabhängig von der Gesamtanzahl der Messungen konstant zu halten. Die marginale Likelihood lässt sich zunächst ohne weitere Annahmen zerlegen in \begin{IEEEeqnarray}{c} \pdfBr{\given{ \measV_{1:\numS} }{ \paramV }} = \pdfBr{\given{ \measV_1 }{ \paramV }} \prod_{\indexS=2}^{\numS}{\pdfBr{\given{ \measV_\indexS }{ \measV_{1:\indexS-1}, \paramV }}}. \label{eq:tsParamProductOfLikelihoods} \end{IEEEeqnarray} Die Produktterme in \eqref{eq:tsParamProductOfLikelihoods} enthalten immer noch eine Konditionierung auf alle vorangegangenen Messungen. Der nächste wichtige Zwischenschritt ist die Nutzung der Markov-Eigenschaft des Zustandsvektors, in dem beim Zustandsraummodell alle Informationen der vorangegangenen Messungen stecken: \begin{IEEEeqnarray}{rCl} \pdfBr{\given{ \measV_\indexS }{ \measV_{1:\indexS-1}, \paramV }} &=& \int{ \pdfBr{\given{ \measV_\indexS }{ \stateV_\indexS, \measV_{1:\indexS-1}, \paramV }} \: \pdfBr{\given{ \stateV_\indexS }{ \measV_{1:\indexS-1}, \paramV }} \dif \stateV_\indexS} \IEEEnonumber \\ &=& \int{ \pdfBr{\given{ \measV_\indexS }{ \stateV_\indexS, \paramV }} \: \pdfBr{\given{ \stateV_\indexS }{ \measV_{1:\indexS-1}, \paramV }} \dif \stateV_\indexS} . \label{eq:tsParamMarkov} \end{IEEEeqnarray} Die beiden Dichtefunktionen im Integral von \eqref{eq:tsParamMarkov} sind bereits von den rekursiven Filteralgorithmen bekannt, wobei jetzt die Abhängigkeit vom zu bestimmenden Parametervektor $\paramV$ explizit angeführt ist. Die Dichtefunktion $\pdfBr{\given{ \measV_\indexS }{ \stateV_\indexS, \paramV }}$ entspricht der Messgleichung des Zustandsraummodells, und $\pdfBr{\given{ \stateV_\indexS }{ \measV_{1:\indexS-1}, \paramV }} $ ist prädiktive Dichte für den Zustandsvektor. Insgesamt wurde $\pdfBr{\given{ \measV_\indexS }{ \measV_{1:\indexS-1}, \paramV }}$ bereits für die Herleitung der Prüfgrenzen in \eqref{eq:tsGKFPredictMeasurementAgainResult} bzw. \eqref{eq:tsStudentPredictionMarginalResult} genutzt. Die Zerlegung \eqref{eq:tsParamProductOfLikelihoods} in Kombination mit \eqref{eq:tsParamMarkov} ermöglicht also durch Nutzung eines zum Zustandsraummodell passenden rekursiven Filters die ebenfalls rekursive Auswertung der marginalen Likelihood \eqref{eq:tsParamProductOfLikelihoods}. Dieses Vorgehen ist heute als Prediction-Error-Decomposition \cite[][176]{Sarkka2013} oder Innovations-Form \cite[][335]{Shumway2011} der marginalen Likelihood bekannt, eine frühe Referenz findet sich bei Schweppe \cite{Schweppe1965}. Es sei angemerkt, dass die konkrete Form der Dichtefunktionen hier noch nicht eingeflossen ist und die Gleichungen daher auch für die robusten Filter gelten. Für die praktische Auswertung der Integrale sind bei nicht-normalverteilten Größen gegebenfalls wieder Näherungen erforderlich. Für die Umsetzung der Parameterschätzung auf Basis der Prediction-Error-Decomposition wird die \ac{MAP}-Lösung \begin{subequations} \label{eq:tsParamDirectMAP} \begin{IEEEeqnarray}{rCl} \est{\paramV}_{\text{\acs{MAP}}} &=& \arg \max_{\paramV} \costFunSubBr{\text{\acs{MAP}}}{\paramV} , \IEEEyessubnumber* \\ \costFunSubBr{\text{\acs{MAP}}}{\paramV} &=& \pdfBr{\given{ \measV_{1:\numS} }{ \paramV }} \: \pdfBr{\paramV} = \pdfBr{\paramV} \: \pdfBr{\given{ \measV_1 }{ \paramV }} \prod_{\indexS=2}^{\numS}{\pdfBr{\given{ \measV_\indexS }{ \measV_{1:\indexS-1}, \paramV }}} \end{IEEEeqnarray} \end{subequations} durch einen numerischen Optimierer ermittelt. Durch formelles Setzen von $\pdfBr{\paramV} = 1$ ergibt sich die \ac{ML}-Lösung \begin{subequations} \label{eq:tsParamDirectML} \begin{IEEEeqnarray}{rCl} \est{\paramV}_{\text{\acs{ML}}} &=& \arg \max_{\paramV} \costFunSubBr{\text{\acs{ML}}}{\paramV} , \IEEEyessubnumber* \\ \costFunSubBr{\text{\acs{ML}}}{\paramV} &=& \pdfBr{\given{ \measV_{1:\numS} }{ \paramV }} = \pdfBr{\given{ \measV_1 }{ \paramV }} \prod_{\indexS=2}^{\numS}{\pdfBr{\given{ \measV_\indexS }{ \measV_{1:\indexS-1}, \paramV }}} . \label{eq:tsParamDirectMLCostFun} \end{IEEEeqnarray} \end{subequations} Eine Variante der \ac{MAP}- bzw. \ac{ML}-Schätzung liefert die \ac{EM}-Methode \cite{Dempster1977}. Sie lässt sich allgemein zur Schätzung bei Modellen mit nicht messbaren (latenten) Variablen einsetzten. Für die Anwendung zur Parameterschätzung bei Zustandsraummodellen benötigt die \ac{EM} einen Glättungs-Algorithmus, während die direkte numerische Optimierung von \eqref{eq:tsParamDirectMAP} bzw. \eqref{eq:tsParamDirectML} mit den Ergebnissen eines Filters arbeitet. Im Gegenzug kommt die \ac{EM} ohne numerischen Optimierer aus. Details für die Implementierung bei Zustandsraummodellen finden sich bei Särkkä \cite[][182-185]{Sarkka2013}. % [Barber2013, Ch. 24.5] EM mit teilweiser Herleitung, komplette Update-Formeln für Lineare Dynamische Systeme %----------------------------------------------------------------------------- \subsubsection{Approximative Bayes-Schätzer} %----------------------------------------------------------------------------- Während die oben Maximum-Likelihood- und Maximum-a-posteriori-Verfahren eine Punktschätzung für die Modellparameter liefern, ist auch eine vollwertige Auswertung der A-posteriori-Verteilung \eqref{eq:tsParamPosterior} zur Parameterschätzung möglich. Eine analytische Auswertung scheitert jedoch spätestens an der Integration des Normalisierungsfaktors im Nenner von \eqref{eq:tsParamPosterior}. Gängige Approximationsverfahren für die A-posteriori-Verteilung stammen aus der Klasse der \ac{MCMC}-Algorithmen. Sie erzeugen Samples aus der A-posteriori-Verteilung, die diese asymptotisch nachbilden \cite{Andrieu2010}. Wie bereits bei den Filtern für nicht linear-gaußsche Modelle können die Monte-Carlo-basierten Methoden zwar beliebig genau gemacht werden, der numerische Rechenaufwand steigt dafür allerdings unter Umständen über das praktisch vertretbare Maß hinaus an. Eine Alternative zum numerischen Sampling bietet wiederum die \ac{BV} mit ihren analytischen Näherungsansätzen für die A-posteriori-Verteilung. Für die Parameterschätzung von Zustandsraummodellen wurden \ac{BV}-Methoden u.a. von Beal \cite[][Kapitel 5]{Beal2003} und Barber \cite{Barber2007} vorgestellt. %----------------------------------------------------------------------------- \subsubsection{Korrelations-basierte und kovarianz-basierte Methoden} %----------------------------------------------------------------------------- Die stochastischen Filter liefern eine Reihe von Statistiken, die nur erfüllt sind wenn die Modellannahmen (und die zugehörigen Parameterwerte) dem realen Systemverhalten entsprechen. Dazu zählt beim Kalman-Filter die zeitliche Unkorreliertheit der Residuen und die ebenfalls vom Filter gelieferte Kovarianz der Schätzgrößen. Korrelations-basierte Methoden nutzen die Autokorrelationsfunktion des Systemausgangs oder der Residuen, um auf die Modellparameter zu schließen \cite{Mehra1972, Odelson2006}. Kovarianz-basierte Methoden werten hingegen in einem zeitlichen Fenster die Sample-Kovarianzmatrizen der Innovationen und Zustandsschätzungen aus \cite{Mehra1972,Hashlamon2016} %----------------------------------------------------------------------------- \subsubsection{Multiple-Model-Approach} %----------------------------------------------------------------------------- Anstatt die Modellparameter einzustellen besteht auch die Möglichkeit, mehrere Filter mit verschiedenen Parametersätzen parallel laufen zu lassen. Basierend auf den Residuen der einzelnen Filter und deren Likelihood wird jenes Filter ausgewählt, das am besten das realen Systemverhalten wiedergibt. Auch gewichtete Kombinationen der einzelnen Filterergebnisse sind möglich \cite[][441-466]{BarShalom2001}. %----------------------------------------------------------------------------- \subsubsection{Subspace-Methode} %----------------------------------------------------------------------------- In einem Zustandsraummodell ohne Rauschen und externe Eingänge führen die System- und Ausgangsmatrix in jedem Zeitschritt eine lineare Projektion des akutellen Zustandsvektors $\stateV_\indexS$ aus, der wiederum aus einer mehrfachen linearen Projektionen des initialen Zustandsvektors $\stateV_1$ hervorging. Die Zustände $\stateV_\indexS$ und Messungen $\measV_\indexS$ können sich also nur innerhalb von Unterräumen befinden, die durch die Modellparameter festgelegt sind. Über diese Unterräume schließt die Subspace-Methode nun durch eine Singulärwertzerlegung auf die Modellparameter \cite[][Abschnitt 24.5.3]{Barber2013}. %------------------ % Reste Parameterschätzung %------------------ %[Barber2013, Ch. 24.5] Kapitel zu Learning Linear Dynamical Systems %EM mit teilweiser Herleitung, komplette Update-Formeln für Lineare Dynamische Systeme => "The statistics required therefore include smoothed means, covariances, and cross moments" % %Vergleich EM und Gradienten-basierte Optimierung in [Cappe2005, p. 358] %• "EM approach is more generally known" %• Optimization method (line-search etc.) needed by gradient-based methods %• "The EM algorithm often deals with parameter constraints implicitly"; "For gradient-based methods this is not the case, and parameter constraints have to be dealt with explicitly" %• "The EM algorithm is parameterization independent"; Reparametrizations change convergence behavior of gradient-based methods %• "Gradient-based methods do not require the M-step. Thus they may be applied to models for which the M-step does not lead to simple closed-form solutions." %• "Gradient-based methods converge faster. As discussed above, gradientbased methods can reach quadratic convergence whereas EM usually converges only linearly." % %zahlreiche Quellen zu State Space / Probabilistic Models: %[McGoff2015] Statistical inference for dynamical systems: A review. Statistics Surveys, 9:209–252, 2015. % %[Durbin2012, Ch. 7] Maximum likelihood estimation of parameters %Übliche Berechnung der Log-Likelihood als Prediction Error Decomposition (7.2) (p. 171) mit Quellen Schweppe1965 und Harvey1989 %Likelihood with diffuse filter initialization (7.3) (p. 171) und auch Kommentar p. 174 oben % %[Schon2011b] Lösungsmöglichkeiten: Numerische Ableitungsberechnung oder ableitungsfreie Verfahren, inkl. Referenzen [44, 48]; Alternative: EM -> Ableitungsfrei %Motivation: %S. Sangsuk-Iam and T. Bullock, “Analysis of discrete-time Kalman filtering under incorrect noise covariances,” IEEE Trans. Automat. Contr., vol. 35, no. 12, pp. 1304–1309, Dec. 1990. %(Referenz gefunden in [Ardeshiri2015]) % %\texttodo{Auswirkungen falscher Parameterwerte aufzeigen (Kovarianz der Parameter zu klein, Anpassung an Änderungen zu langsam); Beispiel RKF auf Random-Walk mit Q \num{0.001} vom realen Wert} % %Initial Conditions for KF -> beeinflusst Parameterschätzung; z.B. in [Harvey1990], [Jalles2009], [Durbin2012 Kap. 5], [Murphy2013, p. 646] % %Abbildung "Graphical Model" mit Übersicht zu Filter/Smoother (Inference) <-> Parameter-Schätzung (Learning) -> Prüfen (Predictive) % %Siehe Einleitung Särkkä2013 (Kapitel 1.5 "Parameter Estimation") %"Unfortunately, computing this joint posterior of the states and parameters is even harder than computation of the joint distribution of states alone, and thus this task is intractable." %Ganz fundamentale Sicht auf Marginal Posterior of Parameters ("warum geht das nicht direkt?"): [Särkkä2013, p. 175] => Näherungen ohne Joint Posterior States and Parameters aufzustellen %[Särkkä2013, Ch. 12]: optimization based MAP/ML; EM; MCMC; jeweils mit Kalman Filter/Smoother, Gaussian Filter/Smoother, Particle Filter/Smoother %[Särkkä2009, p. 596]: %"The classical way (see, e.g., [8], [9]) of solving the problem of uncertain parameters is to use adaptive filters where the model parameters or the noise statistics are estimated together with the dynamic state. The classical noise adaptive filtering approaches can be divided into Bayesian, maximum likelihood, correlation and covariance matching methods [8]. The Bayesian approach is the most general of these and the other approaches can often be interpreted as approximations to the Bayesian approach. Examples of Bayesian approaches to noise adaptive filtering are state augmentation based methods [10], [11], multiple model methods such as the interacting multiple models (IMM) algorithm [3], [9] and particle methods [12]–[15]." %* interessant für Vergleiche von Methoden %* Quelle [8]: R. Mehra, “Approaches to adaptive filtering,” IEEE Trans. Automat. Control, vol. AC-17, no. 5, pp. 693–698, Oct. 1972. -> siehe Referenzen KF %* Quelle [9]: X. R. Li and Y. Bar-Shalom, “A recursive multiple model approach to noise identification,” IEEE Trans. Aerosp. Electron. Syst., vol. 30, no. 3, pp. 671–684, Jul. 1994. -> siehe Referenzen KF %* Quelle [10]: P. Maybeck, Stochastic Models, Estimation and Control, Volume 2. New York: Academic Press, 1982. -> siehe Referenzen KF %Quelle [11]: E. A. Wan and A. T. Nelson, “Dual EKF methods,” in Kalman Filtering and Neural Networks, S. Haykin, Ed. New York: Wiley, 2001, ch. 6. -> siehe Referenzen KF %Quelle [12]: A. Doucet, N. de Freitas, and N. Gordon, Sequential Monte Carlo Methods in Practice. New York: Springer, 2001. %Quelle [13]: G. Storvik, “Particle filters in state space models with the presence of unknown static parameters,” IEEE Trans. Signal Processing, vol. 50, no. 2, pp. 281–289, Feb. 2002. %Quelle [14]: P. M. Djuric and J. Miguez, “Sequential particle filtering in the presence of additive Gaussian noise with unknown parameters,” in Proc. IEEE Int. Conf. Acoust., Speech, Signal Processing, Orlando, FL, 2002, pp. 1621–1624. %Quelle [15] S. Särkkä, “On sequential Monte Carlo sampling of discretely observed stochastic differential equations,” in Proc. Nonlin. Stat. Signal Processing Workshop, Sep. 2006, [CD ROM]. %"Although, the VB-AKF method presented in this article is a Bayesian method for noise identification, the algorithm much resembles the covariance matching methods presented in, for example, [8], [10]. This is because the algorithm essentially estimates noise parameters such that the parameters become consistent with the observed residual, which is also the basis of covariance matching methods." %------------------ % zu einzelnen Methoden konkret: State augmentation %------------------ %State Augmentation / Dual Estimation %[BarShalom2001, Kapitel 11.9] Augmentation of the State, EKF for Parameter Estimation %Dissertation [Nelson2000], siehe dazu Notizen KF %Hauptargument gegen State Augmentation: keine statistisch solide Basis, "Dynamik" der Parameter wird mit Dynamik der Zustände vermischt ; kann divergieren [Särkkä2013, p. 186] %------------------ % zu einzelnen Methoden konkret: EM %------------------ %Ergebnis: MAP der Modellparameter, hier: Kovarianzen % %[Murphy2013, p. 647] EM for LG-SSM %If we only observe the output sequence, we can compute ML or MAP estimates of the parameters %using EM. The method is conceptually quite similar to the Baum-Welch algorithm for HMMs %(Section 17.5), except we use Kalman smoothing instead of forwards-backwards in the E step, %and use different calculations in the M step. We leave the details to Exercise 18.1. % %Klassiker A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm,” J. Roy. Statist. Soc. B (Methodological), vol. 39, no. 1, pp. 1–38, 1977. %-> Parameter-Schätzung für State-Space als incomplete / missing data problem %laut [Särkkä2013]: braucht Smoother % %EM für State Space konkret in [Särkkä2013, p. 183]; Berechnung nötigen Erwartungswerte für Q bei LGS in [Särkkä2013, p. 189]; Konkrete Update-Formeln für LGS [Särkkä2013, p. 191]; für nichtlineare Systeme Näherung mit Gaussian assumed density approximation [Särkkä2013, p. 194] % %Knackige Beschreibung EM mit schneller Herleitung für State Space: [Shumway2011, p. 340ff] %Complete data likelihood (6.63) -> iterative method EM; Smoother für E-Step % %EM braucht keine Ableitungsberechnung der Likelihood [Cappe2005b, p. 703] %Aber: M-Schritt enthält eine Optimierung -> Ableitungen der Auxiliary Function Q, vgl. [Särkkä2013, p. 184]; speziell für nichtlineare Systeme kann eine numerische Lösung innerhalb des M-Schrittes nötig sein [Särkkä2013, p. 194] -> direkte Optimierung u.U. sinnvoller % %EM für Local Level Model in [Commandeur2007, p. 150] % %EM für Lineare Dynamische Systeme [Barber2013, Ch. 24.5] s.o. %EM auch in Hill2009 Real-time Bayesian anomaly detection in streaming environmental data S. 3 -> tends to suboptimal solutions % %[Agamennoni2011] siehe dort % %EM für State Space: %Roweis, S. and Ghahramani, Z. 2001. Learning nonlinear dynamical systems using the expectation–maximization algorithm. Chapter 6, pages 175–220 of: Haykin, S. (ed.), Kalman Filtering and Neural Networks. Wiley-Interscience. %G. C. Goodwin and J. C. Aguero. Approximate EM algorithms for parameter and state estimation in nonlinear stochastic models. In Proceedings of the 44th IEEE conference on decision and control (CDC) and the European Control Conference (ECC), pages 368–373, Seville, Spain, December 2005. % %[Schon2011b] %Sec. 4 Erklärung EM über Complete Data Log Likelihood von der Basis weg; Q als Minimum-Variance Estimate = Erwartungswert der Complete Data Log Likelihood; %Schöne Herleitung Anstieg der Log Likelihood bei EM % %EM allgemein: %EM Herleitung allgemein siehe z.B. [Barber2013, p. 253], eigene Notizen %------------------ % zu einzelnen Methoden konkret: VB %------------------ %[Murphy2013, p. 620] (VB-?) EM %The details for the HMM case can be found in [MacKay1997], [Beal2003], but the basic idea is this: The E step uses forwards-backwards, but where (roughly speaking) we %plug in the posterior mean parameters instead of the MAP estimates. The M step updates the %parameters of the conjugate posteriors, instead of updating the parameters themselves. % %[Ardeshiri2015], [Beal2003], [Särkkä2009] %D. Barber and S. Chiappa, “Unified inference for variational Bayesian linear Gaussian state-space models,” in Advances in Neural Information Processing Systems 19, 2007, pp. 81–88, MIT Press. %[Sarkka2013b] Non-Linear Noise Adaptive Kalman Filtering via Variational Bayes %[Särkkä2013c] Variational Bayesian Adaptation of Noise Covariances in Non-Linear Kalman Filtering %[Ardeshiri2015c] Analytical Approximations for Bayesian Inference Thesis % %[Ostwald2014] - A tutorial on variational Bayes for latent linear stochastic time-series models, inkl. Supplement % %Hoffman, M. D., Blei, D. M., Wang, C., and Paisley, J. (2013). Stochastic variational inference. The Journal of Machine Learning Research, 14(1):1303–1347. 27 %------------------ % zu einzelnen Methoden konkret: MCMC %------------------ %[Murphy2013, p. 620] MCMC %An alternative to VBEM is to use MCMC. A particularly appealing algorithm is block Gibbs %sampling, which we discuss in general terms in Section 24.2.8. The details for the HMM case %can be found in [FruhwirthSchnatter2006], but the basic idea is this: we sample z1:T given %the data and parameters using forwards-filtering, backwards-sampling, and we then sample the %parameters from their posteriors, conditional on the sampled latent paths. This is simple to %implement, but one does need to take care of unidentifiability (label switching), just as with %mixture models (see Section 11.3.1). % %Allgemeine, einfache Beschreibung von MCMC in [Shumway2011, p. 388] % %Referenz in [Särkkä2013] zu MCMC für HMM: [Ninness2010] Bayesian system identification via Markov chain Monte Carlo techniques %Allgemeines daraus zu Bayesian System Identification s.o.; Grundidee MCMC für Bayesian System Identification siehe Sec. 4 mit Metropolis, Metropolis-Hastings konkret für Systemidentifikation; Sec. 5-6 Analyse Metropolis-Hastings (stationäre Verteilung, Konvergenzrate etc.); Sec. 7 konkretes Anwendungsbeispiel, Vergleich mit numerischer Integration statt MCMC, Auswirkungen der Identifikations-Unsicherheit (Bayes) auf einen Regelkreis => MCMC-Samples der Regelstrecken-Parameter-Posterior direkt einsetzen um Amplituden- und Phasenrand zu bestimmen %Im Anhang Herleitungen zu Markov chains on uncountable spaces, … %[Wills2012] nochmal unter Beteiligung von Ninness, hier mit Blocked Gibbs Sampling, s.u. % %Metropolis-Hastings über Energy Function: [Särkkä2013, p. 188] %MCMC mit genähertem Filter (Gaussian filtering for nonlinear systems) [Särkkä2013, p. 193] -> Posterior der Parameter ist auch eine Näherung, typischerweise zu schmal (siehe [Särkk2013, Fig. 12.2 p. 197]) % %Particle Filter innerhalb MCMC: Particle MCMC und Particle Marginal Metropolis Hastings : %Andrieu, C., Doucet, A., Singh, S., and Tadic, V. 2004. Particle methods for change detection, system identification, and control. Proceedings of the IEEE, 92(3), 423–438. %Andrieu, C., Doucet, A., and Holenstein, R. 2010. Particle Markov chain Monte Carlo methods. The Royal Statistical Society: Series B (Statistical Methodology), 72(3), 269–342. %------------------ % Blocked Gibbs Sampling %------------------ %-> Perfekt für sampling-basierte Näherungen Filter/Smoother, s.u. %vgl. [Brodersen2015], [Wills2012] %[Shumway2011, p. 387] "This technique was first used by Carlin et al. (1992) in the context of general nonlinear and non-Gaussian state-space models. Frühwirth-Schnatter (1994) and Carter and Kohn (1994) built on these ideas to develop efficient Gibbs sampling schemes for more restrictive models." %Einfache Beschreibung Gibbs Sampling [Shumway2011, p. 389] % %Forward-Filtering Backward-Sampling ala Frühwirth-Schnatter1994: [Shumway2011, p. 391] %Für nichtlineare / nicht-Gaußsche Modelle, speziell nicht-Gauß mit scale-mixture-Rauschtermen, einschließlich Student-t: [Shumway2011, p. 392] %"Many examples of these techniques are given in Carlin et al. (1992), including the problem of model choice." %rel. ausführliche Beispiele zu nonlinear Noise, Beispiel zu t-Prozessrauschen % %Für Forward-Filtering, Backward-Sampling ala Frühwirth-Schnatter1994 auch [Durbin2012, p. 107] "draw random samples of the disturbance vectors…, and the sate vector" %"Subsequently, in Durbin and Koopman (2002), we developed a method which is absed only on mean corrections of unconditional vectors and which is much simpler and computationally more efficient than the de Jon-Shepard and earlier (read: Frühwirth-Schnatter) procedures" % %[Wills2012] %• Für Sampling aus Zuständen: a) Nutzung der Multivariaten Normalverteilung zum Sampling aus $p(x_1:N | y_1:N, \Theta)$ grundsätzlich möglich, aber für große N nicht machbar wegen Matrixinvertierungen Standard-Kalman-Smoother funktioniert auch nicht, da nur marginale $p(x_k | y_1:N, \Theta)$ geliefert werden und nicht gemeinsame Dichte $p(x_1:N | y_1:N, \Theta)$ %$p(x_1:N-1, x_N | y_1:N, \Theta) = p(x_N | x_1:N-1, y_1:N, \Theta) * p(x_1:N-1, y_1:N, \Theta) = \prod_{k=1}^{N}{ p(x_k | x_1:k-1, y_1:N, Theta) } $ hilft hier nicht weiter… %Ist forward filter backward simulation diesbezgl. wirklich anders??? %• Sec. 4: Forward filter backward simulation (sampling), schrittweise Herleitung (Methode von Carter/Kohn 1994) %=> rekursives Ziehen von Samples des Zustandes in umgekehrter Richtung der Zeitreihe %Notwendige Informationen: 1) Auswertung der Zustandsübergangs-Dichte; 2) Filter-Lösung in Vorwärtsrichtung; y und R kommen nicht vor! %Square-Root-Implementierung (Cholesky); Transformation in ein System ohne Korrelation zwischen Prozess- und Messrauschen %• Sec. 5: Sampling der Parameter %(p. 2: Mögliche a-Priori-Verteilungen der Parameter bei Gibbs-Sampling eingeschränkt gegenüber MH (Prior bei MH frei); Begründung dort: bei HMM blocked Gibbs Sampling-Schritt aus p(theta | X, Y) nur mit conjugate prior möglich) %Prior für die Parameter: Matrix-Normal-Inverse-Wishart (MNIW) [West1997] = Conjugate Prior für Gaussian State-space Model, d.h. Parameter-Conditional ist auch MNIW-verteilt %Genauer: Systemparameter (A, B, C, D) haben Matrix-Normal-Verteilung; Kovarianzen haben Inverse-Wishart-Verteilung (32) %Achtung: Dichte für Sampling der Parameter nutzt Messungen und Kovarianzen (Q, R) -> muss für Robustheit angepasst werden ?!? %Möglicher Ausweg wird angedeutet für den Fall nicht-konjugierter Prior: Nutzung eines Metropolis-Hastings-Samplers mit analytischer Dichte als Proposal Density %Initialisierung der Parameter mit EM %------------------ % zu einzelnen Methoden konkret: Direkte Optimierung %------------------ %"The minimum of the energy function can be computed by using various gradient-free or gradient based general optimization algorithms (see, e.g., Luenberger and Ye, 2008)" %"It is possible to find the derivatives in basically two ways (see, e.g., [Cappe2005])": %• "By formally differentiating the energy function recursion equations for a particular method." (Sensitivity Equations) %siehe [Särkkä2013, p. 212f] Theorem A.2: Partielle Ableitung der "Energy Function" direkt über Rekursion, schöne Herleitung %Quelle dort: Gupta, N. and Mehra, R. 1974. Computational aspects of maximum likelihood estimation and reduction in sensitivity function calculations. IEEE Transactions on Automatic Control, 19(6), 774–783. %im Fließtext auch [Cappe2005] %Auswertung erfolgt gemeinsam mit Filter-Rekursion, kein Smoother notwendig %=> für robustes Filter herleiten! %=> für Horizont herleiten! %• Using Fisher’s identity which expresses the gradient of the energy function as an expectation of the derivative of the complete data log likelihood over the smoothing distribution. %-> konkret in [Cappe2005b] (Score = first derivative of log-likelihood; Information = "Opposite" (Inverse?) of the second derivative of log-likelihood) %auch nochmals erwähnt in Bezug auf Ableitungsberechnung als Nebenprodukt von EM mittels Fishers Identity [Särkkä2013, p. 184] %[Särkkä2013, p. 213f] Theorem A.3, Partielle Ableitung der Energy Function über Fisher's identity %Kommentare zu Umweg Ableitungsberechnung über Q aus EM in [Särkkä2013, p. 185]: %* Useful because it is often easier to compute than the direct derivative of $p(y_{1:T} | \Theta)$ $\rightarrow$ esp. for linear models %* Not so useful in non-linear models because- gradient of energy function approximation via Gaussian filter- gradient via Fishers Identity on Q using Gaussian smoothing are often different (nochmals erwähnt auf p. 195) %* With Particle Filters, Fishers Identity is feasible to approximate the gradient of the energy function % % %[Shumway2011, p. 335ff] Knackige Beschreibung State Space Maximum Likelihood Estimation %"The innovations form of the likelihood function, which was first given by Schweppe (1965) ..." %... %"… is a highly nonlinear and complicated function of the unknown parameters." %The usual procedure is to fix x0 and then develop a set of recursions for the log likelihood function and its first two derivatives (for example, Gupta and Mehra, 1974). Then, a Newton-Raphson algorithm (see %Example 3.29) can be used successively to update the parameter values until the negative of the log likelihood is minimized." %[Shumway2011, p. "The standard errors are a byproduct of the estimation procedure, and we will discuss their evaluation later in this section, after Property 6.4." % %[Harvey2006, p. 368]: Prediction Error Decomposition -> numerical optimization % %[Durbin2012, p. 177f] Numerical maximization algorithms (for log likelihood with respect to unknown parameters) -> Newton's method; BFGS konkret mit approximate Hessian % %[Barber2013, Ch. 11.6 p. 269 Optimising the Likelihood by Gradient Methods]: Abgrenzung zu EM; Umformung für Ableitungsberechnung für Latent Variable Models (über Fisher's Identity, nicht erwähnt); Problem vermutlich: es muss über alle Zustände integriert werden, schwierig bei langen Zeitreihen? Wo ist der Unterschied zur Integration bei EM (vgl. p. 508)?? -> siehe Kommentar letzter Absatz vor 11.6.1 [Barber2016, p. 269] % %Direkte Optimierung der Likelihood mit gradientenbasierten Verfahren: [Murphy2013, p. 349]; Komplexität steigt bei Nebenbedingungen, z.B. positive Definitheit von Kovarianz-Matrizen % %[Schon2011b, Sec. 3] %Problem bei direkter Optimierung mit unterlagertem Particle Filter: Ableitungsberechnung der Particle-Lösung nach den Parametern %Direkte Optimierung mit Particle Filter / Smoother: problematisch laut [Särkkä2013, p. 197], da Lösung des Particle Filters zufällig und unstetig in Theta %lt. [Särkkä2013, p. 199] Direkte Optimierung mit Particle-MCMC mit Ableitungsberechnung über Fisher's Identity in [Ninness2010] Bayesian system identification via Markov chain Monte Carlo techniques %Mit Sample-basiertem Zustandsschätzer/Smoother (z.B. Particle Smoother) nur Stochastische Optimierung (vgl. [Kucukelbir2016] ADVI Stan) wegen Varianz der Zustandsschätzung? % %Optimierungsverfahren: %fmincon: interior point BFGS (default) -> BFGS erwähnt in [Harvey1990]; Skript/Diss Endisch erwähnen % %[Bavdekar2011] Direct Optimization for EKF (sec. 2.3): %"In problem P1, the decision variables are the covariance matrices Q and R. If the decision variables are chosen as $Q^(1/2)$ and $R^(1/2)$, then the constraints for maintaining positive definiteness of the matrices can, in principle, be relaxed. However, to avoid illconditioning of the optimisation problem during the intermediate gradient steps, itmaystill be necesary to incorporate suitable upper and lower bounds on the elements of $Q^(1/2)$ and $R^(1/2)$" %War ähnliches auch bei [Särkkä2013] rel. weit hinten -> Ableitungsberechnung für QR-Zerlegung??? % %bzgl. GKF: \glqq{}Furthermore, the log-likelihood is non-smooth, which means that trust region methods (e.g. Gauss-Newton and Levenberg-Marquardt) are not directly applicable.\grqq [Agamennoni2015] -> fmincon = trust region method? %============================================================================== \subsection{Umsetzung der Parameterschätzung} %============================================================================== Für die Umsetzung der Parameterschätzung im Rahmen der \ac{EOL}-Prüfung sind der \ac{ML}- (bzw. \ac{MAP}-) und die approximativen Bayes-Schätzer am erfolgversprechendsten. Beide Ansätze können auf Basis der Prediction-Error-Decomposition die vorgestellten robusten Filteralgorithmen als Basis verwenden, um darauf aufbauend eine Parameterschätzung durchzuführen. Als wesentlichem Unterschied ist hervorzuheben, dass die \ac{ML}-Methode Punktschätzungen für die Modellparameter liefert, während eine (approximative) Bayes-Schätzung eine Wahrscheinlichkeitsverteilung der Modellparameter als Resultat hat. Die Ergebnisse der \ac{ML}-Schätzung können also direkt für die Online-Prüfung im Prüfstand mit einem festen Parametersatz $\est{\paramV}$ genutzt werden. Im Sinne einer schnellen Durchführung des Prüfvorgangs (vgl. \secref{sec:EOLSplit}) kommt daher die \ac{ML}-Methode mit direkter numerischer Optimierung der Parameter zum Einsatz. Eine Erweiterung zur \ac{MAP}-Schätzung ist bei Bedarf durch die Erweiterung der Kostenfunktion möglich. Das verwendete Filter und der Optimierungsalgorithmus können dabei unabhängig voneinander ausgetauscht werden, die Verbindung zwischen den beiden Elementen stellt die Kostenfunktion über die Auswertung der Likelihood her. Eine wichtige Forderung an ein lernfähiges Prüfsystem war in \secref{sec:TrainingiO} das Training mit \ac{iO}-Daten bei gleichzeitiger Robustheit gegen nicht bekannte \ac{niO}-Messungen. Um auch für die Parameterschätzung der Zeitreihenmodelle die nötige Robustheit gegen solche Ausreißer herzustellen, beruht die Schätzung wiederum auf den robusten Filtern aus \secref{sec:tsRobustFilter}. Im Training werden alle vorhandenen Messungen ohne vorheriges Entfernen potentieller \ac{niO}-Messungen genutzt. Bei Ausschluss der \ac{niO}-Messungen aus den nachfolgenden Trainingsvorgängen bestünde ansonsten die Gefahr, dass der Trainingsalgorithmus die Varianz der Rauschterme verringert, was in der Prüfung zu engeren Prüfgrenzen führt und nachfolgend wiederum mehr Messungen aus dem Training ausschließt. Es sind daher auch solche Messungen in das Training eingeschlossen, die zuvor vom Prüfsystem als \ac{niO} bewertet wurden. %Die State-Augmentation-Methode stellt sich hingegen als primär relevant für Anwendungen mit Online-Parameterschätzung dar. Die korrelations-, kovarianz- und subspace-basierten Methoden %Begründung Auswahl Maximum-Likelihood (schnell, flexibel bzgl. verwendetem Filter, 1 Parametersatz für Prüfstände) , evtl. Info Gelman (?) suchen dass Hyperparameter mit ML ok weil viele Daten im Gegensatz zu latenten Variablen %vs. State-Augmentation: erscheint mit Filter bei vor allem bei kurzen Walks problematisch; Integration mit robusten Filtern? %vs. EM: [Barber2013, Ch. 24.5.2] Kapitel zu Learning Linear Dynamical Systems -> "The statistics required therefore include smoothed means, covariances, and cross moments" %----------------------------------------------------------------------------- \subsubsection{Numerische Optimierung} %----------------------------------------------------------------------------- \begin{figure} \centering \includegraphics{schema_numerische_ML} \caption[FIXME]{Schema der direkten numerischen Optimierung für die \ac{ML}-Parameterschätzung} \label{fig:tsParamSchematicOptimization} \end{figure} \Figref{fig:tsParamSchematicOptimization} zeigt ein Schema des numerischen Optimierungsvorgangs zur Parameterschätzung. Die Schätzung startet mit initialen Modellparametern $\paramV^{(1)}$. Eines der robusten Filter berechnet die Zustandsschätzung auf dem Trainingsdatensatz mit den aktuellen Modellparametern $\paramV^{(\iteration)}$. Für die prädiktiven Dichtefunktionen $\pdfBr{\given{\measV_\indexS}{\measV_{1:\indexS-1}, \paramV}}$ sind wie bei der Berechnung der Prüfgrenzen Näherungen $\pdfBr{\given{\measV_\indexS}{\measV_{1:\indexS-1}, \paramV}} \approx \pdfqSubBr{\learn}{\given{\measV_\indexS}{\measV_{1:\indexS-1}, \paramV}}$ erforderlich, die nachfolgend vorgestellt werden. Basierend auf den Prädiktionen des Filters erfolgt dann die Auswertung der Kostenfunktion \eqref{eq:tsParamDirectMLCostFun}. Ein numerischer Optimierungsalgorithmus bestimmt daraus den nächsten Schätzwert $\paramV^{(\iteration+1)}$. Mit diesem neuen Schätzwert startet die nächste Iteration, bis die Optimierung konvergiert. Für die Parameterschätzung kommen zahlreiche numerische Optimierungsalgorithmen für nichtlineare Kostenfunktionen mit mehreren reellen Parametern in Frage. Beispiele finden sich in der allgemeinen Literatur zur Optimierung \cite{Nocedal2006, Papageorgiou2012} und in Arbeiten zur Identifikation nichtlinearer Systeme \cite{Endisch2009}. Bei \figref{fig:tsParamSchematicOptimization} handelt es sich um ein vereinfachtes Schema, je nach Optimierungsalgorithmus kann pro Iteration auch einemehrfache Auswertung der Kostenfunktion notwendig sein. Bei den ableitungsbasierten Optimierungsverfahren sind außerdem der Gradient oder sogar die Hessematrix der Kostenfunktion zu bestimmen. Die Ableitungsberechnung gestaltet sich bei den Zeitreihenmodellen aufwendig, da das Ergebnis der Zustandsschätzung für das Sample $\indexS$ über den Prädiktions- und Update-Algorithmus des eingesetzten Filters von den Zustandsschätzungen aller vorangegangenen Zeitschritte abhängt. Es ist also für die Ableitungsberechnung ein rekursives Nachdifferenzieren von jedem Term in \eqref{eq:tsParamDirectMLCostFun} rückwärts bis zur Filterinitialisierung nötig. Dieser Vorgang ist vergleichbar mit dem Backpropagation-Schritt bei der Ableitungsberechnung neuronaler Netze. Für ein linear-gaußsches Zustandsraummodell mit einem Kalman-Filter als Zustandsschätzer findet sich die Herleitung der Rekursion bei Särkkä \cite[][212-214]{Sarkka2013}. Das \ac{GKF} hingegen weist aufgrund des schaltenden Verhaltens des Kalman-Gains $\kalmM_\indexS$ Unstetigkeiten in der Kostenfunktion auf. Diese Unstetigkeiten treten immer dann auf, wenn eine Messung $\measV_\indexS$ auf der Grenze ihres jeweiligen -- von den Modellparametern abhängigen -- Validation-Gate liegt. Ableitungsbasierte Optimierungsverfahren sind daher in Kombination mit dem \ac{GKF} fehleranfällig. Eine analytische Ableitungsberechnung ist außer für das Kalman-Filter aber auch für das \ac{TKF} und \ac{RME} möglich, da beide Filter einen analytischen Update-Schritt ohne Unstetigkeiten oder innere Iterationen aufweisen. Im Gegensatz dazu führt das \ac{ORKF} innerhalb jedes Update-Schrittes eine vorab nicht bekannte Anzahl an Iterationen des Update-Algorithmus aus, was eine analytische Ableitungsberechnung erschwert. Eine Festlegung der Anzahl der Iterationen je Update-Schritt ermöglicht auch für das \ac{ORKF} eine analytisch Ableitungsberechnung. Als Alternative zur analytischen Ableitungsberechnung stehen auch Optimierungsverfahren zur Verfügung, die numerisch Näherungen für die benötigten Ableitungen berechnen. Weit verbreitet sind Quasi-Newton-Verfahren mit Näherung der Hesse-Matrix durch den Broyden–Fletcher–Goldfarb–Shanno-Algorithmus \cite[][135-144]{Nocedal2006} \cite[][2-32, 3-27 - 3-29, 6-5 - 6-10]{TMW2018}. Ein Versuch auf einem Random-Walk mit \ac{RME} als Filter hat ergeben, dass ein Quasi-Newton-Verfahren mit numerischer Ableitungsberechnung sogar schneller arbeitet als mit analytischer Berechnung des Gradienten. Der Grund dafür ist die relativ aufwendige analytische Ableitungsberechnung durch eine eigene Rekursion zusätzlich zum Filtervorgang. Bezüglich der praktischen Umsetzung der Optimierung sei angemerkt, dass \eqref{eq:tsParamDirectML} die gleiche Lösung $\est{\paramV}_{\text{\acs{ML}}}$ wie das äquivalente Optimierungsproblem \begin{IEEEeqnarray*}{rCl} \est{\paramV}_{\text{\acs{ML}}} &=& \arg \min_{\paramV} \costFunBr{\paramV} , \\ \costFunBr{\paramV} &=& - \log \pdfBr{\given{ \measV_{1:\numS} }{ \paramV }} = - \log \pdfBr{\given{ \measV_1 }{ \paramV }} - \sum_{\indexS=2}^{\numS}{\log \pdfBr{\given{ \measV_\indexS }{ \measV_{1:\indexS-1}, \paramV }}} \end{IEEEeqnarray*} hat. Die logarithmierten Wahrscheinlichkeitsdichtefunktionen sind numerisch besser handhabbar und kommen deshalb bei der praktischen Umsetzung zur Anwendung. %%============================================================================== %\subsection{Ableitungsberechnung für RKF} %%============================================================================== % %Ergebnisse direkte Optimierung siehe KalmanLikelihoodCostBatch Test.m %z.B. Plot Kostenfunktion, Gradient mit Quiver, Vergleich Laufzeit ohne/mit analytischer Gradientenberechnung, evtl. Verlauf Schritte Optimierung vgl. Diss CE %außerdem kompliziert für Erweiterungen Horizont, Begrenzung der Kovarianz im Filter, Schätzung Schiefe etc. => keine Gradientenberechnungen nachfolgend %Plot-Einstellungen siehe Commit ed5858ffec2fdc3b5a115c93b4df86bcc9ec079a %----------------------------------------------------------------------------- \subsubsection{Prädiktive Dichtefunktion für das \ac{GKF}} %----------------------------------------------------------------------------- Für die Auswertung der Kostenfunktion sind die prädiktiven Verteilungen aus \eqref{eq:tsParamMarkov} zu berechnen. Dabei muss aufgrund der unterschiedlichen Verteilungsannahmen wieder zwischen dem \ac{GKF} mit normalverteiltem Messrauschen und dem \ac{TKF}, \ac{RME} und \ac{ORKF} mit Student-t-Messrauschen unterschieden werden. Für das \ac{GKF} gilt \eqref{eq:tsGKFPredictMeasurementAgainResult} bei gegebenem $\paramV$, resultierend in einer Normalverteilung für $\pdfBr{\given{ \measV_\indexS }{ \measV_{1:\indexS-1}, \paramV }}$. Die Kostenfunktion auf Basis der Normalverteilung ist analog zu den in \secref{sec:basics_robust} beschriebenen Zusammenhängen empfindlich gegen Ausreißer, d.h. $\costFunSubBr{\text{\acs{ML}}}{\paramV}$ reagiert stark Ausreißer in $\measV_\indexS$ und verfälscht entsprechend die Parameterschätzung. Um eine robuste Parameterschätzung zu erreichen kann wie bei der Zustandsschätzung das Validation-Gate herangezogen werden. Eine ersatzweise Dichtefunktion die alle $\measV_\indexS$ außerhalb des Validation-Gate aus $\costFunSubBr{\text{\acs{ML}}}{\paramV}$ ausschließt ist \begin{IEEEeqnarray*}{rCl} \pdfqSubBr{\learn}{\given{\measV_\indexS}{\measV_{1:\indexS-1}, \paramV}} &=& \begin{cases} \pdfBr{ \given{\measV_{\indexS}} {\measV_{1:\indexS-1}, \paramV} }, & \text{wenn } \SNPE(\resV_\indexS) \leq \SNPE_\gate \\ 1 , & \text{sonst} . \end{cases} \end{IEEEeqnarray*} Es handelt sich um das gleiche schaltende Verhalten wie bei der Berechnung des Update-Gains \eqref{eq:tsGKFGainCases}. %----------------------------------------------------------------------------- \subsubsection{Prädiktive Dichtefunktion für \ac{TKF}, \ac{RME} und \ac{ORKF}} %----------------------------------------------------------------------------- Für die drei auf Student-t-Messrauschen basierenden Filter wurde die prädiktive Dichte in \eqref{eq:tsStudentPredictionMarginalResult} allgemein aufgestellt, eine analytische Auswertung scheitert an der Kombination aus Normal- und Student-t-Verteilung. Die Näherung mit einer Normalverteilung aus \eqref{eq:tsStudentTestPredictionApprox} ist für die Parameterschätzung aufgrund der Ausreißer-Empfindlichkeit nicht geeignet. Daher ist für die Parameterschätzung eine Näherung der Zustandsschätzung durch eine Student-t-Verteilung $\pdfqSubBr{\learn}{\given{\stateV_{\indexS}}{\measV_{1:\indexS-1}, \paramV}}$ der sinnvollere Weg. Es resultiert daraus die prädiktive Dichte $\pdfqSubBr{\learn}{\given{\measV_\indexS}{\measV_{1:\indexS-1}, \paramV}}$ mit \begin{IEEEeqnarray*}{rCl} \pdfqSubBr{\learn}{\given{\measV_\indexS}{\measV_{1:\indexS-1}, \paramV}} &=& \int \pdfBr{\given{\measV_{\indexS}}{\stateV_{\indexS}, \paramV}} \: \pdfqSubBr{\learn}{\given{\stateV_{\indexS}}{\measV_{1:\indexS-1}, \paramV}} \, \dif \stateV_{\indexS}, \\ \pdfqSubBr{\learn}{\given{\stateV_{\indexS}}{\measV_{1:\indexS-1}, \paramV}} &=& \pdftBr{\given{\stateV_{\indexS}}{ \est{\stateV}_{\given{\indexS}{\indexS-1}}, \scaleDOF_{\dofMeasN,\infty}^2 \cdot \stateCM_{\given{\indexS}{\indexS-1}}, \dofMeasN }} , \\ \Rightarrow \pdfqSubBr{\learn}{\given{\measV_\indexS}{\measV_{1:\indexS-1}, \paramV}} &=& \pdftBr{\given{\measV_{\indexS}}{ \est{\measV}_{\given{\indexS}{\indexS-1}}, \measCM_{\indexS,\learn}, \dofMeasN }}, \\ \apr{\measCM}_{\indexS,\learn} &=& \scaleDOF_{\dofMeasN,\infty}^2 \cdot \outputM \stateCM_{\given{\indexS}{\indexS-1}} \trans{\outputM} + \measNCM. \end{IEEEeqnarray*} Zum Vorgehen bei der Approximation einer Normalverteilung durch eine Student-t-Verteilung und den eingesetzten Skalierungsfaktor $\scaleDOF_{\dofMeasN,\infty}$ siehe \secref{sec:basicsMathStoch}. %----------------------------------------------------------------------------- \subsubsection{Festlegung der Zahl der Freiheitsgrade} %----------------------------------------------------------------------------- Für das Zustandsraummodell mit Student-t-Messrauschen stellt die Anzahl der Freiheitsgrade $\dofMeasN$ einen weiteren Modellparameter dar. Grundsätzlich kann $\dofMeasN$ genauso wie die anderen Parameter anhand eines Trainingsdatensatzes geschätzt werden. Enthält nun aber der Trainingsdatensatz keine besonders auffälligen Messungen, dann weist die resultierende Schätzung von $\dofMeasN$ möglicherweise einen hohen Wert auf, da die Robustheitseigenschaft für den Trainingsdatensatz ja nicht erforderlich war. Bei große Werte von $\dofMeasN$ gehen das \ac{TKF}, der \ac{RME} und das \ac{ORKF} aber in das nicht-robuste Verhalten des Kalman-Filters über. Eine damit ausgeführte Prüfung wäre wiederum nicht robust gegen Ausreißer, sodass $\dofMeasN$ für die Prüfung nach oben zu begrenzen wäre. Eine Schätzung von $\dofMeasN$ ist deshalb nur eingeschränkt praxistauglich. Da die Student-t-Verteilung zur Herstellung einer Ausreißer-Robustheit in das Zustandsraummodell eingeführt wurde, kann $\dofMeasN$ auch direkt manuell festgelegt werden. Werte von $\dofMeasN < 3$ sind dabei unrealistisch \cite[][437]{Gelman2014}. Kleine Werte von $\dofMeasN$ sorgen für eine gute Robustheit der Filter gegen Ausreißer. Große Werte von $\dofMeasN$ sorgen hingegen dafür, dass Abweichungen der Zustandsschätzung nach wenigen Messungen ausgeglichen werden, da große Residuen in den Filter-Updates stärker berücksichtigt werden. Lange et al. \cite{Lange1989} berichten von guten Erfahrungen mit $\dofMeasN = 4$ in Regressionsmodellen. Im Zuge der Untersuchungen für diese Arbeit hat sich ein Wert von $\dofMeasN = \num{20}$ bewährt und kommt nachfolgend zum Einsatz. %Zusammenhang zwischen Häufigkeit der Ausreißer (Training/Test) und optimalem DOF finden? Durch Simulation? Durch Experiment mit sehr vielen Prüfläufen? (Achtung typisch: mehrere lokale Optima) %Beispiel für Umgang mit DOF, Problem lokaler Minima, Heuristik zum Einstellen: [Berger2012, p. 77] %siehe auch Unterschied Serienbetrieb (DOF ca. 20 um schnell genug auf Sprünge zu reagieren) vs. Versuchsaufbau (DOF ca. 3 um ausreichend robust bei geringer Datenmenge zu sein) %Ideen zu Sensitivity Analysis: [Gelman2014, Ch 6, Ch 17.1]; z.B. Student's t DOF verändern => Einfluss auf Schätzung beobachten => Einfluss einer Normalverteilungsannahme wird sichtbar %[Gelman2014, p. 437]: "In applications for which the t is chosen simply as a robust alternative to the normal, the degrees of freedom can be fixed at a small value to allow for outliers, but no smaller than prior understanding dictates. For example, t’s with one or two degrees of freedom have infinite variance and are not usually realistic in the far tails." %laut [Piche2012]: DOF 4 empfohlen in [Lange1989] %[Lange1989]: Lange Einleitung zu möglichen robusten Modellen und deren Darstellung als Mixture (mit latent variable) %DOF einstellen: ML, oder: "we have found that the value $\nu$ = 4 has worked well in many of our applications" (p. 883) %----------------------------------------------------------------------------- \subsubsection{Identifizierbarkeit und Nebenbedingungen} %----------------------------------------------------------------------------- Aus einem gegebenen Zustandsraummodell lassen sich durch Transformation des Zustandsvektors im Allgemeinen beliebig viele neue Modelle erzeugen. Diese erzeugen exakt die gleichen Messungen $\measV_\indexS$ wie das ursprüngliche Modell, eine Unterscheidung anhand der Messungen ist also nicht möglich. Für die Parameterschätzung bedeutet das, dass unter Umständen beliebig viele gleichwertige Lösungen für die Modellparameter existieren \cite[][Abschnitt 24.5.1]{Barber2013}. Einschränkungen können zum Beispiel für die Eigenwerte der Systemmatrix $\stateM$ vorgegeben werden, um die Stabilität des Modells zu garantieren \cite[][646]{Murphy2013}. Die mangelnde eindeutige Identifizierbarkeit der Modellparameter ist für die \ac{EOL}-Prüfung insofern unkritisch, als dass eine direkte Interpretierbarkeit des Zustandsvektors wie z.B. in einem Regelungssystem nicht angestrebt wird. Im Vordergrund steht die Bewertung der Messungen der Merkmalswerte $\measV_\indexS$ gegenüber ihren prädizierten Werten $\measV_{\given{ \indexS }{ \indexS-1 }}$, welche die transformierten Modelle wiederum gleichwertig wiedergeben. Random-Walk-Modelle sind diesbezüglich besonders gut handzuhaben, da die feste Vorgabe von $\stateM = \outputM = \eye$ die Mehrdeutigkeiten beseitigt. Bezüglich der Kovarianz- bzw. Skalierungs-Matrizen $\procNCM$ und $\measNCM$ ist zu garantieren, dass deren Schätzwerte symmetrisch und positiv semidefinit sind. Möglichkeiten hierfür sind einerseits Nebenbedingungen in der Parameterschätzung. Andererseits kann die auch eine Umparametrierung erfolgen, sodass zum Beispiel anstatt $\procNCM$ als Hilfsgröße die Cholesky-Zerlegung \begin{IEEEeqnarray*}{c} \procNCM = \trans{\UTM} \UTM \end{IEEEeqnarray*} in Form einer oberen Dreiecksmatrix $\UTM$ geschätzt wird \cite{Pinheiro1996}. Es reicht nun sicherzustellen dass die Diagonalelemente von $\UTM$ größer gleich Null sind, um zu erreichen dass $\procNCM$ positiv semidefinit wird. Nebenbedingungen zur Sicherstellung der Positivität lassen sich durch Optimierung der logarithmierten Parameter vermeiden \cite[][178]{Durbin2012}. % Transformation von Parameterbegrenzungen mit log, logit, probit: [Gelman2014, p. 22] %----------------------------------------------------------------------------- \subsubsection{Filterinitialisierung} %----------------------------------------------------------------------------- Sowohl für das Training als auch für die Prüfung benötigen die Filteralgorithmen eine Initialisierung der Verteilungsparameter $\est{\stateV}_{\given{ 1 }{ 0 }}$ und $\stateCM_{\given{ 1 }{ 0 }}$ der Zustandsschätzung. Beim Prüfvorgang gestaltet sich die Initialisierung einfach, da zwangsweise zuvor ein Trainingsvorgang für die Parameterschätzung stattgefunden hat. Dessen letzte Zustandsschätzung in Form der Parameter $\est{\stateV}_{\given{ \numS_\learn }{ \numS_\learn }}$ und $\stateCM_{\given{ \numS_\learn }{ \numS_\learn }}$ kann direkt in einen Prädiktions-Schritt übernommen werden, der den Zeitraum zwischen dem letzten Sample des Trainingsdatensatzes und dem nächsten Prüfvorgang nach Übernahme der neuen Modellparameter im Prüfstand überbrückt. Bezüglich der Filterinitialisierung in der Parameterschätzung ist eine solche Übernahme zuvor berechneter Zustandsschätzungen spätestens dann nicht mehr möglich, wenn das lernfähige Prüfsystem bei neuen Produkten oder Produktvarianten erstmals zur Anwendung kommt. Auch jenseits eines solchen Szenarios kann es vorteilhaft sein, eine Parameterschätzung \glqq{}von Null an\glqq{} nur anhand der vorhandenen Datensätze zu starten. Eine gängige Methode zur Initialisierung bei Kalman-Filtern ist, $\stateCM_{\given{ 1 }{ 0 }} = \constA \eye, \constA \rightarrow \infty$ zu setzen (diffuse Initialisierung) \cite[][12, 295]{Harvey2004}. Dadurch startet der Filter-Vorgang mit einem hohen Gain $\kalmM_1$, und der Einfluss des initialen Zustandes $\est{\stateV}_{\given{ \numS_\learn }{ \numS_\learn }}$ verschwindet. Dieser hohe Gain widerspricht allerdings dem Ziel der Robustheit gegen Ausreißer, da Ausreißer nur eine geringe Berücksichtigung im Update-Schritt erhalten sollen. Das langsame Driftverhalten der Prüfmerkmale in der \ac{EOL}-Prüfung erlaubt beim Random-Walk-Modell eine heuristische Initialisierung mit \begin{IEEEeqnarray*}{rCl} \stateV_{\given{ 1 }{ 0 }} &=& \trans{\mvbegin{\est{\meanSym}_{1}, \ldots, \est{\meanSym}_{\nStateV}}}, \\ \stateCM_{\given{ 1 }{ 0 }} &=& \diagBr{\est{\stdSym}_{1}^2, \ldots, \est{\stdSym}_{\nStateV}^2}. \end{IEEEeqnarray*} Als einfacher aber robuster Schätzer für $\est{\meanSym}_{\indexFeat}$ kann der Median von Merkmal $\indexFeat$ über einen Abschnitt am Beginn des Trainingsdatensatzes dienen. Nachfolgend werden die ersten zehn Samples dafür herangezogen. Für $\est{\stdSym}_{\indexFeat}$ wird die Stichproben-Standardabweichung des gesamten Trainingsdatensatzes eingesetzt. Um die nötige Robustheit gegen einzelne Ausreißer herzustellen werden hierbei \SI{5}{\percent} der extremsten Merkmalswerte von der Berechnung ausgeschlossen (vgl. \secref{sec:basicsRobust}). %============================================================================== \subsubsection{Begrenzung der Varianz in der Prädiktion} %============================================================================== Bei einer langen Pause $\Delta t_\indexS$ zwischen zwei Prüfläufen kann die Varianz der Zustandsprädiktion in Form der Kovarianzmatrix $\stateCM_{\given{ \indexS+1 }{ \indexS }}$ theoretisch beliebig weit ansteigen. Beim Random-Walk-Modell lässt sich der Anstieg der Kovarianz direkt mit $\stateCM_{\given{ \indexS+1 }{ \indexS }} = \stateCM_{\given{ \indexS }{ \indexS }} + \Delta t_\indexS \: \procNCM $ angeben. Auf Basis von $\stateCM_{\given{ \indexS+1 }{ \indexS }}$ öffnet sich auch die Prüfgrenze entsprechend weit. Außerdem liegt beim Update-Schritt nach der Pause eine Situation wie bei der diffusen Filterinitialisierung vor, sodass eine robuste Zustandsschätzung nicht mehr garantiert werden kann. Ein vergleichbarer Effekt ist bei Zustandsschätzern mit gleichmäßiger Abtastrate als Kovarianz-Windup bekannt. Er tritt auf wenn durch mangelnde Systemanregung die Beobachtbarkeit der Zustände zeitweise verloren geht \cite[][473-480]{Astrom1995}. Stenlund und Gustafsson geben verschiedene Möglichkeiten an, einen solchen Kovarianz-Anstieg zu begrenzen \cite{Stenlund2002}. Eine solche Möglichkeit ist, die Eigenwerte von $\stateCM_{\given{ \indexS+1 }{ \indexS }}$ zu begrenzen. Die Motivation dahinter ist, dass die Eigenwerte die Varianz in Richtung der Hauptachsen der Ellipsoide konstanter Wahrscheinlichkeitsdichte angeben \cite[][78-81]{Bishop2006}. Im Folgenden gilt die Annahme, dass plausible Prüfgrenzen die Gesamtvarianz der Daten im Trainingsdatensatz nicht dramatisch übersteigen sollen. Für das Random-Walk-Modell wird daher die Varianz der Zustandsprädiktion auf das Zweifache der Gesamtvarianz des jeweiligen Prüfmerkmales im Trainingsdatensatz begrenzt. Die zugrunde liegende Varianzschätzung kann aus der Filterinitialisierung übernommen werden. %============================================================================== \subsection{Parameterschätzung auf simulierten Random-Walks} %============================================================================== Wie bereits beim Vergleich der Filteralgorithmen mit festen Parametern in \secref ist auch in Verbindung mit der Parameterschätzung zunächst eine Untersuchung anhand simulierter Datensätze angebracht. Nur bei simulierten Datensätzen ist aufgrund der exakt bekannten Modellparameter ein Vergleich mit der optimalen Lösung möglich. Die Basis bilden \num{1000} simulierte Random-Walks, deren Modellparameter \tabref{tab:tsParamSimWalk} zeigt. Für jeden Random-Walk findet eine Parameterschätzung nach der \ac{ML}-Methode mit direkter numerischer Optimierung statt. \Figref{fig:tsParamGaussWalkScatter} zeigt die Ergebnisse der Parameterschätzung auf Basis der vier Filteralgorithmen für Random-Walks mit $\numS_\learn = \num{100}$ und $\numS_\learn = \num{1000}$ Samples im Trainingsdatensatz. \begin{table} \caption[FIXME]{Parameter für die Simulation der Random-Walks für die Parameterschätzung} \label{tab:tsParamSimWalk} \centering \footnotesize \begin{tabular}{ll} \toprule Parameter & Wert \\ \midrule $\procNC$ & \num{0.1} \\ $\measNC$ & \num{1} \\ $\state_1$ & \num{0} \\ $\Delta t_{\indexS}$ & gemäß geschätzter Verteilung aus \secref{sec:tsTimeDistUnequal} \\ $\stateG$ & \num{1} \\ $\inputG$ & \num{0} \\ $\outputG$ & \num{1} \\ $\feedG$ & \num{0} \\ \bottomrule \end{tabular} \end{table} % ------------------------------------------------------------------- % EC1122 without outliers, histograms %\begin{figure} %\centering %\footnotesize %\setlength\figurewidth{.40\columnwidth} %\setlength\figureheight{.23\columnwidth} %\renewcommand\figureXLabel{FIXME} %\renewcommand\figureYLabel{FIXME} %\subcaptionbox{\ac{GKF}}[.99\columnwidth]{% %\centering %\footnotesize %%\mytikzexton %\input{"graphics/P1122 01 01/011 Histogram Q.tikz"} %\input{"graphics/P1122 01 01/012 Histogram R.tikz"} %%\mytikzextoff %}% %% --------------- comment magic (tm) %\\% %\vspace{3mm}% %% --------------- comment magic (tm) %\subcaptionbox{\ac{TKF}}[.99\columnwidth]{% %\centering %\footnotesize %%\mytikzexton %\input{"graphics/P1122 04 01/011 Histogram Q.tikz"} %\input{"graphics/P1122 04 01/012 Histogram R.tikz"} %%\mytikzextoff %}% %% --------------- comment magic (tm) %\\% %\vspace{3mm}% %% --------------- comment magic (tm) %\subcaptionbox{\ac{RME}}[.99\columnwidth]{% %\centering %\footnotesize %%\mytikzexton %\input{"graphics/P1122 02 01/011 Histogram Q.tikz"} %\input{"graphics/P1122 02 01/012 Histogram R.tikz"} %%\mytikzextoff %}% %% --------------- comment magic (tm) %\\% %% --------------- comment magic (tm) %\subcaptionbox{\ac{ORKF}}[.99\columnwidth]{% %\centering %\footnotesize %%\mytikzexton %\input{"graphics/P1122 03 01/011 Histogram Q.tikz"} %\input{"graphics/P1122 03 01/012 Histogram R.tikz"} %%\mytikzextoff %}% %\caption[FIXME]{FIXME} %\label{fig:tsParamGaussWalkHistogram} %\end{figure} % ------------------------------------------------------------------- % EC1122 without outliers, scatter plots \begin{figure} \centering \footnotesize \setlength\figurewidth{.40\columnwidth} \setlength\figureheight{.23\columnwidth} \renewcommand\figureXLabel{$\sqrt{\est{\procNC} / \procNC}$} \renewcommand\figureYLabel{$\sqrt{\est{\measNC} / \measNC}$} % line 1: \subcaptionbox{\ac{GKF}, $\numS_\learn = \num{100}$}[.49\columnwidth]{% \centering \footnotesize \mytikzexton \input{"graphics/P1122 11 01/022 Scatter QR.tikz"} \mytikzextoff }% \subcaptionbox{\ac{GKF}, $\numS_\learn = \num{1000}$}[.49\columnwidth]{% \centering \footnotesize \mytikzexton \input{"graphics/P1122 01 01/022 Scatter QR.tikz"} \mytikzextoff \vspace{-0.7mm} % caution: this vspace acts weird, negative space moves the figure down }% % --------------- comment magic (tm) \\% \vspace{3mm}% % --------------- comment magic (tm) % line 2: \renewcommand\figureXLabel{$\sqrt{\est{\procNC} / \procNC}$}% \renewcommand\figureYLabel{$\sqrt{\scaleDOF_{\dofMeasN,\infty}^2 \est{\measNC} / \measNC}$}% \subcaptionbox{\ac{TKF}, $\numS_\learn = \num{100}$}[.49\columnwidth]{% \centering \footnotesize \mytikzexton \input{"graphics/P1122 14 01/022 Scatter QR.tikz"} \mytikzextoff }% \subcaptionbox{\ac{TKF}, $\numS_\learn = \num{1000}$}[.49\columnwidth]{% \centering \footnotesize \mytikzexton \input{"graphics/P1122 04 01/022 Scatter QR.tikz"} \mytikzextoff \vspace{-0.7mm} % caution: this vspace acts weird, negative space moves the figure down }% % --------------- comment magic (tm) \\% \vspace{3mm}% % --------------- comment magic (tm) % line 3: \subcaptionbox{\ac{RME}, $\numS_\learn = \num{100}$}[.49\columnwidth]{% \centering \footnotesize \mytikzexton \input{"graphics/P1122 12 01/022 Scatter QR.tikz"} \mytikzextoff }% % --------------- comment magic (tm) \subcaptionbox{\ac{RME}, $\numS_\learn = \num{1000}$}[.49\columnwidth]{% \centering \footnotesize \mytikzexton \input{"graphics/P1122 02 01/022 Scatter QR.tikz"} \mytikzextoff \vspace{-0.7mm} % caution: this vspace acts weird, negative space moves the figure down }% % --------------- comment magic (tm) \\% \vspace{3mm}% % --------------- comment magic (tm) % line 4: \subcaptionbox{\ac{ORKF}, $\numS_\learn = \num{100}$}[.49\columnwidth]{% \centering \footnotesize \mytikzexton \input{"graphics/P1122 13 01/022 Scatter QR.tikz"} \mytikzextoff }% % --------------- comment magic (tm) \subcaptionbox{\ac{ORKF}, $\numS_\learn = \num{1000}$}[.49\columnwidth]{% \centering \footnotesize \mytikzexton \input{"graphics/P1122 03 01/022 Scatter QR.tikz"} \mytikzextoff \vspace{-0.7mm} % caution: this vspace acts weird, negative space moves the figure down }% \caption[FIXME]{Streudiagramme der Parameterschätzung \circleWithParens{myLineOne} auf \num{1000} linear-gaußschen Random-Walks mit $\numS_\learn = \num{100}$ bzw. $\numS_\learn = \num{1000}$ Samples Länge. Die geschätzten Parameter sind bezogen auf ihren Sollwert dargestellt, der Punkt $\trans{\mvbegin{1, 1}}$ entspricht daher dem Sollwert \crossWithParens{myLineTwo}. Bei den Filtern \ac{TKF}, \ac{RME} und \ac{ORKF} wurde $\est{\measNC}$ mit $\scaleDOF_{\dofMeasN,\infty}^2$ skaliert, um auf das normalverteilte Messrauschen zurückzuschließen.} \label{fig:tsParamGaussWalkScatter} \end{figure} % ------------------------------------------------------------------- % EC1122 with outliers \begin{figure} \centering \footnotesize \setlength\figurewidth{.40\columnwidth} \setlength\figureheight{.23\columnwidth} \renewcommand\figureXLabel{$\sqrt{\est{\procNC} / \procNC}$} \renewcommand\figureYLabel{$\sqrt{\est{\measNC} / \measNC}$} % line 1: \subcaptionbox{\ac{GKF}, $\numS_\learn = \num{1000}$}[.49\columnwidth]{% \centering \footnotesize \mytikzexton \input{"graphics/P1122 21 01/022 Scatter QR.tikz"} \mytikzextoff }% \renewcommand\figureXLabel{$\sqrt{\est{\procNC} / \procNC}$}% \renewcommand\figureYLabel{$\sqrt{\scaleDOF_{\dofMeasN,\infty}^2 \est{\measNC} / \measNC}$}% \subcaptionbox{\ac{TKF}, $\numS_\learn = \num{1000}$}[.49\columnwidth]{% \centering \footnotesize \mytikzexton \input{"graphics/P1122 24 01/022 Scatter QR.tikz"} \mytikzextoff \vspace{-0.7mm} % caution: this vspace acts weird, negative space moves the figure down }% % --------------- comment magic (tm) \\% \vspace{3mm}% % --------------- comment magic (tm) % line 2: \subcaptionbox{\ac{RME}, $\numS_\learn = \num{1000}$}[.49\columnwidth]{% \centering \footnotesize \mytikzexton \input{"graphics/P1122 22 01/022 Scatter QR.tikz"} \mytikzextoff }% \subcaptionbox{\ac{ORKF}, $\numS_\learn = \num{1000}$}[.49\columnwidth]{% \centering \footnotesize \mytikzexton \input{"graphics/P1122 23 01/022 Scatter QR.tikz"} \mytikzextoff \vspace{-0.0mm} % caution: this vspace acts weird, negative space moves the figure down }% \caption[FIXME]{Streudiagramme der Parameterschätzung \circleWithParens{myLineOne} auf \num{1000} linear-gaußschen Random Walks mit $\numS_\learn = \num{100}$ bzw. $\numS_\learn = \num{1000}$ Samples Länge und \SI{0.5}{\percent} Ausreißern. Die geschätzten Parameter sind bezogen auf ihren Sollwert dargestellt, der Punkt $\trans{\mvbegin{1, 1}}$ entspricht daher dem Sollwert \crossWithParens{myLineTwo}. Bei den Filtern \ac{TKF}, \ac{RME} und \ac{ORKF} wurde $\est{\measNC}$ mit $\scaleDOF_{\dofMeasN,\infty}^2$ skaliert, um auf das normalverteilte Messrauschen zurückzuschließen.} \label{fig:tsParamGaussOutlierWalkScatter} \end{figure} \begin{table} \caption[FIXME]{Schätzwerte zu den linear-gaußschen Random Walks mit $\numS_\learn = \num{100}$ Samples ohne Ausreißer} \label{tab:tsParamGaussEstimationResultsShort} \centering \footnotesize \begin{tabular}{lllll} \toprule Wert & \ac{GKF} & \ac{TKF} & \ac{RME} & \ac{ORKF} \\ \midrule Mittelwert $\sqrt{\est{\procNC} / \procNC}$ & (\num{2.939}) & \num{0.889} & \num{0.917} & \num{0.897} \\ Standardabweichung $\sqrt{\est{\procNC} / \procNC}$ & (\num{12.905}) & \num{0.908} & \num{0.931} & \num{0.905} \\ Mittelwert $\sqrt{\est{\measNC} / \measNC}$ bzw. $\sqrt{\scaleDOF_{\dofMeasN,\infty}^2 \est{\measNC} / \measNC}$ & \num{0.9480} & \num{1.002} & \num{1.000} & \num{1.002} \\ Standardabweichung $\sqrt{\est{\measNC} / \measNC}$ bzw. $\sqrt{\scaleDOF_{\dofMeasN,\infty}^2 \est{\measNC} / \measNC}$ & \num{0.146} & \num{0.075} & \num{0.074} & \num{0.074} \\ \bottomrule \end{tabular} \end{table} \begin{table} \caption[FIXME]{Schätzwerte zu den linear-gaußschen Random Walks mit $\numS_\learn = \num{1000}$ Samples ohne Ausreißer} \label{tab:tsParamGaussEstimationResultsLong} \centering \footnotesize \begin{tabular}{lllll} \toprule Wert & \ac{GKF} & \ac{TKF} & \ac{RME} & \ac{ORKF} \\ \midrule Mittelwert $\sqrt{\est{\procNC} / \procNC}$ & \num{1.212} & \num{1.008} & \num{1.033} & \num{1.019} \\ Standardabweichung $\sqrt{\est{\procNC} / \procNC}$ & \num{4.418} & \num{0.216} & \num{0.220} & \num{0.216} \\ Mittelwert $\sqrt{\est{\measNC} / \measNC}$ bzw. $\sqrt{\scaleDOF_{\dofMeasN,\infty}^2 \est{\measNC} / \measNC}$ & \num{0.960} & \num{1.000} & \num{1.000} & \num{1.001} \\ Standardabweichung $\sqrt{\est{\measNC} / \measNC}$ bzw. $\sqrt{\scaleDOF_{\dofMeasN,\infty}^2 \est{\measNC} / \measNC}$ & \num{0.045} & \num{0.024} & \num{0.024} & \num{0.024} \\ \bottomrule \end{tabular} \end{table} \begin{table} \caption[FIXME]{Schätzwerte zu den linear-gaußschen Random Walks mit $\numS_\learn = \num{1000}$ Samples mit Ausreißern} \label{tab:tsParamGaussEstimationResultsLongOutlier} \centering \footnotesize \begin{tabular}{lllll} \toprule Wert & \ac{GKF} & \ac{TKF} & \ac{RME} & \ac{ORKF} \\ \midrule Mittelwert $\sqrt{\est{\procNC} / \procNC}$ & \num{2.119} & \num{0.956} & \num{1.067} & \num{1.052} \\ Standardabweichung $\sqrt{\est{\procNC} / \procNC}$ & \num{9.178} & \num{0.241} & \num{0.249} & \num{0.242} \\ Mittelwert $\sqrt{\est{\measNC} / \measNC}$ bzw. $\sqrt{\scaleDOF_{\dofMeasN,\infty}^2 \est{\measNC} / \measNC}$ & \num{0.972} & \num{1.054} & \num{1.051} & \num{1.051} \\ Standardabweichung $\sqrt{\est{\measNC} / \measNC}$ bzw. $\sqrt{\scaleDOF_{\dofMeasN,\infty}^2 \est{\measNC} / \measNC}$ & \num{0.073} & \num{0.024} & \num{0.024} & \num{0.024} \\ \bottomrule \end{tabular} \end{table} % Results from EC1122: %01: %Runtime for training of 1000 walks: 37.3219 s, 0.037322 s per walk %02: %Runtime for training of 1000 walks: 77.314 s, 0.077314 s per walk %03: %Runtime for training of 1000 walks: 67.6297 s, 0.06763 s per walk %04: %Runtime for training of 1000 walks: 51.6531 s, 0.051653 s per walk %11: %Runtime for training of 1000 walks: 15.7998 s, 0.0158 s per walk %12: %Runtime for training of 1000 walks: 51.2948 s, 0.051295 s per walk %13: %Runtime for training of 1000 walks: 38.5142 s, 0.038514 s per walk %14: %Runtime for training of 1000 walks: 40.2283 s, 0.040228 s per walk %21: %Runtime for training of 1000 walks: 30.1791 s, 0.030179 s per walk %22: %Runtime for training of 1000 walks: 47.1014 s, 0.047101 s per walk %23: %Runtime for training of 1000 walks: 47.4448 s, 0.047445 s per walk %24: %Runtime for training of 1000 walks: 42.3813 s, 0.042381 s per walk %============================================================================== \subsection{Horizont für Robustheit gegen Modellverletzungen} %============================================================================== Ziel: gutmütiges Verhalten bei Verletzungen der Modellannahmen Industrielle Anwendung $\Rightarrow$ keine isolierte Umgebung => Robustheit gegen Verletzung der Modellannahmen nötig Achtung: Horizont nicht nur für Langzeit-Prädiktion, sondern z.B. auch für Merkmale mit mehreren Moden wichtig Gründe für Modellverletzungen: kleinere Sprünge z.B. durch Produktionslose / Tage; Merkmale mit mehreren Moden; Fokus auf Langzeitprädiktion Möglichkeiten siehe Auflistung Erfindungsmeldung Lösung 1: Modell erweitern, z.B. um autokorrelierte Rauschtermen (vgl. [BarShalom2001]) -> mehr Zustände pro Merkmal -> mehr Aufwand Modellierung und Berechnung steigt, Modellierung merkmalsabhängig je nach Autokorrelation? Lösung 2: Local linear trend model (auch mehr Zustände) (Teil des Modells für autokorrelierte Rauschterme?) Lösung 3: Multi-step Prediction, Horizont Umfangreiche Herleitung multi-step prediction siehe West1997 S. 106ff => wichtig für Training mit Horizont Als Aspekt der Robustheit, robust gegen Verletzungen der Modellannahmen... Einige Merkmale springen z.B. tageweise rel. stark, dazwischen geringere Varianz als Sprünge vermuten lassen; auch: mehrere Modi um den Prozess-mean, teilweise mit kurzen Aussetzern $\rightarrow$ Filter springt $\Rightarrow$ normalverteiltes Rauschen passt nicht, Q wird viel zu groß geschätzt Ähnlich: Merkmale bewegen sich langfristig (Wochen) weniger, als es die lokalen Veränderungen von Prüflauf zu Prüflauf vermuten lassen (wie Drift um einen langfristigen Mittelwert) $\Rightarrow$ wird vom Random-Walk-Modell nicht erfasst, Random Walk läuft beliebig weit davon (Varianz nicht beschränkt) ohne Horizont und für verschiedene Horizont-Längen zahlreiche generierte Walks schätzen -> Vergleich der Ergebnisse; Bias/Variance (vgl. Bootstrap, z.B. Diss Gaussian Mixture Regression Kapitel 7) Für Gating gemeinsam mit Horizont: siehe EstimateKalmanWindowGroupBatch.m; Schwierig: richtige Samples aus Likelihood-Berechnung ausschließen; sonst: Optimierung läuft nicht da Kostenfunktion durch Outlier kollabiert -> Vergleiche Verfahren mit/ohne Outlier in den Daten mit/ohne Robustheit bzgl. gelernter Prozessparameter Interessant: 05 Drift AdaptiveTracking KalmanLikelihoodCostBatch Test.m $\rightarrow$ bei gleichmäßigen Zeitabständen, ohne Messausreißer: Kostenfunktion hängt nicht vom Horizont (Länge, Anzahl Elemente) ab Commits 2f7239a80b76e2f7f58f81b6d184656e92a4ae4e und 3ba04b49d697b1cf71fa25f644a227085140e853 liefern brauchbare Plots Interessant: bei kleinerer DOF-Einstellung wird das Maximum der Likelihood ausgeprägter, die Schätzung liegt näher am realen Wert Wichtig: Darstellung richtig zuschneiden, schönes Beispiel: Commit a4fc8ed535fd00caece559474372dd16170b7fae Training Kalman-Filter mit Horizont: vgl. Cross-validation (maximum) likelihood (siehe Bachoc2013, VanDerLaan2004) Referenzen für MCMC +/- Robustheit +/- Switching: [FruhwirthSchnatter2006, Ch. 13.5 (p. 416)] Einfluss der DOF auf das Tracking von Sprüngen Kompromiss DOF niedrig für gute Parameter (Q bleibt klein) vs. DOF groß für schnelles Tracking von Sprüngen? %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Berücksichtigung von Sondereffekten} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %============================================================================== \subsection{Sprünge in Merkmalen berücksichtigen} %============================================================================== Gustafsson - Adaptive Filtering and Change Detection: CUSUM [Stenlund2002]: CUSUM für Erkennung von Unstetigkeiten Idee: Process Noise zu definierten Sprungzeitpunkten als Student's t mit niedrigen DOF freigeben (z.B. mit Filter Roth); Problem Outlier Confusion -> als nicht geeignete Lösung aufzeigen Idee: Eingang zu allen Zeitpunkten zur Optimierung gemeinsam mit Prozessparametern freigeben und a) häufige b) knapp aufeinanderfolgende Nutzung des Eingangs bestrafen; [Pearson2002, p. 62]: "A case in point is intervention analysis, proposed for modeling effects like labor strikes in economic time series [26, ch. 12]. There, the key assumption is that outliers occur at known times and the basic idea is to treat the anomalous event as a binary disturbance input taking the value when the event occurs and otherwise." [BarShalom2001, Kapitel 11.4.2] Maneuver Detection and Model Switching [BarShalom2001, Kapitel 11.6] The multiple model approach evtl interessant: interacting multiple models (IMM) Notizen zu Literatur (Adaptive Control) siehe Anfang Papier-Unterlagen Kalman Auflistung der Möglichkeiten siehe Erfindungsmeldung Dummy-Variablen für "structural breaks" siehe [Harvey2006, p. 354] und etwas detaillierter [Harvey2006, p. 399f] Vergleich Gating-KF und Robuster KF auf Zeitskala und Sample-Skala => Gating-KF schwingt nach Pause etwas schneller auf neue Grenze ein? (Abhängig von Sprunghöhe?) %============================================================================== \subsection{Robuste Filter mit Schiefe} %============================================================================== Interessant für Skewness: idWindow 213 233 319 325 (!) 342, siehe auch Liste aktivierter Fenster in EC1012 Problematisch auch bei aktivierter Skewness-Schätzung: idWindow 213 und 319 (skewness schwankt zeitlich extrem); 325 (immer wieder extreme Ausreißer) [Gelman2007, p. 548] Transformations; For additional reading on transformations, see Atkinson (1985), Mosteller and Tukey (1977), Box and Cox (1964), and Carroll and Ruppert (1981). auch [Gelman2014] Auflistung von 3 Möglichkeiten siehe Erfindungsmeldung! Beispiel für automatischen Umgang mit Schiefe: [Berger2012, Ch. 4.1] log, Auswertung Skewness um Transformation automatisch zu aktivieren Highest Density Interval (of predictive distribution) for skewed distributions: Kruschke DBDA 2E, um Seite 342f; Literatur zu Berechnung / Darstellung HDI auch [Hyndman1996] %============================================================================== \subsection{Zwangs-i.O.} %============================================================================== Auflistung Möglichkeiten siehe Erfindungsmeldung %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Ergebnisse beim Einsatz robuster Zustandsraummodelle in der \ac{EOL}-Prüfung} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Einfluss der DOF auf die gelernten Parameter ca. 3-5 repräsentative Einzelmerkmale auswählen, die alle relevanten (=später beschriebenen (!)) Effekte abdecken • Wie viele Datensamples werden für eine bestimmte Genauigkeit der Schätzung und Prädiktion benötigt? vgl [Kruschke2015, Ch. 11.5.1 (p. 330), Ch. 13]? Aufzeigen anhand simulierter Random-Walks • Welchen Einfluss haben die Einstellungen des Algorithmus (Horizont, Schiefe etc.) auf die Schätzung (Bias, Varianz, Anzahl erforderlicher Samples, …) und auf die Prädiktion? • Welchen Einfluss haben Modellverletzungen durch die Daten, z.B. lokale Sprünge, Schiefe etc. auf die Schätzung und Prädiktion? • Vergleich mit einem sehr einfachen Verfahren, z.B. var(y) oder cusum • i.O.-Walks erzeugen, n.i.O. einfügen, Erkennung auswerten • Vergleich mit manueller Grenze qualitativ • Für Auswertung robuster Modelle siehe auch [Gelman2014, Ch. 17 und 7.5] Einfluss der geschätzten Prozessparameter auf die Grenzen (a la Sensitivitätsanalyse) (evtl. erst im Abschnitt Berechnung Prüfgrenzen) Auswirkung von DOF auf die Prüfgrenzen darstellen Rechenzeit Training in Abhängigkeit von der Anzahl der Merkmale, Trainingssamples, Art des Filters, Korrektur Umgebungseinflüsse ja/nein Rechenzeit Prüfung Zusätzlich Fehlererkennungsrate angeben o.ä. Qualitative Ergebnisse: Analyse der Einblicke in die Daten und Ergebnisse, die durch die Verfahren gewonnen werden Einstellbarkeit der Prüfschärfe gut Auswirkung der Untergruppen Beispiel CWZ, DED idWindow=052, alle Prüfstände, Normalisierung ohne/mit Unterscheidung nach Motorvariante und ohne/mit Unterscheidung nach Prüfstand im RKF -> recht gut erkennbare Auswirkung auf die Prüfgrenze Demonstration mit Beispiel-Fehlern evtl. vom Versuchsaufbau %Idee 10/2017: Prädiktive Quantile mit Messdaten vergleichen -> stimmt z.B. das 1 prozent, 25 prozent, 50 prozent, 75 prozent, 99 prozent-Quantil? -> Das ist der interessante Wert um Prüfgrenzen festzulegen, und lässt sich anhand i.O. auswerten! %auch: "average log-likelihood of the test-set items" für verschiedene Verfahren und verschiedene Test-Datensätze (Motortypen, Merkmale) vergleichen (wie [Ruiz2017, p. 19f, Table 3]) -> "empirical performance"; macht nur für iO-Prüflinge Sinn! %=> Für Diss: log predictive likelihoods of held-out data / simulated future data für die quantitative Bewertung von robusten Filtern und robusten Parameterschätzern?! (vgl. Notizen [Gelman2014b]) %Plot: Q/Q oder direkt Rohdaten von Prädiktions-Wert über Rohdaten, vgl. residual plot? %Residuen-Plots verschiedener Merkmale (nun weniger für den Vergleich von Verfahren, sondern für den Vergleich von Merkmalen), evtl. log-Histogramm gegenüber erwarteter / gefitteter Dichte, Q/Q-Plot %Normal Probability Plot der Innovationen bzw. diff(y) => Normalverteilung prüfen %Autokorrelation Residuen plotten? (Achtung verschiedene Zeitabstände zwischen den Samples) %Bei leichten Modellverletzungen (z.B. Unstetigkeit) sollen nach wie vor Updates erfolgen -> Gewichtung (Student's t) besser als Gating; Vergleich zeigen Optimierung beschleunigen: • log(Q), log(R) schätzen • Transformation: (Q/R, R) schätzen; Parametergrenzen müssen auch transformiert werden -> deutliche Beschleunigung ca. Faktor 2 schneller siehe [Kucukelbir2016, p. 5] für Beispiele zur Transformation einer Variablen mit Support x>0 auf log(x), -Inf < log(x) < Inf Fisher Information of Reparametrization: en.wikipedia Fisher Information, section Reparametrization $=>$ könnte über Jakobi-Matrix der Transformation benutzt werden um Optimale Transformation herzuleiten? %============================================================================== \subsection{Genauigkeitsuntersuchung Parameterschätzung; Laplace und MCMC} %============================================================================== auch: Varianz der Parameterschätzung in Abhängigkeit von der Länge des Trainingsdatensatzes und Horizont: a) über grafische Darstellung/Auswertung der Kostenfunktion (= Energy Function in [Särkkä2013]) b) über MCMC-Ziehen aus Energy Function vgl. [Särkkä2013, p. 188] c) Variational Approximation of Parameter Posterior? alles für verschiedene Kombinationen und Größenordnungen von Q und R; schönere Darstellung bei Umparametrierung auf Q/R etc.?? Plot Likelihood (bzw. a-posteriori-Verteilung) vs. Laplace-Approximation vs. Metropolis-Hastings-Samples Diskussion Nachteile Laplace (schlechte Näherung, Problem an den Parameter-Grenzen) Marginale Verteilung der Zustände bei unbekannten Parametern => für Prüfgrenzen! siehe [Särkkä2013, p. 15, (1.10)] => z.B. MCMC [Särkkä2013, p. 178, p. 184, p. 189] "It is also possible to use a Laplace approximation (Gelman et al., 2004) which uses the second derivative (Hessian) of the energy function to form a Gaussian approximation to the posterior distribution" "However, to implement the Laplace approximation, we need to have a method to compute (or approximate) the second order derivatives of the energy function." => Option: Näherung Hesse durch Gradient vgl. Optimierungsmethoden? D.h. Näherung um die fertige Parameter-MAP-Lösung herum mit Hesse-Matrix -> Gauß-Verteilung Optimierungsmethoden? D.h. Näherung um die fertige Parameter-MAP-Lösung herum mit Hesse-Matrix -> Gauß-Verteilung Notizen zu Laplace-Approximation über Hesse-Matrix siehe Notizen Parameterschätzung II nach Ableitungsberechnung => damit quantitative Angabe der Parameterunsicherheit direkt aus Trainingsergebnis ?! Achtung: fmincon-Hesse-Matrix betrachtet Constraints mit (Lagrange-Multiplikatoren, siehe Doku fmincon) [Ninness2010, p. 40]: Bayes-Ansatz speziell für Systemidentifikation auf kurzen Datensätze sinnvoll, um Unsicherheit zu quantifizieren Üblicher Weg (p. 41 rechte Spalte bis p. 42 oben): Normalverteilungsannahmen, Linearisierungen, Annahme von Unabhängigkeit der Parameter… problematisch bei geringer Datenmenge [Ninness2013, p. 42]: Parameter Prior "sufficiently regular (for example, smooth)" -> gradient based search can be applied Aber: posterior der Parameter für Abschätzung der Schätzgenauigkeit wird in der Literatur nicht regelmäßig genutzt Für Beispiel Metropolis-Hastings mit RKF siehe KalmanLikelihoodCostBatch Test.m, Parameter-Grid hineingezoomed Commit 7999254134b5a1a39b3c3481dabd137b52decb6f Algorithmus mit simulierten Daten zerlegen, z.B. geschätztes Q über realem Q parametriert mit realem R etc, evtl inkl Konfidenzband aus Bootstrap Beispiel für den Vergleich MCMC mit einer analytischen Näherungslösung (VI): [Kucukelbir2016, Sec. 4] Nutzung der Hesse-Matrix der (log-) Likelihood (geschätzt durch fmincon), um die Unsicherheit der Parameterschätzung abzuschätzen -> vgl. Fisher Identity? Zusammenhang zur Laplace-Approximation [Beal2003, p. 34] und Konsequenzen / Einschränkungen daraus, z.B. rel. ungenau bei geringer Anzahl an Datenpunkten! Vergleich der Unsicherheit aus Laplace-Approximation mit der Unsicherheit in MC-Experimenten (viele generierte Walks schätzen lassen und Histogramm plotten) -> evtl. kann die Laplace-Approximation als Abschätzung der Unsicherheit verwendet werden; siehe Experimente in KalmanLikelihoodCostTest, Abweichungen Laplace-Approximation zur MCMC erheblich [Shumway2011, p. 337] Code: SE = sqrt(diag(solve(est hessian))) $->$ standard error of parameter estimates, extracted from BFGS optimization using KF likelihood [Shumway2011, p. 344] Asymptotic Distribution of the Estimators liefert eine Verteilung für den Parameter-Schätzfehler für Gaussche SS-Modelle über "asymptotic information matrix" I mit Nutzung der zweiten Ableitung der Likelihood (Formel) "For a Newton procedure, the Hessian matrix at the time of convergence can be used as an estimate of $n*I(\Theta_0)$ to obtain estimates of standard errors. Problem 1: Näherung stimmt nicht für robuste Modelle Problem 2: Näherung stimmt nicht für endliche Anzahl Samples, doch genau das ist der interessante Fall (vgl. Laplace-Approximation, das ist das gleiche). Problem 3: keine Aussage nahe der Parametergrenzen (Matlab-BFGS: Hessian enthält Lagrange-Mult.; Laplace-Approximation macht am Rand keinen Sinn) [Shumway2011, p. 359] "Several researchers have found evidence that samples must be fairly large before asymptotic results are applicable (Dent and Min, 1978; Ansley and Newbold, 1980). Moreover, as we discussed in Example 3.35, problems occur if the parameters are near the boundary of the parameter space." [Shumway2011, p. 359ff] Bootstrapping State-Space Models -> Schätzung Unsicherheit Parameterschätzung für ML-Parameteroptimierung; nach Durchsicht nicht ganz klar, für Verständnis siehe Code (p. 363!) oder weitere Quellen (Stoffer and Wall (1991), Stoffer and Wall (2004)) Code sieht aus als ob für jede Bootstrap-Ziehung eine neue Parameterschätzung nötig ist (Zeile 25) $->$ gleich zu Blocked Gibbs, MH o.ä. wechseln... "The bootstrap, however, allows us to investigate the small sample distribution of the estimators and, hence, provides more insight into the data analysis." [Lange1989, p. 888]: "There is a variety of methods for estimating standard errors after fitting a t-family model. The observed or expected information matrix could be used, or a bootstrap resampling scheme could be employed."; "In our implementation the observed information matrix is obtained by numerical differentiation of the score vector and matrix inversion of the resulting Hessian." %============================================================================== \subsection{Prüfgrenzen mit Parameterunsicherheit} %============================================================================== Bank aus RKF Gaussian Sum Filter: [Alspach1972] Diskussion Änderung in Abhängigkeit von der Anzahl der Trainings-Samples, weitere Analysen je nach Ergebnis Diskussion Alternativen, z.B. Gibbs-Sampler mit forward filtering backward sampling %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Zusammenfassung} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Vorteile: Aufwand steigt linear mit der Anzahl der Datenpunkte Stärken und Schwächen der verwendeten Methodik Ausblick auf Verbesserungsmöglichkeiten