2026/03/10
Dinamični zakon spremenljivke nam pove, kako se količina spreminja in ga podamo v obliki diferencialne enačbe. Za naš sistem \( S \) bo diferenčna enačba, ki služi kot približek diferencialne enačbe
\[ \dot{x} = A(t) x + C(t). \]
V diskretni sliki lahko majhne \( \dot{x} \) zapišemo kot
\[ \frac{x((n + 1)T) - x(nT)}{T} = A(nT) x (nT) - C (nT), \]
kjer je \( T \) čas vzorčenja. Slednje lahko dodatno preobrazimo v
\[ x_{n + 1} = \left( 1 + A(n T) T \right) x_n + C(n T) T. \]
Dobili smo zvezi med koeficienti diferenčne in diferencialne enačbe
\[ \Phi_n = 1 + A(nT) T \quad \text{ in } \quad c_n = C(nT). \]
Dinamični zakon v diferenčni obliki je potemtakem
\[\begin{equation} \label{eq:1} x_{n + 1} = \Phi_n x_n + c_n \left( + \Gamma_n w_n \right), \end{equation} \]kjer zadnji člen predstavlja dinamični šum. Vzrok za ta šum je, če ne poznamo dobro dinamike sistema. Členi \( w_n \) so normalno porazdeljeni s povprečjem \( 0 \) in varianco dinamičnega šuma \( Q_n \)
\[ \left\langle w_n w_{n'} \right\rangle = Q_n \delta_{n n'}. \]
Ob času \( t = n T \) poznamo optimalno oceno \( \hat{x}_n \) in njeno disperzijo \( \hat{\sigma}_n ^2 \). Ob času \( (n + 1)T \) izmerimo \( z_{n + 1} \). Ocena ob času \( t = (n + 1)T \) se bo premaknila glede na \(\ref{eq:1}\) in bo enaka
\[ \bar{x}_{n + 1} = \Phi_n \hat{x}_n + c_n. \]
Disperzija ob tem času pa bo
\[ \bar{\sigma}_{n + 1} ^2 = \Phi_n ^2 \hat{\sigma}_n ^2 + \Gamma_n ^2 Q_n. \]
Z meritvijo \( z_{n + 1} \) lahko sedaj izostrimo oceno, kakor smo to naredili v poglavju optimalnega združevanja meritev
\[ \hat{x}_{n + 1} = \bar{x}_{n + 1} + \frac{\bar{\sigma}_{n + 1} ^2}{\bar{\sigma}_{n + 1} ^2 + \sigma_{n + 1}} \left( z_{n + 1} - \bar{x}_{n + 1} \right). \]
Njena disperzija je
\[ \bar{\sigma}_{n + 1} ^2 = \frac{\bar{\sigma}_{n + 1} ^2 \sigma_{n + 1} ^2}{\bar{\sigma}_{n + 1} ^2 + \sigma_{n + 1} ^2} = \bar{\sigma}_{n + 1} ^2 - \frac{\bar{\sigma}_{n + 1} ^4}{\bar{\sigma}_{n+1}^2 + \sigma_{n + 1} ^2}. \]
Uvedemo nove spremenljivke in sicer
\[\begin{array}{|c|c|} \hline \sigma_{n + 1} ^2 & R_{n + 1} \\ \bar{\sigma}_{n + 1} ^2 & M_{n + 1} \\ \hat{\sigma}_{n + 1} ^2 & P_{n + 1} \\ \hline \end{array} \]V primeru, da ne merimo, potem ni drugega koraka - ostrenja ocene. Napoved je najboljša ocena
\[ \hat{x}_{n + 1} = \bar{x}_{n + 1} \quad \text{ in } \quad P_{n + 1} = M_{n + 1}. \]
Kroglica se zakotali od vrha petih stopnic do dna. Vse stopnice so natančno enako visoke \( \bar{H} = 30 \mathrm{cm} \) in imajo enako ničelno negotovost višine \( \sigma_H ^2 = 0 \). Kroglica se na vsaki stopnici odbije natančno enkrat in ob vsakem odboju izgubi polovico svoje kinetične energije v vertikalni smeri.
Izračunaj največjo višino zadnjega, 5., odboja in določi njegovo disperzijo. Med odboji po stopnicah ne merimo višine. Začetna višina kroglice je \( 0 \) in je znana samo do variance \( \sigma_0 \). Zanemari dinamični šum.
Za začetek reševanja te naloge moramo najprej poiskati rekurzivno zvezo, ki nam bo določila dinamiko. Želimo imeti rekurzivno zvezo oblike Kalmanovega filtra
\[ h_{n + 1} = \Phi_n h_n + c_n + \Gamma_n w_n, \]
kjer je \( h \) spremenljivka, ki nas zanima - višina. Označimo izgubo energije po vsakem odboju z \( \delta = \frac{1}{2} \). \( h_0 \) bo začetna višina žoge, ki je v našem primeru \( h_0 = 0 \).
Potencialna energije kroglice po prvem odboju mora biti enaka polovici potencialne energije pred odbojem, torej \( U_1 = \delta U_0 \). Ob upoštevanju definicije potencialne energije dobimo zvezo
\[ mg h_1 = \delta mg \left( H + h_0 \right) \implies h_1 = \delta h_0 + \delta H. \]
Zaradi ničelne negotovosti velja zapis \( \bar{H} = H \). Po drugem odboju \( U_2 = \delta U_1 \) pridemo do podobne zveze
\[ h_2 = \delta h_1 + \delta H. \]
Po \( n \)-tem odboju bo torej višina
\[ h_{n + 1} = \delta h_n + \delta H. \]
Sedaj lahko odčitamo vrednosti koeficientov iz diferenčne enačbe. To so \( \Phi_n = \delta \) in \( c_n = \delta H \). Hkrati je \( \Gamma_n = 0 \), saj smo zanemarili dinamični šum.
Sedaj želimo izraziti višin \( n + 1 \) odboja s \( h_0 \). Za \( n = 2 \) to pomeni
\[ \hat{h}_2 = \bar{h}_2 = \delta \hat{h}_1 + \delta H = \delta \left( \delta \hat{h}_0 + \delta H \right) + \delta H = \delta ^2 \hat{h}_0 + \left( \delta ^2 + \delta \right) H. \]
Tukaj smo upoštevali, da ker ni vmesnih meritev, je napoved \( \bar{h} \) najboljša ocena \( \hat{h} \). Nadalje
\[ \hat{h}_3 = \delta \hat{h}_2 + \delta H = \delta ^3 \hat{h}_0 + H \left( \delta + \delta ^2 + \delta ^3 \right). \]
Torej bo splošna zveza
\[ h_n = \delta^n h_0 + H \sum\limits_{i = 1}^n \delta^i. \]
Vsoto lahko tudi poenostavimo
\[ S_n = \sum\limits_{i = 1}^n \delta^i = \delta \left( 1 + S_n - \delta^n \right) \implies S_n = \delta \frac{1 - \delta^n}{1 - \delta}. \]
Izostrena ocena po \( n \)-odbojih bo torej
\[ \hat{h}_n = \delta^n \hat{h}_0 + H \delta \frac{1 - \delta^n}{1 - \delta}. \]
V našem primeru za \( n = 5 \) in \( \delta = \frac{1}{2} \) je to enako
\[ \hat{h}_5 = H \left( 1 - \frac{1}{2^5} \right) = \frac{31}{32} H. \]
Sedaj pa lahko preko Kalmanovega filtra določimo še varianco petega odboja. Ker nimamo vmesnih meritev, pomeni da ne moremo izostriti vmesnih ocen. Posledično bomo pri dinamiki
\[ M_{n + 1} = \Phi_n ^2 P_n + \underbrace{\Gamma_n ^2 Q_n}_{= 0} \]
označili \( P_{n + 1} = M_{n + 1} \). Torej bo veljalo
\[ P_{n + 1} = \Phi_n ^2 P_n = \delta ^2 P_n. \]
Izraženo s \( P_0 \) potem velja
\[ P_n = \left( \delta ^2 \right)^n P_0. \]
Začetna varianca je \( P_0 = \sigma_0 ^2 \), torej bo \( n \)-ti odboj
\[ P_n = \delta^{2n} \hat{\sigma}_0 ^2. \]
Zanima nas disperzija, ki je
\[ \hat{\sigma}_n = \sqrt{P_{n}} = \delta^n \sigma^2_0. \]
Disperzija 5. odboja bo tako
\[ \hat{\sigma}_5 = \delta ^5 \hat{\sigma}_0 = \frac{1}{32} \hat{\sigma}_0. \]
Kapico smo dali disperziji, saj je to naša optimalna ocena. Disperzijo pa bi lahko še znižali, če bi dodatno poznali meritve med posameznimi odboji, vendar se disperzija manjša od začetne vrednost že zaradi dinamike sistema
\[ \lim_{n \to \infty} \hat{\sigma}_n =0. \]
Prejšnji nalogi dodamo dinamični šum \( \delta_H \ne 0 \). To pomeni, da smo višini stopnic dodali disperzijo.
V tem primeru zapišemo višino stopnice kot
\[ H = \bar{H} + \eta, \]
kjer je \( \left\langle \eta \right\rangle = 0 \) in \( \left\langle \eta ^2 \right\rangle = \delta_H ^2 \).
V povprečju so stopnice enako visoke \( \bar{H} \).
Nova dinamika, Kalmanov filter pa je sedaj
\[ h_{n + 1} = \delta h_n + \delta \left( \bar{H} + \eta \right) = \delta h_n + \delta \bar{H} + \delta \eta. \]
Posamezni koeficienti enačbe \(\ref{eq:1}\) so sedaj
\[ \Phi_n = \delta,\qquad c_n = \delta \bar{H} \quad \text{ in } \quad \Gamma_n w_n = \delta \eta \]
Slednji člen predstavlja dinamični šum.
Ocena višine odboja še zmeraj ostane ista zaradi \( \left\langle \eta \right\rangle = 0 \)
\[ \bar{h}_{n + 1} = \hat{h}_{n + 1} = \delta \hat{h}_n + \delta \bar{H}. \]
Do spremembe zaradi dinamičnega šuma pride pri varianci. Po izpeljavi v teoretičnem delu je
\[ P_{n + 1} = M_{n + 1} = \Phi_{n + 1} ^2 P_n + \Gamma_n ^2 Q_n. \]
Spomnimo, da je varianca dinamičnega šuma
\[ \left\langle \left( \Gamma_n w_n \right) ^2 \right\rangle = \Gamma_n ^2 \left\langle w_n ^2 \right\rangle = \Gamma_n ^2 Q_n. \]
Varianco dinamičnega šuma z našo disperzijo zapišemo
\[ \Gamma_n ^2 Q_n = \left\langle \left( \delta \eta \right) ^2 \right\rangle = \delta ^2 \left\langle \eta ^2 \right\rangle = \delta ^2 \sigma_H ^2, \]
kjer smo upoštevali navodila naloge.
Disperzija izostrene ocene po prvem odboju je torej
\[ P_1 = \delta ^2 \hat{\sigma}_0 ^2. \]
Po drugem odboju je disperzija
\[ P_2 = \delta ^2 \left( \delta ^2 P_0 + \delta ^2 \sigma_H ^2 \right) + \sigma_H ^2 = \delta ^4 P_0 + \delta_H ^2 \left( \delta ^2 + \delta ^4 \right). \]
Po \( n \)-tem odboju velja zveza
\[ P_n = \delta^{2n} \hat{\sigma}_0 ^2 + \sigma_H ^2 \sum\limits_{i = 1}^n \delta^{2i}. \]
Podobno kakor prej lahko zgornjo vsoto pretvorimo na zvezo
\[ S_n = \sum\limits_{i = 1}^n \delta^{2i} = \delta ^2 \left( 1 + S_n - \delta^{2n} \right) \implies S_n = \delta ^2 \frac{1 - \delta^{2n}}{1 - \delta ^2}. \]
Disperzija izostrene ocene glede na Kalmanov filter po dveh odbojih je torej
\[ P_2 = \delta ^4 \hat{\sigma}_0 ^2 + \delta_H ^2 \left( \delta ^2 + \delta ^4\right) \]
Napoved višine je po drugem odboju
\[ h_2 = \delta ^2 h_0 + H \left( \delta + \delta ^2 \right) = f(h_0, H), \]
kar je funkcija dveh spremenljivk. Vsaka spremenljivka \( h_0 \) in \( H \) ima svojo napako in negotovost lahko izračunamo s širjenjem napake po Gaussovi porazdelitvi.
\[ \hat{\sigma}_2 ^2 = \delta ^4 + \left( \delta + \delta ^2 \right) ^2 \sigma_H ^2 = \delta ^4 \hat{\sigma}_0 ^2 + \left( \delta ^2 + \delta ^4 + 2 \delta ^3 \right) \sigma_H ^2 \]
Po Kalmanovem filtru dodatnega člena \( 2 \delta ^3 \) ni! Pri širjenju napake po Gaussu smo se pretvarjali, da je vsaka stopnica meritev iste naključne spremenljivke. \( 2 \delta ^3 \) je korelacija med meritvami \( H \).
Širjenje napake po Gaussu bi bilo pravilno, če bi bila vsaka meritev posamezne stopnice svoje spremenljivke \( H_1, H_2, \ldots \).
Višina po \( n \)-tem odboju je potem funkcija \( n + 1 \) spremenljivk
\[ h_n = f \left( h_0, H_1, H_2, \ldots, H_n \right). \]
V tem primeru dobimo pravilno rešitev, saj so meritve med seboj nekorelirane.
Imamo \( k \) spremenljivk \( x_1, x_2, \ldots, x_k \), ki jih zapišemo v vektorski obliki kot
\[ \mathbf{x} = \begin{bmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{k} \end{bmatrix}. \]
Ocena in izostrena ocena pa bosta v vektorski obliki
\[ \bar{\mathbf{x}} = \begin{bmatrix} \bar{x}_{1} \\ \bar{x}_2 \\ \vdots \\ \bar{x}_k \end{bmatrix} \quad \text{ in } \quad \hat{\mathbf{x}} = \begin{bmatrix} \hat{x}_{1} \\ \hat{x}_{2} \\ \vdots \\ \hat{x}_k \end{bmatrix}. \]
Za vektorske spremenljivke moramo tako kot skalarnih zapisati dinamiko, ki bo ponovno oblike
\[ \mathbf{x}_{n + 1} = \Phi_n \mathbf{x}_n + \mathbf{c}_n + \Gamma_n \mathbf{w}_n. \]
Razlika je sedaj v tem, da je \( \Phi_n \) sedaj matrika.
Za vzorčni čas \( T \) in čas \( t = n T \) želimo narediti napoved za čas \( t = (n + 1)T \). Ob času \( nT \) imamo izostreno oceno \( \hat{\mathbf{x}} \) ter pripadajočo kovariančno matriko izostrene ocene \( P_n \). Napoved ocene za novi čas je zaradi dinamike
\[ \bar{\mathbf{x}}_{n + 1} = \Phi_n \hat{\mathbf{x}}_n + \mathbf{c}_n \]
s pripadajočo kovariančno matriko
\[ M_{n + 1} = \left\langle \left( \bar{\mathbf{x}}_{n + 1} - \mathbf{x}_{n + 1} \right) \left( \bar{\mathbf{x}}_{n + 1} - \mathbf{x}_{n + 1} \right)^T \right\rangle. \]
Razliko napovedi ocene in spremenljvke, ki sledi diferenčni enačbi je
\[ \bar{\mathbf{x}}_{n + 1} - \mathbf{x}_{n + 1} = \Phi_n \hat{\mathbf{x}}_n + c_n - \left( \Phi_n \mathbf{x}_n + \mathbf{c}_n + \Gamma_n \mathbf{w}_n \right) = \Phi_n (\hat{\mathbf{x}}_n - \mathbf{x}_n) - \Gamma_n \mathbf{w}_n. \]
Kovariančna matrika npovedi ocene je zato
\[\begin{align*} M_{n + 1} &= \left\langle \left( \Phi_n \left( \hat{\mathbf{x}}_n - \mathbf{x}_n \right) - \Gamma_n \mathbf{w}_n \right) \left( \Phi_n \left( \hat{\mathbf{x}}_n - \mathbf{x}_n \right) - \Gamma_n \mathbf{w}_n\right)^T \right\rangle \\ &= \Phi_n P_n \Phi_n^T + \Gamma_n Q_n \Gamma_n ^T, \end{align*} \]kjer je
\[ P_n = \left\langle \left( \hat{\mathbf{x}}_n - \mathbf{x} \right)\left( \hat{\mathbf{x}}_n - \mathbf{x} \right)^T \right\rangle \quad \text{ in } \quad \left\langle \mathbf{w}_n \mathbf{w}_n^T \right\rangle = Q_n. \]
Matriki \( Q_n \) pravimo kovariančna matrika dinamičnih šumov.
Brez dodatnih meritev velja \( P_{n + 1} = M_{n + 1} \) in \( \hat{\mathbf{x}}_{n + 1} = \bar{\mathbf{x}}_{n + 1} \).
Dve kladi z masama \( m_1 \) in \( m_2 \) s pripadajočima hitrostima \( v_1 \) ter \( v_2 \) prožno trčita. Hitrosti imata vsaka svojo nedoločenost \( \sigma_1 \) in \( \sigma_2 \). Izračunaj kovariančno matriko po trku za vektor \( \mathbf{v} = [v_1 ', v_2 ']^T \), kjer sta to hitrosti teles po trku. Predpostavi, da sta začetni hitrosti nekorelirani in telesni dinamiki popolnoma poznani.
Pred trkom imamo vektor
\[ \mathbf{v} = \begin{bmatrix} v_{1} \\ v_{2} \end{bmatrix} \]
s pripadajočo kovariančno matriko
\[ P = \begin{bmatrix} \sigma_{1} ^2 & 0 \\ 0 & \sigma_2 ^2 \end{bmatrix}. \]
Išcemo dinamiko oblike
\[ \mathbf{v}' = \Phi \mathbf{v} + \mathbf{c} + \Gamma \mathbf{w}. \]
V navodilih naloge piše, da imamo popolnoma prožni trk, zato bomo nalogo reševali v težiščnem sistemu (CMS). Hitrost posamezne klade označimo z
\[ u_1 = v_1 - v_T \quad \text{ in } \quad u_2 = v_2 - v_T, \]
kjer je \( v_T \) hitrost težišča sistema definirana kot
\[ v_T = \frac{m_1 v_1 + m_2 v_2}{m_1 + m_2}. \]
Iz zakona o ohranitvi gibalne količine \( \sum\limits_{}^{} p_i = 0 \) dobimo dve enačbi za pred in po trku
\[\begin{align*} 0 &= m_1 u_1 + m_2 u_2 \\ 0 &= m_1 u_1 ' + m_2 u_2 '. \end{align*} \]Iz danih enačb lahko izrazimo hitrost druge klade glede na hitrost prve klade
\[\begin{align*} u_2 &= - \frac{m_1}{m_2} u_1 \\ u_2 ' &= - \frac{m_1}{m_2} u_1 '. \end{align*} \]Upoštevati moramo tudi zakon o ohranitvi kinetične energije \( \sum\limits_i^{} T_i = \sum\limits_i^{} T_i ' \). Dobimo enačbo
\[ m_1 u_1 ^2 + m_2 u_2 ^2 = m_1 u_1 ' ^2 + m_2 u_2 ' ^2. \]
Upoštevamo identitete dobljene iz zakona o gibalni količini in zapišemo
\[ u_1 ^2 \left( m_1 + \frac{m_1 ^2}{m_2} \right) = u_1' ^2 \left( m_1 + \frac{m_1 ^2}{m_2} \right), \]
iz česar dobimo
\[ u_1 ' = \pm u_1. \]
Za hitrost po trku si bomo izbrali negativen predznak hitrost, saj drugače do trka ni prišlo! Analogno pridemo tudi do \( u_2 ' = - u_2 \).
Po dobljenih identitetah lahko ponovno preidemo iz težiščnega sistema nazaj v normalen sistem. Za \( i = 1, 2 \) dobimo izraza
\[\begin{align*} u_i ' &= - u_i \\ \left( v_i ' - v_T \right) &= - \left( v_i - v_T \right) \\ v_i ' &= - v_i + 2 v_T. \end{align*} \]