Behandlung von Unstetigkeiten

1. Wir wollen uns heute dem Problem der Behandlung von Unstetigkeiten in der Modellbeschreibung zuwenden. Modelle aus den Ingenieurbereichen weisen häufig Unste-tigkeiten auf, die Schaltvorgänge, Begrenzer, trockene Reibung, Stösse oder ähnliche Vorgänge repräsentieren mögen. Die Modellierungsumgebung muss sich dieser Vorgänge speziell annehmen, da sie in starkem Masse die Numerik der Differentialgleichungslöser beeinflussen. Behandlung von Unstetigkeiten

2. Die numerische Lösung von Differentialgleichungen Unstetigkeiten in Zustandsgleichungen Integration über Unstetigkeiten Zustandsereignisse Ereignisbearbeitung Mehrwertige Kennlinien Der elektrische Schalter Die ideale Diode Die Reibungskennlinie Übersicht

3. Alle heute im Einsatz stehenden numerischen Differentialgleichungslöser arbeiten mit Polynom-extrapolation. Der Wert einer Zustandsvariablen x zum Zeitpunkt t+h, wobei h die momentane Schrittlänge ist, wird ermittelt, indem durch bereits bekannte Stützwerte von x und dx/dt der jetzigen Zeit t sowie der Vergangenheit ein Polynom nter Ordnung hindurchgelegt wird. Der Wert dieses Polynoms zum Zeitpunkt t+h repräsentiert die angenäherte Lösung der Differentialgleichung. Bei impliziten Verfahren wird die Zustandsableitung zum Zeitpunkt t+h ebenfalls als Stützwert beigezogen. Numerische Differentialgleichungslöser

4. · x(t+h)  x(t) + h · x(t) · x(t+h)  x(t) + h · x(t+h) Beispiele Explizites Eulerverfahren erster Ordnung: Implizites Eulerverfahren erster Ordnung:

5. Polynome sind immer stetig und stetig differenzierbar. Wenn also die Zustandsraumbeschreibung des Systems: eine Unstetigkeit aufweist, ist die Polynomextrapolation eine sehr schlechte Approximation der Wirklichkeit. Somit weisen Verfahren mit fester Schrittlänge einen grossen Integrationsfehler auf, während Verfahren mit variabler Schrittlänge den Schritt in der Umgebung der Unstetigkeit sehr stark reduzieren. · x(t) = f(x(t),t) Unstetigkeiten in der Zustandsbeschreibung

6. Ein Verfahren variabler Schrittlänge reduziert die Schrittlänge bei jeder Unstetigkeit. Nach der Unstetigkeit wird die Schrittlänge nur langsam wieder vergrössert, da das Verfahren nicht unterscheiden kann zwischen einer Unstetigkeit und einem Punkt grosser lokaler Steifigkeit (mit grossem absolutem Wert der Ableitung). Unstetigkeiten Die Schrittlänge ist dauernd viel zu klein. Das Verfahren ist zu-mindest sehr ineffizient, falls nicht ungenau. Integration über Unstetigkeiten h t

7. Diese Probleme können vermieden werden, wenn dem System explizit mitgeteilt wird, wo es Unstetigkeiten in der Modellbeschreibung gibt. 3 1 2 f(x) fp f = fm xm a f = m·x x xp f = fp fm m = tg(a) f = ifx < xm then fm else if x < xp then m*x elsefp ; Das Zustandsereignis Beispiel: Begrenzerfunktion

8. Iteration x xp Schwellwert t h h f(x) fp Modellumschaltung h x xp xm a x xp t t h fm Schrittlängenreduktion während Iteration Ereignis Die Ereignisbearbeitung I

9. h t h Schrittlänge ohne Ereignis-behandlung Schrittlänge mit Ereignis-behandlung t Die Ereignisbearbeitung II

10. In Modelica wird die Unstetigkeit durch ein if-Statement repräsentiert. Bei der Übersetzung werden diese Statements in korrekte Ereignisbeschreibungen (Sätze von Modellen mit Um-schaltbedingungen) abgebildet. Der Benützer braucht sich um die Ereignisbeschreibungs-mechanismen selbst nicht zu kümmern. Diese bleiben hinter den if-Statements verborgen. f = ifx < xm then fm else if x < xp then m*x elsefp ; Repräsentation von Unstetigkeiten

11. Der Benützer sollte beachten, dass während der Iteration die unstetige Kennlinie temporär verlassen wird. ist gefährlich, da absDp temporär negativ werden kann. löst dieses Problem. Dp = p1 – p2 ; absDp = ifDp > 0 thenDp else –Dp ; q = sqrt(absDp) ; q =  | Dp |  Dp = p1 – p2 ; absDp =noEvent(ifDp > 0 thenDp else –Dp ) ; q = sqrt(absDp) ; Probleme

12. Die noEvent-Klausel sorgt dafür, dass if-Statements oder Boole’sche Ausdrücke, die normalerweise in Simulations-code übersetzt würden, bei welchem die Diskontinuitäten durch Ereignisse abgefangen werden, so stehengelassen werden, wie sie im Modellierungscode formuliert wurden. Dadurch wird die korrekte Simulation dieser Diskontinuitäten der Schrittsteuerung des numerischen Integrationsverfahrens überlassen. Dp = p1 – p2 ; absDp =noEvent(ifDp > 0 thenDp else –Dp ) ; q = sqrt(absDp) ; Die „noEvent“-Klausel

13. Die bisher eingeführten Sprachelemente reichen nicht aus, um damit mehrwertige Kennlinien, wie zum Beispiel Hysteresefunktionen zu beschreiben. Wenn xgrösser wird als xp, muss f von fm nach fp umgeschaltet werden. Wenn xkleiner wird als xm, muss f von fp nach fm umgeschaltet werden. f(x) fp x xm xp fm Mehrwertige Kennlinien I

14. f(x) fp when initial() then reinit(f, fp); end when; when x > xp or x < xm then f = if x > 0 then fp else fm; end when; x xm xp Ausgeführtam Anfangder Simulation. wird grösser } fm Diese Statements werden nur ausgeführt, wenn entwederxgrösser wird als xpoder aber wennxkleiner wird als xm. ist grösser Mehrwertige Kennlinien II

15. i u 0 = ifoffenthenielse u ; Der elektrische Schalter I Wenn der Schalter offen ist, ist der Strom i=0. Wenn der Schalter zu ist, ist die Spannung u=0. Das if-Statement in Modelica ist akausal. Es wird sortiert, genau wie alle anderen Statements.

16.  0 = s · i + ( 1 – s ) · u Schalter offen: SF f = 0 e  s Sw Die Kausalität des Switch Elements ist eine Funktion des Werts des Steuersignals s. Schalter geschlossen: f e = 0 SE Der elektrische Schalter II Mögliche Implementation: Schalter offen: s = 1 Schalter geschlossen: s = 0

17. Wenn u <0ist, ist der Schalter offen. Kein Strom fliesst. i Schalter geschlossen i u u Wenn u >0ist, ist der Schalter geschlossen. Strom kann fliessen. Die ideale Diode verhält sich wie ein Kurzschluss. Schalter offen f D offen = u <0; 0 = ifoffenthenielse u ; e Die ideale Diode I

18. Da Strom, wenn er durch die Diode fliesst, nicht einfach unterbrochen werden kann, muss das Diodenmodell noch leicht modifiziert werden. Die Variable offen muss als Boolean deklariert werden. Ihr wird der Wert des logischen Ausdrucks rechts vom Gleichheitszeichen zugewiesen. offen = u <= 0and noti > 0; 0 = ifoffenthenielse u ; Die ideale Diode II

19. Kompliziertere Vorgänge, wie zum Beispiel die Reibungs-kennlinie müssen sorgfältig analysiert werden. Der Vorgang wird hier an Hand der Reibung erklärt. fB R0 Rm Gleitreibung v Haftreibung -Rm -R0 Die Reibungskennlinie I Wenn v  0 , ist die Rei-bungskraft eine Funktion der Geschwindigkeit. Wenn v = 0 , wird die Reibungskraft so gerech-net, dass die Geschwin-digkeit 0 bleibt.

