Modellbildung und Simulation, SS 2010
Blatt 5: Inverses Pendel
> | restart; |
> | with(DEtools): |
> | with(LinearAlgebra): |
> | with(plots): |
> | with(plottools): |
> |
Beispielgr??en
Zahlenwerte im Beispiel (Wagenmasse, Stabmasse, halbe Stabl?nge, Erdbeschleunigung):
> | subs_zahlen := {M=0.981,m=0.08,L=0.312,g=9.81}; |
> | K_ungeregelt := {Ks=0,Ksp=0,Kp=0,Kpp=0}; |
> | K_geregelt := {Ks=5.088,Ksp=5.258,Kp=35.39,Kpp=7.174}; |
> |
> |
Sinus und Cosinus
Richtig
Exakte trigonometrische Ausdr?cke:
> | SIN := phi->sin(phi); |
> | COS := phi->cos(phi); |
Linearisiert
Linearisierte trigonometrische Ausdr?cke (kleine
):
> | SIN := phi->phi; |
> | COS := phi->1; |
> |
Listen der vorkommenden Gr??en
Werden sp?ter bei der Termvereinfachung hilfreich sein:
Zustandsgr??en
> | vars := [s(t), diff(s(t),t), phi(t), diff(phi(t),t)]; |
> | subs_vars := seq(vars[i] = s||i,i=1..4); |
Die Sequenz subs_vars erlaubt, die Zustandsgr??en durch s1,..,s4 zu ersetzen:
> | diff(s(t),t)=XYZ; |
> | subs([subs_vars],%); |
Zweite Ableitungen
Die vorkommenden zweiten Ableitungen:
> | abls := diff(s(t),t$2); |
> | ablphi := diff(phi(t),t$2); |
> | abl := [abls,ablphi]; |
> |
> |
Bild vom Stab malen
> | bild := proc(s,phi)
local LL,stab,schwerpunkt; global subs_zahlen; LL := subs(subs_zahlen,L); stab := line([s,0],[s+2*LL*SIN(phi),2*LL*COS(phi)], color=RED,thickness=3): schwerpunkt := disk([s+LL*SIN(phi),LL*COS(phi)],LL/10, color=BLUE); display([stab,schwerpunkt],view=[-2*LL..2*LL,0..2*LL], scaling=constrained) end: |
> | bild(0,1); |
> | display([seq(bild(0,T/10),T=-10..10)],insequence=true,scaling=constrained); |
> |
Bahn des Stabschwerpunkts
Bahn des Stabschwerpunkts (horizontal und vertikal):
> | x := t -> s(t) + L*SIN(phi(t)); |
> | y := t -> y0 + L*COS(phi(t)); |
> |
Beschleunigung des Stabs
H,V: Kraft am Stab (horizontal und vertikal).
Da H und V keine Zustandsgr??en sind, also eliminiert werden sollen, schreiben wir sie isoliert auf die linke Seite der Gleichungen. Dann brauchen wir gar keine Gleichung mehr zu speichern, sondern nur noch die rechte Seite, so dass in Zukunft H und V zu Ausdr?cken in Zustandsgr??en expandiert werden.
> | H := m*diff(x(t),t$2); |
> | V := m*diff(y(t),t$2)+m*g; |
> |
Drehmoment
Tr?gheitsmoment des Stabs:
> | theta := m*L^2/3; |
Stab-Drehmoment.
Um Ordnung zu schaffen, verwenden wir die Liste der zweiten Ableitungen und lassen uns den Term danach sortieren.
Das Ergebnis ist eine weitere Gleichung f?r unser System.
> | theta*diff(phi(t),t$2) - L*(-H*COS(phi(t)) + V*SIN(phi(t)))=0; |
> | gl1 := collect(%,abl); |
> |
Kraft auf Wagen
Kraft auf Wagen (nur horizontale Komponente):
> | M*diff(s(t),t$2) - (u(t)-H) =0; |
> | gl2 := collect(%,abl); |
> |
Stellkraft
Die Stellkraft u(t) ist ein linearer Ausdruck in den Zustandsgr??en (linearer PD-Regler):
> | subs_u := u(t) = Kp*phi(t)+Ks*s(t)+Kpp*diff(phi(t),t)+Ksp*diff(s(t),t); |
> |
Umformen in System linearer DGLn erster Ordnung
Da wir H und V bereits eliminiert haben haben wir nur noch gl1 (Drehmoment) und gl2 (Kraft auf Wagen):
> | gl1; |
> | gl2; |
Noch kommen beide zweiten Ableitungen in beiden Gleichungen vor. Das m?ssen wir ?ndern (in der Form
d?rfen die zweiten Ableitungen nur auf der linken Seite Platz nehmen und zwar immer nur eine!). Das ist aber ganz gew?hnliches Eliminieren von
und Einsetzen:
> | gl1 - 4/3*L*gl2; |
> | gl3 := simplify(isolate(%,abls)); |
> | gl4 := simplify(isolate(subs(gl3,gl1),ablphi)); |
Setzen wir die Stellkraft ein:
> | gl3a := subs(subs_u,gl3); |
> | gl4a := subs(subs_u,gl4); |
Und wir bekommen das gesuchte DGL-System
(Fehlen uns nicht zwei Gleichungen? Nein, das sind die beiden trivialen Gleichungen, mit denen gesagt wird, wie
und
sowie
und
zusammenh?ngen):
> | gls := [diff(s(t),t)=diff(s(t),t), gl3a, diff(phi(t),t)=diff(phi(t),t), gl4a]; |
> |
Koeffizientenmatrix aufstellen
Bisher liegt unser DGL-System als Liste von Gleichungen vor. Zur Bestimmung der Eigenwerte brauchen wir aber die Koeffizientenmatix.
Von Hand geht das schnell, aber es ist etwas m?hsam, das Maple beizubringen. Im Prinzip geht das mit coeff und einigen Vorbereitungsschritten:
Rechte Seiten der DGLn:
> | map(rhs,gls); |
Umbenennen in s1,...,s4 (gibt sonst ?rger, weil z.B.
auch in
vorkommt:
> | subs([subs_vars],%); |
Sortieren:
> | collect(%,[s1,s2,s3,s4]); |
Und Koeffizienten extrahieren:
> | coeff(%[2],s2); |
Das ganze als Prozedur (f?ngt auch noch die F?lle ab, wo die Gr??e gar nicht im Ausdruck vorkommt, der Koeffizient also 0 ist):
Es wird der Koeffizient von Zustandsgr??e j in Gleichung i zur?ckgegeben:
> | getcoeff := proc(i,j)
global gls,vars; try coeff(subs([subs_vars],collect(expand(rhs(gls[i])),vars)),subs([subs_vars],vars[j])) catch: 0 end try end: |
Und das ist der Lohn der M?he:
> | A := Matrix(4,(i,j) -> getcoeff(i,j)); |
> |
Stabilit?tsanalyse
Es geht um die Eigenwerte von A, also die Nullstellen seines charakteristischen Polynoms:
> | collect(CharacteristicPolynomial(subs(subs_zahle=1,A),lambda),lambda); |
> |
Ungeregelter Fall
Da sieht die Matrix ?bersichtlich aus:
> | A0 := subs(subs_zahlen,K_ungeregelt,A); |
Daf?r sind die Eigenwerte gar nicht gut (vor allem der mit positivem Realteil i.e. der umkippende Stab)!
> | Eigenvectors(A0); |
> |
Geregelter Fall
Nun ist das A zwar komplizierter...
> | A1 := subs(subs_zahlen,K_geregelt,A); |
... daf?r sind die Realteile aller Eigenwerte deutlich negativ: Stabilit?t!
> | Eigenvectors(A1); |
> |
Und der geregelte Fall mit bewegtem Stab
> | Tend := 5; |
> | DGLs := subs(K_geregelt,subs_zahlen,[gl3a,gl4a]); |
> | Dsol := dsolve({op(DGLs),s(0)=0,D(s)(0)=1,phi(0)=0,D(phi)(0)=0.1},numeric,{s(t),phi(t)},range=0..Tend); |
> | Dsol(1); |
> | display(seq(bild(
op(evalf(subs(subs(Dsol(Tend*T/100),[phi(t),s(t)])))) ),T=0..100),insequence=true,scaling=constrained); |
> |