//Geschlossene Lösung:
tt=0:0.1:2;
A = exp(log(2)*tt);
plot(tt,A);

//Numerische Lösung:
A=1;
dt=0.1;

Aerg = [A];

for t=dt:dt:2
   Aneu = log(2)*A*dt + A;
   A    = Aneu; 

   Aerg = [Aerg,A]; //Sammeln der Ergebnisse
end

plot(tt,Aerg,'re-');

Code 0-1: Vergleich: geschlossene und numerische Lösung zum Seerosenteich

clear;
//Geschlossene Lösung:
dt=0.001;
tt=0:dt:2;
A = exp(log(2)*tt);
//plot(tt,A);

//Numerische Lösung:
A=1;

Aerg = [A];

for t=dt:dt:2
   Aneu = log(2)*A*dt + A;
   A    = Aneu; 

   Aerg = [Aerg,A]; //Sammeln der Ergebnisse
end

Ageschlossen = exp(log(2)*tt);

plot(tt,Aerg,'re',tt,Ageschlossen,'gr--');

//siehe help axes_properties
a = gca();
a.x_label.text = 't/Woche';
a.y_label.text = 'A/m^2';
a.title.text = 'Flächenwachstum eines Seerosenteiches';
hl=legend(['numerisch';'geschlossen'],a=1);

size(tt)
size(Aerg)

Code 0-2: Musterlösung

A=[1,2;3,4];
b=[19;43];

xy = inv(A)*b

Code 0-3: Lösung mit Hilfe von Scilab.

clear;

dt=0.001;
tmax = 10;
tt=0:dt:tmax;

// m*x.. = -C*x
//
// <=>
// x. = v
// v. = (-C/m)*x
//Numerische Lösung:
C=1; // [N/m]
m=1; // [kg]
x=1;
v=0;

Xerg = [x];
Verg = [v];

for t=dt:dt:tmax
   // "Neuer Wert = Alter Wert + aktuelle Steigung * dt"
         xneu     = x          + v                 * dt;
         vneu     = v          + (-C/m)*x          * dt; 

   x = xneu;
   v = vneu;

   Xerg = [Xerg,x]; //Sammeln der Ergebnisse
   Verg = [Verg,v]; //Sammeln der Ergebnisse
end

plot(tt,Xerg,'re',tt,Verg,'gr--');

//siehe help axes_properties
a = gca();
a.x_label.text = 't/s';
a.y_label.text = 'x/m';
a.title.text = 'Linearer Feder-Masse-Schwinger';
hl=legend(['x/m';'v/m*s^-1'],a=1);

size(tt)
size(Xerg)

Code 0-4: Simulation des Feder-Masse-Schwingers

clear;

dt=0.001;
tmax = 10;
tt=0:dt:tmax;

// m*x.. = -C*x
//
// <=>
// x. = v
// v. = (-C/m)*x + (-D/m)*v
//Numerische Lösung:
C=1; // [N/m]
D=1.0; // [Ns/m]
m=1; // [kg]
x=1;
v=0;

Xerg = [x];
Verg = [v];

for t=dt:dt:tmax
   // "Neuer Wert = Alter Wert + aktuelle Steigung * dt"
         xneu     = x          + v                 * dt;
         vneu     = v          + ((-C/m)*x + (-D/m)*v) * dt; 

   x = xneu;
   v = vneu;

   Xerg = [Xerg,x]; //Sammeln der Ergebnisse
   Verg = [Verg,v]; //Sammeln der Ergebnisse
end

plot(tt,Xerg,'re',tt,Verg,'gr--');

//siehe help axes_properties
a = gca();
a.x_label.text = 't/s';
a.y_label.text = 'x/m';
a.title.text = 'Linearer Feder-Masse-Dämpfer-Schwinger';
hl=legend(['x/m';'v/m*s^-1'],a=1);

size(tt)
size(Xerg)

Code 0-5: Lineares Dämpfungselement ergänzt.

m = 1.0;
C = 1.0;
D = 0.1;

function f = rechteSeite(t,y)
   x = y(1);
   v = y(2);
   f(1) = v;
   f(2) = (-C/m)*x + (-D/m)*v;
endfunction

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

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

Code 0-6: Lineare Feder-Masse-Schwinger mit Dämpfung, simuliert mit dem ode-Integrator.

m = 1.0;
C = 1.0;
D = 1.0;
P = 5.0;
xsoll = 3; //Sollwert für die Auslenkung 
vsoll = 0; //Sollwert für die Geschwindigkeit
function f = rechteSeite(t,y)
   x = y(1);
   v = y(2);
   f(1) = v;
   Fregler = P*(xsoll - x) + D*(vsoll - v);
   f(2) = (-C/m)*x + (1/m)*Fregler;
