4.6 Stabilność Lapunowa i sterowanie nieliniowe · Matematyka i fizyka od podstaw

0 Prolog: stabilność bez rozwiązywania

W 1892 r. Aleksandr Lapunow obronił w Moskwie rozprawę „Ogólny problem stateczności ruchu”. Wyszedł od obserwacji z mechaniki: jeśli energia układu ciągle maleje, a jest najmniejsza w punkcie równowagi, to układ musi do tej równowagi zmierzać, jakkolwiek skomplikowane byłyby jego równania. Nie trzeba znać trajektorii, wystarczy znaleźć funkcję podobną do energii i sprawdzić, że maleje.

Przez pół wieku praca była prawie nieznana na Zachodzie. Gdy w latach 50. i 60. automatyka zderzyła się z nieliniowymi problemami rakiet i statków kosmicznych, metoda Lapunowa stała się podstawą całej teorii sterowania nieliniowego.

Pytanie na start. Wahadło z tarciem puszczone z kąta 143°. Wiadomo, że kiedyś się zatrzyma na dole. Czy da się to udowodnić bez rozwiązywania równania $\ddot\theta = -\frac gL\sin\theta - c\dot\theta$ , którego nie da się rozwiązać w funkcjach elementarnych? Odpowiedź w części 1.

1 Funkcja Lapunowa

Twierdzenie Lapunowa (wersja praktyczna)

Punkt równowagi $\mathbf x = 0$ układu $\dot{\mathbf x} = f(\mathbf x)$ . Jeśli istnieje funkcja $V(\mathbf x)$ taka, że $V(0) = 0$ , $V(\mathbf x) > 0$ dla $\mathbf x \ne 0$ (w pewnym otoczeniu) oraz

\dot V = \nabla V\cdot f(\mathbf x) < 0\quad\text{dla } \mathbf x \ne 0,

to równowaga jest asymptotycznie stabilna. Gdy $\dot V \le 0$ , jest stabilna (trajektorie nie uciekają).

Klucz: $\dot V$ liczymy z równania ruchu, bez jego rozwiązywania. Odpowiedź na pytanie z prologu. Dla wahadła z tarciem weźmy energię (na jednostkę $mL^2$ ): $V = \frac12\omega^2 + \frac gL(1 - \cos\theta)$ . Wtedy

\dot V = \omega\dot\omega + \frac gL\sin\theta\,\omega = \omega\left(-\frac gL\sin\theta - c\omega\right) + \frac gL\omega\sin\theta = -c\,\omega^2 \le 0 .

Energia nigdy nie rośnie, więc wahadło nie ucieknie z obszaru $V < V(0)$ . $\dot V$ jest zerowe tylko chwilowo (gdy $\omega = 0$ w punkcie zwrotnym), a jedynym miejscem, gdzie wahadło może zostać na stałe z $\omega = 0$ , jest równowaga. Zasada LaSalle’a formalizuje ten krok i daje asymptotyczną stabilność. Z kąta 143° ( $V = 17{,}7$ ) wahadło nie przejdzie przez górę, bo to wymaga $V > 2g/L = 19{,}6$ : zatrzyma się na dole.

Zadania · część 1

3 zadania · ok. 20 min

Dla $\dot x = -x^3$ pokaż stabilność funkcją $V = x^2$ . Dlaczego linearyzacja tu nic nie mówi?

Odpowiedź
$\dot V = 2x\dot x = -2x^4 < 0$ dla $x \ne 0$ : asymptotycznie stabilny. Linearyzacja w zerze daje $\dot x = 0$ (biegun w zerze), czyli przypadek graniczny, z którego nic nie wynika.
Jaka prędkość początkowa na dole wystarcza wahadłu bez tarcia, żeby przejść przez górę ( $g/L = 9{,}81$ 1/s²)?

Odpowiedź
$\frac12\omega_0^2 > 2\frac gL$ : $\omega_0 > 2\sqrt{g/L} = 6{,}26$ rad/s. Z tarciem trzeba więcej (z $c = 0{,}5$ nawet 6,5 rad/s nie wystarczy).
Dla wahadła bez tarcia $\dot V = 0$ . Co twierdzenie mówi o stabilności dolnej równowagi? Co o górnej?

Odpowiedź
Dolna: $V > 0$ wokół niej i $\dot V = 0$ , więc stabilna, ale nie asymptotycznie (wahadło oscyluje bez końca). Górna: $V$ ma tam punkt siodłowy (rośnie z $\omega$ , maleje z odchyleniem kąta), a nie minimum, więc twierdzenie nie ma zastosowania, a równowaga jest niestabilna (lekcja math 4.6).

