Modellbildung und Simulation, SS 2009

Blatt 9: Fuzzy-Regelung

> restart:
with(plots):

with(plottools):

with(LinearAlgebra):

UseHardwareFloats := true;

UseHardwareFloats := true

Modell des inversen Pendels

Unser Modell vom letzten Blatt mit der Linearisierung f?r kleine Winkel phi :

Linearisierte trigonometrische Ausdr?cke (kleine phi ):

> SIN := t->t;

> COS := t->1;

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

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

Werte f?r Wagen und Stab

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

Zustandsvariablen: Wagenposition und -geschwindigkeit; Auslenkwinkel und Winkelgeschwindigkeit

> zust := <sp, spp, pp, ppp>;

zust := Vector[column]([[sp], [spp], [pp], [ppp]])

> A := <<0,0,0,0>|<1,0,0,0>|<0,-3*m*g/(4*M+m),0,3*g*(M+m)/(L*(4*M+m))>|<0,0,1,0>>;

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

> A_num := subs(subs_zahlen, A);

A_num := Matrix([[0, 1, 0, 0], [0, 0, -.5880119880, 0], [0, 0, 0, 1], [0, 0, 24.99522112, 0]])

> B := <0,-4/(4*M+m),0,3/(L*(4*M+m))>;

B := Vector[column]([[0], [-4/(4*M+m)], [0], [3/(L*(4*M+m))]])

> B_num := subs(subs_zahlen, B);

B_num := Vector[column]([[0], [-.9990009992], [0], [2.401444709]])

> K := <5.088 | 5.258 | 35.39 | 7.174>;

K := Vector[row]([5.088, 5.258, 35.39, 7.174])

Im ungeregelten Fall ist das System instabil:

> Eigenvalues(A_num);

Vector[column]([[0.+0.*I], [0.+0.*I], [4.99952208916012264+0.*I], [-4.99952208916012352+0.*I]])

Im geregelten Fall ist das System stabil:

> Eigenvalues(A_num-B_num.K);

Vector[column]([[-4.94794570689810609+0.*I], [-2.00811768961696080+2.00321411908876268*I], [-2.00811768961696080-2.00321411908876268*I], [-3.01103600244036373+0.*I]])

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

 return display([stab,schwerpunkt],

                view=[-2*LL..2*LL,0..2*LL],scaling=constrained)

end:

Ein Euler-Schritt der L?nge delta_t f?r Zustandsvektor x. Die Funktion ufkt ist die Abbildung x -> u

> euler_schritt := (X, ufkt, delta_t) -> X + delta_t*(A_num.X - B_num.ufkt(X));

euler_schritt := proc (X, ufkt, delta_t) options operator, arrow; X+delta_t*((A_num . X . )-(B_num . ufkt(X) . )) end proc

N Schritte mit Startvektor X0. Ergebnis ist der Endvektor; ein Film wird erzeugt.

> simulation := proc(X0, ufkt, delta_t, N)
 local bilder, i, X;

 bilder := table();

 X := X0;

 bilder[0] := bild(X[1], X[3]);

 for i to N do

   X := euler_schritt(X, ufkt, delta_t);

   bilder[i] := bild(X[1], X[3]);

 end do;

 print(display([seq(bilder[i],i=0..N)],insequence=true,scaling=constrained));

 X;

end proc:

Ungeregelt:

> simulation(<0,0,evalf(-10/360*2*Pi),evalf(40/360*2*Pi)>, x->0, 0.05, 70);

[Plot]

Vector[column]([[2490.97594564788778], [12453.4388304821041], [-105884.168486017210], [-529370.239238192444]])

Geregelt:

> simulation(<0,0,evalf(-10/360*2*Pi),evalf(40/360*2*Pi)>, x->K.x, 0.05, 50);

[Plot]

Vector[column]([[0.149434174006810830e-2], [-0.208440162317456322e-2], [-0.791139200912121198e-3], [0.375461474334099539e-2]])

>

Fuzzy-Mengen

Mit der piecewise-Funktion lassen sich auf einfache Art und Weise st?ckweise lineare Zugeh?rigkeitsfunktionen zu Fuzzy-Mengen darstellen und zeichnen:

> ZeroApprox := piecewise(x<-1, 0, x<0, 1+x, x<1, 1-x, 0);
Positive   := piecewise(x<0, 0, x<1, x, x<2, 2-x, 0);

