Całkowanie numeryczne

Całkowanie numeryczne jest zadaniem dość złożonym. Najbardziej elementarne techniki z tego zakresu to reguły o małej złożoności i ustalonym rozmiarze kroku. W metodzie tej daną transmitancję układu ciągłego ${\displaystyle G(s)}$ zastępuje się przez równanie różniczkowe, a następnie wyprowadza się równania różnicowe będące aproksymacją równań różniczkowych.

Niech dana będzie transmitancja integratora analogowego:

${\displaystyle G(s)={\frac {U(s)}{E(s)}}={\frac {1}{s}},}$

gdzie ${\displaystyle E(s)}$ oraz ${\displaystyle U(s)}$ są odpowiednio transformatami wejścia i wyjścia integratora. Dla integratora tego można określić równoważne równanie różniczkowe

${\displaystyle {\frac {du(t)}{dt}}=e(t),}$

które można zapisać w postaci całkowej:

${\displaystyle u(t)=\int \limits _{0}^{t}e(\tau )d\tau .}$

Wiele reguł opiera się na właściwej sobie metodzie aproksymacji składnika powiększania pola (pod krzywą funkcji, która w powyższym wzorze podlega całkowaniu). Należą do nich:

• reguła prostokąta wprzód,
• reguła prostokąta wstecz,
• reguła trapezu.

W regule prostokątnej wprzód obszar aproksymuje się przez prostokąt wyznaczany wprzód od chwili kT do chwili kT+T i bierze jako amplitudę prostokąta wartość napotkaną w kT. Szerokość takiego prostokąta wynosi T. Można więc zapisać równanie w pierwszej aproksymacji:

${\displaystyle u_{1}(kT+T)=u_{1}(kT)+Te(kT),}$

gdzie wyrażenie ${\displaystyle u_{1}(kT)}$ reprezentuje obszar pod całkowaną krzywą e(t) w przedziale od t = 0 do t = kT. Po zastosowaniu transformaty Z do powyższej zależności otrzymuje się:

${\displaystyle G_{F}(z)={\frac {U_{1}(z)}{E(z)}}={\frac {T}{z-1}}={\frac {1}{\frac {z-1}{T}}}.}$

W regule prostokątnej wstecz obszar aproksymuje się przez prostokąt wyznaczany wstecz od chwili kT do kT-T i bierze jako amplitudę prostokąta wartość napotkaną w kT. Szerokość takiego prostokąta wynosi T. Można więc zapisać równanie w pierwszej aproksymacji:

${\displaystyle u_{2}(kT)=u_{2}(kT-T)+Te(kT).}$

Po zastosowaniu transformaty Z do powyższej zależności otrzymuje się:

${\displaystyle H_{B}(z)={\frac {U_{2}(z)}{E(z)}}={\frac {zT}{z-1}}={\frac {1}{{\frac {1}{T}}{\frac {z-1}{z}}}}.}$

W regule trapezu obszar aproksymuje się przez pole trapezu umieszczonego pod całkowaną krzywą. Równanie aproksymacji ma wówczas postać:

${\displaystyle u_{3}(kT+T)=u_{3}(kT)+{\frac {T}{2}}[e(kT)+e(kT+1)].}$

Po zastosowaniu transformaty Z do powyższej zależności otrzymuje się:

${\displaystyle H_{T}(z)={\frac {U_{3}(z)}{E(z)}}={\frac {T}{2}}{\frac {z+1}{z-1}}={\frac {1}{{\frac {2}{T}}{\frac {z-1}{z+1}}}}.}$

Metoda reguły trapezu jest również znana jako metoda Tustina lub pod nazwą transformacji biliniowej (zob. też płaszczyzna w). Metoda projektowania wykorzystująca tę regułę polega na tym, że daną transmitancję ciągłą, ${\displaystyle G(s),}$ równoważną transmitancja dyskretnej wyznacza się przez podstawienie:

${\displaystyle G_{T}(z)=G(s)|_{s={\frac {2}{T}}{\frac {z-1}{z+1}}}.}$

