Simulation – Fach im Modul „Modellbildung und Simulation ... · Fakult at f ur Informatik und...

Fakultat fur Informatik und AutomatisierungInstitut fur Automatisierungs- und Systemtechnik

Fachgebiet Simulation und Optimale Prozesse

Simulation

– Fach im Modul”Modellbildung und Simulation“ –

– Fach im Modul”Simulationstechnik“

Dr.-Ing. Siegbert Hopfgarten

Tel., E-Mail: +49 3677 69-1418, [email protected]: http://www.tu-ilmenau.de/simulation

Vorlesungsskript: Lehre/Vorlesungen, Seminare, Praktika/Simulation

1 / 172

http://www.tu-ilmenau.de/simulation

Inhalt

1 Einleitung - Begriffe, Motivation, Beispiele

2 Modelle

3 Analoge Simulation

4 Digitale Simulation

5 Simulationssprachen und -systeme

2 / 172

1. Einleitung - Begriffe, Motivation, Beispiele

Was ist Simulation?

< lat., simulare > vortauschen

1.) bewusste Vortauschung von Krankheiten oder ubertriebene Darstellungihrer Symptome

2.) Sammelbegriff fur die Darstellung oder Nachbildung physikalischer,technischer, biologischer, psychologischer oder okonomischer Prozesse oderSysteme durch mathematische oder physikalische Modelle, wobei dieUntersuchung des Modelles einfacher, billiger oder ungefahrlicher als die desOriginals ist und die Erkenntnisse Ruckschlusse auf die Eigenschaften desOriginals erlauben

→ Simulation ist das Experimentieren an einem Modell.

3 / 172


Zwei Aspekte der Simulation:

1. Modellbeschreibung (z. B. x(t) = −3x(t) + 2u(t), x(0) = 0.5)

2. Experimentbeschreibung (z. B. t ∈ [0, 10] )

Oft nicht strikt trennbar; durch Software nicht immer unterstutzt

Gefahr: Experimentdurchfuhrung an Modell, das fur diese Experimente nichtgultig ist

Was ist ein System? (siehe auch DIN 19226)

In betrachtetem Zusammenhang gegebene Anordnung von Gebilden(Teilsystemen); Beziehungen zueinander und zur Umwelt

Gebilde: Gegenstande, Denkmethoden oder deren Ergebnisse

Hullflache (gedacht): Abgrenzung zur/Wechselwirkung mit Umwelt

Betrachteter Zusammenhang: durch Aufgabenstellung gegeben; in derRegelungs- und Systemtechnik: Beziehungen durch Wirkungen

Voraussetzungen fur sinnvolle Untersuchungen: Steuerbarkeit(Beeinflussbarkeit) und Beobachtbarkeit (Messbarkeit von Signalen)

4 / 172


Beispiel: Globaler Wasserkreislauf (Aggregationsgrad,Detailliertheitsgrad)Quelle: http://bildungsserver.hamburg.de/

Schwarze Zahlen in 1.000 km3, blaue Zahlen in 1.000 km3/y5 / 172


Beispiel: Talsperren in ThuringenBildquelle: Thur. Talsperrenverwaltung; jetzt: Thur. Fernwasserversorgung

6 / 172


Beispiel: Nordthuringer TalsperrensystemBildquelle: Thur. Talsperrenverwaltung; jetzt: Thur. Fernwasserversorgung

7 / 172


Beispiel: AbwasserreinigungBildquelle: Wasser- und Abwasserzweckverband Jena

8 / 172


Was ist ein Experiment?

Datengewinnung von einemSystem bei Aufschaltung vonEingangssignalen

Nachteil des Experiments amrealen System: Einfluss vonStorgroßen; gewisse Anzahlnutzlicher Ausgangen nichtdirekt messbar

Vorteil der Simulation: Zugriffauf alle Ein–/Ausgangsgroßen

Problem bei der Simulation:Simulationslaufe außerhalb desExperimentgultigkeitsbereichesdes Experimentes liegen

Was ist ein Modell?

Abbildung eines Systems oder einesProzesses (Gesamtheit derVorgange) in ein anderesbegriffliches oder gegenstandlichesSystem

Bsp.: begriffliches System:

x(t) = f (x(t), u(t), t) , x(0) = x0

Bsp: gegenstandliches System:Hochwasseruberlauf an einerTalsperre im Maßstab getreuenNachbau (s. Abb.)

9 / 172


Gegenstandliches ModellBildquelle: Thuringer Fernwasserversorgung

10 / 172


Warum ist Modellierung wichtig?

Modellierung: Prozess des Ordnens von Wissen uber ein System

Modellierung als Problemlosungswerkzeug

”(Natur–)Wissenschaftler sind glucklich, wenn sie die Welt beobachten und

verstehen konnen.“ (Analyseaspekt)

”Ingenieure mochten die Welt zu ihrem Vorteil verandern.“ (Entwurfsaspekt)

Warum ist Simulation wichtig?

Neben Experimenten mit dem realen System: Simulation einzige verfugbareTechnik zur Analyse komplexer dynamischer Prozesse

Analytische Techniken im Allgemeinen leistungsfahiger, aber durchProzesskomplexitat begrenzt

Simulation: meistgenutzte Problemlosungsstrategie

11 / 172


Grunde zum Einsatz von Simulationstechniken

Entwurf/Veranderung eines reales System; Test in unbekanntenParameterbereichen; gangige Praxis bei technischen Systemen

Experiment am Original gefahrlich, teuer, kompliziert, . . . (Bsp.: Kernreaktor,kostenintensive Wasserqualitatsuntersuchungen)

Zeitkonstanten (Eigenwerte) sind sehr klein/groß (Bsp.: explosionsartigechemische Reaktionen, Talsperrensysteme, Entwicklung von Sternen)

Vermeidung von Storungen am realen System (Bsp.:Wasserqualitatsuntersuchungen in Rohrleitungsnetzen)

Unterdruckung von Storungen am Modell, die am realen Systemunvermeidbar (Bsp.: Freilandwachstumsprozesse)

Storungen, Zustande, Systemparameter nicht zuganglich

12 / 172


Probleme bei der Anwendung der Simulation

Ironischerweise Vorteil der Simulation gleichzeitig Nachteil: allgemeine und

einfache Anwendbarkeit

Leicht: Nutzung von SimulationsprogrammenExperimentieren mit Modell hat keine KonsequenzenAber: Intelligente Nutzung der Simulation setzt voraus, dass manweiß, was man tut!

Modellverliebtheit

Systemverstandnis entsteht oft erst nach vielen, gezielten Simulationslaufen

Luenberger:”Ein kluger Gedanke spart oft Tausende von Rechnerlaufen.“

13 / 172


Kriterien zur Wahl der Simulationsmittel

Aufgabenstellung (Einschatzung von Umfang und Kompliziertheit)

Vorhandene Hard– und Software

Spezielle Kentnisse des Bearbeiters (Einarbeitung, Flexibilitat,Programmierkenntnisse)

Okonomische Kriterien (Rechner, Losungszeitraum, erwarteter Nutzen,Ersparnis teurer Experimente)

Abgrenzung

Schwerpunkt: Modellierung und digitale Simulation zeitkontinuierlicher undzeitdiskreter analoger Systeme

14 / 172


Allgemeine Vorgehensweise

Originalprozess Rechnermodell

Prozessmodell

Modelleignung Modellverifikation

ProgrammierungSystemanalyseModellbildung

Modellvalidierung

15 / 172


Arbeitsschritte

1. Erstellen des (math. oder algor.) Simulations- modells aus der Problemstellung

2. Planung der Rechenexperimente aus Problemsicht

3. Festlegen der Variablen des Simulationsmodells

4. Wahl des Simulationsmittels (Rechner + Sprache)

5. Programmierung der Aufgabe (+ syntaktische Tests)

6. Rechnung mit ausgewählten Testdaten (s. 2.)

7. Programm syntaktisch und algorithm. richtig?

8. Durchführung geplanter Rechenläufe (s. 2. u. 3.)

9. Auswertung der Ergebnisse und Überprüfen der Zielstellung

10. Programmänderungen notwendig?

11. Zusammenstellung und Interpretation der Simulationsergebnisse

A

E

Korr.

Korr.

j

n

n

j

16 / 172

Literatur I

A. Angermann, M. Beuschel, M. Rau, and U. Wohlfahrt.

MATLAB R©-Simulink R©-Stateflow R©.Oldenbourg, 2009.

D.P. Atherton and P. Borne (Eds.).

Concise encyclopedia of modelling and simulation.Pergamon Press, 1992.

P. Beater.

Regelungstechnik und Simulationstechnik mit Scilab und Modelica.Books on Demand, Norderstedt, 2010.

B.S. Bennett.

Simulation Fundamentals.Prentice Hall, 1995.

O. Beucher.

MATLAB und Simulink. Eine kursorientierte Einfuhrung.mitp, 2013.

O. Beucher.

Ubungsbuch Signale und Systeme.Springer, 2015.http://link.springer.com/book/10.1007%2F978-3-662-45955-3.

O. Beucher.

Signale und Systeme: Theorie, Simulation, Anwendung. Eine beispielhafte Einfuhrung mit MATLAB.Springer, 2015.http://link.springer.com/book/10.1007%2F978-3-662-45965-2.

17 / 172

Literatur IIA. Biran and M. Breiner.

MATLAB 5 fur Ingenieure.Addison–Wesley, 2000.

W. Bub and P. Lugner.

Systematik der Modellbildung, Teil 1: Konzeptionelle Modellbildung.VDI-Berichte 925, Modellbildung fur Regelung und Simulation, S. 1-18, VDI-Verlag, Dusseldorf, 1992.

W. Bub and P. Lugner.

Systematik der Modellbildung, Teil 2: Verifikation und Validation.VDI-Berichte 925, Modellbildung fur Regelung und Simulation, S. 19-43, VDI-Verlag, Dusseldorf, 1992.

F. E. Cellier.

Coninuous System Modeling.Springer, 1991.

F.E. Cellier and E. Kofman.

Continuous System Simulation.Springer, 2006.

DIN19226.

Leittechnik: Regelungstechnik und Steuerungstechnik.Deutsches Institut fur Normung, 1994.

P. Fritzson.

Introduction to Modeling and Simulation of Technical and Physical Systems with Modelica.IEEE Press, 2004.

P. Fritzson.

Priciples of Object-Oriented Modeling and Simulation with Modelica 2.1.IEEE Press, 2004.

18 / 172

Literatur III

P. Junglas.

Praxis der Simulationstechnik.Europa-Lehrmittel, 2014.

G. Lappus and G. Schmidt.

Simulation als vielseitiges Werkzeug - Eine Ubersicht.In K. H. Fasol and K. Diekmann, editors, Simulation in der Regelungstechnik, pages 1–29. Springer, 1990.

L Ljung and T. Glad.

Modeling of Dynamic Systems.Prentice Hall, 1994.

G. Ludyk.

CAE von Dynamischen Systemen – Analyse, Simulation und Entwurf von Regelungssystemen.Springer, 1990.

S. E. Mattsson, M. Anderson, and Astrom, K. J.

Object-Oriented Modeling and Simulation.In D. A. Linkens, editor, CAD for Control Systems, pages 31–69. Marcel Dekker, New York, 1993.

M. Otter.

Objektorientierte Modellierung Physikalischer Systeme, Teil 1.at - Automatisierungstechnik, 47(1):A1–A4 und weitere 15 Teile von OTTER, M. als Haupt– bzw. Co–Autor undanderer Autoren), 1999.

K. Reinisch.

Kybernetische Grundlagen und Beschreibung kontinuierlicher Systeme.Verlag Technik, 1974.

19 / 172

Literatur IV

F. Richert, J. Ruckert, and A. Schloßer.

Vergleich von Modelica und Matlab anhand der Modellbildung eines Dieselmotors.at, 51(6):247–254, 2003.

G. Schmidt.

Simulationstechnik.R. Oldenbourg, 1980.

W. Schweitzer.

MATLAB R© kompakt.Oldenbourg, 2013.

H. Schwetlick and H. Kretzschmar.

Numerische Verfahren fur Naturwissenschaftler und Ingenieure.Fachbuchverlag Leipzig, 1991.

U. Stein.

Programmieren mit MATLAB. Programmiersprache, Grafische Benutzeroberlachen, Anwendungen.Hanser, 2015.

Ph. Thomas.

Simulation of Industrial Processes for Control Engineers.Butterworth-Heinemann, 1999.

M. Tiller.

Intoduction to physical modeling with Modelica.Kluwer, 2001.

20 / 172

Literatur V

P. P. J. van den Bosch and van der Klauw, A.

Modeling, Identification and Simulation of Dynamical Systems.CRC Press, 1994.

A. Weinmann.

Computerunterstutzung fu Regelungsaufgaben.Springer, 1999.

J. Wernstedt.

Experimentelle Prozessanalyse.Verlag Technik, 1989.

R.L. Wood and K.L. Lawrence.

Modeling and Simulation of Dynamic Systems.Prentice Hall, 1997.

21 / 172

2. Modelle2.1 Einteilung von Modellen

Einteilung

Prozessmodell = Systemmodell + Signalmodell(e)

Im Folgenden und im Wesentlichen im Rahmen der Lehrveranstaltungbetrachtet: Systemmodell

Physikalische Modelle, Strukturmodelle (”white box“); erhalten aus Bilanz-

und Erhaltungsgleichungen in unterschiedlichen Fachgebieten (Bsp.:Differentialgleichungen (Dgln.), Ubertragungsfunktionen,Zustandsraumdarstellung, linear/nichtlinear, zeitkontinuierlich/zeitdiskret,Modellvereinfachungen, . . .)

Verhaltensmodelle (”black box“); Bsp.: Kennfelder (

”lookup table“),

Polynomansatze, kunstliche neuronale Netze

Mischformen (”grey box“)

Siehe Lehrveranstaltungen: Modellbildung, Regelungs- und Systemtechnik,Elekrotechnik, Mechanik, Physik, . . .

22 / 172

2. Modelle2.2 Zeitkontinuierliche Modelle2.2.1 Systeme mit konzentrierten Parametern

Gewohnliches Differentialgleichungssystem

