Dlaczego numeryczne rozwiązanie ODE odchodzi od niestabilnej równowagi?

15

Chciałbym zasymulować zachowanie systemu podwójnie wahadłowego. System jest manipulatorem robota o 2 stopniach swobody, który nie jest uruchamiany i dlatego będzie zachowywał się głównie jak podwójne wahadło pod wpływem grawitacji. Jedyną główną różnicą w przypadku podwójnego wahadła jest to, że składa się on z dwóch sztywnych ciał o właściwościach masy i bezwładności w ich środkach masy.

Zasadniczo programowałem ode45pod Matlabem, aby rozwiązać system ODE następujących typów:

[10000M110M1200100M120M22][x˙1x˙2)x˙3)x˙4]=[x2)-V.1-sol1x4-V.2)-sol2)]

gdzie x1 jest kątem pierwszego ciała w stosunku do poziomu, x2 jest prędkością kątową pierwszego ciała; x3 to kąt drugiego ciała w stosunku do pierwszego ciała, a x4 to prędkość kątowa drugiego ciała. Wszystkie te współczynniki są określone w kodeksie, przy użyciu rhsi fMassfunkcji stworzyłem.

clear all
opts= odeset('Mass',@fMass,'MStateDependence','strong','MassSingular','no','OutputFcn',@odeplot);
sol = ode45(@(t,x) rhs(t,x),[0 5],[pi/2 0 0 0],opts);

function F=rhs(t,x)
    m=[1 1];
    l=0.5;
    a=[0.25 0.25];
    g=9.81;
    c1=cos(x(1));
    s2=sin(x(3));
    c12=cos(x(1)+x(3));
    n1=m(2)*a(2)*l;
    V1=-n1*s2*x(4)^2-2*n1*s2*x(2)*x(4);
    V2=n1*s2*x(2)^2;
    G1=m(1)*a(1)*g*c1+m(2)*g*(l*c1+a(2)*c12);
    G2=m(2)*g*a(2)*c12;

    F(1)=x(2);
    F(2)=-V1-G1;
    F(3)=x(4);
    F(4)=-V2-G2;
    F=F';     
end

function M=fMass(t,x)
    m=[1 1];
    l=0.5;
    Izz=[0.11 0.11];
    a=[0.25 0.25];
    c2=cos(x(3));
    n1=m(2)*a(2)*l;
    M11=m(1)*a(1)^2+Izz(1)+m(2)*(a(2)^2+l^2)+2*n1*c2+Izz(2);
    M12=m(2)*a(2)^2+n1*c2+Izz(2);
    M22=m(2)*a(2)^2+Izz(2);
    M=[1 0 0 0;0 M11 0 M12;0 0 1 0;0 M12 0 M22];
end

Zauważ, jak ustawiłem warunek początkowy x1 (kąt pierwszego korpusu w stosunku do poziomu), aby system zaczął działać w pozycji całkowicie pionowej. W ten sposób, ponieważ działa tylko grawitacja, oczywistym rezultatem jest to, że system w ogóle nie powinien się poruszać z tej pozycji.

UWAGA: we wszystkich poniższych rysunkach przedstawiłem rozwiązania x1 i x3 w odniesieniu do czasu.

ODE45

Kiedy uruchamiam symulację przez 6 sekund ode45, otrzymuję oczekiwane rozwiązanie bez żadnych problemów, system pozostaje tam, gdzie jest i nie porusza się:

wprowadź opis zdjęcia tutaj

Jednak gdy uruchamiam symulację przez 10 sekund, system zaczyna się poruszać nieracjonalnie:

wprowadź opis zdjęcia tutaj

ODE23

Następnie przeprowadziłem symulację, ode23aby sprawdzić, czy problem nadal występuje. Skończyłem z tym samym zachowaniem, tylko tym razem dywergencja zaczyna się 1 sekundę później:

wprowadź opis zdjęcia tutaj

ODE15s

Następnie uruchomiłem symulację, ode15saby sprawdzić, czy problem nadal występuje i nie, system wydaje się być stabilny nawet przez 100 sekund:

wprowadź opis zdjęcia tutaj

Z drugiej strony, ode15sjest tylko pierwsze zamówienie i zauważ, że jest tylko kilka kroków integrujących. Uruchomiłem więc kolejną symulację w ode15sciągu 10 sekund, ale MaxSteprozmiar 0.01 aby zwiększyć precyzję, i niestety, prowadzi to do tego samego wyniku, co w przypadku obu ode45i ode23.

wprowadź opis zdjęcia tutaj

Zwykle oczywistym rezultatem tych symulacji byłby fakt, że system pozostaje w początkowej pozycji, ponieważ nic go nie zakłóca. Dlaczego występuje ta rozbieżność? Czy ma to coś wspólnego z faktem, że tego rodzaju systemy mają charakter chaotyczny? Czy jest to normalne zachowanie odefunkcji w Matlabie?

