Vaje iz Fizikalnih merjenj, 13. teden v letu

2026/03/24

Teorija: Kalmanov filter za vektorske spremenljivke v zvezni sliki

Dinamiko sistema opisuje enačba

\[ \dot{\mathbf{x}} = A(t) \mathbf{x} + \mathbf{c}(t) + \Gamma \mathbf{w}(t), \]

kjer je \( \mathbf{w} \) vektor dinamičnega šuma. Za tako dinamiko lahko uporabimo Kalmanov filter za kovariančno matriko

\[ \mathbf{P} = AP + PA^T + \Gamma Q \Gamma^T - PH^TR^{-1} HP. \]

Zadnji člen ostri kovariančno matriko, \( H \) je okenska matrika, \( \mathbf{R} \) pa je kovariančna matrika meritev. V večini primerov bomo uporabljali samo

\[ \mathbf{P} = AP + PA^T + \Gamma Q \Gamma^T. \]

Naloga:

Kroglice z maso \( m \) spuščamo skozi dvojna svetlobna vrata. V prečni smeri med padanjem kroglic delujejo naključne sile \( F(t) \), za katere velja \( \left\langle F(t) \right\rangle = 0 \). Imamo dvoja svetlobna vrata enakih dimenzij. Vsaka vrata imajo polmer \( R \) in medsebojna razdalja je \( h \). Iz višine \( h \) nad prvimi svetlobnimi vrati spuščamo kroglice mase \( m \).

Skozi prva vrata pade \( \frac{2}{3} \) vseh kroglic. Kolikšen delež kroglic pade skozi druga vrata? Nedoločenost lege in hitrosti ob času \( t = 0 \) oziroma na izhodišču je \( \sigma_x (t = 0) = 0 \) in \( \sigma_v(t = 0) = 0. \)

Sile bomo opisali kot dinamični šum. Čas ob prehodu skozi prva vrata bomo označili s \( t_1 \), skozi druga vrata pa s \( t_2 \).

Trajektorije v prečni smeri ne moremo izračunati, saj na kroglice delujejo samo naključne sile.

S padanjem kroglice se bo širina Gaussove porazdelitve širila zaradi naključnih sil. Začetna Gaussova porazdelitev je zaradi nedoločenosti lege \( \sigma_x (t = 0) = 0 \) kar \( \delta \) funkcija.

Zanima na torej verjetnost

\[ P \left( - R < x(t_2) < R \right) = P \left( x(t_2) < R \right) - P \left( x(t_2) < -R \right). \]

Središče Gaussovih porazdelitev se ne bo premikalo. Premikalo bi se v primeru, če bi \( \left\langle F(t) \right\rangle = F_0 \). Če bi to držalo, bi naključne sile opisali kot \( F(t) = F_0 + f(t) \), kjer je \( \left\langle f(t) \right\rangle = 0\).

Verjetnost bomo zapisali s kumulativnimi funkcijami

\[ P \left( - R < x(t_2) < R \right) = F \left( \frac{R - 0}{\sigma_x (t_2)} \right) - F\left( \frac{- R + 0}{\sigma_x (t_2)} \right). \]

Upoštevamo lastnost kumulativnih funkcij \( F(-u) = 1- F(u) \). Verjetnost je tako

\[ P \left( - R < x(t_2) < R \right) = 2 F \left( \frac{R}{\sigma_x (t_2)} \right) - 1. \]

Dinamika - gibalna enačba

Poiščimo dinamiko našega sistema, izhajali pa bomo iz drugega Newtonovega zakona.

\[ \sum\limits_i^{} F_i = ma_x = m \ddot{x} = F(t). \]

Pospeške potem označimo kot

\[ \ddot{x} = \frac{F(t)}{m} = a(t) = w(t). \]

Pospešek v prečni smeri bo torej dinamični šum. Ker Kalmanov filter deluje samo za LDE 1. reda, moramo preiti na vektorsko sliko, kjer sta spremenljivki in enačbi \( \dot{x} = v \) in \( \dot{v} = \ddot{x} = w \). V vektorskem zapisu

\[ \mathbf{x} = \begin{bmatrix} \dot{x} \\ \dot{v} \end{bmatrix}. \]

Dinamiko lahko sedaj zapišemo kot

