Modellieren einer Drahtantenne mit Isolation

Ralf Schlatterbeck 2024-09-17 15:45

[Update 2024-08-20: Korrektur der Formel zu L und äquivalentem Radius]

[Update 2024-09-19: Korrektur Literaturangabe Artikel von Stearns, Ergänzung zum Antenna Book [5]]

Für mein Antennen-Modellierungsprogramm pymininec (das eine Reimplementierung des original Mininec Basic-Codes darstellt) habe ich mich damit beschäftigt, wie man isolierte Drähte in Antennen modellieren kann. Ich hatte Roy Lewallen, W7EL, den Autor von EZNEC gefragt was er in EZNEC verwendet und einen Hinweis auf ein Paper von J. H. Richmond [1] erhalten, das das erste Paper war, welches einen Algorithmus für isolierte Drähte vorgeschlagen und auch implementiert hat. Der Algorithmus wurde von Jerry McCormack in seiner Masterarbeit [2] implementiert und später in einem Report [3] (dieser ist fast identisch mit der Masterarbeit bis auf neueren Fortran Code und einen anderen Autor) wiederveröffentlicht. Der original Fortran Code ist nicht nur in den obigen Reporten abgedruckt, sondern auch heute noch verfügbar auf Ray L. Cross' Seite zum "Antenna Scatterers Analysis Program" ASAP. Die Subroutine "DSHELL" implementiert den Algorithmus für isolierte Drähte. Mit einigen Modifikationen konnte ich diesen Code auch heute noch ausführen. Großer Dank gebührt an dieser Stelle Roy, W7EL der mir diese Infos zugänglich machte.

Um meine Implementierung in pymininec gegen gemessene Daten zu testen habe ich ein Paper von David Lamensdorf [4] ausgegraben wo der Effekt von Isolation wirklich mal gemessen wurde. Weil zu dieser Zeit der Einfluss von Isolation nur mit einem unendlichen Zylinder theoretisch modelliert werden konnte verwendete das Experiment Isolation die über das Ende des Drahtes hinausragte bis keine Änderung der gemessenen Parameter mehr beobachtet wurde.

Experimente mit den Formeln von Richmond [1] zeigten mir dass die resultierende Impedanzmatrix schlecht konditioniert ist: Für kurze Drahtsegmente erzeugen die Formeln sehr große absolute Werte die zur Hauptdiagonale addiert und von den beiden Nebendiagonalen (die Matrix ist symmetrisch, die beiden Nebendiagonalen sind für die obere und untere Dreiecksmatrix ident) subtrahiert werden. Weil diese Werte während der Matrix-Inversion subtrahiert werden, kommt es zu katastrophalen Auslöschungseffekten, was die Ergebnisse dominiert.

Also habe ich weiter nachgeforscht und bin auf ein Paper [5] und eine Präsentation [6] von Steve Stearns, K6OIK gestoßen, die Formeln für die verteilte Induktivität und einen äquivalenten Radius angibt mit welchen sich Isolation modellieren lässt. Seine Präsentation hat auch eine gute Übersicht über unterschiedliche Ansätze zur Modellierung von Isolation mit Programmen die dies nicht unterstützen (wie die original NEC-2 Implementierung oder das original Mininec):

\begin{align*} a_e &= a \left(\frac{b}{a}\right)^{\left(1- \frac{1}{\varepsilon_r}\right)} \\ L &= \frac{\mu_0}{2\pi}\left(1-\frac{1}{\varepsilon_r} \right)\log\left(\frac{b}{a}\right) \\ \sigma_e &= \sigma\left(\frac{a}{b}\right)^{2\left(1-\frac{1}{\varepsilon_r}\right)} \\ \end{align*}

Hier ist $a$ der Original-Radius des Drahtes, $b$ der Radius mit Isolation, $\sigma$ ist die Elektrische Leitfähigkeit des Drahtes, $\varepsilon_r$ ist die relative Permittivität der Isolation, $a_e$ ist der äquivalente Radius und $\sigma_e$ die äquivalente Permittivität die verwendet wird um den Effekt des veränderten Radius bei der Berechnung des Skin-Effekts aufzuheben. Die Induktivität $L$ ist die Induktivität pro Länge des isolierten Drahtes (oder Drahtsegments).

Dann habe ich ein Paper gesucht das die Formeln von Steve Stearns bestätigt und habe ein noch älteres Paper von Tai Tsun Wu [7] gefunden das exakt die gleichen Formeln für den äquivalenten Radius $a_e$ und für die verteilte Impedanz $L$ angibt.

Die Formel für die äquivalente Permittivität wird nur verwendet um den Effekt des äquivalenten Radius auf die Berechnung des Skin Effekts aufzuheben. In einem früheren Blog Post [8] hatte ich eine Formel für den Skin Effekt angegeben. In dieser Formel taucht die Elektrische Leitfähigkeit $\sigma$ als ihr Kehrwert $\rho$ auf, der spezifische Widerstand. Die Quadratwurzel aus $\sigma$ wird immer als Produkt mit dem Radius verwendet. Daher können wir den Effekt der Radius-Änderung durch Multiplikation von $\sigma$ mit dem Kehrbruch des Faktors für den äquivalenten Radius multiplizieren. Wir bekommen:

\begin{align*} \sqrt{\sigma}a &= \sqrt{\sigma_e}a\left(\frac{b}{a}\right)^{\left (1-\frac{1}{\varepsilon_r}\right)} \\ \frac{\sigma_e}{\sigma} &= \frac{1}{\left(\frac{b}{a}\right)^{2\left (1-\frac{1}{\varepsilon_r}\right)}} = \left(\frac{a}{b}\right)^{2\left(1-\frac{1}{\varepsilon_r}\right)} \\ \end{align*}

In pymininec brauche ich allerdings den äquivalenten Radius bei der Skin Effekt Berechnung nicht zu korrigieren da ich dort einfach den Original-Radius verwenden kann.

In den folgenden drei Grafiken vergleiche ich diverse Implementierungen von isolierten Drähten mit den Messungen von Lamensdorf [4]. Der Draht hat einen Durchmesser (nicht Radius) von 0.25 Zoll bei einer Frequenz von 600MHz. Ich verwende nur die Messungen mit $\varepsilon_r = 3.2$, die höheren Permittivitäten sind in der Praxis unrealistisch hoch und die Ergebnisse dürften durch das unrealistische experimentelle Setup abweichen. Ich plotte den Real- (elektrische Leitfähigkeit, auch Konduktivität) und Imaginärteil (Suszeptanz) der Admittanz (Kehrbruch der Impedanz). Die drei Grafiken sind:

Ohne Isolation
Durchmesser des Drahtes mit Isolation 0.375 Zoll (Isolation 0.125 Zoll)
Durchmesser des Drahtes mit Isolation 0.5 Zoll (Isolation 0.25 Zoll)

Die unterschiedlichen Kurven in den Plots sind:

ASAP mit zwei, vier, acht und 16 segmenten für eine Dipolhälfte, die Kurven sind mit 'AS02', 'AS04', 'AS08' and 'AS16' markiert, es werden nicht alle Kurven in jedem Plot gezeigt.
4NEC2 das nach Stearns [6] die Formel von Cebik [9] implementiert, diese Kurven sind mit '4N2' markiert
NEC-2 mit Wu's [7] Formel mit äquivalentem Radius und verteilter Induktivität sowie mit Stearns [6] Formel für das äquivalente $\sigma$ markiert als 'NEC'
Meine pymininec Implementierung verwendet Wu's [7] Formel mit äquivalentem Radius und verteilter Induktivität (aber nicht die $\sigma$ Korrektur weil ich für die Skin-Effekt Berechnung den Original-Radius verwenden kann), markiert als 'mini'
Die Messungen von Lamensdorf [4] sind die blauen und orangen Blasen (jeweils für den Real- und Imaginärteil)
In den Grafiken mit Isolation zeige ich auch die Resultate von EZNEC, die Kurven sind ausgegraut und können in der Legende auf der rechten Seite eingeschaltet werden

Die Formeln von Cebik [9] und Wu [7] sind sehr ähnlich, die Kurven '4N2' und 'NEC' sind fast identisch.

Die EZNEC Kurven und einige der ASAP Kurven sind ausgegraut, weil sie numerische Probleme aufgrund der schlecht konditionierte Impedanzmatrix durch die Formeln von Richmond zeigen. Es ist interessant zu sehen, dass EZNEC auch numerische Probleme bei kleinen Segmenten hat. Die ausgegrauten Kurven können mit Klick in die Legende auf der rechten Seite eingeschaltet werden.

Die Segmentierung in allen Beispielen sind 21 Segmente für einen Draht. Das bedeutet recht kurze Segmente für die kleineren Drahtlängen. Dies triggert eine Warnung zu kurzen Segmenten in EZNEC für die ersten drei Punkte (bis zu einschließlich $\frac{4.2}{32} = 0.06558125$ Wellenlängen) aber nicht für den klar numerisch ungültigen 6. Wert in der ersten Grafik mit Isolation und für die numerisch instabilen Werte in der zweiten Grafik mit Isolation.

Es sei darauf hingewiesen dass die Messungen von Lamensdorf [4] an Monopol-Antennen durchgeführt wurden. Ich simuliere Dipole. Die Impedanz eines Monopols ist die Hälfte von der eines Dipols, daher multipliziere ich die Admittanz mit dem Faktor zwei. Alle Admittanzen sind mit milli-Siemens, früher auch milli-mhos (wobei mho ohm rückwärts ist) mit dem Symbol ℧ zur Zeit des Lamensdorf [4] Papers.

Plotting Antenna Pattern

Ralf Schlatterbeck 2024-08-31 19:30

Ich bin der Autor eines Antennen-Plot-Programms namens plot-antenna. Es ist in Python geschrieben und hat Backends für Matplotlib und Plotly. Während ich das hier schreibe kann die Version auf pypi noch nicht unterschiedliche Frequenz-Auflösungen für Gewinn und Impedanz-Plots. Das Matplotlib Backend ist besser für die Erzeugung von Grafiken für gedruckte Dokumentation während Plotly sich eignet, interaktive Grafiken für das Web zu erzeugen. Die Doku von plot-antenna enthält einige Matplotlib Beispiele daher gehe ich hier mehr ins Detail was das Plotly Backend betrifft.

Vor kurzem habe ich Unterstützung für das "Antenna Scatterers Analysis Program" ASAP eingebaut. Das ist ein Antennen-Simulationsprogramm aus einer Zeit vor dem Fortran-77 Standard aus dem Jahr 1974. Es ist heute noch interessant weil es eine Formulierung der Momentenmethode (method of moments) verwendet, die sich sowohl von NEC [3] als auch Mininec [4] unterscheidet. Und ASAP hat zuerst eine Simulation von isolierten Drähten in einer Antenne unterstützt.

Im Folgenden verwende ich eine 3-Elemente Yagi-Uda Antenne die sich in der Beispiel-Sektion von ASAP findet. Die verwendete Eingabedatei exportiert Gewinn-Daten für 280 MHz bis 305 MHz in 5 MHz Schritten. Zusätzlich exportiert es Impedanzwerte in 1 MHz Schritten.

Für die Erstellung von Gewinn-Grafiken unterstützt plot-antenna unterschiedliche Skalierungsmethoden. Der Standard ist ARRL Skalierung [5] die ich auch im folgenden Beispiel verwende.

Im 3D-Plot kann man die Frequenz auf der rechten Seite in der Legende der Grafik umschalten. Ein Klick auf eine Frequenz ändert das Display auf die Gewinn-Grafik für diese Frequenz. Zusätzlich kann man beim Drüberfahren mit der Maus die Gewinnwerte an unterschiedlichen Positionen sehen. Es ist möglich mit dem Maus-Rad zu zoomen und den Plot durch klicken und ziehen zu drehen.

Natürlich kann plot-antenna auch konventionelle 2D-Plots für Azimuth und Elevation erzeugen:

Im 2D-Display kann man wieder in der Legende auf der rechten Seite die Plots für verschiedene Frequenzen einblenden, im 2D-Plot auch mehrere gleichzeitig. Damit ist es einfacher, Gewinnwerte zu vergleichen als im 3D-Plot wo immer nur der Gewinn für eine Frequenz gleichzeitig sichtbar ist. Ähnlich wie im 3D-Plot kann man die Gewinnwerte beim Drüberfahren mit der Maus sehen. Es ist auch hier möglich zu zoomen (obwohl das in einem Polarplot nicht sehr sinnvoll ist) und die Anzeige mit den "Home" Knopf im Menu auf die Standardansicht zurückzusetzen.

Oft ist es nützlich das Stehwellenverhältnis (voltage standing wave ratio, VSWR) und die Antennen-Impedanz anzeigen zu lassen. Das geht mit dem VSWR-Plot. In der folgenden Grafik habe ich die Anzeige der Impedanz mit der Option --swr-show-impedance eingeschaltet. Die vertikale grüne Linie zeigt die Position des minimalen Stehwellenverhältnisses an. Auch hier sehen wir die geplotteten Werte beim Drüberfahren mit der Maus. Wir können heranzoomen indem wir ein Rechteck mit der Maus selektieren oder mit dem +/- Knöpfen im Menu. Auch hier kann die Standardansicht mit dem "Home" Knopf wieder hergestellt werden.

Nicht zuletzt kann die Antennen-Geometrie angezeigt werden. Dies ist wieder ein 3D-Plot wo wir die Position der Antennen-Elemente beim Drüberfahren mit der Maus sehen können. In der Legende auf der rechten Seite kann man z.B. die Sichtbarkeit des Speisepunkts (in Orange) ausschalten. Andere Antennen können z.B. zusätzlich eine Ground-Plane (Boden) anzeigen oder Impedanz-Lasten auf der Antenne. Die Geometrie-Ansicht ist sehr nützlich um zu checken ob eine Antenne korrekt modelliert wurde (z.B. ob der Speisepunkt an der korrekten Position sitzt).

Das Programm plot-antenna unterstützt derzeit Ausgabe von ASAP, vom Original Mininec (das in Basic geschrieben ist), von meiner Python-Version pymininec, und der Ausgabe von NEC-2. Das Format wird automatisch erkannt. Zusätzlich kann mit dem Begleitprogramm plot-eznec der Gewinn einer mit EZNEC simulierten Antenne angezeigt werden. Während ich das hier schreibe gibt es eine Version die auch VSWR-Plots für EZNEC Daten anzeigen kann indem die beim Plotten von VSWR-Daten mit EZNEC erzeugte Datei LastZ.txt verwendet wird.

Skin Effekt Belag und Impedanz: Update