jrojasqu
źródło
Oprócz równań, myślę, że schemat bardzo pomógłby zrozumieć b pytanie.
nicoguaro
Jeśli uważasz, że jest to właściwe, możesz zaakceptować jedną z odpowiedzi (jest zielony przycisk).
Ertxiem
Nie mówisz, ale wydaje się, że planujesz x1i x3. (Wstaw suchy komentarz na temat wykresów bez legend i opisów.) Spróbuj wykreślić logarytmy (wartości bezwzględne) x2i x4.
Eric Towers

Odpowiedzi:

15

Myślę, że dwie główne kwestie zostały już przedstawione przez Briana i Ertxiema: twoja początkowa wartość jest niestabilną równowagą, a fakt, że twoje obliczenia liczbowe nigdy nie są naprawdę dokładne, zapewnia małe zaburzenie, które wywoła niestabilność.

Aby podać nieco więcej szczegółów, jak to się dzieje, rozważ swój problem w postaci ogólnego problemu z wartością początkową

y˙(t)=M1f(t,y(t))
y(t)=(x1(t),x2(t),x3(t),x4(t))

f(t,y(t))=[x2V1G1x4V2G2]

f(0,y0)=0y˙(0)=0f~

W kodzie możesz to sprawdzić, wykonując obliczenia

norm(rhs(0,[pi/2 0 0 0]))

co daje 6,191e-16 - czyli prawie, ale nie dokładnie zero. Jak to wpływa na dynamikę twojego systemu?

fy0y~0

Ponadto, w bardzo krótkim czasie, rozwiązanie twojego systemu wygląda jak rozwiązanie systemu zlinearyzowanego

y˙(t)=f(0,y0)+f(0,y0)(y(t)y0)=f(0,y0)(y(t)y0)

ffrhsy0d(t):=y(t)y0d

d˙(t)=f(0,y0)d(t).

Nie mogłem sobie pozwolić na ręczne obliczenie jakobianu, więc użyłem automatycznego różnicowania, aby uzyskać dobre przybliżenie:

J:=f(0,y0)=[01009.8102.4525000012.452502.45250]

tak aby twoje równanie stało się

d˙(t)=Jd(t),d(0)=y~0y0

Teraz potrzebujemy ostatniego kroku: możemy obliczyć rozkład wartości własnej jakobianu w taki sposób, aby

J=QDQ1

DJQde(t):=Q1d(t)

e˙(t)=De(t),e(0)=Q1d0.

D

e˙i(t)=λiei(t),ei(0)=ith component of Q1d0

i=1,2,3,4λ1=3.2485

e1(t)=e1(0)e3.2485t.

d(0)=0e(0)=Q1d(0)=0e1(0)=0e1(0)

Daniel
źródło
16

π/2π/2

Brian Borchers
źródło
4
Jeśli uważnie monitorujesz zmienne stanu (patrząc na wartości wydrukowane w notacji naukowej), powinieneś być w stanie zobaczyć początkowe bardzo powolne odchodzenie od równowagi.
Brian Borchers,
Ma to sens i rzeczywiście, kiedy uruchamiam system w pozycji pionowej w dół (będąc punktem stabilnej równowagi), system w ogóle się nie porusza, przynajmniej przez symulację 1000 sekund, którą uważam za bardzo długi okres czasu .
jrojasqu
6
x1sin(0)cos(0)sin(pi/2)cos(pi/2)rhs(t,[0,0,0 0] == [0,0,0,0]
π/2
1
θ=0 1016
4

Spójrz na składowe sił obliczonych w twoich funkcjach.

π

1016

a=1.0a=a+1016

alephzero
źródło
4

Początkowe założenie było takie, że pozycja początkowa znajdowała się w stabilnej równowadze (tj. Minimum energii potencjalnej) przy zerowej energii kinetycznej, a układ zaczął się oddalać od równowagi.
Ponieważ fizycznie tak się nie stanie (jeśli weźmiemy pod uwagę mechanikę klasyczną), przyszło mi do głowy dwie rzeczy:

  1. π/2)-π/2)

  2. Po drugie, być może coś jest nie tak z równaniami ruchu (może literówka gdzieś?). Czy możesz napisać równania wprost? Być może mógłbyś wykreślić przyspieszenie kątowe jako funkcję początkowej pozycji każdego wahadła, zakładając zerową prędkość kątową, aby sprawdzić, czy jest coś dziwnego.