Zustandsgroßen x(t) und Zeit t kontinuiertlich veranderbar

Zustand x(t) kann seinen Wert unendlich oft und zu jeder Zeit andern

→ Zustandsbeschreibung:

Gewohnliche Differentialgleichung(Dgl.)/Vektor-Dgl. (engl.: ordinarydifferential equation (ODE))

x(t) = f(x(t), u(t), t), x(0) = x0

Spezialfall: Lineare Systeme

x(t) = A x(t) + B u(t), x(0) = x0

(Anfangswertproblem;Mehrpunkt-Randwertproblem,z. B. x(tf ) = xf ) 23 / 172

2. Modelle2.2 Zeitkontinuierliche Modelle

Differential-algebraische Gleichungssysteme

Engl.: Differential-algebraic equation system (DAE system)

Implizit: f(x(t), x(t), u(t), t) = 0 , x(0) = x0

Semi-explizit: x(t) = f(x(t), u(t), t), x(0) = x0

0 = g(x(t), u(t), t)

Mit konstanter oder zeit- und/oder zustandsabhangiger, singularerMasse-Matrix M(t, x):

M(t, x) x(t) = f(x(t), u(t), t), x(0) = x0

x(t) - Zustandsvariable (differentielle V.), u(t) - algebraische Variable

Konsistente Anfangswerte!

(Differentiations-)Index: minimale Anzahl der analytischen Differentiationender algebraischen Gleichung nach der Zeit, bis explizites Dgl.-system entsteht(−→ Indexreduktion)

24 / 172


Neben Differentiationsindex: Storindex: Maß fur die Empfindlichkeit(Sensitivitat) der Losung gegenuber Stor-Term bzw. Anfangswerten

Z. B. Differentiationsindex 1:

0 =∂g

∂x

∂x

∂t︸︷︷︸f

+∂g

∂u

∂u

∂t︸︷︷︸u

+∂g

∂t

u(t) = −[∂g

∂u

]−1∂g

∂xf(x(t), u(t), t)−

[∂g

∂u

]−1∂g

∂t

Semi-expizites Dgl.-system, wenn[∂g∂u

]−1existiert

Differentiationsindex: 0 (ODE-System), 1, 2, 3, . . .; je hoher, destoschwieriger zu losen

Beim indexreduzierten System: Anfangsbedingungen der abgeleitetenGleichungen mussen auch erfullt sein

Drift: indexreduziertes System erhoht die Losungsmannigfaltigkeit; kannBerechnungsstabilitat beeintrachtigen

25 / 172


Bsp. fur differential-algebraisches Gleichungssystem: Fadenpendel(s. [1])

QF

1m =

1l =

GF

BF

φ

φy

x

26 / 172


Bsp. fur differential-algebraisches Gleichungssystem: Fadenpendel(1a)-(1e)

Variablenzuordnung: Algebraische Variable: λ(t) - Spannkraftstarke,FQ (t) = −λ(t)q = −λ(t)[x1(t) x2(t)]T , - Spannkraft,Zustandsvektor: x(t) = [x1(t) x2(t) x3(t) x4(t)]T = [x(t) y(t) vx (t) vy (t)]T ,Konstanten: g - Erdbeschleunigung, l = 1 - Lange, m = 1 - Masse

Kraftebilanz; kleine Winkel

x1(t) = x3(t) (1a)

x2(t) = x4(t) (1b)

x3(t) = −2λ(t)x1(t) (1c)

x4(t) = −2λ(t)x2(t)− g (1d)

0 = x21 (t) + x2

2 (t)− 1 (1e)

⇓ Umformung in ODE-System

λ(t) = − [2x1(t)x3(t) + 2x2(t)x4(t)]λ(t)− 3

2gx4(t)

Differentiationsindex: 327 / 172

2. Modelle2.2 Zeitkontinuierliche Modelle2.2.2 Systeme mit verteilten Parametern

Zustandsgroßen x(s, t), Ort s und Zeit t kontinuiertlich veranderbar

Zustand x(s, t) kann seinen Wert unendlich oft, zu jeder Zeit andern

Beschreibung durch partielle Dgl. (engl.: partial differential equation (PDE))

Bsp.: Diffusionsgleichung

∂x(s, t)

∂t= σ

∂2x(s, t)

∂s2

mit x - Konzentration, σ - Diffusionskonstante, s - Ort (Weg) und t - Zeit .

-1.5 -1 -0.5 0 0.5 1 1.5-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1 Color: c Vector field: -grad(c)

0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

28 / 172

2. Modelle2.3 Zeitdiskrete Modelle

Zeit t ist diskretisiert, d. h. nur zuausgewahlten Zeiten tk

bzw. diskreten Zeitpunkten(k = 0, 1, . . . ,K) ist dieBeobachtung des Systems und/oderder Eingriff in das System moglich

Ausgehend vonDifferentialgleichung: Einfuhrungeiner Diskretisierung∆t = tk+1 − tk ; Reprasentationzeitdiskreter Modelle durchDifferenzengleichungen

Zeitachse aquidistant geteilt

dx

dt≈ ∆x

∆t=

x(tk+1)− x(tk )

∆t=

xk+1 − xk

∆t= fk (xk , uk , tk ) , x0

(2)

29 / 172

2. Modelle2.3 Zeitdiskrete Modelle

Haufig ist zeitdiskretes Modell diskretisierte Variante des kontinuierlichenModelles (s. Bsp. System 1. Ordnung):

x(t) = −x(t) + u(t) , x(0) = x0

Haufiges Vorkommen in technischen und computergesteuerten Systemen

Falls Zeitachse nicht aquidistant geteilt: i. Allg. ereignis-diskreteSystembeschreibung bevorzugt

30 / 172

2. Modelle2.4 Qualitative Modelle, 2.5 Ereignis-diskrete Modelle2.4 Qualitative Modelle

Zeitachse von Natur aus diskretisiert (nicht notwendigerweise aquidistant)

Abhangige Variable (Amplitude) ist ebenfalls diskretisiert (wertdiskret)

2.5 Ereignis-diskrete Modelle

Paradoxerweise: Zeit- undAmplitude normalerweisekontinuierlich

Unterschied zu zeit- undwertkontinuierlichen Modellen:in endlicher Zeitdauer kann nureine endliche Anzahl vonZustandsveranderungenauftreten (Bsp.:Pumpenforderregimes)

Zeit- und zustandsabhangigeEreignisse

Auch gemischtereignis-diskrete/zeitkontinuierlicheModelle (hybride Modelle) moglich. 31 / 172

2. Modelle2.6 Systembeschreibung chaotischer Systeme

Gewohnliche oder partielle nichtlineare Dgln./Dgl.-Systeme

Bsp.:”Wetter“

Abbildung : Rotierende Konvektionszellen

32 / 172


Idealisierte zweidimensionale Flussigkeit, die konstantemTemperaturgradienten unterliegt

Modell beschreibt naherungsweise die Kopplung von Konvektion undWarmeleitung (Lorenz-Attraktor)

x1 = −σ (x1 − x2)

x2 = ρ x1 − x2 − x1 x3

x3 = −β x3 + x1 x2

mit σ – hydrodynamische Große, ρ ≈ Temperaturunterschied zwischen oberer

Atmosphare und Erdoberflache, β – beschreibt Geometrie des Systems, x1 ≈Geschwindigkeitsprofil mit einer konvektiven Zelle (x1 > 0: rechtsdrehend, x1 < 0:

linksdrehend), x2 – waagerechte Temperaturverteilung, x3 – senkrechte

Temperaturverteilung.

33 / 172


x1 = −σ (x1 − x2)

x2 = ρ x1 − x2 − x1 x3

x3 = −β x3 + x1 x2

Bsp.: x1 > 0 und x2 > 0: erwarmte Flussigkeit steigt auf, gekuhlte fallt

Gleichgewichtspunkte:

1 x1 = x2 = x3 = 0 keine konvektive Bewegung;festkorperahnliche Warmeableitung

2 x1 = x2 =√β (ρ− 1), x3 = ρ− 1 ; Gleichgewicht stabil fur

1 ≤ ρ < 24, aber fur ρ ≥ 24 instabil3 x1 = x2 = −

√β (ρ− 1), x3 = ρ− 1

Letzte beiden Gleichgewichtspunkte: entsprechen regelmaßigem

”Rollen“ der Konvektionszelle

Animation des Lorenz-Attraktors (Funktion lorenz unter MATLAB R©)

34 / 172

2. Modelle2.7 Stochastisches Modell – Beispiel fur Signalmodelle

Signalbeschreibung durchKenngroßen wie Erwartungswert,Streuung, (Auto-, Kreuz-)Korrelationsfunktionen, (Auto-,Kreuz-)Leistungsdichtespektren

35 / 172

2. Modelle2.8 Spektrum der Modellierung und Simulation

Es∆

Systems

BlackBox

WhiteBox

AEs

Analysis

Circuits

Economic Systems

Electric

Speculation Design

ODEs

ODEs

PsychologicalSystems

BiologicalChemicalSystems

MechanicalSystemsSystems

ControlPrediction

PDEs

Social

Zwischen mechanischen und chemischen Systemen weitere: thermische,hydraulische, pneumatische Systeme

v.r.n.l.: wohl-definiert (white box) bis”schlecht“ definiert (black box)

Modellzweck entscheidend (nach [11])

36 / 172

2. Modelle2.9 Systembegriff und Aufgabenstellungen

Problemstellungen im Zusammenhang mit Systemmodell

Direktes ProblemFall 1: E, S bekannt; A unbekannt(Vorhersage)

Inverse ProblemeFall 2: E, A bekannt; S unbekannt(Parameter-,Strukturidentifikation,Zustandsschatzung)Fall 3: S, A bekannt; E unbekannt(Steuerungsproblem)

37 / 172

2. Modelle2.9.1 Fall 1: Direktes Problem

0( , ),

( , )

x f x u x

y f x u

=

=

&

Abbildung : Vorhersage (Anfangswertproblem)

Fall 1: Vorhersage

Berechnung des zukunftigen Systemverhaltens (t ∈ [t0 = 0, t0 + K ]) unterNutzung des Systemmodells S (t0 = 0:

”Gegenwart“, K : Zeithorizont)

. . . ausgehend vom Anfangszustand x(0)

. . . bei bekannter oder angenommener Steuerung u(t) und evtl. Storung z(t)

D. h., eine Simulationsrechnung

38 / 172

2. Modelle2.9.1 Fall 1: Direktes Problem

0( , ),

( , )

x f x u x

y f x u

=

=

&

Abbildung : Vorhersage (Anfangswertproblem)

Fall 1: Vorhersage

Berechnung des zukunftigen Systemverhaltens (t ∈ [t0 = 0, t0 + K ]) unterNutzung des Systemmodells S (t0 = 0:

”Gegenwart“, K : Zeithorizont)

. . . ausgehend vom Anfangszustand x(0)

. . . bei bekannter oder angenommener Steuerung u(t) und evtl. Storung z(t)

D. h., eine Simulationsrechnung38 / 172

2. Modelle2.9.2 Fall 2: Inverses Problem

ˆmin ( , , , )

ˆ ˆ( , )

px x u p

y f x u

J

=

Abbildung : Parameteridentifikation

Fall 2: z. B. Parameteridentifikation

Berechnung der Modellparameter p unter Nutzung des Systemmodells S,gemessener Eingangs- u(t) und Zustands- x(t) bzw. Ausgangsgroßen y(t)

Typisch: Optimierungsproblem, z. B. Methode der kleinsten Fehlerquadrate

(MKQ) minp||xi − xi ||2 , i = 1, 2, . . . , nm

xi - Messwerte, xi - Modellwerte, nm - Anzahl der Stutzstellen

. . . ausgehend von einem Anfangsparameterdatensatz

. . . bei angenommener Modellstruktur

. . . iterative Losung des Problems (schrittweise Verbesserung der Parameter)

D. h., eine Simulationsrechnung pro Iteration zur Berechnung derModellwerte xi

39 / 172


0

min ( , , )

( , ),

ux u

x f x u x

J t

=&

Abbildung : Optimalsteuerungsproblem

Fall 3: Optimalsteuerungsproblem

Berechnung der optimalen Steuerung u(t) unter Nutzung des SystemmodellsS und gewunschter Zustands- x(t) bzw. Ausgangsgroßen y(t)

Optimalsteuerungsproblem

minu(t)

J(x(t), u(t))

J - Zielfunktional

. . . ausgehend von einer Anfangssteuerung u0(t)

. . . unter Nutzung des Systemmodells S

. . . iterative Losung des Optimalsteuerungs-/nichtlinearenOptimierungsproblems (schrittweise Verbesserung der Steuerung)

D. h., eine Simulationsrechnung pro Iteration zur Berechnung der Steuerung

40 / 172


0

min ( , , )

( , ),

ux u

x f x u x

J t

=&

Abbildung : Optimalsteuerungsproblem

Fall 3: Optimalsteuerungsproblem

Berechnung der optimalen Steuerung u(t) unter Nutzung des SystemmodellsS und gewunschter Zustands- x(t) bzw. Ausgangsgroßen y(t)

Optimalsteuerungsproblem

minu(t)

J(x(t), u(t))

J - Zielfunktional

. . . ausgehend von einer Anfangssteuerung u0(t)

. . . unter Nutzung des Systemmodells S

. . . iterative Losung des Optimalsteuerungs-/nichtlinearenOptimierungsproblems (schrittweise Verbesserung der Steuerung)

D. h., eine Simulationsrechnung pro Iteration zur Berechnung der Steuerung40 / 172

3. Analoge SimulationSimulation von zeit- und wertkontinuierlichen Systemen, z. B. zumSystemverstandnis, d. h. Simulation mit verschiedenen Werten derSytemparameter → Schulung des Systemdenkens

Bsp.: System 1. Ordnung

x(t) = −a x(t) + b u(t) , x(0) = x0

y(t) = c x(t)

z. B. mit u(t) = σ(t) , σ(t) =

1 fur t ≥ 00 fur t < 0

Schaltung/Signalflussbild

x0

x(t)u(t) y(t)x(t)

a

cb∫dt

1

41 / 172

3. Analoge Simulation

Rechnerisch: Auswertung der Bewegungsgleichung