Ralf Schlatterbeck 2024-07-28 20:45

In meinem letzten Post [1] habe ich mich gefragt, warum meine Implementierung von Skin Effekt Belag (Drähte mit einem endlichen Widerstand) nicht mit den Resultaten von NEC-2 übereinstimmten. Es stellt sich heraus dass ich falsch lag in der Annahme dass Skin-Effekt Belag nur eine rein ohmsche Komponente hat, durch den Skin Effekt ergibt sich in Wirklichkeit eine relativ große reaktive Komponente.

Hinweis: Dieser Artikel enthält ziemlich viel Mathematik, vielleicht möchten Sie ja ans Ende scrollen und checken wie die neue Implementierung in pymininec sich verglichen mit meinem letzten Artikel verhält wo nur ein rein ohmscher Widerstand angenommen wurde.

Der englische Wikipedia Artikel zum skin effect enthält die korrekte Formel:

\begin{equation*} \newcommand{\Int}{{\mathrm\scriptscriptstyle int}} \newcommand{\ber}{\mathop{\mathrm{ber}}\nolimits} \newcommand{\bei}{\mathop{\mathrm{bei}}\nolimits} \end{equation*}

\begin{equation*} Z_\Int = \frac{k\rho}{2\pi r}\frac{J_0 (kr)}{J_1 (kr)} \end{equation*}

wobei

\begin{equation*} k = \sqrt{\frac{-j\omega\mu}{\rho}} \end{equation*}

und $r$ ist der Draht-Radius, $J_v$ sind die Bessel-Funktionen der ersten Gattung mit der Ordnung $v$. $Z_\Int$ ist die Impedanz pro Längeneinheit des Drahtes.

NEC-2 verwendet eine andere Formel [2] (S.75):

\begin{equation*} Z_i = \frac{j\Delta_i}{a_i} \sqrt{\frac{\omega\mu}{2\pi\sigma}} \left[\frac{\ber (q) + j\bei (q)} {\ber^\prime(q)+ j\bei^\prime(q)}\right] \end{equation*}

mit

\begin{equation*} q = (\omega\mu\sigma)^{1/2} a_i \end{equation*}

und $a_i$ ist der Draht-Radius, $\Delta_i$ die Länge des Drahtsegments, $\sigma$ die Elektrische Leitfähigkeit des Drahtes und $\ber$ und $\bei$ die Kelvin-Funktionen. Im folgenden werde ich wieder $r$ für den Draht-Radius, und nicht $a_i$ verwenden. Nachdem $\ber^\prime$ die Ableitung von $\ber$ ist, können wir $\ber$ als $\ber_0$ (Ordnung 0) und $\ber^\prime$ als $\ber_1$ (Ordnung 1) schreiben und analog für $\bei$.

Im folgenden bin ich etwas ungenau in der Verwendung der magnetischen Permeabilität im Freiraum $\mu_0$ vs. der Permeabilität eines Drahtes $\mu$, für nicht ferromagnetische Drähte sind diese (fast) identisch: $\mu = \mu_r * \mu_0$ mit $\mu_r\approx 1$ für nicht ferromagnetische Materialien.

Der NEC-2 code verwendet die Fortran Funktion ZINT mit zwei Parametern, der Leitfähigkeit multipliziert mit der Wellenlänge $\lambda$ (Parameter SIGL) und der Radius geteilt durch die Wellenlänge $\lambda$ (Parameter ROLAM). Neben der Approximation der Kelvin-Funktionen (der Term mit dem Bruch aus Kelvin-Funktionen ist BR1), berechnet ZINT:

X = SQRT (TPCMU * SIGL) * ROLAM
ZINT = FJ * SQRT (CMOTP / SIGL) * BR1 / ROLAM

mit TPCMU = 2368.705, CMOTP=60.0, und FJ ist $j$. Das berechnete X wird der Berechnung der Kelvin-Funktionen übergeben mit BR1 als Resultat des Bruch-Terms. Das berechnete ZINT enthält nicht die Segment-Länge, ist also äquivalent mit der Wikipedia Definition in pro Längeneinheit. Mit den Reverse-Engineerten magischen Zahlen

CMOTP $\approx \frac{\mu_0 c}{2\pi}$

TPCMU $\approx 2\pi c\mu_0$

($c$ ist die Lichtgeschwindigkeit) bekommen wir für X aus dem Code (was wir jetzt $q$ nennen um mit der obigen Formel aus [2] konsistent zu sein):

\begin{align*} q &= \frac{r}{\lambda}\sqrt{2368.705\sigma\lambda} \\ &= r\sqrt{\frac{2368.705\sigma}{\lambda}} \\ &\approx r\sqrt{\frac{2\pi c\mu_0\sigma}{\lambda}} \\ &= r\sqrt{\frac{2\pi c\mu_0\sigma f}{c}} \\ &= r\sqrt{\frac{2\pi \mu_0\sigma \omega}{2\pi}} \\ &= r\sqrt{\omega\mu_0\sigma} \\ \end{align*}

Was mit der obigen Definition von $q$ in [2] übereinstimmt.

Für $Z_\Int$ der Fortran-Implementierung bekommen wir:

\begin{align*} Z_\Int &= \frac{j}{r} \sqrt{\frac{\mu_0 c}{2\pi\sigma\lambda}} \left[\frac{\ber_0 (q) + j\bei_0 (q)} {\ber_1(q)+ j\bei_1(q)}\right] \\ &= \frac{j}{r} \sqrt{\frac{\mu_0 c f}{2\pi\sigma c}} \left[\frac{\ber_0 (q) + j\bei_0 (q)} {\ber_1(q)+ j\bei_1(q)}\right] \\ &= \frac{j}{r} \sqrt{\frac{\mu_0 f}{2\pi\sigma}} \left[\frac{\ber_0 (q) + j\bei_0 (q)} {\ber_1(q)+ j\bei_1(q)}\right] \\ &= \frac{j}{r} \sqrt{\frac{\mu_0 \omega}{2\pi 2\pi\sigma}} \left[\frac{\ber_0 (q) + j\bei_0 (q)} {\ber_1(q)+ j\bei_1(q)}\right] \\ &= \frac{j}{2\pi r} \sqrt{\frac{\mu_0 \omega}{\sigma}} \left[\frac{\ber_0 (q) + j\bei_0 (q)} {\ber_1(q)+ j\bei_1(q)}\right] \\ \end{align*}

Diese Formel stimmt NICHT mit der Formel für $Z_i$ aus [2] weiter oben überein. Es sieht so aus als ob die Formel aus [2] (S. 75) die Frequenz $f$ statt der Kreisfrequenz $\omega$ unter der Wurzel haben sollte. Die aus der NEC-2 Implementierung abgeleitete Formel ist aber konsistent mit der obigen Formel aus Wikipedia.

Der Term vor dem Bruch mit den Kelvin-Funktionen unterscheidet sich um einen Faktor von der Wikipedia Formel (Ich lasse $\Delta_i$ weg und habe so eine Formel pro Längeneinheit wie in Wikipedia). Ausserdem gilt $\rho=1/\sigma$:

Aus dem Code abgeleitet:

\begin{equation*} \frac{j}{2\pi r}\sqrt{\frac{\omega\mu}{\sigma}} \end{equation*}

Aus Wikipedia:

