Dieser Beitrag bezieht sich auf PGAPack und PGApy oder andere Genetische Algorithmus (GA) Implementierungen die Differential Evolution (vielleicht übersetzbar mit Differentielle Evolution, in der Folge abgekürzt mit DE) unterstützen.

Differential Evolution (DE) [1], [2], [3] ist ein Optimierungsalgorithmus der ähnlich wie andere evolutionäre Algorithmen auf einer Population basiert. Der Algorithmus ist recht mächtig, in PGAPack implementiert und auch in PGAPy nutzbar. Um einige Punkte zu diesem Algorithmus zu veranschaulichen habe ich einen einfachen Parcours in OpenSCAD konstruiert.

Um den Optimierer mit diesem Parcours zu testen ist das Gen (DE nennt es die Parameter) die X- und Y-Koordinate um eine Position im Parcours zu bestimmen. In der Folge bezeichnen wir diese zwei Werte auch als Vektor, wie das in der Literatur zu DE üblich ist. Wir initialisieren die Population in der Gegend beim Start der Rampe. Wir erlauben der Population, den Initialisierungsbereich zu verlassen.

Es ist zu beachten, dass die Ecken des Parcours flach sind, dass dort also alle Punkte die selbe Höhe haben. Solche Bereiche sind generell schwierig für Optimierungsalgorithmen weil der Algorithmus nicht wissen kann, in welche Richtung bessere Werte (wenn es diese überhaupt gibt) zu finden sind. Dies ist auch der Grund, warum wir den Initialisierungsbereich nicht in den flachen Bereich am Boden des Parcours legen. Üblicherweise würde man nämlich die Population über den gesamten zu durchsuchenden Bereich initialisieren. In diesem Fall möchten wir demonstrieren, dass der Algorithmus auch fähig ist, eine Lösung weit vom Initialisierungsbereich zu finden.

Der DE Algorithmus ist recht einfach (siehe z.B. [3] S.37-47): Mit jedem Vektor \(\mathbf x_{i,g}\) der Population der aktuellen Generation \(g\), wobei \(i\) der Populationsindex ist wird ein Genaustausch mit einem Mutanten-Vektor vorgenommen. Der Mutanten-Vektor wird zufällig (ohne Zurücklegen) aus drei Populations-Vektoren bestimmt, indem diese wie folgt kombiniert werden.

Wie man sieht, wird eine gewichtete Differenz von zwei Vektoren der Population gebildet und zu einem dritten Vektor addiert. Alle Indices \(r_0, r_1,\) und \(r_2\) sind voneinander und von \(i\) verschieden. Das Gewicht \(F\) ist ein Konfigurationsparameter der typischerweise \(0.3 \le F < 1\) ist. Der beschriebene Algorithmus ist der klassische DE Algorithmus, oft auch als de-rand-1-bin bezeichnet. Eine Variante nutzt statt einem zufälligen Index \(r_0\) den Index des besten Vektors und wird als de-best-1-bin bezeichnet. Falls mehr als zwei Differenzen addiert werden wäre eine andere Variante z.B. de-rand-2-bin. Die letzte Komponente der Bezeichnung steht für Uniform Crossover [4], eine binomiale Verteilung. In DE werden die Parameter des Mutanten-Vektors mit einer bestimmten Wahrscheinlichkeit \(Cr\) ausgewählt. Man beachte, dass DE dafür sorgt, dass zumindest ein Parameter des Mutanten-Vektors ausgewählt wird (sonst bliebe der betrachtete Original-Vektor \(\mathbf x_{i,g}\) unverändert).

Mehr Details zu DE und er Nutzung mit PGAPack findet sich in den Abschnitten Differential Evolution und Mutation der PGAPack Dokumentation.

Um das OpenSCAD Modell in der Optimierung zu verwenden, wird es mithilfe eines STL zu PNG Konverters in ein Graustufen-Höhenbild umgewandelt (STL ist ein 3D-Format das von OpenSCAD generiert werden kann). Die Bewertungsfunktion des Algorithmus liefert dann einfach den Graustufen-Wert an der gegebenen X/Y-Position zurück (nach Rundung der Gen-Werte und Konvertierung in einen Integer-Wert). Werte ausserhalb des Bildes werden als 0 (schwarz) geliefert. Der resultierende Optimierungslauf ist unten in einem animierten Bild dargestellt. Die Population ist ziemlich klein und enthält nur 6 Vektoren. Wir sehen dass, wenn die Abstände zwischen den Punkten groß werden (auf den geraden Strecken des Parcours), die Optimierung in größeren Schritten voranschreitet. Wir sehen auch, dass der Algorithmus in den flachen Ecken einige Zeit brauchen kann bis er aus diesem Bereich entkommt. Schließlich, in der letzten Ecke, wird der Kegel erklommen und der Algorithmus stoppt wenn das erste Individuum eine Bewertung mit dem größten Graustufen-Wert liefert (entspricht weiss). Bitte beachten, dass ich ein bisschen geschummelt habe: Viele Optimierungläufe dauern deutlich länger (Insbesondere in den flachen Bereichen kann der Algorithmus längere Zeit stecken bleiben). Ich habe einen Lauf ausgewählt der sich für eine Animation gut eignet. Dieser Lauf ist nicht im Bereich des Durchschnitts sondern deutlich kürzer.

Dieser Beitrag beschäftigt sich wieder mit PGAPack und PGApy, sowie anderen Implementierungen von Genetischen Algorithmen (GA).

Bei der Suche von Lösungen mit einem GA passiert es, dass eine sub-optimale Lösung gefunden wird, weil die Population frühzeitig in einem Bereich des Suchraumes konvergiert wo keine besseren Lösungen mehr gefunden werden. Dieser Effekt heisst vorzeitige Konvergenz.

Es ist normalerweise schwierig, vorzeitige Konvergenz zu entdecken. Eine gute Metrik ist die mittlere Hamming Distanz zwischen Individuen. In PGAPack kann man Auswertungen der Hamming Distanz über die Option PGA_REPORT_HAMMING der Funktion PGASetPrintOptions aktivieren. In PGApy geht das mit dem Konstruktor Parameter print_options. Leider ist dies nur für den binären Datentyp implementiert.

Ein möglicher Grund für das Auftreten von vorzeitiger Konvergenz ist die Verwendung von proportionaler Selektion wie in einem früheren Artikel in diesem Blog [1] ausgeführt. Wenn während der Frühphase der Suche ein Individuum gefunden wird das viel besser ist als alles bisherige, so ist die Chance groß, dass dieses Individuum die Population übernimmt wenn proportionale Selektion verwendet wird, wodurch ein weiterer Fortschritt des Algorithmus verhindert wird. Der Grund dafür ist, dass ein Individuum das wesentlich besser ist als alle anderen einen unangemessen hohen Anteil des (gedachten) Roulette-Rades einnimmt wenn Individuen für die nächste Generation ausgewählt werden, was dazu führt, dass alles andere genetische Material nur eine geringe Chance hat, in die nächste Generation zu kommen.

Ein anderer Grund für vorzeitige Konvergenz kann eine zu kleine Populationsgröße sein. Goldberg et. al. [9] geben Abschätzungen der Populationsgröße für einen klassischen GA mit kleinem Alphabet (mögliche verschiedene Allele, also 0 und 1 für einen binären GA) der Kardinalität \(\chi\) an. Dabei wird eine maximale Länge einer Symbolfolge (Building Block) angenommen die irreführende Symbolfolgen (deception) überwindet mit \(k \ll n\) wobei \(k\) die Symbolfolgen-Länge (ein Maß für die Schwierigkeit des Problems) und \(n\) die Gen-Länge ist. Es wird gezeigt, dass die Populationsgröße \(O(m\chi^k)\) sein sollte mit \(m=n/k\) so dass mit einer vorgegebenen Schwierigkeit \(k\) die Populationsgröße proportional zur Gen-Länge \(n\) ist. Dieses Ergebnis kann aber nicht auf Probleme mit einem sehr großen Alphabet wie Fließkomma-Genen wie dem real Datentyp von PGAPack verallgemeinert werden. Für Fließkomma Darstellungen kann man die Schwierigkeit üblicherweise einschätzen wenn man die Multimodalität, also die Anzahl der Gipfel oder Täler (jenachdem ob es sich um ein Maximierungs- oder Minimierungsproblem handelt) der Zielfunktion betrachtet.

Wenn mit einer angemessenen Populationsgröße und einem sinnvollen Selektionsmechanismus noch immer vorzeitige Konvergenz auftritt, können einige andere Mechanismen zur Anwendung kommen.

Einige Leser hier wissen dass ich PGApack, ein Paket für Genetische Algorithmen, zusammen mit PGApy einer Python Bibliothek pflege (maintaine). Vor kurzem kam die Frage auf, warum PGApack wenn eine Selektionsmethode verwendet wird die proportional zur Fitness selektiert (Fitness Proportionate Selection), im englischen auch als Roulette-Rad Selektion bezeichnet, mit einem Fehler aussteigt weil es einen Fitness-Wert nicht normalisieren kann. Im folgenden werde ich kurz von proportionaler Selektion schreiben.

In PGApack definieren wir als Anwenderin eine Reelle Bewertungsfunktion (die für jedes Individuum aufgerufen wird). Ausserdem definieren wir ob diese Bewertung minimiert oder maximiert wird (wir legen die Optimierungsrichtung fest).

Bei proportionaler Selektion muss diese Bewertung auf positive Werte abgebildet werden so dass jeder Bewertung ein Teil des Roulette-Rades zugewiesen wird. Wenn minimiert wird sind kleine Werte besser und müssen einen größeren Bereich des Roulette-Rades bekommen so dass die Wahrscheinlichkeit der Selektion höher ist. Wir müssen also die Bewertungsfunktion auf eine nicht-negative monoton steigende Fitness abbilden. Bei einem Minimierungsproblem berechnen wir die maximale (schlechteste) Bewertung und bilden dann die Differenz dieses Maximums mit der Bewertung des Individuums (nach einer Skalierung des Maximums so dass kein Fitnesswert negativ wird):

wobei \(F\) die Fitness des aktuellen Individuums ist, \(E_{max}\) das Maximum aller Bewertungen der aktuellen Generation und \(E\) die Bewertung des Individuums.

Wenn sich nun die Bewertungen um mehrere Größenordnungen unterscheiden kann es vorkommen, dass die Fitness zu \(E_{max}\) wird, für viele (unterschiedliche) Individuen. Ich nenne das einen Überlauf (overflow), vielleicht nicht der beste Name dafür.

Der Überlauf passiert wenn \(E_{max}\) groß ist verglichen mit der aktuellen Bewertung \(E\) so dass die obige Differenz \(E_{max}\) wird (die Subtraktion von \(E\) also keine Wirkung hat) Im Code checken wir diese Bedingung:

if (cmax == cmax - evalue && evalue != 0)

Die Bedingung ist erfüllt wenn die Subtraktion von \(E\) von \(E_{max}\) \(E_{max}\) nicht verändert obwohl \(E\) ungleich Null ist. \(E\) ist so klein gegenüber \(E_{max}\) dass der Datentyp double den Unterschied nicht repräsentieren kann. Das passiert wenn die Einheit der letzten Stelle (units in the last place von Goldberg (nicht dem Goldberg der das Genetische Algorithmen Buch geschrieben hat soweit ich weiss) auch ulps genannt) der Mantisse größer ist als der Wert der subtrahiert wird [1].

In unserem Beispiel \(E_{max} = 1.077688 * 10^{22}\) und die Bewertung bei der das schiefging war \(E = 10000\). Ein IEEE 754 Fließkomma-Wert doppelter Genauigkeit (double) hat 53 bit Mantisse die damit Zahlen bis \(2^{54} - 1 = 18014398509481983\) darstellen kann, also etwa \(1.8 * 10^{16}\). Man sieht dass die Zahl 10000 gerade unterhalb des oben erwähnten ulps ist. Wir können das in Python ausprobieren (das für Fließkommazahlen double verwendet):

>>> 1.077688e22 - 10000 == 1.077688e22
True

Warum machen wir diesen Check im Programm? Würde die Suche bei einem solchen Überlauf (oder wie auch immer wir das nennen wollen) weiterlaufen, würde der Algorithmus viele verschiedene Bewertungen auf die gleiche Fitness abbilden. Der Genetische Algorithmus könnte also diese Individuen nicht unterscheiden.

Was können wir tun wenn dieser Fehler auftritt?

Die kurze Antwort: Einen anderen Selektionsmechanismus verwenden. Es gibt einen Grund, warum proportionale Selektion nicht die Standardeinstellung in PGApack ist.

Proportionale Selektion (Fitness proportionate selection) hat noch andere Probleme. Sie hat zu großen Selektionsdruck am Anfang zu wenig gegen Ende der Optimierung (auch erwähnt im obigen Wikipedia Artikel, aber Vorsicht, ich habe Teile davon geschrieben).

Blickle and Thiele [2] haben eine mathematische Analyse von unterschiedlichen Selektions-Algorithmen publiziert und gezeigt dass proportionale Selektion typischerweise keine gute Idee ist (proportionale Selektion war historisch der erste Selektionsmechanismus in der Literatur und war ausführlich in Goldbergs (der andere Goldberg) Buch [3] beschrieben und ist vielleicht deshalb heute noch manchmal in Verwendung). Es sei darauf hingewiese dass Blickle und Thiele in einem früheren Report [4] direkter waren was die Beurteilung von proportionaler Selektion betrifft (meine Übersetzung): "Alle unerwünschten Eigenschaften zusammen führten uns zu dem Schluss dass proportionale Selektion ein sehr ungeeigneter Selektionsmechanismus ist. Informell könnten wir sagen dass der einzige Vorteil von proportionaler Selektion ist, dass es so schwierig ist die Nachteile zu beweisen" ([4], S. 42), in der Endfassung im veröffentlichten Artikel waren sie dann nicht mehr so deutlich :-)

Wir sehen in obigem Beispiel: Wir haben sehr große Unterschiede zwischen guten und schlechten Bewertungen (eben so groß dass die Fitness nicht eindeutig berechnet werden kann, s.o.). Wenn man dafür proportionale Selektion verwendet, werden sehr gute Individuen mit zu großer Wahrscheinlichkeit selektiert was zu verfrühter Konvergenz (premature convergence) führt.

Zum Schluss all dieser Ausführungen: Wenn wir Optimierungen mit Fließkomma-Repräsentationen durchführen (diese werden durch double Datentypen in PGApack repräsentiert) sollten wir Differential Evolution ausprobieren [5], [6], [7]. Zumindest in meinen Experimenten mit Antennen-Optimierung [8] sind die Resultate damit deutlich besser als mit standard Genetischen Algorithmen was auch von verschiedenen Praktikern bestätigt wird [7]. Beispiele dazu finden sich in examples/deb/optimize.c oder examples/mgh/testprogde.c in PGApack.

Der Parameter PGASetDECrossoverProb ist kritisch. Für Probleme wo die Dimensionen nicht einzeln optimierbar sind, sollte der Wert nahe bei, aber nicht gleich 1 sein.

Für mein pymininec Projekt (eine Re-Implementierung des Original-Mininec [1] in Python, siehe die Beschreibung dort wie es möglich war, 40 Jahre alten Basic-code zu Reverse-Engineeren) hatte ich eine Möglichkeit gebraucht, die Resultate graphisch darzustellen. Daher hatte ich zuerst ein Visualisierungswerkzeug in pymininec inkludiert.

Nachdem ich dann aber auf die Schnelle einen Parser für die Ausgabe von nec2c geschrieben hatte, stellte sich heraus, dass das Programm auch allein nützlich ist. Daher habe ich es als eigenes Projekt plot-antenna genannt.

NEC (vielleicht übersetzbar mit numerischer elektromagnetischer Code) wurde ursprünglich am Lawrence Livermore National Laboratory in Fortran entwickelt. Bis Version 2 ist das freie Software, spätere Versionen sind nicht frei (inzwischen ist NEC bei Version 5). Neoklis Kyriazis hat die Fortran Version von NEC-2 nach C übersetzt und das Resultat nec2c genannt. Er hat einige Optimierungen vorgenommen, z.B. das ursprünglich separate Programm SOMNEC (um ein Sommerfeld/Norton Erdungsmodell zu berechnen) integriert und einige Limitierungen des ursprünglichen Fortran-Codes beseitigt. Das Programm nec2c ist eines von wenigen Open Source Antennen-Modellierungsprogrammen. Es ist ein Kommandozeilen-Programm das aus der Beschreibung einer Antenne eine Ausgabedatei mit vielen Tabellen berechnet. Es gibt andere Programme wie xnecview die mit dieser Ausgabedatei eine Antennen-Visualisierung erzeugen können.