x(t) = x0 e−at +

t∫0

e−a(t−τ)b u(τ)dτ

Ruckkopplung offen: reiner Integrator (System ohne Ausgleich),d. h. a = 0 . Fur b = 1 und x0 = 0.1:

x(t) = x0 +

t∫0

b dτ

= x0 + b t

= x0 + t

42 / 172


Ruckkopplung geschlossen: Verzogerungsglied (System mitAusgleich; stabil), z. B. a = 1

T1= 10, d.h. T1 = 0.1, b = 1,

x0 = 0.1 :

x(t) = x0 e−at +

t∫0

e−a(t−τ)b dτ

= x0 e−at + b e−at

t∫0

eaτdτ

= x0 e−at +b

ae−at

[eaτ

]t

0

= x0 e−at +b

ae−at (eat − 1)

= x0 e−at +b

a

(1− e−at)

43 / 172


Durch stetiges Verandern des Modellparameters a: ausgehend vom Integratorerhalt man ein T1-Glied (analoge Simulation).

Ebenso kann z. B. u(t) von σ(t) nach δ(t) verandert werden (Spungfunktion→ Impulsfunktion).

An einem Analogrechner (AR) ist wahrend dieser Anderungen gleichzeitg dieBeobachtung des Ausgangssignals moglich.

Analogrechner: siehe nachste Seite

44 / 172

3. Analoge SimulationAnalogrechner (Bildquelle: http://rechentechnik.foerderverein-tsd.de/

45 / 172


Einfache Ubertragungsglieder/Rechenelemente des Analogrechners

Summierer

ue1

uen

ue1 c1−ua

R1

R0

uen

-

Rn

+ −ua cn

1

ua = −n∑

i=1

ci uei ; ci =R0

Ri

46 / 172

3. Analoge SimulationEinfache Ubertragungsglieder/Rechenelemente des Analogrechners

Integrator

ue1

uen

R1

−ua0

ue1 c1k0 ua

C

Rnuen

ua cn+

-

1

ua(t) = −

n∑

i=1

1

CRi

∫ t

0

uei (τ)dτ + ua0

ua(t) = −

n∑

i=1

R0

Ri

1

R0C

∫ t

0

uei (τ)dτ + ua0

ua(t) = −

n∑

i=1

Ci

T

∫ t

0

uei (τ)dτ + ua0

ci =R0

Ri, T =

1

k0= CR0 , R0 - fiktiver Widerstand

47 / 172


Einfache Ubertragungsglieder/Rechenelemente des Analogrechners

Komparator

-ue1

R1

ue2

R2

ue1

x1

0

ue2

ua

A+

1

ue1 + ue2 > 0→ x = 1

ue1 + ue2 < 0→ x = 0

Koeffizientenpotentiometer, Verstarker mit einem oder mehreren Eingangen,Multiplizierer, Dividierer, Funktionsgeber, hydride Elemente (Logik, Speicher,Schalter)

48 / 172


Programmierung des Analogrechners in 4 Schritten

1 Aufstellen des Koppelplanes

2 Normierung

3 Koeffizientenvergleich

4 Entnormierung

3 Betriebsarten: Anfangswertannahme, Rechnen, Halten

Gegenuberstellung der analogen, digitalen und hybriden Simulation(Analogrechner (AR), Digitalrechner (DR), Hybridrechner (AR steuert DR;Verbindung von AR und DR uber Steuer- und Datenkoppelwerk)

49 / 172


Vor- und Nachteile der analogen Simulation im Vergleich zur digitalen Simulation

Vorteile

Berucksichtigung: Wesentliche Elemente dynamischer Systeme habenanalogen Charakter (Prozesse, Messglieder, Umformer)

”Parallelprozessor“-Eigenschaften: Rechengeschwindigkeit der Gesamtaufgabe

identisch mit der der Teilaufgaben an Einzelelementen

AR besitzt um 1 bis 2 Zehnerpotenzen hohere Rechengeschwindigkeit als DR

Interaktivitat

Anschluss von Echtzeitelementen einfach

Digitale Simulatoren haben oft Analogrechner als Vorbild

50 / 172


Vor- und Nachteile der analogen Simulation im Vergleich zur digitalen Simulation

Nachteile

Große statische und dynamische Fehler der Rechenelemente

Schaltfehler bei Betriebsartensteuerung (Anfangswertannahme, Rechnen,Halten)

Skalierung immer notwendig

Koppelplanerstellung manuell

Begrenzte Rechenkapazitat

Hoher Preis

Vorteile der digitalen Simulation: umfangreichere Palette und Moglichkeiten

(Programmierbarkeit, Softwareunterstutzung)

51 / 172

3. Analoge SimulationArbeitsschritte zur Programmierung des Analogrechners -1-

Problem Maschine1. Programmierung und Normierung (T1–Glied )

Problemgleichung Maschinengleichung(abgeleitet aus unnormierter Problemgl.)

q = − 1T

q + 1T

u , t ∈ [0, t1]1

k0

dx

dt= − c1α1x(t) + c2α2u(t) , t ∈ [0, t1]

q∫dt1

T

1T

u q

1

−xu

c1

c2k0

x

k0

α1

α2

1

52 / 172


Arbeitsschritte zur Programmierung des Analogrechners -2-

Problem Maschine2. Normierung

a) Zeit a) Zeit

τ = λt =⇒ t =τ

λτ = k0t =⇒ t =

τ

k0

τ – dimensionslose Variable τ – dimensionslose Variableλ –Normierungsfaktor k0 – Integrationsgeschwindigkeit

b) Amplitude b) Amplitudeq = Qqm ; u = Uum x = XE ; u = UE ; E = ±10 V

qm, um – bekannte oder geschatzteMaximalwerte (Probelauf)

53 / 172



Problem Maschine2. Normierung

c) Problemgleichung c) Maschinengleichungdq

dt=

dQ

d( τλ

)qm = − 1

TQqm +

1

TUum

1

k0

dx

dt=

1

k0

dX

d( τk0

)E = c1α1XE + c2α2UE

Multiplikation mit

∣∣∣∣ · 1λqm

Multiplikation mit

∣∣∣∣ · 1

E

dQ

dτ= − 1

λTQ +

1

λT

um

qmU (∗) dx

dτ= − c1α1X + c2α2U (∗∗)

54 / 172



3. Koeffizientenvergleich (*) mit (**)

c1α1 =1

λT=⇒ α1 =

1

c1λT

c2α2 =1

λT

um

qm=⇒ α2 =

1

c2λT

um

qm

Mit ciλ ist αi ∈ [0, 1] zu sichern!

55 / 172



4. Entnormierung

a) Zeit

Es muss sein: λt = k0t =⇒ t =k0

λt

t – problemabhangig festlegenk0

λ– Zeitmaßstab

k0

λ

< 1 Zeitdehnung

k0

λ=

t

t;

t

t< 1 =⇒ t > t

= 1 Echtzeit> 1 Zeitraffung =⇒ t < t

56 / 172



b) Amplitude

Es muss sein: Q = X =⇒q

qm=

X

E=⇒ q =

qm

EX

u

um=

U

E=⇒ u =

um

EU

qm

Eist Amplitudenmaßstab

57 / 172

4. Digitale Simulation

Historie

Mitte der 1960-er Jahre zugige Entwicklung blockorientierterSimulationssysteme

Motiv:

insbesondere Regelungstechniker forderten Werkzeuge fur Analyseund Entwurfaußerdem: Entwicklung der zustandsorientierten Simulation zurLosung vornehmlich der Zustandsdgl./-dzgl.

58 / 172

4. Digitale Simulation4.1 Blockorientierte Simulation

Zerlegung des zu untersuchenden Systems in geeignete Blocke; daraufaufbauend: moglichst problemnahe Notation

Bei blockorientierter Simulation:

Simulationssprache (Metasprache)Translator ubersetzt Metasprache in BasisspracheTypisch: Reihenfolge der Notation der Blocke unabhangig vonderen AbarbeitungSystem von Unterprogrammen (Funktionen, Prozeduren)Beispiele: MATLAB R©/Simulink R©, Scilab, . . .Installation auf Workstations oder PCsDigitalrechner arbeitet seriell (Einprozessormaschine oder inKernen)Zeitachse diskretisiertFunktionswerte werden nacheinander berechnet

59 / 172


Erlauterung der sequentiellen Vorgehensweise bei der digitalen Simulation an einemSystem 1. Ordnung

x = f (x , u) ≈ xk+1 − xk

h= a xk + b uk

k - diskrete Zeit , h - Schrittweite

xk+1 − xk = h (a xk + b uk )

xk+1 = (h a + 1)xk + h b uk

Mult.h b

u(k)Sprung-funktion

Mult.(1 + h a)

Summ.x(k+1) Mult.

c

x(k) y(k)Integrator

60 / 172


Fazit: Enthalt das zu simulierende System Ruckfuhrungen, so sind diezuruckgefuhrten Signale xk+1 noch nicht bekannt, wenn aus ihnen derEingang des Integrators berechnet werden soll. Es steht nur xk zur Verfugung.

Eingang des Integrators kann nur berechnet werden aus (h a + 1)xk undh b uk , nicht aber aus (h a + 1)xk+1.

Losung des Problems:

a) Reales System besitzt Totzeit Tt ; evtl. kann h = Tt gewahltwerden

b) Fehler umso kleiner, je mehr Zeitschritte (Aufwand!)c) Abtastsysteme fehlerfrei simulierbar (t = Tt)d) Iteration der Ruckfuhrung (meist benutzt), d. h. in jedem

Arbeitsschritt wird gesamtes System (einschließlichRuckfuhrungen) so lange iteriert, bis sich am Ausgang stationarerWert einstellt; dann nachster Zeitschritt

61 / 172

4. Digitale Simulation4.1 Blockorientierte SimulationKernstuck der Simulation – Numerische Integration

Allgemein: Numerische Integration zur Berechnung . . .

. . . der Flache unter einerKurve (hier wenig relevant)

Losen des bestimmten Integrals

I =b∫

a

f (t)dt (Naherung → Fehler)

0 0.5 1 1.5 2

0

20

40

60

80

100Flaeche = 7.1132

t

f(t)

a b

. . . eines Zeitverlaufes,z. B. x(t) (relevant!)

Problemstellung bei derUntersuchung dynamischerSysteme

0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

Zustandsverlauf

t

x(t)

62 / 172

4. Digitale Simulation4.1 Blockorientierte SimulationWiederholung: numerische Berechnung eines bestimmten Integrals

Aufgabenstellung: I =b∫

a

f (t)dt

Newton-Cotes-Formeln, z. B. aquidistant unterteilte Zeitachse

tj = a + j h (j = 0, 1, 2, . . . n) mit h = b−an

fur n = 1: Trapez–Regel:

I1 =h

2(f0 + f1)

fur n = 2: Simpson–Regel (Keplersche Fass–Regel):

I2 =h

3(f0 + 4 f1 + f2)

bei hinreichend glattem Funktionsverlauf arbeitet Simpson–Regelsehr genaufur n = 3: 3

8 –Regel:

I3 =3

8h (f0 + 3 f1 + 3 f2 + f3)

63 / 172

4. Digitale Simulation4.1 Blockorientierte SimulationProblemstellung bei dynamischen Systemen: Schwerpunkt im Rahmen derVeranstaltung

”Simulation“!

- Andere Aufgabenstellung im Vergleich zu vorheriger Folie -

Gegeben: Systembeschreibung x(t) = f (x, u, t) (Dgl. n–ter Ordnung bzw. n

Dgln. 1. Ordnung), p Zeitfunktionen ui (t), u = (u1 u2 . . . up)T , nAnfangswerte xi (t0)

Gesucht: (Numerische) Losung x(t) = g(t, t0, x0) (Approximation)

64 / 172

4. Digitale Simulation4.1 Blockorientierte SimulationRichtungsfeld → Losung x(t) = g(t, t0, x0)

t0 t1

x0

x1

t

x(t)

∂ x / ∂ t = f(x1,t

1)

Voraussetzungen:

Falls f (x, t) stetig ∀x ∈ X , dann geht durch (t0, x0) mindestens eineLosungskurve

Falls f (x, t) die Lipschitz-Bedingung |fk+1 − fk | ≤ K · |xk+1 − xk |, K > 0, Kkonst. ∀x ∈ X erfullt, geht durch den Punkt (t0, x0) genau eine Losungskurve

Praktisch erfullt, wenn ∂f∂x ≤ K , ∀x ∈ X

65 / 172

4. Digitale Simulation4.1 Blockorientierte SimulationGrundgedanke: Losung in diskreten Schritten

x(t0)

f(x, t)

fkx = f(x, t)

(Approximation

∫dt

x x

von xk)

x

x(t0)

xk

1

Abbildung : Zeitliche Approximation von x

66 / 172


Zeitliche Diskretisierung:

Aquidistant (zunachst betrachtet)Nicht aquidistantDurch Kollokation (Zeitdiskretisierung und Werteberechnungdurch Polynomapproximation)

Aufstellen einer numerischen Integrationsformel (Formel furZuwachs/Anderung von x)

Schreibweise: xk = x(tk ), xk+1 = x(tk+1) = x(tk + h)

Naherung durch Taylorreihe, Abbruchfehler:

x(tk + h) ≈ x(tk ) + ∆x

= x(tk ) +dx(t)

dt

∣∣∣∣∣t=tk

·h +1

2!

d2x(t)

dt2

∣∣∣∣∣t=tk

·h2 + . . .

Methode p-ter Ordnung, wenn die mit ihr berechnete Losung bis zum p-tenGlied mit der Taylorreihe ubereinstimmt

67 / 172


Abbruch der Taylorreihenach dem linearen Termliefert folgendeRekursionsformel(EULER-Verfahren):

xk+1 = xk + xk · h = xk + fk · h

Fehler:

εk+1 = |xk+1 − xk+1| =

∣∣∣∣xkh2

2!+ . . .

∣∣∣∣≈ |xk |

h2

2!=

1

2

∣∣∣fk

∣∣∣ h2 = O(h2)

Der Fehler ist von derOrdnung h2.

hopt

Schrittweite h

Fehl

er

Diskretisierungsfehler

Rundungsfehler

Gesamtfehler

Abbildung : Fehlerarten

68 / 172


Wahl der Schrittweite h

h orientiert sich an der kleinsten Zeitkonstante Tmin des Systems;Empfehlung: h = Tmin

10. . . Tmin

4

Numerische Stabilitat

nicht: Systemstabilitat

limt→∞

x(t) = xs , ∀x(t0) = x0

sondern: numerische Stabilitat

Definition:

∣∣∣∣∣δxk+1

δxk

∣∣∣∣∣ ≤ 1 ∀k

69 / 172


Beispiel:

x = −a x , x(0) = x0, a =1

T, t ≥ 0

x(t) = e−at , stabil fur a > 0

Stabilitat des EULER–Verfahrens:

xk+1 = xk + h fk = xk + h xk

xk+1 = xk − h a xk = xk (1− h a)

xk+1 = xk − h a xk = xk (1− h a)

δxk = xk − xk

δxk+1

δxk= 1− h a∣∣∣∣∣δxk+1

δxk

∣∣∣∣∣ = |1− h a| ≤ 1, → 0 ≤ h a ≤ 2

System stabil fur a > 0 und Verfahren numerisch stabil fur:

h ≤ 2

a= 2T (Einfluss der Schrittweite!)

70 / 172


Weitere Verfahren

Trapez-Verfahren (Funktionswerte bei k = 0, 1, . . . k − 1 gegeben)

xk+1 = xk +h

2(fk+1 + fk )

Adams-Bashfort-V. (ahnlich EULER-V., aber gewichtetes Mittel bzgl. fk ,fk−1)

xk+1 = xk +h

2(3fk − fk−1)

Runge-Kutta-Verfahren approximieren Taylorreihenentwicklung

Runge-Kutta-V. 2. Ordnung (p = 2, HEUN-V., verbessertesEULER-V.); verbesserte Vorhersage (Pradiktion (P)) durchEULER-Schritt; Nutzung linearer Naherung, so dass bei xk+ 1

2auch

fk

2 gilt

xPk+ 1

2= xk +

h

2f (xk , uk )

xk+1 = xk + h · f (xPk+ 1

2, uk+ 1

2)

71 / 172


Runge–Kutta–Verfahren 4. Ordnung (p = 4)

72 / 172


Erlauterungen zum Runge-Kutta-Verfahren 4. Ordnung

1. x(t) in P0(x0, t0) durch Tangente genahert;bei t0 + h ist k1 eine erste Naherung (Q1);x ≈ ∆x/h:

∆x ≈ k1 = h · tanα = hx = hf (x0, t0)

h = ti+1 − ti .

Bei R1(x0 + k12, t0 + h

2) gilt naherungsweise der

Anstieg p1 (siehe Verlauf x(t))

p1 = f (x0 + k12, t0 + h

2) .

2. Gerade P0Q2 mit p1 liefert

∆x ≈ k2 = hf (x0 + k12, t0 + h

2) .

3. Auf Gerade P0Q2 entsteht bei t0 + h2

R2(x0 + k22, t0 + h

2) mit p2 = f (x0 + k2

2, t0 + h

2) .

73 / 172


4. Gerade P0Q3 mit p2 liefert

∆x ≈ k3 = hf (x0 + k22, t0 + h

2) .

5. Die in R3(x0 + k3, t0 + h) vorgeschriebene Richtung p3

fuhrt zu P0Q4 und damit bei Q4 zu

∆x ≈ k4 = hf (x0 + k3, t0 + h) .

6. Aus k1 · · · k4 wird gewogenes Mittel gebildet∆x = 1

6(k1 + 2k2 + 2k3 + k4),

womit x(t0 + h) = x(t0) + ∆x zu berechnen ist.

Tabelle : Arbeitsschritte des Runge–Kutta-Verfahrens 4. Ordnung (sieheauch vorhergehende Abb.)

74 / 172


Charakterisierung der Integrationsverfahren

explizite/implizite (bzw. offene/geschlossene),

Einschritt-/Mehrschritt-,

Pradiktor-/Korrektor-,

rekursive/iterative Verfahren

explizit (offen) implizit (geschlossen)

EULER–Verfahren (-vorwarts) EULER-V. (-ruckwarts)- Vorwartsdifferenzenapproximation - - Ruckwartsdifferenzenapproximation -xk+1 = f (xk , xk ) xk+1 = f (xk+1, xk+1, uk+1)

(nur losbar mit Iterationsverfahren,z.B. Sukzessive Approximation)

xk+1 = xk + h xk xk+1 = xk + h xk+1

Trapez–Verfahrenxk+1 = xk + h

2(xk+1 + xk )

75 / 172


Charakterisierung der Integrationsverfahren - Fortsetzung

explizit (offen) implizit (geschlossen)

Runge-Kutta-V. 4. Ordnung Heun-V. (verbessertes EULER-V.)xk+1 = xk + h

6(k1 + 2k2 + 2k3 + k4)

k1 = f (xk , uk )k2 = f (xk + h

2k1, uk+ 1

2)

k3 = f (xk + h2

k2, uk+ 12)

k4 = f (xk + h k3, uk+1)

xk+1 = f (xk+1, uk+1)≈ f (xk + h xk , uk+1)

Vorhersage:xk = f (xk , uk )

2 Pradiktorschritte:xP

k+1 = xk + h xk

xPk+1 = f (xP

k+1, uk+1)Korrektorschritt:

xk+1 = xk + h2

(xPk+1 + xk )

Einschritt-V. Mehrschritt-V.

alle Argumente aus [tk , tk+1], Verwendung zeitlich weiter zuruck-max. einen Zeitschritt

”alt“ liegender Stutzwerte

in [tk−m, tk−m+1, . . . tk , tk+1]

Mittelpunktregel:xk+1−xk−1

2h≈ xk

76 / 172


Schrittweitensteuerung

bisher: feste Integrationsschrittweite

Schrittweitensteuerung (nicht aquidistante zeitliche Diskretisierung)

h xIntegrations-algorithmus

Schrittweiten-steuerung

abschatzungFehler-

κsoll

κist

1

hneu = haltκsoll

κist

77 / 172


Steife Differentialgleichungssysteme

Verhaltnis von großter zukleinster Zeitkonstante

Tmax

Tmin≈ (20) 100 . . . 1000

oder noch großer

Bsp.: Zustandsvektor-Dgl.

x =

[−1 00 −20

]x , x(0) =

[11

]

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

0.2

0.4

0.6

0.8

1

Zustandsverlaeufe

t

x(t)

Tmax

/Tmin

=20

78 / 172


Numerische Stabilitat am vorangegangenen Beispiel

EULER–vorwarts: h1 < 0.1 , h2 < 2

EULER–ruckwarts: stabil ∀h

→ Verwendung impliziter Verfahren (z.B. Heun-Gear-Verfahren);

”Problem“: Anlauf, Pradiktion

Integration bei unstetiger rechter Seite der Dgl.

z.B. Sprungsignal, Begrenzung, unstetige Regler

→ Einschrittverfahren mit automatischer Schrittweitensteuerung (h wird sogewahlt, dass Unstetigkeit an der Intervallgrenze tk , tk+1 , . . . liegt)

79 / 172


Integration von differential-algebraischen Systemen(DAE-Systemen)

Sehr hohe Steifigkeit: DAE-System

DAE-System: algebraische Schleife

Verfahren mit Ruckwarts-Differenzen-Formeln (engl.: backward differentiationformulas (BDF)), d. h. implizite Mehrschritt-Verfahren oderKollokationsverfahren

80 / 172


Algebraische Schleifen

ua

ye

1

e = u − y

y = a e

t = 0:

e = y = 0, u = 1

fur t > 0:

e = 1− 0 = 1

y = a · 1 = a

e = 1− a

y = a(1− a) = a− a2

e = 1− (a− a2)

y = a(1− (a− a2)) = a− a2 + a3 , usw.

t >> : y = a− a2 + a3 − a4 + a5 − . . .+ . . .

81 / 172


Algebraische Schleifen - Fortsetzung

Konvergenz: nur fur −1 < a < 1

y =a

1 + a

Ausweg: zeitdiskrete Realisierung

e(k) = u(k)− y(k − 1)

y(k) = a e(k)

Algebraische Schleifen konnen, wenn sie nicht entdeckt werden, zuInstabilitaten bei einem stabilen Modell fuhren.

Auswege beim Auftreten von algebraischen Schleifen:

Vermeiden (so weit wie moglich; nicht immer moglich)Schleifenverstarkung < 1: Tt (Modellveranderung; akzeptabel)Iterative Losung

82 / 172


Algebraische Schleifen - Fortsetzung

Im Allgemeinen fuhrt Auftreten von algebraischen Schleifen (neben anderenSystemeigenschaften) zu einem DAE-System:

x(t) = f (x, u, t)

0 = g(x, u, t)

Losung mittels BDF- oder Kollokationsverfahren

Stationare Punkte

Bestimmung der Ruhelage/des Gleichgewichtszustandes

x(t) = f (x, u, t), x(0) = x0

Losung der algebraischen Gleichung (Nullstellenbestimmung)

0 = f (xs , u), u = us = const.,Anfangswert: x0S = x0

Losungsverfahren (z.B. Newton–Raphson–Verfahren, Picard–Verfahren)

83 / 172


Empfehlungen zum Einsatz numerischer Integrationsverfahren

1. Versuch: Runge-Kutta-V. 4./5. Ordnung mit Schrittweitensteuerung(Dormand-Prince-V., explizites Einschritt-V.)

Bei geringen Genauigkeitsanforderungen: Runge-Kutta-V. 2./3. Ordnung(Bogacki-Shampine-V., explizites Einschrittverfahren; auch geeignet beimoderater Steifigkeit)

Bei hoheren Genauigkeitsanforderungen: Adams-Bashfort-Moulton-V.(explizites Mehrschrittverfahren)

Fur steife Differentialgleichungssysteme:

Gear-V. bzw. Heun-Gear-V. variabler Ordnung (implizitesMehrschrittverfahren; benutzt numerische(Ruckwarts-)Differenzen)implizites Runge-Kutta-V. 2. Ordnung (implizitesEinschrittverfahren; benutzt Trapezregel und Ruckwartsdifferenz)

Fur DAE-Systeme: BDF- oder Kollokationsverfahren

84 / 172


Numerische Differentiation - Pendant zur Integration

Aufbauend auf Newtonschen Interpolationspolynom:

dx

dt

∣∣∣∣∣t=tk

=xk − xk−1

hSchrittweite h = tk − tk−1

Genauere Ergebnisse liefern Ausdrucke, die symmetrisch zu tk liegen, z.B.:

dx

dt

∣∣∣∣∣t=tk

=xk+1 − xk−1

2h

oder fur die 2. Ableitung

d2x

dt2

∣∣∣∣∣t=tk

=1

h2(xk+1 − 2xk + xk−1)

85 / 172

4. Digitale Simulation4.2 Zustandsorientierte Simulation

Leistungsfahige Modellierungsmethoden + Fortschritte in der Rechentechnik→ verbesserte Anwendbarkeit von Methoden zum optimalen Systementwurf(optimale Steuerung, adaptive on–line Steuerung)

Bei diesen Problemstellungen:

Objektsimulation ist kleiner TeilNicht zweckmaßig: Anwendung von SimulationssystemenBesser: numerische Integration von Differentialgleichungssystemenoder Auswertung der Bewegungsgleichung fur lineare bzw.quasi–lineare Systeme

Zur Erinnerung: Zustandsbeschreibung:

x(t) = A x(t) + B u(t)

y(t) = C x(t) + D u(t)

Bei zeitvarianten Systemen sind die Systemmatrizen zeitabhangig,z.B. A = A(t) .

86 / 172

4. Digitale Simulation4.2 Zustandsorientierte Simulation

Allgemeine Losung der Zustands-Dgl.:

x(t) = ΦΦΦ(t)x(0) +

t∫0

ΦΦΦ(t − τ)Bu(τ)dτ

ΦΦΦ(t) = eAt

Numerische Losung der Zustands-Dgl. mittels numerischer Integration undrekursiver Beziehung xk+1 = f (xk )

Weiterfuhrende Literatur, z.B. Ludyk [19] oder Reinisch [22]

87 / 172

4. Digitale Simulation4.3 Objektorientierte Simulation

Objektorientierte Modellierung: Verallgemeinerung der blockorientiertenModellierung

Geeignet fur multidisziplinare Systemmodellierung und -simulation(Simulationssysteme oft auf Fachgebiete zugeschnitten; Modellierung vonKomponenten oder Disziplinen anderer Fachgebiete nicht oder nur ineingeschranktem Maße moglich)

Praxisanwendungen: Automobilindustrie, Luft- und Raumfahrt, Robotik,Thermodynamik (Gase, Flussigkeiten), Klimatisierung, Kraftwerke, Mechanik(Mehrkorpersysteme), Optik, Abwasserreinigung, . . .

Bisher: meist blockorientierte (kausale) Modellierung vorherrschend

88 / 172


Kausale Modellierung - Vorteile/Nutzlichkeit