\[ \begin{bmatrix} \dot{x} \\ \dot{v} \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} x \\ v \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \end{bmatrix} + \begin{bmatrix} 0 \\ w \end{bmatrix}. \]

Kovariančna matrika bo imela obliko

\[ P = \begin{bmatrix} \sigma_{x} ^{2} & \sigma_{xv} \\ \sigma_{vx} & \sigma_{v} ^2 \end{bmatrix} = \begin{bmatrix} p_{11} & p_{12} \\ p_{21} & p_{22} \end{bmatrix}. \]

Mi pa iščemo z novo oznako \( \sigma_x (t) = \sqrt{p_{11}} \).

Kovariančna matrika

Uporabili bomo enačbo povedano v teoretičnem delu

\[ \dot{P} = AP + PA^T + \Gamma Q \Gamma^T. \]

Vrednosti elementov v matriki \( P \) ne poznamo in to niso začetne vrednosti,

\[ P = \begin{bmatrix} p_{11} & p_{12} \\ p_{21} & p_{22} \end{bmatrix}. \]

Ob času \( t = 0 \) je kovariančna matrika

\[ P(t = 0) = \begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix}. \]

Ni potrebe po računanju členov \( AP \) in \( PA^T \), saj velja, da je matrika simetrična \( \left( AP \right)^T = P^T A^T = PA^T \).

Torej zapišemo

\[ AP = \begin{bmatrix} 0 & 1 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} = \begin{bmatrix} p_{12} & p_{22} \\ 0 & 0 \end{bmatrix}, \]

kjer smo upoštevali simetričnost matrike \( p_{12} = p_{21} \). Nadalje

\[ PA^T = \begin{bmatrix} p_{12} & 0 \\ p_{22} & 0 \end{bmatrix}. \]

Poglejmo si še dinamičen šum, ki je definiran

\[ \Gamma Q \Gamma^T = \left\langle \Gamma \mathbf{w} \left( \Gamma \mathbf{w} \right)^T \right\rangle. \]

Vemo, da spremenljivka \( \dot{x} \) nima dinamičnega šuma in bo element \( q_{11} = 0 \). Poznamo pa dinamičen šum \( q_{22} = Q \), kjer je \( Q \) varianca dinamičnega šuma in je neznana. Hkrati pa so te matrike ponavadi diagonalne, torej

\[ \Gamma Q \Gamma^T = \begin{bmatrix} 0 & 0 \\ 0 & Q \end{bmatrix}. \]

Kovariančna matrika je torej

\[ \dot{P} = \begin{bmatrix} 2p_{12} & p_{22} \\ p_{22} & Q \end{bmatrix}. \]

Sedaj lahko iz matričnega sistema preberemo sistem parcialnih enačb

\[\begin{align*} \dot{p}_{11} &= 2p_{12} \\ \dot{p}_{12} &= p_{22} \\ \dot{p}_{22} &= Q. \end{align*} \]

Privzamemo \( Q \ne Q(t) \). Postopek rešévanja bo sledeč: rešili bomo spodnjo, potem srednjo in na koncu rezultat vstavili v prvo enačbo.

Integriranje spodnje enačbe nam poda rezultat

\[ p_{22} (t) = Q t + p_{22} (0). \]

Integriramo še srednjo enačbo ob upoštevanju rezultata spodnje enačbe

\[ p_{12} (t) = \frac{1}{2} Q ^2 + p_{22} (0)t + p_{12} (0). \]

Podobno, vstavimo enačbo v zadnjo enačbo in integriramo

\[ p_{11} (t) = \frac{1}{3} Q t ^3 + p_{22} (0) t ^2 + 2 p_{12 } (0) t + p_{11} (0). \]

Upoštevamo iz začetnih pogojev \( p_{22} = p_{12} = p_{11} = 0 \). Ostale nam enačba

\[ p_{11} (t) = \frac{1}{3} Q t ^3. \]

Varianca lege je torej

\[\begin{equation} \label{eq:1} \sigma_x ^2 (t_2) = p_{11} (t_2) = \frac{1}{3} Q t_2 ^3. \end{equation} \]