2 Lapunow dla układów liniowych

Dla $\dot{\mathbf x} = A\mathbf x$ naturalna jest funkcja kwadratowa $V = \mathbf x^TP\mathbf x$ z $P$ dodatnio określonym. Wtedy $\dot V = \mathbf x^T(A^TP + PA)\mathbf x$ . Wystarczy, żeby

A^TP + PA = -Q,\qquad Q \succ 0 .

To równanie Lapunowa: liniowe względem $P$ . Dla stabilnego $A$ ma ono zawsze jedno rozwiązanie, i to dodatnio określone; dla niestabilnego zwykle też ma rozwiązanie, ale nie dodatnio określone. Ta sama macierz pojawiła się w LQR: koszt pętli zamkniętej $\int\mathbf x^TQ\mathbf x\,\mathrm dt = \mathbf x_0^TP\mathbf x_0$ (lekcja 4.3 LQR: koszt zamiast biegunów). Funkcja Lapunowa jest tu kosztem „od teraz do nieskończoności”.

Praktyczne zastosowanie: funkcja Lapunowa z linearyzacji pozwala oszacować obszar przyciągania układu nieliniowego, czyli największą poziomicę $V$ , wewnątrz której $\dot V < 0$ także dla pełnego, nieliniowego modelu.

Zadania · część 2

2 zadania · ok. 20 min

Pętla PD wysokości: $A = \begin{bmatrix}0 & 1\\ -4 & -2{,}8\end{bmatrix}$ . Rozwiąż $A^TP + PA = -I$ .

Odpowiedź
Z trzech równań na $p_{11}$ , $p_{12}$ , $p_{22}$ : $-8p_{12} = -1$ , $p_{11} - 2{,}8p_{12} - 4p_{22} = 0$ , $2p_{12} - 5{,}6p_{22} = -1$ . Stąd $p_{12} = 0{,}125$ , $p_{22} = 0{,}223$ , $p_{11} = 1{,}243$ . Wartości własne 1,258 i 0,208: dodatnio określona ✓.
Dlaczego dla niestabilnego $A$ równanie Lapunowa nie może mieć rozwiązania $P \succ 0$ ?

Odpowiedź
Gdyby miało, $V = \mathbf x^TP\mathbf x$ byłoby dodatnie i malejące wzdłuż każdej trajektorii, więc wszystkie trajektorie zmierzałyby do zera. To przeczy istnieniu rosnącego modu własnego.

3 Linearyzacja przez sprzężenie

Linearyzacja w punkcie pracy przybliża model. Linearyzacja przez sprzężenie zwrotne (ang. feedback linearization) robi coś mocniejszego: sterowaniem dokładnie kasuje nieliniowość, tak że z punktu widzenia nowego wejścia $v$ obiekt jest liniowy w całym zakresie pracy.

Przykład z każdego kontrolera lotu: dron przechylony o kąt $\varphi$ ma pionową składową ciągu $T\cos\varphi$ . Pionowo: $m\ddot z = T\cos\varphi - mg$ . Jeśli regulator wysokości liczy żądane przyspieszenie $v$ (np. PID z lekcji 3.1), to ciąg

T = \frac{m\,(g + v)}{\cos\varphi}

daje dokładnie $\ddot z = v$ przy każdym przechyleniu. Bez tego członu przechylony dron opada: przy 30° pionowe przyspieszenie spada o $g(1 - \cos 30^\circ) = 1{,}31$ m/s², a regulator dowiaduje się o tym dopiero z uchybu.

Warunki i ograniczenia: trzeba znać model (masa, $\cos\varphi$ ), a podział przez $\cos\varphi$ wymaga $|\varphi| < 90^\circ$ i zapasu ciągu (przy 45° potrzeba 1,41 razy więcej). Przy niedokładnym modelu kasowanie jest niedokładne, a resztę musi poprawić pętla.

Zadania · część 3

3 zadania · ok. 20 min

Ile razy większy ciąg potrzebny jest do zawisu przy przechyleniu 10°, 30°, 45°?

Odpowiedź
$1/\cos\varphi$ : 1,015; 1,155; 1,414.
Wahadło sterowane momentem: $mL^2\ddot\theta = -mgL\sin\theta + \tau$ . Zaproponuj $\tau$ , które daje $\ddot\theta = v$ .