\begin{equation*} \frac{k\rho}{2\pi r} = \sqrt{\frac{-j\omega\mu}{\rho}}\frac{\rho}{2\pi r} = \frac{1}{2\pi r}\sqrt{-j\omega\mu\rho} = \frac{1}{2\pi r}\sqrt{\frac{-j\omega\mu}{\sigma}} \end{equation*}

Das ist ein Faktor von

\begin{equation*} \frac{\sqrt{-j}}{j} = -\frac{1+j}{\sqrt{2}} \end{equation*}

weil $\sqrt{-j} = \frac{1-j}{\sqrt{2}}$

Für die Kelvin-Funktionen gilt die folgende Äquivalenz zu den Bessel Funktionen [3] (10.61, Kelvin Functions):

\begin{equation*} \ber_v x + j \bei_v x = J_v \left(xe^{3\pi j/4}\right) = e^{v\pi j} J_v \left(xe^{-\pi j / 4}\right) \end{equation*}

Mit den Spezialfällen der Ordnung 0 und 1:

\begin{align*} \ber_0 x + j \bei_0 x &= e^{0} J_0 \left(xe^{-\pi j / 4}\right) \\ &= J_0 \left(\frac{1-j}{\sqrt{2}}x\right) \\ \ber_1 x + j \bei_1 x &= e^{j\pi} J_1 \left(xe^{-\pi j / 4}\right) \\ &= -J_1 \left(\frac{1-j}{\sqrt{2}}x\right) \\ \end{align*}

Die Wikipedia Formel übergibt $kr$ an die Bessel-Funktionen:

\begin{equation*} kr = r\sqrt{\frac{-j\omega\mu}{\rho}} = r\sqrt{-j\omega\mu\sigma} \end{equation*}

während NEC-2 folgendes übergibt:

\begin{equation*} q = r\sqrt{\omega\mu\sigma} \end{equation*}

Diese unterscheiden sich um den Faktor $\sqrt{-j}$. Das bedeutet wir übergeben einen Faktor um $\sqrt{2}$ zu groß an die Bessel-Funktionen. Daraus ergibt sich wegen der Symmetrie der Bessel Funktionen und er Tatsache dass die Funktionen im Nenner die Ableitung der Funktionen im Zähler sind, dass das Ergebnis um einen Faktor von $\sqrt{2}$ zu groß ist was den Faktor im ersten Term ausgleicht.

Nach der Implementierung der neuen Formel in pymininec verhält sich der reaktive Teil der Impedanz ähnlich wie bei NEC-2 und der ohmsche Widerstand ist noch näher an den Ergebnissen von NEC als vorher.

Auch wenn man den Imaginärteil (die Reaktanz) plottet sind die Resultate jetzt wie erwartet. Die Werte von MMANA sind immer noch überraschend, eine mögliche Ursache ist dass ich MMANA unter WINE in Linux verwendet habe, die ohmschen Werte der Impedanz die von MMANA berechnet wurden schauen allerdings korrekt aus.

Skin Effekt Belag und Impedanz

Ralf Schlatterbeck 2024-07-04 19:45

Vor kurzem habe ich für pymininec, meine Re-Implementierung des Mininec Antennen-Modellierungs-Codes Skin Effekt Belag (Drähte mit einem endlichen Widerstand) implementiert. Das Original-Mininec hatte dieses Feature nicht.

Um die Implementierung zu verifizieren entschied ich einer uralten Publikation von R. P. Haviland, W4MB [1] zu folgen wo er unter anderem den Effekt von Draht-Widerstand auf die Speisepunkt-Impedanz von Antennen untersucht. Er hat zwei Grafiken wo der Real- und Imaginärteil der Speisepunkt-Impedanz eines Dipols gegen das Verhältnis von Länge zu Durchmesser des Drahtes geplottet wird. Er erwähnt zwar Mininec im Text aber bei genauerer Analyse ist es nicht ganz klar ob für die Plots der Speisepunkt-Impedanz wirklich Mininec oder etwas anderes verwendet wurde. Ich habe versucht, seine Grafiken mit PyNEC, einem Python wrapper für eine C++ Implementierung von NEC-2, und meiner pymininec Implementierung von Mininec zu reproduzieren. Die Resultate von Haviland haben die gleiche Tendenz wie die folgenden NEC Resultate aber sind nicht gleich.

Die X-Achse dieser Grafiken ist logarithmisch und die Drähte werden nach rechts dicker (die X-Achse ist das Verhältnis von Länge zu Durchmesser des Drahtes). In allen Fällen vergleichen wir Drähte ohne Widerstand mit Kupferdrähten. "Resistance" ist der (ohmsche) Widerstand (Realteil der Impedanz) und "Reactance" ist der Blindwiderstand (Reaktanz), der Imaginärteil der Impedanz.

Der ohmsche Widerstand ist nicht überraschend: Der Realteil der Antennen-Impedanz steigt mit dem Draht-Widerstand (weil der Draht dünner wird).

Der interessante Teil ist der Imaginärteil (die Reaktanz) der Antenne. Das NEC-Resultat besagt dass die Reaktanz mit dem Widerstand des Drahtes steigt (nach links) während meine pymininec Implementierung ein leichtes Sinken der Reaktanz vorraussagt.

Um das genauer zu untersuchen habe ich mich entschlossen den spezifischen Widerstand des Drahtes ohne Modifikation des Durchmessers zu ändern um zu zeigen dass der Effekt unabhängig vom Drahtdurchmesser ist. Der Radius ist im Folgenden 66µm bei einer Wellenlänge von 1m, $\frac{l}{d}\approx 7576$.

Zusätzlich zu PyNEC und pymininec habe ich auch Resultate für MMANA-GAL und EZNEC in die Grafik aufgenommen. MMANA basiert ursprünglich auf Mininec. Die EZNEC Resultate stimmen ziemlich genau mit PyNEC überein (beides verwendet den gleichen Modellierungscode, einmal in Fortran, einmal in C++) aber MMANA hält die Reaktanz komplett konstant bei $42.63\Omega$ über alle spezifischen Widerstände des Drahtes. Die interessantesten Werte des spezifischen Widerstands sind ganz links die ersten drei großen Punkte, der erste ist ein idealer Draht (spezifischer Widerstand 0) gefolgt von Kupfer und Alu.

Ich habe derzeit keine Ahnung welche Resultate der physischen Wirklichkeit am nächsten kommen.

Schwebung

Ralf Schlatterbeck 2024-05-11 11:00

Musiker kennen den Effekt beim Stimmen von Instrumenten: Wenn zwei Instrumente nicht ganz den selben Ton spielen klingt das als ob es eine periodische Variation der Lautstärke gibt. Die Periode wird größer je besser die beiden Instrumente gestimmt sind und stoppt wenn beide den selben Ton spielen. Wir nennen diesen Effekt "Schwebung".

Funkamateure kennen den "umgekehrten" Fall dieses Effekts von der Amplitudenmodulation einer Radio-Frequenz mit einer Sinusschwingung im Audio-Bereich: Das Resultat sind zwei Radio Frequenzen, eine um die Audio-Frequenz niedriger als die Modulationsfrequenz, eine um die Audio Frequenz höher.

Amplitudenmodulation ist eine einfache Multiplikation: Wir multiplizieren die Radiofrequenz mit dem Audio, wir haben:

$$m(t) = A(t) \cdot\cos(2\pi f t)$$

Wobei $t$ die Zeit ist, $m(t)$ das amplitudenmodulierte Signal, $\cos(2\pi f_r t)$ ist das Radiofrequenz-Signal mit der Frequenz $f_r$ und $A(t)$ ist das Audiosignal. Für eine einfache Audio Sinusschwingung bekommen wir:

$$m(t) = \cos(2\pi f_a t)\cdot\cos(2\pi f_r t)$$

wobei $f_a$ die Audio-Frequenz ist.

Wir möchten jetzt zeigen, dass beide Effekte das gleiche Ding sind. Wir brauchen keine Radiofrequenzen, wir könnenn das gleiche einfach im Audio-Spektrum tun.

Aus Wikipedia bekommen wir eine trigonometrische Identität über Produkte von Winkelfunktionen, daraus nehmen wir die für zwei Cosinus Funktionen:

$$\displaystyle\cos\theta \cos\phi = \frac{\cos(\theta-\phi) + \cos(\theta+\phi)}{2}$$

Wenn wir die auf unser Modulation-Beispiel anwenden bekommen wir:

$$\displaystyle \begin{align} \ \cos(2\pi f_a t) \cos(2\pi f_r t) & = \frac{\cos(2\pi f_r t-2\pi f_a t) + \cos(2\pi f_r t+2\pi f_a t)}{2} \\ & = \frac{\cos(2\pi(f_r-f_a)t)) + \cos(2\pi(f_r+f_a)t)}{2} \\ \end{align} $$

Die Multiplikation von zwei Cosinus-Schwingungen mit unterschiedlichen Frequenzen resultiert in der Summe von zwei Cosinus-Schwingungen, eine davon die Differenz der beiden ursprünglichen Frequenzen, die andere die Summe.

Egal wie wir das modulierte Signal erzeugen, durch Multiplikation von zwei Signalen oder durch Addition von zwei Signalen, wir bekommen das selbe Resultat.

Im folgenden Beispiel multiplizieren wir ein 3Hz Signal (die Schwebungsfrequenz) mit einem 440Hz Signal (Das ist die Note "A" in der Musik). Wir bekommen wie erwartet ein moduliertes Signal das die Lautstärke mit der Schwebungsfrequenz ändert.

In [1]:

import numpy as np
# Drei Sekunden with 18000 Abtastungen
t = np.linspace (0, 3, 18000)
# Produkt von einer 3Hz Schwingung und einem 440Hz Ton
y = np.cos (2*np.pi*3*t) * np.cos (2*np.pi*440*t)
# Wir plotten nur 1/9 der drei Sekunden = 1/3s = eine Periode des 3Hz Signals
limit = int (2000)

In [2]:

import plotly.graph_objs as go
import plotly.io as pio
from plotly_fix_fn import setup_plotly
setup_plotly ('content/2024/05/plotly', '11-de')
show_opt = dict (include_plotlyjs = 'directory')

fig = go.Figure ()
fig.add_trace (go.Scatter (x = t [:limit], y = y [:limit]))
fig.update_layout (title = "Moduliertes Signal", xaxis_title = "Zeit (s)", yaxis_title = "Amplitude")
fig.show (** show_opt)

Nun machen wir eine Fourier-Transformation des Resultats um das Frequenz-Spektrum zu bekommen. Wir sehen dass das Resultat der Multiplikation aus zwei Frequenzen besteht: Ein 3Hz unter 440Hz, eine um 3Hz darüber.

In [3]:

from numpy.fft import rfft, fftfreq
fft   = rfft (y, norm = 'forward')
freqs = fftfreq (len (y), d = 3/18000)
fig = go.Figure ()
fig.add_trace (go.Scatter (x = freqs [:5000], y = np.abs (fft [:5000])))
fig.update_layout (title = "Frequenz-Spektrum", xaxis_title = "Frequenz (Hz)")
fig.update_xaxes (range = [430, 450])
fig.show (** show_opt)

Wir können das selbe Resultat erzielen wenn wir einfach die beiden Frequenzen addieren (und das Ganze durch zwei teilen damit wir die selbe Amplitute haben)

In [4]:

y = 0.5 * (np.cos (2*np.pi*(440-3)*t) + np.cos (2*np.pi*(440+3)*t))
fig = go.Figure ()
fig.add_trace (go.Scatter (x = t [:limit], y = y [:limit]))
fig.update_layout (title = "Moduliertes Signal", xaxis_title = "Zeit (s)", yaxis_title = "Amplitude")
fig.show (** show_opt)

Schließlich verifizieren wir dass dieses Vorgehen das selbe Audiospektrum erzeugt wie vorher die Multiplikation, wir machen eine weitere Fourier-Transformation des neuen Signals:

In [5]:

fft   = rfft (y, norm='forward')
freqs = fftfreq (len (y), d = 3/18000)
fig = go.Figure ()
fig.add_trace (go.Scatter (x = freqs [:5000], y = np.abs (fft [:5000])))
fig.update_layout (title = "Frequenz-Spektrum", xaxis_title = "Frequenz (Hz)")
fig.update_xaxes (range = [430, 450])
fig.show (** show_opt)

Was wir also beim Stimmen von Instrumenten hören ist keine Illusion: Die Addition von zwei Schwingungen mit unterschiedlichen Frequenzen ist das selbe wie die Modulation eines Tons mit einer Schwebungsfrequenz.

Chipselect Aussetzer beim Schreiben auf SPI NOR des IMX-28

Ralf Schlatterbeck 2024-02-01 11:30

Für einen Kunden pflege ich eine Linux Umgebung (inclusive Kernel) für ein kundenspezifischen Board das den IMX-28 ARM Prozessor verwendet. Für den User Space verwende ich ELBE, das Embedded Linux Build Environment das es erlaubt, das User Space Dateisystem aus eine Debian Linux Release zu bauen. Bei der Migration von einer früheren Debian Version zum derzeitigen stabilen Release mit dem Code-Namen "Bookworm" brauchte es ein Kernel Upgrade weil Systemd einige Dinge benötigt die nicht in früheren Kernels enthalten sind.

Nach dem Upgrade von Kernel 4.14 auf 6.1 (beide eine longterm Support Version) funktionierte unser SPI NOR Flash nicht mehr. Das JFFS2 Dateisystem produzierte Fehlermeldungen die nahelegten dass wir ein Problem beim Schreiben in das NOR Flash hatten (Einige der Fehlermeldungen zeigten dass eine 0 erwartet aber 0xFF gelesen wurde was zum Schluss führte dass diese Bytes nicht geschrieben wurden).

Nach Überprüfung der Taktrate des SPI Busses beim Schreiben des Flash mit der Beobachtung dass sich diese nicht geändert hatte schaute ich mir die Signale des SPI Busses an.

Bild Aussetzer im Chip Select — Aussetzer im Chip Select bei der Übertragung

Bild Aussetzer im Chip Select höhere zeitliche Auflösung — Aussetzer im Chip Select bei der Übertragung: Höhere zeitliche Auflösung

In den Oszilloscop Bildern (Blau ist Chip Select, Gelb ist der SPI Takt) sieht man dass es einen kurzen Aussetzer im Chip Select gibt. Das hat den Effekt dass das "Page Program" Kommando an das NOR Flash frühzeitig beendet wird. Aus dem Bild mit der höheren zeitlichen Auflösung sieht man dass dieser Aussetzer mitten in einem Bit der Übertragung stattfindet. Die exakte Position variiert leicht von Übertragung zu Übertragung (Der Aussetzer trifft nicht jedesmal das selbe Bit aber er trifft das selbe Byte).

