m = 1.0;
r = 0.1;
l = 1.0;
h = l/2;
J = (1/4)*m*r*r + (1/12)*m*l*l;
g = 9.81;

function f = rechteSeite(t,y)
    phi = y(1,1);
    om  = y(2,1);
    x   = y(3,1);
    vx  = y(4,1);

    f(1,1) = om;

    sp = sin(phi);
    cp = cos(phi);

    F = 0;

    Z = -m*h*h*om*om*cp*sp + m*g*h*sp + F*h*cp;
    N = J + m*h*h*sp*sp;

    f(2,1) = Z/N;

    f(3,1) = vx;
    f(4,1) = F/m;
endfunction

t = linspace(0,30,3000);

y0 = [1,0,0,0]';
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);

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

Code 0-1: Scilab-Programm: invertierendes Pendel ohne Antriebskraft.

Schritt zwei: Zustandsregler für linearisiertes System auslegen, aber auf das Originalsystem anwenden:

m = 1.0;
r = 0.1;
l = 1.0;
h = l/2;
J = (1/4)*m*r*r + (1/12)*m*l*l;
g = 9.81;

//Linearisierung:

// phi.. = (mgh*phi + F*h) / J
// [phi.;om.] = A * [phi;om] + B * F
A = [  0         1
       m*g*h/J   0];

B = [  0
       h/J];

spec(A)

R = ppol(A,B,[-10+%i*10,-10-%i*10])

// - B*R*y => F = -R*y

function f = rechteSeite(t,y)
    phi = y(1,1);
    om  = y(2,1);
    x   = y(3,1);
    vx  = y(4,1);

    f(1,1) = om;

    sp = sin(phi);
    cp = cos(phi);

    //F = 0;

    r1 = R(1);
    r2 = R(2);

    F = -r1*phi - r2*om;

    Z = -m*h*h*om*om*cp*sp + m*g*h*sp + F*h*cp;
    N = J + m*h*h*sp*sp;

    f(2,1) = Z/N;

    f(3,1) = vx;
    f(4,1) = F/m;
endfunction

t = linspace(0,30,3000);

y0 = [1,0,0,0]';
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);

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

Code 0-2: Scilab-Programm: Zustands geregeltes invertierendes Pendel.

Download Scilab-Skripte: 62_Regelungssysteme/01_day_by_day/inversependulum.zip

Bild 0-1: Pendel ohne Kraft.

Bild 0-2: Pendel mit Zustandsregler.

Ergebnisse der Analyse des linearisierten Systems: Eigenwerte bei -7.5594729 und +7.5594729, d.h. instabiles System.

Vorgabe von Zieleigenwerten bei Polvorgabe mit ppol: -10 + 10i und -10 -10i.

8 -- 16|05|2018 -- WEDNESDAY

Ziele heute:

Zustandsregler für Gesamtzustand entwerfen: x==0 UND phi==0 als Ziel.
Invertierendes Pendel mittels Gradientenverfahren optimieren.

Schritte:

Zustandsregler

Fangbereich des Zustandsreglers untersuchen
Zustandsregler für Gesamtzustand mit ppol auslegen
Kraftbegrenzung einführen

Optimierung

Optimierer aus vorlestzter Veranstaltung wieder ins Gedächtnis rufen
Optimierer anpassen
Optimierung mit Kraftbegrenzung für phi_soll==0 implementieren
Optimierung mit Kraftbegrenzung für phi_soll==0 UND x_soll==0 implementieren
Vergleich mit Zustandsregler, insbesondere Fangbereich
Regler auf möglichst großen Fangbereich hin optimieren.

Offene Fragen

Warum funtioniert der Regler nur in einem begrenzten Winkelbereich?
Wie kann ein Regler entwickelt werden, der auch bei einem Startwert von phi0==PI (hängendes Pendel) funktioniert?

Aufgabenliste

