Modelle von Satelliten-Komponenten an reale Testergebnisse angleichen

Satelliten kommunizieren untereinander per Laser. Die Elektronik muss dabei extreme Temperaturen und Temperaturwechsel verkraften.courtesy of ESA - P. Carril

Satelliten unterliegen extremen Temperaturschwankungen: Auf der einen Seite erhitzt die Sonne ihre Oberfläche auf über 200 °C, während die andere Seite fast den absoluten Nullpunkt von -273,15 °C erreicht. Eine Dreiviertelstunde später kann es umgekehrt sein. Auch die robustesten weltraumqualifizierten Bauteile widerstehen solchen Schwankungen nicht dauerhaft. Daher ist einer der entscheidenden Designtreiber für Satelliten die Thermalkontrolle der Elektronik. Eine Erprobung unter realen Bedingungen ist vor Flug nicht möglich und Tests, die eine Weltraumumgebung teilweise nachbilden, sind teuer und aufwändig. Folglich reduzieren Satellitenhersteller diese Tests auf ein Minimum und nutzen numerische Modelle für den Großteil der Auslegungen, Nachweise und Vorhersagen.

Da viele wichtige Entscheidungen zum Design und dem Betrieb des Systems auf diesen Modellen beruhen, müssen sie eine ausreichende Genauigkeit und Zuverlässigkeit aufweisen. Um dies sicherzustellen werden standardmäßig diese Modelle mit den Umwelttests verglichen. Hier zeigen sich so gut wie immer Abweichungen zwischen den Vorhersagen eines Modells und den zugehörigen Messungen. Mit feineren Netzen, einer genaueren physikalischen Modellierung und besseren Algorithmen kann man die Abweichungen zwischen Modell und Realität nur teilweise reduzieren. Ein bedeutender Teil der Abweichung ist darauf zurückzuführen, dass bestimmte Parameter nicht hinreichend genau bekannt sind.

Ein typisches Beispiel für solch einen Parameter ist die Wärmeleitfähigkeit der Wärmeisolation des Satelliten. Die verwendete Superisolation (MLI) besteht aus mehreren aluminumbedampften Polyamidfolien, welche die Infrarotstrahlung reflektieren. Das gelbbraune Polyamid der äußeresten Folie gibt den meisten Satelliten ihre goldene Farbe. Unter idealen Laborbedingungen ist die Wärmeleitfähigkeit durch diese Isolation recht gut bekannt, aber in der Realität hat die Folie Knicke, Falten, Löcher und Ränder, welche Wärmebrücken darstellen. Somit variiert die Wärmeleitfähigkeit durch diese Isolationen um Größenordnungen von Aufbau zu Aufbau.

Modellkorrelation statt Einzelmessungen

Idealerweise würde man jeden Modellparameter einzeln vermessen, dies ist allerdings meistens viel zu aufwändig oder schlicht unmöglich. Daher ist es üblich, Modelle mit den Umwelttests zu korrelieren. Dafür werden die Modellparameter so angepasst, dass die Ergebnisse mit der Messung möglichst gut übereinstimmen.

Um diese Korrelation zu automatisieren, haben verschiedene Forschungsgruppen Methoden entworfen und untersucht (Amundsen 2003, Jouffroy 2007, De Palo 2011, Momayez 2009, Harvatine 1994, Roscher 2006, van Zijl 2013, WenLong 2011, Mareschi 2005). Die meisten dieser Methoden basieren auf modernen stochastischen Optimierungsalgorithmen und brauchen mehrere hunderte bis zehntausende Iterationen bis zur Konvergenz. In diesem Artikel wird eine Klasse von Algorithmen untersucht, welche Charles George Broyden schon 1965 in den Anfängen der Informatik definierte. Ziel ist die Korrelation von thermalen Modellen aus der Raumfahrt mit realen Messungen; die Methode eignet sich aber auch für viele andere Modelle.

Mathematische Formulierung

Aus mathematischer Sicht ist ein thermales Modell ein Temperaturvektor, welcher eine Funktion eines Parametervektors ist:

• t_mdl = F(p)

Dabei gilt:

• F ist das mathematische Modell. Es ist eine Funktion F: R^k ⟶ R^m (k ist die Anzahl an Parametern, m ist die Anzahl an Temperaturen)
• t_mdl ist ein Vektor mit m Temperaturen. Die Temperaturen sind durch eine Position an der Geometrie und einem Zeitpunkt definiert.
• p ist der Vektor mit den Parametern für das Modell wie Wärmeleitfähigkeiten, Wärmekapazitäten, Emissionsgrade und solare Absorptionkoeffizienten.