endfunction

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

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

Code 0-7: Lineare Feder-Masse-Schwinger mit Dämpfung und PD-Regler.

m = 1.0;
C = 1.0;
D = 1.0; //liefert FEHLER von 7.6701286
P = 5.0;
//D = 10.0; // liefert FEHLER von 3.9120228
//P = 50.0;
xsoll = 0; //Sollwert für die Auslenkung 
vsoll = 0; //Sollwert für die Geschwindigkeit
function f = rechteSeite(t,y)
   x = y(1);
   v = y(2);
   f(1) = v;
   Fregler = P*(xsoll - x) + D*(vsoll - v);
   f(2) = (-C/m)*x + (1/m)*Fregler;
endfunction

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

FEHLER = sqrt(y(1,:)*(y(1,:)')) //Skalaprodukt des Vektors aller Auslenkungswerte mit sich selbst, davon die Wurzel (Fehlerquadrat)

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

Code 0-8: Hinarbeiten auf eine automatische Optimierung: Festlegen einer s.g. Fehlerfunktion. Hier: Summe der Quadrate der Auslenkung als Kriterium, wie schlecht der Regler arbeitet.

//Geschlossene Lösung:
tt=0:0.1:2;
AG = exp(log(2)*tt);
plot(tt,AG);

A0=1;
dt=0.1;

function f = rechteSeite(t,y)
   A = y(1);
   f(1) = log(2)*A;
endfunction

t = 0:dt:2;
y0 = A0;
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);

plot(t,y(1,:)','re--');

Code 0-9: Übung: Seerosenteich mit ode integrieren.

clear;

dt=0.001;
tmax = 10;
tt=0:dt:tmax;

// m*x.. = -C*x + Freg
//
// <=>
// x. = v
// v. = (-C/m)*x + Freg/m
//Numerische Lösung:
C=1; // [N/m]
m=1; // [kg]
x=1;
v=0;

Xerg = [x];
Verg = [v];
e_alt = 0; // alte Regeldifferenz
for t=dt:dt:tmax
   x_soll = 0.2;
   P = 10.0;
   e = x_soll - x;

   eD = (e - e_alt)/dt;
   D  = 1.0;

   Freg = P*e + D*eD;

   // "Neuer Wert = Alter Wert + aktuelle Steigung * dt"
         xneu     = x          + v                 * dt;
         vneu     = v          + (  (-C/m)*x  + Freg/m  ) * dt; 

   x = xneu;
   v = vneu;

   Xerg = [Xerg,x]; //Sammeln der Ergebnisse
   Verg = [Verg,v]; //Sammeln der Ergebnisse

   e_alt = e;
end

plot(tt,Xerg,'re',tt,Verg,'gr--');

//siehe help axes_properties
a = gca();
a.x_label.text = 't/s';
a.y_label.text = 'x/m';
a.title.text = 'Linearer Feder-Masse-Schwinger';
hl=legend(['x/m';'v/m*s^-1'],a=1);

size(tt)
size(Xerg)

Code 0-10: PD-Regler beim Linearen Schwinger.

clear;

dt=0.001;
tmax = 30;
tt=0:dt:tmax;

// m*x.. = -C*x + Freg
//
// <=>
// x. = v
// v. = (-C/m)*x + Freg/m
//Numerische Lösung:
C=1; // [N/m]
m=1; // [kg]
x=1;
v=0;

Xerg = [x];
Verg = [v];
e_alt = 0; // alte Regeldifferenz
e_summe = 0; //aufsummierte Regelabweichung, um die numerische Integration von e zu bilden
for t=dt:dt:tmax
   x_soll = 0.2;
   P = 10.0;
   e = x_soll - x;

   eD = (e - e_alt)/dt;
   D  = 1.0;

   eI = e_summe * dt;
   I  = 5.0;

   Freg = P*e + I*eI + D*eD; //vollständiger PID-Regler

   // "Neuer Wert = Alter Wert + aktuelle Steigung * dt"
         xneu     = x          + v                 * dt;
         vneu     = v          + (  (-C/m)*x  + Freg/m  ) * dt; 

   x = xneu;
   v = vneu;

   Xerg = [Xerg,x]; //Sammeln der Ergebnisse
   Verg = [Verg,v]; //Sammeln der Ergebnisse

   e_alt = e;
   e_summe = e_summe + e; //Aufsummieren der Regelabweichung
end

plot(tt,Xerg,'re',tt,Verg,'gr--');

//siehe help axes_properties
a = gca();
a.x_label.text = 't/s';
a.y_label.text = 'x/m';
a.title.text = 'Linearer Feder-Masse-Schwinger';
hl=legend(['x/m';'v/m*s^-1'],a=1);

size(tt)
size(Xerg)

Code 0-11: PID-Regler beim Linearen Schwinger.

//Herausgezogene wichtige Parameter des Systems
boolean TAUSCHE_ANTRIEB_LINKS_RECHTS = false;
float VORTRIEB = 0.7;
float PROPORTIONALE_VERSTAERKUNG = 0.58;
float ASYMMETRIE = 1.01; // 1.0==voll symmetrisch, >1, LINKS STAERKER, <1 RECHTS STAERKER

Code 0-12: Anzupassende Parameter (Auszug des Processing-Sketches)

function f = rechteSeite(t,y)
    A = [ 0 , 1 ; -1 , 0 ];
    f = A*y; 
endfunction

t = 0:0.01:10;
y0 = [1;0];
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);

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

Code 0-13: Simulationsgleichung in Matrixschreibweise

m = 1; // kg
C = 1; // N/m
D = 1; // Ns/m
function f = rechteSeite(t,y)
    A = [ 0 , 1 ; -C/m , -D/m ];
    f = A*y; 
endfunction

t = 0:0.01:10;
y0 = [1;0];
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);

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