Systemverstandnis (Pole, Nullstellen, Ubertragungsfunktion, -verhalten,lineares/nichtlineares/linearisiertes System, . . .

Beschreibung der Steuerungs-/Regelungsstruktur

Kausale Modellierung - Nachteile

Aufwandig zu erzeugen/umzuformen; manuelle Ableitung aus Zustands- undErhaltungsgleichungen; fehlertrachtig

Weniger gut wiederverwendbar

Schwierigkeiten bei komplexen Systemen, konsistente Anfangsbedinungenherzustellen (bei differential-algebraischen Gleichungssystemen(DAE-Systemen))

89 / 172


Kausale Modellierung - Nachteile (Fortsetzung)

A priori (Im Vorhinein): Was ist bekannt (Eingang), was ist unbekannt(Ausgang)?

Annahmen uber Kausalitat auf Komponentenebene (nicht: Gesamtsystem)

Robustheit, Wartung: Modellparameter”verteilt“ (nicht im Sinne von

Systemen mit verteilten Parametern, sondern im Sinne”verstreut“ im

gesamten Modell (Signalflussbild))

Diagramm beinhaltet nicht genugend Informationen, umAnfangsbedingungen zu berechnen

90 / 172


Akausale Modellierung (physikalisch, objektorientiert) - Merkmale

Beschreibung eines Systems/einer Komponente anhand von Zustands- undErhaltungsgleichungen (i. d. R. Kombination; −→ DAE-System);physikalische Gesetzmaßgikeiten

Keine explizite Festlegung uber Ein- und Ausgang

Zwei Arten von Variablen: Potenzial- und Flussvariablen; (Potenzialvariablen:Triebkrafte); Einheit des Produkts von Fluss- und Potenzialvariablentypischerweise Einheit der Leistung

Objektdiagramm: hierarchischer Aufbau, Verallgemeinerung desBlockschaltbilds

Beispiele (s. auch nachste beide Folien):

Block- und objektorientierte Modellierung eines elektrischenNetzwerks (Blockschaltbild, Objektdiagramm)Objektdiagramm eines AntriebsstrangsObjektdiagramm eines Motors

91 / 172


Blockorientierte vs. physikalische objektorientierte Modellierung

Blockorientiert (Ursache → Wirkung, Ein-/Ausgang, kausal)

Objektorientiert (physikalisch, akausal)

92 / 172


Verbindungen

Signalflusse

Starre mechanische

Regler Getriebe LastMotor

1

Abbildung : Objektdiagramm eines Antriebsstrangs

93 / 172


richterStrom-

Motor-tragheit

Strom

EMK

Drehzahl

FlanschMechanischer

R L

1

Abbildung : Objektdiagramm eines Motors

Bestandteile des Objektdiagramms

Grafische Darstellung der physikalischen Komponenten

Schnittstellen, mit denen die Komponente mit anderen Bauteilen verbundenwerden kann

94 / 172


Bestandteile des Objektdiagramms (Fortsetzung)

Gerichtete oder ungerichtete Verbindungslinien (elektrische, hydraulischeLeitungen, mechanisch starre Verbindungen, aber auch Signalflusse)

Komponenten werden unabhangig von Umgebung definiert (lokale undSchnittstellenvariablen)

Komponente ist hierarchisch aufgebaut oder durch algebraische bzw. Dgln.beschrieben

Arten von Verbindungsgleichungen zwischen Komponenten

Verbundene Variablen haben denselben Wert: Potenzialvariablen, Bsp.:Spannung, Weg, Geschwindigkeit, Druck

Die Summe der verbundenen Variablen verschwindet: Flussvariablen, Bsp.:Strom, Kraft, Moment, Volumen–, Massenstrom; an allen Elementen gleichepositive Flussrichtung, z.B. in Element hinein

95 / 172


Beispiele fur elektrische Einzelkomponenten (Objektdiagramm undBeschreibungsgleichungen)

Widerstand

v1 v2

i2i1

u

R+

1

0 = i1 + i2

u = v1 − v2

u = R i1

Kapazitat

v1 v2

i2i1 C

u

+

1

0 = i1 + i2

u = v1 − v2

i1 = C u

Elektrische Zweipole haben Gemeinsamkeiten (Stromsumme = Null,Spannung = Potentialdifferenz) ⇒ Definition einer (noch nicht vollstandigen)Klasse (eines Teil-Modells) Zweipol

96 / 172


Beispiele fur elektrische Einzelkomponenten (Objektdiagramm undBeschreibungsgleichungen) - Fortsetzung

Induktivitat

v1 v2

i2i1

u

L+

1

0 = i1 + i2

u = v1 − v2

u1 = Ldi1dt

Spannungsquelle

∼v2

i2i1

v1

u

u(t)

+

1

0 = i1 + i2

u = v1 − v2

u = A sinωt

97 / 172


Beispiele furelektrische Einzel-komponenten(Objektdiagrammund Beschrei-bungsgleichungen)- Fortsetzung

Erdung

v i

1

v = 0

Vorgehensweise bei der Modellierung

a)”buttom-up“-Entwurf

Definition einer (noch nichtvollstandigen) Klasse (einesTeil-Modells) Zweipol mitallgemeinen Eigenschaften sowiePotential- und FlussvariablenDefinition einer Klasse(z. B. Widerstand), dass dieEigenschaften von Zweipol erbtInstanziierung einer Klasse (=Objekt) (z. B. R vom TypWiderstand mit dem Wert 1 Ω)Verbindung zum Gesamtmodell

b)”top-down“-Entwurf (umgekehrt zu a) )

98 / 172


Zusammenschaltung o. g. Komponenten

ElektrischerSchaltkreis(Objektdiagramm)

u(t)

g

3

1

2

R

C+

+

A

+

∼

1

Anmerkungen

17 Gleichungen, 17 Unbekannte;R.i1 R.i2 R.v1 R.v2 R.uC .i1 C .i2 C .v1 C .v2 C du

dt

A.i1 A.i2 A.v1 A.v2 A.ug .i g .v

1 Zustand (C .u); Beschreibung lauftauf Zustandsbeschreibung hinaus, d.h.x(t) = f (x , u, t); Anfangszustandbekannt, Parameter bekannt

99 / 172


Ablauf der objektorientierten Simulation

Modellierungs-system

Objektor.Transformations-Symbolische

algorithmen

numerischeUbliche

Integrations-verfahren

Losung

diagramm

Zustands-Objekt-

form

besetzteMatrizen)

DAEs

