inv_pendel.html

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};

$subs_zahlen := {M = .981, m = 0.8e-1, L = .312, g = 9.81}$

> K_ungeregelt := {Ks=0,Ksp=0,Kp=0,Kpp=0};

$K_ungeregelt := {Ks = 0, Ksp = 0, Kpp = 0, Kp = 0}$

> K_geregelt := {Ks=5.088,Ksp=5.258,Kp=35.39,Kpp=7.174};

$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: