Seminar: "Simulation dynamischer Systeme mit Matlab und Simulink"

15. Februar 2004

Marco Fütterer

Gebäude Uni 07/112

Tel.: 0391/ 67-11320

e-mail : marco.fuetterer@ovgu.de

Das folgende Skript soll eine selbständige Einarbeitung in Matlab Simulink aus der Sicht des Regelungstechnikers ermöglichen. In fünf Seminarstunden ist eine vollständige Einführung in dieses Softwarepaket nicht möglich. Deswegen werde ich hier nur bestimmte Techniken vorstellen, die für die Simulation und Synthese von Systemen notwendig sind. Vorschläge sowie Hinweise auf Fehler, nehme ich gerne entgegen, um dieses Skript weiter zu verbessern.

Verwendet den Internet Explorer zur korrekten Anzeige der Formeln!

Inhaltsverzeichnis

1 Seminar
    1.1 Einführung in die Matlab- Befehlszeile
    1.2 Der Aufruf von Simulink
    1.3 Simulation eines RC- Gliedes in Simulink
    1.4 Simulationsparameter in Matlab m-File organisieren
    1.5 Simulation eines aktiven Tiefpasses in Simulink
    1.6 Einstellung der Simulationsparameter
    1.7 Simulation eines Simulink- Modells von der Kommandozeile
    1.8 Analyse eines Simulink- Modells von der Kommandozeile
    1.9 Hausaufgabe
    1.10 Auflösung der Hausaufgabe
2 Seminar
    2.1 Satellitenbewegung
        2.1.1 M- function zur DGL- Beschreibung
        2.1.2 Simulationsfehler
    2.2 Simulation eines RLC- Gliedes
        2.2.1 M- function zur DGL- Beschreibung
        2.2.2 Parameterübergabe
    2.3 Nichtlinearer Oszillator
        2.3.1 Darstellung der nichtlinearen Kennlinie
        2.3.2 Simulation des van der Pol'schen Oszillators
        2.3.3 Phasenebene des van der Pol'schen Oszillators
        2.3.4 Simulation mit m = 0
        2.3.5 Simulation mit unterschiedlichen m
        2.3.6 Simulation mit m=100
        2.3.7 Theoretische Analyse
3 Seminar
    3.1 Simulation eines Gleichstrommotors
        3.1.1 Herleitung der Differentialgleichung
        3.1.2 Statische Beziehungen
        3.1.3 Berechnung der Kennwerte
        3.1.4 Die Drehzahl- Drehmoment- Kennlinie
        3.1.5 Der Signalflussplan in Simulink
        3.1.6 Blockmaskierung in Simulink
        3.1.7 Bode- Diagramm
        3.1.8 Arbeitspunkteinstellung
        3.1.9 Sprungantworten
        3.1.10 Drehzahlregelung mit P- Regler
        3.1.11 Drehzahlregelung mit PI- Regler
        3.1.12 Unterlagerte Stromregelung
    3.2 Gleichstrommotor mit Erregerfeld
        3.2.1 Herleitung der Differentialgleichung
        3.2.2 Statische Beziehungen
        3.2.3 Reibkennlinie
        3.2.4 Erregerflusskennlinie
        3.2.5 Der Signalflussplan in Simulink
        3.2.6 Linearisierung im Arbeitspunkt
    3.3 Hausaufgabe
4 Seminar
    4.1 Die Symbolic Math Toolbox
        4.1.1 Definition von symbolischen Variablen
        4.1.2 Linearisierung im Arbeitspunkt
        4.1.3 Symbolische Lösung der algebraischen Gleichungen
        4.1.4 Numerische Lösung
    4.2 Die Control- Toolbox
        4.2.1 Definition von Transferfunktionen
        4.2.2 Definition von Zustandsraumdarstellungen
        4.2.3 Übergangsfunktion und Gewichtsfunktion
        4.2.4 Bode- Diagramm, Ortskurve und Pole- Zero- Darstellung
        4.2.5 Diskretisierung von kontinuierlichen Systemen
        4.2.6 Simulation von Systemen
5 Seminar
    5.1 S- Function
        5.1.1 Aufruf Syntax
        5.1.2 S- Function am Beispiel des linearen Gleichstrommotors
        5.1.3 Parameterübergabe und Blockmaskierung
    5.2 Klassische Reglersynthese mit Hilfe des SISO- Tools
        5.2.1 5.2.1 Einführung in das Wurzelortsverfahren
        5.2.2 5.2.2 Konstruktionsregeln für das Wurzelortsverfahren
        5.2.3 Aufruf des SISO- Tools
        5.2.4 P- Regler
        5.2.5 P- Regler mit Verzögerungsglied
        5.2.6 I- Regler
        5.2.7 I- Regler mit Frequenzkorrektur

1 Seminar

1.1 Einführung in die Matlab- Befehlszeile

Nach dem Starten von Matlab erscheint das Befehlszeilenfenster, welches das Kernstück von Matlab ist. Das Prompt > > signalisiert die Eingabebereitschaft. Wie bei einem symbolischen Taschenrechner, kann man mit der Tastatur einfache mathematische Ausdrücke eingeben, und diese sich mit der Betätigung der ,,Enter" Taste berechnen lassen. Dabei gilt für die Berechnung ,,Punkt vor Strich". Klammerausdrücke werden wie gewohnt zuerst berechnet.

Beispiel:

>> (22+15*44)^(2/5)

13.5992

Die Ausgabegenauigkeit ist auf 4 Stellen nach dem Punkt voreingestellt. Mit dem Befehl format, kann diese auf 14 Stellen nach dem Punkt eingestellt werden.

Beispiel:

>> format long

>> (22+15*44)^(2/5)

13.59921218019641

Auch die Ausgabeform kann mit dem Befehl format geändert werden, z.B. in die Exponentialdarstellung.

Beispiel:

>> format short e

>> (22+15*44)^(2/5)

1.3599e+001

Für alle Befehle existiert eine Kurzhilfe, die man durch den Befehl help mit dem nachfolgenden eigentlichen Befehl aufrufen kann. Damit lässt sich sofort der richtige Syntax eines Befehls ermitteln.

Beispiel:

>> help format

FORMAT Set output format.

All computations in MATLAB are done in double precision.

.....

Umfangreichere Hilfen lassen sich durch die Eingabe von helpwin, helpdesk oder lookfor aufrufen. Grosse Zahlen können auch jederzeit in Exponentialdarstellung eingegeben werden.

Beispiel:

>> format short

>> 1e3

1000

Vorteilhaft lässt sich der Befehlszeilenspeicher verwenden. Durch Verwendung der Cursortasten für hoch und runter, kann man sich durch die letzten Eingaben navigieren und eine bestimmte Eingabe durch drücken der ,,Enter" -Taste wiederholen. Gibt man die ersten Buchstaben des Befehls bereits vor, so werden nur alle diejenigen Befehle angezeigt die ebenfalls mit dieser Buchstabenkombination anfangen. Es lassen sich auch die vorhergehenden Befehlszeileneingaben durch die rechts und links Cursortasten bearbeiten.

Das Beenden der Matlab- Sitzung geschieht durch die Eingabe von quit oder exit.

1.2 Der Aufruf von Simulink

Matlab besteht aus einer Ansammlung von Toolboxen. Die wohl populärste Toolbox ist Simulink, welche eine Blockorientierte Eingabe von Simulationsmodellen erlaubt. Der Aufruf dieser Toolbox erfolgt durch die Eingabe von Simulink in der Matlab Befehlszeile.

Beispiel:

>> simulink

Anschließend öffnet sich der Simulink Block Library Browser. Dieser besitzt eine Baumstrukturansicht, in der die bereitgestellten Simulink Blöcke angeordnet sind. Ein neues leeres Blockdiagramm lässt sich erzeugen, indem man im Menu des Simulink Block Library Browser den Befehl File® New® Model anwählt. Ein neues leeres Fenster mit der Bezeichnung ,,Untitled" öffnet sich. Aus dem Simulink Block Library Browser lässt sich nun ein gewünschter Block auswählen und durch drücken der linken Maustaste in das leere Blockdiagramm kopieren. Der Name eines jeden Blockes im Blockdiagramm muss zur Identifizierung eindeutig sein. Zunächst wird von Simulink der eigentliche Blockname zugewiesen. Bei Auswahl des gleichen Blocktyps wird lediglich noch eine Zahl angefügt. Durch doppelklicken auf dem Blocknamen lässt sich über die Tastatur jedem Block einen benutzerdefinierten Namen zuweisen.

Bereits im Blockdiagramm vorhandene Blöcke können einfach durch anwählen, drücken der rechten Maustaste, verschieben der Maus und loslassen der rechten Maustaste kopiert werden.

Mitunter ist es wünschenswert Blöcke zu drehen oder zu spiegeln. Solche Manipulationen können durch Aufruf des Kontextmenüs des gewünschten Blockes gemacht werden. Der Aufruf des Kontextmenüs geschieht, indem der Mauszeiger auf dem Block positioniert wird und anschließend die rechte Maustaste gedrückt. Aus dem Untermenü Format lassen sich die gewünschten Befehle auswählen. Das Kopieren und Einfügen von Blöcken ist ebenfalls durch das Kontextmenü möglich.

Das Löschen eines Blockes erreicht man durch das Markieren dieses Blockes durch kurze Anwahl mit einem linken Mausklick und anschließender Betätigung der ,,Delete" Taste auf der Tastatur.

Mehrere Blöcke lassen sich gleichzeitig markieren, indem die linke Maustaste gedrückt und über den entsprechenden rechteckigen Ausschnitt im Blockdiagramm gezogen wird.

Die Ein- und Ausgänge der Blöcke lassen sich untereinander verbinden, indem der Mauszeiger über den Ausgang eines Blockes gehalten wird, anschließend die linke Maustaste gedrückt, den Mauszeiger zum Eingang des nächsten Blockes gezogen und losgelassen wird.

1.3 Simulation eines RC- Gliedes in Simulink

Der reale Prozess muss in eine geeignete Form mathematisch beschrieben werden. Für die Modellbildung des RC- Gliedes verwenden wir hier die Methode der Operatorimpedanzen. Gesucht ist das elektrische Verhalten zwischen Ein- und Ausgang. Bei der Methode der Operatorimpedanzen werden die Bauelemente als komplexe Widerstände angesehen. Durch die Anwendung der Laplace- Transformation lässt sich dieses erreichen. Unter Vernachlässigung der Anfangsbedingung, entsteht der folgende Zusammenhang zwischen Zeit- und Bildbereich.

Zeitbereich	Bildbereich
R	R
L·[d/dt]( ... )	s·L
[1/C]·ò₀^t ...dt	[1/(s·C)]

Die Übertragungsfunktion eines Systems ist definiert durch das Verhältnis zwischen Ausgangssignal und Eingangssignal im Bildbereich. Betrachtet man das RC- Glied nun als komplexen Spannungsteiler lässt sich das Modell eines RC- Gliedes wie folgt herleiten:

G(s)=

U_a (s)

U_e (s)

[1/(s·C)]

R+

s·C

s·R·C+1

Für die Transferfunktion findet man im Simulink Block Library Browser unter Continuous den ,,Transfer Fcn"- Block. Dieser wird aus der Library entnommen und in ein neues leeres Blockdiagramm eingefügt. Durch doppelklicken auf diesen Block wird das Blockparametermenü aufgerufen. In den Feldern ,,Numerator" und ,,Denominator" lassen sich die Koeffizienten der Zähler- und Nennerpolynome in absteigender Ordnung als Zeilenvektor eingeben. Für unser RC- Glied geben wir im Feld Numerator [ 1] und im Feld Denominator [R·C 1] ein. Anschließend wird das Menü wieder geschlossen. Nun soll noch ein Einheitssprungsignal am Eingang des RC- Gliedes angeschlossen werden. Dazu wird aus dem Simulink Block Library Browser von Simulink® Sources der Block Step entnommen und vor dem Transferfunktionsblock platziert. Der Ausgang des Step Blockes wird mit dem Eingang des RC- Gliedes verbunden. Zur Anzeige des Ausgangssignales wird aus Simulink® Sinks der Block Scope entnommen. Der Ausgang des Transferfunktionsblocks wird mit dem Eingang des Scope Blocks verbunden. Das Ergebnis ist untenstehend abgebildet.

Bevor die Simulation gestartet werden kann müssen noch die numerischen Werte für R und C eingegeben werden. Dies geschieht in der Matlab Kommandozeile.

>> R=1000

>> C=100e-6

Dadurch werden zwei Variablen R und C definiert. Diese Variablen werden beim Start der Simulation vom Transferfunktionsblock übernommen, da wir R und C dort als Koeffizient definiert hatten.

Die Simulation kann nun durch den Menüpunkt Simulation® Start gestartet werden. Durch doppelklicken auf dem Scope- Block, lässt sich das Ergebnis grafisch ersehen. Nach zwei Sekunden ist der Übergang bereits abgeschlossen. Unter dem Menüpunkt Simulation® Simulation parameters... lässt sich ein Menü mit mehreren Reitern öffnen. Unter dem Reiter Solver kann die Start- und Stoppzeit angegeben werden. Auch hier kann man anstatt eines numerischen Wertes eine Variable z.B. tend für Stopp time angeben. In der Matlab Kommandozeile wird dann wieder die Variable tend definiert.

>> tend=2

Bei erneutem Start der Simulation, wird diese nach zwei Sekunden beendet.

1.4 Simulationsparameter in Matlab m-File organisieren

Oft werden Simulationen mehrfach wiederholt, so dass es sinnvoll ist alle Simulationsparameter in eine Script- Datei zu organisieren. Hierfür wird in der Matlab Befehlszeile der Befehl edit eingegeben.

>> edit

Danach öffnet sich ein Matlab eigener Texteditor. In diesem werden nun alle Kommandos geschrieben die zuvor in die Befehlszeile direkt eingegeben wurden.

R=1000;

C=100e-6;

tend=2;

Jede Befehlszeile wird nun durch ein Semikolon abgeschlossen. Damit wird verhindert, dass später bei Ausführung des Skriptes in der Matlab Kommandozeile keine Bestätigungen für die Variablenzuweisungen angezeigt werden. Anschließend wird das m-File durch den Menüpunkt File® Save unter den Namen RC_param abgespeichert. Wird nun in der Kommandozeile RC_param eingegeben werden alle Befehle die im m- File stehen ausgeführt.

Hinweis: Das Simulink- Modell wird mit einer Dateierweiterung *.mdl abgespeichert, während die Scriptdatei mit der Dateiendung *.m abgespeichert wird. Gibt man für beides den gleichen Dateinamen an, wird Matlab immer nur das m- File ausführen. Sind beide Namen unterschiedlich, lässt sich das Simulink-Modell auch direkt von der Kommandozeile aufrufen.

1.5 Simulation eines aktiven Tiefpasses in Simulink

Untenstehend ist der Schaltplan eines aktiven Tiefpasses angegeben. Es soll wieder mit der Methode der Operatorimpedanzen das mathematische Modell hergeleitet werden.

Bei der Modellbildung von elektrischen Schaltungen mit Operationsverstärker geht man häufig von einem idealisierten OPV aus. Das heißt man nimmt an, dass deren Verstärkung unendlich ist. Diese Annahme erlaubt es uns weiterhin anzunehmen, dass zwischen invertierten und nicht invertierten Eingang des OPV's kein Spannungsabfall entsteht. Eine zusätzliche Annahme ist, dass die Eingangsströme null sind. Bei der Analyse wird hier ein Hilfspunkt H mit der Hilfsspannung U_h eingeführt. Als erstes wird die Übertragungsfunktion der umrahmten Teilschaltung analysiert.

U_a (s)

U_h (s)

[1/(s·C₂ )]

R₂

s·R₂ ·C₂

U_h (s)=-s·R₂ ·C₂ ·U_a (s)

Nun wird die Kirchhoffsche Knotenpunktregel im Punkt H angewandt.

U_e -U_h

R₁

U_h

[1/(s·C₁)]

U_h

R₂

U_h -U_a

R₃

Wird die Gleichung zwei in drei eingesetzt, U_e und U_a ausgeklammert sowie die Gleichung nach dem Verhältnis U_a zu U_e umgestellt, erhält man folgende Übertragungsfunktion:

G(s)=

U_a (s)

U_e (s)

R₃

R₁

s²·R₂·R₃ ·C₁ ·C₂ +s·C₂ ·

æ
è

R₂ ·R₃

R₁

+R₂ +R₃

ö
ø

Die Koeffizienten des Zähler- und Nennerpolynoms lassen sich wieder als Funktionen der Nennwerte der Bauelemente bestimmen. Es ist hier günstig im Transferfunktionsblock von Simulink die allgemeinen Koeffizienten b₀, a₀, a₁ und a₂ als Variablen einzuführen. Die Zuordnung der Kennwerte der Bauelemente zu den Koeffizienten, wird in einer Script- Datei vorgenommen.

G(s)=

b₀

a₂ ·s²+a₁ ·s·+a₀

Für die Simulation des aktiven Tiefpasses soll nun ein Pulse Generator mit einer Periode von einer Sekunde und einer Pulsbreite von 50% verwendet werden. Dieser Block befindet sich ebenfalls in Simulink® Sources. Die Parameter werden im Parameterdialog des Blockes eingestellt.

Schreiben Sie die notwendige Parameterdatei und simulieren Sie die Schaltung des aktiven Tiefpasses.

Beim doppelklicken auf dem Scope- Block ist zu ersehen, dass das Ergebnis nur noch von 9,2s bis 10s dargestellt wurde. Dies liegt an der Voreinstellung des Scope- Blockes, der nur die letzen 5000 Simulationswerte speichert. In der Symbolleiste befindet sich das Icon ,,Parameters". Nach dem anwählen öffnet sich die Parameterdialogbox von Scope. Klicken Sie auf den zweiten Reiter Namens ,,Data history". Deaktivieren Sie die Check- Box ,,Limit Data points to last". Simulieren Sie erneut den aktiven Tiefpass.

Es ist ein starkes Rauschen auf das Ausgangssignal wieder zu finden. Woher kommt dieses zustande?

1.6 Einstellung der Simulationsparameter

Unter dem Menüpunkt Simulation® Simulationsparameter, kann man diverse Einstellungen für die Simulation vornehmen. Die Start- und Stoppzeit der Simulation kann man unter dem Reiter Solver angeben. Weiterhin kann man den Solver auswählen. Matlab hat verschiedene numerische Methoden zur Lösung von Differentialgleichungen implementiert. Voreingestellt, ist ein Runge- Kutta Verfahren von Dormand und Prince. Für die meisten gängigen Simulationen ist dieses Einschrittverfahren ausreichend. Für steife Anfangswertbedingungen, welche bei Systemen auftreten deren Eigenwerte weit auseinander liegen, sollte man den impliziten Solver ode23s verwenden. Alle Verfahren haben eine automatische Schrittweitenregelung. Über einen Integrationsschritt wird die Integrationsmethode zweimal mit unterschiedlicher Fehlerordnung angewandt. Die Differenz der beiden Lösungen wird verwendet, um den Fehler des Integrationsschrittes zu schätzen. Dabei wird vorausgesetzt, dass die niederwertigen Glieder der Taylorreihenentwicklung um die Lösungstrajektorie dominant sind. Die maximale und minimale Schrittweite, kann man im Dialogfeld vorgeben, ebenso wie die Anfangsschrittweite. Der Fehler eines Integrationsschrittes wird mit Hilfe der relativen und absoluten Toleranz geregelt. Für jeden Zustand des Systems wird versucht innerhalb eines Integrationschrittes den numerischen Fehler nach folgender Formel zu begrenzen.

e(i) £ max(RelTol·abs(x(i)),AbsTol(i))

Da der Übergangsvorgang im Falle des aktiven Tiefpasses sehr schnell vonstatten geht, muss eine kleine Schrittweite gewählt werden. Wählt man die absolute Toleranz zu AbsTol=1e-12, wird das korrekte Ergebnis angezeigt. Dieses Beispiel soll aufzeigen, dass die numerische Simulation immer fehlerbehaftet ist, und Simulationsergebnisse immer kritisch begutachtet werden sollte.

Definition lokaler Fehler

1.7 Simulation eines Simulink- Modells von der Kommandozeile

Möchte man den aktiven Tiefpass mit eine Reihe unterschiedlicher Parameter simulieren, ist es zweckmäßig das Modell von der Kommandozeile aus zu simulieren. Diese Vorgehensweise eröffnet uns die Möglichkeit das Modell im Frequenzbereich zu analysieren. Um ein Simulink- Modell von der Kommandozeile aus zu simulieren, müssen die Ein- und Ausgänge des Modells nach außen geführt werden. Dies Geschieht durch die Simulink- Blöcke In und Out, welche in der Simulink Block Library unter Simulink® Sources und Simulink® Sink zu finden sind. Diese Blöcke werden zusätzlich eingefügt. Das zusätzliche Eingangssignal wird einfach durch einen Summationsblock, der sich in Simulink® Math befindet, zum vorhandenen addiert. Man erhält folgendes Blockschaltbild. Hinweis: Durch drücken der Strg- bzw. Ctrl- Taste lässt sich an einer vorhandenen Signalleitung eine zweite anschließen.

Speichern Sie das Simulink- Modell des aktiven Tiefpasses unter den Namen ATPSim.mdl ab. Wechseln Sie in die Kommandozeile von Matlab und führen sie folgenden Befehl aus.

>> [t,x,y] = sim('ATPSim');

Nachdem das Modell simuliert wurde, erhält man drei Variable. Die Variable t ist der Zeitvektor. Die Variable x ist eine Matrix in der die Spalten jeweils einem internen Zustand des Models entsprechen. Die Variable y ist der Ausgangsvektor. Das Simulationsergebnis kann man sich mit dem Befehl plot anzeigen lassen.

>> plot(t,y);

Es wird das gleiche Ergebnis, wie in Simulink angezeigt. Mit Hilfe des Kommandos grid lassen sich Gitternetzlinien ins plot- Fenster einzeichnen.

>> grid on;

Achsenbezeichnungen sowie einen Titel können ebenfalls angegeben werden mit dem folgenden Fragment.

>> title(`Simulation eines aktiven Tiefpasses');

>> xlabel(`Zeit in Sekunden');

>> ylabel('Ausfangssignal');

>> axis([ 0 tend -1.5 0.5]);

Es sei noch einmal erinnert, dass man für alle Befehle mit dem help Kommando eine kurze Erklärung bekommen kann. Probiert es bitte aus!

>> help title

>> help axis

1.8 Analyse eines Simulink- Modells von der Kommandozeile

Man kann sich auch den Frequenzgang des Simulink- Modells anzeigen lassen. Da der Frequenzgang nur für lineare Systeme Gültigkeit besitzt, muss das Simulink- Modell zuerst linearisiert werden. Auch wenn in unserem Beispiel bereits ein lineares System vorliegt, muss dieser Schritt ausgeführt werden. Die Linearisierung erfolgt durch den Befehl linmod.

>> [A,B,C,D] = linmod('ATPSim');

Dieser Befehl liefert die vier Matrizen zur Beschreibung eines Differentialgleichungssystems zurück. Jedes lineare System kann man als Differentialgleichungssystem erster Ordnung in folgender Form darstellen.

=A·x+B·u

y=C·x+D·u

Aus diesen vier Matrizen wird nun mit dem Befehl ss ein System- Objekt mit den Namen atp_sys erzeugt. Dies ist eine Matlab spezifische Datenstruktur zur allgemeinen Darstellung von Systemen.

>> atp_sys=ss(A,B,C,D);

Mit dem Befehl bode, kann man nun das Bode- Diagramm des Systems zeichnen.

>> bode(atp_sys);

Bode- Diagramm des aktiven Tiefpassfilters

1.9 Hausaufgabe

1) Es ist die nachstehende aktive Schaltung angegeben. Es soll wieder mit der Methode der Operatorimpedanzen das mathematische Modell hergeleitet werden. Ermitteln Sie den Frequenzgang der Schaltung. Um was für eine Schaltung handelt es sich?

1. aktive Schaltung

2) Es ist eine weitere aktive Schaltung gegeben. Auch hier soll wieder mit der Methode der Operatorimpedanzen das mathematische Modell hergeleitet werden. Hinweis: Die Ströme sind hier vorgegeben! Ermitteln Sie den Frequenzgang der Schaltung. Um was für eine Schaltung handelt es sich hierbei?

2. aktive Schaltung

1.10 Auflösung der Hausaufgabe

1) Die Übertragungsfunktion des Gliedes lautet:

G(s)=

U_a (s)

U_e (s)

s·R₂ ·C₁

s²·R₁·R₂ ·C₁ ·C₂ +s·C₂ ·( R₁ +R₂ )+1

Dies ergibt ein Bandpassfilter.

Bandpassfilter

2) Die Übertragungsfunktion des Gliedes lautet:

G(s)=

U_a (s)

U_e (s)

R₂ -s·R₃ ·R₁ ·C₁

R₂·( 1+s·R₁ ·C₁ )

Für R₁=R₂=R₃ ergibt es ein Allpass.

Allpassfilter

2 Seminar

2.1 Satellitenbewegung

Ein Satellit bewegt sich, auf einer kreisförmigen Umlaufbahn, um die Erde. Seine Bewegung soll mit Hilfe von Matlab simuliert werden.

Satellitenbewegung

Für die Herleitung der dynamischen Gleichungen, verwenden wir hier Newtons Grundgesetz der Mechanik.

m₂ ·

××

å
i

F_{x_i} , m₂ ·

××

å
i

F_{y_i}

In unserem Beispiel soll nur die Gravitationskraft berücksichtigt werden, die zwischen Satellit und Erde wirkt. Diese wurde von Newton wie folgt angegeben:

F=-g·

m₁ ·m₂

r²

Um das Kraftgesetz in vektorieller Form zu erhalten, muss die Kraft F mit dem Einheitsvektor multipliziert werden, der vom Anziehungszentrum m₁ zum Satteliten m₂ gerichtet ist. Als Koordinatenursprung wird der Ort der Masse m₁ festgelegt. Den Ortsvektor [r\vec]=x·e_x +y·e_y , kann man mit seinem Betrag | [r\vec] |=r=Ö{x²+y²} in die nachstehende Gleichung einsetzen und erhält die vektorielle Beschreibung in kartesischen Koordinaten.

=-g·

m₁ ·m₂

r²

ê
ê

æ
ç
ç
è

F_x

F_y

ö
÷
÷
ø

=-g·

m₁ ·m₂

( x²+y²)^[3/2]

æ
ç
ç
è

ö
÷
÷
ø

Werden diese Kräfte in ( 1) eingesetzt, erhält man ein DGL- System zweiter Ordnung.

××

=-g·

m₁

( x²+y²)^[3/2]

·x

××

=-g·

m₁

( x²+y²)^[3/2]

·y

Fast alle Software- Pakete für die Simulation unterstützen nur DGL- Systeme erster Ordnung. Da jedes DGL- System höherer Ordnung in ein System erster Ordnungen transformiert werden kann, stellt dies auch kein Problem dar. Hierzu werden folgende Variablen eingeführt.

;

××

Die endgültigen Gleichungen für die Simulation der Satellitenbewegung sind nachfolgend dargestellt.

=-g·

m₁

( x²+y²)^[3/2]

·x

=-g·

m₁

( x²+y²)^[3/2]

·y

2.1.1 M- function zur DGL- Beschreibung

Erzeugen Sie eine neue M- Datei mit den Namen satellite durch das nachfolgende Kommando.

>> edit satellite

Schreiben Sie in dieses M- File die Bewegungsgleichungen für den Satelliten, wie nachfolgend angegeben.

function dqdt=satellite( t, q)

gamma=6.672e-11; % [Nm^2/kg^2]

m1=5.99e24; % [kg]

dqdt=zeros(4,1);

dqdt(1)=-gamma*m1*q(3)/(q(3)^2+q(4)^2)^(3/2);

dqdt(2)=-gamma*m1*q(4)/(q(3)^2+q(4)^2)^(3/2);

dqdt(3)=q(1);

dqdt(4)=q(2);

Es wird hier eine Funktion satellite definiert, der die Zeit t und den Zustandsvektor q übergeben werden. Der Zustandsvektor besteht aus vier Elementen, auf die mit den Klammeroperator ( ) zugegriffen werden kann. Das M- File berechnet die Ableitungen des DGL- Systems und speichert es in den Vektor dqdt ab, welches vom M- File zurückgegeben wird.

Um das DGL- System zu lösen, braucht man einen Differentialgleichungslöser. Hier soll zuerst das Runge- Kutta- Verfahren ode45 verwendet werden. Verwenden Sie folgende Anweisungen, um die DGL's zu lösen. Sie können dafür auch eine Script- Datei schreiben.

gamma=6.672e-11; % Gravitationskonstante [Nm^2/kg^2]

m1=5.99e24; % Erdmasse [kg]

re=6370000; % Erdradius [m]

hs=200000; % Satellitenhöhe [m]

vs=sqrt(gamma*m1/(re+hs)); % Bahngeschwindigkeit [m/s]

tend=2*pi*sqrt((re+hs)^3/(gamma*m1));

% eine Periode um die Erde

options = odeset('RelTol',1e-12);

[t,q]=ode45(@satellite,[0 tend],[ 0 vs re+hs 0], options);

plot(q(:,3),q(:,4));

grid on;

title('Satellitenbewegung');

xlabel('x /[m]');

ylabel('y /[m]');

Nach der Berechnung der gewünschten Anfangsbedingungen, wird mit odeset die relative Toleranz auf eine hohe Genauigkeit eingestellt. Dem Befehl ode45 wird die Funktionsadresse der M- Function satellite übergeben. Dies geschieht mit dem Adressoperator @. Das zweite Funktionsargument ist ein Zeilenvektor, indem Start und Endzeit definiert sind. Das dritte Argument ist die Anfangsbedingung des Zustandsvektors. Er ist ebenfalls ein Zeilenvektor. Das vierte Funktionsargument ist die Struktur options. In diesem stehen alle Einstellungen für den Gleichungslöser. Diese Struktur lässt sich einfach einsehen, indem in der Matlab Kommandozeile der Strukturname options eingegeben wird.

