Zustandsraumdarstellung

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ausarbeitung zum Thema Prozessanalyse Modellierung Matlab 7.1
% Vorlesung: Entwurfmethoden für Regler
% Copyright 2007 Thomas Ludwig, Slawa Braun
% This file is distributed under the terms of the GNU General Public License,
% see http://www.gnu.org/copyleft/gpl.txt
%
% Beschreibung: Modellieren eines Systems in die Zustandsraumdarstellung
%               am Beispiel eines Gleichstrommotors.
%
%               1. Systemzerlegung - Zerlege das System in Komponenten
%               2. Beschreibe das Verhalten der Komponenten
%                  (physikalische, chemische... Grundlagen)
%               3. Beschreibe die Beziehungen zwischen den Komponenten
%               4. Fasse Gleichungen zu einer DGL zusammen
%               5. Normieren
%               6. DGL in Matrixschreibweise der Zustandsraumdarstellung bringen
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function gleichstrommotor();

    % Allgemeine Optionen
    format long;
    close all;
    clear all;
    clc

    help gleichstrommotor
   
    % Daten laden
    obj= daten();

    % Kovertiere Daten nach SS, TF, FRD und ZPK
    obj= convert( obj );   
   
    % System abfragen
    obj= identify( obj );
   
    % Ausgabe im Command Window
    ausgabe( obj );   
    
    % Grafiken ausgeben
    switch input('  >>  Diagramme mit LTIView [L], SISOTOOL [S] oder einzeln [E] generieren?  ','s');
        case 'L'
            ltiview( obj.ss.sys );
        case 'S'
            sisotool( obj.ss.sys );
        case 'E'
            disp('Die Grafiken werden nun generiert....');
            grafik( obj );           
    end
   
    % Regler Vergleich
    %regler( obj );
   
    % Objekt an Workspace uebergeben
    assignin( 'base', 'obj', obj );
   
end

%% Daten - Fuer Zustandsraum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function obj= daten();
    %Parameter der Aufgabenstellung
    R = 0.5;    %Ohms
    L = 0.05;   %mHenrys
    Km = 0.03;    %Vs/rad
    Kt = 30*10^(-3);    %Nm/A
    Kr = 30*10^(-6);    %Nms
    J = 10*10^(-4);     %Nms^2
    % Zustandsdifferentialgleichung: d/dt*x(t)= A*x(t) + b*u(t)
    % Ausgangsgleichung: y(t)= c'*x(t) + d*u(t)
    % x(t) - Zustandsvektor
    % u(t) - Eingangsvariable
    % y(t) - Ausgangsvariable
   
    A=[ -R/L     -Km/L;
         Kt/J    -Kr/J];   % Systemmatrix

    b=[ 1/L;
         0 ];               % Eingangsvektor

    c=[ 0 1 ];              % Ausgangsvektor

    d=[ 0 ];                % Durchgangsfaktor
   
    % Ts = 0 - Continuous-Time State-Space Model
    % Ts > 0 - Discrete-Time Model
    Ts= 0;                  % Abtastzeit: 0.01s -> Ts= 0.01
   
    t = 0:Ts:10;            % Simulationszeit: 0~10s
   
    % Daten an Objekt uebergeben
    obj.A= A;
    obj.b= b;
    obj.c= c;
    obj.d= d;
    obj.t= t;
    obj.Ts= Ts;

end

%% Konvertierung
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function obj= convert( obj );

    % Zustandsraum
    obj.ss.sys = ss(obj.A, obj.b, obj.c, obj.d, obj.Ts);
    obj.ss.tfsys= tf( obj.ss.sys );

    % Frequenzgang-Daten Modell
    obj.frd.sys= frd(obj.ss.sys, obj.Ts,'Units','Hz');
   
    % Uebertragungsfunktion
    [obj.tf.zaehler, obj.tf.nenner]= ss2tf(obj.A, obj.b, obj.c, obj.d);     % Zaehler / Nenner
    obj.tf.GS=  tf(obj.tf.zaehler, obj.tf.nenner, obj.Ts);

    % Pol-Nullstellenmodell
    [obj.zpk.z, obj.zpk.p, obj.zpk.k] = ss2zp(obj.A, obj.b, obj.c, obj.d);  % Zero / Pole / Gain
    obj.zpk.GS= zpk(obj.zpk.z, obj.zpk.p, obj.zpk.k, obj.Ts);
   
end