Code 0-14: Simulationsgleichung parametrisiert und mit Dämpfung in Matrixschreibweise

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

// G(s) = 1 / (s^2 + 2s + 1)
// x/u = 1 / (s^2 + 2s + 1)
// x*(s^2 + 2s + 1) = u
// x.. + 2x. + x = u
// x.. = -x -2x. + u

// x.  = v
// v.  = -x -2v + u

A = [0,1 ; -1, -2];

function f = rechteSeite(t,y)
    u = 1;
    f = A*y + [0;1]*u; 
endfunction

t = 0:0.01:10;
y0 = [0;0];
t0 = 0;
y2  = ode(y0,t0,t,rechteSeite);

plot(t,y(1,:)','gre',t,y2(1,:)','blu--');

Code 0-15: Simuilation des gleichen Systems erst mit csim (Laplace-Transformierte), dann mit ode.

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


A = [0,1 ; 0, -(1/3)];

function f = rechteSeite(t,y)
    u = 1;
    f = A*y + [0;(1/3)]*u; 
endfunction

t = 0:0.01:10;
y0 = [0;0];
t0 = 0;
y2  = ode(y0,t0,t,rechteSeite);

plot(t,y(1,:)','gre',t,y2(1,:)','blu--');

Code 0-16: Ergebnis zur Aufgabe 5.

//s als Argument für ein Polynom definieren:
s = poly(0,"s"); 
//Übertragungsfunktion definieren:
G = syslin('c',[1],[1+2*s+s^2]);

//PID-Regler nach Ziegler und Nichols Methode 2 bilden:
K  = 12.8; //(16.8)
TN = 1.4;  //(0.4)
TV = 0.1; 

R = syslin('c',[K*TN*s + K + K*TN*TV*s^2],[TN*s]);

H = R*G;

Q = H / (1 + H);

t=[0:0.01:10]; 
w=ones(1,1001); //Sollwert w ist Input für Gesamtübertragungsverhalten
y=csim(u,t,Q);  //Gesamtübertragungsfunktion des geschlossenen Systems simulieren.
//plot2d(t,y);
//plzr(G);

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

Code 0-17: Beispielumsetzung des zweiten Verfahrens nach Ziegler und Nichols mit Scilab

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

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

Code 0-18: Variante, Entwerfen Sie einen passenden PI-Regler.

-1.4655712 + 0.i      
   0.2327856 + 0.792552i
   0.2327856 - 0.792552i

Code 0-19: Eigenwerte für b), System ist instabil, also auch nicht technisch stabil.

clear;

m = 2.0;
C = 2.0;
D = 0.5;

function f = rechteSeite(t,y)
   x = y(1);
   v = y(2);
   F=10;
   f(1) = v;
   f(2) = (-C/m)*x + (-D/m)*v + (1/m)*F;
endfunction

t = 0:0.01:10;
y0 = [0,0]';
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);

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


//siehe help axes_properties
a = gca();
a.x_label.text = 't/s';
a.y_label.text = 'x/m';
a.title.text = 'Linearer Feder-Masse-Schwinger';
hl=legend(['x/m'],a=1);