>> options

Die mit odeset vorgenommenen Einstellungen, sind sofort sichtbar. Als Lösung erhält man einen Spaltenvektor t und eine Matrix q zurück. In t sind die Zeitpunkte der einzelnen Simulationsschritte enthalten. In q ist für jedem Zeitpunkt eine Zeile mit den dazugehörigen Zustandswerten enthalten. Da nur die Position des Satelliten von Interesse ist, wird dem plot- Kommando nur die dritte und vierte Spalte von q übergeben. Der Doppelpunkt gibt hier an, dass alle Elemente der Zeilen verwendet werden sollen. Die Lösung der DGL's ist untenstehend dargestellt.

Lösung

2.1.2 Simulationsfehler

Wiederholen Sie die Simulation mit der 10-fachen Simulationszeit. Das bedeutet, dass sich der Satellit nun zehnmal um die Erde dreht.

% 10- Perioden um die Erde

options = odeset('RelTol',1e-12);

[t,q]=ode45(@satellite,[0 10*tend],[ 0 vs re+hs 0], options);

plot(q(:,3),q(:,4));

grid on;

title('Satellitenbewegung');

xlabel('x /[m]');

ylabel('y /[m]');

Wie Sie aus der Lösung entnehmen können, bewegt sich der Satellit immer auf der Kreisbahn. Wiederholen Sie die Simulation mit einer geringeren Toleranz, wie unten angegeben. Um das vorhergehende Plot- Fenster nicht zu überschreiben, öffnet man mit dem Befehl figure ein neues.

% 10- Perioden um die Erde mit geringer Toleranz

options = odeset('RelTol',1e-3);

[t,q]=ode45(@satellite,[0 10*tend],[ 0 vs re+hs 0], options);

figure;

plot(q(:,3),q(:,4));

grid on;

title('Satellitenbewegung');

xlabel('x /[m]');

ylabel('y /[m]');

Es ist deutlich der Fehler des Simulationsergebnisses zu ersehen. Mit jedem Simulationsschritt wird ein Fehler begangen, welcher sich aufsummiert. Die Steuerung des Fehlers in einem Schritt kann durch die RelTol und die AbsTol vorgegeben werden.

2.2 Simulation eines RLC- Gliedes

Nachstehend ist der Schaltplan eines RLC- Gliedes abgebildet. Die Herleitung der Modellgleichungen soll nun im Zeitbereich erfolgen.

RLC- Glied

Die Spannungsabfälle in der Masche lassen sich sofort aus dem Bild entnehmen.

R·i(t)+L·

di(t)

t
ó
õ
0

i(t) dt=U₀

Um das Integral zu eliminieren, wird die gesamte Gleichung differenziert. Die Spannung U ist dabei als konstante Betriebsspannung aufzufassen.

R·

di(t)

+L·

d²i(t)

dt²

·i(t)=0

Die Differentialgleichung des Stromes kann man schließlich in die folgende Standardform bringen.

R·C·

di(t)

+L·C·

d²i(t)

dt²

+i(t)=0

Die Eigenfrequenz des RLC- Gliedes ist w_n = [1/(Ö{L·C})] und die Dämpfung D=[1/2]·R·Ö{[C/L]} .

Die Differentialgleichung zweiter Ordnung muss wieder in zwei Differentialgleichungen erster Ordnung transformiert werden. Es werden die üblichen Phasenvariablen eingeführt.

x₁ = i

x₂ =

Die endgültigen Gleichungen für die Simulation des RLC- Gliedes sind nachfolgend dargestellt.

=x₂

L·C

·x₁ -

·x₂

2.2.1 M- function zur DGL- Beschreibung

Erzeugen Sie eine neue M- Datei mit den Namen rlc.

>> edit rlc

Schreiben Sie in dieses M- File die Differentialgleichungen für das RLC- Glied, wie nachfolgend angegeben.

function dxdt=rlc(t,x)

R=100; % Widerstand [Ohm]

L=10e-3; % Induktivität [H]

C=100e-9; % Kapazität [F]

dxdt=zeros(2,1);

dxdt(1)=x(2);

dxdt(2)=-1/(L*C)*x(1)-R/L*x(2);

Simulieren Sie das RLC- Glied mit einem Anfangsstrom von ein Ampere. Benutzen Sie folgende Zeilen dazu.

% simulation eines RLC- Gliedes

tend=1e-3;

[t,x]=ode45(@rlc,[0 tend],[ 1 0]);

plot(t,x(:,1));

grid on;

title('Stromverlauf am RLC- Glied');

xlabel('t /[s]');

ylabel('i /[A]');

Sie sollten folgenden Graph erhalten. Es ist deutlich die Dämpfung der Schwingung zu ersehen.

Lösung

2.2.2 Parameterübergabe

Angenommen Sie wollten für verschiedene Parameter für R, L und C die Lösung numerisch berechnen lassen, dann müssten Sie immer die M- function rlc abändern. Um dieses zu vermeiden, kann man die Werte der Bauelemente als zusätzliche Parameter übergeben. Ändern Sie die M- function, wie nachfolgend dargestellt.

function dxdt=rlc(t,x,R,L,C)

dxdt=zeros(2,1);

dxdt(1)=x(2);

dxdt(2)=-1/(L*C)*x(1)-R/L*x(2);

Alle zusätzlich vereinbarte Parameter müssen nun in der gleichen Reihenfolge beim Aufruf von ode45 übergeben werden. Da die Struktur options nicht verwendet werden soll, wird für sie eine leere Matrix als Platzhalter übergeben.

% simulation eines RLC- Gliedes

R=100; % Widerstand [Ohm]

L=10e-3; % Induktivität [H]

C=100e-9; % Kapazität [F]

tend=1e-3;

[t,x]=ode45(@rlc,[0 tend],[ 1 0],[],R,L,C);

plot(t,x(:,1));

grid on;

title('Stromverlauf am RLC- Glied');

xlabel('t /[s]');

ylabel('i /[A]');

Wenn nun die Temperaturabhängigkeit des Widerstandes mit dem Einfluss auf die Lösung untersucht werden soll, kann das einfach durch folgende paar Zeilen geschehen.

% simulation eines RLC- Gliedes

L=10e-3; % Induktivität [H]

C=100e-9; % Kapazität [F]

tend=1e-3;

R=95; % Widerstand [Ohm]

[t1,x1]=ode45(@rlc,[0 tend],[ 1 0],[],R,L,C);

R=100; % Widerstand [Ohm]

[t2,x2]=ode45(@rlc,[0 tend],[ 1 0],[],R,L,C);

R=105; % Widerstand [Ohm]

[t3,x3]=ode45(@rlc,[0 tend],[ 1 0],[],R,L,C);

plot(t1,x1(:,1),':b',t2,x2(:,1),'-r',t3,x3(:,1),'-m');

grid on;

title('Stromverlauf am RLC- Glied');

xlabel('t /[s]');

ylabel('i /[A]');

legend('R=95 \Omega','R=100 \Omega','R=105 \Omega');

Dem plot- Kommando wurden nun alle drei Lösungen übergeben. Um die Lösungen besser unterscheiden zu können, werden diese jeweils andersfarbig dargestellt. Auch die Liniendarstellung ist unterschiedlich. Geben Sie den Hilfe- Befehl für dieses Kommando ein, damit Sie alle Möglichkeiten einsehen können.

>> help plot

Zusätzlich wurde auch eine Legende angelegt, in der der Widerstand der Lösung zugeordnet ist. Dabei wurde für den griechischen Buchstaben der entsprechende Latex- Befehl angegeben.

Ist man nur an einem kleinen Ausschnitt der Lösung interessiert, dann lässt sich diese mit dem axis- Befehl anzeigen.

>> axis([0.2e-3 0.4e-3 -0.25 0.4]);

Berechnen Sie nun die Dämpfung und die Eigenfrequenz in der Kommandozeile.

>> f=1/(2*pi*sqrt(L*C))

>> D=0.5*R*sqrt(C/L)

2.3 Nichtlinearer Oszillator

In der frühen Rundfunk- und Fernsehtechnik war es notwendig Oszillatorschaltungen zu entwickeln, die auf einer festen Frequenz stabile Schwingungen ausführen. Der RLC- Oszillator musste dabei so erweitert werden, dass seine gedämpften Schwingungen nicht auslöschen. Die Schlüsselidee war es, ein Oszillator zu konstruieren der für kleine Strombeträge instabil ( aufschwingend) und für große Strombeträge stabil ( gedämpft) sein soll. Der passive Widerstand des RLC- Oszillators wird ersetzt durch ein aktives Element. Balthazar van der Pol, ein Holländischer Ingenieur, benutzte in den zwanziger Jahren des 20. Jahrhunderts dafür mehrere Trioden- Röhren. Heute würde eine Tunnel- Diode ausreichen. Die nichtlineare Charakteristik dieses Elements besitzt einen negativen Widerstand in einem Bereich nahe null und einen positiven Widerstand bei höheren Strömen. Dadurch wird erreicht, dass immer wieder neue Energie in die Schaltung zugeführt wird, sobald die Schwingungen vom Betrag kleiner werden.

Van der Pol'sche Oszillator

Eine solche Schaltung ist obenstehend dargestellt. Ähnlich, wie beim RLC- Glied werden die Spannungsabfälle in der Masche dem Bild entnommen.

f( i(t) )+L·

di(t)

t
ó
õ
0

i(t) dt=U₀

Auch hier wird die gesamte Gleichung differenziert, um das Integral zu eliminieren.

f¢( i(t) )·

di(t)

+L·

d²i(t)

dt²

·i(t)=0

Die nichtlineare Funktion des aktiven Elements, kann man in ein Polynom dritten Grades approximieren. Van der Pol benutzte dabei folgenden Ansatz:

f( i )=-a·i+b·i³.

Damit vereinfacht sich die DGL zu:

( -a+3·b·i² )·

di(t)

+L·

d²i(t)

dt²

·i(t)=0.

Um allgemeine Aussagen über das nichtlineare Verhalten des Oszillators zu machen, ist es zweckmäßig die Differentialgleichung zu normieren. Hierfür wird eine neue dimensionslose Zeitbasis t = w·t eingeführt, sowie eine noch festzulegende lineare Transformation der Stromstärkex=¹/_p·i. Damit folgt der Übergang zu den neuen Variablen mit

i(t)=p·x(t(t)),

di(t)

=p·

dx(t)

=p·w·

dx(t)

d²i(t)

dt²

=p·

d²x(t)

dt²

æ
è

ö
ø

+p·

dx(t)

d²t

dt²

=p·w²·

d²x(t)

dt²

zu:

3·b·p²·

æ
Ö

æ
è

3·b·p²

+x²

ö
ø

dx(t)

d²x(t)

dt²

+x(t)=0.

Die Wahl von p=Ö{[a/(3·b)]} vereinfacht die Gleichung zu:

d²x(t)

dt²

+a·

æ
Ö

( x²-1)·

dx(t)

+x(t)=0.

Wird nun auch noch der Koeffizient zu m = a·Ö{[C/L]} gewählt, erhält man die Standardform der van der Pol'schen Differentialgleichung.

x"+m·( x²-1 )·x¢+x=0.

2.3.1 Darstellung der nichtlinearen Kennlinie

Als erstes soll die nichtlineare Kennlinie des aktiven Elementes dargestellt werden. Dafür wird ein Vektor mit 101 Elemente erzeugt, deren Werte von -0.5 bis +0.5 gleich verteilt sind. Anschließend werden für diese Argumente die Funktionswerte berechnet. Benutzen Sie folgende Anweisungen.

a=2.2;

b=18.1;

i=linspace(-0.5,0.5, 101);

u=-a*i+b*i.^3;

plot(i,u);

title('nichtlineare Kennlinie des aktiven Elementes');

xlabel('i/ [A]');

ylabel('u/ [V]');

grid on;

Damit erzeugen Sie den nachfolgenden Graphen. Wenn eines der Befehle unklar ist, rufen Sie die dazugehörige Kommandozeilen- Hilfe auf!

z.B.

>> help linspace

Nichtlineares Element

2.3.2 Simulation des van der Pol'schen Oszillators

Nachdem die normierten Gleichungen des van der Pol'schen Oszillators hergeleitet wurden, muss diese wieder in ein DGL System erster Ordnung überführt werden. Durch Einführung der Phasenvariablen:

x₁ = x,

x₂ = x¢,

erhält man die für die Simulation geeignete mathematische Beschreibung,

x₁ ¢=x₂ ,

x₂ ¢=-m·( x₁ ²-1 )·x₂ -x₁ .

und das korrespondierende M- File,

function dxdt=Pol( t, x, mu)

dxdt=zeros(2,1);

dxdt(1)=x(2);

dxdt(2)=-mu*(x(1)^2-1)*x(2)-x(1);

Für die Simulation des DGL- Systems, können Sie folgende Kommandos benutzen.

% Simulation des van der Pol'schen Oszillators

mu=1; % Parameter

tend=40;

[t,x]=ode45(@pol,[0 tend],[ 2 0],[],mu);

plot(t,x(:,1));

grid on;

title('normierter Stromverlauf des van der Pol`schen Oszillators');

xlabel('\tau');

ylabel('x_1');

Signalverlauf am van der Pol'schen Oszillator

2.3.3 Phasenebene des van der Pol'schen Oszillators

Ähnlich wie bei der Simulation der Planetenbewegung ist das periodische Verhalten optimal in der Phasenebene zu ersehen. Folgende Anweisungen machen diese sichtbar.

plot(x(:,2),x(:,1));

grid on;

title('Phasenebene des van der Pol`schen Oszillators');

xlabel('x_2');

ylabel('x_1');

Phasenebene des van der Pol'schen Oszillator

Benutzen Sie nun unterschiedliche Anfangsbedingungen, um zu zeigen, dass alle in eine stabile periodische Schwingung münden. Wiederholen Sie alle vorhergehenden Anweisungen und ändern sie folgende Zeile ab.

[t,x]=ode45(@pol,[0 tend],[ 0.1 0],[],mu);

bzw.

[t,x]=ode45(@pol,[0 tend],[ 10 0],[],mu); ( Hinweis: Vergrößern Sie die Endzeit, tend.)

x=[0.1 0]

x=[10 0]

2.3.4 Simulation mit m = 0

Setzen Sie nun den Parameter m zu null, und starten Sie wieder die Simulation des van der Pol'schen Oszillators mit der Anfangsbedingung x=[2 0]. Was stellen Sie fest? Schauen Sie sich die Trajektorien auch im Zeitbereich an. Warum erhalten Sie dieses Ergebnis? Analysieren Sie dazu auch die Differentialgleichung.

2.3.5 Simulation mit unterschiedlichen m

Vergleichen Sie nun die Simulationsergebnisse des van der Pol'schen Oszillators bei unterschiedlichen m mit der Anfangsbedingung x=[0.01 0]! Erhöhen Sie dabei schrittweise m von 0.1 nach 5. Zwei solcher Ergebnisse sind untenstehend dargestellt. Welche Schlussfolgerung ziehen Sie, wenn m größer als zwei wird?

x=[0.01 0]; m=0.1

x=[0.01 0], m=5

2.3.6 Simulation mit m=100

Wiederholen Sie die Simulation mit m=100, tend=400 und x=[2 0]. Was fällt Ihnen auf?

x=[2 0]; m=100

Es gibt Bereiche in der Phasenebene, wo die Anfangswertaufgabe steif wird. Das bedeutet, dass die Schrittweite des Standard Runge- Kutta Gleichungslösers hier stark verringert werden muss, um Stabilität zu gewährleisten. Das klassische Runge- Kutta- Verfahren ist ein explizites Verfahren. Für steife Anfangswertaufgaben sind implizite Verfahren zu bevorzugen. Wiederholen Sie die Simulation mit dem modifizierten Rosenbrock- Verfahren ode23s.

mu=100; % Parameter

tend=400;

[t,x]=ode23s(@pol,[0 tend],[ 2 0],[],mu);

plot(x(:,2),x(:,1));

grid on;

title('Phasenebene des van der Pol`schen Oszillators');

xlabel('x_2');

ylabel('x_1');

Man erkennt deutlich den Rechengeschwindigkeitsvorteil. Nun können Sie sich auch die Simulationsergebnisse für größere m, z.B. m =1000, anzeigen lassen.

2.3.7 Theoretische Analyse

Zur Vollständigkeit soll hier noch eine kurze theoretische Analyse der DGL gemacht werden. Die möglichen Ruhepunkte des Systems erhält man, indem die Ableitungen des DGL- Systems zu null gesetzt werden.

0=x₂

0=-m·( x₁ ²-1 )·x₂ -x₁

Man erkennt sofort, dass nur ein Ruhepunkt vorhanden ist, nämlich x₁=0 und x₂=0. Um Aussagen zur Stabilität der Ruhelage zu machen, muss die Jacobimatrix des Systems berechnet werden. Im Nullpunkt ist diese:

J(0,0)=

æ
ç
ç
è

-1

ö
÷
÷
ø

Da die Determinante verschieden von null ist, ist das lineare Glied der Taylorreihe dominant und eine Linearisierung in einem hinreichend engen Bereich um den Ruhepunkt erlaubt. Die Eigenwerte der linearen Systemmatrix ergeben sich zu:

ê
ê
ê
ê

-1

l-m

ê
ê
ê
ê

=l²-l·m+1,

l_1,2 = +

æ
Ö

m²

-1

Für positive m kleiner als zwei, ergibt dies ein konjugiert komplexes Eigenwertpaar in der rechten Halbebene. Es entsteht eine spiralförmig aufstrebende Schwingung, welche durch die Simulation nachgewiesen werden konnte. Für m größer gleich zwei, entstehen immer zwei reelle Pole in der rechten Halbebene. Ist m gleich null, entsteht ein konjugiert komplexes Eigenwertpaar auf der imaginären Achse, welche eine harmonische Schwingung im Zeitbereich verursacht.

3 Seminar

3.1 Simulation eines Gleichstrommotors

Mittels Simulation soll das dynamische Verhalten eines Gleichstrommotors einschließlich der mit dieser gekoppelten Arbeitsmaschine untersucht werden. Das System wird in der Umgebung des Arbeitspunktes durch das nachfolgende Ersatzschaltbild beschrieben. Dabei sind die Werte der Variablen die Abweichungen von den Werten im Arbeitspunkt.

Ersatzschaltbild eines Gleichstrommotors

Parameter: R_AAnkerwiderstand L_AAnkerinduktivität J Trägheitsmoment Motor und Arbeitsmaschine k_M·F_eMotorkonstante M_R= b·w_M Reibmoment b Viskositätskonstante	Eingangsgrößen: U_AAnkerspannung M_LLastmoment Ausgangsgrößen: n_MMotordrehzahl i_A Ankerstrom
sonstige Größen: U_iinduzierte Spannung ( erzeugt durch den Kommutator) w_MWinkelgeschwindigkeit der Motorwelle
Technische Daten des Gleichstrommotors Typ: GFQUJG 080-22 P_mech =4,6 kW n =2500 min^-1 U_A =420 V I_A =13,3 A U_E =360 V I_E =0,8 A

3.1.1 Herleitung der Differentialgleichung

Zur Herleitung der DGL soll hier die Bilanzmethode verwendet werden. Diese umfassende Methode stützt sich auf den Energiegehalt des Systems.

Es gibt zwei dominierende Speicher:

elektromagnetische Energie in der Ankerinduktivität,
kinetische Energie als Rotationsenergie des Rotors und Arbeitsmaschine.

Für jeden dieser Speicher muss eine Bilanzgleichung aufgestellt werden.

Bilanz der elektromagnetischen Energie:

é
ë

·L_A ·i_A ( t)²

ù
û

gespeicherte elektro-

magnetische Energie

U_A ·i_A ( t)

zugeführte

Leistung

U_i ·i_A ( t)

abgegebene Leistung

des Motors

R_A ·i_A ( t)²

Abgegebene Leistung

in Wärme der Windungen

L_A ·[d/dt]i_A ( t )=U_A -U_i -R_A ·i_A ( t)

mit U_i = k_M ·f·w_mech und w_mech = 2·p·n.

Bilanz der kinetischen Energie:

é
ë

·J·w_mech( t )²

ù
û

gespeicherte Rotationsenergie

der Motorwelle und Arbeits-

maschine

M_Motor ·w_mech ( t)

zugeführte mechanische

Leistung des Motors

M_Reib ·w_mech ( t)

abgegebene Reib-

leistung des Motors

M_Last ·w_mech ( t)

abgegebene mechanische

Leistung des Motors

J·[d/dt]w_mech ( t )=M_Motor -M_Reib-M_Last

mit M_Motor = k_M ·f·i_A ( t ) und M_Reib = b·w_mech ( t )

Damit sind schon die zwei Differentialgleichungen des Motors gefunden. Diese beiden Differentialgleichungen werden nun in die Zustandsraumdarstellung überführt. Als Zustandsgrößen werden der Ankerstrom und die Winkelgeschwindigkeit der Motorwelle verwendet. Eingangsgrößen sollen die Ankerspannung sowie das Lastmoment sein. Für die Ausgangsgrößen werden hier Ankerstrom und Motordrehzahl festgelegt. Damit ergibt sich folgendes Modell um den Arbeitspunkt.

æ
ç
ç
è

Di_A ( t )

Dw_mech ( t )

ö
÷
÷
ø

é
ê
ê
ë

-[(R_A )/(L_A )]

-[(k_M ·f_e )/(L_A )]

[(k_M ·f_e )/J]

-[b/J]

ù
ú
ú
û

æ
ç
ç
è

Di_A ( t )

Dw_mech ( t )

ö
÷
÷
ø

é
ê
ê
ë

[1/(L_A )]

-[1/J]

ù
ú
ú
û

æ
ç
ç
è

DU_A ( t )

DM_Last ( t )

ö
÷
÷
ø

æ
ç
ç
è

Di_A ( t )

Dn( t )

ö
÷
÷
ø

é
ê
ê
ë

[60/(2·p)]

ù
ú
ú
û

æ
ç
ç
è

Di_A ( t )

Dw_mech ( t )

ö
÷
÷
ø

x=A·x+B·u

y=C·x

3.1.2 Statische Beziehungen

Eine gute Herangehensweise ist es, aus der Differentialgleichung die statischen Beziehungen zwischen Ein- und Ausgangsgrößen herzuleiten. Diese Beziehungen sind nützlich für die Identifikation von Parametern, sowie um Aussagen über das allgemeine Verhalten des Systems zu treffen. Die statischen Beziehungen erhält man, indem die Ableitungen der Zustandsgrößen zu null gesetzt werden.

DU_A ( t )=R_A ·Di_A ( t)+

k_M ·f_e ·

2·p

·Dn( t )

DU_i

DM_Last ( t )=k_M ·f_e ·Di_A (t )-b·

2·p

·Dn( t )

Man erkennt sofort, dass man durch Messung der statischen Kennlinien einige der Parameter berechnen kann. Da in der Praxis Messwerte rar sind, soll versucht werden, möglichst viele Parameter aus den Kenndaten des Typenschildes zu berechnen. Dabei wollen wir annehmen, dass die Gleichungen ( 9) global gilt, und nicht nur im Arbeitspunkt. Wird die erste Gleichung mit dem Ankerstrom durchmultipliziert und die zweite mit der Winkelgeschwindigkeit, erhält man zwei Leistungsgleichungen.

U_A ( t )·I_A ( t )

P_e

R_A ·I_A ( t )²

P_v

k_A ·f_e ·

2·p

·n( t)

U_i

·I_A ( t )

P_mech

M_Last ( t )·

2·p

·n( t )

P_mech

k_M ·f_e ·

2·p

·n( t )

U_i

·I_A ( t)-b·

æ
è

2·p

·n( t )

ö
ø

Aus diesen Gleichungen lassen sich nun einige der gesuchten Parameter aus den Kenndaten des Typenschildes berechnen.

3.1.3 Berechnung der Kennwerte

Schreiben Sie eine Skript- Datei indem Sie die Kennwerte des Gleichstrommotors definieren und die fehlenden Parameter berechnen. Da die Viskosität b unbekannt ist, wird sie vorläufig zu null gesetzt. Sie können folgende Vorlage benutzen.

% Kennwerte des Gleichstrommotors aus Typenschild

Pm=4600; % mechanische Nennleistung in [W]

Ua=420; % Ankerspannung in [V]

Ia=13.3; % Ankerstrom bei Nennleistung in [A]

Ie=0.8; % Erregerstrom in [A]

Ue=360; % Erregerspannung in [V]

nn=2500; % Nenndrehzahl in [1/min]

% Berechnung der Kennwerte und Parameter

Mn=60*Pm/(2*pi*nn); % Berechnung des Nennmomentes in [Nm]

Ui=Pm/Ia; % Berechnung der induzierten Spannung in [V]

Ra=(Ua-Ui)/Ia; % Berechnung des Ankerwiderstandes in [Ohm]

KPhi=60*Ui/(2*pi*nn); % Berechnung der Motorkonstante in [Vs]

Eine Berechnung des Trägheitsmomentes J sowie der Ankerinduktivität L_A, kann aus den statischen Beziehungen nicht gewonnen werden. Diese Größen sind nur aus der Dynamik errechenbar. Verwenden Sie die folgenden Werte und fügen Sie diese der Skript- Datei zu.

% Angenommene Werte

La=0.0402; % Ankerinduktivität [H]

J=15; % Trägheitsmoment Motor [kg*m^2]

b=0; % Viskosität [Nms]

3.1.4 Die Drehzahl- Drehmoment- Kennlinie

Aus der statischen Beziehung ( 9) lässt sich auch die ideelle Drehzahl- Drehmomenten- Gleichung herleiten, indem die Gleichung nach der Drehzahl umgestellt wird. Auch hier wird angenommen, dass die statische Beziehung global gilt sowie die Viskosität null ist.

n( t )=

U_A

2·p·k_M ·f_e

R_A

2·p·( k_M ·f_e )²

·M( t )

n( t )=n₀ -m·M( t )

Lassen Sie sich die ideelle Drehzahl- Drehmoment- Kennlinie mit Hilfe von Matlab anzeigen.

% Drehzahl- Drehmoment- Kennlinie

M=linspace(0,Mn,11);

n=60*Ua/(2*pi*KPhi) - 60*Ra/(2*pi*(KPhi)^2)*M;

figure;

plot(M,n);

title('Drehzahl- Drehmoment- Kennlinie');

xlabel('M / [Nm]');

ylabel('n / [1/min]');

axis([0 20 0 3500]);

grid on;

An diesem Gleichstrommotor wurde durch ein Experiment die Drehzahl- Drehmoment- Kennlinie aufgenommen. Ihre Werte sind untenstehend angegeben. Fügen Sie diese der ideellen Drehzahl- Drehmoment- Kennlinie zu. Dies kann durch das Kommando hold on geschehen. Lassen Sie sich die Matlab- Hilfe zu diesem Befehl anzeigen!

% Gemessene Drehzahl- Drehmomentkennlinie

Mm=[0 5 10 15]';

nm=[2710 2660 2620 2630 ]';

hold on;

plot(Mm,nm,'-r');

hold off;

legend('ideelle n( M)', 'reale n( M)');

Zwischen der ideellen und der realen Kennlinie ist eine deutliche Abweichung zu erkennen. Durch das Vernachlässigen der Viskosität und der Annahme eines linearen Modells im gesamten Betriebsbereich, ist dieses auch zu erwarten. Trotz dieser Modellierungsfehler spiegelt das Modell das wesentliche Verhalten wieder. Mit zunehmendem Lastmoment sinkt die Drehzahl.

Ideelle und reale Drehzahl- Drehmoment- Kennlinie

3.1.5 Der Signalflussplan in Simulink

Um das dynamische Verhalten zu untersuchen, muss der Signalflussplan entwickelt werden. Zum Zeichnen des Signalflussplanes geht man von den einzelnen Integratoren der Zustandsgleichung ( 7) aus. Die Zustandsgrößen sind die Ausgänge der Integratoren. Vor jeden Eingang eines Integrators wird eine Mischstelle gesetzt. In ihr fließen unter Beachtung des Vorzeichens alle Terme die auf der rechten Seite der DGL stehen. Dabei werden die Eingangsgrößen sowie die Zustandsgrößen mit ihren Koeffizienten wieder benötigt. Es ergibt sich der folgende Signalflussplan in Simulink.

Signalflussplan Gleichstrommotor

Öffnen Sie ein neues Simulink- Blockschaltbild, indem Sie den Simulinkbrowser öffnen und File® New® Model wählen. Speichern Sie das Modell unter motor.mdl ab. Holen Sie sich aus der Block-Library die notwendigen Blöcke, um das Signalflussbild zu zeichnen. Die Eingänge sind InPort- Blöcke aus Sources und die Ausgänge sind OutPort- Blöcke aus Sinks. Beachten Sie die Numerierung!

Durch Aufruf des Kontextmenüs eines Blockes mit Hilfe der rechten Maustaste, lässt sich unter Format® Rotate die Anweisung für das Drehen des Blockes finden. Mit Format® Hide Name, kann der Blockname ausgeblendet werden. Ein Knotenpunkt bzw. Signalabzweig bekommt man, indem an der entsprechenden Stelle die linke Maustaste und gleichzeitig auch die Strg- Taste gedrückt wird. Die Parameter der Differentialgleichung sollen alle durch Variable definiert werden.

3.1.6 Blockmaskierung in Simulink

Um komplexe Systeme in kleine Systeme aufzuspalten, ist es unerlässlich, Teilsysteme als ein Block zusammenzufassen. Hier soll nun der komplette Motor als ein Block dargestellt werden. Gehen Sie auf den Menüpunkt Edit® Select all. Damit sind nun alle Blöcke und Signalleitungen ausgewählt und als markiert angezeigt. Gehen Sie wieder ins Menü und wählen Sie Edit® Create Subsystem. Dadurch entsteht ein einzelner Block mit zwei Ein- und Ausgängen. Durch Doppelklick auf den Blocknamen, können Sie den Blocknamen auf Motor umändern. Die Größe des Blockes können Sie verändern, durch das Anwählen des Blockes mit einem Mausklick, und das Ziehen der Maus von einer Ecke des Blockes. Durch Doppelklick auf den Block, öffnet sich ein neues Fenster, indem der Inhalt des Blockes angezeigt wird. Schließen Sie dieses Fenster wieder.

