//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: Verwendung von Scilab zur Handhabung von Übertragungsfunktionen (Beispiel)

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

t = linspace(0,30,3000);
y0 = [0,0]';
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);

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

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

// Laplace-Darstellung:


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

u=10*ones(1,3000);  //entspricht F
y=csim(u,t,G);
plot(t,y,'blu--');
//plzr(G);

Code 0-2: Feder-Masse-Schwinger mit äußerer Kraft F im Zeit- und Laplacebereich simulieren.

m = 1.0;
C = 1.0;
D = 1.0;
F = 10.0;
K = 10.0;   //erreicht 9,5 (näher an Sollwert, stärkeres Schwingen zu Beginn)
//K = 1.0;  //erreicht 5 (halb so groß wie Sollwert)
// Laplace-Darstellung:


//s als Argument für ein Polynom definieren:
s = poly(0,"s"); 
//Übertragungsfunktion definieren:
G = syslin('c',[1],[m*s*s + D*s + C]);
R = syslin('c',[K],[1]);
H = R*G;
Q = R*G/(1+R*G);

t = linspace(0,30,3000);
u=10*ones(1,3000);  //entspricht xsoll
y=csim(u,t,Q);
plot(t,y,'blu--');
//plzr(Q);  //Polstellen des Gesamtsystems anzeigen.

Code 0-3: Linearer Schwinger mit Dämpfung mit P-Regler, Realisierung im Laplace-Bereich mit Scilab.

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

t = linspace(0,30,3000);
y0 = [0,0]';
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);

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

Code 0-4: Linearer Schwinger mit Dämpfung mit P-Regler, Realisierung im Zeit-Bereich mit Scilab.

clear();
disp("Schritt 1: Eigenwerte des linearisierten Systems:");

//  phi.. = phi * mgh/J + FA * h/J

//  =>   phi. = om
//       om.  = phi * mgh/J + FA * h/J

//  =>

m = 1;  //kg
g = 9.81; //m/s^2
h = 0.5; //m
r = 0.1; //m
l = 2*h;
J = 0.25*m*r^2 + (1/12)*m*l; // kg*m^2

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

B = [0;h/J]

lambda = spec(A);

disp("EW ohne Regler:");
disp(lambda);  //   lambda1/2=  +/- 7.5594729  (instabil)

disp("Schritt 2: Bestimmung der Parameter des Zustandsreglers mit ppol:");

//R=ppol(A,B,EW),   [phi.; om.] = A*y - B*R*y
//(vergl. Kapitel 8 unten)
//Ergebnis überprüfen:
EW = [-1+%i,-1-%i]
R=ppol(A,B,EW)  // == 10.153333    0.3433333  
spec(A-B*R)

disp("Testen des Refglers im Originalsystem: siehe inverspendel_test.sce");

//  FA = -10.153333*phi -0.3433333*om

Code 0-5: inverspendel.sce

clear();
disp("Verwenden der Regelparameter im Originalsystem:");

m = 1;  //kg
g = 9.81; //m/s^2
h = 0.5; //m
r = 0.1; //m
l = 2*h; //m
J = 0.25*m*r^2 + (1/12)*m*l; // kg*m^2

function f = rechteSeite(t,y)
    phi = y(1,1);
    om  = y(2,1);
    f(1,1) = om;
    FA = -10.153333*phi -0.3433333*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 = [0.2,0.1]';
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);

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

Code 0-6: inverspendel_test.sce

R1=10.153333;
R2=0.3433333;

m = 1;  //kg
g = 9.81; //m/s^2
h = 0.5; //m
r = 0.1; //m
l = 2*h; //m
J = 0.25*m*r^2 + (1/12)*m*l; // kg*m^2

function f = rechteSeite(t,y)
    phi = y(1,1);
    om  = y(2,1);
    f(1,1) = om;
    FA = -R1*phi -R2*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

Code 0-7: modell.sce

function z = berechneFehler(RR1,RR2)
    R1 = RR1;
    R2 = RR2;

    t = linspace(0,20,3000);
    y0 = [0.2,0.1]';
    t0 = 0;
    yy  = ode(y0,t0,t,rechteSeite);

    fehler = sum(abs(yy(1,:)));    
    z = fehler;    
endfunction

//Wie oben, aber mit Ausgaben:
function z = berechneFehlerTest(RR1,RR2)
    R1 = RR1;
    R2 = RR2;

    t = linspace(0,20,3000);
    y0 = [0.2,0.1]';
    t0 = 0;
    yy  = ode(y0,t0,t,rechteSeite);

    fehler = sum(abs(yy(1,:)));  
    disp(fehler);
    plot(t,yy(1,:)',t,yy(2,:)');  
    z = fehler;    
endfunction

Code 0-8: fehlerfunktion.sce

clear();

exec modell.sce
exec fehlerfunktion.sce


RRR1_start = 10;
RRR2_start = 0.3;

RRR1 = RRR1_start;
RRR2 = RRR2_start;

ALFA1 = 0.01;
ALFA2 = 0.01;

FAKTOR = 1.01;

aktueller_fehler = berechneFehler(RRR1,RRR2);

for i=1:1000  //innere Befehle iterativ ausführen
    dx = floor(grand(1, 1, "unf", 0, 2.999999999)) - 1;  //liefert -1 oder +1 zufällig
    dy = floor(grand(1, 1, "unf", 0, 2.999999999)) - 1;  //liefert -1 oder +1 zufällig

    RRR1_neu = RRR1 + ALFA1*dx;
    RRR2_neu = RRR2 + ALFA2*dy;

    neuer_fehler = berechneFehler(RRR1_neu,RRR2_neu);

    //Prüfen, ob neuer Fehler kleiner, wenn ja, zu neuem Punkt gehen, sonst bei altem bleiben    
    if neuer_fehler<=aktueller_fehler then
        RRR1 = RRR1_neu;
        RRR2 = RRR2_neu;
        if neuer_fehler<=aktueller_fehler then
            disp('fehler='+string(aktueller_fehler));     
            disp('RRR1='+string(RRR1));     
            disp('RRR2='+string(RRR2));     
        end
        aktueller_fehler = neuer_fehler;
        if dx>0 then
            ALFA1 = ALFA1 * FAKTOR;
        end
        if dx<0 then
            ALFA1 = ALFA1 * FAKTOR;
        end
        if dy>0 then
            ALFA2 = ALFA2 * FAKTOR;
        end
        if dy<0 then
            ALFA2 = ALFA2 * FAKTOR;
        end
    else
        if dx>0 then
            ALFA1 = ALFA1 * (1/FAKTOR);
        end
        if dx<0 then
            ALFA1 = ALFA1 * (1/FAKTOR);
        end
        if dy>0 then
            ALFA2 = ALFA2 * (1/FAKTOR);
        end
        if dy<0 then
            ALFA2 = ALFA2 * (1/FAKTOR);
        end         
    end
end

berechneFehlerTest(RRR1,RRR2);

Code 0-9: optimierer.sce