ZeroApprox := PIECEWISE([0, x < -1], [1+x, x < 0], [1-x, x < 1], [0, otherwise])

Positive := PIECEWISE([0, x < 0], [x, x < 1], [2-x, x < 2], [0, otherwise])

> plot([ZeroApprox, Positive], x=-3..3, thickness=3);

[Plot]

Die Verwendung von piecewise-Funktionen erm?glicht eine sehr einfache Verkn?pfung von Fuzzy-Mengen ?ber ihre Zugeh?rigkeitsfunktionen:

> erg := max(min(0.5, ZeroApprox), Positive);

erg := max(min(.5, PIECEWISE([0, x < -1], [1+x, x < 0], [1-x, x < 1], [0, otherwise])), PIECEWISE([0, x < 0], [x, x < 1], [2-x, x < 2], [0, otherwise]))

> plot(erg, x=-3..3, thickness=3);

[Plot]

>

Erstellen der Fuzzy-Mengen f?r phi , d*phi/dt und u .

> fuzzySet := proc()
 local l;

 l := (x<args[2], args[1],

       seq(op([x<args[2*i+2],

              args[2*i+1]*(args[2*i+2]-x)/(args[2*i+2]-args[2*i])

              +args[2*i+3]*(args[2*i]-x)/(args[2*i]-args[2*i+2])]),i=1..floor(nargs/2.0)-1),

       args[nargs]);

 return piecewise(l);

end proc:

> fuzzySets := proc()
 return (fuzzySet(1, args[1], 1, args[2], 0),

         seq(fuzzySet(0, args[i], 0, args[i+1], 1, args[i+2], 0), i=1..nargs-2),

         fuzzySet(0, args[nargs-1], 0, args[nargs], 1));

end proc:

> plot([fuzzySets(1,2,3)],x=-1..4,color=[RED,BLUE,GREEN], thickness=3);

[Plot]

Wir betrachten nur den aufgerichteten, linearisierbaren Fall, d.h. wir w?hlen Winkel zwischen -10 und +10 Grad:

Winkel m?ssen wir allerdings im Bogenma? und nicht im Gradma? angeben:

> grad2rad := proc(g)
 g/360*2*Pi;

end proc:

F?r unseren Winkel phi gibt dies dann:

> (pNL, pN, pZ, pP, pPL) := fuzzySets(grad2rad(-10),grad2rad(-5),0,grad2rad(5),grad2rad(10));

