Das Modell läuft, aber die Simulation ist langsam? Dymola warnt vor Steifheit? Dieser Blog erklärt, was steife Gleichungssysteme sind, wie man sie erkennt und mit welchen Ansätzen man die Simulation stabiler und effizienter macht.
chatGPT
Das Simulationsmodell ist korrekt aufgebaut, die Simulation läuft, doch sie rechnet viel langsamer als erwartet? In Dymola wird die Warning „Probably the communication interval is too large or the system is stiff.“ angezeigt?
Ein Grund dafür kann ein steifes Gleichungssystem sein. Was bedeutet das? Wie finde ich die Gründe für die Steifheit? Und wie schaffe ich es, dass meine Simulation effizienter rechnet? All das erfahrt ihr in diesem Blogbeitrag!
Finden in einem dynamischen Simulationsmodell gleichzeitig Änderungen auf sehr schnellen und auf sehr langsamen Zeitskalen statt, spricht man von Steifheit oder steifen Systemen. Solche Systeme stellen Lösungsalgorithmen vor Probleme, verlangsamen die Simulation und führen sogar zu falschen Ergebnissen oder zum Abbruch der Simulation.
Ein Maß für die Änderungsgeschwindigkeiten von Zuständen in unserem Simulationsmodell sind die Zeitkonstanten. Alles rund um das Thema Zeitkonstanten haben wir in unserem Blogartikel zur Zeitkonstante behandelt. Im Folgenden knüpfen wir daran an und erläutern, wie man von Zeitkonstanten auf die Steifheit von System schließen kann und die dadurch entstandenen Probleme löst.
Beispielhaft betrachten wir ein Simulationsmodell aus zwei thermischen RC-Gliedern. Eine Temperaturquelle mit periodisch schwankender Temperatur ist über je einen Wärmewiderstand R mit je einer Wärmekapazität C verbunden. Dabei sind der Widerstand und die Kapazität des ersten RC-Gliedes im Vergleich zum zweiten RC-Glied groß. Das zweite RC-Glied hat also eine viel kleinere Kapazität mit einem kleineren Widerstand.
Ein vergleichbares System wäre zum Beispiel eine Sauna, die wir mit einem Ofen aufheizen wollen. Die kleine Wärmekapazität stellt eine Ofeneinfassung aus Metall dar, die über einen kleinen Wärmewiderstand direkt mit der Temperaturquelle – also dem Ofen – verbunden ist. Die große Kapazität ist der Beton-Boden, der sich langsam erhitzt und einen deutlich größeren Wärmewiderstand hat, da er nicht diekt mit dem Ofen verbunden ist. Im Laufe des Artikels werden wir immer wieder auf dieses Beispiel zurückgreifen.
Wie aus dem Blogbeitrag zur Zeitkonstante bekannt, kann in diesem einfachen System die Zeitkonstante als Produkt der Kapazität C mit dem Widerstand R berechnet werden (Zeitkonstante: \(\tau\) = RC).
Haben wir also einen sehr kleinen und einen sehr großen Widerstand, führt das zu einer sehr kleinen und einer sehr großen Zeitkonstante und damit zu gleichzeitig schneller und langsamer Dynamik im System.
Was passiert, wenn wir dieses System in Dymola simulieren?
Ziel der Simulation ist es herauszufinden, wie lange es dauert, bis sich ein stationärer Systemzustand eingestellt hat. Das steife System wird also so lange simuliert, bis beide Wärmekapazitäten stabile Temperaturen erreicht haben. Da die Quelltemperatur zwischen 75 °C und 80 °C schwankt, simulieren wir so lange, bis beide Kapazitäten, die bei 25 °C starten, erstmals eine Temperatur von über 77 °C erreichen.
Die Simulation mit dem CVode Solver in der Simulationssoftware Dymola zeigt, dass die Kapazitäten die Grenztemperatur von 77 °C nach rund 646 simulierten Stunden erreicht haben. Der Solver benötigt etwa 2,9s für diese Simulation und macht dafür 1.578.812 Simulationsschritte.
Dies zeigt uns einen der großen Nachteile von steifen Systemen: Sie sind sehr rechenintensiv und dadurch sehr langsam.
Bevor wir zum Debugging von steifen Systemen kommen, müssen zunächst einige Hintergründe beleuchtet werden. In diesem Abschnit schauen wir uns nur lineare Gleichungsysteme and. Die Basis hierfür bildet die folgende Zustandsraumdarstellung eines linearen Differentialgleichungssystems:
$$\frac{dx}{dt} = Ax + Bu$$
Dabei ist für uns vor allem die Systemmatrix A interessant. Sie gibt an, wie sich die Zustandsvariablen x ohne den Einfluss von externen Einflüssen u zeitlich verändern. Daher spricht man auch von der Eigendynamik.
Die Systemmatrix A besteht in linearen Gleichungssystemen ausschließlich aus konstanten Einträgen.
Die Zustandsvariablen sind die Unbekannten eines Gleichungssystems, aus denen alle anderen Outputs berechnet werden können.
Die Hauptaufgabe eines Solver bei einer Simulation ist die Zustandsvariablen x im jeweils nächsten Zeitschritt zu ermitteln.
Mithilfe der Systemmatrix A können wertvolle Informationen über das Systemverhalten abgeleitet werden. Dazu ist es notwendig, sich die Eigenwerte dieser Matrix anzusehen. Eigenwerte sind eine skalare Charakteristik von Matrizen, die in der komplexen Zahlenebene angeordnet sind.
Grundsätzlich können aus den Eigenwerten der Systemmatrix A folgende drei zentralen Informationen zum Systemverhalten abgeleitet werden:
Insgesamt spricht man hier von der Eigenbewegung des Systems.
Ein Gleichungssystem gilt als stabil, wenn bei fortschreitender Zeit alle Zustandsvariablen einem stabilen Endwert zustreben. Wächst eine Zustandsvariable unendlich an, spricht man von Instabilität.
Anhand des Real-Teils eines Eigenwerts können wir auf die Stabilität schließen: Sind die Real-Teile aller Eigenwerte negativ, dann ist das System stabil. Ist der Real-Teil von mindestens einem Eigenwert positiv, dann ist das System instabil.
Anhand des imaginären Teils der Eigenwerte, können Informationen zum Schwingungsverhalten des Systems abgeleitet werden. Gibt es keinen Imaginär-Teil (die Eigenwerte liegen also auf der Re-Achse), dann bedeutet dies, dass das System nicht schwingt. Liegen Eigenwerte hingegen im imaginären Bereich, so bedeutet dies, dass es im System zu Schwingungen kommt.
Zuletzt sind durch die Eigenwerte Informationen über die Zeitskalen im System bekannt. Die Zeitkonstante ist ein Maß für die Änderungsgeschwindigkeiten von Zuständen.
Bei einem nicht schwingenden System entspricht die Zeitkonstante dem Kehrwert des negativen Real-Teils der Eigenwerte.
Bei einem schwingenden System (also einem System mit Imaginär-Teil) wird der negative Kehrwert des Real-Teils als Dämpfung bezeichnet und gibt Auskunft, über die Zeitskala des Schwingungs-Abklingens. Die Kreisfrequenz der Schwingung entspricht dem Imaginär-Teil des Eigenwerts und lässt sich einfach durch Teilung durch 2π in die Schwingungsfrequenz umformen.
Wichtig ist zu beachten, dass diese Ableitungen von Systemverhalten nur für lineare Systeme gelten, da nur dann die Systemmatrix über die gesamte Simulationsdauer konstant ist.
Bei nicht-linearen Gleichungssystemen ist die Systemmatrix A zeitlich veränderlich, was dementsprechend auch zu veränderlichen Zeitkonstanten, Schwingungen und Stabilitäten führen kann.
Möchte man das Systemverhalten von nicht-linearen Systemen über die Eigenwerte analysieren, ist dies trotzdem möglich: Für einzelne Zeitpunkte kann das jeweilige System linearisiert und dann entsprechend in diesem Moment analysiert werden.
Doch wie schließen wir jetzt von diesen Informationen auf Steifheit?
Für Steifheit gibt es keine einheitliche Definition. In diesem Artikel werden zwei unterschiedliche Definitionen von Steifheit beleuchtet: Die eine bezieht sich nur auf das Gleichungssystem selbst und die andere auf das Zusammenspiel von Gleichungssystem und Lösungsalgorithmus.
Die erste Definition haben wir bereits anhand des einführenden Beispiels und der mathematischen Hintergründe erläutert. Ein System gilt als steif, wenn die Größenordnung der minimalen und maximalen Zeitkonstante \(\tau\) stark variiert. Grundsätzlich kann man sagen:
Je stärker die Größenordnung der Zeitkonstanten variieren, desto steifer ist das System.
Eine klare Abgrenzung, ab wann ein System steif ist, gibt es in dieser Definition nicht. Mathematisch ausgedrückt bedeutet das:
$$\frac{\tau_{Max}}{\tau_{Min}} >>1$$
Die Steifheit ist also allein vom Gleichungssystem abhängig.
Wir wissen nun, wie sich Aussagen über das Systemverhalten durch die Interpretation der Systemmatrix A treffen lassen. Kehren wir zurück zu unserem Beispiel und vergleichen die Größenordnungen der Zeitkonstanten, zeigt sich Folgendes:
$$\frac{\tau_{Max}}{\tau_{Min}}=\frac{500.000}{1}>>1 1$$
Die Größenordnung der Zeitkonstanten im Beispiel variiert also sehr stark. Nach der ersten Definition liegt ein steifes System vor. Es lässt sich grundsätzlich folgendes sagen:
Die maximale Zeitkonstante bestimmt, wie lange die langsame Kapazität zum Aufheizen benötigt und damit die Simulationsdauer. Die minimale Zeitkonstante legt die erforderliche Schrittweite des Solvers für korrekte Ergebnisse fest. Ein großer Quotient der beiden verlängert die Simulationszeit.
Im Falle eines so einfachen Systems ist es noch möglich die Zeitkonstanten und damit die Steifheit selbst zu prüfen – doch was macht man bei deutlich komplexeren Simulationsmodellen mit vielen Zuständen und Zeitkonstanten?
Dymola bietet ein praktisches Tool zum Debuggen: Full-Linear-Analysis
Zunächst führen wir die Analyse mit den Default-Einstellung für unser Beispielsystem durch. In einem Pop-Up werden nun die Eigenwerte unseres Systems im komplexen Zustandsraum dargestellt:
Eigenwert-Analyse unseres Beispielsystems:
Zusätzlich werden im Commands-Window weitere Informationen zur Verfügung gestellt. Angezeigt werden die Aufschlüsselung unseres Systems in die bereits erwähnte Zustandsraumdarstellung, eine Aussage über die Systemstabilität und eine Eigenwertanalyse – also all die Informationen, die wir aus der Systemmatrix A ziehen können. Für unser Beispielsystem sieht dieEigenwert Analyse wie folgt aus:
Wir haben zwei Eigenwerte, die jeweils zu 100% den jeweiligen Temperatur-Zuständen der Kapazitäten zuzuordnen sind. Beide Eigenwerte haben einen negativen Real-Teil (das System ist also stabil) und keinen Imaginär-Teil (das System schwingt also nicht). Außerdem haben die Zeitkonstanten die Werte 1s und 500.000s, was einem Quotienten von 500.000 entspricht. Die schnelle Kapazität ändert ihre Temperatur also im Bereich von Sekunden, während die langsame Kapazität eine Dauer von mehreren Tagen für eine vergleichbare Änderung benötigt – dies ist also ein gutes Beispiel für die Steifheit eines Systems.
Nun stellt sich die Frage, wie man die Recheneffizienz eines solchen Systems steigern kann. Ein Ansatz ist die Reduktion der Steifheit – wir müssen also den Quotienten der Zeitkonstanten verringern.
Dazu gibt es zwei Optionen:
Beide Anpassungen würden zu einer Veränderung des physikalischen Systemverhaltens führen, aber:
Eine Erhöhung des großen Widerstands (und damit der großen Zeitkonstante) beeinflusst die Simulationsdauer bis zum stationären Zustand erheblich. Beim kleinen Widerstand ist das anders: Aufgrund seiner schnellen Dynamik bleibt die Temperatur der schnellen Kapazität nahezu identisch zur Quelltemperatur (z.B. Ofen). Eine moderate Erhöhung des kleinen Widerstands verlangsamt zwar die Kopplung, verändert jedoch das Gesamtverhalten nur minimal. Daher eignet sich diese Anpassung gut zur Reduzierung der Steifheit.
Lösung: Der kleine Widerstand wird um den Faktor 20 auf \(R_{Nachher}=0,2 K/W\) angehoben.
Die Simulation erfolgt mit den gleichen Einstellungen und führt zu folgendem Ergebnis: Nach erneut 646h erreichen beide Kapazitäten die Zieltemperatur. Der Solver benötigt dafür rund 1,4s und 585.293 Schritte. Die Simulationsdauer konnte also etwa halbiert werden, während die Schrittanzahl weniger als ein Drittel beträgt.
Die Simulation ist also bei nahezu gleichbleibendem Ergebnis deutlich effizienter geworden.
Ein Blick in die Full-Linear-Analysis zeigt, dass der Quotient der Zeitkonstanten in diesem neuen Fall nur noch bei 25.000 liegt und somit ebenfalls um den Faktor 20 reduziert wurde. Während diese Zusammenhänge in einem so simplen System noch recht klar sind, können die Eigenwerte in komplexen Simulationsmodellen deutlich schwerer zu berechnen sein, was die Full-Linear-Analysis zu einem sinnvollen Tool beim Debugging von Simulationsmodellen macht.
Wie wir wissen, ist die Analyse der Eigenwerte nur für lineare Systeme möglich, da die Systemmatrix in nicht-linearen Systemen zeitlich veränderlich ist. Doch die Full-Linear-Analysis von Dymola bietet auch die Option, nicht-lineare Systeme um einen bestimmten Zeitpunkt zu linearisieren und dann zu analysieren. Dafür wird im Tab „simulationSetup“ der Zeitpunkt für die Linearisierung ausgewählt. Anschließend wird die Analyse durchgeführt und die Eigenwerte können interpretiert werden.
Doch das Gleichungssystem ist nicht der einzige Faktor, der die Steifheit einer Simulation beeinflussen kann. Das Zusammenspiel von Solver und Gleichungssystem kann ebenfalls zu steifen Systemen führen. Das führt uns zur zweiten Steifheitsdefinition. Doch zunächst schauen wir uns einmal genauer an, was die Solver (Lösungsalgorithmen) überhaupt tun.
In numerischen Simulationen werden die Zustandsgrößen eines Systems nicht direkt berechnet, sondern über die Zeit integriert. Dabei kommen spezielle Lösungsalgorithmen, sogenannte numerische Solver, zum Einsatz. Sie ermitteln die Zustände im nächsten Zeitschritt, indem sie die zugrunde liegenden Differentialgleichungen lösen, was oft über iterative Lösungen der Gleichungssysteme passiert.
Eine der ersten grundlegenden Unterscheidungen zwischen verschiedenen Solvern betrifft die Schrittweite: Einige Algorithmen arbeiten mit einer festen Schrittweite, während andere diese dynamisch anpassen, um Genauigkeit und Effizienz zu optimieren. Besonders Solver mit fester Schrittweite können in steifen Systemen schnell an ihre Grenzen stoßen – warum das so ist, zeigt sich gut an unserem nächsten Beispiel.
Wir betrachten unser Beispielmodell mit zwei parallelen thermischen RC-Gliedern. Anhand der bisherigen Erkenntnisse des Blogartikels gilt dieses System als steif, aber stabil, da die Eigenwerte nur negativen Real-Teil haben. Nun simulieren wir das System mit dem expliziten Euler-Verfahren mit einer festen Schrittweite von 3s.
Entgegen der Erwartung ist das System scheinbar instabil. Die Temperatur der ersten Kapazität schwingt sich stark zu völlig unrealistischen Werten auf, obwohl wir laut Eigenwerten gar keine Schwingung erwarten, da es keinen Imaginär-Teil gibt.
Das Problem entsteht durch das Zusammenspiel von Solver und Zeitkonstanten, wodurch das System numerisch instabil wird. Die feste Schrittweite des Solvers von 3s überschreitet die kleinste Zeitkonstante von 1s. Dadurch kann der Solver die schnelle Dynamik des Systems nicht präzise erfassen, was zu ungenauen oder sogar völlig fehlerhaften Ergebnissen führt.
Um dieses Problem in den Griff zu bekommen, ändern wir die Solver-Schrittweite auf beispielsweise 0.5s, also einen Wert kleiner als die kleinste Zeitkonstante. Durch diese Anpassung erhalten wir das erwartete stabile Ergebnis.
Die Stabilität eines Lösungsalgorithmus ist dementsprechend vom Zusammenspiel von Schrittweite und Eigenwerten abhängig und unterscheidet sich von Solver zu Solver.
Um die Stabilität von Solver grafisch darzustellen, werden sogenannte Stabilitätsregionen definiert.
Stabilitätsregionen sind in der komplexen Ebene angeordnet und zeigen den Bereich der Stabilität eines Solvers in Abhängigkeit der Eigenwerte und Schrittweiten. Dafür werden die jeweiligen Eigenwerte λ jeweils mit der Schrittweite h multipliziert und in ein Stabilitätsdiagramm eingetragen. Ein Solver gilt für ein Gleichungssystem als stabil, wenn sich alle Eigenwerte λ nach Skalierung mit h innerhalb der Stabilitätsregion des Solvers befinden. Schauen wir auf das Beispiel von vorhin. Die Stabilitätsregion des expliziten Euler-Verfahrens ist ein Kreis mit dem Radius 1 um den Punkt -1 auf der Real-Achse.
Im ersten Fall beträgt die Schrittweite h = 3s und die Real-Teile der Eigenwerte liegen bei -1 und -2e-06, während die Imaginär-Teile beide 0 sind. Die beiden im Diagramm einzutragenden Punkte liegen dementsprechend bei -3 und -6e-6 auf der Real-Achse. Damit liegt der Punkt bei -3 auf der x-Achse außerhalb der Stabilitätsregion, was sich im instabilen Verhalten zeigt.
Nach der Anpassung liegt der kritische Eigenwert nach Multiplikation mit der neuen Schrittweite bei \(-0.5\cdot 1=-0.5\) und damit innerhalb der Stabilitätsregion. Das Ergebnis ist folglich stabil.
Der Solver muss also die Integrationsschritte mit einer geringeren Schrittweite wiederholen, wenn er die vorgegebene Toleranz nicht einhalten kann. Der lokale Fehler bremst die Simulation aus. Aber auch die Bedingung zur Einhaltung der Stabilitätsregionen kann den Solver ausbremsen und dazu zwingen, andere Schrittweiten zu wählen. Und damit kommen wir zur zweiten Definition der Steifheit:
Ein System gilt als steif, wenn der Solver für die Gewährleistung der Stabilität zur Anpassung der Schrittweite gezwungen ist. Die Steifheit ist also (im Gegensatz zur ersten Definition) vom Zusammenspiel aus Gleichungssystem und Solver abhängig.
Moderne Solver mit Schrittweitenanpassung
Durch große Stabilitätsregionen und die Möglichkeit der modernen Solver ihre Schrittweite zur Einhaltung der Stabilitätsgrenzen anzupassen, kommt es seltener zu Problemen mit der Stabilität (solange das Gleichungssystem selbst auch stabil ist). Doch auch wenn Steifheit nicht zwingend zu Instabilität führt, kann sie auch bei modernen Solvern zu großen Effizienzproblemen führen.
Dymola bietet ein praktisches Tool zum Debuggen: Simulation Analysis
Neben der Full-Linear-Analysis stellen wir noch ein weiteres Tool vor, das das Debugging steifer Systeme erleichtert: Die Simulationsanalyse mit der Ausgabe, welche Zustände die Fehler dominieren und zur Begrenzung der Schrittweite führen.
Über den Tab „Simulation Analysis“ und dem gesetzten Haken bei „Which states that dominate error“ erhalten wir eine hilfreiche Übersicht.
Für alle im System enthaltenen Zustandsvariablen wird angegeben:
Mithilfe dieses Tools erhalten wir Informationen über den Zustand, der den Solver am ehesten ausbremst und können dann die entsprechenden Parameter (und damit die Eigenwerte) anpassen, um die Effizienz zu steigern.
Betrachten wir das in einem Beispiel: Wir führen die Simulationsanalyse für unser Beispielsystem mit den zwei parallelen RC-Gliedern und den stark unterschiedlichen Zeitkonstanten durch. Das Ergebnis für die komplette Simulationsdauer von 650h sieht wie folgt aus:
Der Temperaturzustand der schnellen Kapazität limitiert die Solver-Schrittweite während der Simulation insgesamt über 77.000-mal. Aus der vollständigen linearen Analyse ist uns bekannt, dass die kleine Zeitkonstante von 1s zu dieser Zustandsvariable gehört. Das schnelle Verhalten der Variable sorgt dafür, dass der Solver die Schrittweite immer wieder nach unten korrigieren muss – also bereits getätigte Schritte abbricht und neu integrieren muss. Das kostet die Simulation viel Zeit.
Erhöhen wir die Zeitkonstante – beispielsweise um den Faktor 50 – durch die Vergrößerung des Wärmewiderstands und simulieren erneut, zeigt uns die Analyse folgendes Ergebnis:
Die Schrittweite wird nur noch rund 15.000-mal reduziert, was immer noch viel ist, aber diesen Punkt trotzdem auf etwa 20% reduziert.
Diese Anpassung führt, wie wir im ersten Beispiel bereits gesehen haben, zu einer großen Reduktion der Simulationszeit und damit zu einer starken Effizienzsteigerung.
Die Simulationsanalyse zur Schrittweitenlimitierung stellt somit ein gutes ergänzendes Tool dar, um in umfangreichen Modellen die problematischen Zustände zu identifizieren. Zusammen mit der vollständigen linearen Analyse ist dann eine Anpassung der relevanten Parameter (und damit der Eigenwerte) möglich, um die Systemeffizienz zu steigern
Die beiden gezeigten Definitionen der Steifheit haben ihre Daseinsberechtigung. Die Steifheit eines Simulationsmodells kann allein vom Gleichungssystem verursacht sein. Ebenso ist eine ungünstige Kombination von Eigenwerten des Gleichungssystems und dem gewählten Solver eine Möglichkeit, durch die steife Systeme entstehen und ineffizient oder sogar instabil werden.
Die in Dymola vorgestellten Tools ermöglichen ein besseres Verständnis der Systeme. Dadurch kann ein optimales Zusammenspiel von großen und kleinen Zeitkonstanten sowie dem passenden Solver gefunden werden, um ein effizient rechnendes, stabiles System zu erzielen, das die Anforderungen an die Simulationsergebnisse erfüllt.
© 2025 TLK Energy. Alle Rechte vorbehalten.