Odpowiedź
$\tau = mgL\sin\theta + mL^2v$ . Z $v = -k_1\theta - k_2\dot\theta$ każda pozycja (także odwrócona) jest stabilizowana jak liniowy oscylator, w całym zakresie kątów (o ile moment wystarcza).
Dlaczego linearyzacja przez sprzężenie jest ryzykowna, gdy obiekt ma zero w prawej półpłaszczyźnie (lekcja 2.2)?

Odpowiedź
Kasowanie dynamiki to w istocie odwrócenie modelu. Odwrócenie prawego zera daje niestabilny biegun w regulatorze: wewnętrzna dynamika (ang. zero dynamics) eksploduje, choć wyjście wygląda dobrze. Metoda wymaga, żeby ta wewnętrzna dynamika była stabilna.

4 Sterowanie ślizgowe

Co, jeśli modelu nie znamy dokładnie, ale wiemy, że zakłócenie jest ograniczone: $|d| \le D$ ? Sterowanie ślizgowe (ang. sliding mode control) robi dwie rzeczy. Wybiera powierzchnię ślizgu, na której ruch jest taki, jak chcemy, np. dla $\ddot x = u + d$ :

s = \dot x + \lambda x = 0\quad\Rightarrow\quad x(t) = x(t_r)\,e^{-\lambda(t - t_r)}\quad(t \ge t_r,\ t_r\text{ — chwila dojścia do powierzchni}) .

Potem buduje sterowanie, które sprowadza stan na tę powierzchnię i tam trzyma, cokolwiek robi zakłócenie: $u = -\lambda\dot x - k\operatorname{sign}(s)$ z $k > D$ . Uzasadnienie przez Lapunowa: dla $V = \frac12s^2$ jest $\dot V = s\dot s = s(d - k\operatorname{sign}s) \le -(k - D)|s| < 0$ . Stąd też oszacowanie czasu dojścia: co najwyżej $|s_0|/(k - D)$ .

Na powierzchni sterowanie przełącza się nieskończenie szybko między $\pm k$ i w średniej dokładnie kasuje zakłócenie. W praktyce przełączanie z częstotliwością próbkowania to drgania sterowania (ang. chattering): niszczą siłowniki i pobudzają niemodelowane rezonanse. Typowe lekarstwo to warstwa przejściowa: $\operatorname{sign}(s) \to \operatorname{sat}(s/\varepsilon)$ . Drgania znikają, a błąd zostaje ograniczony, ale już nie zerowy (rzędu $\varepsilon/\lambda$ ).

Zadania · część 4

3 zadania · ok. 25 min

Start $x = 1$ , $\dot x = 0$ , $\lambda = 2$ , $k = 1{,}5$ , $D = 0{,}8$ . Oszacuj czas dojścia do powierzchni ślizgu.

Odpowiedź
$s_0 = 0 + 2\cdot 1 = 2$ ; $t \le 2/(1{,}5 - 0{,}8) = 2{,}86$ s. Symulacja: ok. 1,8 s (oszacowanie jest pesymistyczne, bo zakłócenie nie zawsze działa przeciwko).
PD $u = -4x - 2{,}8\dot x$ z tym samym zakłóceniem: jaka amplituda drgań ustalonych? (Policz $\left|\frac{1}{(j\omega)^2 + 2{,}8j\omega + 4}\right|$ dla $\omega = 2$ .)

Odpowiedź
$(j2)^2 + 4 = 0$ , zostaje $|2{,}8\cdot 2j| = 5{,}6$ ; amplituda $0{,}8/5{,}6 = 0{,}143$ . Ślizgowy z $\operatorname{sign}$ : ok. 0,0006, z warstwą $\varepsilon = 0{,}05$ : ok. 0,009.
Dlaczego ruch na powierzchni ślizgu nie zależy od zakłócenia, a ruch PD tak?

Odpowiedź
Na powierzchni obowiązuje $s = 0$ , czyli $\dot x = -\lambda x$ : równanie pierwszego rzędu, w którym $d$ nie występuje, bo sterowanie przełączające w średniej dokładnie je kompensuje. PD ma skończone wzmocnienie i zakłócenie przechodzi przez $S$ (lekcja 2.5).

5 Metoda i przykłady rozwiązane

Metoda: kiedy liniowe już nie wystarcza

Sprawdź zakres: jeśli punkt pracy się zmienia, ale powoli i mierzalnie, wystarczy gain scheduling (lekcja 3.3).
Znana nieliniowość (grawitacja, $\cos\varphi$ , aerodynamika z modelu): linearyzacja przez sprzężenie plus liniowa pętla na resztę.
Nieznane, ale ograniczone zakłócenia: sterowanie ślizgowe z warstwą przejściową.
Dowód stabilności: funkcja Lapunowa (energia, $\frac12s^2$ , $\mathbf x^TP\mathbf x$ ).
Sprawdzenie: symulacja nieliniowa z błędami modelu, opóźnieniem i próbkowaniem.