Weiter Experimente zeigten dass das Problem nur bei Übertragungen mit direktem Speicherzugriff (DMA) auftritt, nicht wenn kein DMA (sondern programmierter Input-Output, PIO) verwendet wurde. Der SPI Treiber optimiert Übertragungen: Nur längere Übertragungen werden mit DMA durchgeführt, kürzere Übertragungen verwenden PIO weil der DMA Overhead dann den Zeit-Gewinn überwiegt. Bei genauerer Inspektion der DMA Implementierung für den IMX-28 in drivers/dma/mxs-dma.c stellte sich heraus dass es eine Änderung bei den Flags des DMA Transfers gab, die alte Version hat:

if (flags & DMA_CTRL_ACK)
        ccw->bits |= CCW_WAIT4END;

Die neue Version hat:

if (flags & MXS_DMA_CTRL_WAIT4END)
        ccw->bits |= CCW_WAIT4END;

Das Bit CCW_WAIT4END (Übersetzt aus der IMX-28 Prozessor Referenz): "Der Wert eins gibt an dass der Kanal auf das Ende des Kommandosignalsvom APBH Device zum DMA wartet bevor ein neues DMA Kommando gestartet wird". Dieses Flag wird für den DMA Transfer zum NOR Flash gesetzt. Eine Besonderheit des IMX-28 ist dass auch Register Einstellungen über DMA übertragen werden können. Vor dem eigentlichen Speichertransfer finden ein Register Transfer statt der das IGNORE_CRC Bit im Prozessor SSP Register setzt. Aus der Prozessor Referenz zum IGNORE_CRC Bit (übersetzt): "Im SPI/SSI Modus: Wenn auf 1 gesetzt, setze den Chip Select (SSn) Pin zurück nachdem das Kommando ausgeführt wurde".

Was also passiert wenn das CCW_WAIT4END nicht gesetzt wird, ist dass der Prozessor das Chip-Select zu früh ausschaltet.

Die DMA Übertragung die vom SPI Treiber in drivers/spi/spi-mxs.c initiiert wird setzt das neue MXS_DMA_CTRL_WAIT4END nicht, so dass das flag CCW_WAIT4END für die Übertragung nicht gesetzt wird.

Der Fix ist damit offensichtlich, wir müssen nur dieses Flag bei der Initiierung der DMA-Übertragung im SPI Treiber in drivers/spi/spi-mxs.c setzen.

Linux Kernel route_localnet Einstellung

Ralf Schlatterbeck 2024-01-18 12:00

Der Linux Kernel hat Einstellungen für Netzwerk Schnittstellen die es erlauben, das lokale Netz (127.0.0.0/8 für IP Version 4, auch oft Loopback Adresse bezeichnet weil sie typischerweise der lo Loopback Schnittstelle zugewiesen ist) zu routen. Diese Einstellung kann man über die Datei

/proc/sys/net/ipv4/conf/eth0/route_localnet

für die typische eth0 Netzwerk Schnittstelle erreichen. Die Einstellung kann man ändern, indem man dort entweder 0 (Standardeinstellung) für kein Routing oder 1 für Routing hineinschreibt. Einstellungen kann man bei den üblichen Linux Distributionen persistent machen, indem man einen Eintrag in /etc/sysctl.conf vornimmt. Dabei wird das Präfix /proc/sys/ vom obigen Pfad entfernt und alle / Zeichen durch einen Punkt ersetzt, das Resultat für obiges Beispiel wäre dann:

net.ipv4.conf.eth0.route_localnet = 1

Was kann man mit dieser Einstellung tun?

Im Folgenden verwende ich Beispiele für die IP Version 4, das gleiche kann man aber für die IP Version 6 local host Adresse ::1 erreichen.

Im Linux Netzwerk-Stack werden IP Pakete die entweder von einer Adresse aus dem lokalen Netz 127.0.0.0/8 kommen oder eine Zieladresse dort haben verworfen, wenn sie über eine Schnittstelle kommen die nicht für diesen Adressbereich konfiguriert ist. Wir verifizieren das im Folgenden zuerst mal.

Für einen einfachen UDP Server in Python der UDP Pakete empfängt und deren IP Adresse und Port (und Länge) ausdruckt können wir das folgende (gekürzte) Script verwenden:

import socket
[...]
sock = socket.socket (socket.AF_INET, socket.SOCK_DGRAM,
socket.IPPROTO_UDP)
sock.bind ((args.address, args.port))
while True:
    buf, adr = sock.recvfrom (1024)
    l = len (buf)
    print ('Received %(l)s bytes from %(adr)s' % locals ())

Alle Werte in der Variable args kommen von der Kommandozeile mit der Klasse ArgumentParser. Die Adresse auf die gebunden wird (bind) ist im folgenden standardmässig 0.0.0.0 was bedeutet dass Pakete von allen Netzwerk-Schnittstellen empfangen werden.

Wir testen das erstmal mit einem einfachen Client-Programm. Im folgenden haben wir zwei Maschinen, .9 sendet und .5 empfängt. Zum Versenden von UDP Paketen können wir den folgenden einfachen Code verwenden:

import socket
sock = socket.socket (socket.AF_INET, socket.SOCK_DGRAM, socket.IPPROTO_UDP)
sock.sendto (b'This is a test', (args.address, args.port))

Auch hier werden die Parameter auf der Kommandozeile angegeben. Wenn wir von .9 an .5 senden bekommen wir vom obigen Server

Received 14 bytes from ('XXXXXXX.9', 38232)

wobei der der Quell-Port 38232 eine zufällige vom Linux Kernel erzeugte Portnummer ist und XXXXXXX die oberen Bytes der IP-Adresse des Senders repräsentiert. Wenn wir auf der empfangenden Maschine an 127.0.0.1 senden bekommen wir:

Received 14 bytes from ('127.0.0.1', 58196)

Wir sehen dass der Server sowohl auf der Ethernet Schnittstelle als auch auf der Loopback Schnittstelle empfängt.

Soweit nichts neues. Schauen wir mal was passiert wenn wir Adressen fälschen. Um beliebige IP Pakete (oder sogar Nicht-IP Pakete) zu erzeugen verwende ich im Folgenden das Python Paket Scapy.

Ein Ausschnitt aus dem Script:

p = Ether (dst = args.mac_address, src = args.mac_address) \
    / IP (dst = args.address, src = args.source_address) \
    / UDP (dport = args.port, sport = args.source_port) \
    / Raw (b'This is a test')
sendp (p, iface = args.interface)

Das baut ein Ethernet Paket (Ether) mit IP Inhalt. Das IP Paket wiederum enthält ein UDP Paket was wiederum einen Raw Zeichenfolge enthält. Das Paket wird über einen Netzwerk Layer 2 Socket versendet. Wiederum kommen die Parameter von der Kommandozeile. Standard für alle Ports ist 22222, die Quell-IP ist XXXXXXX.9 und die Ziel-IP ist XXXXXXX.5. Ziel- und Quell MAC Adresse werden beide auf die MAC-Adresse der Zielmaschine gesetzt. Wir wiederholen was wir mit dem einfachen Client ausprobiert haben:

sudo scapyclient.py

Was in der folgenden Ausgabe resultiert:

Received 14 bytes from ('XXXXXXX.9', 22222)

