Day by Day
(EN google-translate)
(PL google-translate)
Verzeichnis der im Verlauf des Semesters behandelten Themen

Dienstag 24.11.2020
|
Einführung in Scilab


Einführung in Regelungstechnik



//s als Argument für ein Polynom definieren: s = poly(0,"s"); //Übertragungsfunktion definieren: G = syslin('c',[1],[1+2*s+s^2]); t=[0:0.01:10]; u=ones(1,1001); y=csim(u,t,G); plot2d(t,y); //plzr(G);
Code 0-1: Übertragungsfunktionen beschreiben und auswerten mit Scilab.



Übung 1
|
$ \left[\begin{array}{cc}3 & 2 \\ 5 & 1\end{array}\right] \cdot \left[\begin{array}{cc}x \\ v\end{array}\right] = \left[\begin{array}{cc}-1 \\ 3\end{array}\right] $
Formel 0-1: Aufgabe 2: Lösen Sie obiges Gleichungssystem erst von hand und dann mit Hilfe von Scilab.
Aufgabe 3: Ein linearer Feder-Masse-Schwinger habe die Masse 1kg und die Federkonstante C = 1 N/m. Wie muß die Dämpfung gewählt werden, damit die Eingefrequenz 1Hz beträgt?
Beachte: Die Eigen-Kreisfrequenz(en) entspricht/entsprechen dem imaginären Anteil der Eigenwerte in einem linearen Differentialgleichungssystem (LDGLS).
clear; function f = rechteSeite(t,y) x = y(1,1); v = y(2,1); f(1,1) = v; f(2,1) = -x-1.56*v; endfunction t = 0:0.1:100; y0 = [1,0]'; t0 = 0; y = ode(y0,t0,t,rechteSeite); plot(t,y(1,:)',t,y(2,:)');
Code 0-2: Scilab-Code zur Simulation des Systems von Übung1, Aufgabe 3.
clear; s = poly(0,"s"); //Übertragungsfunktion definieren: G = syslin('c',[1],[1+1.56*s+s^2]); t=[0:0.01:100]; u=ones(1,10001); y=csim(u,t,G); plot2d(t,y);
Code 0-3: Scilab-Code zur Simulation des Systems von Übung1, Aufgabe 3 im Laplacebereich.
Übung 2
|

|

Aufgabe 3: Ein linearer Schwinger (m=1kg, C=1N/m) ohne Dämpfungselement wird mit einem aktiven Dämpfer ausgerüstet, also einem D-Element. Die Regeldifferenz zum Sollwert Null m/s wird mit Faktor 10 als äußere Kraft auf das System gegeben.
|

Während der Vorlesung entstandene Videos:


Themen im letzten Teil:
|




