kramann.info
© Guido Kramann

DAY BY DAY zu Simulations- u. Regelungstechnik 2 für 5-MT, WS 2023/24

(EN google-translate)

(PL google-translate)

AKTUELL: Ab KW48, ab 30.11. beide Blöcke Donnerstags, erst Mechatroniklabor ab 09:30Uhr, dann D.5.13.

Übersicht

Die vorliegende Seite stellt den Einstiegspunkt für diese Lehrveranstaltung dar und verzeichnet chronologisch die behandelten Inhalte.
Aber ein großer Teil der Inhalte findet sich nicht direkt hier, sondern die Seite hier verlinkt auf andere Bereiche von kramann.info.
Die Prüfung findet Semester begleitend in elektronischer Form statt (E-Test).

Beim Kurs Regelungstechnik 2 kommen zu den im ersten Kurs behandelten klassinschen Verfahren der Regelungstechnik Optimierungsverfahren (statisch / kleinste Quadrate Methode, dynamisch / modifiziertes Gradientenverfahren) hinzu, sowie Beispiele aus dem Bereich Softcomputing (Evolutionäre Algorithmen, Fuzzy-Regelung, mit Einschränkungen Neuronale Netze). Die im ersten Kurs behandelten Verfahren bilden eine Art Referenz im Hintergrund, gegen die die neuen Verfahren getestet und diskutiert werden. Deshalb sollen die alten Inhalte in den ersten Stunden auch erst einmal wieder wach gerufen werden.

Die Inhalte der vorangehenden Vorlesung können hier nachgeschlagen werden:

04_SoSe2023/02_RTS_day_by_day

Chronologisches Verzeichnis der im Verlauf des Semesters behandelten Themen

Rückblick
Statische Optimierung: Die Methode der Kleinsten Quadrate
Dynamische Optimierung: Das Gradientenverfahren
Evolutionäre Algorithmen
Fuzzy-Logik und Fuzzy-Regler
Reglerauslegung für nicht lineare dynamische Systeme

LV#1 Mittwoch 22.11.2023

1. Rückblick

Rückmeldung zu den Praktika -- Wie hat sich die Sicht auf Beruf und Studium verändert?

Wiederholung der Inhalte aus Regelungstechnik 1

Scilab (ode, csim)
Laplacetransformation
Auslegungsmethoden nach Ziegler und Nichols
Eigenwertberechnung
Polvorgabe
Zustandsregler

zu 4. Eigenwertberechnung

s. Tafel im Unterricht

zu 5. Polvorgabe

62_Regelungssysteme/08_Polvorgabe

zu 6. Zustandsregler

Theorie Zustandsregler mit Beobachter und Blockbezeichnungen A,B,C,D: 62_Regelungssysteme/09_Beobachter

Übungen zu Laplace als Wiederholung

Übung 1 -- Wandeln Sie jeweils um in den Zeit- bzw. Laplacebereich

1a)

$ \ddot x + 3 \dot x + x = 4 u $

Formel 0-1: Transformieren Sie die lineare Differentialgleichung in den Laplacebereich.

1b)

$ X s^3 + X s - 4 U s^2 = U $

Formel 0-2: Transformieren Sie in den Zeitbereich.

1c)

$ G\left(s\right) = \frac {s^2+1}{s^3+s} $

Formel 0-3: Transformieren Sie in den Zeitbereich. Eingang sei U & Ausgang sei X.

1d)

$ G\left(s\right) = \frac {s^2+1}{s^3+s} $

Formel 0-4: Transformieren Sie in den Zeitbereich. Eingang sei U & Ausgang sei X.

1e) Nachfolgend sind die Laplace-Transformierten einer Regelstrecke G(s) und eines Reglers R(s) dargestellt.

Bilden Sie das Gesamtübertragungsverhalten des offenen und des geschlossenen Regelkreises.
Transformieren Sie beide danach in den Zeitbereich.
Welcher Typ von Regler liegt hier vor?

$ G\left(s\right) = \frac {1}{s^2+5s} $

Formel 0-5: Regelstrecke

$ R\left(s\right) = \frac {1}{s} + 5 + s $

Formel 0-6: Regler