Es ist wünschenswert, die Parameter des Blockes in einen Blockdialog anzuzeigen, wie man es von den Standard Simulink- Blöcken gewohnt ist. Dies wird erreicht, indem mit der rechten Maustaste das Kontextmenü, für den soeben erstellten Block des Gleichstrommotors, aufgerufen wird. Wählen Sie Mask subsystem. Es öffnet sich eine mehrreitige Dialogbox. Wählen Sie mit der Maus den zweiten Reiter ,,Initialization". Betätigen Sie die Schaltfläche Add. Geben Sie nach Prompt: ,,Ankerwiderstand:" ein. Bei Variable: geben Sie ,,Ra" ein. Betätigen Sie erneut die Add Schaltfläche. Achten Sie dabei, dass der blaue Balken auf < < end of parameter list > > steht. Geben Sie nun unter Prompt: ,,Ankerinduktivität:" ein. Bei Variable geben Sie ,,La" ein. Betätigen Sie wieder die Add Schaltfläche. Bei Prompt: geben Sie bitte ,,Motorkonstante:" und bei Variable ,,KPhi" ein. Drücken Sie erneut die Add Taste, um anschließend bei Prompt: ,,Trägheitsmoment:" und bei Variable ,,J" einzugeben. Wiederholen Sie das ganze für die Viskosität, indem Sie die Add Taste betätigen und bei Prompt: ,,Viskosität" und bei Variable ,,b" eingeben. Schließen Sie das Dialogfeld mit OK.

Erneutes Doppelklicken auf den Block lässt den folgenden Parameterdialog erscheinen.

Blockschaltbild Motor

Dazugehöriger Parameterdialog

3.1.7 Bode- Diagramm

Ähnlich wie im Seminar 1, lassen sich durch den Befehl linmod die Systemmatrizen bestimmen. Anschließend kann das Bode- Diagramm angezeigt werden. Laden Sie die Parameter des Gleichstrommotors und erzeugen Sie danach das Bode- Diagramm.

>> [A,B,C,D]=linmod('motor');

>> figure;

>> bode(A,B,C,D);

>> grid on;

Bode- Diagramm des Gleichstrommotors

Versuchen Sie mit Hilfe des Blockschaltbildes die Übertragungsglieder zu deuten!

3.1.8 Arbeitspunkteinstellung

Angenommen der Gleichstrommotor treibt eine Arbeitsmaschine an, die eine Drehzahl von n0 = 2500 u/min verlangt. Das erforderliche Drehmoment ist im Normalfall M0=10. Eine Drehzahlregelung soll entwickelt werden, die auch bei Abweichungen des Drehmomentes, die Drehzahl konstant hält. Anhand der statischen Gleichungen, muss man zuerst den notwendigen Ankerstrom sowie die notwendige Ankerspannung berechnen. Fügen Sie zur Parameterdatei die Arbeitspunktvariablen zu.

% Werte im Arbeitspunkt

M0=10; % Lastmoment [Nm]

n0=2500; % Drehzahl [1/min]

Ia0=M0/KPhi; % Ankerstrom [A]

Ua0=Ra*Ia0+KPhi*2*pi/60*n0; % Ankerspannung [V]

Es ist nun wünschenswert mit den absoluten Werten anstatt mit den Abweichungen zu arbeiten. Hierfür müssen nun die Arbeitspunktwerte unter Beachtung des Vorzeichens zu den Abweichungen addiert werden. Ändern Sie das Modell wie folgt ab.

Simulink- Modell

3.1.9 Sprungantworten

Es soll nun untersucht werden, wie sich der Gleichstrommotor bei einer sprungförmigen Änderung des Lastmomentes verhält. Erzeugen Sie nach fünf Sekunden einen Sprung im Lastmoment von M0 nach M0+0.1*M0. Dazu müssen Sie in dem Block für den Lastmomentsprung die entsprechenden Parameter einstellen.

Setzen Sie die beiden Sprungamplituden für den Ankerspannungssprung auf Ua0. Stellen Sie die Simulationszeit auf tend=200s, und starten Sie die Simulation. Sie sollten folgendes Ergebnis erhalten.

Simulationsergebnis: Ankerstrom

Simulationsergebnis: Drehzahl

Da eine höhere mechanische Leistung gefordert wird, muss sich der Ankerstrom erhöhen. Die Drehzahl bleibt jedoch nicht konstant und fällt langsam auf einen geringeren Wert.

Setzen Sie nun beide Sprunghöhen des Lastmomentes auf M0, und erzeugen Sie nach fünf Sekunden einen Sprung von Ua0 nach Ua0+0.05*Ua0. Starten Sie die Simulation erneut.

Simulationsergebnis: Ankerstrom

Simulationsergebnis: Drehzahl

Der Ankerstrom erhöht sich fast sprunghaft, und die damit erhöhte Antriebsenergie wird zur Beschleunigung des Ankers verwendet. Die Drehzahl steigt, während der Ankerstrom wieder sinkt.

3.1.10 Drehzahlregelung mit P- Regler

Damit die Drehzahl nicht absinkt, wenn das Lastmoment steigt, soll diese geregelt werden. Hierzu wird die Drehzahl mit der Solldrehzahl verglichen. Die Differenz wird verstärkt und als Ankerspannung zurückgeführt.

Ändern Sie das Blockschaltbild, wie nachfolgend dargestellt, ab. Die Abweichung der Solldrehzahl vom Arbeitspunkt ist selbstverständlich null. Als Verstärkungsfaktor des P- Reglers soll fünfzig gewählt werden.

Drehzahlregelung mit P- Regler

Stellen Sie einen Lastsprung bei 0,5 s von M0 nach M0*(1+0.1) im Parameterdialog des Sprungblockes ein. Setzen Sie die Simulationsendzeit auf 2 s und starten Sie die Simulation.

Simulationsergebnis: Ankerspannung

Simulationsergebnis: Drehzahl

Die Drehzahl wird nun auf 2500 u/min geregelt. Dazu erhöhte der P- Regler die Ankerspannung. Da für die Erhöhung der Ankerspannung eine bleibende Regelabweichung notwendig ist, fällt die Drehzahl geringfügig.

3.1.11 Drehzahlregelung mit PI- Regler

Um eine Regelabweichung von null zu bekommen, kann ein PI- Regler verwendet werden. Holen Sie aus der Simulink- Block- Library aus Simulink Extras ® Additional Linear den Block PID- Controller, und ersetzen Sie damit den P- Regler.

Drehzahlregelung mit PI- Regler

Stellen Sie einen Proportionalwert von fünfzig und ein Integralwert von fünfzig ein. Erhöhen Sie die Simulationszeit auf 10 s. Eine erneute Simulation sollte folgendes Ergebnis liefern.

Simulationsergebnis: Ankerspannung

Simulationsergebnis: Drehzahl

Die Regelabweichung in der Drehzahl verschwindet nun.

3.1.12 Unterlagerte Stromregelung

Für eine schnelle Ausregelung einer Störung im Lastmoment, ist es wichtig diese Störgröße frühzeitig messtechnisch zu erfassen. Aus der statischen Beziehung ( 9) ist zu entnehmen, dass der Strom linear vom Lastmoment abhängt. Die messtechnische Erfassung des Lastmomentes ist wesentlich schwieriger als der des Ankerstromes. Deswegen soll ein einfacher P- Regler als unterlagerte Stromregelung verwendet werden. Eine solche Anordnung nennt man Kaskadenregelung.

Verändern Sie das Blockschaltbild, wie untenstehend, ab. Der P- Regler soll den Verstärkungswert 500 besitzen. Starten Sie erneut die Simulation.

Drehzahlregelung durch Kaskade von PI- und P- Regler

Das Simulationsergebnis ist unten abgebildet.

Simulationsergebnis: Ankerspannung

Simulationsergebnis: Drehzahl

Da die Ankerspannung frühzeitiger durch die unterlagerte Stromregelung nachgeregelt wird, ist die maximale Regelabweichung in der Drehzahl kleiner.

3.2 Gleichstrommotor mit Erregerfeld

Der Gleichstrommotor soll nun mit Erregerfeld simuliert werden. Das Ersatzschaltbild wird erweitert durch die Erregerwicklung. Der magnetische Fluss soll nicht mehr als konstant angenommen werden.

Ersatzschaltbild eines Gleichstrommotors mit Erregerfeld

Parameter: R_AAnkerwiderstand L_AAnkerinduktivität J Trägheitsmoment Motor und Arbeitsmaschine k_Mgeometrische Konstante M_R(w_M) Reibmoment	Eingangsgrößen: U_AAnkerspannung M_LLastmoment Ausgangsgrößen: n_MMotordrehzahl i_A Ankerstrom
sonstige Größen: U_iinduzierte Spannung w_MWinkelgeschwindigkeit der Motorwelle F_eErregerfluss
Technische Daten des Gleichstrommotors Typ: GFQUJG 080-22 P_mech =4,6 kW n =2500 min^-1 U_A =420 V I_A =13,3 A U_E =360 V I_E =0,8 A

Die geometrische Konstante K_M hängt von der Anzahl der Pole und der Ankerzweige ab. Eine genauere Modellierung würde auch eine Rückwirkung des magnetischen Flusses auf den Erregerstromkreis in Form einer induzierten Spannung berücksichtigen. Diese induzierte Spannung ist jedoch weitaus geringer und wird daher meist vernachlässigt.

3.2.1 Herleitung der Differentialgleichung

Um eine exaktere Modellierung des nichtlinearen Verhaltens vorzunehmen, soll hier nicht der gewohnte lineare Ansatz für die elektromagnetische Energie angesetzt werden, da dieser streng genommen nur im Arbeitspunkt Gültigkeit besitzt.

Die gespeicherte Elektromagnetische Energie ist allgemein:

E_magn. ( t )=

t
ó
õ
0

u_L ( t )·i(t) dt=

t
ó
õ
0

( t )·i(t) dt.

Es gibt einen nichtlinearen Zusammenhang zwischen den Fluss F und dem Strom i,

f = f( i( t ) ).

Das Differential kann im Arbeitspunkt linearisiert werden,

» L·

Wird dieser Ausdruck integriert, erhält man den bekannten linearen Ansatz der elektro- magnetischen Energie.

E_magn. =

t
ó
õ
0

·i dt »

t
ó
õ
0

L·

·i dt=

·L·i²

Es gibt drei dominierende Speicher:

elektromagnetische Energie in der Ankerinduktivität,
elektromagnetische Energie in der Erregerinduktivität
kinetische Energie als Rotationsenergie des Rotors und Arbeitsmaschine.

Für jeden dieser Speicher muss wieder eine Bilanzgleichung aufgestellt werden.

Bilanz der elektromagnetischen Energie im Anker:

é
ë

t
ó
õ
0

df_A

(t )·i_A (t) dt

ù
û

gespeicherte elektro-

magnetische Energie

U_A ·i_A ( t)

zugeführte

Leistung

U_i ·i_A ( t)

abgegebene Leistung

des Motors

R_A ·i_A ( t)²

Abgegebene Leistung

in Wärme der Windungen

L_A ·[d/dt]i_A ( t )=U_A -U_i -R_A ·i_A ( t)

mit U_{i_A} = k_M ·f_E ( t )·w_mech (t ) und w_mech ( t )=2·p·n( t).

Beachten Sie, dass der Erregerfluss nicht mehr als konstant angenommen wird. Dadurch entsteht ein nichtlinearer Term.

Bilanz der elektromagnetischen Energie im Erreger:

é
ë

t
ó
õ
0

df_E

(t )·i_E (t) dt

ù
û

gespeicherte elektro-

magnetische Energie

U_E ·i_E ( t)

zugeführte

Leistung

R_E ·i_E ( t)²

Abgegebene Leistung

in Wärme der Windungen

Da für die erste Differentialgleichung der Erregerfluss gesucht ist, muss die Differentialgleichung durch den Erregerfluss ausgedrückt werden.

[d/dt]f_E ( t )=U_E -R_E ·f( f_E (t)), da i_E (t)=f( f_E (t) ).

Der Erregerfluss ist eine nichtlineare Funktion des Erregerstroms. Die Inverse dieser Funktion soll hier als f( f_E (t) ) bezeichnet werden.

Bilanz der kinetischen Energie:

é
ë

t
ó
õ
0

(t)·w_mech ( t ) dt

ù
û

gespeicherte Rotationsenergie

der Motorwelle und Arbeits-

maschine

M_Motor ·w_mech ( t)

zugeführte mechanische

Leistung des Motors

M_Reib ·w_mech ( t)

abgegebene Reib-

leistung des Motors

M_Last ·w_mech ( t)

abgegebene mechanische

Leistung des Motors

J·[d/dt]w_mech ( t )=M_Motor -M_Reib-M_Last ,

mit [dL/dt]( t ) » J·[(dw_mech )/dt]( t ),

und M_Motor = k_M ·f_E ( t )·i_A ( t), M_Reib = M_Reib ( w_mech (t) ).

Damit ergibt sich das folgende nichtlineare Modell.

di_A ( t )

R_A

L_A

·i_A ( t)-

k_M

L_A

·f_e ( t )·w_mech( t )+

L_A

·U_A ( t )

dw_mech ( t )

k_M

·f_e( t )·i_A ( t )-

·M_Reib (w_mech (t) )-

·M_Last ( t )

df_E ( t )

=-R_E ·f( f_e ( t) )+U_E ( t )

3.2.2 Statische Beziehungen

Die statischen Beziehungen sind ähnlich denen des linearen Models. Die zusätzliche Gleichung liefert mit U_e eine weitere Eingangsgröße und somit eine weitere Eingriffsmöglichkeit für die Drehzahlregelung.

U_A ( t )=R_A ·i_A ( t )+k_M ·f_e( t )·

2·p

·n( t )

M_Last ( t )=k_M ·f_e ( t )·i_A( t )-M_Reib

æ
è

2·p

·n( t)

ö
ø

f_e ( t )=f^-1

æ
è

U_E ( t )

R_E

ö
ø

3.2.3 Reibkennlinie

In linearen Modellen wird eine viskose Reibung angenommen. Für einen nichtlinearen Ansatz, wird meist die Striebeck- Kennlinie verwendet. Man kann diese mathematisch, wie folgt, beschreiben.

M_Reib ( w( t ) )=( b₁ ·e^{-[(| w( t ) |)/(b₂ )]}+b₃ ·|w( t ) |+b₄ )·

w( t)

| w( t ) |

Lassen Sie sich diese Reibkennlinie in Matlab anzeigen. Benutzen Sie die folgenden Anweisungen.

% Striebeck- Kennlinie

b1=1.0; % [Nm]

b2=5.0; % [s/rad]

b3=0.002; % [Nm*s/rad]