Przykład 5.1 Wahadło z 143°

$g/L = 9{,}81$ 1/s², $c = 0{,}5$ 1/s, start $\theta_0 = 2{,}5$ rad, $\omega_0 = 0$ . Co wiemy bez rozwiązywania?

$V(0) = 9{,}81(1 - \cos 2{,}5) = 17{,}67 < 2g/L = 19{,}62$ : wahadło nie przejdzie przez górę.
$\dot V = -0{,}5\omega^2 \le 0$ , a LaSalle daje zbieżność do $\theta = 0$ .
Symulacja: V w kolejnych sekundach 17,7; 11,5; 8,0; 5,2; 3,2; 2,0, po 20 s prawie zero.

Wynik: wahadło zatrzyma się na dole, bez żadnej trajektorii w ręku.

Sprawdzenie: symulacja 1 z tymi ustawieniami ✓.

Przykład 5.2 Kompensacja przechylenia w dronie

Dron 1 kg leci w bok z przechyleniem 30°, regulator wysokości żąda $v = 0$ . Porównaj ciąg $mg$ i ciąg z linearyzacji przez sprzężenie.

Ciąg $mg = 9{,}81$ N: pionowa składowa $9{,}81\cos 30^\circ = 8{,}50$ N, dron przyspiesza w dół 1,31 m/s².
Ciąg $mg/\cos 30^\circ = 11{,}33$ N: pionowa składowa dokładnie 9,81 N, wysokość stała.

Wynik: jedna linijka (dzielenie przez $\cos\varphi$ ) usuwa zakłócenie, które regulator liniowy musiałby zbijać uchybem.

Sprawdzenie: przy $\varphi = 0$ wzór daje zwykły ciąg zawisu ✓.

6 Błędne intuicje, których trzeba się oduczyć

„Jeśli nie umiem znaleźć funkcji Lapunowa, układ jest niestabilny.”: Twierdzenie daje warunek wystarczający. Brak funkcji to brak dowodu, nie dowód niestabilności.
„Linearyzacja w punkcie pracy i linearyzacja przez sprzężenie to to samo.”: Pierwsza przybliża model, druga kasuje nieliniowość sterowaniem i działa w całym zakresie, o ile model jest dokładny.
„Sterowanie ślizgowe jest odporne na wszystko.”: Tylko na zakłócenia dopasowane do wejścia (wchodzące tam, gdzie sterowanie) i ograniczone przez $k$ . Płaci się drganiami sterowania albo warstwą z niezerowym błędem.
„Układ nieliniowy zawsze wymaga regulatora nieliniowego.”: Większość dronów i samolotów lata na regulatorach liniowych z harmonogramem nastaw. Nieliniowe metody stosuje się tam, gdzie te nie wystarczają.

7 Zadania do przemyślenia

Te zadania wymagają więcej myślenia niż liczenia.

Zadania do przemyślenia

4 zadania · ok. 1 godzina