Übung 2 -- Umsetzung regelungstechnischer Systeme mit Scilab

Nehmen Sie sich erneut das System von 1e) vor.

2a)

Übertragen Sie R(s) und G(s) mit Hilfe von syslin() nach Scilab.
Ergänzten Sie Ihr Scilab-Skript so, dass aus R(s) und G(s) das Übertragungsverhalten H(s) des offenen und des geschlossenen Regelkreises Q(s) gebildet werden.

2b)

Vervollständigen Sie zwei Varianten des Scilab-Skripts aus 2a) so, dass ein
Einheitssprung auf den offenen und den geschlossenen Regelkreis simuliert wird.

2c)

Verdeutlichen Sie sich durch Zeichnen eines Blockschaltbildes:

Was ist jeweils Ein- bzw. Ausgang beim offenen und geschlossenen Regelkreis?

Lineare Regression
Kleinste Quadrate Methode

50_Simulationstechnik/05_Parameterindentifikation/01_KleinsteQuadrate/01_Linearregression

50_Simulationstechnik/05_Parameterindentifikation/01_KleinsteQuadrate

50_Simulationstechnik/05_Parameterindentifikation/01_KleinsteQuadrate/02_Beispiel

Allgemeines zu Optimierungsverfahren

50_Simulationstechnik/06_Optimierung

B = [ 1  0
      1  16
      1  25 ];

d = [  -100
       -70
       -50];

A = B'*B

b = B'*d

//  Ac+b=0
//   c  = inv(A)*(-b)

c=inv(A)*(-b)

//...oder direkt mit der Scilab-Funktion lsq:
cc=lsq(B,-d)

Code 0-1: Aufgabe aus dem Unterricht mit Scilab gelöst.

Übung

VoltCm = [ 2.5 20
           2.0 30
           1.5 40
           1.05 60
           0.75 90
           0.5  130
           0.45 150
         ];

// Aufgabe:
// a) Messwerte möglichst gut durch quadratische Funktion annähern
// b) Messwerte möglichst gut durch kubische Funktion annähern
// Verläufe aus a) , b) und den Messwerten übereinander als Plot darstellen

// a) und b) sind mit Scilab mittels kleinster Quadrate-Methode zu lösen.

Code 0-2: Näherung der Meßwerte zu einem Entfernungssensor über Funktionen

Quelle / Datenblatt des Sensors: https://www.sparkfun.com/datasheets/Sensors/Infrared/gp2y0a02yk_e.pdf

LV#2 Donnerstag 23.11.2023

ACHTUNG: ab 30.11. beginnt die LV Donnerstags immer um 12:20 im Raum D.5.13!

Themen:

Besprechung der gestrigen Übung zur Kleinsten Quadrate Methode
Beispiel zur Kleinsten Quadrate Methode mit einer transzendenten Funktion
Heuristiken
Gradientenverfahren

1. Besprechung der gestrigen Übung zur Kleinsten Quadrate Methode

y=[20; 30; 40; 60; 90; 130; 150] //Entfernungen in cm

x=[2.5; 2; 1.5; 1.05; 0.75; 0.5; 0.45 ] //gemessenen Spannungen am Sensor

//Ansatzfunktion:
// y = c0 + c1*x + c2*x^2
// B = [1, xi , xi.*xi]

zs = size(x);
zeilen = zs(1);

B = [ones(zeilen,1),x,(x.*x)]

//B= [1,2.5,6.25; 1,2,4; 1,1.5,2.25; 1,1.05,1.025; 1,0.75,0.5629; 1,0.5,0.25; 1,0.45,0.2025] 

d=-y;

A=B'*B

b =B'*d
c=inv(A)*(-b)

// oder
cc=lsq(B, -d)
plot(x,y,'* ') //Plot der Messpunkte

xx = 0.3:0.1:2.5;
yy = c(1) + c(2).*xx + c(3).*xx.*xx;
yy2 = cc(1) + cc(2).*xx + cc(3).*xx.*xx;

plot(xx,yy,'red'); //Lösung mit Matrizen
plot(xx,yy2,'blu--');//Lösung mit lsq (identisches Ergebnis)