b4=0.5; % [Nm]

w=linspace(-2*pi/60*3000,2*pi/60*3000);

figure;

plot(60/(2*pi)*w,(b1*exp(-abs(w)/b2)+b3*abs(w)+b4).*sign(w));

xlabel('n / [1/min]');

ylabel('M_{Reib} / [ Nm]');

title('Striebeck- Kennlinie');

axis([-3000 3000 -1.5 1.5]);

grid on;

Statische Reibkennlinie

Fügen Sie die Parameter der Striebeckkennlinie der Parameterdatei zu!

3.2.4 Erregerflusskennlinie

Für den Erregerfluss in Abhängigkeit des Erregerstromes sind folgende Messwerte bekannt.

i_e / A	-0,85	-0,64	-0,52	-0,36	-0,23	-0,09	0	0,09	0,24	0,36	0,5	0,69	0,85
F / Vs	-12,9	-12,1	-11,5	-10,1	-8,2	-3,9	0	4,1	8,3	10,2	11,4	12,3	12,8

Lassen Sie sich auch diese Kennlinie anzeigen mit Hilfe der folgenden Anweisungen.

% Erregerfluss- Kennlinie

i=[-0.85;-0.64;-0.52;-0.36;-0.23;-0.09;0;0.09;0.24;0.36;0.5;0.69;0.85];

phi=[-12.9;-12.1;-11.5;-10.1;-8.2;-3.9;0;4.1;8.3;10.2;11.4;12.3;12.8];

figure;

plot(i,phi);

ylabel('\phi( i_e) / [ Vs]');

xlabel('i_e / [ A]');

title('Erregerfluss- Kennlinie');

axis([-0.8 0.8 -15 15]);

grid on;

Erregerfluss- Kennlinie

3.2.5 Der Signalflussplan in Simulink

Um nun das dynamische Verhalten des nichtlinearen Modells zu untersuchen, muss der Signalflussplan entwickelt werden. Man geht wieder von den einzelnen Integratoren der Zustandsgleichung ( 13) aus. Vor jeden Eingang eines Integrators wird eine Mischstelle gesetzt. Unter Beachtung des Vorzeichens, fließen in ihr alle Terme die auf der rechten Seite der DGL stehen. Dies ergibt den folgenden Signalflussplan.

Signalflussplan des nichtlinearen Gleichstrommotors

Öffnen Sie ein neues Simulink- Blockschaltbild, indem Sie den Simulinkbrowser öffnen und File® New® Model wählen. Speichern Sie das Modell unter motor_nl.mdl ab. Holen Sie sich aus der Blocklibrary die notwendigen Blöcke, um das Signalflussbild zu zeichnen. Die Eingänge sind wiederum ,,InPort"- Blöcke aus Sources und die Ausgänge sind ,,OutPort"- Blöcke aus Sinks. Beachten Sie die Numerierung!

Für die nichtlineare Reib- Kennlinie holen Sie sich bitte aus der Block- Library unter Functions & Tables den Block Fcn. In diesem können Sie direkt eine mathematische Funktion in Abhängigkeit des Einganges eingeben. Benutzen Sie folgende Formel:

(b1*exp(-abs(u(1))/b2)+b3*abs(u(1))+b4)*sign(u(1)).

Die Erregerflusskennlinie ist nur durch einige Messwerte bekannt. Sie muss daher interpoliert werden. Holen Sie sich aus der Block- Library unter Functions & Tables den Block ,,Look- up Table". In ,,Vector of input variable" geben Sie die Variable phi ein und in ,,Vector of output variable" die Variable i. Sie sollten beide Variablen in der Parameterdatei definieren. Für die Multiplikation findet sich der passende Block ,,Product" in ,,Math". Fügen Sie der Parameterdatei noch die folgenden Variablen hinzu.

Km=0.104; % geometrische Konstante

Re=Ue/Ie; % Erregerwiderstand [Ohm]

3.2.6 Linearisierung im Arbeitspunkt

Damit im Frequenzbereich Aussagen über die Systemeigenschaften in einem Arbeitspunkt getroffen werden können, muss das System linearisiert werden. Um den Ruhepunkt für eine bestimmte Eingangssignalkombination zu finden, steht die Funktion trim zur Verfügung. Dieser Ruhepunkt kann anschließend der Funktion linmod übergeben werden, um die Systemmatrizen zu finden.

Es soll im Arbeitspunkt ein Lastmoment von 10 Nm bei einer Drehzahl von 2500 u/min auftreten. Mit Hilfe der statischen Beziehungen des linearen Modells, werden der Ankerstrom, die Ankerspannung und die Erregerspannung geschätzt. Diese stehen als Startwerte der Funktion trim zur Verfügung. Nutzen Sie folgende Anweisungen.

% Werte im Arbeitspunkt

M0=10; % Lastmoment [Nm]

n0=2500; % Drehzahl [1/min]

Ia0=M0/KPhi; % Ankerstrom [A]

Ua0=Ra*Ia0+KPhi*2*pi/60*n0; % Ankerspannung [V]

Ue0=Ue; % Erregerspannung [V]

phi0=KPhi/Km; % magnetischer Fluss [Vs]

% Berechne Arbeitspunkt

[X0,U0,Y0,DX0]=trim('motor_nl',[Ia0;2*pi/60*n0;phi0],[Ua0;M0;Ue0],[Ia0;n0;phi0],2,[2; 3],2);

% Linearisiere Modell im Arbeitspunkt

[A,B,C,D]=linmod('motor_nl',X0,U0,Y0);

% Bode- Diagramm

figure;

bodemag(ss(A,B,C,D));

grid on;

Man erhält für diesen Arbeitspunkt die untenstehenden Amplitudengänge. Es ist sinnvoll die Vektoren der Arbeitspunktvariablen sich anzuschauen und zu deuten.

>> U0

>> X0

>> Y0

Amplitudengänge des Gleichstrommotors in verschiedenen Arbeitspunkten

Um die Änderung des linearen Verhaltens in Abhängigkeit des Arbeitspunktes zu visualisieren, können mit den nachstehenden Anweisungen dem Amplitudengang zwei weitere Arbeitspunkte hinzugefügt werden. Diese sind im obigen Bild schon eingetragen.

% Berechne Arbeitspunkt bei Nennwerten

[X0,U0,Y0,DX0]=trim('motor_nl',[Ia;2*pi/60*nn;phi0],[Ua;Mn;Ue],[Ia;nn;phi0],2,[2; 3],2);

% Linearisiere Modell bei Nennwerten

[A,B,C,D]=linmod('motor6',X0,U0,Y0);

% füge Bode- Diagramm hinzu

hold on;

bodemag(ss(A,B,C,D));

% Berechne Arbeitspunkt bei Leerlauf

[X0,U0,Y0,DX0]=trim('motor_nl',[0;2*pi/60*nn;phi0],[Ua;0;Ue],[Ia;nn;phi0],[],[2]);

% Linearisiere Modell bei Leerlauf

[A,B,C,D]=linmod('motor6',X0,U0,Y0);

% füge Bode- Diagramm hinzu

bodemag(ss(A,B,C,D));

hold off;

Die Amplitudengänge zeigen auf, dass der Gleichstrommotor sehr gut durch ein lineares Modell beschrieben werden kann. Lediglich im Leerlauf ist der Einfluss der Erregerspannung auf dem Ankerstrom in den unteren Frequenzen geringer. Je stärker die Belastung wird, desto größer wird ihr Einfluss.

3.3 Hausaufgabe

Entwickeln Sie selbständig eine Blockmaskierung für den nichtlinearen Gleichstrommotor. Die Blockmaskierung kann wie folgt aussehen.

Maskierung des nichtlinearen Gleichstrommotors

Simulieren Sie den Anlaufvorgang in den Arbeitspunkt. Folgendes Blockschaltbild kann verwendet werden.

Blockschaltbild

Sie sollten folgendes Ergebnis erhalten.

Ankerstrom

Drehzahl

Erregerstrom

Der hohe Anlaufstrom ist eine große Herausforderung. Die Steuerelektronik muss für den Maximalwert ausgelegt werden, welches hohe Kosten mit sich zieht. Darum ist es günstiger, während der Anlaufphase den Ankerstrom zu begrenzen. Dies erhöht zwar die Anlaufzeit, sinkt jedoch die Kosten der Steuerelektronik erheblich.

Fügen Sie diesen Block anstelle des linearen Gleichstrommotors in die Kaskadenregelung ein. Vergleichen Sie das Simulationsergebnis mit dem des linearen Gleichstrommotors.

4 Seminar

In diesem Seminar sollen einige für den Regelungstechniker brauchbare Techniken vorgestellt werden.

4.1 Die Symbolic Math Toolbox

Viele Reglersynthese- Verfahren setzen ein lineares Prozessmodell voraus. Um diese anzuwenden, kann man versuchen das nichtlineare System im Arbeitspunkt zu Linearisieren. Für Systeme hoher Ordnung kann die Linearisierung per Hand sehr zeitaufwendig sein. Aus diesem Grunde, soll dies mit Hilfe der Symbolic Math Toolbox geschehen.

Der nichtlineare Gleichstrommotor lässt sich durch die folgenden Differentialgleichungen beschreiben.

di_A ( t )

R_A

L_A

·i_A ( t)-

k_M

L_A

·f_e ( t )·w_mech( t )+

L_A

·U_A ( t )

dw_mech ( t )

k_M

·f_e( t )·i_A ( t )-

·M_Reib (w_mech (t) )-

·M_Last ( t )

df_E ( t )

=-R_E ·f( f_e ( t) )+U_E ( t )

Mit den nichtlinearen statischen Funktionen

M_Reib ( w( t ) )=( b₁ ·e^{-[(| w( t ) |)/(b₂ )]}+b₃ ·|w( t ) |+b₄ )·

w( t)

| w( t ) |

i_e ( f( t ) )=

·tan

æ
è

f( t )

f_max

ö
ø

4.1.1 Definition von symbolischen Variablen

Bevor man symbolische Berechnungen durchführen kann, müssen alle verwendeten Variablen als Symbol definiert werden. Dies kann durch die Anweisung syms geschehen. Definieren Sie symbolisch die Variablen des Gleichstrommotors.

% Definiere symbolische Variablen

syms ia wm phi Ua Mlast Ue; % Signale des Gleichstrommotors

syms Ra Re La Km J; % Parameter des Gleichstrommotors

syms b1 b2 b3 b4; % Parameter der Reibkennlinie

syms phi_max m; % Parameter der Erregerflusskennlinie

Hat man die Variablen definiert, kann man diese in symbolischen Berechnungen verwenden. Geben Sie nun die Gleichungen der nichtlinearen statischen Beziehungen für das Reibmoment sowie den Erregerfluss ein.

% Definiere Striebeck- Kennlinie

Mreib=(b1*exp(-abs(wm)/b2)+b3*abs(wm)+b4)*wm/abs(wm);

% Definiere inverse Magnetisierungskennlinie

ie=1/m*tan(pi/(2*phi_max)*phi)

Nachdem die Kennlinien definiert sind, kann man sie bei der Definition der Differentialgleichungen weiterverwenden. Definieren Sie die Systemfunktion.

% Definiere Systemfunktion f(x,u)

f1=-Ra/La*ia-Km/La*phi*wm+1/La*Ua;

f2=Km/J*phi*ia-1/J*Mreib-1/J*Mlast;

f3=-Re*ie+Ue;

% Fasse die Funktionen zu einer Vektorfunktion zusammen

f=[f1;f2;f3];

4.1.2 Linearisierung im Arbeitspunkt

Nun soll die Linearisierung im Arbeitspunkt vorgenommen werden. Es ist zweckmäßig Arbeitspunktvariablen und Abweichungsvariablen zu definieren. Nehmen Sie dieses vor.

% Definiere Variable im Arbeitspunkt

syms ia0 wm0 phi0 Ua0 Mlast0 Ue0;

% Definiere Abweichungsvariable vom Arbeitspunkt

syms dia dwm dphi dUa dMlast dUe;

Um diese Variablen benutzen zu können, muss noch der Zusammenhang zwischen ihnen hergestellt werden. Anschließend werden diese Zusammenhänge in der Systemfunktion substituiert.

% Definiere Zusammenhang

ia=ia0+dia;

wm=wm0+dwm;

phi=phi0+dphi;

Ua=Ua0+dUa;

Mlast=Mlast0+dMlast;

Ue=Ue0+dUe;

% Substituiere in Gleichung

f=subs(f)

Nachdem die vollständige Systemgleichung definiert wurden ist, können die Jakobi- Matrizen berechnet werden. Die Ableitungen müssen nach den Abweichungsvariablen vorgenommen werden.

% Berechne symbolisch Jacobi Matrix

A = jacobian(f, [dia dwm dphi])

B = jacobian(f, [dUa dMlast dUe])

Diese Jakobi- Matrizen entsprechen den Systemmatrizen des linearen Systems im Arbeitspunkt. Die Abweichungsvariablen im Arbeitspunkt sind natürlich null. Deswegen werden diese nun numerisch als null definiert und deren Werte anschließend in den Matrizen substituiert.

% Setze für Matrizen die Abweichungen zu null

dia=0;

dwm=0;

dphi=0;

dUa=0;

dMlast=0;

dUe=0;

% Substituiere

A = subs(A)

B = subs(B)

Damit sind die Systemmatrizen A und B schon gefunden. Die Linearisierbarkeit des nichtlinearen Systems im Arbeitspunkt kann durch die Determinante überprüft werden. Diese muss verschieden von null sein.

% Berechne symbolisch die Determinante

detA = simple(det(A));

4.1.3 Symbolische Lösung der algebraischen Gleichungen

Angenommen im Arbeitspunkt ist Drehzahl, Lastmoment und Erregerspannung vorgegeben, dann kann man die Systemfunktion nach Ankerspannung, Ankerstrom und Fluss auflösen. Hierzu muss die Systemfunktion im Arbeitspunkt betrachtet werden, die als f0 bezeichnet wird.

f0 = subs(f)

Durch den Befehl solve, kann nun die Systemfunktion im Arbeitspunkt nach den gesuchten Größen aufgelöst werden.

[Ua0,ia0,phi0]=solve(f0(1),f0(2),f0(3),Ua0,ia0,phi0)

Mit dem Befehl simple, kann man versuchen die Gleichungen kompakter umzuformen.

z.B.

>> simple( Ua0)

4.1.4 Numerische Lösung

Sind zum Schluss auch noch numerische Lösungen gesucht, lassen sich die Parameter, wie gewohnt, numerisch definieren. Mit dem subs Befehl kann man dann diese in den symbolischen Gleichungen einsetzen.

% Definiere numerisch die Parameter Ruhepunkt

La=0.0402; % Ankerinduktivität [H]

J=15; % Trägheitsmoment Motor [kg*m^2]

Km=0.104; % geometrische Konstante

Ra=5.6; % Ankerwiderstandes in [Ohm]

Re=450; % Erregerwiderstand [Ohm]

% Parameter der Erregerfluss- Kennlinie

phi_max=15; % max. Fluss [Vs]

m=5; % Anstieg im Nullpunkt

% Parameter der Striebeck- Kennlinie