Dienstag 01.12.2020
Themen
TEIL 1
|
TEIL 2
|
Zu Übung 2, Aufgabe 3 von letzter Woche:
Um eine exakte Entsprechung zwischen der Simulation der Variante im Zeitbereich und der im Laplacebereich zu erzielen, müssen w und dw/dt in beiden Varianten gleich sein. Bei Vorgabe eines Sprungs ist das nicht gewährleistet, da dw/dt nicht eindeutig bestimmt ist. Darum wird die Aufgabe etwas angepaßt und w=sin(t)*exp(-t), bzw. dw/dt=cos(t)*exp(-t)-sin(t)*exp(-t) gesetzt:
$ \ddot x + x = u $
Formel 0-2: System im Zeitbereich.
$ u = P \cdot \left(w-x\right) + D \cdot \left(\dot w - \dot x\right) $
Formel 0-3: Regler im Zeitbereich.
$ G\left(s\right)= \frac {1}{s^2+1} $
Formel 0-4: System im Laplace-Bereich \left(y/u\right).
$ R\left(s\right)= \frac {P+s \cdot D}{1} $
Formel 0-5: Regler im Laplace-Bereich \left(u/e\right).
$ Q\left(s\right)= \frac {P+s \cdot D}{s^2+1} $
Formel 0-6: Offener Regelkreis im Laplace-Bereich \left(y/e\right).
$ H\left(s\right)= \frac {P+s \cdot D}{s^2+s \cdot D+P+1} $
Formel 0-7: \left(Geschlossener\right) Regelkreis im Laplace-Bereich \left(y/w\right).
clear; P = 10.0; D = 1.0; s = poly(0,"s"); //Übertragungsfunktion definieren: G = syslin('c',[1],[1+s^2]); // Regelstrecke R = syslin('c',[P+s*D],[1]); // PD-Regler Q = R*G; //Übertragungsv. offener Rk. H = Q/(1+Q); //Übertrv. geschl. Rk. t=[0:0.01:15]; u=sin(t).*exp(-t); y=csim(u,t,H); plot(t,y,"r"); //plzr(G);
Code 0-4: Simulation des PD-geregelten linearen Schwingers bei Darstellung des Modells im Laplacebereich.
clear; P = 10.0; D = 1.0; function f = rechteSeite(t,y) x = y(1,1); v = y(2,1); w = sin(t)*exp(-t); derw = cos(t)*exp(-t) - sin(t)*exp(-t); u = P*(w-x) + D*(derw-v); f(1,1) = v; f(2,1) = -x+u; endfunction t = 0:0.01:15; y0 = [0,0]'; t0 = 0; y = ode(y0,t0,t,rechteSeite); plot(t,y(1,:),"b--");
Code 0-5: Simulation des PD-geregelten linearen Schwingers bei Darstellung des Modells im Zeitbereich.
Nachtrag: Wie man einen I-Anteil ergänzen kann:
clear; P = 10.0; D = 1.0; I = 2.0; function f = rechteSeite(t,y) x = y(1,1); v = y(2,1); e_integral = y(3,1); w = sin(t)*exp(-t); derw = cos(t)*exp(-t) - sin(t)*exp(-t); u = P*(w-x) + D*(derw-v) + I*(e_integral); f(1,1) = v; f(2,1) = -x+u; f(3,1) = w-x; endfunction t = 0:0.01:15; y0 = [0,0,0]'; t0 = 0; y = ode(y0,t0,t,rechteSeite); plot(t,y(1,:),"b--");
Code 0-6: I-Anteil als Dummy-Systemzustand ergänzt.

Bild 0-1: Aus den beiden vorangehenden Simulationen resultierende übereinanderliegende Plots.
Links zu den Inhalten der Lehrveranstaltung


Übung 1
Welche der Aufgaben (Übung 5, ganz unten auf Webseite, Link nachfolgend) läßt sich nach der 1. Methode nach Ziegler und Nichols lösen?. Setzen Sie diese Lösungen für P, PI und PID-Regler um.

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-7: Studentische Lösung für die erste Übertragungsfunktion G(s) = exp(-0.1s)/(s^2+2s+1) ergibt: K_krit=19,75; T_Krit=1,45s
|
Teil 2 entfällt. Aufgabe bis 14Uhr: Bei obigem G(s) = exp(-0.1s)/(s^2+2s+1) und bei dem zweiten System G(s) = exp(-s)/((s+1)*(s+2)) den PID-Regler nach Ziegler-Nichols auslegen und insbesondere den I-Anteil so wie weiter oben gezeigt ergänzen. Das heißt: Das System mit dem fertigen PID-Regler soll simuliert werden.

|
Übung 2


Dienstag 08.12.2020
|
Teil 1: Invertierendes Pendel