20. Wir unterscheiden fünf Fälle: Haftung: Die Reibungskraft kompensiert die Summe der angrei-fenden Kräfte, es sei denn |Sf | > R0 . v = 0 a = 0 Vorwärtsbewegung: v > 0 Die Reibungskraft berechnet sich zu: fB = Rv·v + Rm. Rückwärtsbewegung: v < 0 Die Reibungskraft berechnet sich zu: fB = Rv·v -Rm. Einsetzen der Vorwärtsbewegung: v = 0 a > 0 Die Reibungskraft berechnet sich zu: fB = Rm. Einsetzen der Rückwärtsbewegung: v = 0 a < 0 Die Reibungskraft berechnet sich zu: fB =-Rm. Die Reibungskennlinie II

21. Der genaue Vorgang kann in einem Zustandsübergangs-diagramm dargestellt werden. Anfang v < 0 v > 0 v = 0 Sf < -R0 Sf > +R0 Rückwärts-beschleunigung (a < 0) Vorwärts-beschleunigung (a > 0) Haftung (a = 0) Rückwärts-bewegung (v < 0) Vorwärts-bewegung (v > 0) a  0 and not v < 0 a  0 and not v > 0 v  0 v  0 v < 0 v > 0 Das Zustandsübergangsdiagramm

22. Das Reibungsmodell I model Reibung; parameter Real R0, Rm, Rv; parameter Boolean ic=false; Real fB, fc; Boolean Haft(final start = ic); Boolean Vor(final start = ic), Rück(final start = ic); Boolean StartVor(final start = ic), StartRück(final start = ic); fB = if Vor then Rv*v + Rm else if Rück then Rv*v - Rm else if StartVor then Rm else if StartRück then -Rm else fc; 0 = if Haft or initial() then a else fc;

23. when Haft and not initial() then reinit(v,0); end when; Vor = initial() and v > 0or pre(StartVor)and v > 0or pre(Vor) and not v <= 0; Rück = initial() and v < 0or pre(StartRück)and v < 0or pre(Rück) and not v >= 0; Das Reibungsmodell II

24. StartVor = pre(Haft) and fc > R0or pre(StartVor)and not (v > 0or a <= 0and not v > 0); StartRück = pre(Haft) and fc < -R0or pre(StartRück)and not (v < 0or a >= 0and not v < 0); Haft = not (Vor or Rück or StartVor or StartRück); end Reibung; Das Reibungsmodell III

25. Cellier, F.E. (1979), Combined Continuous/Discrete System Simulation by Use of Digital Computers: Techniques and Tools, Swiss Federal Institute of Technology, ETH Zürich, Switzerland. Elmqvist, H., F.E. Cellier, and M. Otter (1993), “Object-oriented Modeling of Hybrid Systems,” Proc. ESS'93, SCS European Simulation Symposium, Delft, The Netherlands, pp.xxxi-xli. Cellier, F.E., M. Otter, and H. Elmqvist (1995), “Bond Graph Modeling of Variable Structure Systems,” Proc. ICBGM'95, 2nd SCS Intl. Conf. on Bond Graph Modeling and Simulation, Las Vegas, NV, pp.49-55. Referenzen I

26. Elmqvist, H., F.E. Cellier, and M. Otter (1994), “Object-oriented Modeling of Power-electronic Circuits Using Dymola,” Proc. CISS'94, First Joint Conference of International Simulation Societies, Zurich, Switzerland, pp. 156-161. Glaser, J.S., F.E. Cellier, and A.F. Witulski (1995), “Object-oriented Switching Power Converter Modeling Using Dymola with Event-handling,” Proc. OOS'95, SCS Object-Oriented Simulation Conference, Las Vegas, NV, pp. 141-146. Referenzen II