---
# System prepended metadata

title: Projekt PPSK - magnetyczna stabilizacja satelity

---

# Projekt PPSK - magnetyczna stabilizacja satelity
*Krzysztof Haładyn, Michał Szaknis*

[toc]

## Wstęp
W ramach projektu zajmowaliśmy się analizą problemu stabilizacji satelity na orbicie okołoziemskiej z wykorzystaniem do tego celu jej pola magnetycznego. Na satelicie umieszczany jest zestaw 3 cewek na ortogonalnych ściankach, których odpowiednie zasilanie pozwala stabilizować obrót satelity w przestrzeni trójwymiarowej.

W ramach projektu przygotowaliśmy symulację, pokazującą jak satelita wyposażony w taki system zachowuje się w przestrzeni kosmicznej. Ze względu na ograniczenia czasowe, ograniczyliśmy się do rozpatrzenia problemu dla układu z jedną cewką.

Symulacja została napisana w języku Python. Dodatkowo, przygotowaliśmy jej wizualizację w Unity.

## Symulacja zachowania na orbicie

Poniższe wzory bazowane są na https://louis.uah.edu/cgi/viewcontent.cgi?article=1219&context=honors-capstones

### Ruch bryły sztywnej
$\vec{\tau} = \matrix{I} \cdot \vec{\alpha} = \matrix{I} \cdot \ddot{\theta}$
**Gdzie**:
* $\vec{\alpha}$ - Wektor przyspieszenia kątowego
* $\matrix{I}$ - Macierz momentu bezwładności satelity

Jednakże powyższy wzór działa tylko w inercyjnym układzie odniesienia. Dla układów odniesienia z występującą rotacją nalezy zastosować dodatkowo twierdzenie tzw. "rate of change transport theorem". Niestety nie znaleźliśmy jego polskiej nazwy. Zgodnie z tą zasadą do wyliczenia pochodnych wektorów w takim układzie odniesienia nalezy zastosować operator:

$\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{f} = \left[ \left(\frac{\mathrm{d}}{\mathrm{d}t}\right)_{\mathrm{r}} + \boldsymbol{\Omega} \times \right] \boldsymbol{f}$

Zatem:

$\matrix{I} \cdot \vec{\alpha} = \matrix{I} \cdot \dot{\omega} = \vec{\tau} - \vec{\omega} \times \matrix{I} \vec{\omega}$

Po przekształceniu:

$\vec{\alpha} = \ddot{\vec{\theta}} = \matrix{I}^{-1} \cdot (\vec{\tau} - \vec{\omega} \times \matrix{I} \vec{\omega})$

### Cewka magnetyczna

$\vec{\tau} = \vec{\mu} \times \vec{B}$
**Gdzie**:
* $\vec{B}$ - Zewnętrzny wektor indukcji pola magnetycznego (pole magnetyczne Ziemi)
* $\vec{\mu}$ - Moment magnetyczny cewki stabilizującej (na satelicie)

$\vec{\mu} = i \cdot \vec{A} = NiS \cdot \vec{n}$
**Gdzie**:
* $\vec{A}$ - Wektor powierzchni cewki
* $\vec{n}$ - Wektor normalny powierzchni cewki
* $S$ - Powierzchnia cewki

### Reprezentacja obrotu

Do reprezentacji obrotów wykorzystywane są kwaterniony. Zatem musimy wyznaczyć pochodną kwaterniona obrotu naszego satelity po czasie.

Mamy kwaternion $q = \{a, b, c, d\}$, w zapisie wykładniczym $q = e^{\frac{\theta}{2}(b\boldsymbol{i} + c\boldsymbol{j} + d\boldsymbol{k})}$.
Podnoszenie kwaterniona obrotu do potęgi można interpretować jako składanie obrotów. Podobnie jak w przypadku jednostkowej liczby zespolonej jej potęgowanie jest obrotem przeciwnym do ruchu wskazówek zegara wzgledem początku układu współrzednych.

$q^t = e^{t\frac{\theta}{2}(b\boldsymbol{i} + c\boldsymbol{j} + d\boldsymbol{k})}$

Zatem obrót kwaterniona w czasie:

$q(t) = q_0e^{t\frac{\theta}{2}\vec{a}}$
**Gdzie**:
* $\vec{a} = b\boldsymbol{i} + c\boldsymbol{j} + d\boldsymbol{k}$
* $q_0$ - początkowy obrót

Powyższe pozwala nam wyznaczyć pochodną $q$ po $t$.

$\frac{dq}{dt} = \frac{d q_0e^{t\frac{\theta}{2}\vec{a}}}{dt} = \frac{\theta}{2}\vec{a} \cdot q_0e^{t\frac{\theta}{2}\vec{a}} = \frac{\theta}{2}\vec{a} \cdot q(t)$

Pełna rotacja o kąt $\theta$ względem osi $\vec{a}$ w czasie jednostkowym jest prędkością kątową $\omega$. Zatem:

$\theta \vec{a} = \vec{\omega}$

Ostatecznie

$\frac{dq(t)}{dt} = \frac{1}{2}\omega q(t)$

### Podsumowanie - równanie stanu

Stanem naszego satelity w symulacji jest:
* $\vec{r}$ - Wektor pozycja względem Ziemi
* $\vec{v}$ - Wektor prędkości orbitalnej względem Ziemi
* $\vec{\omega}$ - Wektor prędkości kątowej satelity względem Ziemi
* $q$ - Kwaternion obrotu satelity śledzący obrót względem Ziemi

Stan jest aktualizowany poprzez numeryczne rozwiązanie następujących równaniń różniczkowych:

$\left\{ \begin{aligned} 
& \vec{r} \\ 
& \vec{v} \\ 
& \vec{\omega} \\ 
& q \\
\end{aligned} \right\} \left\{ \begin{aligned} 
& \dot{\vec{r}} = \vec{v} \\ 
& \dot{\vec{v}} = \frac{GM}{||\vec{r}^3||}\vec{r} \\ 
& \dot{\vec{\omega}} =  \matrix{I}^{-1} \cdot (\vec{\tau} - \vec{\omega} \times \matrix{I} \vec{\omega}) \\
& \dot{q} = \frac{1}{2}\vec{\omega} q \\
\end{aligned} \right\}$

### Układy współrzędnych

Wszystkie powyższe rozważania zapisane są przy założeniu układu współrzędnych ECI (Earth-centered inertial). Jest to układ współrzędnych związany z Ziemią (umieszczamy ją na 0,0,0), w którym Ziemia nadal wykonuje obrót wokół własnej osi. Bardzo wygidnie opisuje się w tym układzie równania ruchu, ale należy zwrócić szczególną uwagę na rozróżnienie między ECI a ECEF (Earth-centered, Earth-fixed), związanym z nieobracającą się Ziemią. Czasami będziemy potrzebowali przejść między jednym układem a drugim, np. w celu wyznaczania współrzędnych geograficznych.

