kramann.info
© Guido Kramann

Day by Day

(EN google-translate)

(PL google-translate)

Nachfolgend wird chronologisch verzeichnet, welche Inhalte in den Einzellehrveranstaltungen im Simulations- und Regelungstechnik für den Bachelor Ingenieurwissenschaften der Studienrichtung Mechatronik7 im Wintersemester 2021/22 behandelt wurden.

#1 Mittwoch 24.11.2021

Übersicht und Vereinbarungen zur Ausgestaltung des Unterrichts

64_Regelungssysteme

Themen

Ergänzungen zur klassischen linearen Regelungstechnik
Behandlung nicht linearer Probleme
Optimierungsprobleme statisch / dynamisch
Evolutionäre Algorithmen
Fuzzy-Logik / Fuzzy-Systeme

Statische Optimierung

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

studentische Lösung

clear;
mode(0)

//Messwerte eines Sensors 
VjeCm=[2.5 20;2.0 30;1.5 40;1.05 60;0.75 90;0.5 130;0.45 150]

//a)Messwerte möglichst gut durch quadratische Funktion annähern
//b)Messwerte möglichst gut durch kubische Funktion annähern
//c)Verläufe aus a) und b) zusammen in einem Plot darstellen

//Spannungsmesswerte zu den jeweiligen Entfernungen
a=2.5;
b=2;
c=1.5;
d=1.05;
e=0.75;
f=0.5;
g=0.45;
//a) y=c1*x^0 + c2*x + c3*x^2

B1=[1 a a^2
    1 b b^2
    1 c c^2
    1 d d^2
    1 e e^2
    1 f f^2
    1 g g^2]
d1=[-20;-30;-40;-60;-90;-130;-150]
c1=lsq(B1,-d1)



clf;
V= [0.5:0.1:2.5];
plot(c1(1)+c1(2)*V+(c1(3)*V^2),V);
xlabel("Weg [cm]")
ylabel("Spannung [V]")


//b
B2=[1 a a^2 a^3
    1 b b^2 b^3
    1 c c^2 c^3
    1 d d^2 d^3
    1 e e^2 e^3
    1 f f^2 f^3
    1 g g^2 g^3]
d2=[-20;-30;-40;-60;-90;-130;-150]
c2=lsq(B2,-d2)
plot(c2(1)+c2(2)*V+(c2(3)*V^2)+(c2(4)*V^3),V,'r');


plot(VjeCm(:,2),VjeCm(:,1),'g');
hl=legend(['Quadratisch';'Kubisch';'Messwerte']);

Code 0-3: studentische Lösung

Lösungsvariante

clear;

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

zeilen_spalten = size(VoltCm);
zeilen = zeilen_spalten(1);
spalten = zeilen_spalten(2);

VC = VoltCm(zeilen:-1:1,:)

// a) Annäherung durch quadratische Gleichung

Volt = VC(:,1);
Cm   = VC(:,2);

d = -Cm

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


c = lsq(B,-d)

v = [0.45:0.1:2.5]; //Ordinate für quadratische Funktion
q = c(1)+c(2)*v+c(3)*(v.*v);

plot(Volt,Cm,'gre*')
plot(v,q)

// b) Annäherung durch kubische Gleichung


d = -Cm

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

c = lsq(B,-d)

v = [0.45:0.1:2.5]; //Ordinate für quadratische Funktion
q = c(1)+c(2)*v+c(3)*(v.*v)+c(4)
*(v.*v.*v);

plot(v,q,'red')

Code 0-4: Lösungsvariante

Übung 2

Implementierung von LSQ und Polynomannäherung an Messwerte mit Java / Processing

Ausgangspunkt: Implementierung von Matrizenmultiplikation:

LSQ001.zip

LSQ002.zip -- Ergänzen Sie die Klasse LSQ in diesem Projekt so, dass sie Least Square Fit implementiert.

LSQ002b_LGSloesen.zip -- mit LGS Löser für 2x2 und 3x3

LSQ002c_weitereMethoden.zip -- weitere Matrizen-Methoden

LSQ003_LSQ2x2.zip -- Musterlösung, in der LSQ für A=B'*B, mit A: 2x2-Matrix implementiert ist.

LSQ004_LSQ3x3.zip -- ... für A: 3x3 Matrix

Quelltext zu LSQ004_LSQ3x3

/**
   Kleinste Quadrate-Methode für A: 2x2-Matrix 
*/
LSQ la = new LSQ();

public void setup()
{
     double[][] A = {
                       {1,2,3},
                       {4,5,6}
                    };
     double[][] B = {
                       {1,2},
                       {3,4},
                       {5,6}
                    };
                    
    double[][] C = la.mult(A,B); 
    la.showMatrix(C);
    
    println("lgs2 testen:");
    
    /*
          x=2, y=1
          3*x+4*y=10
          1*x+5*y=7
          
          A*x=b
          
          A=[3 4
             1 5]
          b=[10
             7]
    */
    
    double[] xy = la.lgs2(new double[][] {{3,4},{1,5}},new double[] {10,7});
    la.showVektor(xy);
    
    double[] xyz = la.lgs3(new double[][] {{-1,2,3},{4,5,6},{7,8,-9}},new double[] {4,12,-36});
    la.showVektor(xyz);
    
    println("Testen von LSQ für A: 2x2 Matrix: (Beispiel aus LV)");
    
    double[][] BBB = {
                       {1,0},
                       {1,16},
                       {1,25}
                   };
    double[] ddd = {-100,-70,-50};
    
    double[] c = la.lsq2(BBB,ddd);
    println("Parameter c0 und c1:");
    la.showVektor(c);
    
    //######################################
    //Sensor-Beispiel bei quadratischer Näherungskurve:
    double[][] BBBB ={
   {1.0,   0.45,   0.2025},
   {1.0,   0.5,    0.25  },
   {1.0,   0.75,   0.5625},
   {1.0,   1.05,   1.1025},
   {1.0,   1.5 ,   2.25  },
   {1.0,   2.0,     4.    },
   {1.0,   2.5 ,   6.25  }};
   
   
    double[] dddd = {  -150.0,
  -130.0,
  -90.0 ,
  -60.0 ,
  -40.0 ,
  -30.0 ,
  -20.0 };
    
    double[] cc = la.lsq3(BBBB,dddd);
    println("Parameter c0, c1 und c2: Centimenter = c0+c1*Volt+c2*Volt^2");
    la.showVektor(cc);
    
}

public void draw()
{
}

Code 0-5: Haupt-Tab im Sketch

public class LSQ extends LineareAlgebra
{
    /**
        Beschränkte Implementierung:
        A = B'*B mit A:2x2 Matrix als Bedingung
        entsprechend:
        b = B'*d mit b: 2x1 Vektor
    */
    public double[] lsq2(double[][] B, double[] d)
    {
        double[][]  A = mult(trans(B),B);
        double[][] dd = new double[d.length][1];
        for(int i=0;i<dd.length;i++)
            dd[i][0]=d[i];
        double[][]  b = mult(trans(B),dd);
        
        double[] bb = {-b[0][0],-b[1][0]}; // wegen: A*c+b=0, gilt: A*c=-b
        
        showVektor(size(A));
        
        return lgs2(A,bb);
    }
    
    public double[] lsq3(double[][] B, double[] d)
    {
        double[][]  A = mult(trans(B),B);
        double[][] dd = new double[d.length][1];
        for(int i=0;i<dd.length;i++)
            dd[i][0]=d[i];
        double[][]  b = mult(trans(B),dd);
        
        double[] bb = {-b[0][0],-b[1][0],-b[2][0]}; // wegen: A*c+b=0, gilt: A*c=-b
        
        showVektor(size(A));
        
        return lgs3(A,bb);
    }
    
}

Code 0-6: Klasse LSQ