An dieser Stelle ist zu bemerken, dass ich keine offizielle Webseite von nec2c gefunden habe. Das Programm ist in Debian und Ubuntu (Ubuntu sucht auch nach der offiziellen Upstream-Version) enthalten und es scheint die ursprüngliche Version war mal auf Google Code gehostet. Es gibt eine Version auf Github die auf eine andere Version auf Github verweist. Die Software scheint heutzutage nicht gewartet zu werden – kein großes Problem weil auch NEC in Version 2 seit ein paar Dekaden keine Updates mehr bekommt.

Zurück zu plot-antenna: Ursprünglich hatte das Programm nur einen Parser für die Ausgabe von pymininec (und für das Original Mininec [1] Basic-Programm). Später habe ich dann einen Parser für NEC hinzugefügt. Die derzeitigen Parser sind wahrscheinlich nicht sehr robust, funktionieren aber für meine Zwecke.

Die erste Version verwendet Matplotlib für die Grafikausgabe und diese Variante ist immer noch unterstützt. Urprünglich hatte ich Azimuth- und Elevationsdiagramme, sowie eine 3D-Anzeige implementiert.

Kürzlich, nachdem ich einige Zeit vorher Plotly entdeckt hatte, habe ich HTML/Javascript Ausgabe mithilfe von Plotly hinzugefügt. Die 3D Möglichkeiten von Plotly sind deutlich besser als die von Matplotlib – sogar mit 1 Grad Auflösung bei den Antennen-Winkeln (Azimuth/Elevation), also 180 * 360 3D Punkten, funktioniert das Display im Browser (Firefox). Mit der 3D Variante von Matplotlib muss 5 Grad Auflösung verwendet werden (36 * 72 3D Punkte), sonst ist die Anzeige sehr ruckelig. Mit der folgenden Grafik kann selbst damit experimentiert werden (kann auch in einem vollen Browserfenster dargestellt werden):

Nicht zu reden von Zoom: Das funktioniert in Plotly ausgezeichnet (mit dem Scroll-Rad der Maus in obiger Grafik) und scheint in Matplotlib garnicht implementiert zu sein.

In der Plotly Version zeichnen wir die Visualisierung für alle Frequenzen in eine Grafik. Für die 3D-Ansicht sind alle Visualisierungen bis auf die der gerade selektierten Frequenz versteckt. Über die Frequenz-Legende auf der rechten Seite kann die jeweils anzuzeigende Frequenz ausgewählt werden.

Eine Besonderheit der Plotly Version ist die Anzeige der 3D-Koordinaten auf der Oberfläche der 3D-Visualisierung beim Drüberfahren mit der Maus. In der derzeitigen Implementierung werden dabei Azimuth und Elevations-Winkel sowie absoluter und relativer Gewinn für den gerade ausgewählten Punkt angezeigt.

Allerdings habe ich mit Plotly auch einige Workarounds gebraucht:

• Es gibt einen Bug im Koordinatendisplay

• Plotly unterstützt kein Menü (Legende) bei 3D Oberflächen, Ich musste eine zusätzliche (leere mit einem weissen Pixel) Grafik in die Anzeige aufnehmen um Menü (Legende) anzuzeigen. Das Menü wird verwendet um auf der rechten Seite die Anzeige für eine bestimmte Frequenz auszuwählen.

• In Plotly gibt es keine Möglichkeit, nur einen bestimmten Menüeintrag gleichzeitig anzuzeigen. Ich wollte dass sich die Menü-Einträge wie Radio-Knöpfe (Radio-Buttons) verhalten: Beim Klicken auf eine Frequenz sollte die Visualisierung von nur dieser Frequenz erfolgen. Für 3D Grafiken macht es nicht viel Sinn, mehrere 3D Visualisierungen in einer Anzeige zu haben. Dies wird mit ein paar Zeilen individuellem Javascript erreicht.

Das Antennendiagramm in obigem Beispiel ist aus einem Artikel von L.B. Cebik [2]. Beim Vergleich des Diagramms im Artikel mit obiger 3D-Ansicht fällt auf (auch ein Vergleich der Azimuth/Elevationsdiagramme in der Beschreibung von plot-antenna ist möglich) dass die Plotly 3D Ansicht leicht verzerrt ist (die Keule ist zu lang). Ich habe noch nicht rausgefunden ob das ein Problem von Plotly oder meiner Verwendung davon ist. Ich versuche das gleiche Bildseitenverhältnis für alle Richtungen einzustellen indem ich den Wertebereich (range) für xaxis, yaxis und zaxis einstelle aber es schaut nicht korrekt aus.

Beim Vergleich mit dem Azimuth Diagramm welches mit Plotly erzeugt wurde:

sehen die Größenverhältnisse richtig aus. Auch im 2D-Diagramm wird der Winkel (Azimuth in diesem Fall) beim Überfahren mit der Maus sichtbar. Auch mit Matplotlib ist dies möglich, jedoch rastet der Cursor nicht auf den jeweiligen Punkt auf der Anzeige ein. Damit kann mit Matplotlib nur geraten werden ob sich der Cursor an der korrekten Stelle befindet. Weiters kann wie bei der 3D-Grafik die Frequenz für die Anzeige rechts in der Legende ausgewählt werden. Diesmal ist die gleichzeitige Anzeige von mehreren Frequenzen möglich, was bei einer opaken 3D-Ansicht nicht sehr sinnvoll wäre.

[Änderung 2022-10-18: Ersetze epsilon durch delta in der Beschreibung von Beispiel 7 in pgapack]

Viele Optimierungsprobleme beinhalten Randbedingungen die von gültigen Lösungen erfüllt sein müssen. Ein Optimierungsproblem mit Randbedingungen wird typischerweise als nichtlineares Programmierungs-Problem formuliert [1].

In diesem Problem gibt es \(n\) Variablen (der Vektor \(\vec{x}\) hat die Länge \(n\)), \(J\) Ungleich-Bedingungen, \(K\) Gleichheits-Bedingungen und die Variable \(x_m\) muss im Bereich \(|x_m^l, x_m^u|\) sein (als Box-Bedingung bezeichnet). Die Funktionen \(f_i\) heissen Zielfunktionen. Probleme mit mehreren Zielfunktionen habe ich schon früher in diesem Blog [2] beschrieben. Im folgenden werde ich einige Begriffe ohne weitere Erklärung verwenden die in diesem früheren Blogbeitrag eingeführt wurden.

Die Zielfunktionen werden nicht notwendigerweise minimiert (wie in der Formel angegeben) sondern können auch maximiert werden wenn das Problem dies verlangt. Die Ungleich-Bedingungen werden auch oft mit einem \(\ge\) Zeichen formuliert, die Formal kann einfach umgestellt werden (z.B. durch Multiplikation mit -1) um ein \(\le\) Zeichen zu verwenden.

Nachdem es recht schwer ist, Gleichheits-Bedingungen zu erfüllen, besonders wenn diese in nichtlinearen Funktionen der Eingabevariablen auftreten, werden Gleichheits-Bedingungen oft in Ungleich-Bedingungen umgewandelt unter Verwendung einer δ‑Umgebung:

Wobei δ so gewählt wird dass die Lösung gut genug für das zu lösende Problem ist.

Eine sehr erfolgreiche Methode zur Lösung von Problemen mit Randbedingungen verwendet eine lexicographische Ordnung von Randbedingungen und Zielfunktion(en). Lösungskandidaten des Problems werden zuerst nach verletzten Randbedingungen sortiert (typischerweise der Summe der Verletzung von Randbedingungen) und dann nach dem Wert der Zielfunktion(en) [1]. Beim Vergleich von zwei Individuen während der Selektionsphase des Genetischen Algorithmus gibt es drei Fälle: Wenn beide Individuen die Randbedingungen verletzen gewinnt das Individuum mit der kleineren Verletzung der Bedingungen. Wenn ein Individuum keine Bedingung verletzt, das andere aber schon, gewinnt das Individuum ohne Verletzung von Randbedingungen. Im letzten Fall dass keines der Individuen die Bedingungen verletzt, findet ein normaler Vergleich der Zielfunktionen statt (was vom verwendeten Algorithmus abhängt und ob maximiert oder minimiert wird). Diese Methode, ursprünglich von Deb [1] vorgeschlagen ist in PGAPack und im Python Wrapper PGAPy implementiert.

