Software: FEM - Tutorial - Magnetfeld - Probabilistik - Kennfeld-Identifikation

↑

← →

Identifikation von Ersatzmodellen

Die Kraft- und Koppelfluss-Kennfelder eines E-Magneten kann man in Modellen der System-Simulation als Ersatzmodell für eine konkrete Wandler-Geometrie benutzen:

1. 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.
2. 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 nun konstante Größen. Diese Streuungsmittenwerte sind so zu wählen, dass über den gesamten Streubereich numerisch instabile Werte nahe Null vermieden werden:
  • i.Wert=5.01 A (i≥0.01 A)
  • s.Wert=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
  • i_.Nennwert und s_.Nennwert Null (Toleranzmittenabstand=0)
• 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 x₁ und x₂ im Abstand (x₁−x₂):

1. Square Exponential
   
2. Exponential
   
3. Gamma-Exponential
   
4. Matern Class 3/2
   
5. Matern Class 5/2
   
6. Periodic
   
7. 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, für dessen Konfiguration die Standardparameter=True nutzen. Die in der Versuchsplanung zu wählende Optimierungsmethode für die Robustoptimierung 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 besten Polynom-Ordnung für jede Streuungsrichtung.
• Bei der Option "Best Covariance" versucht OptiY automatisch die beste Covarianz-Funktion zu ermitteln, was inzwischen recht gute Ergebnisse liefert:

===>>> ab hier erfolgt zur Zeit eine Überarbeitung der Anleitung!!!

• Nach dem Neuberechnen der Antwortflächen ergibt sich eine Approximation nur für das F-Kennfeld, während als Psi-Kennfeld nur das Polynom der gewählten Ordnung=1 erscheint:
  
• Optimierungsmethode=Evolutionsstrategien in der Versuchsplanung führt nach dem Neuberechnen der Antwortflächen zu einem alternativen Ergebnis:
  

Hier existiert also ein grundsätzliches Konvergenz-Problem bei der Suche nach der optimalen Antwortfläche auf Basis eines Gauß-Prozesses:

• Versuche haben für das Gradienten-Verfahren gezeigt, dass der Wertebereich der Antwortfläche eine entscheidende Rolle spielt
  • F    = 0 ... 100 (Verfahren konvergiert)
  • Psi = 0 ... 0,1 (Verfahren konvergiert nicht!)
• Es wird vermutet, dass hier das Rauschen des Gauß-Prozesses bei kleinen Absolutwerten eine Rolle spielt.
• Workaround (bis ein OptiY-Update dieses Problem mindert):
  • Das FEMM-Programm berechnet den Koppelfluss Psi in [Wb]
  • Im Workflow multiplizieren wir diesen Wert mit dem 1000 und erhalten damit Psi in [mWb] → Einheit für Gütekriterium korrigieren:
    
  • Danach sollte das Gradientenverfahren problemlos funktionieren:
    
• Positiv im Vergleich mit dem Polynom-Ansatz fällt beim Gauß-Prozess auf, dass:

1. physikalisch korrekt keine negativen Kräfte oder Koppelflüsse berechnet werden und
2. die Residuen zeigen, dass sämtliche Punkte der echten Stichprobe exakt (max. Abweichung ca. 1E-13) in die identifizierte Ersatzfunktionen F=f(s,i) bzw. Psi=f(s,i) eingebettet sind:
   

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 identifizierten Kennfelder besitzen auf Grund der Polynomordnung=0 für die globale Anpassung schwache Knicke an den Stützstellen. Hier würde man physikalisch eine stetige Ableitung der Kennfelder erwarten.
• Die 95% Vertrauensintervalle zwischen den Stützstellen sind sehr eng:
  

Insbesondere den Aspekt des Vertrauensintervalls sollten wir näher betrachten:

• Um die obige Darstellung der Vertrauensintervalle zu erhalten, müssen wir:

1. die Anzeige der Stützpunkte für die Schnittdiagramme aktivieren.
2. 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).

• An den Stützstellen selbst ist die Breite des Vertrauensintervalls immer Null, da der Wert bekannt ist.

Hinweis:

• 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.

Um ein Gefühl für die Wirkung der Polynomordnung und der unterschiedlichen Covarianz-Funktionen auf die Qualität der Antwortflächen zu erhalten, werden diese Parameter im Folgenden systematisch verändert:

1. Polynomordnung

• Wir erhöhen die Polynomordnung für beide Ergebnisgrößen schrittweise um 1.
• Nach dem Neuberechnen der Antwortflächen ergibt sich im Beispiel eine Verschlechterung des Ersatzmodells in Hinblick auf das physikalische Verständnis:
  

1. Für die Magnetkraft entstehen wie beim Polynomansatz negative Werte.
2. Die Welligkeit beider Antwortflächen nimmt mit steigender Ordnungszahl zu.

Man sollte die kleinste Polynomordnung verwenden, welche zu physikalisch plausiblen Ersatzmodellen führt! Das wird meist die Polynomordnung=0 sein.

2. Covarianz-Funktion

Wir setzen die Polynomordnung auf den mit dem einfachen Exponential-Ansatz ermittelten optimalen Wert, bevor wir mit den unterschiedlichen Covarianz-Funktionen experimentieren:

• Höherwertige Covarianz-Funktionen führen bei unseren monotonen, stetigen Kennfeldern zu minderwertigen Ergebnissen im Vergleich zum einfachen Exponential-Ansatz:
  

1. Es wird zwar ein schmaleres Vertrauensintervall berechnet, was aber nicht mit dem physikalischen Verständnis des Magnetmodells korreliert!
2. Die Kennfelder erhalten eine verstärkte Welligkeit, welche z.B. bei der Magnetkraft zu physikalisch sinnlosen negativen Kraftwerten führen kann.

• Je komplizierter die gewählte Covarianz-Funktion ist, desto größer ist die Wahrscheinlichkeit, dass die globale Approximation für den Gauß-Prozess misslingt.

Hinweise

Die beschriebene Vorgehensweise für die Identifikation von Ersatzmodellen als Antwortflächen mittels Gauß-Prozess ist nur unter den folgenden Randbedingungen durchführbar:

1. Es existieren bereits hinreichend viele Stützstellen für das zu identifizierende Modell innerhalb des vorgesehenen Wertebereiches der freien Parameter.
2. Die Anzahl der freien Parameter ist überschaubar (z.B. maximal drei).
3. 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:
  1. Knicke bewirken Sprünge der Ableitung
  2. 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!

← →