Še zmeraj moramo določiti varianco dinamičnega šuma, kar bomo naredili s pomočjo zadnjega neupoštevanega dejstva \( P \left( - R < x(t_1) < R \right) = \frac{2}{3} \).

Pomagali si bomo z dejstvom, da

\[ P \left( \sigma_1 (t_1) < x (t_1) < \sigma_x (t_1) \right) = 68 \%. \]

Ko primerjamo parametre, razberemo \( \sigma_x (t_1) = R \).

Pri prehodu skozi prva svetlobna vrata torej velja enačba

\[ \sigma_x ^2 (t_1) \cdot R ^2 = \frac{1}{3} Q ^3 t_1 ^3 \implies Q = \frac{3 R ^2}{t_1 ^3}. \]

Vstavimo vrednost \( Q \) v enačbo \(\ref{eq:1}\) in dobimo, da je varianca lege torej

\[ \sigma_x ^2 (t_2) = R ^2 \left( \frac{t_2}{t_1} \right) ^3 \]

Sedaj izračunajmo še časa \( t_1 \) in \( t_2 \)

\[\begin{align*} h &= \frac{1}{2} g t_1 ^2 2h &= \frac{1}{2} g t_2 ^2, \end{align*} \]

kar dobimo

\[ t_1 = \sqrt{\frac{2 h}{g}} \quad \text{ in } \quad t_2 = \sqrt{\frac{4h}{g}}. \]

Variance lege ob času \( t_2 \) je potem

\[ \sigma_x ^2 (t_2) = R ^2 \cdot *^{\frac{1}{2}} \implies \sigma_x (t_2) = R \cdot 8^{\frac{1}{4}}. \]

Vrnimo se nazaj k verjetnosti

\[ P \left( - R < x(t_2) < R \right) = 2 F \left( \frac{R}{\sigma_x (t_2)} \right) - 1 = 2 F \left( 8^{- \frac{1}{4}} \right) -1 = 2 F \left( 0.595 \right) = 45\%. \]

Vrednost \( F(0.595) \) smo prebrali iz tabele.

Note

Po prehodu skozi prva svetlobna vrata lahko opravimo meritev \( z_1 \) in opravimo izostreno oceno \( \hat{x} \). Na vratih torej uporabimo diskreten Kalmanov filter. Dobljene rezultate lahko potem propagiramo preko zveznega Kalmanovega filtra do drugih vrat.

Naloga: jašek

Imamo globok jašek napolnjen z vodo. Po sredini jaška spuščamo kroglice z maso \( m \).

Kroglico spustimo v vodo. Po dolgem času je varianca hitrosti enaka \( \frac{1}{4} \) začetne variance hitrosti.

a) Zapiši Kalmanov filter za vektor \( \mathbf{x} = [y, v]^T \) b) Poišči varianco hitrosti ob času \( \tilde{t} \), kjer je \( \beta \tilde{t} = 1 \). \( \beta \) je definiran s silo upora

\[ \mathbf{F}_u = - m \beta \mathbf{v}. \]

Upoštevaj tudi vpliv naključnih sil (aka dinamičen šum).

Uskladimo za začetek koordinatni sistem. Naj bo os \( y \) usmerjena navzdol v pozitivni smeri, vzporedno z jaškom. Kroglice spuščamo iz koordinate \( y = 0 \).

Dinamika - gibalna enačba

Podobno kot v prejšnji nalogi zapišemo 2. Newtonov zakon

\[ \sum\limits_i^{} F_i = m a_y = m \ddot{y}. \]

Na kroglico deluje sila teže v vzporedni smeri \( \mathbf{F} = m \mathbf{g} \), v negativni smeri pa še sila vzgona \( \mathbf{F}_{vzg} \) in sila upora \( \mathbf{F}_u \). V poljubni navpični smeri deluje pa še dinamičen šum \( F(t) \).

Drugi Newtonov zakon je za naš specifičen primer

\[ m \ddot{y} = mg - F_{vzg} - m \beta v + F(t). \]

Vsoto sile teže in sile vzgona bomo definirali kot efektivni težnosti pospešek

\[ F_{eff} = F_g - F_{vzg} = mg - F_{vzg} = m g_{eff}. \]

Preostane nam LDE 2. reda

\[ \ddot{y} = g_{eff} - \beta v + \frac{F(t)}{m}. \]