Mit diesem Algorithmus zur Behandlung von Randbedingungen, werden zuerst Lösungskandidaten gefunden die keine Randbedingung verletzen bevor der Algorithmus die Zielfunktionen überhaupt "anschaut". Daher passiert es dann oft, dass die Suche in einer Region endet wo es keine guten Lösungen gibt (aber dafür keine Randbedingungen verletzt sind). Schwierige Probleme, bei welchen dies auftritt, sind oft Probleme mit Gleichheits-Bedingungen aber es gibt auch andere "schwierige" Randbedingungen. In einem früheren Blog Post [2] zur Antennen-Optimierung habe ich geschrieben: "dass der Optimierungs-Algorithmus Schwierigkeiten hat, überhaupt sinnvolle Lösungen für die Direktor-Variante zu finden. Nur bei einer handvoll Experimente war es überhaupt möglich, die obige Pareto-Front zu finden."

In diesem Experiment habe ich 50 Optimierungsläufe durchgeführt und nur 5 davon sind nicht in einem lokalen Optimum steckengeblieben. Etwas ähnliches passiert bei Problem 7 in Deb's Artikel [1] wo Gleichheits-Bedingungen verwendet werden. Ich habe dieses als Beispiel 7 in PGAPack implementiert. Es wird nur eine Lösung nahe dem (bekannten) Optimum gefunden wenn \(\delta \ge 10^{-2}\) für alle Gleichheits-Bedingungen ist (ich habe nicht mit verschiedenen Anfangswerten für den Zufallszahlengenerator experimentiert, es könnte sein dass man mit einem anderen Startwert bessere Lösungen findet). Im Artikel [1] verwendet Deb \(\delta = 10^{-3}\) aus dem selben Grund.

Eine Methode, um dieses Problem zu verbessern, war attraktiv, weil sie sowohl einfach zu verstehen als auch zu implementieren ist: Takahama und Sakai haben zuerst mit einer Methode experimentiert um unter gelockerten Randbedingungen in der Frühphase der Optimierung zu arbeiten, sie nannten dies einen Genetischen Algorithmus mit α‑Randbedingungen [3]. Die Methode wurde später vereinfacht und als Optimierung mit ε‑Randbedingungen bezeichnet. Sie kann für verschiedene Optimierungsalgorithmen verwendet werden, nicht nur für Genetische Algorithmen und Varianten [4]. Von speziellem Interesse in diesem Zusammenhang ist die Anwendung der Methode für Differential Evolution [5], [6], aber natürlich kann sie auch für andere Formen von Genetischen Algorithmen verwendet werden.

Das ε im Namen dieser Methode kann für das δ bei der Konversion von Gleichheits-Bedingungen in Ungleich-Bedingungen verwendet werden, ist aber nicht auf diesen Anwendungsfall eingeschränkt.

Während des Optimimierungslaufes wird in jeder Generation ein neuer Wert für ε berechnet. Der Vergleich von Individuen wie oben skizziert wird so modifiziert, dass ein Individuum so behandelt wird als verletze es keine Randbedingung wenn die Verletzung der Randbedingungen kleiner ist als ε. Wenn also beide Individuen die Randbedingungen um mehr als ε verletzen gewinnt das Individuum mit der kleineren Verletzung. Wenn ein Individuum die Bedingungen um weniger als ε verletzt, das andere aber mehr, gewinnt das erste. Und schließlich wenn die Verletzung der Randbedingungen bei beiden Individuen kleiner ist als ε findet der normale Zielfunktions-Vergleich statt.

Dieser letzte Fall ist der Schlüssel zum Erfolg des Algorithmus: Obwohl die Suche in eine Richtung zu weniger Verletzung der Randbedingungen fortschreitet, wird gleichzeitig eine gute Lösung bezogen auf die Zielfunktionenen gefunden.

Der Algorithmus startet mit der Initialisierung von \(\varepsilon_0\) mit der Summe der Verletzung der Randbedingungen des Individuums mit dem Index \(\theta\) der nach Verletzung der Randbedingungen sortierten initialen Population. Dabei ist \(\theta\) ein Parameter des Algorithmus zwischen 1 und der Populationsgröße, ein guter Wert ist etwa 20% der Populationsgröße, was auch der Standardwert in PGAPack ist. In jeder Generation \(t\) wird \(\varepsilon_t\) berechnet zu:

bis zur Generation \(T_c\). Nach dieser Generation wird ε auf 0 gesetzt. Der Exponent \(cp\) ist zwischen 2 und 10. Der Artikel von 2010 [6] empfiehlt zur Generation \(T_\lambda = 0.95 T_c\) den Wert \(cp = 0.3 cp + 0.7 cp_\min\) zu setzen, wobei \(cp_\min\) der fixe Wert 3 ist. Der Anfangswert von \(cp\) wird so gewählt dass zur Generation \(T_\lambda\) \(\varepsilon_\lambda=10^{-5}\) ist ausser \(cp\) wäre kleiner, dann wird es auf \(cp_\min\) gesetzt. PGAPack implementiert diese Empfehlung für \(cp\) als Standard-Einstellung, erlaubt aber das Setzen von \(cp\) zu Beginn und während der Laufzeit der Optimierung. Es ist also einfach möglich, andere Änderungen von \(cp\) zu implementieren – die Standard-Einstellung funktioniert jedoch recht gut.

Durch Optimierung mit ε Randbedingungen konnten in meinen Experimenten die Gleichheits-Bedingungen in Beispiel 7 von Deb [1] mit einer Genauigkeit von \(10^{-6}\) approximiert werden, siehe die Variable epsilon_generation im optimizer Beispiel.

Wird der antenna-optimizer mit ε‑generation 50 gestartet (das ist der \(T_c\) aus obigem Algorithmus) bleibt der Algorithmus nur in einem einzigen Fall im lokalen Optimum stecken, alle anderen Fälle finden gute Lösungen:

In dieser Grafik sind alle Lösungen von einem Optimierungslauf die von den Lösungen eines anderen Laufs dominiert werden in Schwarz eingezeichnet. Wir sehen dass die Daten von Lauf 16 keine nicht-dominierte Lösung beigetragen haben (auf der rechten Seite in der Legende fehlt die Nummer 16). In der Grafik kann man die dominierten Lösungen aus der Anzeige nehmen indem man auf den schwarzen Kreis in der Legende klickt.

Wenn der Parameter ε‑generation für den Lauf der im lokalen Optimum geendet hat auf 60 erhöht wird, findet auch der Lauf mit Zufallszahlen-Startwert 16 eine Lösung:

Es ist auch sichtbar dass die Lösungen für alle Experimente ziemlich gut (weil nahe an der Pareto-Front) sind. Der schwarze "Schatten" der dominierten Lösungen ist ziemlich nahe an der wirklichen Pareto-Front und ist gut genug um mit einem einzigen Optimierungslauf ausreichend gute Lösungen zu finden.

Vor kurzem waren in der deutschen Zeitschrift "Funkamateur" zwei Artikel wo der Autor behauptet hat, es sei möglich, ein besseres Stehwellenverhältnis (SWV) mit Hilfe eines 50Ω Kabels zu erreichen.

Der erste Artikel beschrieb eine Duoband-Antenne für die 70cm und 2m Amateurfunkbänder [1]. Der Artikel erwähnte nebenbei dass das SWV durch eine spezifische Länge von 50Ω Kabel verbessert werden könne und enthielt einen Link zu einer Software des selben Autors um diese Länge zu berechnen (es wurde behauptet dass verschiedene Längen von Kabel für den zu erreichenden Zweck möglich seien). Die Behauptung war, es fände eine Leitungstransformation durch das 50Ω Kabel statt.

Im zweiten Artikel – ausgelöst von mehreren Lesern die Leserbriefe geschrieben hatten – lieferte der Autor Messergebnisse die klar zeigten dass, in der Tat, das SWV durch eine bestimmte Kabellänge verbessert wurde [2].