Hier wird standardmässig immer Port 22222 verwendet es sei denn wir setzen den Port über eine Kommandozeilen-Option des Scripts. Wir müssen sudo verwenden weil nur der Superuser beliebige Pakete auf Layer-2 des Netzwerk-Stacks senden darf.

Was passiert nun wenn wir Pakete fälschen? Damit wir sehen was gesendet (und empfangen) wird starten wir tcpdump wie folgt:

sudo tcpdump -n -i eth0 udp port 22222

Natürlich muss das Argument der -i Option geändert werden wenn die Netzwerk-Schnittstelle anders heisst als eth0. Die Option -n sagt tcpdump dass keine DNS-Auflösung von IP Adressen erfolgen soll. Wir senden ein Paket mit einer localhost Adresse als Quell-Adresse:

sudo scapyclient.py --source-address=127.0.0.1

Auf beiden Maschinen sehen wir in der tcpdump Ausgabe

TT:TT:TT.TTTTTT IP 127.0.0.1.22222 > 10.23.5.5.22222: UDP, length 14

wobei TT:TT:TT.TTTTTT ein Zeitstempel ist. Das Paket geht also auf der Ethernet-Schittstelle der sendenden Maschine raus und wird von der empfangenden Maschine empfangen (weil die MAC-Adresse richtig ist) aber von unserem Server wird nichts ausgegeben.

Das gleich passiert wenn wir die Zieladresse ändern:

sudo scapyclient.py --address=127.0.0.1

Die tcpdump Ausgabe ist jetzt:

TT:TT:TT.TTTTTT IP 10.23.5.9.22222 > 127.0.0.1.22222: UDP, length 14

Wieder wird das von beiden Maschinen von tcpdump ausgegeben aber nichts wird von unserem Server ausgegeben.

Schauen wir nun was passiert wenn wir route_localnet für unsere Ethernet Schnittstelle einschalten (das folgende funktioniert natürlich nur als root):

echo 1 > /proc/sys/net/ipv4/conf/eth0/route_localnet

Wir fälschen wieder die Quelladresse:

sudo python3.11 scapyclient.py --source-address=127.0.0.1

Wieder wird das Paket von beiden laufenden tcpdump Prozessen ausgegeben aber wieder nicht von unserem Server. Aber wenn wir die Zieladresse fälschen:

sudo python3.11 scapyclient.py --address=127.0.0.1

Bekommen wir zusätzlich zur tcpdump Ausgabe auch eine Antwort von unserem Server:

Received 14 bytes from ('10.23.5.9', 22222)

Und bemerkenswert ist, dass das auch funktioniert wenn wir unserem Server sagen dass er nur auf der localhost Adresse hören soll:

python3 server.py --address=127.0.0.1

Der Server empfängt trotzdem das gefälschte Paket obwohl es nicht über die Loopback Schnittstelle empfangen wurde:

Received 14 bytes from ('10.23.5.9', 22222)

Sicherheit

Manche Anwendungen verwenden lokale Server die nur auf die Loopback Schnittstelle hören. Es wird erwartet dass diese nur Pakete über die Loopback Schnittstelle empfangen – gesendet von Prozessen die auf der selben Maschine laufen wie der Server.

Dies ist immer eine gefährliche Annahme, es sollten nie gefährliche Dinge getan werden die sich auf diese (Nicht-) Authentifizierung verlassen. Zumindest unter Linux gibt es eine gewisse Sicherheit solange die Einstellung route_localnet auf 0 gesetzt ist. Auf anderen Betriebssystemen ist das Glück vielleicht nicht so hold. Und auch unter Linux kann dieses Verhalten ausgeschaltet werden.

Warum wollen wir diese Einstellung vornehmen?

Manchmal sollen Verbindungen die auf der lokalen Maschine zu einem Port auf der lokalen Maschine aufgebaut werden zu einer anderen Maschine verbunden werden, wenn z.B. /etc/hosts folgenden Eintrag hat

127.0.0.1 some.remote.example.com

und eine Anwendung zu dieser Maschine verbinden soll können wir Firewall Regeln konfigurieren die die Pakete so umschreiben dass das funktioniert:

iptables -t nat -A OUTPUT -d 127.0.0.1 -p tcp --dport 443 \
    -j DNAT --to-destination XXXXXXX.5:1443
iptables -t nat -A POSTROUTING -d XXXXXXX.5 -p tcp --dport 1443 \
    -j SNAT --to-source XXXXXXX.9

Die erste Regel schreibt die Zieladresse um (und den Zielport) während die zweite Regel sicherstellt dass Antwortpakete ankommen indem die Quelladresse umgeschrieben wird. Aber ohne die Einstellung route_localnet für die ausgehende Ethernet-Schnittstelle wird das nicht funktionieren, die Pakete werden verworfen.

Aber Vorsicht: Weil diese Einstellung auch auf die empfangenen Pakete einen Einfluss hat ist zu beachten was im Abschnitt Sicherheit beschrieben wurde.

Optimierung von Fliesskomma-Problemen mit Evolutionären Algorithmen

Ralf Schlatterbeck 2024-01-07 18:00

[Dies ist ein erweiterter Abstract für einen geplanten Vortrag]

Evolutionäre Algorithmen (EA) sind ein Oberbegriff für Genetische Algorithmen (GA) und Verwandte Algorithmen. GA arbeiten üblicherweise mit Bits oder kleinen Integer-Zahlen.

Es gibt andere EA die direkt mit Fliesskomma-Zahlen arbeiten können, unter ihnen Differential Evolution (DE) [1] [2] [3].

Der Vortrag gibt eine Einführung in die Optimierung von Fließkomma Problemen anhand von Beispielen aus der Elektrotechnik sowie der Optimierung von Kurvenformen zur Ansteuerung von piezoelektrischen Inkjet Druckern. Bei diesen Druckern hängt die Form des gejetteten Tropfens (neben anderen Parametern wie der gejetteten Flüssigkeit) von der zur Ansteuerung verwendeten Kurvenform ab. Diese bestimmt die Qualität der gejetteten Tropfen.

Für die Software verwende ich die Python Bindings PGAPy [4] für das ursprünglich an den Argonne National Laboratories entwickelte "Parallel Genetic Algorithm Package" PGAPack [5]. Beide Open Source Pakete maintaine ich seit einigen Jahren. Unter anderen wurde diverse Algorithmen wie Differential Evolution und Strategien zur Optimierung von Multi-Objective Optimization (also Problemen mit mehreren Zielfunktionen mit NSGA-II [6]) neu implementiert.

Differential Evolution veranschaulicht

Ralf Schlatterbeck 2023-04-07 22:00

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.

\begin{equation*} \mathbf v_{i,g} = \mathbf x_{r_0,g} + F \cdot (\mathbf x_{r_1,g} - \mathbf x_{r_2,g}) \end{equation*}

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.

Notizen zur Vorzeitigen Konvergenz in Genetischen Algorithmen

Ralf Schlatterbeck 2023-01-06 18:45

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.

Vermeidung von Duplikaten

PGAPack implementiert einen Mechanismus zur Vermeidung von doppelt vorkommenden Genen (Individuen). In früheren Implementierungen war der Aufwand dafür quadratisch in der Populationsgröße $N$, also $O(N^2)$ (es wurde jedes neue Individuum mit allen Mitgliedern der aktuellen Population verglichen). In der aktuellen Version wird eine Hashfunktion verwendet um Duplikate zu entdecken, was den Aufwand auf $O(N)$ reduziert, also einen kleinen konstanten Aufwand für jedes neue Individuum.