a) Prüfen, ob Drift bei x, bei Zustandsregler für phi!=0
b) Erweitern des Zustandsreglers auf Gesamtzustand, d.h. x!=0 & phi!=0.
c) Maximalkraft (Nachbilden F=-r1*y(1,:)-r2*y(2,:)...) von b) ermitteln und als Obergrenze für Optimierung nehmen
d) Optimierer nochmal ansehen
e) Fehlerfunktion formulieren
f) Optimierer für phi!=0 MIT Kraftbegrenzung aus c)
g) Optimierer auch für phi!=0 & x!=0 MIT Kraftbegrenzung aus c)
h) Vergleich mit Zustandsregler, insbesondere Fangbereich.
i) Regler auf möglichst großen Fangbereich hin optimieren.
j) I-Anteil einführen

Studentische Teillösungen

62_Regelungssysteme/01_day_by_day/Teil_a_b_Zwingmann.zip

9 -- 23|05|2018 -- WEDNESDAY

62_Regelungssysteme/01_day_by_day/optipendel.zip

Aufgabe: Steuerung der Schrittweite alpha ergänzen.

Musterlösung, auch mit Kraftbegrenzung:

clear;

//Simulationsmodell

m = 1.0;
r = 0.1;
l = 1.0;
h = l/2;
J = (1/4)*m*r*r + (1/12)*m*l*l;
g = 9.81;

//Linearisierung:

// phi.. = (mgh*phi + F*h) / J
// [phi.;om.] = A * [phi;om] + B * F
A = [  0         1
       m*g*h/J   0];

B = [  0
       h/J];

spec(A)

R = ppol(A,B,[-10+%i*10,-10-%i*10])

// - B*R*y => F = -R*y

function f = rechteSeite(t,y)
    phi = y(1,1);
    om  = y(2,1);
    x   = y(3,1);
    vx  = y(4,1);

    f(1,1) = om;

    sp = sin(phi);
    cp = cos(phi);

    //F = 0;

    r1 = R(1);
    r2 = R(2);

    F = -r1*phi - r2*om;

    if F>100 then
        F = 100;
    end
    if F<-100 then
        F = -100;
    end

    Z = -m*h*h*om*om*cp*sp + m*g*h*sp + F*h*cp;
    N = J + m*h*h*sp*sp;

    f(2,1) = Z/N;

    f(3,1) = vx;
    f(4,1) = F/m;
endfunction

t = linspace(0,1.5,300);

y0 = [1,0,0,0]';
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);

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

//----------------------------------------------------------

//Fehlerfunktion

RR = [44.143333,    3.4333333  ];

function err = fehlerfunktion(RRR)
    R = RRR;
    y  = ode(y0,t0,t,rechteSeite);
    err = sum( abs(y(1,:)) );
endfunction

//Ausgangspunkt und dessen Fehler:

err_aktuell = fehlerfunktion(RR);
disp("Start-Fehler bei Verwendung der ppol-Ergebnisse:");
disp(err_aktuell);

startfehler = err_aktuell;
RRstart     = RR;

//----------------------------------------------------------

//Optimierer:

// Optimierer formulieren:
SCHRITTE = 2000;
alfamat = [0.1,0.1;0.1,0.1]; //Zeile==Paramternummer, Spalte==Richtung

for i=1:1:SCHRITTE

    RRneu = RR;
    //Zufällig wählen, welcher Parameter variiert wird.
    index = grand(1, 1, "uin", 1, 2);    
    //Zufällig Richtung wählen, in die ich den Parameter ändere
    index_richtung = grand(1, 1, "uin", 1, 2);    // [1 oder 2]
    richtung = 2*(index_richtung-1) - 1;    // [-1 oder +1]

    alfa = alfamat(index,index_richtung);

    RRneu(index) = RRneu(index) + richtung*alfa;
    err_neu = fehlerfunktion(RRneu);
    if err_neu < err_aktuell then
            disp(err_neu);

    end

    if    err_neu <= err_aktuell then
        RR = RRneu;
        err_aktuell = err_neu;
        alfamat(index,index_richtung) = alfa *1.01;
    else
        alfamat(index,index_richtung) = alfa *0.99;
    end