(groß,schwach

1

100 / 172


Anmerkungen

Manuelle Umformung aufwandig (nach der interessierendenZustandsableitung); → Algorithmen sortieren DAEs (rekursiveBerechnungsvorschrift, nicht vereinfacht; symbolische Verarbeitung)

Vorgehen systematisch und einfach

Rechnerunterstutzung unabdingbar

Seit den 1990-er Jahren: Entwicklung von Softwaresystemen und Sprachen(unvollstandige Liste: ABACUS, Omola, Dymola R©, Modelica R©,OpenModelica, SimulationX, Wolfram System Modeler R©, MapleSimTM . . . )

Weiterfuhrende Literatur: Otter, et. al. [21], Fritzson [15]

101 / 172

5. Simulationssprachen und -systeme5.1 Ubersicht – kleine Auswahl

Bsp.: Softwarebibliotheken

The Numerical Algorithm Group (NAG) / Kommerzielle numerischeAlgorithmen / nag.co.uk

Netlib / Freie mathematische Software (Numerik, wissenschaftlichesRechnen) / netlib.org

Bsp.: Simulationssysteme mit grafischer Oberflache

LabVIEW / Grafische Programmierumgebung (Messdatenerfassung und-verarbeitung, Messgeratesteuerung, eingebettete Systeme, Steuerung,Regelung, . . .) / ni.com/labview/d/

102 / 172

nag.co.uk

netlib.org

ni.com/labview/d/

5.1 Ubersicht – kleine Auswahl

Bsp.: Integrierte Systeme (Modellierungssprache undSimulationsumgebung)

OpenModelica, Dymola / Objektorientiert; Modellbildung und Simulation /openmodelica.org, http://www.3ds.com/products-services/catia/products/dymola

MATHMODELICA / Objektorientiert; Modellbildung und Simulation /http://www.mathcore.com/products/mathmodelica/

MATLAB R©/Simulink R© / MATLAB R©: fuhrende Sprache furwiss.-techn. Berechnungen, Steuerungsentwurf, . . . ; Modellierung,Simulation, Analyse dynamischer Systeme, usw. / mathworks.com

GNU Octave / Freie Software; numerische Berechnungen, . . .; uberwiegendMATLAB R©-kompatibel / octave.org

Scilab / Freie Software; numerische Berechnungen, ahnlich MATLAB R© . . . /scilab.org

Bsp.: Sprachen

MATLAB R© / Hochsprache fur wiss.-techn. Berechnungen / mathworks.com

Modelica R© / Freie Software; Objektorientierte Sprache / modelica.org

103 / 172

openmodelica.org

http://www.3ds.com/products-services/catia/products/dymola

http://www.mathcore.com/products/mathmodelica/

mathworks.com

octave.org

scilab.org

mathworks.com

modelica.org


Bsp.: Simulation industrieller Prozesse

gPROMS / Modellierung, Simulation und Optimierung kont. dyn. Systeme;Simulator einzubetten in Kundenanwendungen /http://psenterprise.com/gproms.html

WITNESS / Modellierung, Simulation und Optimierung derArbeitsumgebung einschl. Geschaftsprozesse / lanner.com/en/witness.cfm

Bsp.: Web-basierte und verteilte Simulation

RT-LAB / Echzeit-Modellierungs- und Steuerungswerkzeug (HiL) /opal-rt.com

JSIM / Java-basiert; Schwerpunkt: Physiologie, Biomedizin /http://www.physiome.org/jsim/

104 / 172

http://psenterprise.com/gproms.html

lanner.com/en/witness.cfm

opal-rt.com

http://www.physiome.org/jsim/


Bsp.: Sonstige

Hybride Systeme

Netzwerke

Anwendungsspezifische Systeme (Welt, Umwelt, Wirtschaft, Produktion,CAD, Ressourcen, u.v.a.m.)

Ubersicht

en.wikipedia.org/wiki/List_of_computer_simulation_software

105 / 172

en.wikipedia.org/wiki/List_of_computer_simulation_software

5.2 MATLAB R©/Simulink R© - Einfuhrung5.2.1 Einleitung

MATLAB R©: Hochsprache fur wissenschaftlich-technische Berechnungen

(Berechnungen, Visualisierung, Programmierung) in leicht handhabbarer

Umgebung

Mathematisch-numerische Berechnungen (auch symbolische)Modellierung, Simulation, Prototyping, Algorithmen, . . .Datenanalyse, -auswertung und -visualisierungWissenschaftlich-technische grafische DarstellungenAnwenderprogrammentwicklung (einschl. grafischesBenutzerinterface)Physikalische objekt-orientierte Modellierung

Simulink R©: Block- (tw. objekt-)orientiertes grafisches Simulationssystem

Fachgebietsspezifische Toolboxen

Abkurzung MATLAB R©: matrix laboratory

Grundlegendes Datenelement: Feld (array, matrix)Interaktiv (Kommando–Interpreter)Ursprungliche Entwicklung: leichter Zugang zu LINPACK undEISPACK (zwei Bibliotheken,

”state-of-the-art“ der

Matrizenrechnung) und grafische Veranschaulichung106 / 172

5.2 MATLAB R©/Simulink R© - Einfuhrung

Uberblick (Quelle: www.de.mathworks.com/products)

107 / 172

www.de.mathworks.com/products

5.2 MATLAB R©/Simulink R© - Einfuhrung5.2.1 Einleitung

Entwicklung von MATLAB R© im Verlaufe der Zeit: Toolboxen,Objektorientierung (Klassendefinitionen, . . .)

Verfugbarkeit: PCs und Workstations, verschiedene Betriebssysteme(WINDOWS, Unix, Linux, . . .), TU Ilmenau: im Uni-Rechnernetz (ThuringerLandeslizenz der Unis; zeitparallele Zugriffe)

Kommerzielles Produkt; Studentenversion (82,00 e)(https://www.academic-center.de/cgi-bin/program/S1672) mit(kleinen) Einschrankungen (MATLAB R©, Simulink R© (1000 Blocke), 12Toolboxen; auch: Prototyping, Test auf Arduino, LEGO MINDSTORMSNXT, Raspberry Pi)

Informationen

http://www.mathworks.com

http://www.mathworks.de

Literaturliste (Bibliothek), Copy-Shop, Internet

Alternativen zu MATLAB R©:

Scilab/Xcos unter Linux, Windows, usw. (Freie, quelloffeneSoftware, GPL-kompatibel), https://www.scilab.orgGNU Octave unter Linux, Windows, usw. (Freie, quelloffeneSoftware), https://www.gnu.org/software/octave 108 / 172

https://www.academic-center.de/cgi-bin/program/S1672

http://www.mathworks.com

http://www.mathworks.de

https://www.scilab.org

https://www.gnu.org/software/octave

5.2 MATLAB R©/Simulink R© - Einfuhrung5.2.2 Arbeitsumgebung

Kommando-, Workspace-, Verzeichnis-, Historiefenster

Reiter Home, Plots, Apps

Kommandomanipulation; -folgen in ASCII-Skript-Dateien; (*.m)

Hilfe, Informationen, Ubersichten, usw.

Befehl Erlauterung

demo Beispiele

doc schlusselwort Hypertext-Dokumentation (auch PDF-Format)

help, help help,help fcn | m-file | thema Textuelle Online-Hilfe

exist Test, ob Variable oder Funktione definiert

type Auflisten vom M-Files

edit Editieren von M-Files

109 / 172

5.2 MATLAB R©/Simulink R© - Einfuhrung5.2.2 Arbeitsumgebung

Wichtige”Betriebssystemkommandos“

cd Verzeichniswechsel

pwd Ausgabe des gegenwartigen Arbeitsverzeichnisses

ls, dir Inhaltsverzeichnis

delete Loschen von Dateien

!befehl Betriebssystemskommandos ((WIN)DO(W)S, LINUX, . . .)

Sonstiges

format Zahlenformat (short, long, e, g, bank, . . .)

; / ... Ausgabeunterdruckung/Kommandofortsetzung

Workspace:who, whos Info zu Variablen (Kurz-, Langform); s. auch Workspace

clear Loschen von Variablen oder Funktionen

save, load Speichern bzw. Laden von Variablen in bzw.

von Datei mit angebbarem Format (-mat, -ascii, u.a.)

diary Tagebuchfunktion, d. h. Speichern aller Kommandos undErgebnisausgaben des Kommandofensters in einer Datei

110 / 172

5.2 MATLAB R©/Simulink R© - Einfuhrung5.2.3 MATLAB R©-Kommandosprache

Variable und Datentypen

Variable: erste 31 Zeichen signifikant; Buchstaben, Ziffern, ; beginnend mitBuchstaben; keine explizite Vereinbarung; Speicherung im Workspace(Arbeitsspeicher)

Spezielle Variable und Variablenwerte

Variable Erlauterung

ans Antwortvariable, wenn das Ergebnis keiner Variablenzugewiesen wurde

eps Relative Genauigkeit

realmax, realmin Großte bzw. kleinste darstellbare Gleitkommazahl

pi π = 3.1415926 . . .

i,j Imaginare Einheit

inf Unendlich (∞)

NaN”not a number“ , Ergebnis eines unbestimmten

Ausdruckes, z. B. 00

111 / 172


Datentypen (Klassen) – eine Auswahl

Datentyp Erlauterung

double Vektoren, Matrizen, 3D-Felder;Elemente reelle oder komplexe Zahlen; voll besetzt;Indizierung bei 3-D-F.: (Zeile, Spalte, Seite)

sparse Schwach besetzte Matrizen (sparse matrices)

char Zeichenketten(-felder)

cell Mehrdimensionale”Zellenfelder“; Zellen: Elemente

unterschiedlicher Datentypen, auch Strukturen

struct Elemente unterschiedlicher Datentypen

einer Variablen zusammengefaßt

int8, ..., int64, Vorzeichenbehaftete unduint16, ..., uint64, . . . vorzeichenlose ganze Zahlenfelder

logical Boolesche Felder

. . . . . .

112 / 172

5.2 MATLAB R©/Simulink R© - Einfuhrung5.2.3 MATLAB R©-KommandospracheZahlen, Zeichen(ketten)

Vorzeichenbehaftete reelle oder komplexe Zahlen in Festkomma- oder inwissenschaftlicher Darstellung (Mantisse, Exponent):

Bsp.: 3 9.6 -.99 1e3 -5.34e-4 -3.14j 1i 3-4i

Zeichen(ketten) in Hochkomma eingeschlossen, Bsp. ’xyz’ ,’er hat’’s’

Operatoren (doc ops)

Arithmetische, Vergleichs–, logische, bitweise, Mengen- und spezielleOperatoren

Funktionen (doc elfun)

Mathematische Funktionen (Winkelfunktionen, Exponentialfunktion,Logarithmen, Wurzel, Rundung, Teiler, Rest, Komplexe Zahlen, ...) (docelfun) u.v.a., z.B. Umwandlung von Datentypen (Zeichenketten,numerisch), Zeichenkettenmanipulation

Ausdrucke

Ausdrucke werden aus Variablen, Zahlen, Zeichenketten, Operatoren undFunktionen gebildet.

113 / 172


Ein-/Ausgabe

save Speicherung aller Workspace-Variablenin einer binaren oder ASCII-Datei

Bsp.:save fname Speicherung der Datei fname.matsave matlab.mat

save fname X,Y,Z fname.mat; Inhalt: X, Y und Zsave fname -ascii fname.mat im ASCII-Format

load Laden von Variablen aus einer binarenoder ASCII-Datei in den Workspace

Bsp.:load fname Laden aus dem binaren MAT-File fname.matload dto. matlab.matload test.dat Laden aus der ASCII-Datei test.dat bzw.(oder load test -ascii) test; Abspeichern in Variable/Matrix test

114 / 172


Ein-/Ausgabe (Fortsetzung)

input Tastatureingabe durch NutzerBsp.: n=input(’Wieviele Messwerte ?’);

Eingabe: MATLAB R©-AusdruckBsp.: Name=input(’Wie ist Ihr Name? : ’,’s’);

Eingabe: Zeichenkette; wird nicht ausgewertet

fopen, fclose,fread, fwrite,fscanf, fprintf . . .

”low–level“-E/A-Kommandos

print Ausgabe von Grafiken auf Drucker, Plotter oder inDatei; Formate (PS, EPS, JPG, PNG, TIF, HPGL,druckerspezifisch); Zusammenwirken mit GhostScript

plot(x,y) Grafik mit linearer Achseneinteilung, y = f (x)

Weitere Kommandos zur grafischen Darstellung: siehe Visualisierung

Weitere Sprachelemente: siehe M-File-Programmierung

115 / 172

5.2 MATLAB R©/Simulink R© - Einfuhrung5.2.4 Matrizen, lineare Algebra

Matrizen, Teilmatrizen (doc elmat)

Matrix A =

a11 a12 a13

a21 a22 a23

a31 a32 a33

=

1 2 34 5 67 8 6

A=[1 2 3;4 5 6;7 8 6]

Teilmatrix B =

[a11 a12

a21 a22

]B=A(1:2,1:2)

Generierung von Matrizen

Funktion Erlauterung

zeros, ones, eye Null-, Eins-, Einheitsmatrizen

rand, randn zufallig belegte Matrizen (Gleich-, Normalverteilung)

linspace, logspace Vektor mit aquidistanter/logarithmischer Teilung

diag Diagonalmatrizen/Diagonalen

116 / 172


Matrizenoperationen

Operator Erlauterung

’ Transposition

+, - Addition, Subtraktion

* Matrizenmultiplikation

.* Elementweise Multiplikation

\, /, ./ Links-, Rechts- und elementweise Division

ˆ, .ˆ Potenzierung, elementweise Potenzierung

117 / 172


Lineare Gleichungssysteme (m = n)

1. Bsp.: m, n - Anzahl d. Gleichungen bzw. Anzahl d. Variablen

3 x1 + 2 x2 − 3 x3 = 2

−5 x1 − 4 x2 + 6 x3 = −1

2 x1 − x3 = 1

Matrix-/Vektor-Schreibweise: 3 2 −3−5 −4 6

2 0 −1

x1

x2

x3

=

2−1

1

A x = b

x = A−1 b

Losung: x=inv(A)*b oder x=A\b

118 / 172


Lineare Gleichungssysteme (m = n)

1. Bsp.: m, n - Anzahl d. Gleichungen bzw. Anzahl d. Variablen

3 x1 + 2 x2 − 3 x3 = 2

−5 x1 − 4 x2 + 6 x3 = −1

2 x1 − x3 = 1

Matrix-/Vektor-Schreibweise: 3 2 −3−5 −4 6

2 0 −1

x1

x2

x3

=

2−1

1

A x = b

x = A−1 b

Losung: x=inv(A)*b oder x=A\b 118 / 172


Lineare Gleichungssysteme (uberbestimmt; m > n)

2. Bsp.:

3 x1 + 2 x2 = 3

x1 − 2 x2 = 1

x1 − x2 = 1

4 x1 + x2 = 4

Matrix-/Vektor-Schreibweise:3 21 −21 −14 1

[ x1

x2

]=

3114

A x = b

Losung (i. S. kleinste Fehlerquadrate): hier nur moglich: x=A\b

119 / 172


Lineare Gleichungssysteme (uberbestimmt; m > n)

2. Bsp.:

3 x1 + 2 x2 = 3

x1 − 2 x2 = 1

x1 − x2 = 1

4 x1 + x2 = 4

Matrix-/Vektor-Schreibweise:3 21 −21 −14 1

[ x1

x2

]=

3114

A x = b

Losung (i. S. kleinste Fehlerquadrate): hier nur moglich: x=A\b 119 / 172


Lineare Gleichungssysteme (unterbestimmt; m < n)

3. Bsp.: keine eindeutige Losung

6 x1 + 8 x2 + 7 x3 + 3 x4 = 7

3 x1 + 5 x2 + 4 x3 + x4 = 8

Losung: siehe MATLAB-Online-Hilfe (Linear Equations/Systems of LinearEquations/Underdetermined Systems (Stichworte: Basislosungen; hochstensm Nicht-Null-Komponenten; Linearkombination; Nullraum-Vektor))

120 / 172


Matrizenfunktionen (doc matfun)

Befehl Erlauterung

norm, cond Matrix- oder Vektornorm, Konditionszahl

rank, inv Rang, Inverse einer Matrix

det Determinante einer Matrix

trace Spur einer Matrix (Summe der Diagonalelemente)

chol, lu Cholesky-, LU-Faktorisierung

eig Eigenwerte, -vektoren

poly Koeffizienten des charakteristischen Polynoms

. . .

121 / 172

5.2 MATLAB R©/Simulink R© - Einfuhrung5.2.5 Gewohnliche Differentialgleichungen

- Anfangswertprobleme -

Funktionen fur Dgln. (doc funfun)

Befehl/Datei Erlauterung

odeset, odeget Einstellen/Lesen von Parametern zur Dgl.losung

Num. Integrationsverf. zur Losung von Dgl.systemen . . .ode45, ode23, ode113 . . . nicht-steifeode15s, ode23s, ...

ode23t, ode23tb . . . steifeode15s, (ode23t) . . . DAEs (mit konst. Masse-Matrix)ode15i . . . implizite

odeplot Zeitverlaufgrafik wahrend der Simulation

odephas2, odephas3 2D-, 3D-Phasendiagramm wahrendder Simulation

odeprint Tabellarische Ergebnisse im Kommandofenster

122 / 172


Beispiel 1: Logistisches Wachstum (Wachstums-Dgl.)

y(t) = a y(t)(1− y(t)) , y(0) = y0 = 0.1

mit Wachstumsparameter a = 3.2

MATLAB R©-Funktion (Datei: growth.m)

function ypunkt=growth(t,y)

a=3.2; % Wachstumsparameter

ypunkt=a*y*(1-y);

Aufruf (im Kommandofenster, im Skript, in einer Funktion)

[t,y]=ode45(@growth,[0 5],0.1);plot(t,y) % odersol=ode45(@growth,[0 5],0.1);plot(sol.x,sol.y) % oder—–sol=ode45(@growth,[0 5],0.1);x=linspace(0,5,101); % Zeitbasis neuy=deval(sol,x,1); % Interpolation fuer neue Zeitbasis

plot(x,y)

123 / 172


Beispiel 2: Van-der-Pol-Schwinger (Dgl.-system)

y(t) =

[y1(t)y2(t)

]=

[y2(t)

µ(1− y 2

1 (t))

y2(t)− y1(t)

], y(0) = y0 =

[20

]µ = 1, 1000

MATLAB R©-Funktion (Datei: shvdpmue.m)

function dydt=shvdpmue(t,y,mue) % mit weiterer Parameteruebergabe

dydt=[y(2);mue*(1-y(1)ˆ2)*y(2)-y(1)];


options=odeset(’OutputFcn’,’odeplot’); mue=1;

[t,y]=ode45(@shvdpmue,[0 50],[2;0],options,mue); % oder

[t,y]=ode15s(@shvdpmue,[0 50],[2;0],options,mue);

Eigenschaft fur µ = 1, 1000?

µ = 1: nicht-steif, µ = 1000: steif

124 / 172


Beispiel 2: Van-der-Pol-Schwinger (Dgl.-system)

y(t) =

[y1(t)y2(t)

]=

[y2(t)

µ(1− y 2

1 (t))

y2(t)− y1(t)

], y(0) = y0 =

[20

]µ = 1, 1000

MATLAB R©-Funktion (Datei: shvdpmue.m)

function dydt=shvdpmue(t,y,mue) % mit weiterer Parameteruebergabe

dydt=[y(2);mue*(1-y(1)ˆ2)*y(2)-y(1)];


options=odeset(’OutputFcn’,’odeplot’); mue=1;

[t,y]=ode45(@shvdpmue,[0 50],[2;0],options,mue); % oder


Eigenschaft fur µ = 1, 1000?µ = 1: nicht-steif, µ = 1000: steif

124 / 172


Beispiel 2: Van-der-Pol-Schwinger (Steifes Dgl.-system)

y(t) =

[y1(t)y2(t)

]=

[y2(t)

µ(1− y 2

1 (t))

y2(t)− y1(t)

], y(0) = y0 =

[20

]µ = 1000

Jacobi-Matrix, MATLAB R©-Funktion (Datei: shvdpjacob1000.m)

Jacobi-Matrix: J =∂f

∂y=

[∂f1∂y1

∂f1∂y2

∂f2∂y1

∂f2∂y2

]=

[0 1

−2µy1y2 − 1 µ(1− y 21 )

]function J = shvdpjacob1000(t,y,mue)

J = [0 1; -2*mue*y(1)*y(2)-1 mue*(1-y(1)ˆ2)];


options=odeset(’Jacobian’,@shvdpjacob1000);mue=1000;


125 / 172


Beispiel 3: DAE-Systeme

a) Implizites ODE- oder DAE-System (hochstens Index 1) liegt bereits vor

2 Schritte:

1 Konsistente Anfangsbedingungen (AB) zur Integration mit ode15i[y0mod,yp0mod] =

decic(odefun,t0,y0,fixed y0,yp0,fixed yp0,options)

y0 = [1; 0; 1e-3]; % Nicht konsistente AB

yp0 = [0; 0; 0]; % Nicht konsistente AB

tspan = [0 4*logspace(-6,6)];

M = [1 0 0;0 1 0;0 0 0]; % Masse-Matrixoptions = odeset(’RelTol’,1e-4, ...

’AbsTol’,[1e-6 1e-10 1e-6],’Jacobian’,[],M);[y0,yp0] = decic(@f,0,y0,[1 1 0],yp0,[],options)

126 / 172


Beispiel 3: DAE-Systeme (Fortsetzung)

zu a) Implizites ODE- oder DAE-System (hochstens Index 1) liegt bereits vor

2 Losen mit ode15i

[t,y] = ode15i(@f,tspan,y0,yp0,options)

MATLAB R©-Funktion (Datei: f.m)

function res = f(t,y,yp)

res = [ yp(1) + 0.04*y(1) - 1e4*y(2)*y(3);

yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)*y(2);

y(1) + y(2) + y(3) - 1];

Bsp.: Chemische Reaktion (s. o., ihb1dae.m)

Weiteres Bsp.: Diskretisierte partielle Dgl. mit beweglichem Gitter(Burgers-Gleichung; Mechanik, Nichtlineare Akustik, Fluid-, Gasdynamik,Verkehrsfluss), iburgersode.m)

127 / 172


Beispiel 3: DAE-Systeme (Fortsetzung)

b) DAE-System (Differentiationsindex > 1) liegt vor

6 Schritte (symbolisch (Symbolic Math ToolboxTM

), numerisch):

1 Definition der symbolischen Variablen und des DAE-Systems2 Falls Dgl.-system hoherer Ordnung: → Dgl.-system 1. Ordnung3 Falls Differentiationsindex > 1: Reduzierung4 Symbolisches System → numerisches

”MATLAB R© function

handle“5 Konsistente Anfangsbedingungen (AB)6 Losen mit ode15i, ode15s, ode23t

Bsp.: Planares Pendel (siehe MATLAB-Online-Dokumentation))