Nun sagt die Leitungstheorie dass eine Impedanztransformation mit Änderung des Stehwellenverhältnisses nur mit einer Leitung einer anderen Impedanz als der Systemimpedanz erreicht werden kann, siehe auch meinen kürzlich erschienenen Artikel zu diesem Thema. Es muss also einen anderen Grund geben, warum sich die Impedanz ändert – wie klar von den Messungen des Autors gezeigt wurde. Der Autor meint in seinem Artikel dass es wohl eine Diskrepanz zwischen Theorie und Praxis gäbe – ich zeige im folgenden dass es eine ganz gute Erklärung für die Beobachtungen gibt ohne die existierende elektromagnetische Theorie über Bord zu werfen.

Der Hinweis hier ist, dass der Autor keinen Balun oder einen anderen Mechanismus zur Unterdrückung von Mantelwellen auf der Zuleitung verwendet hat. Also habe ich mir den Spaß gemacht, die Verhältnisse mit dem Antennen-Simulationsprogramm NEC2 zu simulieren. Die Resultate sollten mit einem anderen Modellierungsprogramm reproduzierbar sein, zumindest wenn dieses einen NEC Kern (z.B. EZNEC) enthält der auch ein Bodenmodell enthält und nicht nur im Freiraum simulieren kann. Ich habe die (einfache Zwei-Elemente) Antenne aus dem Original-Artikel modelliert.

Die folgende Grafik zeigt das Stehwellenverhältnis (SWV) und den Real- und Imaginäranteil der Impedanz direkt an der Antenne. Das erste NEC input file kann verwendet werden um die Resultate zu reproduzieren. Ich habe das Sommerfeld/Norton Bodenmodell mit einer mittleren Bodenqualität verwendet. Die Antenne ist 15m über dem Boden in diesem Modell weil sie zur Befestigung auf einem Balkon vorgesehen ist.

Das nächste Bild zeigt die selbe Antenne mit einem Stück Speiseleitung (die Länge ist die von der der Autor behauptet sie würde das SWV verbessern). Die Speiseleitung ist nicht als Draht modelliert, sondern als Übertragungsleitung (für die bietet NEC eine eigene Abstraktion). Wieder zeigen wir Real- und Imaginärteil der Impedanz. Die elektrische Länge (NEC verwendet bei Übertragungsleitungen die elektrische Länge ohne Verkürzungsfaktor) ist etwa 1.81m berechnet mit der Formel aus [1] mit einem Verkürzungsfaktor von 1. Es ist klar ersichtlich, dass sich zwar die Impedanz ändert aber das SWV gleich bleibt. Das korrespondierende zweite NEC input file kann auch hier verwendet werden um die Resultate zu reproduzieren.

Um nun die Verhältnisse am Speisepunkt zu simulieren wenn es keinen Balun gibt habe ich einen Draht zwischen der Antenne und dem Speisepunkt am Ende des Kabels im dritten NEC input file modelliert. Das folgende Bild zeigt das Modell, wir haben einen zusätzlichen Draht der die Aussenseite der Coax-Speiseleitung darstellt. Dieser Draht ist verbunden mit der Antenne auf der einen Seite und dem Ende des Coaxkabels auf der anderen Seite. Die Länge des Kabels ist die physikalische Länge. Ausser einer geringfügigen Erhöhung der elektrischen Länge wegen der Isolation (nicht mit dem üblichen Coax Verkürzungsfaktor, wir simulieren ja nur das äussere Drahtgeflecht!) muss der Draht im Modell die echte physische Länge des Coax wie mit der Formel aus dem ersten Artikel berechnet [1] verwenden. Wir verwenden 1.484m aus [1], Tabelle 3 für ein Aircell5 mit dem Faktor n=1 aus der Tabelle im Artikel. Wir sehen dass sich in der Tat, wie im zweiten Artikel [2] beobachtet, das Stehwellenverhältnis SWV verkleinert hat. Der Grund ist dass, wegen des fehlenden Balun, sich die Antenne geändert hat: Der Aussenmantel der Coax-Speiseleitung ist Teil der Antenne.

Wenn man sich die Antenne mit dem Linux-Programm xnec2c ansieht, bemerkt man die Ströme auf der Aussenseite der Speiseleitung. Wir können klar sehen, dass da ein Strom auf der Speiseleitung fließt (die Speiseleitung in dem Bild ist der lange Teil rechts, parallel dazu liegt die modellierte Übertragungsleitung zur Speisung).

Also ist die elektromagnetische Theorie gerettet, es gibt eine physikalische Erklärung des Phänomens. Trotzdem wäre es noch interessant die verschiedenen Einflüsse von Höhe über Grund oder unterschiedliche Bodencharacteristika zu modellieren – um herauszufinden wie reproduzierbar die Ergebnisse mit unterschiedlichen Parametern sind , weil der Autor ja behauptet das beobachtete Verhalten sei "seit Jahren bekannt" ([2] S. 35).

Seit einiger Zeit optimiere ich Antennen mit Genetischen Algorithmen. Ich verwende dafür das Paket für parallele Genetische Algorithmen, pgapack das ursprünglich von David Levine vom Argonne National Laboratory entwickelt wurde und das ich pflege (wie übersetze ich "maintain"?). Noch länger als die Pflege von pgapack entwickle ich eine Python Bibliothek für pgapack namens pgapy.

Für die Antennensimulation verwende ich Tim Molteno's PyNEC Paket, eine Bibliothek in Python für ein Paket zur elektromagnetischen Simulation namens "Numerical Electromagnetics Code (NEC) Version 2", geschrieben in C++ (auch bekannt als NEC++).

Mit diesen Paketen habe ich ein kleines Open Source Framework für die Antennen-Optimierung namens antenna-optimizer geschrieben. Dieses kann traditionelle Genetische Algorithmen mit Bit-Strings als Genen als auch Fließkomma-Repräsentationen mit dafür geeigneten Operatoren verwenden.

Das "Parallel" in pgapack sagt uns, dass die Evaluierungsfunktion des Genetischen Algorithmus parallelisiert werden kann. Bei der Optimierung von Antennen simulieren wir für jede Paramter-Kombination die einen Antenne darstellt diese mit PyNEC. Antennensimulation ist nach wie vor (der ursprüngliche NEC Code ist aus den 80er Jahren und wurde im Zeitalter von Lochkarten entwickelt) ein CPU-instensives Unterfangen. So ist es eine gute Nachricht, dass wir mit pgapack viele Simulationen parallel mit dem Message Passing Interface (MPI) Standard [1] durchführen können.

Für pgapack – und damit auch für pgapy – habe ich in letzter Zeit einige bewährte klassische Algorithmen implementiert:

• Differential Evolution [2], [3], [4] ist ein erfolgreicher Optimierungs-Algorithmus für Fließkomma-Gene der für elektromagnetische Aufgaben sehr interessant ist.

• Der "elitist Nondominated Sorting Genetic Algorithm" NSGA-II [5] erlaubt es in einem einzigen Optimierungslauf mehrere Zielfunktionen zu optimieren

• Wir können Wertebeschränkungen im Zielbereich einer Funktion definieren die minimiert werden, sogenannte "Constraints" [6]. Damit eine Lösung gültig ist, müssen alle Constraints 0 oder negativ sein.

Traditionell ist mit Genetischen Algorithmen nur eine Evaluierungsfunktion, die sogenannte Zielfunktion möglich. Mit NSGA-II können mehrere Zielfunktionen oder Kriterien verwendet werden. Das Englische nennt solche Algorithmen "Multi-Objective Optimization" (Optimierung mehrerer Ziele).

Für die Antennensimulation heisst dass, wir müssen nicht mehrere Kriterien für eine gute Antenne wie Gewinn, Vorwärts/Rückwärtsverhältnis und Stehwellenverhältnis (SWV) in einer einzigen Evaluierungsfunktion zusammenfassen, was ich bisher im antenna-optimizer gemacht habe, sondern können sie separat definieren und dem Genetischen Algorithmus die Optimierung überlassen.