Ponovno želimo zapisati v vektorski obliki, da bomo lahko uporabili Kalmanov filter

\[ \dot{\mathbf{x}} = \begin{bmatrix} \dot{x} \\ \dot{v} \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 0 & - \beta \end{bmatrix} \begin{bmatrix} y \\ v \end{bmatrix} + \begin{bmatrix} 0 \\ g_{eff} \end{bmatrix} + \begin{bmatrix} 0 \\ w \end{bmatrix}. \].

Matrika \( AP \) je

\[ AP = \begin{bmatrix} 0 & 1 \\ 0 & -\beta \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{21} & p_{22} \end{bmatrix} = \begin{bmatrix} p_{12} & p_{22} \\ -\beta p_{12} & -\beta p_{22} \end{bmatrix}. \]

Transponiramo, da dobimo

\[ PA^T = \begin{bmatrix} p_{12} & - \beta p_{12} \\ p_{22} & - \beta p_{22} \end{bmatrix}. \]

Podobno kakor prej ima matrika dinamičnega šuma obliko

\[ \Gamma Q \Gamma^T = \begin{bmatrix} 0 & 0 \\ 0 & Q \end{bmatrix}. \]

Imamo torej matrično enačbo

\[ \mathbf{P} = \begin{bmatrix} 2 p_{12} & p_{22} - \beta p_{12} \\ p_{22} - \beta p_{12} & Q - 2 \beta p_{22} \end{bmatrix}. \]

Točka b

Iščemo \( \sigma_v ^2 (\tilde{t}) = p_{22} (\tilde{t}) \).

Za razliko od prejšnjega primera tokrat ne rabimo reševati celotnega sistema, ampak rešimo samo

\[ \dot{p}_{22} = Q - 2 \beta p_{22}. \]

V enačbi ne poznamo dinamičnega šuma \( Q \).

Zgornjo nehomogeno LDE lahko rešimo z uvedbo nove spremenljivke

\[ u = Q - 2 \beta p_{22}. \]

Novo spremenljivko odvajamo

\[ \dot{u} = - 2 \beta \dot{p}_{22}. \]

Namesto \( \dot{p}_{22} \) pa vstavimo \( u \) in rešujemo torej

\[ \frac{\dot{u}}{2 \beta} = - u \implies u(t) = C e^{- 2 \beta t} \]

z začetnim pogojem \( u(0) = C \).

Zapišemo

\[ u(t) = u(0) e^{- 2 \beta t}, \]

in upoštevamo definicijo

\[\begin{equation} \label{eq:2} Q - 2 \beta p_{22} (t) = \left( Q - 2 B p_{22} (0) \right) e^{- \beta t}. \end{equation} \]

Po dolgem času \( t \to \infty \), dobimo vrednost

\[ Q - 2 \beta p_{22} (\infty) = 0 \implies Q = 2 \beta p_{22} (\infty). \]

Dobljeno vrednosti vstavimo v enačbo \(\ref{eq:2}\)

\[ 2 \beta p_{22} (\infty) - 2 \beta p_{22} (t) = \left( 2 \beta p_{22}(\infty) - 2 \beta p_{22} (0) \right) e^{- 2 \beta t}. \]

Pokrajšamo člene in nam preostane

\[ p_{22} (t) = p_{22}(\infty) - \left( p_{22} (\infty) - p_{22} (0) \right) e^{- 2 \beta t}. \]

Poznamo podatek \( p_{22} (\infty) = \frac{1}{4} p_{22} (0) \) in ga upoštevamo

\[ p_{22} (t) = \frac{1}{4} p_{22} (0) t \frac{3}{4} p_{22} (0) e^{- 2 \beta t} = \frac{p_{22} (0)}{4} \left( 1 + 3 e^{- 2 \beta t} \right). \]

Ob času \( \tilde{t} \) je potem

\[\begin{align*} p_{22} (\tilde{t}) &= \frac{p_{22} (0)}{4} \left( 1 + 3 e^{- 2 \beta \tilde{t}} \right) \\ &= \frac{p_{22} (0)}{4} \left( 1 + 3 e^{-2} \right) = \frac{p_{22} (0)}{3} = 0.35 \cdot p_{22} (0). \end{align*} \]