128 / 172


- Mehrpunkt-Randwertproblem -

Beispiel 4: Dgl. mit Anfangs- und Endbedingung (s. [1]):Mehrpunkt-Randwertproblem (Kleine Durchbiegungen am

Biegebalken; 1 +(

dydx

)2 ≈ 1)

( ) 0M x >

0 1

y

x

Dgl.: − d2y

dx2=

M(x)

EI (x)= b(x) =

5− 4x

100e−x

mit l = 1, y(0) = 0, y(1) = 0, Biegemomentverlauf M(x), Biegesteifigkeit EI (x)

129 / 172


Beispiel 4: (Fortsetzung)

Umwandlung in Dgl.-system 1. Ordnung: y1 = y , y2 = dydx

MATLAB R©-Routinen (Dgl., Randbedingungen, Losung):

function dydx = ystrich b(x,y)

dydx = [ y(2); (4*x-5)/(100*exp(-x))];

—

function res = bc b(ya,yb)

res=[ya(1);yb(1)]; % d. h. y(0) - 0 und y(1) - 0

—

xmesh=linspace(0,1,100); % Schaetzung fuer Startloesung

yinit=[0;0]; % Anfangswerte

loesinit=bvpinit(xmesh,yinit); % Anfangsloesung

loesung=bvp4c(@ystrich b,@bc b,loesinit); % Loesung

plot(loesung.x,loesung.y(1,:)) % Druck der Biegelinie

130 / 172


- Dgln. mit Totzeit(en) -

Beispiel 5: Dgl.-system mit Totzeiten

x1 = −0.4 x1(t − 1)

x2 = x1(t − 1)− 2 x2(t − 0.5)

MATLAB R©-Routinen (Dgl., Historie, Losung):

function xp = dgl Ttbsp(t,x,z)

xT1=z(:,1);xT2=z(:,2);xp=[-0.4*xT1(1); xT1(1)-2*xT2(2)];

—

function s = history Ttbsp(t)

s=[1;-2];

—

tspan=[0 6]; Tt=[1 0.5]; options=ddeset(’InitialY’,[0;0]);

loesung=dde23(@dgl Ttbsp,Tt,@history Ttbsp,tspan,options);

t interp=linspace(0,6,80); x=deval(loesung,t interp);

plot(t interp,x(1,:),t interp,x(2,:)) % Druck131 / 172


Im Gegensatz: numerische Integration zur Ermittlung desbestimmten Integrals

Flacheninhalt (Bsp.)

I =

5∫0

|x |dx

Befehl Erlauterung

quad, quad8 Numerische Integration (bestimmte Integration)

MATLAB R©-Befehl:

F=quad(’abs’,0,5);

Ergebnis:

F=12.5

132 / 172

5.2 MATLAB R©/Simulink R© - Einfuhrung5.2.6 Partielle Differentialgleichungen

Partial Differential Equation ToolboxTM

Funktionen zur Losung von PDEs (2D-, 3D-, Zeit) mittelsFinite-Differenzen-Methode

Statische, Zeit-, Frequenzbereichs-, Eigenwertprobleme

Nachbearbeitung, Visualisierung

Diffusion, Waremubertragung, Strukturmechanik, Elektrostatik,Magnetostatik, AC-Leistungsmagnetismus; nutzerspezifische, gekoppelteSysteme

Siehe Online-Nutzerdokumentation

133 / 172

5.2 MATLAB R©/Simulink R© - Einfuhrung5.2.7 Polynome, Interpolation

Darstellung von Polynomen: Koeffizienten des Polynoms in absteigenderReihenfolge der Potenzen in vektorieller Form,z. B. steht p=[1 0 -2 -5] fur p(x) = x3 − 2 x − 5

Polynomfunktionen (doc polyfun)

Befehl Erlauterung

roots(p) Wurzeln (Nullstellen) des Polynoms

poly Berechnung der Polynomkoeffizienten aus den Wurzeln(siehe auch poly bei Matrizen)

polyval(p,5) Berechnung des Wertes des Polynoms fur x = 5

conv, deconv Polynommultiplikation, -divisionBsp.: a=[1 2 3], b=[4 5 6], c=conv(a,b)

(Ergebnis: c=[4 13 28 27 18])[q,r]=deconv(c,a)

(Ergebnisse: q=[4 5 6], r=[0 0 0 0 0])

134 / 172

5.2 MATLAB R©/Simulink R© - Einfuhrung5.2.7 Polynome, Interpolation

Polynomfunktionen (Fortsetzung)

Befehl Erlauterung

polyder Ableitung eines Polynoms

polyfit Anpassung der Polynomkoeffizienten an Daten

residue Residuenberechnung (Partialbruchzerlegung)Bsp.: 2 s+5

s4+3 s3+2 s2 = As+2

+ Bs+1

+ Cs

+ Ds2

Ergebnis: A=-0.25, B=3, C=-2.75, D=2.5

Interpolation

Befehl Erlauterung

interp1 1D-Interpolation (linear, spline, kubisch)

interp2 2D–Interpolation

. . .

spline Spline-Interpolation (siehe auch Splines-Toolbox)

135 / 172

5.2 MATLAB R©/Simulink R© - Einfuhrung5.2.8 Schwach besetzte Matrizen

Funktionen fur schwach besetzte Matrizen (doc sparfun)

Befehl Erlauterung

speye Einheitsmatrix

sprand, sprandn Zufallige belegte Matrizen(Gleich-, Normalverteilung)

spdiags Diagonalmatrizen/Diagonalen (jeweils schwach besetzt)

. . .

sparse, full Umwandlung in schwach besetzte Matrix bzw.umgekehrt

spconvert Einlesen von Daten und Speichern in schwachbesetzter Matrix

. . .

136 / 172


Beispiel

D =

8 1 6.003 5 7.004 9 2.009 9 0

(Variante 1)

bzw. eine ASCII-Datei mydata.dat enthalt die o.a. Werte der Matrix D (Variante 2)

MATLAB R©-BefehleS=spconvert(D) (Variante 1)load mydata.dat (Variante 2)

137 / 172


Weitere Funktionen fur schwach besetzte Matrizen

Befehl Erlauterung

nnz Anzahl von Null verschiedener Elemente

find Auffinden von Null verschiedener Elemente

spfun Anwendung von Funktionen auf von Nullverschiedene Elemente

spy Visualisierung der Besetzheitsstruktur

BeispielZ=sprandn(100,100,0.3);

[L,U,P]=lu(Z,1e-3);

[L,U,P]=luinc(Z,’0’);

spy(Z), spy(L), spy(U), spy(P)

Algorithmen zur Umordnung

Funktionen zur lineraren Algebra

Iterative Methoden zur Losung linearer Gleichungen, usw.

138 / 172

5.2 MATLAB R©/Simulink R© - Einfuhrung5.2.9 M-File-Programmierung

– Erganzung zu Abschnitt 5.2.3; weitere Sprachelemente (doc lang) –

Skripte, Funktionen (doc function)

Skripte, Funktionen und Unterfunktionen mit variabler Anzahl von Argumentenfunction c=test(a,b)

if nargin== 1,

c = a. 2; % Quadratzahl einer Zahl

elseif nargin==2

c = a + b; % Summe zweier Zahlen

end

139 / 172


Lokale und globale Variable

Variable innerhalb von Funktionen lokal

Wenn Variable vor ersten Auftreten mit global vereinbart, Zugriff in allenFunktionen bzw. im Workspace, in denen Variable auch mit globalvereinbart wurde, auf die gleiche Variable

Beispiel: Rauber-Beute-Modell

y1(t) = y1 − α y1(t) y2(t)

y2(t) = −y2(t) + β y1(t) y2(t)

function [yp]=sh lotka(t,y)

% LOTKA Lotka-Volterra predator-prey model

global ALPHA BETA

yp = [y(1) - ALPHA*y(1)*y(2); -y(2) + BETA*y(1)*y(2)];

Befehle im Workspace:global ALPHA BETA

ALPHA=0.01; BETA=0.02;

[t,y]=ode23(@sh lotka,[0 50],[1;1]); plot(t,y)

140 / 172


Programmablaufsteuerung

if, else, elseif - Verzweigung

Beispiel: Test, ob nicht negative Zahl gerade oder ungeradeif n < 0, % n negativ ? : Fehlermitteilung

disp(’n muss nicht-negativ sein’)

elseif rem(n,2) == 0,

disp(’n gerade’)

else

disp(’n ungerade’)

end

141 / 172


Programmablaufsteuerung (Fortsetzung)

switch - Ausfuhrung bestimmter Anweisungen in Abhangigkeit des Werteseiner Variablen oder eines Ausdruckes

Beispiel:switch zahl

case -1

disp(’minus Eins’)

case 0

disp(’Null’)

case 1

disp(’plus Eins’)

case 2,3,4

disp(’plus Zwei, Drei oder Vier’)

otherwise

disp(’anderer Wert’)

end

142 / 172


Programmablaufsteuerung (Fortsetzung)

while - Ausfuhrung von Anweisungen solange ein Ausdruck wahr ist

Beispiel:n=1;

while prod(1:n) < 1e100,

n = n + 1;

end

Ausstieg aus while-Schleife mit break moglich

for - Laufanweisung

Beispiel:A=zeros(6);

for i=1:6; for j=1:2:6;

A(i,j)=j*3+i;

end; end

143 / 172



Befehl Erlauterung

break Unterbrechung von Schleifen

pause Pause (Benutzereingabe, Zeitverzogerung)

keyboard Eingaben von Tastatur (aus M-File heraus aufgerufen)

return Rucksprung aus Funktionen, Beenden des keyboard-Modus

eval, feval Zeichenkettenauswertung und -berechnungBeispiele:n=3; t=’1/(i+j-1)’;

for i=1:n; for j=1:n;

a(i,j)=eval(t);

end; end

for i=1:10,eval([’P’,int2str(i),’=i. 2’]);end

fun=[’sin’;’cos’;’log’];

k=input(’Waehle Funktions-Nummer: ’);

x=input(’Gib einen Zahlenwert ein: ’);

feval(fun(k,:),x)

144 / 172



Befehl Erlauterung

error, warning Fehlermitteilung (danach Programmstopp),Warnung (danach Programmfortsetzung)

Bsp.: if n <= 1, error(’n muss groesser als 1 sein!’),end

. . .

Debugger Programmtest (Schrittbetrieb, Unterbrechungspunkte, . . .)

Hinweise zur Optimierung des MATLAB R©-Codes

Vektorisierung der Schleifen

Anstelle von

i=0; for t=0:0.01:10; i=i+1; y(i)=sin(t); end

ist die vektorisierte Version besser:

t=0:0.01:10;y=sin(t);

Vorherige Speicherplatzzuweisung fur Matrizen:

y=zeros(1,100); for i=1:100, y(i)=det(x i);end

145 / 172

5.2 MATLAB R©/Simulink R© - Einfuhrung5.2.10 Visualisierung

2D-Darstellung (doc graph2d)

Befehl Erlauterung

plot 2D-Grafik mit linearer AchsenteilungBeispiele: plot(x), plot(A), plot(x,y)

plot(x,y,’m--’), plot(x,y,’o’,x,2*y,’g:’)

plotyy 2D-Grafik mit linearer Achsenteilung,2 Ordinaten mit unterschiedlicher Skalierung

loglog 2D-Grafik mit logarithmischer Achsenteilung

semilogx, semilogy 2D-Grafik mit linearer Teilung der Abszisse undlogarithmischer Teilung der Ordinate u. umgekehrt

146 / 172


Grafik- und Achsenbeschriftung, allgemeine Grafikbefehle (docgraph2d, doc graphics)

Befehl Erlauterung

title Grafiktitel

xlabel, ylabel Abszissen-, Ordinatenbeschriftung

legend Grafiklegende

text Text an Position (xpos,ypos); Bsp. text(2,3,’Min.’)

gtext Platzierung des Textes mittels Maus

axis Achsenskalierung; Bsp. axis([xmin xmax ymin ymax])

grid Gitternetzlinien

hold on Beibehalten der Grafik (weitere Kurven in gleiche Grafik)

hold off Freigabe der Grafik zum Uberschreiben

subplot Unterteilung eines Bildes (figure) in mehrere GrafikenBsp.: subplot(3,2,1) 6 Grafiken, 3x2, Ansprechen der 1. Grafik

147 / 172



Befehl Erlauterung

figure Erzeugung von Bildern (Grafikfenstern)/Bild im Vordergrund

clf Loschen des Bildinhaltes

close [all] Schließen des/aller Fenster(s)

gcf”Handle“ des aktuellen Bildes (get current figure)

gca”Handle“ des aktuellen Achsenkreuzes

(get current axis)

get, set Lesen bzw. Setzen von Eigenschaften

print Drucken von Grafikfensterinhalten

Grafikobjekte: figure, axes, line, text, patch, surface, image, light, uicontrol,uimenu

148 / 172



Befehl Erlauterung

plot3 Linien und Punkte im 3D-Raum

meshgrid, mesh 3D-Werteberechnung und 3D-Netz

surf 3D-Oberflache

colorbar Farbskalierung

colormap Farbeinstellung

view Beleuchtung, Schattierung, Ansicht (Blickwinkel)

rotate, rotate3d interaktive Ansicht

149 / 172


Spezielle grafische Darstellungen

Befehl Erlauterung

bar, barh, Balkendiagrammebar3, bar3h

area ausgefullte Flachen unter Kurven

comet”Kometen-Bahn“

errorbar Funktionsverlauf mit”Fehlerbereich“

ezplot, fplot Grafik von Funktionen; Bsp.: ezplot(’sin(x)’)

hist Histogramm (Haufigkeitsverteilung)

pie, pie3 2D- bzw. 3D-Tortendiagramm

stem, stem3, 2D-, 3D-Verlauf mit Stutzstellen,stairs 2D-Stufenverlauf

150 / 172


Spezielle grafische Darstellungen (Fortsetzung)

contour, contourf, Iso-Liniendarstellung (2D, ausgefullt)contour3 Iso-Liniendarstellung (3D)

pcolor”Farbausdruck“ einer Matrix

meshc, meshz, kombinierte Darstellungen (raumlich, Iso-Linien)waterfall, surfc