<!-- złe osie xd ![](https://upload.wikimedia.org/wikipedia/commons/3/32/Earth_Centered_Inertial_Coordinate_System.png) -->

## Modele magnetyczne
Pole magnetyczne Ziemii jest bardzo niejednorodne i zmienia się w czasie, zatem kluczowe dla symulacji jest dobranie jego odpowiedniego modelu. Nasz symulator obecnie obsługuje dwa modele magnetyczne:
* prosty model dipolowy
* WMM2020, z wykorzystaniem implementacji z biblioteki https://github.com/space-physics/wmm2020

### Prosty model dipolowy
Najprostszym, lecz niezbyt rzeczywistym modelem pola magnetycznego jest model dipolowy

![](https://i.imgur.com/R3kmn7s.png)

Do jego implementacji wykorzystaliśmy bibliotekę `magpylib`, zakładając poziom namagnesowania około 40000 nT.

### Model WMM2020
World Magnetic Model 2020 jest aproksymacją Ziemskiego pola magnetycznego. Działa poprzez wykorzystanie tzw. sferycznych harmonicznych, czyli szeregu funkcji zdefiniowanych na powierzchiach sferycznych tworzących konstrukcje działającą na podobnej zasadzie co szereg Fouriera. Model do dokładnego opisu ziemskiej magnetosfery wykorzystuje pierwsze 12 harmonicznych. Pozwala to na już dość efektywne uwidocznienie niejednorodności i kształtu linii pola.

Ponieważ pole magnetyczne Ziemi zmienia się w czasie, co kilka lat wydawane są nowe wersje tego modelu. WMM2020 ważny jest na lata 2020-2025.

![](https://i.imgur.com/pCL4A29.jpg)

Jako wejścia, model ten przyjmuje:
* długość i szerokość geograficzną punktu
* wysokość nad poziomem morza
* czas, podany w formacie `decimal year` - rok z ułamkiem

Na wyjściu, dostajemy (między innymi) wektor pola magnetycznego w układzie NED (north-east-down) w nanoteslach (nT).

Aby skorzystać z tego modelu, najpierw musimy wyznaczyć długość i szerokość punktu, nad którym znajduje się satelita. W tym celu, najpierw konwertujemy współrzędne satelity z ECI do ECEF, a następnie stosujemy powszechnie znaną konwersję

$R = |Position|$
$g_{lat} = arctan2(-Position[y], R)$
$g_{lon} = arctan2(Position[x], Position[z])$

<!---
:::warning
symulatorowi brakuje ECI->ECEF XD
:::
-->

Wyjście modelu musimy przekonwertować z powrotem na globalny układ współrzędnych naszej symulacji. Ponieważ mieliśmy już zdecydowanie dosyć matematyki, użyliśmy tu gotowej funkcji `enu2uvw` z biblioteki `pymap3d`, która konwertuje wektor wyrażony w układzie ENU (East-North-Up, analogicznym do NED) na zwykłe XYZ w układzie ECEF. Operacja ta wymaga oczywiście podania położenia obserwatora, w formacie długość/szerokość/wysokość n.p.m. (uwaga: jest to inna operacja niż `enu2ecef` - ta zwraca pozycję bezwzględną, więc nie nadaje się do konwersji wektora wyrażonego w nT)

<!---
:::warning
Nie powinnięmy tu czasem jeszcze zrobić ECEF->ECI na koniec...?
:::
-->

![](https://upload.wikimedia.org/wikipedia/commons/7/73/ECEF_ENU_Longitude_Latitude_relationships.svg)

## Kontroler PID
Gdybyśmy nie modulowali w żaden sposób prądu przesyłanego przez cewkę, otrzymalibyśmy oscylator harmoniczny - dlatego prądem tym trzeba odpowiednio sterować. W ramach symulacji zaimplementowany został prosty kontroler PID, który stara się wyrównać satelitę tak, aby jej przód ustawił się w kierunku wektora magnetycznego. Niestety, brakło nam już czasu na odpowiednie dobranie parametrów kontrolera.

### Wyznaczanie odchyłu
Do wyznaczenia odchyłu mierzymy kąt między wektorem indukcji magnetycznej, a wektorem momentu magnetycznego cewki. Niestety, tradycyjne meteody mierzenia kątów w przestrzenicach opierające się na obliczaniu iloczynu wektorowego lub skalarnego nie daje informacji o kierunku odchylenia. Tak uzyskane wartości są zawsze dodatnie przez co całkująca część regulatora PID staje się rozbieżna co uniemożliwia sterowanie. Do rozwiązania tego problemu staraliśmy się znaleźc metodę wyznaczania kąta ze znakiem. Znalezione przez nas metody w literaturze polegały na wyznaczaniu kąta względem pewnej arbitralnej płaszczyzny obserwatora. Niestety, w naszym problemie symetria osiowa uniemożliwia wyznaczenia takiej płaszczyzny obserwacji. Do obejścia tego problemu zastosowaliśmy więc rozwiązanie "studenckie" polegające na mierzeniu pochodnej zmiany kąta i zmianie jego znaku przy zmianie znaku pochodnej. Taka zmiana znaku występuje przy spotkaniu się wektora Ziemskiej induckji magnetycznej z momentem cewki w satelicie. Całość pozwoliła już na uruchomienie regulatora PID.

### Sterowanie cewką
Pozostaje jeszcze problem sterowania prądem w cewkach na podstawie wyjscia regulatora PID. W naszym przypadku wektor momentu magnetycznego zawsze będzie generował taki moment siły, aby ustawić się równolegle do wektora indukcji Ziemskiej magnetosfery. Wyjście regulatora PID interpretuje się natomiast jako kierunek zmiany odchyłu i jego amplitudę. Z tego powodu nie możemy po prostu zinterpretować wyjścia regulatora jako prądu. Spowodowało by to odwrócenie kierunku wektora momentu magnetycznego cewki, czego efektem będzie dalsze rozpędzanie się satelity po przejściu odchyłu przez 0. W celu rozwiązania tego problemu wraz ze zmianą znaku odchyłu zmieniamy znak wyjścia z regulatora. Dzięki temu prąd w cewce zawsze płynie w odpowiednią stronę.

### Strojenie
Istnieje wiele metod strojenia regulatorów PID, ktore możemy podzielić na:
* metody manualne - polegające na ręcznym wyznaczeniu kazdego z 3 parametrów obserwując zachowanie się sterowanego obiektu,
* metody półautomatyczne lub automatyczne - polgają na wstępnym wyznaczeniu parametrów na podstawie zadanego algorytmu, a następnie ich dostrojeniu.

Niestety, żadne z metod nie zadziałała w naszym przypadku. Ręczne strojenie okazało się bardzo czasochłonne i nie dawało żadnych, nawet częściowych rezultatów. W przypadku metod półautomatycznych np. metody Zieglera-Nicholsa nie byliśmy w stanie poprawnie wyznaczyć potrzebnych danych do algorytmu. W obu przypadkach narpiew należało doprowadzić układ na granicę stabilności regulacji co w naszym przypadku okazało się niemożliwe ze względu na brak oporów ruchu w naszej symulacji. W przypadku układów na Ziemi zawsze występują opory ruchu co powoduje, że nawet oscylatory harmoniczne po pewnym czasie stabilizują się w pewnym położeniu poprzez utratę energii drgań. W przypadku satelity na oribcie taka utrata energi albo nie występuje albo jest znikoma. Efektem tego są niegasnące drgania, których nie udało się nam wygasić używając komponentu proporcjonalnego regulatora PID, co jest konieczne do zastosowania tej metody. 

## Opis interfejsu symulatora

Sam symulator składa się z dwóch komponentów - symulatora właściwego, napisanego w Pythonie oraz opcjonalnej wizualizacji 3D, zrealizowanej w Unity. Symulator właściwy otwiera okno pozwalające wybrać i dostosowywać w czasie rzeczywistym parametry symulacji:
* wykorzystywany model magnetyczny
* parametry kontrolera PID
* datę i godzinę symulacji (istotna dla modelu WMM2020)
* zwiększenie predkości symulacji (umożliwia oglądanie symulacji z prędkością większą niż rzeczywista)

Dodatkowo, widzimy wykres odchylenia i podawanego sterowania (wejścia i wyjścia) kontrolera PID w czasie.

Brakło nam czasu na dorobienie wyboru początkowych parametrów orbity (trzeba je ustawić ręcznie w kodzie, podając początkową pozycję i prędkość jako wektory 3D). Nie da się też używając GUI zmienić masy satelity (1kg), rozmiaru (10cm x 10cm x 10cm), parametrów cewki (100 zwojów, prąd max 1A) i jej orientacji (+X w układzie współrzędnych satelity).

## Screenshoty

<!-- to jest to samo co niżej tylko nieustawione, meh
![](https://i.imgur.com/hBliidl.png)-->
![](https://i.imgur.com/NbTglnc.png)
![](https://i.imgur.com/OkPB0f5.png)
![](https://i.imgur.com/64aXuZc.png)

## Podsumowanie

Udało nam się napisać w pełni działający symulator zachowania się satelity na orbicie, wraz z symulacją zachowania cewki w polu magnetycznym Ziemi. Problemem okazało się być zestrojenie kontrolera, tak, aby odpowiednio sterować satelitą z jej wykorzystaniem. Uzyskane rezultaty pokazują jednak, że wykonanie takiego sterowania w celu stabilizacji satelity w polu magnetycznym Ziemi jest jak najbardziej możliwe.
