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

SIN := proc (phi) options operator, arrow; sin(phi) end proc

COS := proc (phi) options operator, arrow; cos(phi) end proc

Linearisiert

Linearisierte trigonometrische Ausdr?cke (kleine phi ):

> SIN := phi->phi;

> COS := phi->1;

SIN := proc (phi) options operator, arrow; phi end proc

COS := proc (phi) options operator, arrow; 1 end proc

>

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

vars := [s(t), diff(s(t), t), phi(t), diff(phi(t), t)]

subs_vars := s(t) = s1, diff(s(t), t) = s2, phi(t) = s3, diff(phi(t), t) = s4

Die Sequenz subs_vars erlaubt, die Zustandsgr??en durch s1,..,s4 zu ersetzen:

> diff(s(t),t)=XYZ;

> subs([subs_vars],%);

diff(s(t), t) = XYZ

s2 = XYZ

Zweite Ableitungen

Die vorkommenden zweiten Ableitungen:

> abls := diff(s(t),t$2);

> ablphi := diff(phi(t),t$2);

> abl := [abls,ablphi];

>

abls := diff(s(t), `$`(t, 2))

ablphi := diff(phi(t), `$`(t, 2))

abl := [diff(s(t), `$`(t, 2)), diff(phi(t), `$`(t, 2))]

>

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

[Plot]

> display([seq(bild(0,T/10),T=-10..10)],insequence=true,scaling=constrained);

[Plot]

>

Bahn des Stabschwerpunkts

Bahn des Stabschwerpunkts (horizontal und vertikal):

> x := t -> s(t) + L*SIN(phi(t));

> y := t -> y0 + L*COS(phi(t));

x := proc (t) options operator, arrow; s(t)+L*SIN(phi(t)) end proc

y := proc (t) options operator, arrow; y0+L*COS(phi(t)) end proc

>

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

H := m*((diff(s(t), `$`(t, 2)))+L*(diff(phi(t), `$`(t, 2))))

> V := m*diff(y(t),t$2)+m*g;

V := m*g

>

Drehmoment

Tr?gheitsmoment des Stabs:

> theta := m*L^2/3;

theta := 1/3*m*L^2

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

1/3*m*L^2*(diff(phi(t), `$`(t, 2)))-L*(-m*((diff(s(t), `$`(t, 2)))+L*(diff(phi(t), `$`(t, 2))))+m*g*phi(t)) = 0

gl1 := L*m*(diff(s(t), `$`(t, 2)))+4/3*m*L^2*(diff(phi(t), `$`(t, 2)))-L*m*g*phi(t) = 0

>

Kraft auf Wagen

Kraft auf Wagen (nur horizontale Komponente):

> M*diff(s(t),t$2) -  (u(t)-H) =0;

> gl2 := collect(%,abl);

M*(diff(s(t), `$`(t, 2)))-u(t)+m*((diff(s(t), `$`(t, 2)))+L*(diff(phi(t), `$`(t, 2)))) = 0

gl2 := (M+m)*(diff(s(t), `$`(t, 2)))-u(t)+m*L*(diff(phi(t), `$`(t, 2))) = 0

>

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

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;

L*m*(diff(s(t), `$`(t, 2)))+4/3*m*L^2*(diff(phi(t), `$`(t, 2)))-L*m*g*phi(t) = 0

(M+m)*(diff(s(t), `$`(t, 2)))-u(t)+m*L*(diff(phi(t), `$`(t, 2))) = 0

