Software: FEM - Tutorial - Magnetfeld - Probabilistik - Kennfeld-Identifikation
Die Kraft- und Koppelfluss-Kennfelder eines E-Magneten kann man in Modellen der System-Simulation als Ersatzmodell für eine konkrete Wandler-Geometrie benutzen:
- Die Ergebnisse einer FEM-Rastersuche kann man als Daten exportieren und die berechneten Abtaststellen des FEM-Modells als Stützstellen von 3D-Funktionsflächen in der System-Simulation verwenden. Aus diesen Stützstellen werden in der System-Simulation dann die Funktionswerte der Zwischenräume interpoliert.
- Die für die Durchführung der virtuellen Stichprobe von OptiY generierten mathematischen Funktionen der Antwortflächen kann man aber auch direkt als Modell-Code in Modelle der System-Simulation implementieren. Die Gewinnung von Ersatzmodellen nach diesem Prinzip wollen wir abschließend in diesem Übungskomplex durchführen.
Ausgehend vom Raster-Experiment Kennfeld-Berechnung konfigurieren wir nach dem Duplizieren ein Experiment Kennfeld-Identifikation.
Workflow
- Wir werden die probabilistische Simulation benutzen, um über den gesamten Arbeitsbereich unseres E-Magneten die Antwortflächen für F(i,s) und Psi(i,s) bilden zu lassen:
- Die Nennwerte von Strom und Arbeitsluftspalt sind durch Streuungen zu ersetzen. Die Streuungsmittenwerte sind so zu wählen, dass über den gesamten Streubereich numerisch instabile Werte nahe Null vermieden werden:
- i.Nennwert=5.01 A (i≥0.01 A)
- s.Nennwert=2.03 mm (s≥30 µm)
- Die Streuungen von Strom und Arbeitsluftspalt umfassen den gesamten Arbeitsbereich des E-Magneten. Die Gleichverteilung innerhalb der Streu-Bereiche gewährleistet eine gleichmäßige Abtastung, falls man eine Sample Methode verwendet:
- i.Toleranz=10 A
- s.Toleranz= 4 mm
- Die beiden anderen Luftspalte setzen wir manuell in der Input-Datei auf ihre Nennwerte:
- sDeckel = 20 µm
- sGleit = 80 µm
Versuchsplanung
Wir werden dieses Simulationsexperiment (Optimierungsverfahren=Simulation) so konfigurieren, dass wir möglichst anschaulich erkennen, wie genau die approximierten Antwortflächen die wirklichen Übertragungsfunktionen des Modells abbilden:
- Grundlage für die Berechnung der Antwortfläche soll das Full Factorial Design sein, welches praktisch einer Rastersuche entspricht. Wir wählen dafür 6 Stufen, um die Anzahl der Modellberechnungen mit 36 gering zu halten. Die Anzahl 36 ergibt sich aus dem Raster 6x6.
- Wir setzen den virtuellen Stichprobenumfang=0, da wir in diesem Experiment nur die Antwortflächen gewinnen.
Approximation mit Polynom-Ansatz
- Wir verwenden zuerst die bereits für die Sample Methode genutzte polynomiale Approximation für die Gütekriterien und beginnen mit der einheitlichen Polynomordnung=1.
- Nach Abschluss der Simulation liegen die Antwortflächen als Funktionen im (i,s)-Raum vor.
- In die Diagramme der 3D-Antwortflächen blenden wir die Punkte der berechneten Stützstellen ein.
- Schrittweises Erhöhen der Polynomordnung mit anschließendem Neuberechnen der Antwortflächen führt zu einer optimalen Anpassung der Antwortflächen. Aus einem zu großem Wert für die Polynomordnung (z.B. ab Ordnung=6) resultiert dann wieder eine Verschlechterung der Approximation durch extrem starke Welligkeit der Fläche. Für unsere beiden Kennfelder führt eine manuelle, individuelle Anpassung der Polynomordnung an die Abszissen-Richtungen (Strom und Luftspalt) zu keiner weiteren Verbesserung.
Achtung: Ein Rücksetzen des Experiments mit erneuter Berechnung der realen Stichprobe ist beim Ermitteln der optimalen Polynomordnung nicht erforderlich!
- Ein Drehen des Diagramms der 3D-Antwortfläche zeigt deutlich die Lage der realen Modellwerte abseits der Antwortfläche. Besonders stark sind die Abweichungen bei der Kombination von minimalem Luftspalt mit minimalem Strom.
- Das Residuum-Plot zeigt quantitativ die Abweichungen der realen Modellberechnungen von den identifizierten Antwortflächen. Bei der Magnetkraft liegt diese Abweichung überwiegend im Bereich von ca. 10%, was für ein Ersatzmodell häufig noch akzeptiert werden kann.
- Nicht akzeptabel ist jedoch, dass die Ersatzfunktion im Bereich kleiner Ströme zu negativen Kraftwerten führt, was physikalisch nicht korrekt ist!
- Ähnlich ungenau sehen die Ergebnisse für das Psi-Kennfeld aus:
Schlussfolgerung:
Polynom-Ansätze sind für die Bildung von globalen Ersatzmodellen bei nichtlinearem Übertragungsverhalten nicht besonders gut geeignet.
Approximation mit Gauss-Prozess
Der Gauß-Prozess, angewandt in der Geostatistik auch als Kriging bekannt, ist ein statistisches Verfahren, mit dem man Werte an Orten, für die keine Probe vorliegt, durch umliegende Messwerte interpolieren oder auch annähern kann.
Der Gauß-Prozess besteht aus einem globalen Modell f(x) und einem stochastischen Prozess Z(x), welcher die mögliche Abweichung von dem globalen Modell beschreibt:
- x ist ein m-dimensionaler Parametervektor.
- Y(x) ist der Ergebnisvektor für den Punkt x im Parameterraum.
- f(x) sind Polynome beliebiger Ordnung, welche zusammen mit den unbekannten Regressionskoeffizienten βi die Regressionsfunktion bilden.
- Z(x) ist ein stationärer stochastischer Prozess mit dem Mittelwert Null, der Varianz σ und der Covarianz R. Dieser Anteil des Gauß-Prozesses beschreibt das 95% Erwartungsintervall für jeden Punkt x des Parameterraumes.
In dem englischen Wikipedia-Artikel zum Kriging wird das für die Interpolation einer eindimensionalen Funktion sehr anschaulich dargestellt:
Die berechneten Grenzverläufe des Erwartungsintervalls werden wesentlich bestimmt durch das Erfahrungswissen in Hinblick auf den erwartenden Verlauf der zu interpolierenden Funktion zwischen den bekannten Werten der Stichproben-Exemplare. Diese Erwartung wird durch die Wahl einer geeigneten Covarianz-Funktion R beschrieben.
- In OptiY sind unterschiedlichste Covarianz-Funktionen R implementiert. Diese beschreiben den Verlauf des 95% Erwartungsintervalls zwischen den Stützstellen in Abhängigkeit von der Stützstellendichte. Die folgende Notation bezieht sich auf zwei Stützstellen x1 und x2 im Abstand (x1−x2):
- Square Exponential
- Exponential
- Gamma-Exponential
- Matern Class 3/2
- Matern Class 5/2
- Periodic
- Rational Quadratic
- Der allgemeine Fall für Exponential-Funktionen ist Gamma-Exponential, die anderen beiden Exponentialfunktionen sind die Spezialfälle für γ=1 bzw. γ=2.
- Die Matern Class Funktionen sind Erweiterungen der Exponential-Funktion.
- Die Hyper-Parameter w, γ , λ und α werden mittels der Maximierung der Likelihood-Funktion der multivariaten Normalverteilung ermittelt.
Die Ermittlung dieser Hyper-Parameter des Gauß-Prozesses erfolgt in der statistischen Versuchsplanung mittels eines implementierten Optimierungsverfahrens, in dessen Konfiguration wir die Standardparameter nutzen. Die in der Versuchsplanung zu wählende Optimierungsmethode für das "Robust Design" ist davon völlig unabhängig:
Über die Polynomordnung des globalen Modells kann man einen Kompromiss finden zwischen bester Anpassung der interpolierten Regressionsfunktion an vorhandene Stützstellen und optimalem Verlauf zwischen diesen Stützstellen:
- Wie im vorherigen Abschnitt "Approximation mit Polynom-Ansatz" beschrieben, bestimmt die Polynomordnung des globalen Modells f(x) die allgemeine Richtung (globale Anpassung) der Regressionsfunktion. Die dabei verbleibenden Residuen sind noch sehr groß.
- Der stochastische Prozess Z(x) hat dann die Aufgabe, diese verbleibenden Residuen mittels der Covarianz-Funktion (Normal-Verteilung) zu eliminieren (lokale Anpassung).
- Wenn die Anzahl der Stützstellen bzw. der Daten hinreichend groß ist und man die verbleibenden Residuen statistisch auswertet, entsteht eine Normal-Verteilung (auch Gauss-Verteilung genannt) mit dem Mittelwert=0. Das ist die ursprüngliche Idee des Gauß-Prozesses.
- Wenn man die globale Anpassung mit einer zu hohen Ordnungen der Polynome durchführt, besitzt die Kurve mehr Freiheitsgrade als nötig. Das führt dann zu Welligkeiten zwischen den Stützstellen, weil dies durch keine Zwangsbedingungen verhindert wird.
Dieses Prinzip der Bildung von Ersatzmodellen (Antwortflächen) werden wir nun auf unser Kennfeld-Problem anwenden. Wir nutzen dafür das bereits konfigurierte Experiment Kennfeld-Identifikation ohne erneute Berechnung der Stichprobe:
- Wir wählen für die beiden Ergebnisgrößen F und Psi zur Approximation den Gauss Prozess und überlassen dem Programm für den Polynom-Typ die Ermittlung der Beste Ordnung für die Polynome in jede Streuungsrichtung.
- Mit der Option Best Covariance versucht OptiY automatisch die beste Covarianz-Funktion zu ermitteln, was inzwischen recht gute Ergebnisse liefert:
- Nach dem Neuberechnen der Antwortflächen ergeben sich folgende Approximationen für das F-Kennfeld und das Psi-Kennfeld:
- Die Residuen zeigen, dass die Punkte der echten Stichprobe etwas genauer in die identifizierte Ersatzfunktionen F=f(s,i) bzw. Psi=f(s,i) eingebettet sind, als beim reinen Polynom-Ansatz:
- Hinweis: In Versionen vor OptiY 4.6 kam es zu Problemen bei der Residuenberechnung, welche eine höhere Genauigkeit darstellten, als mit den Antwortflächen eigentlich erreicht wurde.
- Die jetzt korrekt berechneten Residuen decken sich in etwa mit der 3D-Darstellung der Antwortflächen. Trotzdem ist zu beachten::
- Die dargestellten Flächen beruhen auf einer Approximation der wirklichen 3D-Funktionen.
- Die Genauigkeit der Darstellung wird über die Zahl der Rasterpunkte gesteuert (Standardwert=20). Dabei kann jedoch die Genauigkeit der Darstellung nicht besser werden, als die zugrundeliegende Genauigkeit der Approximation.
- Nur zur Info - muss in der Übung nicht ausprobiert werden:
- Um den Effekt hier in dieser Anleitung zu verdeutlichen, wurde innerhalb des Workflows der Psi-Wert mit dem Faktor 1000 multipliziert (Psi-Werte also in mWb).
- Nach erneuter Berechnung der Stichprobe ergibt sich dann eine wesentlich bessere Approximation der darzustellenden Fläche (hier dargestellt mit dem Raster 100x100).
- Dabei wird deutlich, dass die Funktionsfläche des Gauss-Prozesses so "zurecht gezupft" wurde, dass die Stützstellen wirklich exakt auf dieser Fläche liegen.
- Dadurch entstehen "Unstetigkeiten" (hier in Form schmaler Rinnen bzw. Grate), welche in der Realität nicht vorhanden sind:
Schlussfolgerung:
- Man muss die mittels automatisiertem Gauss-Prozess erzeugten Antwortflächen einer gründlichen Analyse in Hinblick auf unzulässige Welligkeiten und Unstetigkeiten untersuchen.
- Bei Bedarf muss man den Gauss-Prozess dann doch manuell optimieren:
- Verwendung einer möglichst kleinen Polynomordnung.
- Wahl einer geeigneten Covarianzfunktion zur Erzeugung realistischer Übergänge zwischen den Stützstellen der realen Stichprobe.
Zusätzlich zu den Residuen und den 3D-Antwortflächen öffnen wir für jede Ergebnisgröße auch die Schnittdiagramme. Innerhalb dieser Schnittdiagramme aktivieren wir die Darstellung des Vertrauensintervalls:
- Die automatisch identifizierten Kennfelder besitzen Welligkeiten zwischen den Stützstellen. Hier würde man physikalisch monotone Verläufe der Kennfelder erwarten.
- Die 95% Vertrauensintervalle zwischen den Stützstellen sind recht eng:
Insbesondere den Aspekt des Vertrauensintervalls sollten wir näher betrachten:
- Um die obige Darstellung der Vertrauensintervalle zu erhalten, müssen wir:
- die Anzeige der Stützpunkte für die Schnittdiagramme aktivieren.
- die aktuellen Werte der Streugrößen (Virtueller Entwurf - Nennwert) auf eine Stützstellen-Reihe setzen (z.B. _i=6.01 A und _s=2.43 mm).
Wichtig:
- An den Stützstellen selbst ist die Breite des Vertrauensintervalls immer Null, da der Wert bekannt ist.
- Die berechneten Vertrauensintervalle resultieren nur aus den benutzten mathematischen Ansätzen. Entscheidend für uns sind jedoch die interpolierten Funktionsverläufe von Kraft und Koppelfluss! Wir wissen (im Unterschied zum mathematischen Formalismus), dass sich diese physikalischen Größen zwischen den Stützstellen monoton ändern (ohne "Welligkeit").
- Die Vertrauensintervalle zeigen den Worst Case der Unsicherheit bei fehlenden Kenntnissen zu den abgebildeten physikalischen Zusammenhängen.
Hinweise
Die beschriebene Vorgehensweise für die Identifikation von Ersatzmodellen als Antwortflächen mittels Gauß-Prozess ist nur unter den folgenden Randbedingungen durchführbar:
- Es existieren bereits hinreichend viele Stützstellen für das zu identifizierende Modell innerhalb des vorgesehenen Wertebereiches der freien Parameter.
- Die Anzahl der freien Parameter ist überschaubar (z.B. maximal drei).
- Man besitzt bereits qualitative Vorstellungen zu den Eigenschaften der Antwortflächen.
Die Identifikation von Kennfeldern für einen elektro-magnetischen Wandler ist nicht ganz einfach:
- Entscheidend für die Genauigkeit des damit aufzubauenden Dynamik-Modells ist die hinreichend genaue zeitliche Ableitung des Psi-Kennfeldes. Die damit berechnete Induktionsspannung bestimmt die zeitliche Entwicklung des Stromes in der Spule.
- Der Gauss-Prozess berücksichtigt nur, wie exakt die berechneten Stützstellen in die Antwortfläche passen und wie breit das Vertrauensintervall der identifizierten Antwortfläche zwischen den berechneten Stützstellen ist. Knicke und Welligkeiten der Antwortfläche fließen nicht in die Bewertung der Antwortflächen-Güte ein, sind aber entscheidend für die Güte der benötigten Ableitung!
- Die Güte der Ableitung eines Kennfeldes kann man mittels OptiY zurzeit nur iterativ durch die Wahl einer geeigneten Covarianz-Funktion in Kombination mit einer günstigen Polynomordnung beeinflussen.
Qualität der Ableitung:
- Benötigt man innerhalb der Systemsimulation partielle Ableitungen der identifizierten Funktionen, so muss man dies als zusätzliches Qualitätsmerkmal berücksichtigen. Im Beispiel trifft dies für Psi(i,s) zu, denn die zeitliche Ableitung des Koppelflusses entspricht im Netzwerkmodell der induzierten Spannung in der Spule.
- Leider kann man die Qualität partieller Ableitungen der identifizierten Antwortflächen innerhalb des OptiY nur sehr grob beurteilen:
- Knicke bewirken Sprünge der Ableitung
- Welligkeit kann zu Vorzeichenwechsel der Ableitung führen
- Um die physikalische Sinnfälligkeit der Ableitung zu beurteilen, muss man zurzeit diese Ableitung auf Basis des exportierten Modell-Codes berechnen und darstellen!