Das Ziel einer Korrelation lautet, einen Parametersatz p_corr zu finden, der den Betrag des Differenzvektors zwischen dem Vektor der gemessenen Temperaturen t_mes und dem entsprechenden Ergebnisvektor F(p_corr) minimiert. Dieser Differenzvektor (F(p_corr) – t_mes) heißt im Folgenden Abweichungsvektor. Mathematisch ausgedrückt ergibt sich:

• ‖(F(p_corr) – t_mes)‖ = min{‖F(p) – t_mes ‖ | p∈P}

mit P als Lösungsraum für den Parametervektor p.

Die verwendeten Methoden

Im Gegensatz zu den meisten bisher entwickelten Ansätzen minimiert diese Methode nicht direkt den Betrag des Abweichungsvektors ‖F(p) – t_mes‖ , sondern versucht das folgende Gleichungssystem zu lösen:

• F(p_corr) – t_mes = 0

Der Vorteil dieses Ansatzes ist, dass er nicht nur einen Skalar-Wert beurteilt sondern einen ganzen Vektor an Informationen verwendet. Dies führt zu einem Gleichungssystem mit relativ einfachen Funktionen statt einer einzigen, aber sehr komplexen Funktion. Eine der besten Methoden um solch ein Problem zu lösen ist das multidimensionale Newton-Verfahren. Dieses löst interaktiv folgendes Gleichungssystem:

• p_(n+1) = p_n – J(p_n)^-1 × F(p_n)

Dabei bedeuten:

• p_n ist der Parametervektor.
• J(p_n) ist die Jacobimatrix der (Modell-)Funktion F(p_n) an der Stelle p_n.

Der Nachteil dieser Methode ist, dass man die Jacobimatrix oft nur numerisch berechnen kann: Für jeden Parameter p_n ist eine Berechnung von F(p_n) nötig. Daher entwickelte Broyden eine Klasse von Methoden, um die Änderung der Jacobimatrix anhand der Ergebnisse der vorhergehenden Iteration abzuschätzen (Broyden 1965). Diese Methoden sind in der Literatur und in Broydens Veröffentlichung beschrieben. Für diese Untersuchung an Satelliten eignet sich die erste Methode nach Broyden, welche auch als „Broydens gute Methode“ bekannt ist. Zusätzlich wird eine vom Autor dieses Artikels entwickelte Methode verwendet, welche die Elemente (b_n^ij) der approximierten Jacobimatrix (B≈J) mit dem größten Einfluss in der letzten Iteration am stärksten verändert. Die Aktualisierung erfolgt anhand folgender Formel:

• b_n+1^ij = (1 + k_nⁱ · b_n^ij · s_n^j) · b_n^ij

mit:

• s_n := p_n – p_n-1
• y_ex,n := B_n · s_n
• y_is,n := F(p_n) – F(p_n-1)
• kⁱ = (y_is,nⁱ – y_ex,nⁱ) / (∑_j=1^m (b_n^ij · s_n^j)²)

Monotonie, Differenzierbarkeit

Jede mathematische Methode hat bestimmte Grenzen und Randbedingung, unter denen sie funktioniert und effizient arbeitet. Nachfolgend werden einige Randbedingungen erwähnt, die sich im praktischen Einsatz als relevant erwiesen. Diese Randbedingungen sollten Anwender schon bei der Wahl der Methode, der Parameter und der Zielgrößen (Sensorpositionen und Zeiten) beachten.

Eine wichtige Einschränkung von Broydens guter Methode ist, dass die Funktionen im betrachteten Parameterraum monoton und differenzierbar sein müssen. Bei thermalen Modellen ist das meist der Fall; Ausnahmen wären zum Beispiel geregelte Systeme oder wenn einer der Parameter die Position im Orbit ist. Auch durch numerisches Rauschen und Runden sind numerische Modelle bei kleinen Parameteränderungen weder monoton noch differenzierbar. In diesen Fällen konvergiert der Algorithmus oft nicht.

Beobachtbarkeit und Beeinflussbarkeit