Code 0-3: Studentische Lösung

Bild 0-1: Plot der Lösung.

x=[0,1,2]
y=[2,3,4]

//Ansatzfunktion: y = c0 + c1*2^x

B = [ones(3,1),(2^x)']
d = -y'

A=B'*B
b=B'*d

c=inv(A)*(-b)

plot(x,y,'* ')

xx = -1:0.1:5;
yy = c(1)+c(2)*(2^xx);
plot(xx,yy);

Code 0-4: Kleines Beispiel mit transzendenter Ansatzfunktion, statt Polynom.

Bild 0-2: Plot dazu.

3. Heuristiken

"Heuristik ... bezeichnet Methoden, die mit begrenztem Wissen (unvollständigen Informationen) und wenig Zeit dennoch zu wahrscheinlichen Aussagen oder praktikablen Lösungen kommen. "

Quelle: https://de.wikipedia.org/wiki/Heuristik, Aufruf 23.11.2023.

62_Regelungssysteme/05_Regleroptimierung/02_Heuristiken

62_Regelungssysteme/05_Regleroptimierung/01_Guetefunktion

4. Gradientenverfahren

62_Regelungssysteme/05_Regleroptimierung/04_Gradientenverfahren

62_Regelungssysteme/05_Regleroptimierung/05_ModifizierteG

function F = fehlerfunktion(u)
   x = u(1);
   y = u(2);
   F = x^3-4*x+y^3-16*y
endfunction

ustart = [1,1];
uopt   = [2/sqrt(3),4/sqrt(3) ]

Fstart = fehlerfunktion(ustart)
Fopt = fehlerfunktion(uopt)

alfa = 0.01;
schritte = 1000
u=ustart;
Fakt=fehlerfunktion(u)
for i=1:schritte
   r = floor(grand(2, 1, "unf", 0, 2.999999999))-1;
   uneu = u + alfa*(r');
   Fneu = fehlerfunktion(uneu);
   if Fneu<Fakt
      Fakt = Fneu;
      u=uneu;
   end
end

disp("Endergebnis der Optimierung:")
disp(u)
disp("Fehler am Ende:")
disp(fehlerfunktion(u))

Code 0-5: Implementierung eines modifizierten Gradientenverfahrens bei bekannter Fehlerfunktion F.

//Alfa steuern
function F = fehlerfunktion(u)
   x = u(1);
   y = u(2);
   F = x^3-4*x+y^3-16*y
endfunction

ustart = [1,1];
uopt   = [2/sqrt(3),4/sqrt(3) ]

Fstart = fehlerfunktion(ustart)
Fopt = fehlerfunktion(uopt)

alfa = 0.1;
alfaunterschwelle=0.003;
schritte = 100
u=ustart;
Fakt=fehlerfunktion(u)
//for i=1:schritte
while alfa>alfaunterschwelle
   r = floor(grand(2, 1, "unf", 0, 2.999999999))-1;
   uneu = u + alfa*(r');
   Fneu = fehlerfunktion(uneu);
   if Fneu<=Fakt //wenn sich nichts ändert, trotzdem wandern!!!
      Fakt = Fneu;
      u=uneu;
      alfa=alfa*1.1;
   else
      alfa=alfa*0.9;
   end
end

disp("Endergebnis der Optimierung:")
disp(u)
disp("Fehler am Ende:")
disp(fehlerfunktion(u))
disp("Alfa am Ende:")
disp(alfa);

Code 0-6: Schrittweitensteuerung und Abbruchkriterium ergänzt.

AKTUELL: Ab KW48, ab 30.11. beide Blöcke Donnerstags, erst Mechatroniklabor ab 09:30Uhr, dann D.5.13.

LV#3 Donnerstag 30.11.2023

Themen:

Wiederholung Invertierendes Pendel
Wiederholung Zustandsregler Invertierendes Pendel
Formulierung einer Fehlerfunktion für die Beurteilung der Regelung des Invertierenden Pendels
Optimierung eines Reglers für das invertierende Pendel mittels Gradientenverfahren
Übung: Implementierung Parameterstudie als Alternative Auslegungsmethode
Zustandsregelung des Gesamtzustandes / Hinzunahme der Position in x-Richtung
Übung: Gradientenverfahren anpassen um Regelparameter zur Ausregelung des Gesamtzustandes zu finden
Einführung in Evolutionäre Algorithmen
Genetische Optimierung zum Erraten eines Binärpatterns
Genetische Optimierung des Reglers für das Invertierende Pendel

Material zum Invertierenden Pendel mit Zustandsregler: 62_Regelungssysteme/98_day_by_day_WS2021_SoSe21

Genetische Optimierung: 50_Simulationstechnik/07_Genalgorithmus

Optimierung der Regelparameter beim Invertierenden Pendel mittels modifiziertem Gradientenverfahren

//Auslegung des Reglers für das invertierende Pendel
//mit Hilfe eines Optimierers
clear();
m = 1;     //kg
g = 9.81;  // m/s^2
l = 1;     // m
h = l/2;   // m
r = 0.1;   // m

J = 0.25*m*r*r + (1/12)*m*l*l; // kg*m^2 

R = [10.0,   0.5];
Ropt = [11.63,   0.7904851];
Fgrenz=1.73;
function f = rechteSeite(t,y)
    phi = y(1,1);
    om  = y(2,1);
    f(1,1) = om;
    FA = -R * [phi;om];
    if FA>Fgrenz then
        FA=Fgrenz;
    end
    if FA<-Fgrenz then
        FA=-Fgrenz;
    end
    N = J + m*h*h*sin(phi)*sin(phi);
    f(2,1) = (-m*h*h*om*om*sin(phi)*cos(phi) + h*m*sin(phi)*g + h*cos(phi)*FA )/N;
endfunction

t = 0:0.01:8; 
y0 = [10.0*%pi/180,0.0]';  //Start in der Nähe der instabilen Ruhelage
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);


FFA = -R(1)*y(1,:)-R(2)*y(2,:);

//plot(t,y(1,:)',t,y(2,:)',t,FFA);


function err = fehlerfunktion(RR)
    R = RR;
    y  = ode(y0,t0,t,rechteSeite);
    err = sqrt(y(1,:)*(y(1,:)'));
endfunction

err = fehlerfunktion(Ropt)
R1 = [10.0,0.5];
err = fehlerfunktion(R1)
R2 = [11.0,0.7];
err = fehlerfunktion(R2)



ustart = R1; //[1,1];
//uopt   = [2/sqrt(3),4/sqrt(3) ]

Fstart = fehlerfunktion(ustart)
//Fopt = fehlerfunktion(uopt)

alfa = 0.1;
alfaunterschwelle=0.003;
schritte = 100
u=ustart;
Fakt=fehlerfunktion(u)
//for i=1:schritte
while alfa>alfaunterschwelle
   r = floor(grand(2, 1, "unf", 0, 2.999999999))-1;
   uneu = u + alfa*(r');
   Fneu = fehlerfunktion(uneu);
   if Fneu<=Fakt //wenn sich nichts ändert, trotzdem wandern!!!
      Fakt = Fneu;
      u=uneu;
      alfa=alfa*1.1;
   else
      alfa=alfa*0.9;
   end
end

disp("Endergebnis der Optimierung:")
disp(u)
disp("Fehler am Ende:")
disp(fehlerfunktion(u))
disp("Alfa am Ende:")
disp(alfa);

R = u;
y  = ode(y0,t0,t,rechteSeite);


FFA = -R(1)*y(1,:)-R(2)*y(2,:);

plot(t,y(1,:)',t,y(2,:)',t,FFA);

Code 0-7: Scilab Script zur Optimierung.

Parameterstudie als alternative Möglichkeit, sinnvolle Parameter zu finden

//Auslegung des Reglers für das invertierende Pendel
//mit Hilfe eines Optimierers
clear();
m = 1;     //kg
g = 9.81;  // m/s^2
l = 1;     // m
h = l/2;   // m
r = 0.1;   // m

// siehe: http://www.kramann.info/54_Kinetik/02_NewtonEuler/01_Traegheitsmomente
J = 0.25*m*r*r + (1/12)*m*l*l; // kg*m^2 

R = [10.0,   0.5];
Ropt = [11.63,   0.7904851];
Fgrenz=1.73;
function f = rechteSeite(t,y)
    phi = y(1,1);
    om  = y(2,1);
    f(1,1) = om;
    FA = -R * [phi;om];
    if FA>Fgrenz then
        FA=Fgrenz;
    end
    if FA<-Fgrenz then
        FA=-Fgrenz;
    end
    N = J + m*h*h*sin(phi)*sin(phi);
    f(2,1) = (-m*h*h*om*om*sin(phi)*cos(phi) + h*m*sin(phi)*g + h*cos(phi)*FA )/N;
endfunction

t = 0:0.01:8; 
y0 = [10.0*%pi/180,0.0]';  //Start in der Nähe der instabilen Ruhelage
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);


FFA = -R(1)*y(1,:)-R(2)*y(2,:);

//plot(t,y(1,:)',t,y(2,:)',t,FFA);


function err = fehlerfunktion(RR)
    R = RR;
    y  = ode(y0,t0,t,rechteSeite);
    err = sqrt(y(1,:)*(y(1,:)'));
endfunction

err = fehlerfunktion(Ropt)
R1 = [10.0,0.5];
err = fehlerfunktion(R1)
R2 = [11.0,0.7];
err = fehlerfunktion(R2)



ustart = R1; //[1,1];
//uopt   = [2/sqrt(3),4/sqrt(3) ]

Fstart = fehlerfunktion(ustart)
//Fopt = fehlerfunktion(uopt)

//  11.428411   0.53358
//  r1          r2
//  Parameterstudie:
//  r1 = [9..13] r2=[0.4..0.9]
err_alt = 99999.99;
r1best=0;
r2best=0;
for r1=20:0.1:30
   for r2=1.0:0.01:2.0
      u=[r1,r2];
      err_neu = fehlerfunktion(u);
      if err_neu<err_alt  
          err_alt = err_neu;
          r1best = r1;
          r2best = r2;
      end
   end
end

disp("Beste Kombination aus der Parameterstudie:");

R = [r1best,r2best]
errbest = err_alt

y  = ode(y0,t0,t,rechteSeite);


FFA = -R(1)*y(1,:)-R(2)*y(2,:);

plot(t,y(1,:)',t,y(2,:)',t,FFA);

Code 0-8: Parameterstudie.

Ausregelung aller Zustände: phi, omega, x, vx

clear();
m = 1;     //kg
g = 9.81;  // m/s^2
l = 1;     // m
h = l/2;   // m
r = 0.1;   // m

// siehe: http://www.kramann.info/54_Kinetik/02_NewtonEuler/01_Traegheitsmomente
J = 0.25*m*r*r + (1/12)*m*l*l; // kg*m^2 

R =  [27.834942,   3.2444592,  -6.9996602,  -4.8997621];

function f = rechteSeite(t,y)
    phi = y(1,1);
    om  = y(2,1);
    x   = y(3,1);
    v   = y(4,1);
    FA = -R * [phi;om;x;v];
    f(1,1) = om;
    N = J + m*h*h*sin(phi)*sin(phi);
    f(2,1) = (-m*h*h*om*om*sin(phi)*cos(phi) + h*m*sin(phi)*g + h*cos(phi)*FA )/N;
    f(3,1) = v;
    f(4,1) = FA/m;
endfunction


t = 0:0.01:10;
y0 = [10.0*%pi/180,0.0,0.0,0.0]';  //Start in der Nähe der instabilen Ruhelage
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);

plot(t,y(1,:),t,y(2,:),t,y(3,:),t,y(4,:));

Code 0-9: Zustandsregelung Gesamtzustand.

4 Zustände optimieren

clear();
m = 1;     //kg
g = 9.81;  // m/s^2
l = 1;     // m
h = l/2;   // m
r = 0.1;   // m

// siehe: http://www.kramann.info/54_Kinetik/02_NewtonEuler/01_Traegheitsmomente
J = 0.25*m*r*r + (1/12)*m*l*l; // kg*m^2 

R =  [27.834942,   3.2444592,  -6.9996602,  -4.8997621];

function f = rechteSeite(t,y)
    phi = y(1,1);
    om  = y(2,1);
    x   = y(3,1);
    v   = y(4,1);
    FA = -R * [phi;om;x;v];
    f(1,1) = om;
    N = J + m*h*h*sin(phi)*sin(phi);
    f(2,1) = (-m*h*h*om*om*sin(phi)*cos(phi) + h*m*sin(phi)*g + h*cos(phi)*FA )/N;
    f(3,1) = v;
    f(4,1) = FA/m;
endfunction

t = 0:0.01:10;
y0 = [10.0*%pi/180,0.0,0.0,0.0]';  //Start in der Nähe der instabilen Ruhelage
t0 = 0;
//y  = ode(y0,t0,t,rechteSeite);

///////////


//FFA = -R(1)*y(1,:)-R(2)*y(2,:);

//plot(t,y(1,:)',t,y(2,:)',t,FFA);


function err = fehlerfunktion(RR)
    R = RR;
    y  = ode(y0,t0,t,rechteSeite);
    err = sqrt(y(1,:)*(y(1,:)'));
endfunction

//err = fehlerfunktion(Ropt)
R1 = [28, 3, -6, -5];
err = fehlerfunktion(R1)
R2 = [26, 2, -5, -4];
err = fehlerfunktion(R2)

//#########################

ustart = R1; //[1,1];
//uopt   = [2/sqrt(3),4/sqrt(3) ]

Fstart = fehlerfunktion(ustart)
//Fopt = fehlerfunktion(uopt)

alfa = 0.1;
alfaunterschwelle=0.003;
schritte = 100
u=ustart;
Fakt=fehlerfunktion(u)
//for i=1:schritte
while alfa>alfaunterschwelle
   r = floor(grand(4, 1, "unf", 0, 2.999999999))-1;
   uneu = u + alfa*(r');
   Fneu = fehlerfunktion(uneu);
   if Fneu<=Fakt //wenn sich nichts ändert, trotzdem wandern!!!
      Fakt = Fneu;
      u=uneu;
      alfa=alfa*1.1;
   else
      alfa=alfa*0.9;
   end
end

disp("Endergebnis der Optimierung:")
disp(u)
disp("Fehler am Ende:")
disp(fehlerfunktion(u))
disp("Alfa am Ende:")
disp(alfa);

R = u;
y  = ode(y0,t0,t,rechteSeite);


//FFA = -R(1)*y(1,:)-R(2)*y(2,:);
///


plot(t,y(1,:),t,y(2,:),t,y(3,:),t,y(4,:));

Code 0-10: 4 Zustände optimieren.

DAY BY DAY zu Simulations- u. Regelungstechnik 2 für 5-MT, WS 2023/24

(EN google-translate)

(PL google-translate)

Übersicht

Chronologisches Verzeichnis der im Verlauf des Semesters behandelten Themen

LV#1 Mittwoch 22.11.2023

1. Rückblick

Rückmeldung zu den Praktika -- Wie hat sich die Sicht auf Beruf und Studium verändert?

Wiederholung der Inhalte aus Regelungstechnik 1

zu 1.: Scilab (ode, csim)

zu 2.: Laplacetransformation

zu 3.: Auslegungsmethoden nach Ziegler und Nichols

zu 4. Eigenwertberechnung

zu 5. Polvorgabe

zu 6. Zustandsregler

Übungen zu Laplace als Wiederholung

Übung 1 -- Wandeln Sie jeweils um in den Zeit- bzw. Laplacebereich

Übung 2 -- Umsetzung regelungstechnischer Systeme mit Scilab

Allgemeines zu Optimierungsverfahren

Übung

LV#2 Donnerstag 23.11.2023

Themen:

1. Besprechung der gestrigen Übung zur Kleinsten Quadrate Methode

3. Heuristiken

4. Gradientenverfahren

LV#3 Donnerstag 30.11.2023

Themen:

Optimierung der Regelparameter beim Invertierenden Pendel mittels modifiziertem Gradientenverfahren

Parameterstudie als alternative Möglichkeit, sinnvolle Parameter zu finden

Ausregelung aller Zustände: phi, omega, x, vx

4 Zustände optimieren