Noch kommen beide zweiten Ableitungen in beiden Gleichungen vor. Das m?ssen wir ?ndern (in der Form diff(x, t) = A*x 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 diff(phi(t), `$`(t, 2)) und Einsetzen:

> gl1 - 4/3*L*gl2;

> gl3 := simplify(isolate(%,abls));

L*m*(diff(s(t), `$`(t, 2)))+4/3*m*L^2*(diff(phi(t), `$`(t, 2)))-L*m*g*phi(t)-4/3*L*((M+m)*(diff(s(t), `$`(t, 2)))-u(t)+m*L*(diff(phi(t), `$`(t, 2)))) = 0

gl3 := diff(s(t), `$`(t, 2)) = -(3*m*g*phi(t)-4*u(t))/(m+4*M)

> gl4 := simplify(isolate(subs(gl3,gl1),ablphi));

gl4 := diff(phi(t), `$`(t, 2)) = 3*(m*g*phi(t)-u(t)+g*phi(t)*M)/(L*(m+4*M))

Setzen wir die Stellkraft ein:

> gl3a := subs(subs_u,gl3);

> gl4a := subs(subs_u,gl4);

gl3a := diff(s(t), `$`(t, 2)) = -(3*m*g*phi(t)-4*Kp*phi(t)-4*Ks*s(t)-4*Kpp*(diff(phi(t), t))-4*Ksp*(diff(s(t), t)))/(m+4*M)

gl4a := diff(phi(t), `$`(t, 2)) = 3*(m*g*phi(t)-Kp*phi(t)-Ks*s(t)-Kpp*(diff(phi(t), t))-Ksp*(diff(s(t), t))+g*phi(t)*M)/(L*(m+4*M))

Und wir bekommen das gesuchte DGL-System  diff(x, t) = A*x (Fehlen uns nicht zwei Gleichungen? Nein, das sind die beiden trivialen Gleichungen, mit denen gesagt wird, wie s(t) und diff(s(t), t) sowie phi(t) und diff(phi(t), t) zusammenh?ngen):

> gls := [diff(s(t),t)=diff(s(t),t), gl3a, diff(phi(t),t)=diff(phi(t),t), gl4a];

gls := [diff(s(t), t) = diff(s(t), t), diff(s(t), `$`(t, 2)) = -(3*m*g*phi(t)-4*Kp*phi(t)-4*Ks*s(t)-4*Kpp*(diff(phi(t), t))-4*Ksp*(diff(s(t), t)))/(m+4*M), diff(phi(t), t) = diff(phi(t), t), diff(phi(...gls := [diff(s(t), t) = diff(s(t), t), diff(s(t), `$`(t, 2)) = -(3*m*g*phi(t)-4*Kp*phi(t)-4*Ks*s(t)-4*Kpp*(diff(phi(t), t))-4*Ksp*(diff(s(t), t)))/(m+4*M), diff(phi(t), t) = diff(phi(t), t), diff(phi(...

>

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

[diff(s(t), t), -(3*m*g*phi(t)-4*Kp*phi(t)-4*Ks*s(t)-4*Kpp*(diff(phi(t), t))-4*Ksp*(diff(s(t), t)))/(m+4*M), diff(phi(t), t), 3*(m*g*phi(t)-Kp*phi(t)-Ks*s(t)-Kpp*(diff(phi(t), t))-Ksp*(diff(s(t), t))+...

Umbenennen in s1,...,s4 (gibt sonst ?rger, weil z.B. s(t) auch in diff(s(t), t) vorkommt:

> subs([subs_vars],%);

[s2, -(3*m*g*s3-4*Kp*s3-4*Ks*s1-4*Kpp*s4-4*Ksp*s2)/(m+4*M), s4, 3*(m*g*s3-Kp*s3-Ks*s1-Kpp*s4-Ksp*s2+g*s3*M)/(L*(m+4*M))]

Sortieren:

> collect(%,[s1,s2,s3,s4]);

[s2, 4*Ks*s1/(m+4*M)+4*Ksp*s2/(m+4*M)-(3*m*g-4*Kp)*s3/(m+4*M)+4*Kpp*s4/(m+4*M), s4, -3*Ks*s1/(L*(m+4*M))-3*Ksp*s2/(L*(m+4*M))+3*(m*g-Kp+g*M)*s3/(L*(m+4*M))-3*Kpp*s4/(L*(m+4*M))]

Und Koeffizienten extrahieren:

> coeff(%[2],s2);

4*Ksp/(m+4*M)

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

A := Matrix([[0, 1, 0, 0], [4*Ks/(m+4*M), 4*Ksp/(m+4*M), -3*m*g/(m+4*M)+4*Kp/(m+4*M), 4*Kpp/(m+4*M)], [0, 0, 0, 1], [-3*Ks/(L*(m+4*M)), -3*Ksp/(L*(m+4*M)), 3*m*g/(L*(m+4*M))-3*Kp/(L*(m+4*M))+3*g*M/(L*...

>

Stabilit?tsanalyse

Es geht um die Eigenwerte von A, also die Nullstellen seines charakteristischen Polynoms:

> collect(CharacteristicPolynomial(subs(subs_zahle=1,A),lambda),lambda);

lambda^4-(-3*Kpp+4*Ksp*L)*lambda^3/(L*(m+4*M))-(3*g*M-3*Kp+3*m*g+4*Ks*L)*lambda^2/(L*(m+4*M))+3*Ksp*g*lambda/(L*(m+4*M))+3*Ks*g/(L*(m+4*M))

>

Ungeregelter Fall

Da sieht die Matrix ?bersichtlich aus:

> A0 := subs(subs_zahlen,K_ungeregelt,A);

A0 := Matrix([[0, 1, 0, 0], [0., 0., -.5880119880, 0.], [0, 0, 0, 1], [0., 0., 24.99522113, 0.]])

Daf?r sind die Eigenwerte gar nicht gut (vor allem der mit positivem Realteil i.e. der umkippende Stab)!

> Eigenvectors(A0);

Vector[column]([[0.+0.*I], [0.+0.*I], [4.99952209016021864+0.*I], [-4.99952209016021864+0.*I]]), Matrix([[1.+0.*I, -1.+0.*I, -0.461277527921204970e-2+0.*I, 0.461277527921204970e-2+0.*I], [0.+0.*I, 0.2...

>

Geregelter Fall

Nun ist das A zwar komplizierter...

> A1 := subs(subs_zahlen,K_geregelt,A);

A1 := Matrix([[0, 1, 0, 0], [5.082917084, 5.252747254, 34.76663337, 7.166833168], [0, 0, 0, 1], [-12.21855068, -12.62679628, -59.99190712, -17.22796434]])

... daf?r sind die Realteile aller Eigenwerte deutlich negativ: Stabilit?t!

> Eigenvectors(A1);

Vector[column]([[-4.94794570165402003+0.*I], [-2.00811769082704528+2.00321412034486590*I], [-2.00811769082704528-2.00321412034486590*I], [-3.01103600269189275+0.*I]]), Matrix([[0.303059067825318462e-2...Vector[column]([[-4.94794570165402003+0.*I], [-2.00811769082704528+2.00321412034486590*I], [-2.00811769082704528-2.00321412034486590*I], [-3.01103600269189275+0.*I]]), Matrix([[0.303059067825318462e-2...Vector[column]([[-4.94794570165402003+0.*I], [-2.00811769082704528+2.00321412034486590*I], [-2.00811769082704528-2.00321412034486590*I], [-3.01103600269189275+0.*I]]), Matrix([[0.303059067825318462e-2...Vector[column]([[-4.94794570165402003+0.*I], [-2.00811769082704528+2.00321412034486590*I], [-2.00811769082704528-2.00321412034486590*I], [-3.01103600269189275+0.*I]]), Matrix([[0.303059067825318462e-2...Vector[column]([[-4.94794570165402003+0.*I], [-2.00811769082704528+2.00321412034486590*I], [-2.00811769082704528-2.00321412034486590*I], [-3.01103600269189275+0.*I]]), Matrix([[0.303059067825318462e-2...Vector[column]([[-4.94794570165402003+0.*I], [-2.00811769082704528+2.00321412034486590*I], [-2.00811769082704528-2.00321412034486590*I], [-3.01103600269189275+0.*I]]), Matrix([[0.303059067825318462e-2...

>

Und der geregelte Fall mit bewegtem Stab

> Tend := 5;

Tend := 5

> DGLs := subs(K_geregelt,subs_zahlen,[gl3a,gl4a]);

DGLs := [diff(s(t), `$`(t, 2)) = 34.76663337*phi(t)+5.082917084*s(t)+7.166833168*(diff(phi(t), t))+5.252747254*(diff(s(t), t)), diff(phi(t), `$`(t, 2)) = -59.99190714*phi(t)-12.21855068*s(t)-17.227964...DGLs := [diff(s(t), `$`(t, 2)) = 34.76663337*phi(t)+5.082917084*s(t)+7.166833168*(diff(phi(t), t))+5.252747254*(diff(s(t), t)), diff(phi(t), `$`(t, 2)) = -59.99190714*phi(t)-12.21855068*s(t)-17.227964...

> 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 := proc (x_rkf45) local res, data, vars, solnproc, outpoint, ndsol, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if...

> Dsol(1);

[t = 1., phi(t) = -0.595123997048194325e-1, diff(phi(t), t) = .379833963077229164, s(t) = .386088140850482875, diff(s(t), t) = -.566231610006889929]

> display(seq(bild(
 op(evalf(subs(subs(Dsol(Tend*T/100),[phi(t),s(t)]))))

),T=0..100),insequence=true,scaling=constrained);

[Plot]

>