end

disp(RR)
err_best = fehlerfunktion(RR)
R = RR;
y  = ode(y0,t0,t,rechteSeite);
plot(t,y(1,:)',t,y(2,:)');

disp("Fehler nach Optimierung:");
disp(err_best);

disp("Fehler VOR Optimierung:");
disp(startfehler);

Code 0-3: Optimierung mit Schrittweitensteuerung.

10 -- 30|05|2018 -- WEDNESDAY

E-Test 8:30

Einfache Modellierung eines Fahrzeugs / Regelung mit N.N. -- nur Kinematik!!!

Cognimem: Neurocar.

Cognimem: Tutorials / Demonstrations

11 -- 06|06|2018 -- WEDNESDAY

siehe auch: Neuronales Netz auf FPGA-Basis ganz unten in day_by_day 69_FPGA/30_day_by_day/ :-)

Ergänzungen

Im bisherigen Verlauf der LV wurde konzentriert auf den Zustandsregler "hin gearbeitet". Einige Inhalte wurden dabei übersprungen. Ohne bei diesen zu sehr in die Tiefe zu gehen, sollen diese im Nachhinein angesprochen werden.

Gutes Stör- / gutes Führungsverhalten:

62_Regelungssysteme/05_Regleroptimierung/07_Stoerverhalten

Dynamische / statische Stabilität:

62_Regelungssysteme/11_Stabilitaet

Ortskurve, Frequenzgang, Nyquist-Kriterium:

62_Regelungssysteme/11_Stabilitaet/02_Nyqusitkriterium

Bode-Diagramm und vereinfachtes Nyquist-Kriterium:

62_Regelungssysteme/11_Stabilitaet/04_Bode

Methode 2 nach Ziegler und Nichols:

62_Regelungssysteme/07_Einstellregeln/03_Methode2

Windup-Effekt

62_Regelungssysteme/11_Stabilitaet/03_Windup

Dämpfungsgrad

62_Regelungssysteme/07_Einstellregeln/05_Daempfungsgrad

Matrizen A,B,C,D
Zustandsregler mit Beobachter

62_Regelungssysteme/09_Beobachter

Adaptiver Regler

62_Regelungssysteme/12_Adaptiv/01_Definition

Methode der kleinsten Quadrate zur Modell-Identifikation

50_Simulationstechnik/05_Parameterindentifikation/01_KleinsteQuadrate

Genetische Algorithmen zur Optimierung

50_Simulationstechnik/07_Genalgorithmus

Regelung nicht-linearer Systeme

62_Regelungssysteme/13_Analyse/02_Nichtlinear

Das Kalman-Filter

62_Regelungssysteme/14_Kalmanfilter

12 -- 13|06|2018 -- WEDNESDAY

Kleinste Quadrate Methode

Gemeinsamer Entwurf einer Aufgabe und dessen Lösung.

//postulierter funktionaler Zusammenhang:
//  a= c0 + c1 * u^2 + c2 * 2^u
//Messwerte:
//   u | a
//   0 | 4
//   1 | 3
//   2 | 2
//   3 | 1

B = [1 0 1;1 1 2;1 4 4;1 9 8]
d = [-4;-3;-2;-1]
A = B'*B
b = B'*d
c = -inv(A)*b
c0=4
c1=0
c2=-0.4
u = [0 1 2 3]
a = c0 + c1*u + c2*2^u
plot(u,a)
a2 = [4 3 2 1];
plot(u,a2,'*')
cc = lsq(B,-d)

Code 0-4: Beispiel aus dem Unterricht (Handrechnung fehlt hier).

Anderes Beispiel, siehe:

50_Simulationstechnik/05_Parameterindentifikation/01_KleinsteQuadrate/02_Beispiel

Beispiel mit Realitätsbezug: Parameter eines Elektromotors bestimmen:

62_Regelungssysteme/08_Polvorgabe