Ertxiem - przywróć Monikę
źródło
1
Rzeczywiście, uruchomiłem system w pozycji pionowej w górę. Dlatego jest to punkt niestabilnej równowagi. Komentarz Briana Borchera uzupełnia twoją odpowiedź, wyjaśniając problemπzbliżenie, które prowadzi do przesunięcia się systemu w końcu z tej pozycji.
jrojasqu
2
Nawiasem mówiąc, dla zabawy, jeśli chcesz utrzymać układ w niestabilnej pozycji pionowej, możesz zmienić początek współrzędnych, aby kąt równy zero był skierowany w górę.
Ertxiem
@Ertxiem inną opcją jest wprowadzenie małego tarcia w pinach, który zjadłby błędy numeryczne.
svavil
@Ertxiem Dla zabawy próbowałem zmienić układ współrzędnych, aby kąt zerowy powodował, że układ jest skierowany w górę. To naprawdę najlepsza parametryzacja tutaj. Oczywiście system pozostaje w pozycji do góry przez czas nieokreślony. Jednak nadal pojawiają się oscylacje (minimalne przez 1000 sekund symulacji, ale tam) w stabilnej pozycji równowagi (prosto w dół), ponieważ wtedy istniejegrzech(π)do obliczenia w sile pochodzącej z energii potencjalnej. Dlatego nalegam na to, że jeśli symuluję to wystarczająco długo, system zacznie odchylać się także od tej pozycji.
jrojasqu
Ponieważ fizycznie nie może się to zdarzyć - biorąc pod uwagę, że jesteśmy w niestabilnej równowadze, nieco to podważam. Układy fizyczne (bez zbyt dużego tarcia) nie pozostają w niestabilnej równowadze. Mówiąc bardziej ogólnie, jeśli symulujesz prawdziwe systemy, chciałbyś uniknąć przypadkowego utknięcia w niestabilnej równowadze (jakkolwiek się tam dostało) - to cecha, a nie błąd. (Istnieją rzadkie wyjątki od tego, takie jak stan niezakażony w immunologii, który jest niestabilną równowagą, którą można utrzymać.)
Wrzlprmft
0

Powinieneś poszukać więcej o podwójnych wahadłach: są to tak zwane „systemy chaotyczne”. Mimo że zachowują się według prostych zasad, zaczynając od nieco innych warunków początkowych, rozwiązania różnią się dość szybko. Wykonywanie symulacji numerycznych dla tego rodzaju systemów nie jest łatwe. Obejrzyj poniższy film, aby uzyskać lepszy wgląd w problem.

Dla prostego lub podwójnego wahadła możesz napisać wzór na całkowitą energię układu. Zakładając, że siły tarcia są zaniedbywane, całkowita energia jest zachowywana przez system analityczny. Liczbowo jest to zupełnie inna kwestia.

Przed wypróbowaniem podwójnego wahadła wypróbuj prosty wahadło. Zauważysz, że dla metod małego rzędu Runge-Kutty energia systemu będzie rosła w symulacji numerycznych, zamiast pozostać stała (tak dzieje się w twoich symulacjach: ruch jest niczym z niczego). Aby temu zapobiec, można zastosować metody RK wyższego rzędu (ode45 jest rzędu 4; RK rzędu 8 działałby lepiej). Istnieją inne metody zwane „metodami symplektycznymi”, które są zaprojektowane tak, aby symulacje numeryczne oszczędzały energię. Zasadniczo należy przerwać symulację, gdy tylko energia znacznie wzrośnie w porównaniu do inicjalizacji.

I spróbuj zrozumieć proste wahadło, zanim przejdziesz do podwójnego.

Beni Bogosel
źródło
2
Nie jest to jednak kwestia chaosu systemu. Możesz mieć niestabilną równowagę również w układach nie chaotycznych, np. Pojedynczy wahadło jest „na głowie” i będzie wykazywał takie samo zachowanie opisane w pytaniu.
Daniel
1
Nie jest również prawdą, że energia rośnie dla RKM małego rzędu: niejawny Euler jest pierwszym rzędem i wykazuje dokładnie odwrotne zachowanie.
Daniel
@BeniBogosel Wspominasz metody symplektyczne, które zwróciły moją uwagę, ponieważ oczywiście w moim przykładzie energia nie jest zachowana. Czy mógłbyś jednak wskazać konkretną metodę symplektyczną, którą można tu zastosować?
jrojasqu
@jrojasqu, dlaczego według ciebie energia nie jest oszczędzana w twoim systemie?
Ertxiem
@ Ertxiem Kiedy obliczam całkowitą energię mechaniczną układu (energie kinetyczne + potencjalne) z wyjściami dostarczonymi przez ode45, otrzymuję wartość, która zaczyna się od zera, a następnie rośnie z czasem. Wartość jest bardzo mała w pierwszych sekundach, ale wciąż konsekwentnie rośnie od zera. Uważam, że dzieje się tak z powodu problemów, które zostały poruszone w powyższych odpowiedziach (przybliżenieπitp.).
jrojasqu