Code 0-20: Lösung zu Übung 3a)

clear;

m = 2.0;
C = 2.0;
D = 0.5;

//s als Argument für ein Polynom definieren:
s = poly(0,"s"); 
//Übertragungsfunktion definieren:
G = syslin('c',[1],[m*s^2+D*s+C]);
t=[0:0.01:10]; 
F=10*ones(1,1001);
y=csim(F,t,G);

plot(t,y(1,:),'blu--');

//siehe help axes_properties
a = gca();
a.x_label.text = 't/s';
a.y_label.text = 'x und v';
a.title.text = 'Linearer Feder-Masse-Schwinger';
hl=legend(['x/m';'v/m*s^-1'],a=1);

Code 0-21: Lösung zu Übung 3b)

clear;

dt=0.01;
tmax = 10;
tt=0:dt:tmax;

m=2; // [kg]
C=2; // [N/m]
D=0.5; // [Ns/m]
x=0;
v=0;

Xerg = [x];
Verg = [v];

for t=dt:dt:tmax
   F=10.0;
   // "Neuer Wert = Alter Wert + aktuelle Steigung * dt"
   xneu     = x          + v                 * dt;
   vneu     = v          + ((-C/m)*x + (-D/m)*v + (1/m)*F) * dt; 

   x = xneu;
   v = vneu;

   Xerg = [Xerg,x]; //Sammeln der Ergebnisse
   Verg = [Verg,v]; //Sammeln der Ergebnisse
end

plot(tt,Xerg,'blu--',tt,Verg,'blk--');

//siehe help axes_properties
a = gca();
a.x_label.text = 't/s';
a.y_label.text = 'x/m';
a.title.text = 'Linearer Feder-Masse-Dämpfer-Schwinger';
hl=legend(['x/m';'v/m*s^-1'],a=1);

Code 0-22: Lösung zu Übung 3c)

clear;

m = 2.0;
C = 2.0;
D = 0.5;

function f = rechteSeite(t,y)
   f(1) = sin(t);
endfunction

t = 0:0.01:6;
y0 = 0;
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);

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

a = gca();
a.x_label.text = 't/s';
a.y_label.text = 'x';
a.title.text = 'dx/dt=sin(t)';
hl=legend(['x'],a=1);

Code 0-23: Lösung zu Übung 3d)

clear;

m = 2.0;
C = 2.0;
D = 0.5;

function f = rechteSeite(t,y)
   f(1) = sin(t);
endfunction

t = 0:0.01:6;
y0 = -1;
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);

plot(t,y(1,:),'gr',t,sin(t),'re');

a = gca();
a.x_label.text = 't/s';
a.y_label.text = 'x';
a.title.text = 'dx/dt=sin(t)';
hl=legend(['x','sin(t)'],a=1);

Code 0-24: Variante zu Lösung zu Übung 3d) mit x0=-1

clear();

xtot=0;

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

       Kkrit = 9.7;       
       xsoll = 1.0;
       //e     = xsoll - x;
       e     = xsoll - xtot;
       u=Kkrit*e;
       
       f(1,1)=v;
       f(2,1)=-0.5*x-v+u; //angepasst!!!
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 = 99.0;
dt   = 0.01;
schritte = ceil(tmax/dt);

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

Ttot = 0.1; //angepasst!!!
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-25: Angepasstes Scilab-Skript.

P-Regler:    K=0.5*Kkrit=4.85
PI-Regler:   K=0.45*Kkrit=4.365 und TN=0.85*Tkrit=1.717
PID-Regler:  K=0.6*Kkrit=5.82   und TN=0.5 *Tkrit=1.01  und TV=0.12*Tkrit=0.2424

Code 0-26: Parameter für die verschiedenen Regler

R=ppol([0,1;-0.5,-1],[0;1],[-2+0.5*%i,-2-0.5*%i])

Code 0-27: Lösung mit Scilab, liefert R=[3.75,3]

R=ppol([0,1;-0.5,-1],[0;1],[-1+0.5*%i,-1-0.5*%i])

Code 0-28: Lösung mit Scilab, Variante, liefert R=[0.75,1]

det|A* - lambda*E|=0

Code 0-29: ...gilt es jetzt zu lösen.

-0.5-0.5r2=-1 => r2=1
in Wurzel eingesetzt ergibt:
1/4 = 1/2 +r1 -1 => r1=3/4

Code 0-30: Bestimmung der Werte von r1 und r2.