Każda z powyższych aproksymacji może być potraktowana jako przekształcenie płaszczyzny s na płaszczyznę z.

Porównując transmitancje operatorowe z trzema aproksymacjami dyskretnymi można zauważyć, że transmitancję dyskretną można uzyskać bezpośrednio z transformaty operatorowej podstawiając za zmienną zespoloną „s” jej aproksymatę.

W przypadku reguły prostokąta wprzód jest to podstawienie ${\displaystyle s\leftarrow {\frac {z-1}{T}}.}$

W przypadku reguły prostokąta wstecz jest to podstawienie ${\displaystyle s\leftarrow {\frac {z-1}{Tz}}.}$

W przypadku reguły trapezu jest to podstawienie ${\displaystyle s\leftarrow {\frac {2}{T}}{\frac {z-1}{z+1}}.}$

Szczególnie interesujące jest to, że reguła bilinearna odwzorowuje stabilną półpłaszczyznę ${\displaystyle s}$ dokładnie na stabilny obszar płaszczyzny z, przy tym cała oś ${\displaystyle j\omega }$ płaszczyzny s jest skompresowana na długości obwodu okręgu jednostkowego.

Dyskretyzacja szumu procesu

Numeryczna ewaluacja ${\displaystyle \mathbf {Q} _{d}}$ jest nieco bardziej złożona z uwagi na całkę eksponenty macierzy. Można ją, jednakże, wyliczyć poprzez skonstruowanie najpierw macierzy, a następnie wyliczenie na komputerze jej eksponenty:

${\displaystyle \mathbf {F} ={\begin{bmatrix}-\mathbf {A} &\mathbf {Q} \\\mathbf {0} &\mathbf {A} ^{T}\end{bmatrix}}T}$

${\displaystyle \mathbf {G} =e^{\mathbf {F} }={\begin{bmatrix}\dots &\mathbf {A} _{d}^{-1}\mathbf {Q} _{d}\\\mathbf {0} &\mathbf {A} _{d}^{T}\end{bmatrix}}.}$

Zdyskretyzowany szum procesu jest wówczas wyliczany poprzez przemnożenie transponowanej dolnej, prawej partycji macierzy G z górną, prawą partycją macierzy G:

${\displaystyle \mathbf {Q} _{d}=(\mathbf {A} _{d}^{T})^{T}(\mathbf {A} _{d}^{-1}\mathbf {Q} _{d}).}$

Wyprowadzenie

Rozpoczynając z modelem ciągłym

${\displaystyle \mathbf {\dot {x}} (t)=\mathbf {A} \mathbf {x} (t)+\mathbf {B} \mathbf {u} (t)}$

wiadomo, że eksponenta macierzy jest następująca:

${\displaystyle {\frac {d}{dt}}e^{\mathbf {A} t}=\mathbf {A} e^{\mathbf {A} t}=e^{\mathbf {A} t}\mathbf {A} }$

i przez wcześniejsze przemnożenie modelu uzyskuje się:

${\displaystyle e^{-\mathbf {A} t}\mathbf {\dot {x}} (t)=e^{-\mathbf {A} t}\mathbf {A} \mathbf {x} (t)+e^{-\mathbf {A} t}\mathbf {B} \mathbf {u} (t),}$

co zapisać można jako

${\displaystyle {\frac {d}{dt}}(e^{-\mathbf {A} t}\mathbf {x} (t))=e^{-\mathbf {A} t}\mathbf {B} \mathbf {u} (t),}$

a następnie całkując:

${\displaystyle e^{-\mathbf {A} t}\mathbf {x} (t)-e^{0}\mathbf {x} (0)=\int _{0}^{t}e^{-\mathbf {A} \tau }\mathbf {B} \mathbf {u} (\tau )d\tau ,}$

${\displaystyle \mathbf {x} (t)=e^{\mathbf {A} t}\mathbf {x} (0)+\int _{0}^{t}e^{\mathbf {A} (t-\tau )}\mathbf {B} \mathbf {u} (\tau )d\tau ,}$