slice Scheiben in einer bestimmten Ebene

Bildverarbeitungsbefehle, Image Processing Toolbox

151 / 172

5.2 MATLAB R©/Simulink R© - Einfuhrung5.2.11 Simulink R© – Simulation dynamischer Systeme

Modellierung, Simulation und Analyse dynamischer Systeme

Erzeugen von Simulationsmodellen mit Hilfe hinterlegter Elemente inBibliotheken: Continuous, Discrete, Nonlinear, Sources, Sinks, MathOperations, Signal Routing, Ports & Subsystems, User-Defined Functions,Simulink Extras

Verbinden und Parametrieren der Elemente (Blocke); auch selbstprogrammierte Inhalte

Einstellen der Experimentierbedingungen

Simulation

Auswertung

Merkmale

Zusammenfassung von Blocken zu Untermodellen und Bibliotheken

Kopplung zwischen MATLAB R© und Simulink R©, MATLAB R©/Simulink R© undanderen Programmiersprachen (z. B. C, C++, FORTRAN)

Simulink R©-Debugger

152 / 172

5.2 MATLAB R©/Simulink R©

5.2.11 Simulink – Simulation dynamischer Systeme

Beispiel: Regelkreis

SRw u y

Z. B. Ubertragungsfunktionsblocke fur Strecke GS (s) = 1(T1s+1)(T2s+1)

,

T1 = 5, T2 = 1 und Regler GR (s) = KR = 5 (Bibl. Continuous)

Summationsstelle (Math Operations), Quelle/Signal = Fuhrungsgroße(Sources), Senke/Signal = Regelgroße (Sinks), Verbindungslinien

153 / 172


5.2.11 Simulink – Simulation dynamischer Systeme

Beispiel: Tanksystem

Ausfuhrliche Behandlung in der Ubung

Systembeschreibung

Aufbau eines Simulink-Simulationsdiagramms: eintank.mdl, eintank1.mdl

- Modell mit interner/externer Eingangsgroße (bzw. *.slx)

154 / 172


5.2.11 Simulink R© – Simulation dynamischer Systeme

Beispiel: Tanksystem (Fortsetzung)

Simulation

Arbeitspunktbestimmung (et l.mdl, externe Ein- und Ausgangsgroße) undLinearisierung

Sprungantwort des linearisierten Modells (et l ss.mdl,Zustandsraumdarstellung)

Verhalten im Regelkreis mit PI-Regler (et vlim.mdl)

Verhalten im Regelkreis mit Zwei-Punkt-Regler (et zpr.mdl)

155 / 172



Kopplung MATLAB R©-Simulink R©, Simulink R©-MATLAB R©

MATLAB R©-Workspace (From/To Workspace-), MAT-Files (From/ToFile-Blocke)

MATLAB R© unter Simulink R© nutzen (Interpreted MATLAB Fcn, MATLABFunction, Fcn, S-function (system function))

Simulink unter MATLAB nutzen (MDL-File), Bsp:

[t,x,y]=sim(’eintank1’,[0 500],[],[[0:500]’ ...

8e-5*ones(501,1)]);

156 / 172



Warum S-functions? (Level-2 S-function als M-File)

Hinzufugen neuer algorithmischer Blocke

Einbeziehung z.B. existierenden C-Codes

Systembeschreibung als mathematisches Gleichungssystem

Grafische Animationen

Arbeitsweise von S-functions

Simulink-Block mit Eingang, Ausgang und Zustand x = [xTc xT

d ]T mit xc -kontinuerlicher Zustandsvektor und xd - diskontinuierlicher Zustandsvektor;Achtung! : Direkter Durchgriff → Gefahr algebraischer Schleifen

Simulationsschritte in Simulink (wiederholter Aufruf von Blocken)

157 / 172



Aufbau von M-File S-functions

M-File S-Function Template (Vorlage):

matlabroot\toolbox\simulink\blocks\sfuntmpl.m (doc sfuntmpl, typesfuntmpl)

Empfehlung: Orientierung an Beispielen

Simulationsschritt Funktion

Initialisierung mdlInitializeSizes

Berechnung nachster Abtastzeitpunkt mdlGetTimeOfNextVarHit(optional bei variabler Abtastzeit)

Berechnung der Ausgange mdlOutputs

Aktualisierung der diskreten Zustande mdlUpdate

Berechnung der Ableitungen mdlDerivatives(rechte Seite der Dgl.)

Beenden der Simulation mdlTerminate

158 / 172



Anpassung von M-File S-functions am Bsp. Tanksystem;Simulink-Modell: et m.mdl, M-FIle S-function: et mfsf.m

Initialisierung: (mdlInitializeSizes)Sizes.NumContStates=1; % Anzahl kontinuierl. Zustande

Sizes.Outputs=1; % Anzahl Ausgange

Sizes.Inputs=1; % Anzahl Eingange

Sizes.DirFeedthrough=0; % Durchgriff

x0=0.2; % Anfangszustand

Ausgangsgleichung: (mdlOutputs)sys=x; % Ausgang(1)=Zustand(1)

Ableitung (rechte Seite der Dgl.):sys=64.9612*(u(1)-3e-5*sqrt(2*9.81*x(1)));

159 / 172



Aufbau von C-MEX S-functions

C-MEX S-Function Template (Vorlage):

matlabroot\simulink\src\sfuntmpl\basic.c(type matlabroot\simulink\src\sfuntmpl\basic.c)

Empfehlung: Orientierung an Beispielen

Simulationsschritt Funktion

Initialisierung/Dimensionen mdlInitializeSizes

Initialisierung/Anfangszustand mdlInitializeConditions

Berechnung der Ausgange mdlOutputs

Aktualisierung der diskreten Zustande mdlUpdate

Berechnung der Ableitungen mdlDerivatives(rechte Seite der Dgl.)

Beenden der Simulation mdlTerminate

160 / 172



Anpassung von CMEX S-functions am Bsp. Tanksystem;Simulink-Modell: et c.mdl, M-FIle S-function: et cmsf.c

Typ und Name der S-function:# define S FUNCTION LEVEL 2

# define S FUNCTION NAME et cmsf

Initialisierung/Dimensionen: (mdlInitializeSizes)ssSetNumContStates(S, 1); /* Anzahl kontinuierl. Zustande */

ssSetOutputPortWidth(S, 0, 1); /* Anzahl Ausgange */

ssSetInputPortWidth(S, 0, 1); /* Anzahl Eingange */

ssSetInputPortDirectFeedThrough(S, 0, 1); /* Durchgriff */

Initialisierung/Anfangszustand: (MDL INITIALIZE CONDITIONS)*x0=0.2; % Anfangszustand

Ausgangsgleichung: (mdlOutputs)y[0]=x[0];

Ableitung (rechte Seite der Dgl.):dx[0]=64.9612*(u(0)-3e-5*sqrt(2*9.81*x[0]));

161 / 172



Einbeziehung der CMEX S-function et cmsf.c

Ubersetzung: mex et cmsf.c; Erzeugung einer dynamischen Bibliothek(WINDOWS: et cmsf.mexw32, et cmsf.mexw64, Linux: et cmsf.mexglx,et cmsf.mexa64)

Simulationsoptionen (Menu oder simset) und Simulation

Hinweise zur effizienten Simulation

Vermeidung von”Interpreted MATLAB R©-Fcn-Blocken“ (besser: Fcn mit

elementaren mathematischen Funktionen)

Vermeidung von M-File S-functions (etwas besser: CMEX S-function)

Nachteilig: Memory-Blocke

Maximale Schrittweite zu klein, Zeitbasis zu lang?

Zu große Genauigkeit? (Standard: 0.1 % meist ausreichend)

Liegt steifes Dgl.-system vor? (ode15s)

Verschiedene Abtastzeiten?

Algebraische Schleifen? 162 / 172


5.2.12 SimscapeTM

und Objektorientierung

Objektorientierte Erweiterungen

SimscapeTM

: Modellierungs- und Simulationsumgebung fur physikalischeSysteme (domaneubergreifend); Bibliotheken; Erweiterbarkeit;Modellaustausch; HiL; C-Code-Export

SimPowerSystemsTM

(fruher: SimPowerTM

): ElektrischeEnergieversorgungssysteme

SimMecanicsTM

: Mehrkorper-Systeme (3D-Mechanik)

SimElectronicsTM

: elektronische und mechatronische Systeme

163 / 172

5.3 Modelica R© und OpenModelica

Modelica R©

Freie, objektorientierte, gleichungsbasierte Sprache zur zweckmaßigen Modellierungphysikalischer Systeme (https://modelica.org)

OpenModelica

Quelloffene Modelica R©-basierte Modellierungs- und Simulationsumgebung(https://openmodelica.org) mit OMEdit (OpenModelica Connection Editor:quelloffener benutzerfreundliche grafische Benutzerschnittstelle, Simulation,Ergebnisdarstellung, Modellaustausch (engl.: FMI = Functional Mockup Interface))

Eine kommerzielle Alternative

DYMOLA R©: DYnamic MOdeling LAboratory http://www.3ds.com/de/produkte-und-services/

catia/funktionsumfang/catia-systementwicklung/modelica-systemsimulation/dymola/

164 / 172

https://modelica.org

https://openmodelica.org

http://www.3ds.com/de/produkte-und-services/catia/funktionsumfang/catia-systementwicklung/modelica-systemsimulation/dymola/

http://www.3ds.com/de/produkte-und-services/catia/funktionsumfang/catia-systementwicklung/modelica-systemsimulation/dymola/


Weitere Werkzeuge

OMC: OpenModelica-Compiler

OMShell: OpenModelica-Shell

OMNotebook: Mathematica-ahnliches Notebook fur Lehrzwecke

DrControl under OMNotebook: Interaktives Kursmaterial furRegelungstechnik

OMOptim: Optimierung

OMPython: Schnittstelle zur Python-Umgebung

. . .

165 / 172


Struktur der OpenModelica-Umgebung

Bildquelle: https://openmodelica.org

166 / 172

https://openmodelica.org


Modellierungs- und Simulationsumgebung (allg.)

Bibliothek Modell

Übersetzung

ÜbersetztesModell

p =

x =0

t , t0 1

AndereProgramme

Simulations-programm

AnimationGrafik- und

Numerikausgabe

Modelleditor

Simulationssteuerung

(Parameter,Anfangswerte,Zeithorizont,Simulation,Skripte)

Visualisierung

167 / 172


Vorgehen bei der Modellerstellung

Klassenkonzept (ahnlich objektorienierter Programmierung)

Hierarchische Modellerstellung (Objektorientierung und Gleichungen)

Festlegung von Potenzial- und Flussvariablen

Definition von Konnektoren mit Potenzial- und Flussvariablen

Verallgemeinerbare Sachverhalte extrahieren (Teilmodelle (Klassen)),z. B. elektrischer Zweipol)

Spezifizierung mit Vererbung (z. B. Ohmscher Widerstand)

Verbinden von Objekten

→ Gesamtmodell (”Bottom-Up“-, (

”Top-Down“)-Entwurf)

168 / 172


OpenModelica Connection Editor

Verschiedene Ansichten/Fenster (Modell (Designer-, Hauptfenster), Suche,Bibliotheken-, Modell-Navigator (Hierarchische Struktur; Bibliotheken:Elektrik, Magnetismus, Mechanik, Hydraulik, Flussigkeiten, Medien (Stoffe),Thermodynamik, Mathematik, Blocke, Zustandsgraphen, SI-Einheiten),Dokumentation, Variablen, Meldungen)

Designer-Fenster: Hauptfenster (Symbol (Icon), Diagramm, Modelica-Text,Dokumentation)

169 / 172


Bearbeiten von Aufgabenstellungen (OMEdit)

1.) Simulationsdiagramm aus Bibliothekselementen (Bsp.: RC-Netz mitSpannungsquelle und Erdung, s. Bsp., Abschnitt 4.3)

2.) Problemstellung bearbeiten anhand math.-physikal. Beschreibung;Programmieren in Modelica (Bsp.: planares, reibungsfreies Pendel, s. o.)

Math.-physikal. Beschreibung

J ω = −m g l sin(ϕ)

ϕ = ω

mit m - Masse, l - Lange, g - Erdbeschleu-nigung, J - Tragheitsmoment (J = m l2), ϕ -Auslenkwinkel, ω - Winkelgeschwindigkeit,ω - Winkelbeschleunigung.

Illustration

170 / 172


Bearbeiten von Aufgabenstellungen (Fortsetzung)

Programmtechnische Umsetzung; aufeinander folgende Wahl der

Reiter/Menupunkte/Fenster:

Im Hauptfenster: Deklaration von Variablen und Parametern;Werte; Modellgleichungen

parameter Real m=1;

parameter Real l=1;

parameter real g=9.81;

parameter Real J=m*l 2;

Real phi(start=1),w;

der(phi)=w;

J*der(w)=-m*g*l*sin(phi);

Datei/SpeichernSimulation/Simulieren - Ubersetzen und Simulation (mitSimulationsparametereinstellungen: Start-, Endzeit,Integrationsverfahren (Standard: dassl, (Euler-, Runge-Kutta-,Radau-, Lobatto-, Inline-Verfahren))Im Variablen-Fenster: Variablen auswahlen furAnzeige/Speichern/Druck (→ Reiter: Plotten)

171 / 172


Bearbeiten von Aufgabenstellungen (Fortsetzung),Integrationsverfahren, Skripte

3.) Teilweise Nutzung von Bibliothekselementen, teilweise Erstellung neuer

Modelle, Aufbau eines Gesamtmodells (s. o., ausfuhrlich: Ubung)

Nutzung des vordefinierten Konnektors”Pin“

Definition eines Teilmodells”Elektrisches Bauelement mit zwei

Anschlussen“Nutzung des Teilmodells und Spezifizierung eigener Modelle furKonstantspannungquelle, Widerstand und InduktivitatAufbau des Stromkreises

Numerische Integrationsverfahren

Skripte (nutzbar unter OMShell)172 / 172

Simulation – Fach im Modul „Modellbildung und Simulation ... · Fakult at f ur Informatik und...

Documents

Transcript of Simulation – Fach im Modul „Modellbildung und Simulation ... · Fakult at f ur Informatik und...