Jeder Parameter muss einen relevanten Einfluss auf mindestens ein Ergebnis ausüben. Diese Bedingung gilt nicht nur für diese Methode, sie ist allgemein nötig für eine Korrelation. Auch wenn diese Bedingung trivial erscheint, ist sie doch oft der Grund dafür, dass ein Parameter nicht konvergiert. Spezifisch für die beschriebene Methode können schlecht beobachtbare Parameter mit numerischem Rauschen zur Instabilität führen.

Die Relevanz eines Parameters lässt sich oft aus dem Wissen über das Modell abschätzen. In komplexeren Systemen ist das aber oft schwierig und aufwändig. Hier kann die erstellte Jacobimatrix Hinweise geben: Wenn eine Spalte der Jacobimatrix nur aus Nullen besteht, hat der entsprechende Parameter keine Relevanz. In praktischen Anwendungen kommen eher sehr kleine Werte statt Nullen vor. Sind diese Werte kleiner als die Messgenauigkeit des entsprechenden Tests, können sie als nicht beobachtbar gelten.

Jedes Ergebnis muss von mindestens einem Parameter beeinflusst werden. Dies ist ebenfalls eine allgemeine Bedingung für Korrelationen. Kann ein Ergebnis von keinem Parameter beeinflusst werden, dann kann dies durch numerisches Rauschen zu Instabilitäten führen. Auch hierfür gibt die Jacobimatrix wertvolle Hinweise: Besteht eine Zeile nur aus Nullen oder sehr kleinen Werten, dann ist das entsprechende Ergebnis nicht oder kaum beeinflussbar.

Durch Ungenauigkeiten in der Messung und in der Modellierung wird der Abweichungsvektor so gut wie nie gegen Null konvergieren. Sind diese Ungenauigkeiten zu groß, kann es vorkommen, dass das Modell nicht innerhalb des physikalisch sinnvollen Parametraums zu sinnvollen Ergebnissen konvergiert.

Grenzen des Parameterraums

Fast jeder Parameter hat sinnvolle physikalische Grenzen. Beispielsweise wird eine Wärmeleitfähigkeit nie negativ sein und auch nie einen bestimmten Wert überschreiten. Da Broydens Methode nur für unbeschränkte Probleme anwendbar ist, wird folgender Algorithmus benutzt um dies zu umgehen:

• Berechne p_n+1 ohne die Grenzen des Parameterraums zu berücksichtigen.
• Liegt ein Parameter von p_n+1 außerhalb des erlaubten Parameterraums, dann skaliere die Änderung der Parameter so, dass alle Parameter innerhalb des Parameterraums liegen.
• Wenn einige der Parameter bereits an den Grenzen des Parameterraums waren, dann berechne p_n+1 für die Bedingung, dass diese Parameter auf deren Grenze fixiert sind.

Broydens Algorithmus wurde für diesen Beitrag an zwei Modellen ausprobiert: an einem einfachen imaginären Beispiel und an einem komplexen Modell eines realen Systems. Diese Modelle wurden in der Berechnungssoftware Systema/Thermica simuliert. Diese raumfahrtspezifische Software ähnelt den bekannteren Produkten Esatan und Sinda. Die Beschreibung der Modelle orientiert sich an den üblichen Definitionen dieser Software.

Bild 1: Das Schema des einfachen Modells zeigt vier Knoten und als Randbedingung eine Umgebung mit 0 °C.Tesat-Spacecom

Das einfache Modell

Das einfache Modell besteht aus einer stationären Simulation von vier veränderbaren Knoten und einem Umgebungsknoten bei 0 °C als Randbedingung. Jeder Knoten hat eine Fläche von 0,1 m² mit einer Emissivität von 1 und strahlt mit einem Sichtfaktor von 1 gegen den Umgebungsknoten. 10 W werden in jedem Knoten dissipiert. Knoten 1 ist mit dem Umgebungsknoten mit einem Leitwert von 1 W/K verbunden. Knoten 1 bis 4 sind untereinander mit sechs Leitwerten GL₁ bis GL₆ verbunden. Bild 1 zeigt den Aufbau.

Die sechs Leitwerte sind die Eingangsparameter des Modells, welche so verändert werden sollen, dass die Temperaturen der Knoten einer angenommenen Messung entsprechen. Die angenommenen Messergebnisse entsprechen den Temperaturen für die Leitwerte: GL_i = (0,1 + i/100) W/K. Drei Fälle werden mit diesem Modell untersucht:

• Ein unterbestimmtes System: Alle Leitwerte sind variabel. Es gibt daher sechs Parameter (Leitwerte) für vier  Ergebnisse (Temperaturen)
• Ein bestimmtes System: Leitwerte 3 und 6 sind festgesetzt (GL₃ = 0,13 W/K und GL₆ = 0,16 W/K). Es gibt daher vier Parameter (Leitwerte) für vier Ergebnisse (Temperaturen)
• Ein überbestimmtes System: Leitwerte 3, 5 und 6 sind festgesetzt (GL₃ = 0,13 W/K, GL₅ = 0,5 W/K und GL₆ = 0,16 W/K). Es gibt daher drei Parameter (Leitwerte) für vier Ergebnisse (Temperaturen). Der Leitwert GL₅ ist so gewählt, dass es unmöglich ist das Gleichungssystem zu lösen.

Für alle Leitwerte, die nicht festgesetzt sind, beträgt der Anfangswert 0,5 W/K.

Ergebnisse des unterbestimmten Falls

Bild 2 zeigt die Entwicklung des Betrags des Abweichungsvektors (RSS) für beide Algorithmen. Ein paar relevante Punkte dieses Verlaufs sind in Tabelle 1 zusammengetragen. Der Algorithmus konvergierte schnell gegen Null für beide Algorithmen. Pro Größenordnung benötigte der Algorithmus zwei bis drei Iterationen bis als Grenze die numerische Auflösung der Ergebnisse von zirka 10^-5 K erreicht wurde.

Tabelle 1: Ausgewählte Punkte der Konvergenz des Algorithmus für das unterbestimmte Modell.Tesat-Spacecom

Bild 2: Betrag des Abweichungsvektors über den Iterationen des unterbestimmten Models.Tesat-Spacecom

Bis zur vierten Iteration verhielten sich die beiden Algorithmen fast identisch. Bei der fünften und sechsten Iteration war der selbst entwickelte Algorithmus schneller, wurde dann aber vom Broyden-Algorithmus eingeholt. Es lässt sich folgern, dass sich beide Algorithmen für diesen Fall eignen.

Ergebnisse des bestimmten Falls

Auf gleiche Art und Weise, wie schon für das unterbestimmte Modell gezeigt, zeigt Bild 3 die Entwicklung des Betrags des Abweichungsvektors und Tabelle 2 die Werte ausgewählter Punkte des Diagramms. Im Mittel brauchte die Methode nach Broyden zirka elf Iterationen pro Größenordnung wobei der selbstentwickelte Algorithmus im Mittel nur fünf benötigte.

Tabelle 2: Ausgewählte Punkte der Konvergenz des Algorithmus für das bestimmte Modell.Tesat-Spacecom

Bild 3: Betrag des Abweichungsvektors über den Iterationen des bestimmten Modells.Tesat-Spacecom

Ergebnisse des überbestimmten Falls

Das überbestimmte Modell konvergierte wie erwartet nicht gegen Null sondern gegen einen endlichen Wert (0,0375 K). Dies kann man in Bild 4 und Tabelle 3 sehen. Für diesen Fall konvergierte die selbst entwickelte Methode mit zirka vier Iterationen pro Größenordnung deutlich schneller als Broydens Methode mit zirka acht Iterationen.

Tabelle 3: Ausgewählte Punkte der Konvergenz des Algorithmus für das überbestimmte Modell.Tesat-Spacecom

Bild 4: Betrag des Abweichungsvektors über den Iterationen des überbestimmten Modells.Tesat-Spacecom

Das komplexe Modell

Das komplexe Modell ist eine transiente Simulation eines realen Raumfahrtsystems in einer Simulation der thermalen Weltraumumgebung in einer Thermovakuumkammer. Insgesamt wurden 13 Parameter für die Korrelation ausgesucht. Dies waren thermale Randbedingungen, Leitwerte und Emissivitäten. Die Temperaturen von 26 Sensoren für sechs Zeitschritte wurden ausgewertet. Damit ergeben sich 156 verwendete Messwerte.

Bild 5 und Tabelle 4 zeigen in einer leicht abgewandelten Form die Konvergenz des komplexen Modells. Da es sich hier um ein überbestimmtes Modell handelt, konvergiert es nicht zu Null. Daher wurde die auf den Anfangswert normierte Differenz des Abweichungvektorsbetrags zum niedrigsten erreichten Abweichungsvektorbetrag dargestellt: (RSS – RSS_min) / (RSS_initial – RSS_min). Der niedrigste erreichte Wert ergab sich nach einem Neustart des Algorithmus’ nahe dem Optimum.