|
Übungen wurden in den Chat geschrieben.
clear; dt=0.1; t = 0:dt:100; zs = size(t); zeilen = zs(1); spalten = zs(2); // Punkt A: [x;y] zeitlicher Verlauf == [0;0] pA = zeros(2,spalten); // Punkt B: [x;y] zeitlicher Verlauf == [sin(...);1] pB = [sin(2*%pi*0.2*t);ones(1,spalten)]; //Eine Animation soll eine Verbindungslinie von pA mit pB im zeitlichen Verlauf darstellen: isoview; clf; // ABBRUCH MIT CTR C !! for i=1:spalten clf; plot([-2,-1.999,1.999,2],[0,2,2,0]); //Rahmen, um Skalierung konstant zu halten plot([pA(1,i),pB(1,i)],[pA(2,i),pB(2,i)]); sleep(1000*dt); // "Echtzeit" end
Code 0-8: Einfache Möglichkeit, mit Scilab eine Animation zu erzeugen: Mehrfach plotten mit sleep() dazwischen.
Übung 2 im Chat: Es sollte das inverse Pendel mit Scilab simuliert werden und auch eine Version umgesetzt werden, um Ziegler und Nichols anzuwenden.
Dabei hat sich gezeigt, dass nur durch Hinzufügen eines Dämpfungsterms die Methode nach Ziegler und Nichols anwendbar ist. Die während der LV entstandenen Skripte finden sich nachfolgend als zip-File:

Teil 2: Polvorgabe zur Reglerauslegung

Dienstag 15.12.2020
Themen:
|
1. Zustandsregler, 2. Polvorgabe, 4. Zustandstransformation in den Nullpunkt




3. Dämpfungsgrad / "Technisch stabil"

4. Zustandstransformation in den Nullpunkt

Übung 1: Polvorgabe für x.. = -x
Übung 1b: Polvorgabe für x.. = -x, aber Regelung zum Punkt x=1, v=0.
ÜBUNG 2: 5. Polvorgabe für das inverse Pendel
|
Videos aus der LV (Zustandsregler)





Beispiellösung / Anderer Sollwert, als Nullvektor:



Dienstag 05.01.2021
Themen
|
Musterlösung zu Übung 2 vom 15.12.2020
|
Modell des invertierenden Pendels:

Bild 0-2: Herleitung zum invertierenden Pendel, Teil 1.

Bild 0-3: Herleitung zum invertierenden Pendel, Teil 2.

$ \ddot \phi = \frac {-m \cdot h^2 \cdot \dot \phi ^2 \cdot \sin \phi \cdot \cos \phi + h \cdot m \cdot g \cdot \sin \phi + h \cdot \cos \phi \cdot F_A}{J+m \cdot h^2 \cdot \left( \sin \phi \right)^2} $
Formel 0-8: Modellteil des invertierenden Pendels für die Verkippung \phi .
$ \ddot x = \frac {F_A}{m} $
Formel 0-9: Modellteil des invertierenden Pendels für die Position x auf der Schiene.
|
$ \ddot \phi = \frac {h \cdot m \cdot g \cdot \phi + h \cdot F_A}{J} $
Formel 0-10: Erster Modellteil nach der Linearisierung.
Hier tritt deutlich die In- bzw. Grenzstabilität des Ausgangssystems zutage: kein negatives Vorzeichen vor dem Termin mit phi.
Quadratische Terme entfallen, phi statt sin(phi), 1 statt cos(phi) gemäß Taylor-Entwicklung.
$ \ddot x = \frac {F_A}{m} $
Formel 0-11: Der zweite Modellteil ist bereits linear.
|
$ \dot \vec y = A \cdot \vec y + B \cdot \vec u $
Formel 0-12: Allgemeine Form eines linearen Systems mit Zus \tan dsmatrix A und Eingriffsmatrix B.
$ \dot \vec y = A \cdot \vec y - B \cdot R \cdot \vec y $
Formel 0-13: Allgemeine Form des Zus \tan dsreglers.
u wird ersetzt durch -Ry.
$ \left[\begin{array}{cc} \dot \phi \\ \dot \omega \end{array}\right] = \left[\begin{array}{cc}0 & 1 \\ \frac {h \cdot m \cdot g}{J} & 0\end{array}\right] \cdot \left[\begin{array}{cc} \phi \\ \omega \end{array}\right] + \left[\begin{array}{cc}0 \\ \frac {h}{J}\end{array}\right] \cdot F_A $
Formel 0-14: Übertragung auf den linearisierten ersten Modellteil. Der Antriebskraft FA entspricht u.
|
clear(); m = 1; //kg g = 0.91; // 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 A = [ 0 , 1 ; m*g*h/J , 0 ]; B = [ 0 ; h/J ]; spec(A) //Eigenwerte der ursprünglichen Matrix ansehen // ergibt: // 2.3023837 + 0.i // -2.3023837 + 0.i // technisch stabil wählen, an vorhandenen Werten orientieren: EW = [-2.3023837 + %i*2.3023837 , -2.3023837 - %i*2.3023837 ]; R = ppol(A,B,EW) // ergibt (korrigiert 5.1. wg. verändertem B): // 2.7299999 0.7904851 //testweise diesen Regler in das nicht linearisierte Modell einbauen und testen: function f = rechteSeite(t,y) phi = y(1,1); om = y(2,1); f(1,1) = om; FA = -R * [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,:)');
Code 0-9: Scilab-Skript zur Bestimmung eines Zustandsreglers durch Polvorgabe und Simulation des nicht linearen Systems unter Verwendung dieses Reglers.