public class LineareAlgebra
{
     /**
          A = [ 1 2 3
                4 5 6]
          B = [ 1
                2
                3 ]
          (Zeilenanzahl von B == Spaltenanzahl von A)
          C = A*B
     */
     public double[][] mult(double[][] A, double[][] B)
     {
          if(B.length == A[0].length)
          {
               //korrekt, Berechnung ausführen:
               double[][] C = new double[A.length][B[0].length];
               
               for(int i=0;i<C.length;i++)
               {
                     for(int k=0;k<C[i].length;k++)
                     {
                           double h=0.0;
                           for(int p=0;p<A[0].length;p++)
                           {
                                h += A[i][p]*B[p][k];
                           }
                           C[i][k] = h;
                     }
               }
               
               return C;
          }
          else
          {
               println("ERROR in LineareAlgebra.mult: Zeilenanzahl von B != Spaltenanzahl von A");
               return null;
          }
     }
     
     /**
         2x2 Determinante
         A: 2x2 Matrix
     */
     public double det2(double[][] A)
     {
           return A[0][0]*A[1][1] - A[1][0]*A[0][1];
     }
     
     /**
         Löse 2x2 LGS:
         A*x=b
         
         A: 2x2 Matrix
         b: 2x1 Vektor
     */
     public double[] lgs2(double[][] A, double[] b)
     {
          double D = det2(A);
          if(D==0.0)
          {
              println("ERROR in LineareAlgebra.lgs2(..): Determinante ist Null");
              return null;
          }
          double[][] Dx = {
                             {b[0],A[0][1]},
                             {b[1],A[1][1]}
                          };
          double[][] Dy = {
                             {A[0][0],b[0]},
                             {A[1][0],b[1]}
                          };
          return new double[] {  det2(Dx)/det2(A),det2(Dy)/det2(A)  };                
     }

     //########################################
     
     /**
         3x3 Determinante
         A: 3x3 Matrix
     */
     public double det3(double[][] A)
     {
           return A[0][0]*A[1][1]*A[2][2] 
                 +A[0][1]*A[1][2]*A[2][0]
                 +A[0][2]*A[1][0]*A[2][1]
                 
                 -A[2][0]*A[1][1]*A[0][2]
                 -A[2][1]*A[1][2]*A[0][0]
                 -A[2][2]*A[1][0]*A[0][1];
     }
     
     /**
         Löse 3x3 LGS:
         A*x=b
         
         A: 3x3 Matrix
         b: 3x1 Vektor
     */
     public double[] lgs3(double[][] A, double[] b)
     {
          double D = det3(A);
          if(D==0.0)
          {
              println("ERROR in LineareAlgebra.lgs2(..): Determinante ist Null");
              return null;
          }
          double[][] Dx = {
                             {b[0],A[0][1],A[0][2]},
                             {b[1],A[1][1],A[1][2]},
                             {b[2],A[2][1],A[2][2]}
                          };
          double[][] Dy = {
                             {A[0][0],b[0],A[0][2]},
                             {A[1][0],b[1],A[1][2]},
                             {A[2][0],b[2],A[2][2]}
                          };
          double[][] Dz = {
                             {A[0][0],A[0][1],b[0]},
                             {A[1][0],A[1][1],b[1]},
                             {A[2][0],A[2][1],b[2]}
                          };
          return new double[] {  det3(Dx)/det3(A),det3(Dy)/det3(A),det3(Dz)/det3(A)  };                
     }
     
     
     //########################################
     
     //weitere Methoden:
     public double[] plus(double[] a, double[] b)
     {
          double[] c = new double[a.length];
          for(int i=0;i<c.length;i++)
              c[i] = a[i]+b[i];
          return c;    
     }
     public double[][] plus(double[][] a, double[][] b)
     {
          double[][] c = new double[a.length][a[0].length];
          for(int i=0;i<c.length;i++)
            for(int k=0;k<c[0].length;k++)
              c[i][k] = a[i][k]+b[i][k];
          return c;    
     }
     
     /**
         Element für Element, statt Matrizenmultiplikation!
     */
     public double[] dotMult(double[] a, double[] b)
     {
          double[] c = new double[a.length];
          for(int i=0;i<c.length;i++)
              c[i] = a[i]*b[i];
          return c;    
     }
     public double[][] dotMult(double[][] a, double[][] b)
     {
          double[][] c = new double[a.length][a[0].length];
          for(int i=0;i<c.length;i++)
            for(int k=0;k<c[0].length;k++)
              c[i][k] = a[i][k]*b[i][k];
          return c;    
     }
     
     /**
        Matrix transponieren
     */
     public double[][] trans(double[][] m)
     {
          double[][] n = new double[m[0].length][m.length];
          for(int i=0;i<m.length;i++)
              for(int k=0;k<m[0].length;k++)
                   n[k][i] = m[i][k];
          return n;         
     }
     
     public void showMatrix(double[][] m)
     {
          for(int i=0;i<m.length;i++)
          {
               for(int k=0;k<m[i].length;k++)
               {
                    print(m[i][k]+" ");
               }
               println();
          }
     }
     public void showMatrix(int[][] m)
     {
          for(int i=0;i<m.length;i++)
          {
               for(int k=0;k<m[i].length;k++)
               {
                    print(m[i][k]+" ");
               }
               println();
          }
     }
     
     public int[] size(double[][] m)
     {
           return new int[] {m.length,m[0].length};
     }
     public int size(double[] v)
     {
           return v.length;
     }
     
     public void showVektor(double[] m)
     {
               for(int k=0;k<m.length;k++)
               {
                    print(m[k]+" ");
               }
               println();
     }
     public void showVektor(int[] m)
     {
               for(int k=0;k<m.length;k++)
               {
                    print(m[k]+" ");
               }
               println();
     }
}

Code 0-7: Klasse LineareAlgebra

#2 Mittwoch 01.12.2021

Thema: Fortsetzung mit Optimierung, jetzt evolutionäre Algorithmen

Zwar würden klassischerweise nun Gradientenverfahren besprochen werden. Diese sollen angerissen, aber später ausführlich ausgeführt werden zugunsten genetischer Algorithmen, da diese im Zusammenhang mit Neuronalen Netzen morgen in Informatik-Ergänzungen gebraucht werden.

Übung

Das Prinzip der genetischen Optimierung soll an einem einfachen Anwendungsfall durch praktische Implementierung verständlicher werden.

Zudem soll die Anwendung von Objektorientierung erübt werden.

"Lottozahlen" 7 - 14 - 15 - 29 - 36 - 37 sollen durch einen Genalgorithmus gefunden werden. Jedes Gen besitzt als Attribut ein Array aus sechs Integer-Zahlen, die zu Beginn zufällig gesetzt werden [0..100]. Folgende Fehlerfunktion wird genutzt, um zu bestimmen, wie schlecht ein Gen ist:


public double bestimmeFehler(Gen gen)
{
    int[] lotto = {7,14,15,29,36,37};

    double fehler = 0.0;

    for(int index=0;index<lotto.length;index++)
    {
        double diff = (double)(lotto[index] - gen.get(index));
        fehler += diff*diff;
    }
    return Math.sqrt(fehler); //Euklidische Vektornorm als Maß für den Fehler.
}

Code 0-8: Fehlerfunktion.

Hinweis: In diesem Anwendungsbeispiel sind Geno- und Phänotyp indentisch. Das ist normalerweise nicht der Fall, vergl. Vorlesung.

Die Aufgabe:

Arbeiten Sie am besten in Zweiergruppen.
Entwickeln Sie ein objektorientiertes System für genetische Optimierung.
Gen-Objekte beinhalten Integer-Arrays als "Gene", deren Länge und Belegungs-Intervall durch Aufruf des Konstruktors festgelegt wird.
Obige Fehlerfunktion soll benutzt werden, um als Anwendungsfall mit den Genen möglichst nahe an die Lottozahlen zu kommen.
(Der Anwendungsfall ist nicht sehr sinnvoll, aber dafür einfach.)
Sorgen Sie dafür, dass der Verlauf der Optimierung sichtbar wird.

Gehen Sie die Übungsaufgabe in Ruhe durch und sammeln Sie Fragen, die wir zunächst klären bevor Sie mit der Bearbeitung beginnen.

Vorlage für die prozedurale Umsetzung:

Genetisch001.zip

import java.util.Random;

Random zufall = new Random(System.currentTimeMillis());

