A ballisztikus pálya állapottér (state-space) modellje


Néhány napja a Tiborublog és Lemilblog szerzőinek és olvasóinak Facebook csoportjában kisebb kioktatást kaptam a ballisztikus pálya számítási módjából. Nem probléma, mivel nem értek hozzá. Az eset azonban felpiszkált egy kicsit és utána olvastam. Feltűnt, hogy több tudományos cikk is state-space modellt1 használja a pályaadatok számítására. Ezzel szemben egyszerűbb esetekben – mint például ferde hajítás – nem ez a gyakorlat. Sőt ezekre egyáltalán nem találtam state-space modelleket. Ez azért érdekes, mert oktatási célból létezniük kellene. Ezért úgy gondoltam, ha már így alakult megírom.


Frissítés:

2021-11-23 — Akos Gabrial kommentjei alapján magyarítása a “state-space”-nek, és kód javítás

A állapottér modell

Ez a modell típus vektorként ábrázolja egy zárt rendszer aktuális állapotát az euklideszi térben. A jó dolog benne, hogy egyrészt nagyon tömör formában megjeleníti az egész rendszer adott időpillanatban felvett állapotát. Másrészt lehetővé teszi a jövőben felvett állapot gyors és egyszerű számítását. Ez így talán homályos, de a következő egyszerű példák remélhetőleg megvilágítják az ilyen típusú modellezés előnyét. Annyit még érdekességként, hogy Kálmán Rudolf Emil magyar származású mérnök-matematikus úttörő munkát végzet a kidolgozásában.

A ferde hajítás 

Kezdjük a legegyszerűbb esettel, a ferde hajítással. A modellt eredetileg Galileo dolgozta ki 1638-ban. Azzal a lényeges egyszerűsítéssel élt, hogy az eldobott testre ható erők közül csak a gravitációs erőt vette számításba. A linkelt Wikipédia lap részletesen foglalkozik a számítással, így én csak arra a kérdésre koncentrálok, hogy miként számítjuk ki az elhajított test aktuális pozícióját. Az egyszerűség kedvéért csak két dimenzióban fogom bemutatni a rendszert.

Van egy ferdén elhajított testünk, aminek a dobás helyétől számított vízszintes távolságát jelöljük x-el, a függőleges magasságát pedig z-vel.2 Ebben az esetben a t időpillanatban felvett x értéket a következő módon számíthatjuk:

(1)  x_t = v_0 \cdot t \cdot cos(\theta)

Az ugyanekkor a felvett magasság pedig:

(2)  y_t = v_0 \cdot t \cdot sin(\theta)

Ahol: 

  • v_0  —  a ferde hajítás kezdősebbegése
  • θ  —  a hajítás szöge
  • t  —  eltelt idő a hajítás kezdete óta
  • x_t — az aktuális vízszintes távolság a kiindulási ponttól
  • y_t— az aktuális magasság

Ez nagyon egyszerű. Például így valósítható meg Pythonban:

import numpy as np
from math import radians, cos, sin
from scipy.constants import g

# kezdő sebesség
v = 300
# theta (kilövési szög)
theta = radians(50)

# pozíciók az adott időpillanádban
x_t = []
z_t = []
# mérési időpillanatok
t = np.linspace(0, 70, num=10000).tolist()

# az adott időpillanathoz tartozó pozíció számítása
for k in t:
    x = ((v*k)*np.cos(theta)) 
    z = ((v*k)*np.sin(theta))-((0.5*g)*(k**2))
    
    # mentjük a számított pozíciót
    x_t.append(x)
    z_t.append(z)

# szűrjük azokat az időpillanatokat amelyek 0 magasság alatt vannak 
p = [i for i, j in enumerate(z_t) if j < 0] 
for i in sorted(p, reverse = True):
    del x_t[i]
    del z_t[i]
    del t[i]

Az eredmény pedig:

A állapottér modell kulcsa