Bild 0-4: Ergebnis-Plot aus obigem Skript.
|
$ \left[\begin{array}{cc} \dot \phi \\ \dot \omega \\ \dot x \\ \dot v \end{array}\right] = \left[\begin{array}{cc}0 & 1 & 0 & 0 \\ \frac {h \cdot m \cdot g}{J} & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0\end{array}\right] \cdot \left[\begin{array}{cc} \phi \\ \omega \\ x \\ v\end{array}\right] + \left[\begin{array}{cc}0 \\ \frac {h}{J} \\ 0 \\ \frac {1}{m}\end{array}\right] \cdot F_A $
Formel 0-15: Linearisiertes Gesamtmodell in Matrixschreibweise.
clear(); m = 1; //kg g = 0.91; // 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 //Vorgehen für Gesamtsystem analog wie beim Teilmodell: A = [ 0,1,0,0 ; m*g*h/J,0,0,0 ; 0,0,0,1 ; 0,0,0,0 ]; B = [ 0 ; h/J ; 0 ; 1/m ]; spec(A) //Eigenwerte der ursprünglichen Matrix ansehen // ergibt: // 2.3023837 + 0.i // -2.3023837 + 0.i // 0. + 0.i // 0. + 0.i // technisch stabil wählen, an vorhandenen Werten orientieren: // Regleranteil für x langsamer machen: u=2.3023837; v=1.0; EW = [-u+%i*u, -u-%i*u, -v+v*%i, -v-v*%i ]; R = ppol(A,B,EW) // ergibt (korrigiert wg. korrigiertem B 5.1.): // 5.3409701 2.1187267 -3.9999999 -5.7373297 //testweise diesen Regler in das nicht linearisierte Modell einbauen und testen: 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 = linspace(0,20,3000); 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-10: Scilab-Skript, jetzt um das Gesamtsystem zu regeln.