Rozhuśtanie wahadła. Chcemy rozhuśtać wahadło z dołu do góry ograniczonym momentem. Zaproponuj prawo sterowania oparte na energii. Wskazówka: $\dot V = \omega\tau/(mL^2)$ (bez tarcia).
Rozwiązanie
1. Cel: $V = V_\text{góra} = 2g/L$ . Sterowanie $\tau = k\,\operatorname{sign}(\omega)\,(V_\text{góra} - V)$ zwiększa energię, gdy jej brakuje (pchamy zgodnie z ruchem), i zmniejsza, gdy jest za duża.
2. Po dojściu w okolicę góry przełączamy na regulator liniowy (LQR) stabilizujący odwrócone wahadło. To klasyczne rozwiązanie Åströma i Furuty (2000).
Lapunow dla orientacji satelity. Satelita obraca się z prędkością $\boldsymbol\omega$ ; energia kinetyczna obrotu $\frac12\boldsymbol\omega^TJ\boldsymbol\omega$ . Zaproponuj moment sterujący, który na pewno wytraca obrót.
Rozwiązanie
1. Z równań Eulera $\dot V = \boldsymbol\omega^T\boldsymbol\tau$ (człon żyroskopowy nie wykonuje pracy). Wybór $\boldsymbol\tau = -k\boldsymbol\omega$ daje $\dot V = -k|\boldsymbol\omega|^2 < 0$ .
2. To tłumienie prędkości obrotowej, używane do wytracania obrotu po odłączeniu od rakiety; jego przybliżona wersja dla cewek magnetycznych (magnetorquerów) to prawo „B-dot”. Temat kursu astronautyki (sterowanie orientacją).
Ślizg z próbkowaniem. Sterowanie ślizgowe liczone co $T = 1$ ms. Ile wynosi w przybliżeniu amplituda drgań $s$ wokół zera? (Między próbkami $\dot s \approx \pm k$ .)
Rozwiązanie
1. W jednym okresie $s$ zmienia się o ok. $kT$ , a znak przełącza się co próbkę: $|s| \lesssim kT$ (dla $k = 1{,}5$ : ok. 1,5 mm/s).
2. Stąd błąd $x$ rzędu $kT/\lambda$ . Wolniejsze próbkowanie wprost zwiększa drgania. Warstwa przejściowa powinna być co najmniej rzędu $kT$ .
Obszar przyciągania. Odwrócone wahadło z regulatorem liniowym zaprojektowanym w pobliżu pionu. Opisz, jak funkcją Lapunowa z równania Lapunowa oszacować, z jak dużego wychylenia regulator na pewno wróci.
Rozwiązanie
1. Z linearyzacji: $A_\text{cl}^TP + PA_\text{cl} = -Q$ . Dla pełnego modelu nieliniowego policz $\dot V = 2\mathbf x^TPf(\mathbf x)$ i znajdź największe $c$ , dla którego $\dot V < 0$ w całym obszarze $V(\mathbf x) < c$ (numerycznie, na siatce).
2. Obszar $V < c$ to gwarantowany (zwykle zachowawczy) obszar przyciągania. Wraz z ograniczeniem momentu daje odpowiedź „z ilu stopni regulator na pewno złapie wahadło”.

8 Sprawdź się

Odpowiedz w myślach, zanim rozwiniesz odpowiedź.

1 Co mówi twierdzenie Lapunowa?

Jeśli istnieje $V > 0$ z $\dot V < 0$ poza równowagą, równowaga jest asymptotycznie stabilna.

2 Jaka funkcja Lapunowa pasuje do wahadła z tarciem?

Energia: $V = \frac12\omega^2 + \frac gL(1 - \cos\theta)$ , $\dot V = -c\omega^2$ .

3 Co to równanie Lapunowa?

$A^TP + PA = -Q$ ; dodatnio określone $P$ istnieje wtedy i tylko wtedy, gdy $A$ jest stabilna.

4 Na czym polega linearyzacja przez sprzężenie?

Na skasowaniu znanej nieliniowości sterowaniem, np. $T = m(g + v)/\cos\varphi$ .

5 Jak działa sterowanie ślizgowe?

Sprowadza stan na powierzchnię $s = 0$ i trzyma go tam przełączaniem $\pm k$ z $k > |d|$ .

6 Co to chattering i jak go ograniczyć?

Szybkie przełączanie sterowania; warstwa przejściowa $\operatorname{sat}(s/\varepsilon)$ .

9 Karta wzorów i co dalej

Lapunow

$V > 0$ , $\dot V = \nabla V\cdot f < 0$ → asymptotycznie stabilny

wahadło: $V = \frac12\omega^2 + \frac gL(1 - \cos\theta)$ , $\dot V = -c\omega^2$

liniowy: $A^TP + PA = -Q$

Sterowanie nieliniowe

dron: $T = m(g + v)/\cos\varphi$

ślizg: $s = \dot x + \lambda x$ , $u = -\lambda\dot x - k\operatorname{sign}s$ , $k > D$

czas dojścia $\le |s_0|/(k - D)$ ; warstwa: błąd $\sim\varepsilon/\lambda$

Do przeczytania

Obowiązkowo (dla chętnych): J.-J. Slotine, W. Li, Applied Nonlinear Control, rozdz. 3 (Lapunow) i 7 (sterowanie ślizgowe), bardzo przystępne.
Dalej: H. K. Khalil, Nonlinear Systems (standardowy podręcznik); A. Tewari, Advanced Control of Aircraft, Spacecraft and Rockets, dodatek B (stabilność Lapunowa w zwięzłej formie).

To koniec części IV. Część V jest mostem do podręcznika Tewariego: sterowanie optymalne w ogólnej postaci (zasada maksimum Pontriagina, która daje LQR jako przypadek szczególny i sterowanie czasooptymalne jako inny) oraz pytanie, które wisi nad całą częścią IV: co się dzieje z odpornością, gdy LQR połączymy z filtrem Kalmana.