Most pedig nézzük meg ugyanezt másképp. Képzeljük el, hogy egy fényképet készítünk a dobásunkról. Mit látunk a képen? Az eldobott tárgyat abban a távolsági és magassági helyzetben, ahol a fotózás pillanatában látszott. Ha elég kis záridővel dolgoztunk, a tárgy szép élesen látszik és meg tudjuk mérni, hogy hol volt pontosan az adott pillanatban. Ez a fénykép a State Mátrix, vagyis a rendszer pillanatképe, az adott t időpillanatban. Megtehetnénk, hogy minden egyes pixelt megőrzünk a képről, de mivel takarékos és okos emberek vagyunk, rájövünk, hogy teljesen fölösleges számon tartani minden egyes helyet, ahol a tárgy nem volt az adott időpillanatban. Elég ha csak azt tartjuk számon ahol volt, így jelentősen csökkentjük a képünk nagyságát. 

Mi tehetünk, ha tudni szeretnénk, hogy egy bizonyos idő eltelte után – nevezzük ezt az időt Δt-nek – hogy fog kinézni ez a fénykép? Két lehetőség van: készíthetünk egy másik képet, vagy ha ismerjük a rendszer működését, megrajzolhatjuk magunk is. Szerencsére ez esetben ismerjük a rendszer működését, ezért a második lehetőséget választjuk. Miért ? Mert semmi se gátol meg minket abban, hogy az így rajzolt képet ne módosítsuk újra, és rajzoljunk egy másikat, amely azt ábrázolja, hogy 2 x Δt időtávolságra mit láthatunk. Aztán ezt az új ábrát megint módosíthatjuk … És mindehhez csak egyszer kell kiszámolni a Δt-hez kapcsolódó változások értékét, utána már csak ugyanezt a változást alkalmazzuk minden egyes képen. 

Amint az látható a state-space modell kulcsa a Δt időegység alatti rendszer változása. A ferde hajítás esetén ez pedig így írható fel:

(3)  x_t = x_{t-1} + \dot{x}_{t-1} \cdot \Delta t

(4)  z_t = z_{t-1} + \dot{z}_{t-1} \cdot \Delta t  - \frac{1}{2}\cdot g \cdot \Delta t^2

Ahol:

  • x  — a vízszintes távolság, ugyanúgy mint korábban
  • z  —  a magasság
  • \Delta t  —  az időkülönbség
  • \dot{x}  — a vízszintes sebesség
  • \dot{z}  —  a függőleges sebesség
  • g  —  a gravitációs állandó

A (3) és (4) és képletből láthatjuk, amikor leírjuk a rendszer pillanatnyi állapotát (elkészítjük a fényképünket) négy értéket kell megfigyelnünk: az x és z irányú pozíciókat és sebességeket.  A State Mátrix; szokásos jelölése X; a következőképpen fog kinézni:

(5)  X= \begin{bmatrix} x \\ z \\ \dot{x} \\ \dot{z} \end{bmatrix}

A következő lépésben a Δt-hez tartozó rendszer változást kell kifejeznünk.  Hogyan is néz ki egy state-space modell?  Az általános gyakorlat szerint a következő rendszer állapotot így szokták számolni:

(6)  X_{t} = A \cdot X_{t-1} + B \cdot U + w

Ez így elég homályos, ezért nézzük mit is jelentenek az egyes elemek:

  • X  —  a már ismert State Mátrix
  • U  — a Kontrol Vektor
  • A  —  Állapotátmeneti mátrix (State Transition Mátrix3)
  • B  — Kontrol Input Mátrix
  • w  —zaj, amitől most eltekinthetünk, értékét nullának vesszük.

Az A egy mátrix, ami leírja, hogy az X mátrix elemei milyen módon módosítják a következő állapot helyzetét. Esetünkben ez így fog kinézni:

(7)  A = \begin{bmatrix} 1  & 0 & \Delta t & 0 \\ 0 & 1 & 0 & \Delta t \\ 0 &0&1&0 \\ 0 &0&0&1 \\ \end{bmatrix}

Az U vektorba jönnek azok a változók, amelyek befolyásolják a következő állapotot, de nem szerepelnek az X mátrixunkban, ebben az esetben ez a gravitációs állandó:

(8)  U = \begin{bmatrix} g \end{bmatrix}