%% System Identifizieren
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function obj= identify( obj );

    % Nullstellen
    obj.nullstellen= zero(obj.ss.sys);

    % Polstellen
    obj.polstellen= pole(obj.ss.sys);

    % Eigenwerte
    [obj.eigenvektor, obj.eigenwerte]= eig(obj.ss.sys.a);

    % Dämpfung
    [obj.daempfung.Wn, obj.daempfung.Z, obj.daempfung.P]= damp(obj.ss.sys); % Spuckt mehr aus als bei obj.ss.sys
    obj.daempfung.matrix= [obj.daempfung.Z, obj.daempfung.P, obj.daempfung.Wn];

    % Verstärkungsfaktor
    obj.kp= dcgain(obj.ss.sys);
   
    % Amplitudenrand: GainMargin und GMFrequency
    % Phasenrand: PhaseMargin und PMFrequency
    % Totzeitrand: DelayMargin und DMFrequency
    % Stabilitat: Stable, 1 bei Stabilitat, 0 sonst.       
    obj.allmargin= allmargin(obj.ss.sys);
   
    % Bode Daten
    [obj.bode.mag, obj.bode.phase, obj.bode.w]= bode(obj.ss.sys);

    % Nyquist Daten
    [obj.nyquist.mag, obj.nyquist.phase, obj.nyquist.w]= bode(obj.ss.sys);
   
    % Wurzelortskurve Daten
    [obj.Wok.r, obj.Wok.k]= rlocus(obj.ss.sys);

end

%% Ausgabe im Command Window
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ausgabe( obj );

    % System ausgeben
    disp('Das gewählte Zustandsmodell:')
    printsys(obj.A, obj.b, obj.c, obj.d)

    % Abtastzeit
    disp(sprintf('Die gewählte Abtastzeit Ts=%ds', obj.Ts))
   
    % Single Input, Single Output (SISO)
    disp(sprintf('Ist dies ein Single Input Single Ouput System? [J/N]: %i', issiso(obj.ss.sys)))

    % Verstärkungsfaktor
    disp('-> Verstärkungsfaktor'), disp(sprintf('   Kp= %f', obj.kp));

    % Uebertragungsfunktion - TF
    disp('-> Übertragungsfunktion - TF'), obj.tf.GS
    disp('********************************************************************')

    % Übertragungsfunktion - SS
    disp('-> Übertragungsfunktion - SS - Zustandsdarstellung'), obj.ss.tfsys
    disp('********************************************************************')
   
    % Übertragungsfunktion - ZPK
    disp('-> Übertragungsfunktion - ZPK - Nullstellen-Polstellen-Darstellung'), obj.zpk.GS 
    disp('********************************************************************')

    % Frequenz
    disp('-> Frequenz Daten Modelle - FRD - Ausgang'), obj.frd.sys
    disp('********************************************************************')   
   
    % Nullstellen
    disp('-> Nullstellen'), disp(obj.nullstellen)
   
    % Polstellen
    disp('-> Polstellen - Nullstellen der Charakterlichen Gleichung'), disp(obj.polstellen)
   
    % Eigenwerte
    disp('-> Eigenwerte'), disp(obj.eigenwerte)
   
    % Eigenvektor
    disp('-> Eigenvektor'), disp(obj.eigenvektor)

    % Daempfung
    disp('-> Dämpfung           Eigenwerte          Frequenz (rad/s)'), disp(obj.daempfung.matrix)
  
    % Stabilitaet
    disp('-> Amplitudenrand/Phasenrand/Totzeitrand/Stabilität'), disp(obj.allmargin)
   
end

%% Grafik
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function grafik( obj );

    % Bode-diagramm
    figure('Name','Bode Diagramm','NumberTitle','off'), ch_bode= bodeplot(obj.ss.sys,'r');
    setoptions(ch_bode,'FreqUnits','Hz');
   
    % Nyquist-Kurve
    figure('Name','Nyquist','NumberTitle','off'), ch_nyquist= nyquistplot(obj.ss.sys);      
    setoptions(ch_nyquist,'FreqUnits','Hz', 'PhaseUnits','deg', 'MagUnits','dB', 'ShowFullContour','on', 'Grid','on');
   
    % Wurzelortskurve
    figure('Name','Wurzelortskurve','NumberTitle','off'), rlocusplot(obj.ss.sys), sgrid;       

    % Nichols
    %figure('Name','Nichols Kurve','NumberTitle','off'), nicholsplot(obj.ss.sys);
   
    % Null-Polstellen Verteilung
    figure('Name','Null-Polstellen Verteilung','NumberTitle','off'), ch_pzplot= pzplot(obj.ss.sys);      
    setoptions(ch_pzplot,'FreqUnits','Hz'), sgrid;;
   
    % Sprungantwort
    figure('Name','Sprungantwort','NumberTitle','off'), stepplot(obj.ss.sys);

    % Impulsantwort
    figure('Name','Impuslantwort','NumberTitle','off'), impulseplot(obj.ss.sys);

    % Sigma
    figure('Name','Sigma','NumberTitle','off'), sigmaplot(obj.ss.sys);

    % Initial
    %figure('Name','Initial [1; 0]','NumberTitle','off'), initialplot(obj.ss.sys, [1; 0]);

    % Hankel singular values
    %figure('Name','Hankel singular values','NumberTitle','off'), hsvplot(obj.ss.sys);

end

% Ueberpruefung der Stabilitätskriterien
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function obj= stable_check( obj );

    % Dämpfung überprüfen
    if (D > 1)
        status= 'Aperiodischer Fall';
    elseif (D == 1)
        status= 'Aperiodischer Grenzfall';
    elseif ( (D > 0) && (D < 1) )
        status= 'Periodischer Grenzfall - Schwingfall';
    end

end