//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.

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.

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.

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.

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

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.

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.

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.

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

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.