A B mátrix leírja, milyen módon kapcsolódnak az U mátrix és az X mátrix elemei. Ferde hajítás esetén a vízszintes pozíciót és a sebességet nem befolyásolják: így értékük 0. A magasságot viszont csökkenti \frac{1}{2} \cdot g \cdot  \Delta t^2 értékkel, a sebességet pedig g\cdot \Delta t -el. Így:

(9)  B = \begin{bmatrix} 0 \\ -\frac{1}{2} \cdot  \Delta t^2 \\ 0 \\ -\Delta t \end{bmatrix}

Ezzel kész is. Ezt így lehetne Pythonban megvalósítani:

# kezdő sebbeség
v = 300

# pozíciók az adott időpillanádban
x_t = []
y_t = []

# mérési időpillanatok
t = np.linspace(0, 70, num=1000) # Set time as 'continous' parameter.
t = t.tolist()
delta_t = t[1]-t[0]

# állandó mátrixok
A = np.array( [[1,0,delta_t,0], 
            [0,1,0,delta_t], 
            [0,0,1,0], 
            [0,0,0,1]] )
U = np.array([[g]])
B = np.array(  [ [0], [-0.5* (delta_t**2) ], [0],[-1*delta_t] ] )

# kezdő pont
X = np.array( [ [0], [0], [ v * np.cos(kezdoszog) ], [ v * np.sin(kezdoszog) ]  ] )

# az adott időpillanathoz tartozó pozíció számítása
for k in t:
    X = A.dot(X) + B.dot(U)
    
    # mentjük a számított pozíciót
    x_t.append( X[0][0] )
    y_t.append( X[1][0] )

# szűrjük azokat az időpillanatokat amelyek 0 magasság alatt vannak 
p = [i for i, j in enumerate(y_t) if j < 0]                         
for i in sorted(p, reverse = True):
    del x_t[i]
    del y_t[i]
    del t[i]

Nagyon jó! De mire jó ez? Ki tudjuk másképp számolni ugyanazt, és? Vegyük észre, hogy a két számítási mód abban tér el egymástól, hogy mi szükséges a pálya meghatározásához. Míg az első esetben a kezdő sebesség és a kilövési szög , addig a state-space modell esetén egy bármilyen időpillanatban készített rendszerállapot szükségleltetik. Miért érdekes ez? Azért mert a második számítási módhoz könnyebb beszerezni az adatokat.

A ferde hajítás légellenállással

Most nézzünk egy kicsit bonyolultabb esetet. Hogyan alakul a state-space modell, ha van légellenállás? Ekkor az adott magasságváltozást így számolhatjuk:

(10)  z_t = z_{t-1} + \dot{z}_{t-1} \cdot \Delta t  - \frac{1}{2}\cdot \ddot{z} \cdot \Delta t^2

Csupán kicseréltük a gravitációs együtthatót,  a z irányba mutató gyorsulásal4

A másik különbség, hogy ebben az esetben vízszintes irányban is hat a légellenállás, így azt az x számításába is be kell építenünk:

(11)  x_t = x_{t-1} + \dot{x}_{t-1} \cdot \Delta t  - \frac{1}{2}\cdot \ddot{x} \cdot \Delta t^2

Ez említett gyorsulások pedig:

(12)   \ddot{z} = -g-\frac{ \rho \cdot Cd \cdot A}{2\cdot m} \cdot \dot{z} \sqrt{ \dot{x}^2 +\dot{z}^2  }

(13)  \ddot{x} = -\frac{ \rho \cdot Cd \cdot A}{2\cdot m} \cdot \dot{x} \sqrt{ \dot{x}^2 +\dot{z}^2  }

Ahol

  • g  — a már ismert gravitációs együttható
  • \rho  —légsűrűség
  • Cd  —  az alakra jellemző légellenállási együttható
  • A  — a mozgás irányába mutató keresztmetszet területe
  • m  —a tárgy tömege

Innen már egyszerű összerakni a modellt. Az X és az A ugyanaz marad:

(14)     X= \begin{bmatrix} x \\ z \\ \dot{x} \\ \dot{z} \end{bmatrix}