co jest rozwiązaniem analitycznym dla modelu ciągłego.

Teraz należy zdyskretyzować powyższe wyrażenie. Można przyjąć, że ${\displaystyle u}$ jest stała podczas każdego kroku czasowego

${\displaystyle \mathbf {x} [k]\ {\stackrel {\mathrm {def} }{=}}\ \mathbf {x} (kT),}$

${\displaystyle \mathbf {x} [k]=e^{\mathbf {A} kT}\mathbf {x} (0)+\int _{0}^{kT}e^{\mathbf {A} (kT-\tau )}\mathbf {B} \mathbf {u} (\tau )d\tau ,}$

${\displaystyle \mathbf {x} [k+1]=e^{\mathbf {A} (k+1)T}\mathbf {x} (0)+\int _{0}^{(k+1)T}e^{\mathbf {A} ((k+1)T-\tau )}\mathbf {B} \mathbf {u} (\tau )d\tau ,}$

${\displaystyle \mathbf {x} [k+1]=e^{\mathbf {A} T}\left[e^{\mathbf {A} kT}\mathbf {x} (0)+\int _{0}^{kT}e^{\mathbf {A} (kT-\tau )}\mathbf {B} \mathbf {u} (\tau )d\tau \right]+\int _{kT}^{(k+1)T}e^{\mathbf {A} (kT+T-\tau )}\mathbf {B} \mathbf {u} (\tau )d\tau .}$

Wyrażenie w nawiasie można zapisać jako ${\displaystyle \mathbf {x} [k]}$ a drugie wyrażenie można uprościć przez podstawienie ${\displaystyle v=kT+T-\tau .}$ Ponadto można przyjąć, że ${\displaystyle \mathbf {u} }$ jest stałe podczas całkowania, co z kolei daje:

${\displaystyle \mathbf {x} [k+1]=e^{\mathbf {A} T}\mathbf {x} [k]+\left(\int _{0}^{T}e^{\mathbf {A} v}dv\right)\mathbf {B} \mathbf {u} [k]=e^{\mathbf {A} T}\mathbf {x} [k]+A^{-1}\left(e^{\mathbf {A} T}-I\right)\mathbf {B} \mathbf {u} [k],}$

co stanowi dokładne rozwiązanie dyskretyzowanego problemu.

Aproksymacje

Dokładna dyskretyzacja czasami może być trudna z uwagi na dużą eksponentę macierzy i związane z tym operacje całkowania. Znacznie łatwiej wyliczyć, w oparciu o nią, przybliżony model dyskretny dla małych kroków czasowych ${\displaystyle e^{\mathbf {A} T}\approx \mathbf {I} +\mathbf {A} T.}$ Przybliżone rozwiązanie przyjmuje wówczas postać:

${\displaystyle \mathbf {x} [k+1]\approx (\mathbf {I} +\mathbf {A} T)\mathbf {x} [k]+(\mathbf {I} T+{\frac {1}{2}}\mathbf {A} T^{2})\mathbf {B} \mathbf {u} [k],}$

co można dalej aproksymować jeśli ${\displaystyle {\frac {1}{2}}\mathbf {A} T^{2}}$ jest małe; co daje:

${\displaystyle \mathbf {x} [k+1]\approx (\mathbf {I} +\mathbf {A} T)\mathbf {x} [k]+T\mathbf {B} \mathbf {u} [k].}$

Inne możliwe aproksymacje to: ${\displaystyle e^{\mathbf {A} T}\approx \left(\mathbf {I} -\mathbf {A} T\right)^{-1}}$ i ${\displaystyle e^{\mathbf {A} T}\approx \left(\mathbf {I} +{\frac {1}{2}}\mathbf {A} T\right)\left(\mathbf {I} -{\frac {1}{2}}\mathbf {A} T\right)^{-1}.}$ Każda z nich ma inne własności związane ze stabilnością. Ostatnia znana jest jako transformacja Tustina (transformacja bilinearna) i zachowuje stabilność lub odpowiednio niestabilność układu czasu ciągłego.