Bild 0-5: Ergebnis-Plot aus obigem Skript (Gesamtmodell).
clear(); 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 R = [5.3409701, 2.1187267, -3.9999999, -5.7373297]; 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 = linspace(0,20,3000); y0 = [10.0*%pi/180,0.0,0.2,0.0]'; //Start in der Nähe der instabilen Ruhelage t0 = 0; y = ode(y0,t0,t,rechteSeite); x = y(3,:); phi = y(1,:); y = h.*cos(phi); zs = size(t); spalten = zs(2); pA = [x + h.*sin(phi);zeros(1,spalten)]; pB = [x - h.*sin(phi);y + h.*cos(phi)]; //Eine Animation soll eine Verbindungslinie von pA mit pB im zeitlichen Verlauf darstellen: isoview; clf; // ABBRUCH MIT CTR C !! for i=1:spalten clf; plot([-1.2,-1.1999,1.1999,1.2],[-1.2,1.2,1.2,-1.2]); //Rahmen, um Skalierung konstant zu halten plot([pA(1,i),pB(1,i)],[pA(2,i),pB(2,i)]); sleep(10); // "Echtzeit" end
Code 0-11: Animation der Ergenisse
ÜBUNG 5.1.: Hinzufügen eines integralen Zustandes, um einen I-Anteil im Zustandsregler zu haben:
$ \left[\begin{array}{cc} \dot \phi \\ \dot \omega \\ \dot \alpha \end{array}\right] = \left[\begin{array}{cc}0 & 1 & 0 \\ \frac {h \cdot m \cdot g}{J} & 0 & 0 \\ 1 & 0 & 0\end{array}\right] \cdot \left[\begin{array}{cc} \phi \\ \omega \\ \alpha \end{array}\right] + \left[\begin{array}{cc}0 \\ \frac {h}{J} \\ 0\end{array}\right] \cdot F_A $
Formel 0-16: Übertragung auf den linearisierten ersten Modellteil. Der Antriebskraft FA entspricht u.
...jetzt mit integrativem Zustand alpha = Integral von phi, bzw. d alpha / dt = phi.
clear(); m = 1; //kg g = 0.91; // 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 //Einführen eines fiktiven Zustandes phi_integ = Integral phi dt, //oder: d phi_integ / dt = phi (Die Ableitung dieses Zustandes ergibt phi) A = [ 0 , 1 , 0 ; m*g*h/J , 0 , 0 ; 1 , 0 , 0 ]; B = [ 0 ; h/J ; 0 ]; spec(A) //Eigenwerte der ursprünglichen Matrix ansehen // ergibt: // 0. + 0.i // 2.3023837 + 0.i // -2.3023837 + 0.i // technisch stabil wählen, an vorhandenen Werten orientieren: // Der dritte Eigenwert kann nur real sein! EW = [-2.3023837 + %i*2.3023837 , -2.3023837 - %i*2.3023837 , -2.3023837]; R = ppol(A,B,EW) // 4.5499999 1.1857276 4.1903382 //testweise diesen Regler in das nicht linearisierte Modell einbauen und testen: function f = rechteSeite(t,y) phi = y(1,1); om = y(2,1); phi_integ = y(3,1); f(1,1) = om; FA = -R * [phi;om;phi_integ]; 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) = phi; endfunction t = linspace(0,20,3000); y0 = [10.0*%pi/180,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,:)');
Code 0-12: Musterlösung für das Hinzufügen eines integralen Zustandes phi_integ, der das Integral von phi ist.

Bild 0-6: Ergebnis-Plot aus obigem Skript (Teilmodell mit zusätzlichem Zustand phi_integ).
Hinweise zur Vorbereitung und zum Ablauf der Prüfung
|
|
|
Bitte beachten Sie weitere Hinweise während der Lehrveranstaltung.

Daraus: § 19 Weitere Bildungs- sowie Aus-, Fort- und Weiterbildungseinrichtungen: (1) Präsenzangebote in Bildungs- sowie Aus-, Fort- und Weiterbildungseinrichtungen, insbesondere in Hochschulen, Musikschulen, Kunstschulen, Volkshochschulen, Fahr-, Flug- und Segelschulen sind nur mit jeweils bis zu fünf Teilnehmerinnen und Teilnehmern zulässig. Die Personengrenze nach Satz 1 gilt nicht für die Durchführung und Vorbereitung von Prüfungen sowie die Abnahme von Prüfungsleistungen.
Theoretische Behandlung des Zustandsreglers mit Beobachter

Wiederholung einzelner Themen gemäß Ihren Wünschen
**** NEW ****