Für benutzerdefinierte Datentypen bedeutet das, dass eine Hash Funktion für den Datentyp gebraucht wird, zusätzlich zu einer Gen-Vergleichsfunktion. Für die eingebauten Datentypen (binary, integer, real) sind diese Funktionen automatisch verfügbar.

Die Vermeidung von Duplikaten funktioniert gut für Integer und Binäre Datentypen, besonders wenn die verwendeten genetischen Operatoren mit hoher Wahrscheinlichkeit Duplikate erzeugen können. Es funktioniert nicht so gut mit dem real Datentyp weil Individuen normalerweise unterschiedlich sind, auch wenn sie genetisch sehr nah beieinanderliegen.

Restricted Tournament Replacement

Ein Algorithmus der ursprünglich von Harik [2], [3] Restricted Tournament Selection genannt wurde und später von Pelikan [4] mit Restricted Tournament Replacement (RTR) bezeichnet wurde, verwendet die alte und die neue Population um zu entscheiden, ob ein Kandidat in die neue Population aufgenommen wird. Der Algorithmus wählt eine Anzahl von Individuen der alten Population (das sogenannte Selektionsfenster) zufällig aus und vergleicht die Bewertung des Indivuums welches genetisch dem Kandidaten-Individuum am nächsten ist. Nur wenn das Kandidaten-Individuum besser ist als das genetisch ähnlichste aus dem Selektionsfenster wird es in die neue Generation übernommen.

Die Standardeinstellung für die Anzahl der Individuen im Selektionsfenster ist $\min (n, \frac{N}{20})$ [4] wobei $n$ die Genlänge und $N$ die Populationsgröße ist. Dies kann vom Benutzer geändert werden durch Aufruf der Funktion PGASetRTRWindowSize für PGAPack bzw. mit dem Konstruktor-Parameter rtr_window_size für PGApy.

Der RTR-Algorithmus braucht eine Ähnlichkeitsmetrik um zu entscheiden, wie ähnlich ein Individuum einem anderen ist. Die Vorgabe ist die Manhattan Distanz (für binäre Gene identisch mit der Hamming Distanz), also die Summe aller einzelnen Allele-Differenzen. Dies kann auf die Euclidsche Distanz geändert werden durch Setzen einer benutzerdefinierten Funktion für die genetische Distanz. Die Verwendung einer Euclidschen Ähnlichkeitsmetrik sieht man im Beispiel: Magisches Quadrat in PGApy wo man die Euclidsche Metrik über eine Kommandozeilenoption auswählen kann.

Restricted Tournament Replacement funktioniert nicht nur für binäre und Integer Gene gut sondern auch für den real Datentyp. Man kann diesen Ersetzungsmechanismus mit anderen Einstellungen des GA kombinieren.

Die Auswirkung von RTR auf ein Problem das zu vorzeitiger Konvergenz tendiert kann man im Test-Programm examples/mgh/testprog.f sehen. Dieses Programm implementiert einige Testfunktionen aus einem alten Benchmark für nichtlineare Optimierungsalgorithmen [5]. Die Testfunktion mit Tendenz zu vorzeitiger Konvergenz ist was die Autoren eine "gausssche Funktion" nennen, beschrieben im Beispiel (9) im Artikel [5] und implementiert als Funktion 3 in examples/mgh/objfcn77.f. Diese Funktion wird angegeben als

\begin{equation*} f(x) = x_1 \cdot e^{\frac{x_2(t_i - x_3)^2}{2}} - y_i \end{equation*}

mit

\begin{equation*} t_i = (8 - i) / 2 \end{equation*}

Und tabellierten Werten für $y_i$ im Artikel bzw. der Implementierung in examples/mgh/objfcn77.f. Das Minimierungsproblem aus diesen Gleichungen ist

\begin{equation*} f (x) = \sum_{i=1}^{m} f_i^2(x) \end{equation*}

mit $m=15$ für dieses Testproblem. Die Autoren [5] geben den Vektor $x_0 = (0.4, 1, 0)$ für das Minimum $f(x_0) = 1.12793 \cdot 10^{-8}$ an. Die Original Fortran-Implementierung in examples/mgh/testprog.f verwendet eine Populationsgröße von 10000 mit den Standardeinstellungen für den real Datentyp von PGAPack. Die große Populationsgröße wird gewählt weil das Problem sonst vorzeitig konvergiert. Nach 100 Generationen wird die Lösung $x_0=(0.3983801, 0.9978369, 0.009918243)$ mit einem Evaluationswert von $f(x_0)=2.849966\cdot 10^{-5}$ gefunden. Die Anzahl der Funktionaufrufe der Zielfunktion waren 105459 (etwas weniger als $10000 + 100 \cdot 1000$ also Bewertung der initialen Population plus etwa 10% der Population in jeder Generation, die Wahrscheinlichkeit für Crossover und Mutation ist nicht 100%, daher kommt es vor dass keiner der Operatoren angewandt wird und die Evaluationsfunktion daher nicht aufgerufen wird)

Ich habe das gleiche Problem mit Differential Evolution [6], [7], [8] in examples/mgh/testprogde.c implementiert (das Treiberprogramm ist in C geschrieben weil ich nicht wirklich Fortran spreche, verwendet aber die gleichen Funktionen wie das Fortran Treiberprogramm und linkt Fortran und C Code in einem gemeinsamen ausführbaren Programm). Es verwendet eine Populationsgröße von 2000 (was für Differential Evolution groß ist wieder aus Gründen von vorzeitiger Konvergenz) und findet die Lösung $x_0=(0.3992372, 0.9990791, -0.0007697923)$ mit einer Bewertungsfunktion von $f(x_0)=7.393123\cdot 10^{-7}$ in nur 30 Generationen. Dies entspricht 62000 Funktionsaufrufen der Zielfunktion (Differential Evolution berechnet eine komplette Generation und entscheidet dann welche Individuen in die neue Generation übernommen werden).

Wenn man RTR für dieses Problem verwendet wie in examples/mgh/testprogdertr.c, kann die Populationsgröße auf 250 reduziert werden und sogar nach 100 Generationen konvergiert die Suche nicht vorzeitig auf eine suboptimale Lösung. Nach 100 Generationen finden wir $x_0=(0.398975, 1.000074, -6.719886 \cdot 10^{-5})$ und $f(x_0)=1.339766\cdot 10^{-8}$ (allerdings mit ein paar anderen Einstellungen für den Differential Evolution Algorithmus). Dies entspricht 25250 Funktionsaufrufen der Zielfunktion.

Restart

Das letzte Mittel wenn die obigen Mechanismen versagen ist ein regelmässiger Restart des GA wann immer die Population zu stark konvergiert. Der Restart-Mechanismus in PGAPack verwendet das beste Individuum der Population um eine neue Population mit durch Mutation erzeugten Varianten dieses Individuums zu initialisieren. Restarts kann man mit PGASetRestartFlag in PGAPack oder dem restart Konstruktor Parameter in PGApy einschalten. Die Frequenz der Restarts (Standard ist alle 50 Generationen) kann mit PGASetRestartFrequencyValue in PGAPack und dem restart_frequency Konstruktor Parameter in PGApy gesetzt werden.