Tabelle 4: Ausgewählte Punkte der Konvergenz des Algorithmus für das komplexe Modell.Tesat-Spacecom

Bild 5: Betrag des Abweichungsvektors über den Iterationen des komplexen Modells.Tesat-Spacecom

In Bild 5 zeigt sich, dass der Algorithmus schon nach wenigen Schritten ein Optimum fast erreicht. Mit nur fünf Iterationen reduzierte Broydens Algorithmus den Zielwert um fast drei Größenordnungen. Der selbst entwickelte Algorithmus benötigte sogar eine Iteration weniger. Bei zirka 0,22 % blieb Broydens Algorithmus stabil, während der selbstentwickelte Algorithmus sich instabiler verhielt aber letztendlich zu einem geringeren Wert (0,06 %) konvergierte.

Zusammenfassung der Ergebnisse

Sowohl Broydens als auch der selbst entwickelte Algorithmus konvergierte für alle untersuchten Fälle. Der Wert gegen den die Algorithmen konvergierten war entweder die numerische Auflösung der Ergebnisse oder eine Grenze infolge der Überbestimmtheit des Systems. Die Anzahl der Iterationen kann mit der folgenden Formel approximiert werden:

• m = 1 + k + r · c

mit:

• m ist die Gesamtanzahl an Iterationen inklusive der vorangehenden Berechnung mit den Anfangswerten.
• k ist die Anzahl der Parameter (Iterationen für die Berechnung der ersten Jacobimatrix).
• r ist ein problemspezifischer Wert welcher hauptsächlich von den Nichtlinearitäten und der Interaktion der Parameter abhängt. Die Anzahl an Parametern hat keinen Einfluss auf diesen Parameter. Für die untersuchten Fälle lag dieser Wert zwischen 2 und 10 bis ein Grenzwert erreicht wurde. Für ähnliche Problemstellungen ist nicht zu erwarten, dass dieser Wert deutlich größer wird.
• c ist die anvisierte Reduktion der Abweichung zwischen Modell und der Messung.

Der Wert c ist definiert über die folgende Formel:

• c = log₁₀(R₀ – R_min) – log₁₀(R_final – R_min)

mit:

• R₀ ist der Anfangsbetrag des Abweichungsvektors (F(p₀) – t_mes)
• R_final ist der Endbetrag des Abweichungsvektors (F(p_final) – t_mes).
• R_min ist der kleinste mögliche Betrag des Abweichungsvektors (F(p_ideal) – t_mes).

Besonders für Modelle mit vielen Parametern ist es interessant, dass die Anzahl an Parametern, im Gegensatz zu den stochastischen Methoden, höchstens linear in diesen Zusammenhang eingeht.

Bild 6: Vergleich der Anzahl an benötigten Iterationen für eine Modellkorrelation (n = Anzahl der Modelknoten, p = Anzahl der korrelierten Parameter).Tesat-Spacecom

Schneller am Ziel

Wie in Bild 6 zu sehen ist, genügen deutlich weniger Iterationen (10 bis 30) für diese Untersuchung als bei den Untersuchungen mit genetischen Algorithmen (25600: van Zijl 2013, 843 bis 33072: Jouffroy 2007) oder Schwarmalgorithmen (6000: van Zijl 2013). Leider sind weder die Modelle noch die Algorithmen aus den fremden Untersuchungen verfügbar, sodass ein direkter Vergleich nicht möglich ist. Die Problemstellung ist aber ähnlich und die Anzahl an Parametern ist vergleichbar mit dem komplexen Modell. Die Anzahl an Knoten ist deutlich kleiner, da die Autoren üblicherweise eine Korrelation großer Modelle als zu aufwändig oder gar unmöglich ansehen.

Im Beitrag wurde nur die Anzahl an Iterationen verglichen, da die benötigte Rechenzeit um einen neuen Satz an Parametern zu generieren mit wenigen Millisekunden vernachlässigbar klein ist gegenüber den Minuten oder Stunden für das Lösen des Modells. Daher ist die Rechenzeit fast proportional zur Anzahl an Iterationen, unabhängig davon welcher Algorithmus zum Einsatz kommt. Und diese Anzahl liegt für die untersuchte Problemstellung bei numerischen Verfahren aus Broydens Klasse typischerweise 20- bis 1000-mal niedriger als bei stochastischen Methoden.

(lei)