b1=1.0; % [Nm]

b2=5.0; % [s/rad]

b3=0.002; % [Nm*s/rad]

b4=0.5; % [Nm]

% Vorgegebene Werte im Ruhepunkt

wm0=2*pi/60*2500; % [rad/s]

Mlast0=10; % [Nm]

Ue0=360; % [V]

% Berechne die anderen Werte

Ua0=subs(Ua0)

ia0=subs(ia0)

phi0=subs(phi0)

A=subs(A)

B=subs(B)

detA=subs(detA);

4.2 Die Control- Toolbox

Die Control- Toolbox erlaubt das komfortable Definieren von linearen zeitinvarianten Systemen ( kurz LTI in der Matlab- Umgebung). Nach ihrer Definition stehen dem Anwender verschiedene Werkzeuge für die Manipulation und der visuellen Darstellung zur Verfügung. Einige sollen kurz vorgestellt werden.

4.2.1 Definition von Transferfunktionen

Eine oft benutzte Darstellung von LTI- Objekten ist die Polynomdarstellung. Ähnlich wie in Simulink, kann man sie durch zwei Spaltenvektoren mit den gewünschten Polynomkoeffizienten definieren. Dafür verwendet man den Befehl tf. Geben Sie folgende Anweisung im Befehlszeilenfenster von Matlab ein.

>> Gs = tf([1],[2 1])

Sofort wird die Transferfunktion in der üblichen Darstellung angezeigt. Auf dieser Weise lassen sich nun diverse Transferfunktionen definieren. Ein statischer Übertragungsfaktor als Transferfunktion kann, wie folgt, definiert werden.

>> Kp = tf([5],[1])

Für den Regelungstechniker ist es oft erforderlich die Transferfunktion des geschlossenen Regelkreises zu berechnen. Mit LTI- Objekten ist dieses sehr einfach. Angenommen nach der Definition einer Regelstrecke und eines P- Reglers soll die Transferfunktion des geschlossenen Regelkreises berechnet werden, dann erfolgt dies unter Anwendung der Mason- Formel durch folgende Anweisung.

>> Gc=Kp*Gs/(1+Kp*Gs)

Das Ergebnis entspricht nicht den Erwartungen. Es entsteht eine Pol- Nullstellen Kompensation. Um das korrekte Ergebnis zu erhalten ist der Befehl minreal erforderlich.

>> Gc=minreal(Kp*Gs/(1+Kp*Gs))

Es können auch Totzeiten angegeben werden. Leider kann man solche Systeme nicht mehr allen mathematischen Operationen unterziehen. Eine Strecke mit Totzeit kann folgendermaßen angegeben werden.

>> Gs= tf([5],[2 1],'iodelay',3)

4.2.2 Definition von Zustandsraumdarstellungen

LTI- Objekte können auch in Zustandsraumdarstellung definiert werden. Hierfür müssen dem Befehl ss die Zustandsraummatrizen übergeben werden. Benutzen Sie folgende Anweisungen.

>> A=-0.5;

>> B=1;

>> C=0.5;

>> D=0;

>> Gs = ss(A,B,C,D)

Auch hier sind Zusammenschaltungen möglich. Eine Reihenschaltung zweier Systeme ergibt folgendes Gesamtsystem.

>> Gs*Gs

Geben Sie nun den Regler in Zustandsraumdarstellung an.

>> A=0;

>> B=0;

>> C=0;

>> D=5;

>> Kp=ss(A,B,C,D)

Der geschlossene Regelkreis in Zustandsraumdarstellung ergibt.

>> Gc=Kp*Gs/(1+Kp*Gs)

Welches wiederum nicht in der minimalen Realisation vorliegt. Das korrekte System lautet, wie folgt.

>> Gc=minreal(Kp*Gs/(1+Kp*Gs))

Auch in der Zustandsraumdarstellung können Totzeiten angegeben werden. Hier existieren zwei Möglichkeiten. Eingangstotzeiten und Ausgangstotzeiten. Das Ein- Ausgangsverhalten ist bei beiden das Gleiche, die Zustandstrajektorie ist jedoch verschieden.

>> A=-0.5;

>> B=1;

>> C=0.5;

>> D=0;

>> Gs1 = ss(A,B,C,D,'inputdelay',3)

>> Gs2 = ss(A,B,C,D,'outputdelay',3)

4.2.3 Übergangsfunktion und Gewichtsfunktion

Nachdem der geschlossene Regelkreis als LTI- System vorliegt, kann seine Übergangs-funktion und Gewichtsfunktion einfach angezeigt werden.

>> step(Gc)

>> impulse(Gc)

4.2.4 Bode- Diagramm, Ortskurve und Pole- Zero- Darstellung

Auch für die Anzeige des Bode- Diagramms und der Ortskurve stehen Befehle zur Verfügung. Probieren Sie die nachfolgenden Anweisungen aus.

>> bode(Gc)

>> grid on

>> nyquist(Gc)

>> pzmap(Gc)

4.2.5 Diskretisierung von kontinuierlichen Systemen

Heutzutage bedient man sich digitaler Rechner zur Regelung kontinuierlicher Prozesse. Deswegen muss das Prozessmodel diskretisiert werden. Zuvor soll hier auch noch eine andere komfortable Eingabemöglichkeit von Transferfunktionen aufgezeigt werden.

>> s=tf('s');

>> Ks=5;

>> T=2;

>> Gs=Ks/(1+s*T)

Diese definierte kontinuierliche Transferfunktion soll nun mit dem Befehl c2d diskretisiert werden.

>> Ts=0.25;

>> Gz=c2d(Gs,Ts,'zoh');

Alle zuvor aufgezeigte Visualisierungsmöglichkeiten sind auch mit der diskreten Transferfunktion möglich. Probieren Sie es aus.

4.2.6 Simulation von Systemen

Möchte man das Ausgangssignal zu einem beliebigen Eingangssignal berechnen, so steht der Befehl lsim zur Verfügung. Die Antwort auf ein sinusförmiges Eingangssignal für das vorherige System, erfolgt durch folgende Anweisungen.

>> s=tf('s');

>> Ks=5;

>> T=2;

>> Gs=Ks/(1+s*T)

>> t=[0:0.25:20]';

>> u=sin(2*pi/5*t);

>> [y,t]=lsim(Gs,u,t);

>> plot(t,y);

>> xlabel('t / [ s]');

>> ylabel('y( t)');

>> title('Ausgangssignal');

>> grid on

5 Seminar

5.1 S- Function

Mit Simulink lassen sich auch komplexe Systeme simulieren. Ein Nachteil entsteht jedoch in der Ausführungszeit, wenn diese aus einzelnen Blöcken aufgebaut werden. Um dennoch eine schnelle Simulation zu ermöglichen, kann man komplexe Teilsysteme in einer S- Function definieren und diese als einzelnen Block darstellen. Weiterhin besteht die Möglichkeit S- Function neben der Matlab- eigene Programmiersprache auch in C oder Fortran zu erstellen.

5.1.1 Aufruf Syntax

S- Function stellen ein spezielles Interface bereit, über welches die DGL- Löser die notwendigen Informationen über die Systemgleichungen ermitteln. Das Interface wurde dabei sehr allgemein ausgelegt, um stetige, diskrete und hybride Systeme zu definieren. Damit lassen sich eine breite Klasse von Systemen in Simulink darstellen. Die S- Function wird mit Hilfe des Blockes S- Function aus der Nonlinear Library in Simulink eingefügt.

Unten ist ein allgemeiner Simulink- Block dargestellt. Er besteht aus einem Eingangsvektor und einem Ausgangsvektor. Innerhalb des Blockes befindet sich auch ein Zustandsvektor.

Der Zustandsvektor kann aus stetigen oder diskreten Zuständen bestehen, sowie auch aus einer Kombination von beiden. Die mathematischen Beziehungen zwischen Ein- und Ausgang lassen sich, wie folgt, beschreiben.

x=[ x_C ,x_D ]¢ Zustandsvektor aus stetigen und diskreten Zuständen

[x\dot]_C = f(t,x,u) Stetige Zustände, Derivative- Funktion

x_D (k+1)=f( t(k),x(k),u(k) ) Diskrete Zustände, Update- Funktion

y=g(t,x,u) Ausgangsfunktion

Die einzelnen Simulationsschritte werden nach dem folgenden Ablaufplan abgearbeitet. Der S- Function wird dabei ein Flag übergeben, welches den aktuellen Simulationsschritt angibt. Der Nutzer muss anschließend die entsprechende Information über das System zurückgeben.

Abarbeitung einer S- Function

5.1.2 S- Function am Beispiel des linearen Gleichstrommotors

Es soll hier noch mal das Beispiel des linearen Gleichstrommotors von Seminar 3 aufgegriffen werden. Die beschreibende Differentialgleichung ist noch einmal unten angegeben.

æ
ç
ç
è

Di_A ( t )

Dw_mech ( t )

ö
÷
÷
ø

é
ê
ê
ë

-[(R_A )/(L_A )]

-[(k_M ·f_e )/(L_A )]

[(k_M ·f_e )/J]

-[b/J]

ù
ú
ú
û

æ
ç
ç
è

Di_A ( t )

Dw_mech ( t )

ö
÷
÷
ø

é
ê
ê
ë

[1/(L_A )]

-[1/J]

ù
ú
ú
û

æ
ç
ç
è

DU_A ( t )

DM_Last ( t )

ö
÷
÷
ø

æ
ç
ç
è

Di_A ( t )

Dn( t )

ö
÷
÷
ø

é
ê
ê
ë

[60/(2·p)]

ù
ú
ú
û

æ
ç
ç
è

Di_A ( t )

Dw_mech ( t )

ö
÷
÷
ø

x=A·x+B·u

y=C·x

Erzeugen Sie eine neue M- Datei mit dem Namen sfun_motor.m. Schreiben Sie anschließend den folgenden Programmtext hinein.

>> edit sfun_motor

function [sys,x0,str,ts] = sfun_motor(t,x,u,flag)

% sfun_motor S-Funktion welche ein Gleichstrommotor definiert.

switch flag,

%%%%%%%%%%%%%%%%%%

% Initialization %

%%%%%%%%%%%%%%%%%%

case 0,

[sys,x0,str,ts]=mdlInitializeSizes;

%%%%%%%%%%%%%%%

% Derivatives %

%%%%%%%%%%%%%%%

case 1,

sys=mdlDerivatives(t,x,u);

%%%%%%%%%%

% Update %

%%%%%%%%%%

case 2,

sys=mdlUpdate(t,x,u);

%%%%%%%%%%%

% Outputs %

%%%%%%%%%%%

case 3,

sys=mdlOutputs(t,x,u);

%%%%%%%%%%%%%%%%%%%%%%%

% GetTimeOfNextVarHit %

%%%%%%%%%%%%%%%%%%%%%%%

case 4,

sys=mdlGetTimeOfNextVarHit(t,x,u);

%%%%%%%%%%%%%

% Terminate %

%%%%%%%%%%%%%

case 9,

sys=mdlTerminate(t,x,u);

%%%%%%%%%%%%%%%%%%%%

% Unexpected flags %

%%%%%%%%%%%%%%%%%%%%

otherwise

error(['Unhandled flag = ',num2str(flag)]);

end

%================================================================

% mdlInitializeSizes

% Gibt die Anzahl der Zustände, Anfangswerte, und Abtastzeit der S-function zurück.

%================================================================

function [sys,x0,str,ts]=mdlInitializeSizes

sizes = simsizes;

sizes.NumContStates = 2; % Anzahl der stetigen Zustände

sizes.NumDiscStates = 0; % Anzahl der diskreten Zustände

sizes.NumOutputs = 2; % Anzahl der Ausgänge

sizes.NumInputs = 2; % Anzahl der Eingänge

sizes.DirFeedthrough = 0; % FLAG ob D Matrix existiert

sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

% Initialisierung der Anfangswerte

x0=zeros(sizes.NumContStates,1);

% str ist eine leere Matrix

str = [];

% Initialisierung des Vektors der gewünschten Abtastzeiten

ts = [0 0]; % Abtastzeit: [Periode, Offset]

% end mdlInitializeSizes

%================================================================

% mdlDerivatives

% Berechnet die Ableitungen der stetigen Zustände.

%================================================================

function dx=mdlDerivatives(t,x,u)

% Definition der Parameter

Ra=5.5741; % Ankerwiderstand [Ohm]

La=0.0402; % Ankerinduktivität [H]

J=15; % Trägheitsmoment Motor [kg*m^2]

b=0; % Viskosität [Nms]

KPhi = 1.3211; % Motorkonstante [Vs]

% Definition der Differentialgleichung

dx=zeros(2,1);

dx(1)=-Ra/La*x(1) - KPhi/La*x(2) + 1/La*u(1);

dx(2)= KPhi/J*x(1) - b/J*x(2) - 1/J*u(2);

% end mdlDerivatives

%================================================================

% mdlUpdate

% Berechne den nächsten diskreten Zustand, sowie Abtastzeitpunkt

%================================================================

function sys=mdlUpdate(t,x,u)

sys = [];

% end mdlUpdate

%================================================================

% mdlOutputs

% Berechne den Ausgang des Blockes.

%================================================================

function y=mdlOutputs(t,x,u)

y = [1,0;0,60/(2*pi)]*x;

% end mdlOutputs

%================================================================

% mdlGetTimeOfNextVarHit

% Berechnet nächsten Abtastzeitpunkt der diskreten Zustände

%================================================================

function sys=mdlGetTimeOfNextVarHit(t,x,u)

sys = [];

% end mdlGetTimeOfNextVarHit

%================================================================

% mdlTerminate

% Räume vor dem Ende der Simulation auf.

%================================================================

function sys=mdlTerminate(t,x,u)

sys = [];

% end mdlTerminate

Nachdem Sie den Programmtext für die S- Function eingegeben haben, soll dieser nun als ein Simulink- Block verwendet werden. Dazu öffnen Sie ein neues leeres Simulink- Fenster. Holen Sie aus der Blocklibrary unter Functions & Tables den Block S- Function. Geben Sie im Parameterdialog des Blockes unter ,,S- function name" den Namen der S- Function ( sfun_motor) ein. Da keine zusätzlichen Parameter definiert wurden sind, bleibt das Parameterfeld leer. Holen Sie aus der Blocklibrary unter Signals&Systems die Blöcke Mux und Demux. Mit ihnen ist es möglich die Signale zusammenzufassen bzw. zu extrahieren. Anschließend kann man die gleiche Beschaltung vornehmen, wie im Seminar 3. Mit einer Simulationszeit von 200s kann die Simulation gestartet werden. Es sollten die gleichen Ergebnisse, wie im Seminar 3, erreicht werden.

% Werte im Arbeitspunkt

M0=10; % Lastmoment [Nm]

n0=2500; % Drehzahl [1/min]

Ia0=7.57; % Ankerstrom [A]

Ua0=388; % Ankerspannung [V]

Simulink Blockschaltbild

5.1.3 Parameterübergabe und Blockmaskierung

Auch an S- Functions können Parameter übergeben werden. Diese werden als zusätzliche Argumente der S- Funktion deklariert. Ändern Sie in der S- Function die folgenden Zeilen ab.