int[][] gen = new int[100][6];
int[][] bestgen = new int[10][6];
double[] fehler = new double[gen.length];
public void setup()
{
    for(int i=0;i<gen.length;i++)
        for(int k=0;k<gen[i].length;k++)
            gen[i][k] = zufall.nextInt(100); // Vergabe von Zufallszahlen von 0...99
    //...
    for(int DURCHGANG=0; DURCHGANG<100;DURCHGANG++)
    {
          //Evolutionärer Vorgang:
          
          // 1. Bewerten, wie nah die 
          //    zufälligen Zahlen in
          //    gen[i] an den Lottozahlen
          //    dran sind:
          for(int i=0;i<gen.length;i++)
             fehler[i] = bestimmeFehler(gen[i]);
          
          // 2. Die besten 10 müßten aus den 100
          //    gefunden werden, also die Gene
          //    mit den 10 kleinsten Fehlern
          //
          //    Diese würde man dann in bestgen kopieren.
          
          // ....
          
          // 3. Immer zwei Beste zufällig auswählen und
          //    immer im Wechsel mal eine Zahl von dem
          //    einen, mal von dem anderen eine Zahl nehmen.
          // Allerdings: 10 Besten behalte ich.
          // 0..9: sind die Besten dann
          // 10..99: sind dann die kombinierten neu gebildeten.
          
          // 4. Zufällige Mutationen bei den hinteren 90
          //    durchführen
          
          // JETZT geht es von Vorne los,
          // mit der neu gebildeten Generation
          
    }
    
    //...
}

public void draw()
{
}

Code 0-9: Vorlage für die prozedurale Umsetzung

public double bestimmeFehler(int[] gen)
{
    int[] lotto = {7,14,15,29,36,37};

    double fehler = 0.0;

    for(int index=0;index<lotto.length;index++)
    {
        double diff = (double)(lotto[index] - gen[index]);
        fehler += diff*diff;
    }
    return Math.sqrt(fehler); //Euklidische Vektornorm als Maß für den Fehler.
}

Code 0-10: Vorlage für die prozedurale Umsetzung - Tab mit Fehlerfunktion

Genetisch002.zip -- Studentische Lösung für das Finden der 10 besten.

Genetisch003.zip -- Lösung für prozedurale Umsetzung des genetischen Algorizthmus'.

Genetisch004.zip -- OOP Lösung, unvollständig

Genetisch005.zip -- OOP Lösung fertig

#3 Mittwoch 08.12.2021

Themen:

Gradientenverfahren
Modifiziertes Gradientenverfahren
Scilab

Heute lernen Sie den klassischen Zugang zur Optimierung in Form von Gradientenverfahren kennen. Ergänzend werden wir die Fertigkeiten bei der Erstellung von Skripten in Scilab vertiefen.

Heuristiken als Grundlage für die Implementierung von Optimierungsverfahren

50_Simulationstechnik/06_Optimierung/01_Gradientenverfahren -- Theorie zu Gradientenverfahren

50_Simulationstechnik/06_Optimierung/02_Heuristiken -- Voraussetzungen dafür, dass ein Optimierer besser wird als bei "trial and error".

Prinzip des modifizierten Gradientenverfahrens

50_Simulationstechnik/06_Optimierung/03_ModifizierteG -- Beispiel für ein modifiziertes Gradientenverfahren

50_Simulationstechnik/06_Optimierung/04_optim -- Verwendung von optim (Scilab-Funktion)

62_Regelungssysteme/01_day_by_day -- verbesserte Variante eines selbst geschriebenen Algorithmus für ein modifiziertes Gradientenverfahren (siehe Code 0-9: optimierer.sce )

Übung

Ein linearer Schwinger ohne Dämpfung, m=1kg, C=1N/m, soll mittels eines PID-Reglers aktiv beruhigt werden.

Bilden Sie Zweiergruppen, innerhalb derer Sie die nachfolgenden Teilaufgaben bearbeiten.
Lesen Sie zunächst die nachfolgenden Aufgabenteile genau duch und fragen bei der Lehrperson nach, wenn es Verständnisprobleme gibt, oder Zweifel an der Machbarkeit.
Formulieren Sie ein entsprechendes Simulations-System mit Scilab.
Überlegen Sie sich eine Fehlerfunktion, die durch eine einzige Zahl quantifiziert, wie schlecht der Regler ist (integrales Kriterium).
Formulieren Sie in Scilab ein modifiziertes Gradientenverfahren, das den Regler optimiert.
Vergleichen Sie die Qualität der gefundenen Regler Gruppen übergreifend untereiander.

Während des Unterrichts entstandene Lösungen

clear;