(15)      A = \begin{bmatrix} 1  & 0 & \Delta t & 0 \\ 0 & 1 & 0 & \Delta t \\ 0 &0&1&0 \\ 0 &0&0&1 \\ \end{bmatrix}

Az U-nak most két eleme lesz, a függőleges és vízszintes gyorsulás. 

 (16)     U = \begin{bmatrix} \ddot{x} \\ \ddot{z} \end{bmatrix} =  \begin{bmatrix} -\frac{D}{m}\cdot \dot{x} \cdot \sqrt{ \dot{x}^2 + \dot{z}^2 } \\ -g-\frac{D}{m}\cdot \dot{z} \cdot \sqrt{ \dot{x}^2 + \dot{z}^2 } \\ \end{bmatrix}

A B-ben pedig összekapcsoljuk a megfelelő dimenziókat.

(17)     B =  \begin{bmatrix} \frac{1}{2} \cdot \Delta t^2 & 0 \\ 0 &\frac{1}{2} \cdot \Delta t^2 \\ \Delta t & 0 \\ 0 & \Delta t \\ \end{bmatrix}

Ennek Python megvalósítása:

from math import sqrt, exp, radians, cos, sin, pi

# pozíciók az adott időpillanádban
x_t = []
z_t = []

# állandók
#   tömeg kg
m = 46
#   átmérő, m
d = .122
# alak ellenállása (Drag Coefficient)
Cd = .2
# légsürüsség kg/m^3
rho= 1.2
# gravitáció
g = 9.8
# kezdő sebbeség
v = 1500
# theta (kilövési szög)
theta = radians(50)

# Számított állandok
#   keresztmeszet területe
A = pi*(d**2)
# ellenállásra jelemző állandó 
D = rho*Cd*A/2

# mérési időpillanatok
t = np.linspace(0, 100, num=100000) # Set time as 'continous' parameter.
t = t.tolist()
delta_t = t[1]-t[0]


X = np.array( [ [0], [0], [ v * np.cos(theta) ], [ v * np.sin(theta) ]  ] )
A = np.array( [ [1,0,delta_t,0], [0,1,0,delta_t], [0,0,1,0], [0,0,0,1]] )
B = np.array( [ [ 0.5*delta_t**2, 0], [ 0, 0.5*delta_t**2], [delta_t, 0], [0, delta_t] ] )

for k in t:
    U = np.array([
        [ -1.*(D/m) * X[2][0] * sqrt( X[2][0]**2 + X[3][0]**2 ) ], 
        [ -g-(D/m) * X[3][0] * sqrt( X[2][0]**2 + X[3][0]**2 ) ]]
    )
    X = A.dot(X) + B.dot(U)
    
    # mentjük a szimulált poziciót
    x_t.append( X[0][0] )
    z_t.append( X[1][0] )

    
# szűrjük azokat az időpillanatokat amely 0 magasság alatt vannak
p = [i for i, j in enumerate(z_t) if j < 0]                          
for i in sorted(p, reverse = True):
    del x_t[i]
    del z_t[i]
    del t[i]

Ez a következő pályát rajzolja ki:



Lábjegyzet

  1. Tudja valaki a „state-space” magyar megfelelőjét? Köszönet Akos Gabrelnek ez is megvan: állapottér.
  2. A szokásos jelölés két dimenzióban y szokott lenni, de, hogy ne legyen bonyolult a modell későbbi kiterjesztése három dimenzióra ezért eltértem ettől.
  3. Erre mi a magyar megfelelő? Szintén Akos Gabriel nyomán: állapotátmeneti mátrix
  4. A két pont a z betű fölött a gyorsulást jelöli


Hírdetés

A ballisztikus pálya állapottér (state-space) modellje” bejegyzéshez 8 hozzászólás

Vélemény, hozzászólás?

Adatok megadása vagy bejelentkezés valamelyik ikonnal:

WordPress.com Logo

Hozzászólhat a WordPress.com felhasználói fiók használatával. Kilépés /  Módosítás )

Twitter kép

Hozzászólhat a Twitter felhasználói fiók használatával. Kilépés /  Módosítás )

Facebook kép

Hozzászólhat a Facebook felhasználói fiók használatával. Kilépés /  Módosítás )

Kapcsolódás: %s