function [sys,x0,str,ts] = sfun_motor2(t,x,u,flag,Ra,La,KPhi,J,b)

% sfun_motor S-Funktion welche ein Gleichstrommotor definiert.

...

case 1,

sys=mdlDerivatives(t,x,u,Ra,La,KPhi,J,b);

...

function dx=mdlDerivatives(t,x,u,Ra,La,KPhi,J,b)

% Definition der Differentialgleichung

dx=zeros(2,1);

dx(1)=-Ra/La*x(1) - KPhi/La*x(2) + 1/La*u(1);

dx(2)= KPhi/J*x(1) - b/J*x(2) - 1/J*u(2);

% end mdlDerivatives

...

Anschließend muss im Parameterdialog des S- Function Blocks unter ,,S- functions Parameter" die definierte Reihenfolge der Parameter ( Ra, La, KPhi, J, b) angegeben werden. Nachdem die numerischen Werte in der Matlab- Umgebung definiert werden, lässt sich das Simulink- Modell wieder verwenden.

>> Ra=5.5741; % Ankerwiderstand [Ohm]

>> La=0.0402; % Ankerinduktivität [H]

>> J=15; % Trägheitsmoment Motor [kg*m^2]

>> b=0; % Viskosität [Nms]

>> KPhi = 1.3211; % Motorkonstante [Vs]

Damit S- Function auch komfortabel parametrisiert werden können, ist es günstig, eine Blockmaskierung vorzunehmen. Markieren Sie durch Ziehen mit der Maus den Mux, Demux und den S- Function Block. Gehen Sie ins Menu auf Edit® Create Subsystem. Ändern Sie den Text für die Portbezeichnungen und Blockbezeichnung entsprechend ab.

Simulink Blockschaltbild

Sub- System

Anschließend lässt sich für den neuen Block Motor, wie im Seminar 3 eine Blockmaskierung vornehmen.

5.2 Klassische Reglersynthese mit Hilfe des SISO- Tools

Nachdem Möglichkeiten zur Darstellung von Systemen sowie deren Linearisierung aufgezeigt wurden, soll hier eine kurze Einführung in die klassische Reglersynthese erfolgen.

5.2.1 5.2.1 Einführung in das Wurzelortsverfahren

Das Wurzelortsverfahren ist ein grafisches Verfahren, mit dem einfache lineare Regelkreise analysiert und dimensioniert werden können. Mit Hilfe des Wurzelortsverfahrens lassen sich aus den bekannten Eigenschaften des offenen Regelkreises Aussagen über die Lage der Wurzeln der charakteristischen Gleichung des geschlossenen Regelkreises in Abhängigkeit eines Parameters, meist der Regelverstärkung, machen.

Einfacher linearer Regelkreis

Unter der Anwendung der Formel von Mason, lässt sich die Übertragungsfunktion des geschlossenen Regelkreises, wie folgt, angeben.

G_c ( s )=

K·G₀ ( s )

1+K·G₀ (s )

Die charakteristische Gleichung des geschlossenen Regelkreises lautet somit:

P( s )=0=1+K·G₀ ( s ).

In Abhängigkeit des Parameters K ändert sich die Lage der Wurzeln dieser charakteristischen Gleichung. Lässt man K von null bis unendlich laufen und verbindet alle Änderungen der Wurzelorte, erhält man die Wurzelortskurve.

Beispiel:

Offener Regelkreis ohne K: G₀ ( s )=[1/(1+s)]

Charakteristische Gleichung: P( s )=0=1+[K/(1+s)]=1+s+K

Wurzeln in Abhängigkeit von K: s( K )=-1-K

Wurzelortskurve

In diesem Beispiel verbleibt die Wurzelortskurve auf der negativen reellen Achse. Sie beginnt in der Polstelle des offenen Regelkreises und wandert zu minus unendlich. Der Regelungstechniker kann durch die Auswahl eines bestimmten Verstärkungsfaktors K den Pol des geschlossenen Regelkreises auf einem Punkt dieser Kurve festlegen.

5.2.2 5.2.2 Konstruktionsregeln für das Wurzelortsverfahren

Obwohl man in heutiger Zeit die Wurzelortskurve nicht mehr per Hand konstruiert, sondern mit der Hilfe eines Computers, ist es für das Verständnis hilfreich die von Evans eingeführten Konstruktionsregeln zu kennen. Aus diesem Grunde sollen hier die wichtigsten kurz aufgeführt werden.

Das charakteristische Polynom muss so umgeformt werden, dass die höchste Potenz in s einen positiven Faktor von eins besitzt.
Da alle Wurzeln reell oder konjugiert komplex sind, verläuft die Wurzelortskurve symmetrisch zur reellen Achse.
Die Zahl der Äste der Wurzelortskurve ist gleich der Zahl n der Pole der Übertragungsfunktion des aufgeschnittenen Regelkreises G₀.
Die Äste der Wurzelortskurve sind kontinuierliche Kurven. Sie beginnen mit K=0 in den Polen von G₀, da für s=sk die Betragsbeziehung gleich null ist. Mit K® ¥enden m Äste in den Nullstellen von G₀ und n-m Äste im unendlichen.
Wurzelorte auf der reellen Achse liegen für K > 0 links von einer geraden Anzahl von Polen und Nullstellen.
Der Schnittpunkt der Wurzelortskurve mit der imaginären Achse kennzeichnet die Stabilitätsgrenze mit ( K_krit, w_krit). Bei diesem K_krit führt der geschlossene Regelkreis ungedämpfte Dauerschwingungen mit der Frequenz w_krit aus.

5.2.3 Aufruf des SISO- Tools

Mit dem Befehl sisotool kann man eine interaktive Benutzeroberfläche öffnen, mit der man einfache Regelkreise dimensionieren kann. Zuerst muss die Regelstrecke als ein lineares Modell beschrieben werden. Dies kann durch Linearisierung im Arbeitspunkt geschehen. Definieren Sie folgendes lineares Modell in der Matlab Kommandozeile.

>> s=tf('s');

>> Ks=5;

>> T=1;

>> D=0.5;

>> Td=0.6;

>> Gs=Ks*(1-Td*s)/(1+s*2*D*T+s^2*T^2);

Rufen Sie anschließend das SISO- Tool, wie folgt, auf.

>> sisotool(Gs)

In der Standardkonfiguration sollte folgendes Fenster erscheinen.

Grafisches Interface der SISO- Tools

In der Wurzelortskurve ist deutlich zu ersehen, dass die Wurzelorte des geschlossenen Regelkreises in den Polen des aufgeschnittenen Regelkreises beginnen. Mit zunehmenden K wandert ein Pol in die Nullstelle des aufgeschnittenen Regelkreises und ein Pol ins unendliche. Die gegebene Strecke besitzt eine Nullstelle in der rechten Halbebene. Solche nichtminimalphasige Systeme zeichnen sich immer dadurch aus, dass mindestens ein Pol mit zunehmender Verstärkung in die rechte Halbebene wandert. Ab einer bestimmten Verstärkung wird der Regelkreis instabil. Dadurch ist die Kreisverstärkung bei einem nichtminimalphasigen System limitiert. Somit ist auch die erreichbare Regelgüte des Regelkreises bei einfachen Reglerstrukturen limitiert.

5.2.4 P- Regler

Bestimmen Sie als erstes den Verstärkungsfaktor, ab dem der geschlossene Regelkreis instabil wird. Die roten Punkte in der Wurzelortskurve entsprechen den aktuellen Polstellen des geschlossenen Regelkreises. Ziehen Sie mit der Maus einen der roten Punkte zur imaginären Achse. Sind beide Polstellen auf der imaginären Achse, dann ist der Regelkreis grenzstabil, und es entstehen am Ausgang ungedämpfte Schwingungen. Durch das Verschieben der Pole des geschlossennen Regelkreises, haben sie den Verstärkungsfaktor des P- Reglers verändert. Die kritische Schleifenverstärkung ist 0,333. Nach den Einstellregeln von Ziegler und Nichols sollte die Dimensionierung des P- Reglers die Hälfte dieser Verstärkung sein. Geben Sie mit der Tastatur eine Verstärkung von 0,166 für den P- Regler ein. Im Amplitudengang und im Phasengang können Sie die aktuelle Amplitudenreserve von ca. 6dB und die aktuelle Phasenreserve von 66 ablesen. Gehen Sie auf den Menüpunkt Analysis- > Response to Step Command. Es öffnet sich ein neues Fenster, in dem als blaue Kurve die Sprungantwort des geschlossenen Regelkreises dargestellt wird. Die grüne Kurve ist die resultierende Stellamplitude und soll uns vorerst nicht interessieren. Rufen Sie mit der rechten Maustaste durch das Klicken im Graph das Kontextmenü auf. Wählen Sie Characteristik- > Seetling time. Dadurch wird das Einregelband mit der Einregelzeit angezeigt. Sie können mit der Maus auf die Sprungantwort klicken und sich die gewünschten Koordinaten anzeigen lassen. So lässt sich, wie untenstehend abgebildet, die Überschwingweite bestimmen.

Sprungantwort im LTI- Viewer

5.2.5 P- Regler mit Verzögerungsglied

Was sofort auffällt ist, dass die Regelabweichung mehr als 50% beträgt. Dieses ist natürlich nicht zufrieden stellend. Für die Regelabweichung lässt sich mit Hilfe des Grenzwertsatzes der Laplace- Transformation folgende Formel finden.

e( t® ¥ )=[1/(1+K₀ )] bei Regelkreisen mit proportionaler Schleifenverstärkung.

Damit ergibt sich die notwendige Schleifenverstärkung zu:

K₀ =

e_max

-1

Die maximale Regelabweichung soll nicht größer als 5% vom Sollwert sein. Somit ergibt sich eine Schleifenverstärkung von 19. Da die Verstärkung der Strecke bereits 5 ist, ist eine Reglerverstärkung von 3,8 erforderlich.

Geben Sie diese Reglerverstärkung für den Regler C(s) ein. Die Pole des geschlossenen Regelkreises liegen in der rechten Halbebene. Der Regelkreis ist instabil. Um den Regelkreis zu stabilisieren, kann der Amplitudengang abgesenkt werden. Dies erreicht man durch einfügen eines Tiefpassfilters erster Ordnung. Gehen Sie mit der Maus in die Symbolleiste und klicken Sie auf das einfache rote Kreuz, welches für Add Pol steht. Gehen Sie mit der Maus in den Amplitudengang des Bodeplots, und fügen Sie den Pol irgendwo ein. Anschließend können Sie den Pol so lange nach links verschieben, bis im Phasengang eine Phasenreserve von ca. 60 erscheint. Beachten Sie, dass sich die Reglerverstärkung dabei nicht ändert. Korrigieren Sie diese notfalls.

SISO- Tool Fenster

Sprungantwort

Erst durch diese Frequenzkorrektur, ist es möglich mit einem P- Regler eine zufrieden stellende Regelung aufzubauen.

5.2.6 I- Regler

Mit Hilfe eines I- Reglers soll die Regelabweichung auf ein Sprung zu null gemacht werden. Um diesen Integrator einzufügen, gehen Sie in den Menüpunkt Compensators- > Edit- > C. Löschen Sie zuerst den eingefügten Pol, indem Sie die Auswahlbox delete markieren und anschließend auf apply drücken. Danach sollte wieder die Anfangskonfiguration erscheinen. Drücken Sie nun auf add real pole. Es erscheint eine Edit- Box in der Sie eine null eingeben für den Integrator. Schließen Sie die Dialogbox. Verschieben Sie in der Wurzelortskurve das konjugiert komplexe Poolpaar erneut auf die imaginäre Achse, um die kritische Reglerverstärkung zu ermitteln. Sie ist 0,125. Wählen Sie wieder die Hälfte als Reglerverstärkung.

SISO- Tool Fenster

Sprungantwort

5.2.7 I- Regler mit Frequenzkorrektur

Wenn Sie die Sprungantwort mit der Sprungantwort der prediktiven Regelung vergleichen, stellen Sie fest, dass die Einregelzeit wesentlich länger dauert. Um die Bandbreite des geschlossenen Regelkreises zu erhöhen, muss die Durchtrittsfrequenz des offenen Regelkreises erhöht werden. Dies erfordert eine Frequenzkorrektur im höheren Frequenzbereich.

Problematisch ist hier, dass nahe der Durchtrittsfrequenz ein konjugiert komplexes Polpaar existiert. Eine einfache Möglichkeit besteht darin dieses Poolpaar durch ein konjugiert komplexes Nullstellenpaar zu kompensieren. Klicken Sie in der Symbolleiste auf das Icon ,,Add complex Zero". Das Icon besitzt ein Bindestrich mit zwei Kreisen. Gehen Sie mit dem Mauszeiger in die Wurzelortskurve nahe dem konjugiert komplexen Polpaar der Strecke und fügen Sie die Nullstellen durch einen Mausklick dort ein. Sie können auch nachträglich die Nullstellen durch ziehen mit der Maus auf die gewünschte Position bringen.

Dadurch verändert die Wurzelortskurve ihr aussehen. Wie in den Konstruktionsregeln zu ersehen, wandert das konjugiert komplexe Poolpaar in das konjugiert komplexe Nullstellenpaar. Verschieben Sie die Amplitude des Amplitudengangs so, dass eine Phasenreserve von 60 entsteht. Schauen Sie sich im LTI- Viewer die Sprungantwort an. Die Einregelzeit sollte bedeutend geringer sein. Untenstehend ist eine Möglichkeit abgebildet. Die Einregelzeit dauert ca. 2,5 s.

SISO- Tool Fenster

Sprungantwort

Der entwickelte Regler ist improper und somit nicht realisierbar. Dies erkennt man daran, dass das Zählerpolynom größer als das Nennerpolynom ist. Es muss mindestens eine Polstelle im höheren Frequenzbereich eingefügt werden, damit der Regler realisierbar wird. Klicken Sie in der Symbolleiste auf das Icon Add real pole, welches durch ein einfaches x dargestellt ist. Fügen Sie diesen Pol weit links auf der reellen Achse in der Wurzelortskurve ein. Es ist auch möglich den Pol weit rechts in das Bodediagramm einzufügen. Ein guter Ort liegt bei ca. 10 rad/s. Damit ist ein einfacher linearer Regler gefunden mit den folgenden Daten.

C( s )=0.15·

1+1.1·s+s²

s·(1+0.1·s )

Nachfolgend sind die Wurzelortskurve und die Sprungantwort angegeben.

SISO- Tool Fenster

Sprungantwort

Damit soll die Einführung in das SISO- Tool, als auch das Seminar beendet werden. Es wurden hier einige Möglichkeiten zur Simulation dynamischer Systeme mit Hilfe von Matlab vorgestellt. Eine vollständige Einarbeitung ist in diesem, mittlerweile großen Softwarepaket, innerhalb 5 Seminare nicht möglich. Die vorgestellten Techniken sollen als Orientierung dienen, welche Möglichkeiten Matlab bietet. Eine selbständige Vertiefung in den Techniken, ist durch die Umfangreiche Hilfe leicht möglich.