function r = abstieg(x,y)
    grad = [3*x*x - 4;3*y*y - 16];
    normGrad = sqrt(grad'*grad);

    r = - grad ./ normGrad

endfunction

xsoll = 2/sqrt(3);
ysoll = 4/sqrt(3);

x=4
y=4

alpha = 0.1 //Güte der Lösung hängt von alpha ab.
for i=1:100 //Güte der Lösung hängt von Iterationsschritten ab.
    r = abstieg(x,y);
    x = x + alpha*r(1)
    y = y + alpha*r(2)
    err = sqrt(  (x-xsoll)*(x-xsoll) + (y-ysoll)*(y-ysoll)  )
end

Code 0-11: Modifiziertes Gradientenverfahren am Beispiel von F(x,y) = x^3-4x+y^3-16y, iterative Methode des größten Abstiegs

Erschwerte aber realistische Situation: Kein Gradient von F(x,y) verfügbar!!!

clear;
//erschwerte aber realistische Situation:
//Kein Gradient von F(x,y) verfügbar!!!
function F = berechneFehler(x,y)
    F = x^3-4*x+y^3-16*y;
endfunction

xsoll = 2/sqrt(3);
ysoll = 4/sqrt(3);

x=4
y=4



alpha = 0.01 //Güte der Lösung hängt von alpha ab.
for i=1:1000 //Güte der Lösung hängt von Iterationsschritten ab.
    F = berechneFehler(x,y);
    Fxp = berechneFehler(x+alpha,y);
    Fxn = berechneFehler(x-alpha,y);
    Fyp = berechneFehler(x,y+alpha);
    Fyn = berechneFehler(x,y-alpha);

    if (Fxp<=Fxn & Fxp<=Fyp & Fxp<=Fyn) then
         x=x+alpha;
    elseif (Fxn<=Fxp & Fxn<=Fyp & Fxn<=Fyn) then
         x=x-alpha;
    elseif (Fyp<=Fxp & Fyp<=Fxn & Fyp<=Fyn) then
         y=y+alpha;
    else    
         y=y-alpha;
    end

   err = sqrt(  (x-xsoll)*(x-xsoll) + (y-ysoll)*(y-ysoll)  )
end

Code 0-12: Erschwerte aber realistische Situation: Kein Gradient von F(x,y) verfügbar!!!

Studentische Lösung Schrittweitensteuerung UND Abbruchkriterium

clear;
//erschwerte aber realistische Situation:
//Kein Gradient von F(x,y) verfügbar!!!
function F = berechneFehler(x,y)
    F = x^3-4*x+y^3-16*y;
endfunction

xsoll = 2/sqrt(3);
ysoll = 4/sqrt(3);

x=4
y=4
/*
    Wenn alle Fehler Fxn,Fxp,Fyp,Fyn <=F sind, dann soll alpha verkleinert werden
    SONST soll alpha vergroessert werden
*/
/*
    1. Idee fuer Abbruchkriterium: Unterschreiten einer Fehlerschwelle von F
    2. Idee fuer Abbruchkriteirum: Unterschreiten eines minimalen Alpha-Wertes
*/


alpha = 0.01 //Güte der Lösung hängt von alpha ab.
for i=1:1000 //Güte der Lösung hängt von Iterationsschritten ab.
    F = berechneFehler(x,y);
    //In 4 richtungen laufern und kleinsten Fehler suchen
    Fxp = berechneFehler(x+alpha,y);
    Fxn = berechneFehler(x-alpha,y);
    Fyp = berechneFehler(x,y+alpha);
    Fyn = berechneFehler(x,y-alpha);
    //dann in die jeweilige Richtung weiter laufen
    if (Fxp<=Fxn & Fxp<=Fyp & Fxp<=Fyn) then
         x=x+alpha;
    elseif (Fxn<=Fxp & Fxn<=Fyp & Fxn<=Fyn) then
         x=x-alpha;
    elseif (Fyp<=Fxp & Fyp<=Fxn & Fyp<=Fyn) then
         y=y+alpha;
    else    
         y=y-alpha;
    end

    if (Fxp >= F & Fxn >= F & Fyp >= F & Fyn >= F) then
        alpha = alpha/2;
    //else
    //    alpha=alpha*2;
    end
    
    F_alt(i)=F;
    if(i >= 2) then
          Vergleich=F_alt(i-1)/F
             if (Vergleich == 1.0) then
//                i=999;
                  break;
             end
    end

    
    
   err(i) = sqrt(  (x-xsoll)*(x-xsoll) + (y-ysoll)*(y-ysoll)  );
end

Code 0-13: Studentische Lösung Schrittweitensteuerung UND Abbruchkriterium

Neuer Ausgangspunkt mit zufälliger Richtungswahl

clear;
//erschwerte aber realistische Situation:
//Kein Gradient von F(x,y) verfügbar!!!
function F = berechneFehler(x,y)
    F = x^3-4*x+y^3-16*y;
endfunction

xsoll = 2/sqrt(3);
ysoll = 4/sqrt(3);

x=4
y=4

alpha = 0.01 //Güte der Lösung hängt von alpha ab.
alpha_min = 0.0001 //Güte der Lösung hängt von alpha ab.

RICHTUNGEN = [ 1, 0;  -1, 0;  0,1;  0,-1];

for i=1:10000
 
    F = berechneFehler(x,y);

    index = 1 + floor(grand(1, 1, "unf", 0, 3.999999999));  //liefert 1..4

    aktuelleRichtung = RICHTUNGEN(index,:);

    Fneu = berechneFehler(x+alpha*aktuelleRichtung(1),y+alpha*aktuelleRichtung(2));

    if (Fneu<=F) then
        x = x + alpha*aktuelleRichtung(1);
        y = y + alpha*aktuelleRichtung(2);
    end

    if (Fneu>F) then
        alpha = alpha*0.9999;
    else
        alpha = alpha*1.0001;
    end

    err = sqrt(  (x-xsoll)*(x-xsoll) + (y-ysoll)*(y-ysoll)  )    
end

Code 0-14: Neuer Ausgangspunkt mit zufälliger Richtungswahl

Nachtrag: Nachfolgendes Skript besitzt ein funktionierendes Abbruchkriterium

clear;
//erschwerte aber realistische Situation:
//Kein Gradient von F(x,y) verfügbar!!!
function F = berechneFehler(x,y)
    F = x^3-4*x+y^3-16*y;
endfunction

xsoll = 2/sqrt(3);
ysoll = 4/sqrt(3);

x=4
y=4

alpha = 0.01 //Güte der Lösung hängt von alpha ab.
alpha_min = 0.005 //Güte der Lösung hängt von alpha ab.

RICHTUNGEN = [ 1, 0;  -1, 0;  0,1;  0,-1];
i=1
while %T
   
    F = berechneFehler(x,y);

    index = 1 + floor(grand(1, 1, "unf", 0, 3.999999999));  //liefert 1..4
    disp("Durchlauf:");
    disp(i);
    disp("alpha:");
    disp(alpha);
    aktuelleRichtung = RICHTUNGEN(index,:);

    Fneu = berechneFehler(x+alpha*aktuelleRichtung(1),y+alpha*aktuelleRichtung(2));

    if (Fneu<=F) then
        x = x + alpha*aktuelleRichtung(1);
        y = y + alpha*aktuelleRichtung(2);
    end

    if (Fneu>F) then
        alpha = alpha*0.9999;
    else
        alpha = alpha*1.0001;
    end

    err = sqrt(  (x-xsoll)*(x-xsoll) + (y-ysoll)*(y-ysoll)  )    

    if alpha<=alpha_min then
        break;
    end

    i=i+1;
end

Code 0-15: Skript besitzt ein funktionierendes Abbruchkriterium

Nachtrag 2: Streuung der letzten N Werte für F als Abbruchkriterium benutzen

clear;
//erschwerte aber realistische Situation:
//Kein Gradient von F(x,y) verfügbar!!!
function F = berechneFehler(x,y)
    F = x^3-4*x+y^3-16*y;
endfunction

xsoll = 2/sqrt(3);
ysoll = 4/sqrt(3);

x=4
y=4

alpha = 0.01 //Güte der Lösung hängt von alpha ab.
alpha_min = 0.005 //Güte der Lösung hängt von alpha ab.

RICHTUNGEN = [ 1, 0;  -1, 0;  0,1;  0,-1];
i=1

letzteN = 300;  //Streuung der letzten N Werte benutzen
FEHLER = grand(1, letzteN, "unf", 0, 100.0)

while %T
   
    F = berechneFehler(x,y);

    index = 1 + floor(grand(1, 1, "unf", 0, 3.999999999));  //liefert 1..4
    disp("Durchlauf:");
    disp(i);
    disp("alpha:");
    disp(alpha);
    aktuelleRichtung = RICHTUNGEN(index,:);

    Fneu = berechneFehler(x+alpha*aktuelleRichtung(1),y+alpha*aktuelleRichtung(2));

    if (Fneu<=F) then
        x = x + alpha*aktuelleRichtung(1);
        y = y + alpha*aktuelleRichtung(2);
    end

    if (Fneu>F) then
        alpha = alpha*0.999;
    else
        alpha = alpha*1.001;
    end    

    FEHLER(1+modulo(i,letzteN))=F;
    STREUUNG = stdev(FEHLER);
    disp("STREUUNG=");
    disp(STREUUNG);

    err = sqrt(  (x-xsoll)*(x-xsoll) + (y-ysoll)*(y-ysoll)  )    

    

    if STREUUNG<0.00000001 then
        break;
    end

    i=i+1;
end

Code 0-16: Streuung der letzten N Werte für F als Abbruchkriterium benutzen

Musterlösung zur Übung

clear;

//Regelstrecke

//  m x.. = -Cx + F
//    x.. =  -x + F
//    xs^2+x = F

//               1
//  G(s) =  ---------- 
//           s^2 +  1

//Regler

//  klassisch: R(s) = K*(1+1/(Tn*s) + Tv*s)
//  für Optimierung verworfen, da K auf alle drei Terme wirkt!, deshalb:

//  R(s) = PP + II/s  + DD*s

//          PP*s + II + DD*s^2
//  R(s) = --------------------
//                s

PP = 10.0;
II = 0.1;
DD = 0.5;

s = poly(0,"s"); 
//Übertragungsfunktion der Regelstrecke definieren:
G = syslin('c',[1],[s^2 +  1]);
//Übertragungsfunktion des Reglers definieren:
R = syslin('c',[PP*s + II + DD*s^2],[s]);

H = R*G;
Q = H/(1+H);

t=0:0.01:15;
zs = size(t);
w = ones(1,zs(2));

//y=csim(w,t,Q);
//plot(t,y,'blu--');


function F = berechneFehler(PARAM)
    PP = PARAM(1);
    II = PARAM(2);
    DD = PARAM(3);
    R = syslin('c',[PP*s + II + DD*s^2],[s]);
    H = R*G;
    Q = H/(1+H);
    y=csim(w,t,Q);
    v = w - y;
    summeQuadrate = v * v';   //  [1,2,3]  * [1;2;3] = 1*1 + 2*2 +3*3
    F = sqrt( summeQuadrate );
endfunction


P=[10,0.1,0.1]
ERR = berechneFehler(P)

//------ Optierer ------

alpha = 0.01 //Güte der Lösung hängt von alpha ab.
alpha_min = 0.005 //Güte der Lösung hängt von alpha ab.

RICHTUNGEN = [ 1, 0, 0;  -1, 0, 0;  0, 1, 0;  0,-1, 0;  0, 0, 1; 0, 0, -1];
i=1

letzteN = 30;  //Streuung der letzten N Werte benutzen
FEHLER = grand(1, letzteN, "unf", 0, 100.0)

while %T
   
    F = berechneFehler(P);
    index = 1 + floor(grand(1, 1, "unf", 0, 5.999999999));  //liefert 1..4
    disp("Durchlauf:");
    disp(i);
    disp("alpha:");
    disp(alpha);
    aktuelleRichtung = RICHTUNGEN(index,:);

    Pneu = P + alpha*aktuelleRichtung;
    Fneu = berechneFehler(Pneu);

    if (Fneu<=F) then
        P = Pneu;
    end

    if (Fneu>F) then
        alpha = alpha*0.999;
    else
        alpha = alpha*1.001;
    end    

//++++++
    FEHLER(1+modulo(i,letzteN))=F;

    STREUUNG = stdev(FEHLER);
    disp("STREUUNG=");
    disp(STREUUNG);
    disp("PID=");
    disp(P);
    
    if STREUUNG<0.00000001 then
        break;
    end

    i=i+1;
end

Code 0-17: Musterlösung

#3 Mittwoch 15.12.2021

Überarbeitete Lösung zu einem Optimierer eines PID-Reglers

Die nachfolgende überarbeitete Variante des Optimierers läuft 1000 Schritte fest eingestellt durch.
Die zuvor als Abbruchkriterium eingestellte Verringerung der Streuung der auftretenden verbesserten letzten N Parameter müßte um zu einer solchen hohen Anzahl an Achritten zu gelangen noch wesentlich verkleinert werden.
Mit der besten Parameterkombination wird am Ende eine letzte Simulation durchgeführt und in rot über den Verlauf mit den Startparamatern (blau) gelegt:

Plot der beiden simulierten regelungstechnischen Systeme. blau: Startparameter für PID, rot: Simulation mit den Parametern nach der Optimierung.

Bild 0-1: Plot der beiden simulierten regelungstechnischen Systeme. blau: Startparameter für PID, rot: Simulation mit den Parametern nach der Optimierung.

Die Verbesserung ist deutlich zu sehen.
Das System ist wesentlich schneller auf den neuen Sollwert eins eingeregelt.
Startparameter: PP = 10.0; II = 0.1; DD = 0.5;
Parameter nach der Optimierung: Pbest=10.989295; Ibest=1.1802558; Dbest=1.6549533;

clear;

//Regelstrecke

//  m x.. = -Cx + F
//    x.. =  -x + F
//    xs^2+x = F

//               1
//  G(s) =  ---------- 
//           s^2 +  1

//Regler

//  klassisch: R(s) = K*(1+1/(Tn*s) + Tv*s)
//  für Optimierung verworfen, da K auf alle drei Terme wirkt!, deshalb:

//  R(s) = PP + II/s  + DD*s

//          PP*s + II + DD*s^2
//  R(s) = --------------------
//                s

PP = 10.0;
II = 0.1;
DD = 0.5;

s = poly(0,"s"); 
//Übertragungsfunktion der Regelstrecke definieren:
G = syslin('c',[1],[s^2 +  1]);
//Übertragungsfunktion des Reglers definieren:
R = syslin('c',[PP*s + II + DD*s^2],[s]);

H = R*G;
Q = H/(1+H);

t=0:0.01:15;
zs = size(t);
w = ones(1,zs(2));

y=csim(w,t,Q);
plot(t,y,'blu--');

Pbest=0;
Ibest=0;
Dbest=0;

function F = berechneFehler(PARAM)
    PP = PARAM(1);
    II = PARAM(2);
    DD = PARAM(3);
    R = syslin('c',[PP*s + II + DD*s^2],[s]);
    H = R*G;
    Q = H/(1+H);
    y=csim(w,t,Q);
    v = w - y;
    summeQuadrate = v * v';   //  [1,2,3]  * [1;2;3] = 1*1 + 2*2 +3*3
    F = sqrt( summeQuadrate );
endfunction


P=[10,0.1,0.1]
ERR = berechneFehler(P)

//------ Optierer ------

alpha = 0.01 //Güte der Lösung hängt von alpha ab.
alpha_min = 0.005 //Güte der Lösung hängt von alpha ab.

RICHTUNGEN = [ 1, 0, 0;  -1, 0, 0;  0, 1, 0;  0,-1, 0;  0, 0, 1; 0, 0, -1];
i=1

letzteN = 30;  //Streuung der letzten N Werte benutzen
FEHLER = grand(1, letzteN, "unf", 0, 100.0)

for i=1:1000
//while %T
   
    F = berechneFehler(P);
    index = 1 + floor(grand(1, 1, "unf", 0, 5.999999999));  //liefert 1..4
    aktuelleRichtung = RICHTUNGEN(index,:);

    Pneu = P + alpha*aktuelleRichtung;
    Fneu = berechneFehler(Pneu);

    if (Fneu<=F) then
        P = Pneu;
        Pbest = P(1);
        Ibest = P(2);
        Dbest = P(3);
    end

    if (Fneu<F) then
        disp("Durchlauf:");
        disp(i);
        disp("alpha:");
        disp(alpha);       
        disp("Fneu:");
        disp(Fneu);       
    end

    if (Fneu>F) then
        alpha = alpha*0.999;
    else
        alpha = alpha*1.001;
    end    

//++++++
    FEHLER(1+modulo(i,letzteN))=F;

    STREUUNG = stdev(FEHLER);
    disp("STREUUNG=");
    disp(STREUUNG);
    disp("PID=");
    disp(P);
    
    if STREUUNG<0.000000001 then
        break;
    end

    //i=i+1;
end

R = syslin('c',[Pbest*s + Ibest + Dbest*s^2],[s]);
H = R*G;
Q = H/(1+H);
y=csim(w,t,Q);
plot(t,y,'red');

Code 0-18: Überarbeitete Lösung zu einem Optimierer eines PID-Reglers

Übung bis kommende Woche

Statt einer über mehr als eine Stunde währenden Übung mit anschließender Besprechung innerhalb eines acht stündigen Kurses, soll ab nun eine aufwändigere Übungsaufgabe bis zur kommenden Woche gelöst werden und wird dann zu Beginn besprochen. Dabei wird Stoff zugrunde gelegt, der bereits in kleineren Vorübungen vertieft behandelt wurde. Allerdings wird aber oft noch Bezug auf Unterrichtsstoff genommen, der schon weiter zurück liegt.

die Aufgabe:

Obiges System aus Regelstrecke, Regler und Optimierer nach dem Prinzip des modifizierten Gradientenverfahrens soll auf die Regelung eines Einachsers übertragen werden, wie er im letzten Semester hergeleitet wurde. Nehmen Sie sich dazu zunächst folgende Inhalte aus dem letzten Semester vor:

62_Regelungssysteme/98_day_by_day_WS2021_SoSe21 -- Schauen Sie ab: "Musterlösung zu Übung 2 vom 15.12.2020"

Es wurden hier schon recht gute Lösungen durch Vorgabe der Polstellen gefunden.
Teilweise wurde nur der Kippwinkel φ ausgeregelt, teilweise auch die Position des Fahrzeugs x.
Setzen Sie einen Optimierer von φ für das nichtlinearisierte System um.
Versuchen Sie danach auch den Verlauf von φ und x zu optimieren.
Bleiben Sie ruhig bei der Optimierung eines Zustandsreglers.
Erweitern Sie am Ende den zu optimierenden Regler um einen I-Anteil.

Beachten Sie, dass hier die Darstellung der Regelstrecke im Zeitbereich erfolgt. Dementsprechend erfolgt auch die Formulierung des Reglers im Zeitbereich.

Teil 2: Fuzzy-Regelung

Als neues Thema wird im folgenden die Theorie und die Verwendung von Fuzzy-Reglern eingeführt. Eine ausführliche Übung hierzu wird es dann erst nach den Weihnachtsferien geben.

62_Regelungssysteme/18_Fuzzy -- Motivation durch ein konkretes Beispiel

62_Regelungssysteme/18_Fuzzy/01_Fuzzylogik -- Was ist Fuzzy-Logik?

62_Regelungssysteme/18_Fuzzy/02_FuzzyRegler -- Wie ist ein Fuzzy-Regler aufgebaut?

62_Regelungssysteme/18_Fuzzy/03_Uebung9 -- Kleines Beispiel mit Lösung

ÜBUNG: Fuzzy-System mit Scilab formulieren

Fuzzy-Regler gemäß "Übung 2016" hier:

62_Regelungssysteme/18_Fuzzy/03_Uebung9

1. Funktionen zur Abbildung der Typen von Fuzzygrößen formulieren

function Z = berechneFuzzygroesseDreieck(linksNull, rechtsNull, sensorwert)
    mitte = 0.5*(linksNull + rechtsNull);
    steigung = 1 / (mitte - linksNull);
    if sensorwert < linksNull then
        Z = 0;
    elseif sensorwert > rechtsNull then
        Z = 0;
    elseif sensorwert < mitte then
        Z = steigung*(sensorwert - linksNull);
    else
        Z = 1 - steigung*(sensorwert - mitte);
    end
endfunction

function Z = berechneFuzzygroesseLinksOffen(linksEins, rechtsNull, sensorwert)
    steigung = -1 / (rechtsNull - linksEins);
    if sensorwert < linksEins then
        Z = 1;
    elseif sensorwert > rechtsNull then
        Z = 0;
    else
        Z = 1 + steigung * (sensorwert - linksEins);
    end
endfunction

function Z = berechneFuzzygroesseRechtsOffen(linksNull, rechtsEins, sensorwert)
    steigung = 1 / (rechtsEins - linksNull);
    if sensorwert < linksNull then
        Z = 0;
    elseif sensorwert > rechtsEins then
        Z = 1;
    else
        Z = steigung * (sensorwert - linksNull);
    end
endfunction

function A = berechneTrapezflaeche(linksNull, rechtsNull, hoehe)
    breiteUnten = rechtsNull - linksNull;
    breiteOben = breiteUnten - hoehe*breiteUnten;
    A = hoehe * 0.5 * ( breiteUnten + breiteOben );
endfunction

//Implementierung des Fuzzy-Reglers als Scilab-Funktion:
function U = berechneSpannung(T,V)
    // 1) Bestimmung der Zugehörigkeitswerte der Fuzzygroessen der Eingangsfuzzysets
    TsehrNiedrig = berechneFuzzygroesseLinksOffen(-20, -10, T);
    Tniedrig     = berechneFuzzygroesseRechtsOffen(-20, -10, T);
    Vschwach     = berechneFuzzygroesseLinksOffen(20, 40, V);
    Vstark       = berechneFuzzygroesseRechtsOffen(20, 40, V);
    // 2) Aktivierung der einzelnen Regeln
    R1 = min(TsehrNiedrig,Vstark);
    R2 = min(TsehrNiedrig,Vschwach);
    R3 = min(Tniedrig,Vstark);
    R4 = min(Tniedrig,Vschwach);
    // 3) Schwerpunktbildung im Ausgangs-Fuzzyset
    A1 = berechneTrapezflaeche(10, 30, R1);
    A2 = berechneTrapezflaeche(0, 20, R2);
    A3 = berechneTrapezflaeche(10, 30, R3);
    A4 = berechneTrapezflaeche(0, 20, R4);

    //Schwerpunkte der Trapeze, die dnn gewichtet werden:
    s1 = 20;
    s2 = 10;
    s3 = 20;
    s4 = 10;

    U = (s1*A1+s2*A2+s3*A3+s4*A4)/(A1+A2+A3+A4);
endfunction

Code 0-19: Funktionen zur Abbildung der Typen von Fuzzygrößen formulieren

#5 Mittwoch 05.01.2022

Die Prüfung wird als E-Test gemeinsam mit Informatik-Ergänzungen am Donnerstag 20.01. ab 8:30 in IWZ140 durchgeführt. Siehe auch: Infos auf moodle.

Themen in Simulations- und Regelungstechnik 2 werden sein (alles auf der Basis von Scilab):

Kleinste Quadrate Methode
Modifiziertes Gradientenverfahren
Fuzzy-Regler
(Evolutionäre Algorithmen) -- wird ev. Informatik-Erg. zugeschlagen, falls dort nicht noch Deep Learning behandelt wird.

Im Vergleich die Themen in Informatik-Ergänzungen:

Java (OOP Modifikatoren, Schnittstellen, Vererbung, Kapselung, Threads)
Neuronale Netze und Backpropagation
Genetische Algorithmen

Themen heute:

Besprechung der Hausaufgabe (Modifiziertes Gradientenverfahren zur Optimierung eines invertierenden Pendels)
Wiederhol-Übung
Entwurf eines Fuzzy-Reglers für das invertierende Pendel
Optimierung des Fuzzy-Reglers mittels eines modifizierten Gradientenverfahrens

Thema 1: Besprechung der Hausaufgabe (Modifiziertes Gradientenverfahren zur Optimierung eines invertierenden Pendels)

clear();
clf;
m = 1; //kg
g = 0.91; // 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

A = [ 0 , 1 ; m*g*h/J , 0 ];
B = [ 0 ; h/J ];
P=[1.5, 1.3];
Palt=P;
function f = rechteSeite(t,y)
    //PAR=xxx; //Probelm: er kennt das PAR nicht, deswegen laueft hier gar nichts
    phi = y(1,1);
    om = y(2,1);
    f(1,1) = om;
    FA = -P * [phi ; 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;
endfunction

t = linspace(0,20,3000);
y0 = [10.0*%pi/180,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,:)');

w=0;

function F = berechneFehler()
    rechteSeite(t,y);
    y = ode(y0,t0,t,rechteSeite);
    diff = w(1,:) - y(1,:);
    summeQuadrate = diff * diff';
    F = sqrt( summeQuadrate );
endfunction


ERR = berechneFehler();


alpha = 0.05 //Güte der Lösung hängt von alpha ab.
alpha_min = 0.01 //Güte der Lösung hängt von alpha ab.

RICHTUNGEN = [ 1, 0; -1, 0; 0, 1; 0,-1];
i=1;

letzteN = 30; //Streuung der letzten N Werte benutzen
FEHLER = grand(1, letzteN, "unf", 0, 10.0)

PP=P;

for i=1:20//0
    //while %T
    P=PP;
    F = berechneFehler();
        index = 1 + floor(grand(1, 1, "unf", 0, 3.999999999)); //liefert 1..4
        aktuelleRichtung = RICHTUNGEN(index,:);
    
    Pneu = PP + alpha*aktuelleRichtung;
    P=Pneu;
    Fneu = berechneFehler();
    //liefer komplexe werte --> real part nutzen
    if (real(Fneu)<=real(F)) then
        PP = Pneu;
    end
    
    if (real(Fneu)<=real(F)) then
        disp("Durchlauf:");
        disp(i);
        disp("alpha:");
        disp(alpha);
        disp("Fneu:");
        disp(Fneu);
    end
    
    if (real(Fneu)<=real(F)) then
        alpha = alpha*0.999;
      else
        alpha = alpha*1.001;
    end
    /*
    //++++++
    FEHLER(1+modulo(i,letzteN))=real(F);
    
    STREUUNG = stdev(FEHLER);
    disp("STREUUNG=");
    disp(STREUUNG);
    disp("PID=");
    disp(P);
    
    if STREUUNG<0.000000001 then
    break;
    end
    
    //i=i+1;
    */
end

Palt
P
y = ode(y0,t0,t,rechteSeite);
plot(t,y,'red');

Code 0-20: Studentische Lösung, Aufgabe siehe:"Übung bis kommende Woche" bei #3 Mittwoch 15.12.2021.

Thema 2: Wiederhol-Übung

Nehmen Sie sich das Dokument RT2_OHNE.pdf auf Moodle vor und lösen die Aufgaben zur Methode der kleinsten Quadrate und zu Fuzzy.
Überprüfen Sie Ihre Lösung selber mittels RT2_MIT_LOESUNG.pdf und fragen im Zweifelsfall im Plenum nach.

Nachfolgende Aufgabe noch nicht bearbeiten, erst einmal nur die beiden E-Test-Aufgaben.

VERÄNDERT!: Nehmen Sie sich die studentische Lösung aus Teil 1 vor:

Sorgen Sie dafür, dass die Stellkraft extrahiert und geplottet werden kann.
Finden Sie eine günstige Lösung mittels Polvorgabe, wenn nur ein Zustandsregler (PD-Regler) realisiert wird.
Begrenzen Sie die Stellkraft auf einen Maximalwert, der 2/3 des maximalen Wertes darstellt, der bei der Lösung mit Polvorgabe entsteht.
Führen Sie eine automatische Optimierung (für eine PD-Regler) unter diesen Bedingungen durch und vergleichen Sie die Güte der anfänglichen Lösung aus der Polvorgabe mit der aus der Optimierung. Statt des Abbruchkriteriums, optimieren Sie einfach eine begrenzte Anzahl von Schritten, z.B. 100.

Hinweis: Der ODE-Integrator könnte versagen, wenn F begrenzt wird, da dann ein "Knick" (nicht differenzierbar an der Stelle) im Verlauf auftritt, darum sollten Sie den Runge-Kutta-Integrator aus folgendem Beispiel benutzen:

clear();

xtot=0;

function f=modell(y,t,dt)
x=y(1,1);
v=y(2,1);

Kkrit = 19.75;
xsoll = 1.0;
//e = xsoll - x;
e = xsoll - xtot;
u=Kkrit*e;

f(1,1)=v;
f(2,1)=-x-2*v+u;
endfunction

function yneu=ruku(y,t,dt)
k1=modell(y,t,dt);
k2=modell(y+0.5.*dt.*k1,t+0.5*dt);
k3=modell(y+0.5.*dt.*k2,t+0.5*dt);
k4=modell(y+dt.*k3,t+dt);
yneu=y+dt.*(k1+2.*k2+2.*k3+k4)./6;
endfunction

tmax = 20.0;
dt = 0.01;
schritte = ceil(tmax/dt);

yalt = [0,0]';
ysim = yalt;
t=0.0;
tt=t;

Ttot = 0.1;
anztot = round(Ttot/dt)
xtotarr = zeros(anztot,1);

for i=1:1:schritte
yneu=ruku(yalt,t,dt);
yalt=yneu;
ysim=[ysim,yalt];
tt =[tt,t];
t=t+dt;

xtot = xtotarr(modulo((i-1),anztot)+1);
xtotarr(modulo((i-1),anztot)+1)=yneu(1);
end

plot(tt,ysim(1,:));

Code 0-21: Beispiel mit selbst implementiertem Runge-Kutta-Verfahren (Einschritt-Verfahren 3. Ordnung)

begrenzteStellkraft.zip -- Musterlösung (in mehreren Schritten, Variante 004 ist unbrauchbar, 004b ist aktuell)

Nachtrag:

clear();
clf;
m = 1; //kg
g = 0.91; // 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

A = [ 0 , 1 ; m*g*h/J , 0 ];
B = [ 0 ; h/J ];

P=[1.5, 1.3];
Palt=P;

FMAX=0.167;

//MODELL, invertierende Pendel
function f = modell(y,t)
    phi = y(1,1);
    om = y(2,1);
    f(1,1) = om;
    FA = -P * [phi ; om];
    if FA>FMAX then 
        FA=FMAX;
    elseif FA<-FMAX then
        FA=-FMAX;
    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

function yneu=ruku(y,t,dt)
   k1=modell(y,t);
   k2=modell(y+0.5.*dt.*k1,t+0.5*dt);
   k3=modell(y+0.5.*dt.*k2,t+0.5*dt);
   k4=modell(y+dt.*k3,t+dt);
   yneu=y+dt.*(k1+2.*k2+2.*k3+k4)./6;
endfunction

t = linspace(0,20,3000);
y0 = [10.0*%pi/180,0.0]'; //Start in der Nähe der instabilen Ruhelage
t0 = 0;

//Simulation mit den Startparametern:
//y = ode(y0,t0,t,rechteSeite);
Fv = 0;
Faktuell = 0;
yalt=y0;
t=0;
dt=0.01;
schritte=2000;
tt=0;
ysim=yalt;
for i=1:1:schritte
   yneu=ruku(yalt,t,dt);
   yalt=yneu;
   ysim=[ysim,yalt];
   tt =[tt,t];
   t=t+dt;
end

Fv = 0;
for i=1:1:schritte
   FA = -P(1) * ysim(1,i) - P(2) * ysim(2,i);
    if FA>FMAX then 
        FA=FMAX;
    elseif FA<-FMAX then
        FA=-FMAX;
    end
   Fv = [Fv,FA];
end

plot(tt,ysim(1,:)','red');

//Fverlauf = -P(1) * y(1,:) - P(2) * y(2,:);
//plot(t,Fverlauf,'red');
//plot(tt,Fv,'red');



function F = berechneFehler()
    //y = ode(y0,t0,t,rechteSeite);
    yalt=y0;
    t=0;
    dt=0.01;
    schritte=2000;
    tt=0;
    ysim=yalt;
    for i=1:1:schritte
        yneu=ruku(yalt,t,dt);
        yalt=yneu;
        ysim=[ysim,yalt];
        tt =[tt,t];
        t=t+dt;
    end
    summeQuadrate = ysim(1,:) * ysim(1,:)';
    F = sqrt( summeQuadrate );
endfunction


ERRvorher = berechneFehler()


alpha = 0.05 //Güte der Lösung hängt von alpha ab.
alpha_min = 0.01 //Güte der Lösung hängt von alpha ab.

RICHTUNGEN = [ 1, 0; -1, 0; 0, 1; 0,-1];
i=1;

PP=P;

for i=1:100
    P=PP;
    F = berechneFehler();
    index = floor(grand(1, 1, "unf", 1, 4.999999999)); //liefert 1..4
    aktuelleRichtung = RICHTUNGEN(index,:);
    
    Pneu = PP + alpha*aktuelleRichtung;
    P=Pneu;
    Fneu = berechneFehler();
    //liefer komplexe werte --> real part nutzen
    if (real(Fneu)<=real(F)) then
        PP = Pneu;
    end
    
    if (real(Fneu)<=real(F)) then
        disp("Durchlauf:");
        disp(i);
        disp("alpha:");
        disp(alpha);
        disp("Fneu:");
        disp(Fneu);
    end
    
    if (real(Fneu)<=real(F)) then
        alpha = alpha*0.999;
      else
        alpha = alpha*1.001;
    end
end

Palt
P
//y = ode(y0,t0,t,rechteSeite);
yalt=y0;
t=0;
dt=0.01;
schritte=2000;
tt=0;
ysim=yalt;
for i=1:1:schritte
   yneu=ruku(yalt,t,dt);
   yalt=yneu;
   ysim=[ysim,yalt];
   tt =[tt,t];
   t=t+dt;
end

//plot(t,y(1,:),'gre');
plot(tt,ysim(1,:)','gre');
//FverlaufNEU = -P(1) * y(1,:) - P(2) * y(2,:);
Fv = 0;
for i=1:1:schritte
   FA = -P(1) * ysim(1,i) - P(2) * ysim(2,i);
    if FA>FMAX then 
        FA=FMAX;
    elseif FA<-FMAX then
        FA=-FMAX;
    end
   Fv = [Fv,FA];
end


//plot(tt,Fv,'gre');


ERRvorher
ERR = berechneFehler()

Code 0-22: FMAX=0.167N führt zu einem Fehler, der am Anfang bei 2.27 liegt und am Ende bei 1.92.

Bild 0-2: Auslenkung phi, rot: vor Optimierung, grün: nach Optimierung.

Bild 0-3: Kraft FA, rot: vor Optimierung, grün: nach Optimierung.

Thema 3: Entwurf eines Fuzzy-Reglers für das invertierende Pendel

03_Archiv/02_WS2020_21/01_RTS/01_day_by_day --- "Dienstag 05.01.2021, Formel 0-8, 0-9, 0-10" dienen als hier zu verwendendes Modell. Es wird nur phi geregelt. Parameter, siehe Code 0-9 ebenda.

Es macht Sinn, die entwickelten Scilab-Funktionen in verschiedenen Files zu halten und diese dann zur Bereitstellung einer Funktionsgruppe jeweils einmal auszuführen.

Teilaufgaben:

Machen Sie einen Vorschlag für ein Fuzzy-System (Eingangs-Fuzzy-Set / Ausgangs-Fuzzy-Set / Regeln).
Formulieren Sie Hilfsfunktionen, die die Basis für alle Berechnungen bilden, insbesondere für eine nicht äquidistante Anordnung aller Stützstellen aller Fuzzy-Sets.
Vervollständigen und testen Sie das Gesamtsystem.

Thema 4: Optimierung des Fuzzy-Reglers mittels eines modifizierten Gradientenverfahrens

Jede Stützstelle eines Fuzzy-Sets ist ein Parameter für den Optimierer. Überlegen Sie, wie Verschiebungen dieser Stützstellen vom Optimierer vorgenommen werden können, ohne dass sich deren Reihenfolge ändern kann.
Implementieren und testen Sie den Optimierer.

#6 Mittwoch 12.01.2022

Themen:

Beantwortung von Fragen zu den Übungsaufgaben und zu der geplanten Klausur.
Umsetzung und Optimierung eines Fuzzy-Reglers für das invertierende Pendel, je nach Verfügbarkeit als Besprechung der "Hausaufgabe", oder als Übung.

Konkrete Umsetzung des Fuzzy-Reglers und dessen Optimierung Mit Scilab

Wie sieht noch gleich ein Fuzzy-Regler aus? -- Siehe noch einmal folgendes Beispiel:

62_Regelungssysteme/18_Fuzzy/03_Uebung9

Wichtige "Erkenntnis" vorab: globale Variablen lassen sich in Scilab-Funktionen lesen, aber nicht ändern.

Versuch dazu:

disp("Experiment zur Sichtbarkeit und Veränderbarkeit globaler Variablen in Scilab");

global = 7;

function ergebnis = meineFunktion()
    disp("In der Funktion:");
    disp(global);
    global = 8;
    disp("In der Funktion nach Änderung von global:");
    disp(global);

    ergebnis = global;
endfunction

disp("Außerhalb der Funktion VOR dessen Aufruf:");
disp(global);

x = meineFunktion()

disp("Außerhalb der Funktion NACH dessen Aufruf:");
disp(global);

Code 0-23: Experiment zu Funktionen und globalen Variablen.

 "Außerhalb der Funktion VOR dessen Aufruf:"
   7.


  "In der Funktion:"

   7.

  "In der Funktion nach Änderung von global:"

   8.
 x  = 

   8.


  "Außerhalb der Funktion NACH dessen Aufruf:"
   7.

Code 0-24: Ausgabe des Experiments.

Damit ein Optimierer das Fuzzy-System verändern kann, werden Stützstellen der Fuzzy-Sets global definiert, im Fuzzy-Regler benutzt und im Optimierer verändert.
Der Optimierer darf dann keine Funktion sein.
Jedes der drei Fuzzy-Sets soll aus fünf Fuzzy-Größen bestehen.
Alle Fuzzy-Größen werden mit Indizes angesprochen.
Die Regeln werden durch eine Matrix beschrieben: Zeile==Auslenkungs-Fuzzy-Größe PHI, Spalte==Winkelgeschwindigkeits-Fuzzy-Größe OMEGA, Matrixeintrag: Index der Ausgangs-Fuzzy-Größe.
Jedes Eingangs-Fuzzyset kann durch 5 Stützstellen auf der Ordinate ("X-Achse") eindeutig definiert werden.
Jedes Ausgangs-Fuzzyset kann durch 7 Stützstellen auf der Ordinate ("X-Achse") eindeutig definiert werden.
Optimieren heißt hier Stützstellen verschieben und behalten, wenn sich das Verhalten bessert, oder die Änderung andernfalls zu verwerfen.
Der Fuzzy-Regler fuzzyRegler(PHI,OMEGA) hat zwei Input-Werte: PHI und OMEGA und einen Ausgangswert: Die Antriebskraft FA.
Ansonsten soll eine Hierarchische Struktur von Hilfsfunktionen entwickelt werden, die die Arbeit von fuzzyRegler(...) unterstützen.
fuzzyRegler(PHI,OMEGA) ersetzt den Zustandsregler beim invertierenden Pendel.

Nachfolgende Schaubilder illustrieren die aufgelisteten Punkte:

Bild 0-4: Aufbau des umzusetzenden Fuzzy-Systems

Bild 0-5: Programmstruktur.

Teilimplementierung der geforderten Funktionen:


function x = fuzzyGroesseDreieck(WERT,a,b,c)
    if ((WERT > a) & (WERT <= b)) then
        m = 1/(b-a);
        x = m*(WERT-a);
    elseif ((WERT > b) & (WERT < c)) then
        m = 1/(c-b);
        x = 1 - m*(WERT-b);
    else
        x = 0;
    end  
endfunction


function s = schwerpunktDreieckLinks(a,b)
    s = a + (2/3)*(b-a);
endfunction

function b = bestimmeB(hoehe,a,q,d)
    m = 1/(q-a);
    // hoehe = m*(b-a)
    b = (hoehe/m) + a;
endfunction

Code 0-25: Teilimplementierung der geforderten Funktionen.

Übung

Entwickeln und optimieren Sie den Fuzzy-Regler für das invertierende Pendel

Vorgehen:

Legen Sie sinnvolle Ausgangswerte für die Stützpunkte der Fuzzy-Sets und für die Regelmatrix fest.
Implementieren Sie alle geplanten Scilab-Funktionen, die das Fuzzy-Regler repräsentieren.
Simulieren Sie das invertierende Pendel unter Verwendung des Fuzzy-Reglers.
Implementieren Sie den Optimierer für den Fuzzy-Regler.
Optimieren Sie den Fuzzy-Regler.

StudentischeLoesung_Pendel_FuzzyRegler.zip -- studentische Teillösung (Stützstellen noch nicht parametrisiert / ohne Optimierung)

Bild 0-6: Schaubild zur studentischen Teillösung

musterloesung_fuzzyRegler_Optimierung_invPendel.zip -- Musterlösung.

Day by Day

(EN google-translate)

(PL google-translate)

#1 Mittwoch 24.11.2021

Übersicht und Vereinbarungen zur Ausgestaltung des Unterrichts

Themen

Statische Optimierung

Allgemeines zu Optimierungsverfahren

Übung

studentische Lösung

Lösungsvariante

Übung 2

Quelltext zu LSQ004_LSQ3x3

#2 Mittwoch 01.12.2021

Thema: Fortsetzung mit Optimierung, jetzt evolutionäre Algorithmen

Heuristiken als Grundlage für die Implementierung von Optimierungsverfahren

Prinzip des Gradientenverfahrens

Genetische Algorithmen

Übung

Vorlage für die prozedurale Umsetzung:

#3 Mittwoch 08.12.2021

Heuristiken als Grundlage für die Implementierung von Optimierungsverfahren

Prinzip des modifizierten Gradientenverfahrens

Übung

Während des Unterrichts entstandene Lösungen

Erschwerte aber realistische Situation: Kein Gradient von F(x,y) verfügbar!!!

Studentische Lösung Schrittweitensteuerung UND Abbruchkriterium

Neuer Ausgangspunkt mit zufälliger Richtungswahl

Musterlösung zur Übung

#3 Mittwoch 15.12.2021

Überarbeitete Lösung zu einem Optimierer eines PID-Reglers

Übung bis kommende Woche

Teil 2: Fuzzy-Regelung

ÜBUNG: Fuzzy-System mit Scilab formulieren

1. Funktionen zur Abbildung der Typen von Fuzzygrößen formulieren

#5 Mittwoch 05.01.2022

Thema 1: Besprechung der Hausaufgabe (Modifiziertes Gradientenverfahren zur Optimierung eines invertierenden Pendels)

Thema 2: Wiederhol-Übung

Thema 3: Entwurf eines Fuzzy-Reglers für das invertierende Pendel

Thema 4: Optimierung des Fuzzy-Reglers mittels eines modifizierten Gradientenverfahrens

#6 Mittwoch 12.01.2022

Konkrete Umsetzung des Fuzzy-Reglers und dessen Optimierung Mit Scilab

Übung

Entwickeln und optimieren Sie den Fuzzy-Regler für das invertierende Pendel