Mit mehreren Zielfunktionen ist es allerdings typischerweise so, dass die Verbesserung bei einem Ziel eine Verschlechterung bei einem anderen Ziel zur Folge hat und umgekehrt. Wir suchen also Lösungen die strikt besser sind als andere Lösungen. Eine Lösung dominiert eine andere Lösung wenn sie in einem Kriterium besser aber in keinem der anderen Kriterien schlechter ist als die andere Lösung. Alle Lösungen die dieser Definition entsprechen heissen Pareto-Optimal nach dem Italienischen Wissenschaftler Vilfredo Pareto der das Konzept von Pareto-Optimalität zuerst definiert hat. Alle Lösungen die dieses Optimalitäts-Kriterium erfüllen, liegen auf einer sogenannten Pareto-Front. Wenn wir nur zwei Zielfunktionen haben, können wir die Pareto-Front gut mit einem Scatterplot (Streudiagramm) darstellen, wie wir weiter unten sehen werden.

Weil pgapack einen Ansatz verfolgt, der die freie Kombination von Algorithmen unterstützt, können wir verschiedene erfolgreiche Strategien für unterschiedliche Teile des Genetischen Algorithmus kombinieren:

• Für den Mutation/Crossover Teil können wir Differential Evolution (DE) verwenden

• DE wiederum kann für mehrere Zielfunktionen mit der (Generations-) Ersetzungsstrategie von NSGA-II kombiniert werden

• Wir können einige der Zielfunktionen als Beschränkungen (Constraints) definieren. Für unser Problem ist es sinnvoll nur Antennen zu erlauben die ein bestimmtes Stehwellenverhältnis (SWV) nicht überschreiten. Wir erlauben also keine Antenne mit einem SWV > 1.8. Die zugehörige Zielfunktion ist \(S - 1.8 \le 0\) wobei \(S\) das SWV ist.

Mit dieser Kombination können wir erfolgreich Antennen für das 70cm Amaterufunkband (430 MHz - 440 MHz) berechnen. Die Antenne verwendet einen sogenannten Faltdipol (das Ding mit den Rundungen in der Grafik) und ein gerades Element. Die Messlinien in der Grafik repräsentieren die vom Genetischen Algorithmus optimierten Längen. Die zwei Punkte in der Mitte des Faltdipols zeigen den Punkt wo die Speiseleitung angeschlossen wird.

Ein erstes Beispiel simuliert Antennenparameter für die höchste, niedrigste und mittlere Frequenz. Der Gewinn und das Vorwärts/Rückwärtsverhältnis (Forward/Backward ratio in der Grafik) werden nur für die mittlere Frequenz berechnet:

In dieser Grafik (ein Scatterplot) wird die erste Zielfunktion (der Gewinn) gegen die zweite Zielfunktion, das Vorwärts/Rückwärtsverhältnis aufgetragen. Alle Zahlen sind für die mittlere Frequenz. Jeder Punkt in der Grafik repräsentiert eine Antenne. Alle Antennen haben ein SWV kleiner als 1.8 auf der kleinsten, mittleren und höchsten Frequenz.

Nach diesem Erfolg habe ich dann angefangen, mit unterschiedlichen Einstellungen für die Parameter des Differential Evolution (DE) Algorithmus zu experimentieren. Es ist in der Literatur zu DE bekannt dass für "zerlegbare" (decomposable) Probleme eine niedrige Crossover-Rate besser ist, für nicht-zerlegbare eine höhere. Ein zerlegbares Problem zeichnet sich dadurch aus, dass die unterschiedlichen Dimensionen getrennt voneinander optimierbar sind, dies wurde zuerst von Salomon 1996 [7] beobachtet. Ursprünglich hatte ich eine Crossover-Rate von 0.2 eingestellt und meine Hoffnung war, dass die Optimierung mit einer höheren Rate besser und schneller funktionieren würde. Das Experiment unten verwendet eine Crossover-Rate von 0.9.

Zusätzlich habe ich mit "Dither" für den Skalierungsfaktor \(F\) von Differential Evolution (könnte vielleicht mit einem leichten "Zittern" Skalierung übersetzen) experimentiert. Mit diesem Faktor wird die Differenz von zwei Vektoren der Input-Parameter bei der Anwendung von DE multipliziert. In der ersten Implementierung hatte ich Dither auf 0 gesetzt, jetzt hatte ich 0.2 eingestellt. Ich war dann sehr überrascht dass ich mit diesen Einstellungen eine neue Pareto-Front als Lösung gefunden habe:

Damit man besser sieht, dass die zweite Front komplett die zuerst gefundene Front dominiert habe ich beide in einer Grafik geplottet:

Weil diese zweite Front (über den gesamten Frequenzbereich) für eine Zwei-Elemente Antenne zu gut ausschaut um wahr zu sein schauen wir uns das im Folgenden genauer an. Zuerst schauen wir auf die Orientierung der Antenne und des berechneten Gewinn-Diagramms für eine der Antennen in der Mitte der unteren Front:

Die Antenne hat – wie schon aus der ersten Pareto-Front Grafik zu entnehmen – einen Gewinn von etwa 6.6 dBi und ein Vorwärts/Rückwärtsverhältnis von etwa 11 dB in der Mitte des Bandes bei 435 MHz. Die Einfärbung der Antenne deutet die Ströme in der Antennenstruktur an. Wer sich das selbst anschauen möchte: Hier ist ein Link zur NEC Input Datei für Antenne 1.

Vergleichen wir diese Antenne mit einer aus der Mitte der "orangen Front" wo wir deutlich bessere Ergebnisse bekommen:

Diese Antenne ist aus der Mitte der oberen Pareto-Front und hat einen Gewinn von etwa 6.7 dBi sowie ein Vorwärts/Rückwärtsverhältnis von etwa 16 dB in der Mitte des Bandes bei 435 MHz. Wer entdeckt den Unterschied zur vorherigen Antenne? Ja: der maximale Gewinn liegt in der Gegenrichtung der ersten Antenne. Wir sagen dass bei der ersten Antenne das gerade Element ein Reflektor ist, während bei der zweiten Antenne das gerade Element als Direktor arbeitet. Auch hier ist ein Link zur NEC Input Datei für Antenne 2.

Nun schauen wir uns über die ganze Frequenz den Gewinn und das Vorwärts/Rückwärtsverhältnis an. Die linke Grafik ist für die erste Antenne (die mit dem Reflektor-Element), die rechte für die zweite Antenne wo das gerade Element als Direktor arbeitet.

Wir sehen dass sich das Vorwärts/Rückwärtsverhältnis bei der Direktor Antenne von etwas mehr als 10 dB bis zu über 25 dB bewegt, während es für das Reflektor-Design von 9.3 dB bis 11.75 dB geht. Der Minimalgewinn ist beim Reflektor-Design etwas besser (6.35-6.85 dBi bzw. 6.3-7.05 dBi). Also brauchen wir weitere Experimente. Wenn man die Lösungen auf ein Reflektor-Design beschränkt und fordert dass der Minimalgewinn und das minimale Vorwärts/Rückwärtsverhältnis an den drei Stellen (untere, mittlere, obere Frequenz) genommen wird bekommen wir:

Für das gleiche bei einem Direktor Design (auch mit minimalem Gewinn und Vorwärts/Rückwärtsverhältnis bei allen drei Frequenzen (untere, mittlere, obere Frequenz) bekommen wir:

Mit diesen Resultaten ist der Sweet Spot für eine Antenne die man wirklich bauen will wahrscheinlich bei etwa 10 dB Vorwärts/Rückwärtsverhältnis oder darüber und bei einem Gewinn von etwa 6.2 dBi. Ein paar zehntel-dB mehr Gewinn rauszuholen und dabei mehrere dB Vorwärts/Rückwärtsverhältnis zu verlieren scheint nicht sehr sinnvoll. Beim Vergleich des Direktor mit dem Reflektor-Design fällt auf (was zumindest meiner Intuition widerspricht) dass das Direktor-Design ein besseres Vorwärts/Rückwärtsverhältnis über den gesamten Frequenzbereich hat. Wenn allerdings die Antenne für Relais-Kommunikation Verwendung finden soll wo die Sendefrequenz (die Relais-Eingabe) im unteren Teil des Frequenzbandes liegt und die Relais-Ausgabe (unsere Empfangsfrequenz) in der oberen Hälfte, würden wir vermutlich ein Reflektor-Design realisieren weil der Gewinn beim Senden besser ist und das Vorwärts/Rückwärtsverhältnis beim Empfang besser ist (vergleiche die beiden vorherigen Grafiken zu Gewinn und Vorwärts/Rückwärtsverhältnis).

Man beachte auch, dass der Optimierungs-Algorithmus Schwierigkeiten hat, überhaupt sinnvolle Lösungen für die Direktor-Variante zu finden. Nur bei einer handvoll Experimente war es überhaupt möglich, die obige Pareto-Front zu finden. Das Direktor-Design ist schmalbandiger als das Reflektor-Design und da das Stehwellenverhältnis ein Beschränkung (Constraint) ist, werden oft nur lokale Optima gefunden. Der größere Unterschied im Wertebereich von Gewinn- und Vorwärts/Rückwärtsverhältnis für das Direktor Design sagt uns ausserdem dass es schwieriger zu realisieren sein wird: Wenn die Dimensionen nicht exakt richtig getroffen werden wird die Antenne wohl weiter von den vorhergesagten Simulationsergebnissen abweichen. Das Reflektor-Design ist etwas toleranter in dieser Beziehung.

Als Funkamateur ist man hin und wieder mit der Herausforderung konfrontiert, eine Antenne an einen Transceiver über eine Übertragungsleitung anzuschließen. Oft ist dabei die Impedanz der Antenne unterschiedlich von der Impedanz der Übertragungsleitung und des Transceivers. Die Impedanz von Transceiver und einer typischen Coax Übertragungsleitung ist üblicherweise 50Ω im Amateurfunkbereich. Die Impendanz der Antenne ist je nach technischen und Gegebenheiten und Umgebungseinflüssen unterschiedlich.

Eine Übertragungsleitung transformiert je nach Länge die Impedanz der Antenne. Eine gut bekannte Formel wird z.B. von Chipman [1] S.134 Formel 7.15 angegeben:

In dieser Formel ist \(Z_d\) die Impedanz in einem Abstand von der Last \(d\) zur Antenne, \(Z_l\) ist die Impedanz an der Last (also der Antenne) und \(Z_0\) ist die characteristische Impedanz der Übertragungsleitung und des Transceivers, typischerweise 50Ω. \(\gamma\) ist der komplexe Übertragungskoeffizient, diesen kann man in den Real- und den Imaginärteil aufteilen, wobei dann \(\alpha\) die Dämpfung in Neper/m und \(\beta\) die Phasenkonstante ist, \(j\) ist die imaginäre Einheit (oft auch als \(i\) geschrieben aber in der Elektrotechnik meist als \(j\)):

Es wird manchmal behauptet, dass ein Kabel das Stehwellenverhältnis am Transceiver verbessern kann, weil es die Impedanz transformiert. Dies stimmt aber nur dann, wenn die characteristische Impedanz des Kabels unterschiedlich von der Impedanz des Transceivers ist (wenn wir annehmen dass das Kabel verlustfrei ist was in vielen Fällen von kurzen Kabeln gilt) was ich im folgenden zeigen werde.

Im verlustfreien Fall ist \(\alpha\) in der obigen Formel 0. Wir können \(\beta\) über die Frequenz und schließlich über die Wellenlänge \(\lambda\) ausdrücken:

und:

Die Frequenz \(f\) ist in Hz, \(c\) ist die Lichtgeschwindigkeit und \(\mathbb{VF}\) der Verkürzungsfaktor (im Englischen velocity factor) der Übertragungsleitung. Wenn wir die Entfernung von der Last \(d\) in Chipman's Formal als Vielfaches von \(\lambda\), \(l_\lambda\) ausdrücken, kürzen sich \(f\), \(c\) und \(\mathbb{VF}\), und wir bekommen für den verlustfreien Fall:

Der komplexe Reflexionskoeffizient \(\rho\) ist (siehe z.B. bei Chipman [1] 7.9 S.128):

wobei \(Z\) die Impedanz am Punkt der Leitung ist wo wir den Reflexionskoeffizienten wissen wollen. Daraus können wir das Stehwellenverhältnis berechnen (vgl. z.B. Chipman [1] 8.21 S.165) das ich hier mit \(S\) bezeichne:

Ein einfaches Python Programm (Python ist angenehm für solche Berechnungen weil es komplexe Zahlen unterstützt und freie Software ist) um obiges zu berechnen könnte wie folgt aussehen:

from math import e, pi
def impedance (z_load, l_lambda, z0 = 50.0):
    zl = z_load / z0
    ex = 2j * pi * l_lambda
    lp = e **  (ex) * (zl + 1)
    lm = e ** (-ex) * (zl - 1)
    return z0 * (lp + lm) / (lp - lm)

def vswr (z, z0 = 50.0):
    absrho = abs ((z - z0) / (z + z0))
    return (1 + absrho) / (1 - absrho)

Wenn wir das nun für einige Punkte berechnen, z.B. für eine Antenne mit 72Ω und für einige Vielfache von \(\lambda\) bekommen wir immer den selben Wert für das Stehwellenverhältnis:

Zu beweisen dass das Stehwellenverhältnis immer gleich ist, egal wie lang die Übertragungsleitung ist, wird an dieser Stelle als Hausaufgabe für den Leser gelassen. Ein Hinweise: Es genügt zu zeigen dass der Absolutbetrag von \(\rho\) gleich bleibt.

Der interessantere Fall ist aber nun, wenn wir Kabelverluste mit einbeziehen. Ich habe ein Stück Software geschrieben, das das Verhalten einer Coax Leitung aus den Herstellerangaben berechnen kann, eine Idee die schon vor langer Zeit von Frank Witt, AI1H [2] publiziert wurde. Die Implementierung ist Teil meines Open Source antenna-optimizer Projekts. und beinhaltet ein kleines Programm namens coaxmodel. Mit diesem kann die Impedanz (und Stehwellenverhältnis) für ein echtes Kabel berechnet werden. Die Implementierung enthält schon die Modelle einiger Kabel aber es ist recht einfach, neue dazuzugeben (siehe am Ende von coaxmodel.py). Zusätzlich kann man damit auch Anpassungen mit Hilfe einer Stichleitung berechnen: Wenn man ein Stück Übertragungsleitung parallel zur Speiseleitung in einem bestimmten Abstand zur Last anschließt (dieses Stück heisst dann Stichleitung) transformiert sich die Impedanz auf eine Weise dass der Generator (der Transceiver) eine Last von 50Ω sieht. Für obiges Beispiel würde es z.B. berechnen (Nur ein Teil der Ausgabe ist hier wiedergegeben, probieren Sie es selbst):

% coaxmodel -z 72 -f 435e6 -l .057 match
0.06 m at 435.00 MHz with 100 W applied
           Load impedance 72.000 +0.000j Ω
          Input impedance 46.857 -17.433j Ω
             VSWR at load 1.440
            VSWR at input 1.439
Inductive stub with open circuit at end:
            Stub attached 0.06246 m from load
              Stub length 0.20224 m
      Resulting impedance 50.00 -0.00j

Das verrät uns dass ein 5.7cm Stück Übertragungsleitung die 72Ω Impedanz an der Last (Antenne) in eine 46.86-17.43jΩ Impedanz am Eingang der Übertragungsleitung transformiert. Es wird auch klar dass sich das Stehwellenverhältnis durch Kabelverluste ganz leicht verbessert hat und dass man diesen Unterschied bei dieser Länge des Kabels vernachlässigen kann.

Wenn man nun eine 20.2cm Stichleitung von 50Ω mit einem offenen Ende parallel zur Speiseleitung im Abstand von 6.2cm von der Last anschließt, wird die Impedanz auf 50Ω transformiert mit einem Stehwellenverhältnis von 1:1.

Der Orange-Pi Zero ist ein beliebter Open Hardware Einplatinencomputer der unter Linux laufen kann. Er ist gut von neuen Linux Kerneln und dem U-Boot Bootloader unterstützt.

Neuere boards haben ein NOR-Flash das verwendet werden kann um den U-Boot Bootloader aufzuspielen. U-Boot kann aus den Sourcen mit folgenden Kommandos gebaut werden:

make orangepi_zero_defconfig
make

Dann kann mit den sunxi-tools (Debian Linux hat diese Tools als Paket seit etwa zwei Releases) der Bootloader ins NOR-Flash geladen werden. Zuerst wird der Orange-Pi USB-OTG Port (das ist der kleine µ-USB auf dem Board) mit einem Micro-USB Kabel mit einem Computer verbunden (am besten ein Debian Linux, mit anderen Systemen müssen Sie selbst rausfinden wie das geht). Der Orange-Pi darf keine SD-Karte eingesteckt haben. In meinen Experimenten war die Stromversorgung aus dem USB genug um mit folgendem Kommando den Bootloader ins NOR-Flash zu bekommen:

sunxi-fel -v -p spiflash-write 0 u-boot-sunxi-with-spl.bin

Die Datei u-boot-sunxi-with-spl.bin haben wir gerade mit obigen Kommandos gebaut. Bis hierher ist alles auf der sunxi.org Orange-Pi Zero Seite dokumentiert. Aber sobald der Bootloader mal geflasht ist, bootet der Orange-Pi nicht mehr in den FEL-Modus der für obiges Kommando gebraucht wird. Stattdessen startet der Bootloader aus dem NOR-Flash. Wenn jemals wieder eine neuere Version von U-Boot geflasht werden soll, muss der Orange-Pi wieder im FEL-Modus starten. Einige Methoden dazu sind auch auf der sunxi.org Seite dokumentiert. Leider ist die einfachste, dass man den RECOVERY-Kontakt auf Masse legt, was bei vielen Boards über einen Jumper möglich ist, beim Orange-Pi nicht so einfach weil der Anschluss über einen Pull-Up Widerstand (R 123) ohne Jumper fest angeschlossen ist. Zum Glück ist dieser Widerstand auf dem Board einfach zugänglich. Ich habe im Bild die Seite von R 123 auf dem Board mit dem rechten roten Pfeil markiert.

Diese Seite des Widerstands muss mit Masse verbunden werden. Ich verwende die Masse von der Debug-Schnittstelle (markiert mit dem linken roten Pfeil) weil alle anderen gut zugänglichen Masse-Anschlüsse direkt daneben einen +5V Kontakt haben. Wenn man versehentlich +5V mit irgendwas anderem auf dem Board kurzschließt kann man das Board zerstören. Daher ist der Masseanschluss an der Debug-Schnittstelle die sicherste Lösung. Wir verbinden Masse und den R 123 Kontakt, beide mit roten Pfeilen markiert, und schalten den Strom ein (z.B indem wir den OTG-USB verbinden). Dann ist zu überprüfen ob das Board wirklich im FEL-Modus ist mit dem Kommando:

sunxi-fel ver

Das Kommando druckt Versionsinformation zum Board oder eine Fehlermeldung wenn der FEL-Modus nicht erfolgreich gestartet wurde. Sobald das Board dann mal im FEL-Modus ist, kann ein neuer Bootloader aufgespielt werden.

Erfolgsmeldung vom obigen Kommando:

AWUSBFEX soc=00001680(H3) 00000001 ver=0001 44 08 scratchpad=00007e00 00000000 00000000

Fehlermeldung wenn der FEL-Modus nicht gestartet wurde:

ERROR: Allwinner USB FEL device not found!

[Edit 2021-10-25 Tippfehler gefixt, Danke Peter!]

Als ich vor kurzem Übertragungsleitungen für mein antenna-optimizer Projekt studiert habe, bin ich über eine Formel für den Q-Faktor eines Coax-Resonators von Frank Witt [1] gestolpert:

In dieser Formel ist \(F_0\) die Frequenz des Resonators in MHz, \(A\) is der Verlust in dB pro 100 ft, und VF ist der velocity factor (oder auf Deutsch "Verkürzungsfaktor") des Kabels. Die Formel war für einen \(\frac{\lambda}{4}\) Resonator angegeben. Ich habe mich über die Formel gewundert, weil ich dachte dass ein logarithmischer Wert für den Verlust eine Exponentiation bei der Berechnung des Q-Faktors erfordern sollte.

Wenn wir bei Wikipedia die Definition der Güte eines Schwingkreises nachschlagen bekommen wir:

Wobei \(E_s\) die gespeicherte Energie ist und \(E_d\) die während der Schwingung als Wärme verbrauchte Energie.

Wir wissen, dass der Verlustfaktor eines Kabels in dB Leistungsverluste definiert. Wenn der Verlust in dB pro 100m (wir verwenden natürlich metrische Einheiten) \(a\) ist, haben wir für den Verlust in dB:

Wobei \(l\) die Länge in Meter ist. Für einen \(\frac{\lambda}{4}\) Resonator bekommen wir:

Um die Verlustleistung als Faktor zu bekommen (statt in dB) bekommen wir:

Wir setzen

in die Formel ein, wobei \(c\) die Lichtgeschwindigkeit ist und \(f_0\) die Resonanzfrequenz in Hz:

Zurück zur Wikipedia-Formel die ja Energien benötigt müssten wir hier integrieren – aber weil beim Integrieren eines Verhältnisses wieder das selbe Verhältnis rauskommt machen wir so eine Art vereinfachende Integration wo wir uns nur überlegen müssen wie lang die von der Energie zurückgelegte Distanz ist. Die Wikipedia-Definition behandelt eine volle Periode, also \(\lambda\), nicht \(\frac{\lambda}{4}\). Und der Q-Faktor ist das Verhältnis von gespeicherter zu verlorener Energie. Wir haben also:

Wie wir sehen ist das \(Q\) des Resonators unabhängig vom Resonator-Typ, sei es ein \(\frac{\lambda}{4}\) oder \(\frac{\lambda}{2}\) Resonator.

Wenn wir Q-Faktoren für einige Kurzwellen-Frequenzen (in der Grafik sind die Frequenzen in MHz) gegen Verluste in dB auftragen, sehen wir dass obige Formel recht nahe bei der Näherungsformel von Witt [1] ist.

Plotten wir den Q-Faktor gegen die Frequenz (auch im Kurzwellen-Bereich) für typische Kabel, sehen wir auch dass der Fehler nicht sehr hoch ist, jedenfalls bei höheren Frequenzen.

Den relativen Fehler können wir auch darstellen, wieder für typische Kabel über den gesamten Kurzwellenbereich. Wir sehen, dass der Fehler recht hoch ist für niedrige Frequenzen und verlustreiche Kabel wie RG174.

Wir sehen also dass Witts Formel vermutlich eine Approximation ist, die recht gut für hohe Frequenzen und niedrige Verluste funktioniert. Die Frage ist aber offen: Woher kommt die Formel und woher kommt die magische Konstante die vermutlich alle physischen Konstanten in der Formel zusammenfasst?

Beim Studium von Übertragungsleitungen bin ich auch über ein älteres Buch gestolpert, das mal ein Universitäts-Lehrbuch war [2]. Auf S. 222 leitet Chipman eine Näherungsformel für \(Q\) her, die als vereinfachende Annahme ein niedriges \(\alpha\) (den Dämpfungskoeffizienten in Neper) voraussetzt. Chipman's Näherungsformel ist:

Hier ist \(\alpha\) der Dämpfungskoeffizient in Neper/m und \(\beta\) ist der Phasenfaktor in Radians pro Meter. Das Subscript \(r\) steht für Resonanz. Wir schreiben:

und

und Witts Verlust \(A\) pro 100 ft können wir auch schreiben (wir konvertieren Neper zu dB) als:

Wobei die Konstante 3.2808 der Konversionsfaktor m/ft ist. Lösen wir für \(\alpha_r\) und ersetzen \(\lambda\) und \(\alpha_r\) in Chipman's Formel bekommen wir:

Wobei das \(F_0\) die Frequenz \(f_0\) in MHz ist.

Wir sehen dass die Annahme von niedrigem \(\alpha\) in Chipman's (und Witt's) Formel für höhere Frequenzen und niedrige Verluste zutrifft: \(\alpha\) in Neper pro Meter ist ein konstanter Faktor wenn wir es in dB (pro 100m oder pro 100 ft) umrechnen. Dass die Näherung mit hohen Frequenzen besser wird, liegt daran, dass Kabelverluste typischerweise mit der Wurzel aus \(\lambda\) steigen, während \(\lambda\) invers proportional mit der Frequenz sinkt. Also hat z.B. ein \(\frac{\lambda}{4}\) Resonator höhere Verluste bei niedrigen Frequenzen.