pNL, pN, pZ, pP, pPL := PIECEWISE([1, x < -1/18*Pi], [36*(-1/36*Pi-x)/Pi, x < -1/36*Pi], [0, otherwise]), PIECEWISE([0, x < -1/18*Pi], [-36*(-1/18*Pi-x)/Pi, x < -1/36*Pi], [-36*x/Pi, x < 0], [0, other...

> plot([pNL,pN,pZ,pP,pPL], x=grad2rad(-15)..grad2rad(15), thickness=3);

[Plot]

>

Entsprechend beschr?nken wir uns auf Winkelgeschwindigkeiten d*phi/dt zwischen circa -40 und 40 Grad pro Sekunde
Das entspricht in etwa :

> (ppNL, ppN, ppZ, ppP, ppPL) := fuzzySets(grad2rad(-40),grad2rad(-20),0,grad2rad(20),grad2rad(40));

ppNL, ppN, ppZ, ppP, ppPL := PIECEWISE([1, x < -2/9*Pi], [9*(-1/9*Pi-x)/Pi, x < -1/9*Pi], [0, otherwise]), PIECEWISE([0, x < -2/9*Pi], [-9*(-2/9*Pi-x)/Pi, x < -1/9*Pi], [-9*x/Pi, x < 0], [0, otherwise...

> plot([ppNL,ppN,ppZ,ppP,ppPL], x=grad2rad(-45)..grad2rad(45), thickness=3);

[Plot]

>

F?r die Ausgabe kommen Werte zwischen -15 und +15 Newton in Frage:

> (uNL, uN, uZ, uP, uPL) := fuzzySets(-15,-7.5,0,7.5,15);
u_range := -15..15;

uNL, uN, uZ, uP, uPL := PIECEWISE([1, x < -15], [-.9999999998-.1333333333*x, x < -7.5], [0, otherwise]), PIECEWISE([0, x < -15], [2.000000000+.1333333333*x, x < -7.5], [-.1333333333*x, x < 0], [0, oth...uNL, uN, uZ, uP, uPL := PIECEWISE([1, x < -15], [-.9999999998-.1333333333*x, x < -7.5], [0, otherwise]), PIECEWISE([0, x < -15], [2.000000000+.1333333333*x, x < -7.5], [-.1333333333*x, x < 0], [0, oth...

u_range := -15 .. 15

> plot([uNL,uN,uZ,uP,uPL], x=u_range, thickness=3);

[Plot]

Aufstellen der Regelbasis

Das Aufstellen der Regeln geht in Maple sehr sch?n und ?ber die symbolischen Bezeichner, da die Auswertung unterbunden werden kann.
Wir erhalten eine Matrix (Zeilen: Winkelgeschwindigkeit) x (Spalten: Winkel).

Die beiden Listen Zeilen und Spalten helfen uns nachher, bei der Fuzzifizierung: Konkrete Werte f?r phi m?ssen auf alle Eintr?ge von Spalten angewendet werden, entsprechend Werte f?r d*phi/dt auf Zeilen.

> Zeilen := ['ppNL', 'ppN', 'ppZ', 'ppP', 'ppPL'];

> Spalten := ['pNL', 'pN', 'pZ', 'pP', 'pPL'];

> Regeln := Matrix([[0, 0, 'uNL', 0, 0],
                 [0, 'uNL', 'uN', 'uZ', 0],

                 ['uNL', 'uN', 'uZ', 'uP', 'uPL'],

                 [0, 'uZ', 'uP', 'uPL', 0],

                 [0, 0, 'uPL', 0, 0]]);

Zeilen := [ppNL, ppN, ppZ, ppP, ppPL]

Spalten := [pNL, pN, pZ, pP, pPL]

Regeln := Matrix([[0, 0, uNL, 0, 0], [0, uNL, uN, uZ, 0], [uNL, uN, uZ, uP, uPL], [0, uZ, uP, uPL, 0], [0, 0, uPL, 0, 0]])

Alternative Regeln, ohne 0-Eintr?ge

> #Regeln := Matrix([['uNL', 'uNL', 'uNL', 'uN', 'uZ'],
#                  ['uNL', 'uNL', 'uN', 'uZ', 'uP'],

#                  ['uNL', 'uN', 'uZ', 'uP', 'uPL'],

#                  ['uN', 'uZ', 'uP', 'uPL', 'uPL'],

#                  ['uZ', 'uP', 'uPL', 'uPL', 'uPL']]);

>

Defuzzifizierung

Die Stellgr??e muss defuzzifiziert werden, d.h. die x-Koordinate des Schwerpunkts der Zugeh?rigkeitsfunktion muss bestimmt werden:

> defuzzify := proc(f, r::range)
 return 1/evalf(Int(f, x=r))*evalf(Int(x*f,x=r));

end proc:

>

Fuzzifizierung und Inferenz

Zum Fuzzifizieren m?ssen nur f?r konkrete Werte alle entsprechenden Zugeh?rigkeitsfunktionen ausgewertet werden.
Die Inferenz erfordert einen Durchlauf durch die Regelmatrix.

Aufstellen der Ergebnismenge f?r die Steuergr??e:

> fuzzifyInference := proc(phi_pkt::numeric, Zeilen::list, phi::numeric, Spalten::list, Regeln::Matrix)
 local l1, l2, i1, i2, erg;

 # Fuzzifizieren

 l1 := evalf([seq(unapply(Zeilen[i], x)(phi_pkt), i=1..nops(Zeilen))]);

 l2 := evalf([seq(unapply(Spalten[i], x)(phi), i=1..nops(Spalten))]);

 # Inferenz

 erg := 0;

 for i1 from 1 to nops(Zeilen) do

   for i2 from 1 to nops(Spalten) do

     if l1[i1]<>0 and l2[i2]<>0 then

       erg := max(erg, min(min(l1[i1],l2[i2]),Regeln[i1,i2]));

     end if;

   end do:

 end do:

 return erg;

end proc:

>

Und ein konkretes Beispiel:

> phi := 0.1; phi_pkt := -0.2;

phi := .1

phi_pkt := -.2

> st := time():
erg := fuzzifyInference(phi_pkt, Zeilen, phi, Spalten, Regeln);

time()-st;

st := time():

defuzzify(erg, u_range);

time()-st;

erg := max(0, min(PIECEWISE([0, x < 0], [.1333333333*x, x < 7.5], [2.000000000-.1333333333*x, x < 15], [0, otherwise]), .4270422048), min(.1459155900, PIECEWISE([0, x < 7.5], [-.9999999998+.1333333333...erg := max(0, min(PIECEWISE([0, x < 0], [.1333333333*x, x < 7.5], [2.000000000-.1333333333*x, x < 15], [0, otherwise]), .4270422048), min(.1459155900, PIECEWISE([0, x < 7.5], [-.9999999998+.1333333333...

0.32e-1

3.406420661

2.056

Die Fuzzifizierung und die Inferenz gingen ja ziemlich schnell. Nur das Defuzzifizieren ben?tigt viel zu viel Zeit...

Wir zeichnen die Zugeh?rigkeitsfunktion der Ergebnismenge:

> plot(erg, x=u_range, thickness=3);

[Plot]

>

Fuzzifizierung und Inferenz II (optimiert):

Da das Defuzzifizieren wegen der min- und max-Prozeduren viel zu lange dauert, implementieren wir uns eigene MIN und MAX-Prozeduren, die das f?r st?ckweise lineare Funktionen viel schneller und einfacher machen:

> # MAX:
MAX := proc(pl1::list, pl2::list)

#option trace;

 # lokale Variablen

 local i1, i2, x_l, x_r, x1_r, x2_r, xs, f1, f2, f1_l, f2_l, f1_r, f2_r, pe;


 # Initialisierungen

 i1 := 1; i2 := 1;

 x1_r := `if`(nops(pl1)=1, infinity, pl1[2*i1]);

 x2_r := `if`(nops(pl2)=1, infinity, pl2[2*i2]);

 f1 := pl1[1]; f2 := pl2[1];

 pe := [];


 x_l := -infinity;

 f1_l := subs(x=x_l, f1);

 f2_l := subs(x=x_l, f2);


 # Schleife, bis beide Intervalle abschnittsweise durchlaufen

 x_r := min(x1_r, x2_r);

 while (x_r < infinity) do

   f1_r := subs(x=x_r, f1);

   f2_r := subs(x=x_r, f2);

   if (f1_l <= f2_l) and (f1_r <= f2_r) then

     # f1 <= f2

     pe := [op(pe), f2, x_r];

   elif (f1_l >= f2_l) and (f1_r >= f2_r) then

     # f1 >= f2

     pe := [op(pe), f1, x_r];

   else

     # Schnitt von f1 und f2

     xs := solve(f1=f2, x);

     if (f1_l < f2_l) then

       # linkes Ende: f1 < f2

       pe := [op(pe), f2, xs, f1, x_r];

     else

       # linkes Ende: f1 < f2

       pe := [op(pe), f1, xs, f2, x_r];

     end if;

   end if;

   # Punkte verschieben:

   if (x1_r < x2_r) then

     # f1 ein Intervall weiter:

     i1 := i1 + 1;

     x1_r := `if`(nops(pl1)<2*i1, infinity, pl1[2*i1]);

     f1 := pl1[2*i1-1];

     f1_l := f1_r;

   else

     # f1 ein Intervall weiter:

     i2 := i2 + 1;

     x2_r := `if`(nops(pl2)<2*i2, infinity, pl2[2*i2]);

     f2 := pl2[2*i2-1];

     f2_l := f2_r;

   end if;

   x_r := min(x1_r, x2_r);  

 end do:

 f1_r := subs(x=x_r, f1);

 f2_r := subs(x=x_r, f2);

 return [op(pe), `if`(f1_r < f2_r, f2, f1)];

end proc:

> # MIN:
MIN := proc(pl1::list, pl2::list)

#option trace;

 # lokale Variablen

 local i1, i2, x_l, x_r, x1_r, x2_r, xs, f1, f2, f1_l, f2_l, f1_r, f2_r, pe;


 # Initialisierungen

 i1 := 1; i2 := 1;

 x1_r := `if`(nops(pl1)=1, infinity, pl1[2*i1]);

 x2_r := `if`(nops(pl2)=1, infinity, pl2[2*i2]);

 f1 := pl1[1]; f2 := pl2[1];

 pe := [];


 x_l := -infinity;

 f1_l := subs(x=x_l, f1);

 f2_l := subs(x=x_l, f2);


 # Schleife, bis beide Intervalle abschnittsweise durchlaufen

 x_r := min(x1_r, x2_r);

 while (x_r < infinity) do

   f1_r := subs(x=x_r, f1);

   f2_r := subs(x=x_r, f2);

   if (f1_l <= f2_l) and (f1_r <= f2_r) then

     # f1 <= f2

     pe := [op(pe), f1, x_r];

   elif (f1_l >= f2_l) and (f1_r >= f2_r) then

     # f1 >= f2

     pe := [op(pe), f2, x_r];

   else

     # Schnitt von f1 und f2

     xs := solve(f1=f2, x);

     if (f1_l < f2_l) then

       # linkes Ende: f1 < f2

       pe := [op(pe), f1, xs, f2, x_r];

     else

       # linkes Ende: f1 < f2

       pe := [op(pe), f2, xs, f1, x_r];

     end if;

   end if;

   # Punkte verschieben:

   if (x1_r < x2_r) then

     # f1 ein Intervall weiter:

     i1 := i1 + 1;

     x1_r := `if`(nops(pl1)<2*i1, infinity, pl1[2*i1]);

     f1 := pl1[2*i1-1];

     f1_l := f1_r;

   else

     # f1 ein Intervall weiter:

     i2 := i2 + 1;

     x2_r := `if`(nops(pl2)<2*i2, infinity, pl2[2*i2]);

     f2 := pl2[2*i2-1];

     f2_l := f2_r;

   end if;

   x_r := min(x1_r, x2_r);  

 end do:

 f1_r := subs(x=x_r, f1);

 f2_r := subs(x=x_r, f2);

 return [op(pe), `if`(f1_r < f2_r, f1, f2)];

end proc:

>

Die Funktion  liefert uns eine Liste von St?tzstellen einer piecewise-Funktion  mit zugeh?rigen Funktionen dazwischen

> convert(piecewise(x<-1, 0, x<0, x+1, x<1, 1-x, 1), pwlist, x);

[0, -1, 1+x, 0, 1-x, 1, 1]

Aufstellen der Ergebnismenge f?r die Steuergr??e:

> fuzzifyInference := proc(phi_pkt::numeric, Zeilen::list, phi::numeric, Spalten::list, Regeln::Matrix)
 #option trace;

 local l1, l2, i1, i2, erg;

 l1 := evalf([seq(unapply(Zeilen[i], x)(phi_pkt), i=1..nops(Zeilen))]);

 l2 := evalf([seq(unapply(Spalten[i], x)(phi), i=1..nops(Spalten))]);

 erg := [0];

 for i1 from 1 to nops(Zeilen) do

   for i2 from 1 to nops(Spalten) do

     if l1[i1]<>0 and l2[i2]<>0 then

       erg := MAX(erg, MIN([min(l1[i1],l2[i2])],convert(Regeln[i1,i2], pwlist, x)));

     end if;

   end do:

 end do:

 if nops(erg) = 1 then

   erg := erg[1]; # konstante Funktion

 else

   erg := piecewise(seq(op([x<erg[2*i],erg[2*i-1]]),i=1..(nops(erg)-1)/2),erg[nops(erg)]):

 end if;

 return erg;

end proc:

Und einmal anwenden:

> st := time():
erg := fuzzifyInference(phi_pkt, Zeilen, phi, Spalten, Regeln);

time()-st;

defuzzify(erg, u_range);

time()-st;

erg := PIECEWISE([0, x < -7.5], [.9999999998+.1333333333*x, x < -3.202816537], [.5729577950, x < 0], [.5729577950, x < 0], [.5729577950, x < 3.202816537], [.5729577950, x < 3.202816537], [.9999999998-...

0.56e-1

3.406420661

0.96e-1

Wieder zeichnen wir die Zugeh?rigkeitsfunktion der Ergebnismenge - und

> plot(erg, x=u_range, thickness=3);

[Plot]

>

Fuzzy-Regelung

Und jetzt regeln wir:

> u_fuzzy := proc(X)
 local erg;

 erg := fuzzifyInference(X[4], Zeilen, X[3], Spalten, Regeln);

 if erg <> 0 then

   defuzzify(erg, u_range)

 else

   0

 end if;

end proc;

u_fuzzy := proc (X) local erg; erg := fuzzifyInference(X[4], Zeilen, X[3], Spalten, Regeln); if erg <> 0 then defuzzify(erg, u_range) else 0 end if end proc

> simulation(<0,0,evalf(-10/360*2*Pi),evalf(40/360*2*Pi)>, u_fuzzy, 0.05, 50);

[Plot]

Vector[column]([[-.418191810239929907], [-0.196597871132134228e-1], [0.808693398268046231e-2], [-.406166263453275334]])

>

>

>