Author D. Selzer-McKenzie
YoutubeVideo:https://youtu.be/nSeqBrvBIjs
Zur Lösung von Optimierungsproblemen und zur Extremwertsuche gibt es neben den an-alytischen Lösungsmethoden auch zahlreiche Algorithmen für die numerische Näherung. Viele Verfahren sind auf gewisse Rahmenbedingungen spezialisiert, woraus natürlich eine starke Limitierung auf bestimmte Probleme aber auch das Erlangen einer gewissen Ef¬fizienz folgt. Beispielsweise kann das Newton-Raphson-Verfahren zur Bestimmung eines Extremums nur auf zweimal stetig differenzierbare Funktionen angewendet werden, stellt dafür aber auch eine schnelle, wenn auch nicht allzu robuste Lösungsmethode dar.
Mit der Zeit wurden immer mehr Verfahren entwickelt, die sich auf gewisse Problemk-lassen spezialisieren. Die relative Anzahl an universell einsetzbaren Methoden nahm in der Folge immer weiter ab. Letztere können auf ein breites Feld von Problemen angewen¬det werden und bieten eine gewisse Robustheit mit entsprechenden Einbußen bei der Effizienz. Doch was spricht für den Einsatz universeller und/oder robuster Verfahren?
Ist die „Gutartigkeit“ eines Problems unbekannt oder ungewiss, verflüchtigen sich die Vorteile spezialisierter Verfahren recht schnell. Verfügt ein Verfahren über eine geringe Resistenz gegenüber Störungen, kann der ursprüngliche Zeitvorteil in einen -nachteil umschlagen - robuste Verfahren bieten dann eine bessere Eignung und stellen damit eine Alternative dar. In der Realität sieht man sich häufig Problemen mit Rauschen, der Suche nach globalen Lösungen oder zeitabhängigen Optima konfrontiert. Für solche Probleme eignen sich universelle Verfahren und im Speziellen genetische Algorithmen.
Genetische Algorithmen gehören zu der Klasse der evolutionären Algorithmen, die zu den universellen Verfahren zählen. Ihr Vorteil ist nicht nur der mögliche Einsatz bei den oben genannten Problemstellungen, sondern auch die prinzipielle Anwendung auf Probleme ohne großes Wissen bezüglich deren Eigenschaften. Ihre Suche basiert auf einer Gruppe von Punkten des Definitionsbereichs, auch Population genannt, wodurch sie im Stande sind lokale Optima zu überwinden und dadurch eine Eignung für multimodale, nichtlineare und diskontinuierliche Optimierungsprobleme erlangen. Zudem bestechen sie durch die recht einfache und damit zeitsparende Implementierung.
Die offensichtlichen Qualifikationen scheinen interessant und vielfältig zu sein, doch wie sehen die Einsatzmöglichkeiten für derartige Algorithmen in der Realität aus und
in welchen Branchen werden sie heutzutage verwendet? Existieren Bereiche, in denen einzig genetische Algorithmen Möglichkeiten zur Lösungsfindung bieten und wie können die Einsatzmöglichkeiten auf weitere Problemstellungen ausgeweitet werden?
Diesem Themengebiet und der Beantwortung der aufgeführten Fragen widmet sich diese Arbeit. Sie befasst sich mit genetischen Algorithmen, deren Theorie und Ein-satzmöglichkeiten, gibt anhand numerischer Tests Ratschläge zur Nutzung und gibt einen Ausblick auf die Anwendungen der genetischen Programmierung, beispielsweise zur Lösung von Systemen partieller Differentialgleichungen.
Die ersten beiden Kapitel führen genetische Algorithmen und deren Theorie ein, gefolgt von einer Untersuchung zur Findung von Richtlinien für Anwendung und Wahl der einzelnen Instrumente und Parameter. Ziel ist eine Hilfestellung für eine Implemen¬tierung zur Lösung unterschiedlicher Probleme. Kapitel 5 und 7 geben einen Überblick über Anwendungsmöglichkeiten in der Wirtschaft anhand expliziter Beispiele.
Zu der Klasse der evolutionären Algorithmen zählt auch die genetische Program¬mierung. Sie erlaubt die Anwendung von biologischen Evolutionsstrategien auf Prob¬lemstellungen unterschiedlicher Arten, beispielsweise zur Lösung von Systemen partieller Differentialgleichungen. Kapitel 6 gibt einen Ausblick über diese Anwendungsmöglichkeiten und deren Nutzen.
2. Einführung in genetische
Algorithmen
Genetische Algorithmen fassen die natürliche Evolution als Optimierungsverfahren auf, bei dem eine Population von Individuen an Umweltbedingungen optimal angepasst wird. Dabei basieren sie auf Mechanismen der biologischen Evolution nach der Theo¬rie des britischen Naturforschers Charles Robert Darwin (1809-1882). Ihren Ursprung haben sie in den Sechziger- und Siebzigerjahren des 20. Jahrhunderts - an der Univer¬sität von Michigan war es John H. Hollands Intention, ein artifizielles Softwaresystem zu entwickeln, welches zum einen dabei hilft den Prozess von natürlichen Systemen zu verstehen und gleichzeitig die Robustheit der Natur beibehält.
Genetische Algorithmen gehören der Klasse der evolutionären Algorithmen an. Dieser Begriff beschreibt Optimierungsverfahren, die von biologischen Paradigmen, wie beispiel-sweise bei Tieren beobachtete Verhaltensweisen oder der Genetik, inspiriert sind. Evo-lutionäre Algorithmen lassen sich nach [19] in vier voneinander unabhängig entwickelte Strömungen unterteilen:
· genetische Algorithmen
· evolutionäre Programmierung
· Evolutionsstrategien
· genetische Programmierung
Von der evolutionären Programmierung wurde erstmals in [7] berichtet. Das Ziel be¬stand darin, eine künstliche Intelligenz durch Simulation der Evolution als einen Lern¬prozess zu erschaffen. Die Evolutionsstrategien gehen zurück auf Ingo Rechenberg und Hans-Paul Schwefel [30]. Ursprünglich ging es um die Entwicklung von Verfahren für die experimentelle Optimierung von verschiedenen physikalischen Problemen, auf die keine gradientenbasierten Methoden angewendet werden konnten. Die Genetischen Pro-grammierung handelt von der Erstellung optimaler Programme zur Berechnung von
gegebenen Eingabe-Ausgabe-Paaren.
2.1. Mechanismen der natürlichen Evolution
Wie bereits erwähnt, formalisieren genetische Algorithmen den Prozess der natürlichen Evolution für die Anwendung auf mathematische Optimierungsprobleme. Aus diesem Grund werden im Rahmen derartiger Verfahren Begriffe aus der biologischen Evolu¬tion und Genetik mit Begriffen aus der mathematischen Optimierung identifiziert. Im Verlauf dieser Arbeit werden dieses Begriffe oftmals als Synonyme verwendet, weshalb eine Gegenüberstellung der Begriffe und ihrer Bedeutungen notwendig erscheint. In der Literatur finden sich Teils widersprüchliche Angaben, weshalb nur die nachstehenden Begriffe für diese Arbeit von Relevanz sind. Sie beziehen sich zur Vereinfachung auf die Binärkodierung; die Bedeutungen lassen sich bei Bedarf auch auf andere Kodierungen übertragen.
biologischer Begriff Bedeutung
Phänotyp Erscheinungsbild; dekodierte Variable
Genotyp Erbbild; kodierte Variable
Population Gruppe von Suchpunkten
Individuum Suchpunkt aus der Population
Chromosom als String kodierter Suchpunkt
Gen Bit
Allel Wert eines Bits
Lokus Bitposition
Fitness Zielfunktionswert
Generation Population zu einem Zeitpunkt
Selektion Auswahl
Rekombination Kreuzungsoperator
Mutation Mutationsoperator
Tabelle 2.1.: Begriffserklärung
Zufällige Variationen im genetischen Code können zu vorteilhaften Eigenschaften führen, welche durch natürliche Auslese an künftige Generationen weitergereicht wer¬den. Individuen einer Population passen sich somit an ihre Umweltbedingungen durch Selektion, Rekombination und Mutation an. Organismen, die an die Umweltbedingun¬gen besser angepasst sind oder über günstige Eigenschaften verfügen, werden begün-
2.2. Aufbau genetischer Algorithmen
stigt, ihre Erbinformation an die Nachkommen weiterzureichen, was den Fortbestand der Art sichert. Individuen, die eine geringere Anpassung aufweisen, sterben mit höherer Wahrscheinlichkeit aus. Umwelteinflüsse und Fehler bei der Rekombination lassen Mu-tationen am Erbgut auftreten, welche zu einer verbesserten Anpassung führen können.
Diese Prozesse wenden genetische Algorithmen in Form von Operatoren auf Werte des Definitionsbereichs, mit dem Ziel einer optimalen Anpassung an das zugrundeliegende Problem, an. Die Anpassung von Generation zu Generation, also Iteration zu Iter-ation, erfolgt durch Änderungen der Erbinformation, der genotypen Darstellung der Suchpunkte. Die genetischen Operatoren lassen sich nicht auf Werte in phänotyper Darstellung anwenden. Das Verfahren betrachtet nicht jeweils nur einen Suchpunkt, sondern ein sogenanntes "Multiset" an Suchpunkten, welche eine Population bilden - in der Biologie: eine Gemeinschaft von Individuen. So werden Punkte aus dem Defin-tionsbereich, welche Teil der Population sind, als Individuen, Chromosomen, oder in der Binärkodierung als Strings bezeichnet. Ein Chromosom/String besteht in der geno-typen Darstellung aus Genen/Bits, welche gewisse Ausprägungen (Allele) haben. In der Binärkodierung besteht folglich ein Chromosom/String aus einer festgelegten Anzahl an Genen/Bits, welche die Werte 0 oder 1 (Allele) annehmen können. Allele stellen somit die einzelnen Elemente des zugrundeliegenden Alphabets dar, welches bei der Kodierung der Suchpunkte definiert wird. Durch die wiederholte Anwendung der genetischen Op¬eratoren auf die Population gelangen genetische Algorithmen zu einer Anpassung der Population an die entsprechenden Umweltbedingungen.
2.2. Aufbau genetischer Algorithmen
Der restliche Teil dieses Kapitels widmet sich dem Aufbau genetischer Algorithmen. Es soll eine klare Vorstellung des Ablaufs vermittelt und die unterschiedlichen, zur Verfü¬gung stehenden Methoden erläutert werden (vgl. [33] und [10]).
Um überhaupt den gesamten Optimierungsprozess starten zu können, bedarf es einer Kodierung der Punkte des Definitionsraumes. Nach der Bestimmung einer Kodierung wird das Verfahren durch die Generierung einer Startpopulation initialisiert. Es folgt eine Ausführung der Prozesse Selektion, Rekombination und Mutation für jedem It-erationsschritt, bis ein Abbruchkriterium erreicht wird. Dieses wird oftmals als eine maximale Anzahl an Iterationen definiert. Ausgegeben wird dann ein Chromosom mit der höchsten Fitness aus der Population des letzten Iterationsschrittes. Bei der Suche nach einem Minimum, kann, wie bei Optimierungsproblemen üblich, die Zielfunktion
2.2. Aufbau genetischer Algorithmen
negiert oder Veränderungen am Algorithmus selbst vorgenommen werden.
Bevor auf die einzelnen Prozesse im Detail eingegangen wird, soll der Ablauf anhand einer Grafik 2.1 verdeutlicht werden:
Abbildung 2.1.: Schematische Darstellung des Ablaufs
Quelle: A. Windisch [1]
2.2.1. Kodierung
Die Kodierung des Definitionsraumes ist problemspezifisch zu wählen. Abhängig von dem zugrundeliegenden Optimierungsproblem, kann eine Binärdarstellung oder beispiel-sweise eine Darstellung von Permutationen optimal sein. Letztere bietet sich beispiel-sweise bei der Suche nach kürzesten Wegen an, bei der jeder Wegpermutation durch die Fitnessfunktion die entsprechende Zeit zugeordnet wird.
Populär ist die Binärkodierung, welche u.a. für die Herleitung des Fundamental¬satzes der genetischen Algorithmen und auch in dieser Arbeit für die numerischen Tests verwendet wird. Mit einer festgelegten Bitlänge können recht einfach natürliche Zahlen dargestellt werden. Eine 5-Bit Kodierung ist in der Lage ganze Zahlen zwischen 0 und 31 darzustellen. Mithilfe von arithmetischen Operatoren können durch die Binärkodierung auch reelle Zahlen umgesetzt werden, was die Anwendungsmöglichkeiten deutlich erhöht. Eine genaue Erläuterung findet sich in Kapitel 5.
2.2.2. Initialisierung
Ein Durchlauf wird Initialisiert mit dem Erzeugen der Startpopulation. Dies geschieht zufallsbasiert und ist abhängig von der zugrundeliegenden Kodierung der Chromosomen.
2.2. Aufbau genetischer Algorithmen
Bei einer Populationsgröße von N müssen zwangsläufig N Chromosomen der Länge L in genotyper Darstellung erzeugt werden, auf die dann in jedem Iterationsschritt die nachfolgenden Operatoren angewendet werden.
2.2.3. Selektion
Jeder Iterationsschritt beginnt mit dem Prozess der Selektion. Aus der resultieren¬den Population des vorangegangenen Iterationsschrittes werden unter gewissen Krite¬rien Chromosomen für den nächsten Iterationsschritt ausgewählt. Die Gesamtanzahl der Chromosomen in der Population wird dabei unverändert gelassen, um die Unversehrtheit des grundlegenden Prinzips genetischer Algorithmen, der populationsbasierten Suche, zu gewährleisten. Allerdings ist es auch möglich die Populationsgröße mit jedem It-erationsschritt zu verkleinern. Welches Vorgehen gewählt wird bleibt dem Anwender überlassen.
Die Selektion dient der Verbesserung der durchschnittlichen Fitness der gesamten Population, durch die Auswahl der Chromosomen mit den höchsten Fitnesswerten. Sie ist somit das wichtigste Instrument für die Suche nach einer geeigneten Näherungslösung.
Es gibt mehrere verschiedene Selektionsvarianten - alle haben jedoch das gemeinsame Ziel der Auswahl der fitnessbezogen besten Chromosomen. Zu den bekannteren Vari¬anten zählen die Roulette-Wheel-Selektion, die Tournament-Selektion und die Truncation-Selektion, welche nun näher erläutert werden sollen.
Roulette-Wheel-Selection
Die Roulette-Wheel-Selection, auch fitnessproportionale Selektion genannt, ist das bekan-nteste Verfahren zur Auswahl geeigneter Individuen. Diese Selektionsmethode war die erste, von John Holland zur Herleitung des Fundamentalsatzes genutzte Variante.
Jedem Individuum, also jedem Chromosom ci mit i = 1,...,N der Population der Größe N, wird eine Auswahlwahrscheinlichkeit zwischen 0 und 1 entsprechend dem Anteil an der Gesamtfitness zugeordnet. Dieser Anteil kann berechnet werden mit
f(ci)
PN (2.1)
j=1 f(cj)
Mit einer höheren Wahrscheinlichkeit steigt auch die Chance für die nächste Gener¬ation, demnach den nächsten Iterationsschritt, ausgewählt zu werden. Nach der Bes¬timmung der Auswahlwahrscheinlichkeiten werden auf dem Intervall [0;1] allen Chro¬mosomen Teilintervalle entsprechend ihren Wahrscheinlichkeiten pci zugewiesen. An-
2.2. Aufbau genetischer Algorithmen
schließend werden gleichverteilt N Zufallszahlen auf dem Intervall [0;1] bestimmt, welche dann insgesamt N Chromosomen für die nächste Generation bestimmen. Hierbei kann ein Chromosom mehrfach (mehrere Kopien des Chromosoms) in die nächste Genera¬tion gelangen. Man kann sich diesen Vorgang mithilfe eines Rouletterads anschaulich erklären. Jedem Chromosom wird entsprechend seiner Auswahlwahrscheinlichkeit eine Fläche auf der Rad zugewiesen. Das Drehen des Rades startet das Zufallsexperiment; kommt das Rad zum Stillstand, wird das Chromosom ausgewählt auf dessen Fläche sich der Zeiger befindet.
Abbildung 2.2.: Veranschaulichung der Roulette-Wheel-Selection
Bei der Roulette-Wheel-Selection besteht durchaus die Möglichkeit, wenn auch mit einer geringen Wahrscheinlichkeit, dass ausschließlich Chromosomen die mit den niedrig¬sten Fitnesswerten ausgewählt werden. Somit ist eine Selektion der fittesten Individuen unter dieser Selektionsmethode nicht garantiert.
Tournament Selection
Die Tournament-Selektion ist die nächste sehr bekannte Selektionsvariante. Wie der Name schon vermuten lässt, werden bei dieser Methode „Turniere“ abgehalten um die Individuen für die nächste Generation zu bestimmen. Für ein „Turnier“ werden min¬destens zwei Chromosomen zufällig aus der aktuellen Population ausgewählt. Das Chro¬mosom mit der höheren Fitness geht als Sieger aus dem „Turnier“ hervor und geht über in die neue Population.
Mehrere Implementierungsmöglichkeiten sind denkbar. Zum einen können bei einer Populationsgröße von N und Auswahl von jeweils zwei Chromosomen insgesamt N Turniere abgehalten werden. Hierbei können dann einige Chromosomen mehrfach in
2.2. Aufbau genetischer Algorithmen
die nächste Generation gelangen. Zum anderen können Turnierteilnehmer aus der Liste für nachfolgende Turniere gelöscht werden, was aber bei einer Wahl von zwei Chro¬mosomen pro „Turnier“ in einer Populationsgröße von insgesamt N/2 Chromosomen für die nächste Generation resultiert. Kann hingegen jedes Chromosom zwei mal für „Turniere“ ausgewählt werden, führt dies zu einer gleichbleibenden Populationsgröße in den nachfolgenden Generationen.
Rank-based Selection
Als weitere nennenswerte Selektionsmethode kann die Truncation-Selection aufgeführt werden. Dem Grunde nach werden bei dieser Methode die stärksten Chromosomen für die nächste Generation ausgewählt. Diese Methode klingt einerseits sehr einfach, ist jedoch relativ zeitaufwendig. Jedem Chromosom wird ein Rang entsprechend seiner Fit¬ness zugewiesen, was eine absteigende Sortierung aller Chromosomen nach ihrer Fitness erforderlich macht (Aufwand N log(N)). Darauf folgt eine Definition der Wahrschein-lichkeitsverteilung der Rangskala; je geringer der Rang, umso höher die Wahrschein¬lichkeit. Die Auswahl der Chromosomen für die nächste Generation erfolgt mithilfe eines Rouletterads anhand der Verteilung. Durch dieses Vorgehen kann das Domi-nanzproblem, bei dem der Fitnesswert die Auswahlwahrscheinlichkeit direkt beeinflusst, vermieden werden.
Elitist Selection
Bei dieser Form der Selektion wird eine gewisse, frei wählbare Anzahl von Chromo¬somen unverändert in die nächste Generation übernommen. Dies bedeutet, dass die Variationsoperatoren nicht auf diese, sondern auf die restlichen Chromosomen angewen¬det werden. Ersichtlich ist, dass dieses Verfahren eine deutliche Zeitersparnis gegenüber anderen Verfahren erwirtschaften kann, da sich der Selektionsprozess dem Grunde nach auf eine Sortierung der Chromosomen nach ihrer Fitness reduziert. Lokale Optima kön¬nen aber schlechter überwunden werden.
Auf die selektierten Chromosomen können nun die beiden Variationsoperatoren Rekom-bination und Mutation angewendet werden. Diese gehören zu den sogenannten Zwei-bzw. Ein-Elter-Operatoren. An dieser Stelle sei darauf hingewiesen, dass auch Mehr-Elter-Operatoren, bei denen z.B. drei Chromosomen miteinander gekreuzt werden, möglich sind.
2.2. Aufbau genetischer Algorithmen
2.2.4. Rekombination
Beim Zwei-Elter-Operator Rekombination (engl. Crossover), findet eine Kreuzung zweier zufällig gewählter Chromsomen aus der Population an mindestens einem zufälligen Punkt statt. Die Rekombinationspunkte befinden sich zwischen den einzelnen Genen, demnach zwischen den einzelnen Bits in der Binärkodierung. Hinter dem Kreuzungspunkt erfolgt ein Tausch der nachfolgenden Abschnitte bzw. Bitfolgen der beiden Chromo¬somen. Dieser Operator ist somit stark abhängig von der Diversität der Population. Befinden sich alle Individuen der Startpopulation im Bereich eines lokalen Optimums, ist die Wahrscheinlichkeit eher gering, dass der Algorithmus diesen Bereich verlässt. Wie bei der Selektion gibt es viele unterschiedliche Rekombinationsvarianten. Im Fol¬genden wird auf 1-Punkt-, 2-Punkt- und uniformen Crossover eingegangen. Zu den weiteren Varianten zählen unter anderem der Cut and Splice Crossover und der Shuffle Crossover.
1-Punkt Crossover
Beim 1-Punkt Crossover, oder auch 1X-Crossover, erfolgt eine Rekombination an einem zufällig gewählten Punkt zwischen den einzelnen Bits. Ein Chromosom der Länge L bietet L − 1 mögliche Kreuzungspunkte. Die Wahl der Schnittpunktes ist zufallsbasiert und wird auf die beiden zufällig gewählten Chromosomen aus der Population angewndet. Ab diesem Punkt erfolgt ein Austausch der Gensequenzen. Auf die Auswahl des Rekom-binationspunktes können unterschiedliche Wahrscheinlichkeitsverteilungen angewendet werden (z.B. die Gleichverteilung oder auch die Normalverteilung).
2-Punkt Crossover
Im Gegensatz zum 1-Punkt Crossover werden beim 2-Punkt Crossover zwei Kreuzungspunkte für die Rekombination gewählt. Der Austausch der Gensequenzen erfolgt zwischen den beiden Kreuzungspunkten.
Uniformer Crossover
Die Wahl der Anzahl der Rekombinationspunkte kann weitergeführt werden bis zum L-Punkt Crossover, bei dem jedes zweite Gen getauscht wird. Daraus entstand die Idee des uniformen Crossover, bei dem mit einer Wahrscheinlichkeit Px für jedes Gen entschieden wird, ob es getauscht wird oder nicht.
2.2. Aufbau genetischer Algorithmen
2.2.5. Mutation
Ergänzend zur Rekombination kann der Variationsoperator Mutation betrachtet werden. Mit diesem Operator können geringfügige Änderungen an den Chromosomen vorgenom¬men werden mit dem Ziel einer stichprobenartigen Erkundung des Suchraumes. Dies verleiht die Möglichkeit, unabhängig von der erzeugten Startpopulation, neue Chromo¬somen erzeugen zu können. Abhängig von der gewählten Kopdierung wird bei diesem Prozess ein Gen mit einer festgelegten, frei wählbaren Wahrscheinlichkeit Pm geändert (mutiert). In der Binärkodierung erfolgt ein einfacher Allel- bzw. Bitwechsel von einer 0 zu einer 1 und umgekehrt. In der Praxis wird für Bitstrings der Länge L eine Muta-tionswahrscheinlichkeit von Pm = 1/L gewählt. Zur Verbesserung der Konvergenz gibt es außerdem die Idee der Verwendung von steigenden Mutationswahrscheinlichkeiten für rechtsliegende Bits. Diese Technik soll ein Phänomen umgehen, bei dem rechtsstehende Bits von der Rekombination unbeeinflusst bleiben.
Anstelle der gewöhnlichen Standardmutation können auch die Allele von zwei Genen innerhalb des Chromosoms vertauscht werden (Zweiertausch). Denkbar wäre auch eine Verschiebung, Permutierung oder Invertierung eines zufällig gewählten Teilstücks eines zufällig gewählten Chromosoms. Letztendlich ist es dem Anwender überlassen, welche Möglichkeit er für die Mutation wählt.
2.2.6. Beispiel Variationsoperatoren
Der Ablauf und das Zusammenspiel der Variationsoperatoren sollen nun an einem ein-fachen Beispiel in der Tabelle 2.2 verdeutlicht werden. Zu sehen sind der 1-Punkt Crossover sowie die Standardmutation.
Population nach Rekombination Population nach Population nach
Selektion Rekombination Mutation
1 1 0 1 1 1 1 0 1 1 1 1 0 1 0 1 1 0 1 0
1 1 0 1 0 1 1 0 1 0 1 1 0 1 1 1 1 1 1 1
1 1 0 0 0 1 1 0 0 0 1 1 0 0 1 1 1 0 0 1
0 1 0 0 1 0 1 0 0 1 0 1 0 0 0 1 1 0 0 0
Tabelle 2.2.: Beispiel Variationsoperatoren
Die Populationsgröße beträgt N = 4 und die Bitlänge L = 5. Somit kann die Popu¬lation in jeder Iteration als eine 4x5-Matrix aus Nullen und Einsen dargestellt werden. Für die Rekombination werden zwei Paare gebildet und jeweils ein Kreuzungspunkt für
2.2. Aufbau genetischer Algorithmen
jedes Paar ausgewählt. Anschließend werden mit der Mutationswahrscheinlichkeit Bits durch die Mutation verändert, womit der Iterationsschritt beendet wird. Diese letzte Matrix repräsentiert die Population, auf die im nächsten Iterationsschritt wieder Selek-tion, Rekombination und Mutation angewendet werden.
3. Das Schema-Theorem
Nachdem nun die Grundlagen für das Verständnis des Aufbaus genetischer Algorith¬men gelegt worden sind, stellt sich die Frage nach der Funktionsweise dieser. Wie gelangt ein derartiger Algorithmus zu einem hinreichend guten Ergebnis? Dieser Thematik nähert sich dieses Kapitel unter Zuhilfenahme von Hollands Schema-Theorem an. Es wurde im Jahre 1975 von John Holland [14] (vgl. auch [33] und [10]) und seinen Studen¬ten hergeleitet und legte das mathematische Fundament für genetische Algorithmen . Das Theorem leitet seine Aussagen aus der Betrachtung der Vervielfältigung von Chro-mosomenschemata über Generationen hinweg her, wobei ein Schema als eine Klasse von Chromosomen betrachtet werden kann. Die Zusammenfassung von Chromosomen zu Schemata steht in Verbindung mit der Notwendigkeit einer populationsübergreifenden Betrachtung, da eine Verfolgung der Entwicklung einzelner Chromosomen wenig sinnvoll erscheint. Vereinfachend werden in der Herleitung ausschließlich fitnessproportionale Se¬lektion, 1-Punkt Crossover und bitweise Mutation betrachtet. Zum besseren Verständnis erfolgt ergänzend eine anschauliche Darstellung der Implikationen des Theorems mittels Hyperwürfeln.
3.1. Herleitung
Für die Herleitung des Theorems sind einige grundlegende Definitionen notwendig, welche im Folgenden aufgeführt werden.
Definition 3.1 (Schema)
Ein Schema H ist ein Bitstring der Länge L über das erweiterte binäre Alphabet {0, 1, *}, dessen Bits nicht ausschließlich mit * besetzt sind H E {0, 1, *}L\{*}L.
Dabei dient * als Platzhalter, auch „Don´t Care“-Symbol genannt, für eine 0 oder eine 1 im Zeichen * E {0, 1}, um einen Satz von Lösungen darstellen zu können. Demnach wird beispielsweise das Chromosom 01101 u.a. durch das Schema 01 * 0* vertreten. Ins¬gesamt existieren 3L−1 denkbare Schemata, da jede Bitposition drei Zustände annehmen
3.1. Herleitung
kann und das suchraumaufspannende Schema {*}L per Defintion nicht als solches gew¬ertet wird.
Definition 3.2 (Ordnung)
Die Ordnung o(H) eines Schemas H gibt die Anzahl der von * verschienen Einträge wider
Somit hat das obige Schema H = 01 * 0* die Ordnung o(H) = 3. Mit der Ein¬führung der Ordnung für Schemata lässt sich folgern, dass ein Schema insgesamt 21−o(H) unterschiedliche Chromosomen darstellen kann. Außerdem gilt:
(L )
2i = 3L − 1
i
denn für die Anzahl von Schemata der Ordnung i gilt:
(L )
2i
i
(L )
Dies ist einleuchtend, da es für ein Schema i-ter Ordnung Möglichkeiten gibt i Bits
i
festzusetzen und 2i Möglichkeiten diesen Werte zuzuweisen. Ergo haben 2L der 3L − 1 Schemata volle Ordnung L und stellen die einzelnen Lösungen des Suchraums dar.
Definition 3.3 (definierende Länge)
Die definierende Länge eines Schemas H ist die Differenz zwischen dem letzten und ersten mit 0 oder 1 besetzten Lokus.
Obiges Schema H = 01 * 0* hat folglich die definierende Länge d(H) = 4 − 1 = 3.
Zur Veranschaulichung lassen sich Schemata mithilfe von Hyperebenen in Hyperwür¬feln darstellen. Betrachtet man einen Bitstring der Länge 3, so enthält der dazugehörige Suchraum 23 = 8 denkbare Lösungen. Dieser Suchraum lässt sich mit einem Würfel visualisieren, wobei jeder Ecke eine Lösung zugewiesen wird. Ein Schema stellt dann eine Hyperebene des Würfels dar; so beschreibt das Schema 0 * * zum Beispiel die Fron¬tebene des Würfels (siehe Abbildung 3.1). In einer derartigen Darstellung von Bits der Länge 3 stellen Schemata voller Ordnung Ecken, Schemata der Ordnung 2 Kanten und Schemata der Ordnung 1 Ebenen des Würfels dar. Es ist ersichtlich, dass jede Lösung, folglich jedes Chromosom, zu 2L − 1 Schemata bzw. Hyperebenen passt, denn jede Lö¬sung gehört, bezogen auf das Beispiel, zu drei Ebenen, drei Kanten und einer Ecke des Würfels (23 − 1 = 7 = 3 + 3 + 1). Oder allgemein:
3.1. Herleitung
!
L = 2L − 1 i
Somit tragen die erzeugten Chromosomen in der Startpopulation Informationen zu einer Vielzahl an Hyperebenen. Es ist plausibel, dass wesentlich mehr Hyperebenen als Chromosomen während des Suchprozesses durchlaufen werden, woraus sich die Hy-pothese des impliziten Parallelismus folgern lässt.
Dank der Verwendung eines binären Alphabets verdoppelt sich die Anzahl der Lösun-gen für jedes zusätzlich kodierte Bit im String, was bei einem Bitstring der Länge 4 zu insgesamt 24 = 16, also doppelt so vielen möglichen Lösungen führt. Der resultierende Suchraum lässt sich ohne großen Aufwand durch zwei ineinander geschachtelte Würfel anschaulich wie folgt darstellen:
Abbildung 3.1.: Anschauliche Darstellung des Suchraumes mit Hyperwürfeln Quelle: Whitley [33]
Es erfolgt eine Erweiterung des dreidimensionalen Falls, indem der äußere Würfel das Präbit 1 und der innere das Präbit 0 erhält. Jeder Würfel stellt weiterhin acht Lösungen des Suchraums dar. Die vier zuvor betrachteten Lösungen des Schemas 0 * * werden nun in der 4-Bit Kodierung durch das Schema 10 * * beschrieben, gekennzeichnet durch die rote Fläche in Abbildung 3.1. Im späteren Verlauf dieses Kapitels soll das Schema-Theorem anhand von Hyperebenen im Hyperwürfel klar machen, wie ein genetischer
3.1. Herleitung
Algorithmus den Suchraum durchläuft.
Nachdem nun die grundlegenden Definitionen und auch Schemata anschaulich dargestellt worden sind, widmet sich der restliche Teil dieses Abschnitts der Herleitung des grundle-genden Theorems, beginnend mit der Formalisierung des Selektionsprozesses. Anschließend wird dieser um Crossover und Mutation erweitert, woraus dann das endgültige Theorem resultiert.
Sei f(ci; t) die Fitness des Chromosoms ci E {0;1}L in der Generation t; weiterhin f(t) = N1 EiN=1 f(ci) die mittlere Fitness der Population und m(ci; t) die Anzahl der Auftritte von ci in der Population. Die Auswahlwahrscheinlichkeit eines Chromosoms ci während der fitnessproportionalen Selektion, auch als Überlebenswahrscheinlichkeit bezeichnet, ist definiert durch:
p(ci) =
f(ci) f(ci)
Enj=1 f(cj) = N ¯f(t)
Es lässt sich nun die erwartete Anzahl der Kopien von ci in der Population nach der
Selektion beschreiben:
E[m(ci, t + Δts)] = f(ci; t)
¯f(t) m(ci; t)
t+Ats beschreibt den Zeitpunkt unmittelbar nach der Selektion; analog beziehen sich
im späteren Verlauf Atc auf den Rekombinations– und Atm auf den Mutationsprozess.
Aus der Gleichung wird klar, dass Chromosomen mit einer überdurchschnittlichen Fit-
ness und/oder höherem m(ci; t) häufiger in der nächsten Generation vertreten sein wer-
den. Dies wird durch folgende Proposition bestätigt:
Proposition 3.4 (Vervielfältigung)
Hat ein Chromosom ci eine überdurchschnittliche Fitness, so vervielfältigt es sich in den
nachfolgenden Generationen annähernd exponentiell
Beweis
Sei f(ci; t) = (1 + 8) f(t) mit 8 > 0. Dann folgt:
E[m(ci, t + Ats)] = f(-ci; t)
) m(ci; t) = (1 + 8)m(ci; t)
Ist die Populationsgröße N ausreichend groß, so ist m(ci; t) relativ klein und 8
3.1. Herleitung
approximativ konstant. Dann gilt:
m(ci; t) = (1 + δ)tm(ci; 0)
II
Die letzte Gleichung impliziert ein annähernd exponentielles Wachstum der Anzahl der Kopien von ci. Analog gilt dies auch für die Sterberate von Chromosomen unter-durchschnittlicher Fitness. Die Ausbreitung der Chromosomen ist jedoch begrenzt, da mit fortlaufen des Algorithmus die durchschnittliche Fitness ¯f(t) durch die Dominanz fitterer Chromosomen in der Population zunimmt. Durch diese Entwicklung sinkt der Selektionsdruck, dem diese unterworfen sind. Darunter ist eine abnehmende Aussterbe-wahrscheinlichkeit für (erfolgreichere) Chromosomen und Schemata zu verstehen.
Sei m(H, t) die Anzahl der Vertreter des Schemas H in der Population c1, c2, ..., cn der Generation t. Die mittlere Fitness des Schemas H für alle ci, die durch H vertreten werden, lautet:
X
¯f(H, t) = 1 m(ci; t)f(ci, t)
m(H, t) ci∈H
Damit lässt sich die hergeleitete Formel (3.1) für die Anzahl der Kopien von Chromo-somen auch auf Schemata anwenden:
f(H; t)
E[m(H, t + zts)] = ¯f(t) m(H; t)
Sie beschreibt die erwartete Anzahl der Vertreter von H in der Population nach der Se-lektion. Chromosomen mit überdurchschnittlicher Fitness haben einen positiven Einfluss auf ¯f(H; t), wodurch eine Präsenz des Schemas in der nächsten Generation wahrschein-
¯f(H;t)
licher ist. ¯f(t) ist der Selektionsdruck, welcher auf dem Schema H lastet.
Mit dem Abschluss der Formalisierung des Selektionsprozesses ist der Grundstein für das Verständnis der Funktionsweise genetischer Algorithmen gelegt worden. Es folgt eine Erweiterung um die beiden Variationsoperatoren: Rekombination und Mutation. Es sei darauf hingewiesen, dass eine Verwendung dieser beiden Prozesse für das Erreichen einer Näherungslösung nicht zwingend erforderlich, aber oftmals ratsam ist. Schließlich können nur durch Variationsmethoden neue Punkte im Suchraum erschlossen werden; eine Fokussierung alleine auf die Selektion führt zu einer starken Abhängigkeit von der zufallserzeugten Startpopulation.
3.1. Herleitung
Für die Rekombination oder auch den Crossover wird zuerst die Auswahlwahrschein-lichkeit pc eingeführt je nach Wahl dieses Wertes müssen nicht alle Chromosomen in diesen Prozess miteinbezogen werden. Schemata der Länge L bieten insgesamt L − 1 mögliche Rekombinationspunkte für 1-Punkt Crossover; eine Rekombination an d(H) dieser Stellen kann das zugrundeliegende Schema stören. In der Folge können dann im Schema festgelegte Genmuster nicht mehr an die Nachkommen weitergereicht werden, wodurch ein aussterben des Schemas möglich ist. Zu beachten ist, dass eine Rekom-bination an den kritischen Stellen zu einer Störung führen kann, es jedoch nicht zwin¬gend muss. Betrachtet man beispielsweise das Schema H = 11 * * * * mit L = 6 und o(H) = 2 − 1 = 1, so führt eine Rekombination von 111010 E H mit 010000 E/ H zwischen den ersten beiden Bits zu keinem Verlust für H, obwohl der Kreuzungspunkt innerhalb der kritischen Zone liegt. Eine Rekombination mit 100110 E/ H stört hingegen das Schema und führt zu einer Auslöschung des Genmusters 11 an den Loki 1 und 2.
Es ist auch möglich, dass eine Rekombination von Chromosomen anderer Schemata zu Nachkommen führt, die dem betrachteten Schema H entsprechen. Somit können im gesamten Crossoverprozess Passungen verloren und hinzugewonnen werden. Die folgende Definition hält dies fest:
Definition 3.5 (Passungen während der Rekombination)
ploss beschreibt die Wahrscheinlichkeit, dass durch Rekombination die Passung eines Chromosoms zum Schema H verloren geht. gains hingegen gibt gewinne an zu H passenden Chromosomen wider
Die erwartete Anzahl der Vertreter von H nach Selektion lässt sich nun um die Rekom-bination erweitern:
E[m(H, t + zts + ztc)] = (1 − pc)E[m(H, t + Lts)]
+ pcE[m(H, t + Lts)](1 − ploss) + gains
Die rechte Seite der Gleichung besteht nun aus drei Summanden. Der erste beschreibt Chromosomen, die nicht von Crossover betroffen sind; diese werden in die nächste Gen-eration direkt übernommen. Der zweite Summand beschreibt Chromosomen, die am Crossover teilnehmen, wodurch eventuell die Passung zum Schema verloren gehen kann. Ergänzt wird diese Berechnung um die Zugewinne an passenden Chromosomen.
Mit der Erkenntnis, dass die Passung verloren gehen kann, wenn ein zu H passendes Chromosom mit einem zu H nicht passenden Chromosom im kritischen Bereich d(H)
L−1
3.1. Herleitung
gekreuzt wird, lässt sich ploss unter der Annahme einer strikten störenden Auswirkung abschätzen:
ploss < d(H) ( 1 L − 1 V )
E[m(H, t + Δts)]
N
Durch das einsetzen dieser Abschätzung und ignorieren von gains kann die Anzahl
der Vertreter weiter ausformuliert werden:
E[m(H, t + Δts + Δtc)] > (1 − pc)E[m(H, t + Δts)]
d(H) (1 E[m(H, t + Ats)] ))
+ pcE[m(H, t + Ats)] (1
L − 1 N
L(H)
( 1 E[m(H t + Ats)] )))
= E[m(H, t + Ats)] (1 − pc + pc (1
L(H) V − 1 N
= E[m(H, t + Ats)] (1 d(H) V (1− N E[m(H, t + Ats)] ))
pU − 1
d
f(H, t) m(H,t) (1 (H) (1
Pc L 1
f(t) NH,t
( ¯f(t) ) m(H' t))) 3.1)
Bevor die Herleitung ihren Abschluss findet, muss noch der letzte Prozess, die Muta-tion, in die Formalisierung integriert werden. Es erfolgt eine bitweise Anwendung der Mutation mit der Mutationswahrscheinlichkeit pm auf die Schemata. Wird ein Bit für die Mutation ausgewählt, verursacht dies einen Wechsel des Allels, auch Bit-Flip genannt. Die Wahrscheinlichkeit, dass dabei ein Schema gestört wird beträgt:
(1 − pm)o(H) (3.2)
1−pm gleicht der Wahrscheinlichkeit eines im Schema festgesetzten Bits, die Mutation zu überstehen; o(H) ist die Ordnung des Schemas H und gibt die Anzahl festgesetzten Bits eines Schemas an.
Mit der Mutation findet die Herleitung des Theorems ihren Abschluss:
Theorem 3.6 (Holland)
Die erwartete Anzahl der Vertreter des Schemas H nach Selektion, Rekombination und Mutation in Periode t kann abgeschätzt werden durch
E[m(H, t + 1)] >
f(H, t) m(H, t) (1 − pc d(H) (1
f(t) L − 1
Pf t) m(H,t))) (1 − pm)o(H) (
3.2. Implikationen
Das Theorem folgt aus der Multiplikation der Mutationswahrscheinlichkeit für Schemata (3.2) mit der rechten Seite der Ungleichung aus der Herleitung des Rekombination-sprozesses (3.1), welche Selektion und Rekombination formal darstellt.
3.2. Implikationen
Aus der Herleitung und dem Theorem selbst ergeben sich einige unmittelbare Implikatio¬nen, welche in diesem Abschnitt erläutert werden sollen. Ersichtlich ist, dass ein Schema in der nachfolgenden Generation häufiger vertreten sein wird, wenn es eine überdurch¬schnittliche Fitness aufweist, von geringer Ordnung ist und eine kurze definierende Länge besitzt. Ein derartiges Schema wird im Selektionsprozess begünstigt und bietet wenige Möglichkeiten für Störungen durch Variationsmethoden, was zu einer Ausbreitung in der Population führt. Zu beachten ist aber auch die Möglichkeit einer überdurchschnittlich hohen Vertretung in der zufallserzeugten Startpopulation.
Schemata mit den eben beschriebenen Eigenschaften können als „Bausteine“ (engl. „building blocks“) bezeichnet werden, da sie durch Rekombination ihre charakteristis¬chen Genkombinationen an ihre Nachkommen weiterreichen und somit zur Eingrenzung des Suchraums beitragen. Je größer die Startpopulation gewählt wird, umso wahrschein¬licher ist das Vorhandensein guter Genkombinationen für den raschen Aufbau guter Lö¬sungen. Bezogen auf das Hyperwürfelmodell stellen „Bausteine“ Hypereben dar, welche eine Vielzahl an Ecken abdecken und die Suche durch Selektion, Rekombination und Mutation auf bestimmte Areale des Hyperwürfels verlagern, vorzugsweise auf diejeni¬gen mit den fitnessbezogen besten Lösungen. Im Verlauf des Algorithmus nimmt die Anzahl der Schemata nach und nach ab und die durchschnittliche Fitness nach und nach zu. Die Suche konzentriert sich dann auf einige bestimmte Areale des Suchraums, bzw. des Hyperwürfels. Ob sich in diesen Arealen die optimale Lösung befindet ist zwar unklar, jedoch ist die Chance bei mehreren separaten Suchorten größer als bei einem einzelnen. Unter Einsatz hoher Mutationswahrscheinlichkeiten und aggressiver Rekom-binationsmethoden, insbesondere dem uniformen Crossover, können zwar recht schnell weite Teile des Suchraums erschlossen werden, bei höherer Iterationsanzahl können sich jedoch diese wegen störender Auswirkungen auf Schemata durchaus kontraproduktiv auswirken. Dieser Umstand verlangt eine sorgfältige Abwägung der Selektionsmeth-ode, der Variationsoperatoren und frei wählbaren Parameter, wie Populationsgröße und Iterationsanzahl, in Abhängigkeit des zu Lösenden Optimierungsproblems.
3.3. Kritik am Schema-Theorem
3.3. Kritik am Schema-Theorem
Einige Kritikpunkte schmälern die Aussagekraft des Schema-Theorems und dessen Imp-likationen ein. Zum einen erfolgt eine Einschränkung durch einige Annahmen zur Vere-infachung - welches Verhalten der Algorithmus bei Wahl anderer Methoden offenbart, wird nicht ersichtlich. Zum anderen erfolgt eine Abschätzung von pjo55 und gains wird gänzlich außer Acht gelassen. In der Praxis ist somit eine Abschätzung der Vertreter eines Schemas in einer bestimmten Generation häufig fehlerhaft, wodurch ein realer Ein¬satz nur in den ersten Generationen sinnvoll erscheint. In höheren Generationen wird die Abschätzung Zusehens schlechter. Eine exakte Berechnung von E[m(H, t + 1)] ist zwar möglich, jedoch relativ aufwändig.
Folglich kann die „Baustein“-Hypothese nur als Erklärung für die ersten Iterationen herangezogen werden, denn mit Fortlaufen des Algorithmus erfolgt eine präzisere Ein-grenzung des Suchraums, wodurch o(H) und d(H) ansteigen. Bei höherer Ordnung und definierender Länge steigt allerdings die Wahrscheinlichkeit für Störungen, die Vari-ationsmethoden an Schemata verursachen können. Die Hypothese ist nur bei hoher Populationsgröße und niedriger Iterationsanzahl sinngemäß.
Die durchschnittliche Fitness eines Schemas H kann sich in wenigen Generationen drastisch ändern, falls eine Konzentration auf einen bestimmten Bereich des Suchraums stattfindet. Somit ist auch die Betrachtung der durchschnittlichen Fitness nur in den er¬sten Iterationsschritten relevant. Danach entwickelt sich eine Abhängigkeit von anderen Schemata, welche vom Theorem gänzlich ignoriert wird. Außerdem sinkt der Selektions-druck mit Zunahme der durchschnittlichen Fitness und Abnahme der absoluten Anzahl der Schemata, resultierend in einer Abnahme der impliziten Parallelisierung. Mit der Möglichkeit der Durchführung einer Skalierung der Fitnesswerte, z.B. nach Goldberg [12] oder der Wahl einer alternativen Selektionsmethode, lässt sich dieser Kritikpunkt allerdings ausmerzen.
Schlussendlich bleibt die Erkenntnis, dass der Erklärungsansatz des Schema-Theorems zum völligen Verständnis leider nicht ausreichend ist. In der Realität scheitert das The¬orem an der Komplexität, der dynamischen Natur und der modularen Struktur genetis¬cher Algorithmen. Für die Wahl der einzelnen verfügbaren Instrumente gibt es zudem keine geschlossene Theorie; hier ist man auf Erfahrungswerte aus der Literatur oder eigene Experimente und Tests angewiesen. Dies gilt ebenso für die Konvergenzdauer, die im Allgemeinen nicht fundiert ist. Trotz allem gelingt mit dem Theorem die Vermit¬tlung einer Idee der Funktionsweise, welche als Hilfreich angesehen werden darf.
4. Implementierung und numerische
Tests
Wie im vorangehenden Kapitel bereits dargestellt, gibt es keine geschlossene Theorie für einen optimalen genetischen Algorithmus, weshalb man für Vergleiche zwischen den einzelnen frei wählbaren Operatoren auf Ergebnisse aus der Literatur oder auf eigene Tests angewiesen ist. In diesem Kapitel wird dem Folge geleistet. Es werden mehrere Tests präsentiert mit dem Ziel einer Hilfestellung für die Wahl der verfügbaren Methoden. Zudem sollen Richtlinien für die Wahl der freien Parameter erarbeitet werden.
Vor der Präsentation der Testergebnisse erfolgt eine kurze Erläuterung des zu testen¬den Verfahrens, eine Darstellung der Implementierung in Pseudo–Code sowie eine der getesteten Selektions- und Rekombinationsmethoden. In den ersten Tests werden an¬hand einer einfachen eindimensionalen Funktion die unterschiedlichen Selektions– und Rekombinationsmethoden im Hinblick auf deren Erfolgsquote und zeitlichen Aufwand miteinander verglichen. Darauf folgend werden weitere Tests für mehrdimensionale Funktionen aufgeführt, die einen Überblick über die Zusammenhänge von Iterationsan-zahl (Abbruchkriterium) und Populationsgröße verschaffen sollen, um nachvollziehen zu können, welche Genauigkeit in welcher Zeit unter welcher Kombination erreicht werden kann.
Abschließend wird auf die mögliche Kombination von genetischen Algorithmen mit anderen Verfahren, im speziellen dem Hill-Climbing Verfahren, eingegangen. Welche Vorteile hat eine Verbindung dieser Art zu bieten?
4.1. Algorithmus und Implementierung 4.1.1. Allgemeiner Aufbau
Es folgt eine kurze Darstellung des grundlegenden Aufbaus genetischer Algorithmen. Zudem wird auf die Frage eingegangen, welche Parameter frei wählbar sind, um den
4.1. Algorithmus und Implementierung
Algorithmus 4.1 : Allgemeiner Ablauf eines GA
Data: Parameter
Result : Chromosom mit maximalem Fitnesswert
1 Initialisierung: Zufallsbasierte Erzeugung der ersten Generation
2 while Abbruchkriterium nicht erfüllt (max. Iterationsanzahl) do
3 Evaluation: Jedem Chromosom wird mit Hilfe der Zielfunktion ein
Fitnesswert zugewiesen
4 Selektion: Auswahl der Chromosomen für die Rekombination
5 Rekombination: Erschaffung neuer Chromosomen durch Kreuzung der durch
die Selektion ausgewählten Chromosomen
6 Mutation: zufällige Veränderung der Nachfahren
7 end
Rahmen für die Tests abstecken zu können.
Daraus ergibt sich die freie Wählbarkeit folgender Parameter:
· Bitlänge (Anzahl der Stellen der Binärzahl)
· Populationsgröße
· Anzahl Iterationen (Abbruchkriterium)
· Methode für Selektion
· Methode für Rekombination
· Methode für Mutation
4.1.2. Bitlänge
Spätestens jetzt stellt sich die Frage, wie man mit Hilfe von Binärzahlen, welche, im Prinzip, Zahlen aus N0 darstellen können, ein reelles Intervall abdecken kann.
Sei B eine Binärzahl, l die Bitlänge und I = [a, b] das gewünschte Intervall. Mit folgender Formel erfolgt eine Transformation auf das Intervall I:
B
2l+1 − 1(b − a) + a
Beispiel: Binärzahl B habe die Bitlänge 10, das gewünschte Intervall sei I = [−5, 5].
Binärzahlen mit einer Bitlänge von 10 können Zahlen von 0 bis 1023 darstellen. Normierung
auf das Intervall [0, 1] : x = B
1023.
4.1. Algorithmus und Implementierung
4.1.3. Initialisierung
Für die Initialisierung müssen die Individuen in der Population zufällig erzeugt werden. Erreicht wird dies mittels einer zufallsgenerierten Matrix M E {0, 1}mxn bestehend aus Nullen und Einsen. Hierbei entspricht in (Anzahl der Zeilen) der Populationsgröße bzw. der Anzahl der Individuen, welche als Binärzahlen gegeben sind. n (Anzahl der Spalten) entspricht der Bitlänge der Binärzahlen. Zur Evaluation, also der Berechnung der Fitness, muss jedesmal die Binärzahl in eine Dezimalzahl umgewandelt werden, bevor sie mithilfe der Zielfunktion ausgewertet werden kann.
4.1.4. Selektionsmethoden
Im Rahmen der numerischen Tests fanden die folgenden drei Selektionsmethoden Beach-tung:
Roullette-Wheel Selection 4.2:
Diese Form der Selektion fand Verwendung während der Herleitung des Schema-Theorems.
Algorithmus 4.2 : Roulette-Wheel Selection
1 Berechne Fitness der Chromosomen fi
2 Berechne Anteil an der Gesamtfitness pi = Pfifi
3 Jedes Chromosom bekommt auf dem "Rouletterad" einen Bereich der Größse pi zugewiesen
4 Berechne in Zufallszahlen aus (0, 1) (Rouletterad drehen) und wähle die entsprechenden Chromosomen aus
Tournament Selection 4.3:
Für die Tournament Selection existieren mehrere unterschiedliche Ansätze. Meistens
werden zufällig 1 k ~ in Chromosomen aus der Population gewählt, von denen
das mit der größten Fitness ausgewählt wird. Da es tendenziell sinnvoll erscheint, einige 'schwächere' Individuen auszuwählen, um im späteren Verlauf weiterhin quali¬tativ hochwertige Nachfahren zu erhalten und trotzdem sicher zu stellen, dass die stärk¬sten Individuen nicht ausselektiert werden, fiel die Wahl auf die Implementierung im Algorithmus 4.3.
Truncation Selection/Elitist Selection 4.4:
Die folgende Methode ist eine Mischung aus der Elitist Selection, bei der Individuen 'ungealtert' mehrere Generationen überleben können, und der Truncation Selection, bei welcher die stärksten ausgewählt werden.
4.1. Algorithmus und Implementierung
Algorithmus 4.3 : Tournament Selection
1 Erstelle Teilnehmerfeld in dem jedes Chromosom zwei mal präsent ist
2 Wähle zufällig zwei Chromosomen aus dem Teilnehmerfeld aus
3 Vergleiche die Fitness
4 Das Individuum mit der höheren Fitness wird ausgewählt
5 Lösche diese beiden Chromosomen aus dem Teilnehmerfeld
6 Fahre fort bis Teilnehmerfeld leer ist
Algorithmus 4.4 : Truncation Selection/Elitist Selection
1 Sortiere Chromosomen nach ihrer Fitness
2 Wähle 50 Prozent der 'stärkeren' für die neue Generation aus
3 Die anderen 50 Prozent entstehen aus Rekombinationen der 'stärkeren' Chromosomen
4.1.5. Rekombinationsmethoden
Es wurden die folgenden drei Rekombinationsvarianten getestet:
· 1-Punkt Crossover
· 2-Punkt Crossover
· Uniformer Crossover
4.1.6. Pseudocode
Zur Übersicht werden nun einige Teile des Verfahrens in Pseudocode 4.5, 4.6, 4.7, 4.8 aufgeführt. Der komplette Quellcode, MATLAB und C++, kann im Anhang eingesehen werden.
Algorithmus 4.5 : Pseudocode genetischer Algorithmus
Data: Zielfunktion, Populationsgröße, Bitlänge und Anzahl Iterationen Result : max. Fitnesswert: maxif(M(i,:))
1 Erzeuge zufällig eine Matrix M E {0, 1}mxn
2 while Abbruchkriterium nicht erfüllt do
3 selection()
4 crossover()
5 mutation()
6 end
4.2. Numerische Tests (eindimensional)
Algorithmus 4.6 : Pseudocode Tournament Selection
1 Erzeuge T = {1, 2, ...,m, 1,2,..., m}
2 fork=1,...,mdo
3 Wähle zufällig i, j T, i =6 j
4 Berechne Fitness von i,j : f(M(i,:)) und f(M(j,:))
5 Wähle stärkeres aus (bspw. i ): Mneu(k,:) = M(i,:)
6 Lösche i,j aus T
7 end
8 Gebe Mneu zurück
Algorithmus 4.7 : Pseudocode Unifromer Crossover
1 ErzeugeT={1,...,m}
2 while T nicht leer do
3 Wähle zufällig i, j T
4 rschaffe 2 neue Chromosomen c1, c2 aus i ( M(i,:) ) und j ( M(j,:) ) durch
uniformen Crossover
5 Übertrage c1, c2 in Mneu
6 Lösche i,j aus T
7 end
8 Gebe Mneu zurück
Algorithmus 4.8 : Pseudocode Mutation
1 Gehe jedes Element von M durch
2 Ändere mit vorgegebener Mutationswahrscheinlichkeit von 0 zu 1 bzw. von 1 zu 0
4.2. Numerische Tests (eindimensional)
Die folgenden Tests verschaffen einen Überblick über die Kombinationen von Itera-tionsanzahl, Populationsgröße, Selektionsmethode, Rekombinationsmethode und Mu¬tation. Da bei diesem Algorithmus die Prozesse zufallsabhängig sind und somit ein ein¬maliges Anwenden des Verfahrens auf das Problem keine objektiven Ergebnisse liefern kann, erfolgt eine 100-malige Ausführung zur Bestimmung einer Erfolgsquote (Anzahl "richtiger Lösungen" auf 6 Nachkommastellen; erster Eintrag) sowie der durchschnit¬tlichen Abweichung von der exakten Lösung (zweiter Eintrag). Die Bitlänge wurde auf n = 10 gesetzt. Die Berechnung dieser ersten Tests erfolgte unter Ubuntu 13.10 64-Bit mit GNU Octave 3.6.4 (Intel Pentium Dual-Core 2x3.2GHz, 4GB RAM). Alle Tests in
4.2. Numerische Tests (eindimensional)
diesem Abschnitt beziehen sich auf nachfolgende eindimensionale Funktion. Die oberen
Zahlen in den Tabellen sind die Erfolgsquoten, die unteren die durchschnittlichen
Abweichungen.
0 0.2 0.4 0.6 0.8 1
x
Abbildung 4.1.: f(x) = sin(irx) auf [0, 1]
Test 1: Populationsgröße 10, Anzahl Iterationen 10, Mutationswahrscheinlichkeit 1/50
1-Punkt 2-Punkt Uniform
Tournament 34% 35% 38%
0.0031938 0.0024602 0.0025656
Roulette Wheel 9% 7% 5%
0.024368 0.027152 0.034031
Elitist 8% 14% 7%
0.011757 0.010879 0.02124
Tabelle 4.1.: Quote und durchschnittliche Abweichung, Test 1
Test 2: Populationsgröße 10, Anzahl Iterationen 50, Mutationswahrscheinlichkeit 1/50
1-Punkt 2-Punkt Uniform
Tournament 62% 62% 58%
0.0035248 0.0018656 0.0028361
Roulette Wheel 32% 25% 29%
0.059202 0.047554 0.051855
Elitist 45% 59% 47%
0.011638 0.011455 0.013761
Tabelle 4.2.: Quote und durchschnittliche Abweichung, Test 2
4.2. Numerische Tests (eindimensional)
Test 3: Populationsgröße 10, Anzahl Iterationen 100, Mutationswahrscheinlichkeit
1/50
1-Punkt 2-Punkt Uniform
Tournament 55% 59% 63%
0.0058105 0.00078751 0.0010315
Roulette Wheel 48% 51% 44%
0.06701 0.061118 0.068841
Elitist 44% 57% 59%
0.023241 0.010784 0.010829
Tabelle 4.3.: Quote und durchschnittliche Abweichung, Test 3
Test 4: Populationsgröße 50, Anzahl Iterationen 10, Mutationswahrscheinlichkeit
1/250
1-Punkt 2-Punkt Uniform
Tournament 91% 77% 74%
3.065e−06 4.8568e−06 1.6173e−05
Roulette Wheel 26% 21% 21%
0.00060505 0.00053886 0.00065658
Elitist 20% 21% 7%
0.00049973 0.00060522 0.00072628
Tabelle 4.4.: Quote und durchschnittliche Abweichung, Test 4
Test 5: Populationsgröße 50, Anzahl Iterationen 50, Mutationswahrscheinlichkeit
1/250
1-Punkt 2-Punkt Uniform
Tournament 99% 95% 96%
1.2732e−06 2.4991e−06 4.8567e−06
Roulette Wheel 45% 35% 42%
0.0017457 0.0045101 0.0026197
Elitist 50% 48% 52%
0.00086335 0.00053838 0.00049904
Tabelle 4.5.: Quote und durchschnittliche Abweichung, Test 5
4.2. Numerische Tests (eindimensional)
Test 6: Populationsgröße 50, Anzahl Iterationen 50, Mutationswahrscheinlichkeit
1/50
1-Punkt 2-Punkt Uniform
Tournament 99% 98% 99%
1.4618e−06 1.3675e−06 1.2732e−06
Roulette Wheel 70% 61% 62%
0.0036952 0.0036904 0.0051865
Elitist 62% 53% 65%
0.00039589 0.00082888 0.00038787
Tabelle 4.6.: Quote und durchschnittliche Abweichung, Test 6
Test 7: Populationsgröße 50, Anzahl Iterationen 50, Mutationswahrscheinlichkeit
1/10
1-Punkt 2-Punkt Uniform
Tournament 100% 99% 91%
1.1788e−06 3.8194e−06 0.00024578
Roulette Wheel 95% 95% 82%
0.0066455 0.0036508 0.0089702
Elitist 86% 79% 86%
0.0045138 0.037691 0.042401
Tabelle 4.7.: Quote und durchschnittliche Abweichung, Test 7 Benötigte Zeit: Populationsgröße: 50, Anzahl Iterationen: 50
Zeit in s 1-Punkt 2-Punkt Uniform
Tournament 5.8844 6.3324 7.7333
Roulette Wheel 6.8087 7.2858 8.7748
Elitist 4.4421 4.4833 4.4698
Tabelle 4.8.: Benötigte Zeit, Test 7
Für den eindimensionalen Fall offenbaren die Tests einerseits große Unterschiede zwis-chen den Selektionsmethoden, andererseits eher marginale Unterschiede zwischen den Rekombinationsmethoden. Es lässt sich eindeutig feststellen, dass die Tournament-Selection die anderen Selektionsmethoden klar dominiert; in allen Tests erzielt sie die besten Resultate. Besonders profitiert sie von einer höheren Populationsgröße und Mu-tationswahrscheinlichkeit. Es können zudem recht schnell gute Ergebnisse in Verbindung mit 1-Punkt Crossover unter niedriger Populationsgröße und Iterationsanzahl erreicht
4.3. Numerische Tests (mehrdimensional)
werden, was zu einer Verringerung des Rechenaufwands genutzt werden kann. Über¬steigt aber die Iterationsanzahl einen gewissen Wert, verschlechtern sich die Ergebnisse, wenn man vom uniformen Crossover absieht. Somit spielen die anderen Selektionver-fahren eine eher untergeordnete Rolle. Die Roulette-Wheel Selection verlangt nach ho¬hen Populationsgrößen, Iterationszahlen und Mutationswahrscheinlichkeiten um vergle¬ichbar zu erscheinen; die Elitist Selection erzielt einzig in Test 2 (2-Punkt Crossover) und 3 (2-Punkt und uniformer Crossover) vergleichsweise gute Ergebnisse, aber mit höheren durchschnittlichen Abweichungen. Zu beachten ist jeodch, dass die Elitist Selection den geringsten Zeitaufwand verursacht.
Bei den Rekombinationsvarianten kann sich unter Tournament Selection der 1-Punkt Crossover behaupten (4 bis 7). Auch unter den anderen Selektionsmethoden schneidet er meistens gut ab und kann den 2-Punkt Crossover in den Tests 4 bis 7 dominieren. Der uniforme Crosover brilliert in Relation zu den anderen Methoden vor allem bei kleineren Populationsgrößen; höhere Mutationswahrscheinlichkeiten dämpfen aber die Resultate. In Test 4 leistet er sich einen Ausreißer unter der Elitist Selection. Die Verbindung von uniformen Crossover und Elitist Selection verlangt anscheinend nach einer höheren Iter-ationsanzahl. Erwähnenswert ist zudem, dass eine höhere Mutationswahrscheinlichkeit häufig gute Ergebnisse bei unveränderter durchschnittlicher Abweichung erzielt.
Aus diesen Ergebnissen kann gefolgert werden, dass eine Erhöhung der Populations-größe deutlich bessere Resultate erzielt. Eine Erhöhung der Iterationsanzahl führt in den meisten Fällen ebenfalls zu besseren Ergebnissen, allerdings nicht im gleichen Ausmaße. Da beide Möglichkeiten ungefähr den gleichen Aufwand verursachen, kann gefolgert wer¬den, dass man zur Verbesserung der Ergebnisse tendenziell eher die Populationsgröße anstelle der Iterationsanzahl erhöhen sollte. In einigen Fällen vollzieht sich nach einer gewissen Anzahl an Iterationen eine Verschlechterung der Resultate. Demnach besteht die Notwendigkeit einer dem Problem angemessenen Wahl der Iterationszahl - mehr Iterationsschritte führen nicht zwangsläufig zu einer Verbesserung.
4.3. Numerische Tests (mehrdimensional)
Nun wird der genetische Algorithmus auf mehrdimensionale Funktionen f : D C R2 → R angewendet. Die Bitlänge wird auf n = 30 festgesetzt, als Selektionsvariante Tourna-ment Selection, als Rekombinationsmethode der uniforme Crossover und eine Mutation-swahrscheinlichkeit von 1
10n gewählt. Die folgenden Tests wurden mit GNU C++ unter Ubuntu 13.10 bei unveränderter Hardware kompiliert. Auch hier wird der Algorithmus
4.3. Numerische Tests (mehrdimensional)
jeweils 100 mal durchgeführt. Die nachstehenden Tabellen zeigen die durchschnittliche Abweichung (oben) von der exakten Lösung sowie die benötigte Zeit (unten), um einen Vergleich zwischen Populationsgröße und Iterationsanzahl zu geben. Beschriftungen in der linken Spalte stellen die Populationsgröße und Beschriftungen in der obersten Zeile die Iterationsanzahl dar.
Test 8:
Abbildung 4.2.: f(x, y) = 50 − x2 − y2 auf [−5,5] x [−5,5]
Population / Iterationen 10 50 100 500
10 0.85861933 0.19815222 0.20910056 1.1546998
0.0008s 0.0039s 0.0077s 0.0387s
Tabelle 4.9.: Durchschnittliche Abweichung und die benötigte Zeit, Test 8
4.3. Numerische Tests (mehrdimensional)
Test 9:
0.5
0
2
0
−1 0 1 2 −2
x
Quelle: eigene Darstellung
Abbildung 4.3.: f(x, y) = xexp(−x2 − y2) auf [−2,2] x [−2,2]
Population / Iterationen 10 50 100 500
10 0.00463929 0.03499361 0.03509808 0.04269244
0.0008s 0.0039s 0.0079s 0.0041s
50 0.00337451 0.00079542 0.00051921 0.00038902
0.0035s 0.0178s 0.0352s 0.1743s
100 0.00017032 0.00045263 0.00029578 0.00028805
0.0071s 0.0349s 0.0707s 0.3425s
500 0.00041027 3.76329e − 05 1.41898e − 09 4.63462e − 05
0.0353s 0.1837s 0.3445s 1.753s
Tabelle 4.10.: Durchschnittliche Abweichung und die benötigte Zeit, Test 9 Test 10:
Abbildung 4.4.: f(x, y) = sin(2irx) + y auf [0, 1] x [0, 1]
4.4. Verbindung mit dem Hill-Climbing Verfahren
Population / Iterationen 10 50 100 500
10 0.10239198 0.09229848 0.20369982 0.04902680
0.0008s 0.0038s 0.0078s 0.039s
50 0.00632808 0.00030542 2.75457e − 06 2.0825e − 05
0.0034s 0.0176s 0.0346s 0.1732s
100 0.00197074 4.30026e − 08 3.05991e − 07 9.02246e − 06
0.007s 0.034s 0.0688s 0.3986s
500 0.00070039 1.1235e − 11 0 0
0.035s 0.1746s 0.3446s 1.719s
Tabelle 4.11.: Durchschnittliche Abweichung und die benötigte Zeit, Test 10
Die hier offerierten Tests zeigen, dass die optimale Populationsgröße sowie die Anzahl der Iterationen von der Zielfunktion abhängen. Bei allen Tests ist eine gewisse Kombi¬nation beider Größen für eine bestimmte Genauigkeit erforderlich. Für eine Genauigkeit auf 5 Stellen sind bei Test 8 mindestens 50 Iterationen und eine Populationsgröße von 100 vonnöten, bei Test 9 hingegen eine Populationsgröße von 500. Test 10 verlangt für die gleiche Genauigkeit moderate 100 Iterationen und eine Populationsgröße von 50.
Allgemeine Richtlinien lassen sich infolgedessen nur schwer angeben; feststellen kann man aber, dass zur Verbesserung der Lösung wieder die Populationsgröße anstelle der Iterationsanzahl angehoben werden sollte. Eine „zu hohe“ Anzahl an Iterationen führt auch hier zu einer Zunahme des Fehlers.
4.4. Verbindung mit dem Hill-Climbing Verfahren
In diesem Abnschitt soll darauf eingegangen werden, wie sich der genetische Algorith¬mus sinnvoll mit anderen Verfahren kombinieren lässt und welche Vorteile dies mit sich bringt. Mit genetischen Algorithmen kann relativ schnell und einfach mittels kleiner Population und Iterationsanzahl eine gute Approximation der gesuchten Lösung bes¬timmt werden. Um eine hohe Genauigkeit zu erreichen wird allerdings eine deutlich höhere Populationsgröße und/oder Anzahl an Iterationen benötigt.
Aufgrund dieser Erkenntnis kann es als sinnvoll erachtet werden, eine Näherungslö-sung mit einer kleinen Population und geringer Iterationsanzahl zu bestimmen, welche als Startwert für ein anderes Verfahren, wie beispielsweise das Gradientenverfahren dient. Da aber das Gradientenverfahren, anders als die genetischen Algorithmen, nur auf dif-ferentierbare und stetige Funktionen anwendbar ist, fiel die Wahl auf eine Kombination mit dem Hill-Climbing Verfahren. Damit ist weiterhin ein breiter Anwendungsbereich sichergestellt.
4.4. Verbindung mit dem Hill-Climbing Verfahren
Um die Funktionsweise klar zu machen, folgt eine kurze Darstellung des Hill-Climbing Verfahrens für den eindimensionalen Fall:
Algorithmus 4.9 : Ablauf des Hill-Climbing Verfahrens
Data : Startwert x0, Schrittweite h
1 Berechne v1 = f(x0 + h), v2 = f(x0 − h)
2 if x0 > v1 undx0 > v2 then
3 x0 ist die Lösung
4 else if x0 < v1 und v1 > v2 then
5 x1 = v1
6 else if x0 < v2 und v2 > v1 then
7 x1 = v2
8 etc
Dieses Verfahren lässt sich zwar auch problemlos auf mehrdimensionale Funktionen anwenden, wodurch jedoch der Aufwand deutlich ansteigt. Dies hängt damit zusammen, dass wesentlich mehr Richtungen 'abgetastet' werden müssen. Es folgen zwei Tests zum kombinierten Verfahren:
Test 11: h = 10−4
Population / Iterationen 10 50 100 500
10 4.50583e − 05 8.35800e − 05 6.11396e − 05 4.43963e − 05
0.0015s 0.0039s 0.0114s 0.0397s
50 3.72960e − 05 2.63988e − 05 2.38828e − 05 3.31514e − 05
0.0036s 0.0172s 0.0372s 0.1723s
100 6.31381e − 05 1.52321e − 06 5.9032e − 09 9.5042e − 06
0.007s 0.037s 0.069s 0.3412s
Tabelle 4.12.: Durchschnittliche Abweichung und die benötigte Zeit, Test 11 Test 12: h = 10−6
Population / Iterationen 10 50 100
10 5.88087e − 07 6.3011e − 07 8.0705e − 07
0.0442s 0.0426s 0.0314s
50 2.76183e − 07 4.64412e − 07 1.79602e − 07
0.0078s 0.017s 0.0364s
Tabelle 4.13.: Durchschnittliche Abweichung und die benötigte Zeit, Test 12
Bei Betrachtung der Testergebnisse fällt sofort auf, dass das kombinierte Verfahren
schon bei Wahl niedriger Parameterwerte hohe Genauigkeiten erzielt. Die Wahl niedrigerer
4.5. Zusammenfassung
Werte resultiert in einer kürzeren Prozessdauer. Hierbei spielt zusätzlich auch die Wahl der Schrittweite eine entscheidende Rolle. Eine kleinere Schrittweite führt zu einer höheren Genauigkeit, aber auch zu höherem Rechenaufwand.
Im Vergleich mit den Ergebnissen des vorigen Abschnitts kann die Kombination von genetischen Algorithmen mit dem Hill-Climbing Verfahren die Performance teil¬weise verbessern. Tendenziell ist es sinnvoll, eine sehr kleine Populationsgröße und Iterationsanzahl zu verwerden, um eine Maximierung der Leistungsfährigkeit zu erre¬ichen. Höhere Genauigkeiten verlangen höhere Parameterwerte, da ansonsten dass Hill-Climbing Verfahren zu viele Iterationsschritte benötigt. Bei kluger Wahl aller Parame-terwerte kann dieses Verfahren eine Zeitersparnis gegenüber dem gewöhnlichen genetis¬chen Algorithmus erzielen.
4.5. Zusammenfassung
Die Tests dieses Kapitels geben einige Einblicke in die praktische Anwendung von genetis-chen Algorithmen. Sie zeigen deutlich, dass die Wahl der einzelnen Methoden und Pa-rameter von dem zugrundeliegenden Problem abhängig ist. Die markanteste Erkenntnis dieses Kapitels ist die übergeordnete Bedeutung der Selektionsmethode und der Pop-ulationsgröße für die Erreichung der gesteckten Näherung. Unter Wahl einer großen Population und der Tournament Selection kann recht schnell ein akzeptables Ergebnis erzielt werden - die weiteren Parameter spielen dann eine eher untergeordnete Rolle. Bei den Rekombinationsvarianten sind die Unterschiede nur moderat, woraus sich nicht fol¬gern lässt, welche Methode bevorzugt werden sollte. In Kombination mit der Tounament Selection erzielt jedoch der 1-Punkt Crossover die besten Ergebnisse im gesamten Test-durchlauf. Entgegen der Theorie aus Kapitel 3 trägt eine höhere Mutationswahrschein-lichkeit in den meisten Fällen zu einer Verbesserung die Lösungsquote bei. Auch der uniforme Crossover scheint bei hohen Iterationszahlen nicht so destruktiv zu sein wie prognostiziert.
Unter der Kombination des genetischen Algorithmus mit dem Hill-Climbing Verfahren kann eine Zeitersparnis erreicht werden, was jedoch eine smarte Wahl der Parameter und der Schrittweite voraussetzt.
Für weiterführende Aussagen und zur Festigung der hier erzielten Resultate sind weit¬ere Tests für eine Vielzahl an Funktionen und Problemstellungen erforderlich. Eine de¬rartige Untersuchung würde aber den Rahmen dieser Arbeit sprengen und so bleibt an dieser Stelle nur der Verweis auf die Notwendigkeit weiterer Untersuchungen.
5. Das Problem Economic Dispatch
Der Anwendungsbereich der genetischen Algorithmen in der Wirtschaft ist sehr groß. In diesem Kapitel wird eine sehr wichtiges Optimierungsproblem, das bei allen En-ergieunternehmen auftritt, mithilfe der GA gelöst. Das Problem Economic Dispatch beschreibt Bestimmung der Menge an Strom, die an mehreren Anlagen zu möglichst geringen Kosten produziert werden muss, um die Nachfrage der Verbraucher zu decken. Nach R. Ouiddir, M. Rahli und L. Abdelhakem-Koridak [25] muss die Balance zwischen dem Stromerzeugnis und der variablen Auslastung des Stromnetzes (bedingt durch die variable Nachfrage) gefunden werden. Dabei ist wichtig, welche Anlage wie viel Strom produziert. Diese Entscheidung muss alle wenige Stunden bis alle wenige Minuten getrof¬fen werden.
5.1. Formulierung
Das Problem Economic Dispatch wird nach [25] und P. Chen, H-C. Chang [27] wie folgt formuliert:
min F = Xn i=1 fi (Pi)
s.t. Xn Pi = PD +Ploss, (ED)
i=1
P
min i < Pi < Pmax i .
Dabei bezeichnen F die gesamten Kosten der Stromerzeugung, Pi die Erzeugung an
der Anlage i, Pmin
i und Pmax
i die minimale bzw. die maximale Erzeugung der Anlage i, n die Anzahl der Anlagen, fi (Pi) die Produktionskosten der Anlage i, PD die Auslastung des Systems und Ploss die Übertragungsverluste.
Die Produktionskosten fi (Pi) werden als quadratische Funktion dargestellt:
fi(Pi) = aiPi + biPi+ ci, (5.1)
5.1. Formulierung
wobei ai, bi und ci Konstanten sind.
Die Übertragungsverluste Ploss haben geographische und/oder physikalisch-technische Gründe und können z.B. als abgegebene Wärme aufgefasst werden. Ihre Bestimmung kann auf zwei Wegen erfolgen:
· mit den Penalty-Faktoren (dazu siehe [28]),
· mit den B-Koeffizienten.
Die letzte Methode wird hier verwendet. L. K. Kirchmayer, G. W. Stagg [21], A. F. Glimn, L. K. Kirchmayer, J. J. Skiles [11] und E. E. George [9] unterstellen, dass die Verluste auf gegenseitige Eigenschaften der Stromlinien (Entfernungen zwischen Kraftwerken und Umspannwerken, Material der Stromleitungen) und der Sammelschienen der Umspannwerke zurückzuführen sind. Deswegen werden Übertragungsverluste als quadratische Funktion in den Produktionsmengen betrachtet:
Die Koeffizienten Bij beinhalten mehrere charakteristische Größen der Stromnetze: Spannungen, Widerstände der Linien, Entfernungen usw.
Die Betriebsbereiche der einzelnen Kraftwerke können durch kraftwerkspezifische Gren-zen limitiert werden. Ebenso spezifisch ist der Ausmaß der Anpassung (ramp rate limits) der zu produzierenden Menge an Strom begrenzt. Wenn am Kraftwerk i heute Pi0 MWh erzeugt werden, und morgen soll die Erzeugung auf Pi erhöht werden, dann kann man sie maximal um URi erhöhen. Analog, wenn morgen weniger produziert werden soll, kann man die Produktion höchstens um DRi herabsenken. Also:
· Pi− Pi0 ≤URi, falls die Produktion erhöht werden soll,
· Pi0− Pi
Die Betriebsbereiche der Kraftwerke können auch durch s.g. prohibited zones ver-
botene Zonen eingeschränkt werden. Diese Zonen sind bestimmte Leistungen, die ein Kraftwerk aufgrund seiner bautechnischen Besonderheiten nicht erzeugen kann. Diese können z.B. Vibrationen in den Maschinen sein. Die Produktion in den verbotenen Zo¬nen kann zum Ausfall einzelner Anlagen oder sogar ganzer Kraftwerke führen oder der Verlauf der Input-Output-Kurven ist in diesen Zonen unbekannt, nicht eindeutig oder zu komplex. Nach [27] ist die optimale produzierte Leistung diejenige, die außerhalb
5.2. Implementierung
dieser Zonen liegt. Um diese Zonen umzugehen, werden Heuristiken angewendet. Sollte Pi
min F = Xn i=1 fi (Pi)
s.t. Xn Pi=PD+Ploss, (ED*)
i=1
( ) ( )
max P min
i , P i 0 − DRi < Pi < min P max
i , P i 0 + URi ,
Somit ergibt sich, dass die Lösungsmenge des Problems aus mehreren Intervallen in R besteht (roter Bereich in der Abbildung 5.2) und damit nicht konvex ist. Die konvexen Optimierungsverfahren kommen damit für die Lösung des Problems nicht infrage. Es bietet sich ein genetischer Algorithmus an.
5.2. Implementierung
Die einzige Variable, die kodiert wird, ist λnm mit 0 < λnm < 1, das normierte Zusatzkosten beschreibt, die in diesem Fall mit Grenzkosten der Produktion übereinstimmen1. Damit ist das Problem unabhängig von der Anzahl der Kraftwerke, was die Lösungsmethode auch für große Stromnetze geeignet macht. Es wird behauptet, dass durch diese Wahl
t t + 1
5.2. Implementierung
Ppz− Ppz +
der Kodierung die Lösungszeit linear mit der Anzahl der Kraftwerke steigt.
Die Autoren benutzen 10 Bits für die Kodierung. Sie behaupten, dass die Genauigkeit der Lösung mit der Anzahl der Bits zunimmt, wobei gleichzeitig die Konvergenz langsamer wird. Die Länge von 10 Bits wird hier wie folgt benutzt:
d1 d2 d3 d4 d5 d6 d7 d8 d9 d10
x x x x x x x x x x
2−1 2−2 2−3 2−4 2−5 2−6 2−7 2−8 2−9 2−10
wobei di E {0,1}, i = 1,2, ...10. Die Dekodierung erfolgt gemäß
λnm = 10
X
i=1 di2−i, di E {0, 1}. (5.3)
Die Beziehung zwischen den tatsächlichen und normierten Zusatzkosten ist gegeben
durch
Aact = Arsny'isn + Anm (ArZax − AZ") , (5.4)
Mit der Annahme, dass die Lösungsmenge des Problems konvex ist, kann man (ED*)
5.2. Implementierung
mit dem Lagrangeansatz lösen:
Pfi (2aiPi + bi) = Aact, max (Prin, Pi0 − DRi < Pi < min (Pimax, Pi0 − URi
Pfi (2aiPi + bi) < Aact, Pi min (Pimax, Pi0 URi (*)
Pfi (2aiPi + bi) > = max (Pmini, − DRi ,
wobei Pfi = 1
∂Ploss
1+
Es wird ∂Pi
Pi − PD + Ploss ~~~~~ (5.5)
definiert. Das Abbruchkriterium ist erfüllt, wenn ε unter einem bestimmten Wert liegt. Die normierte Fitnessfunktion wird wie folgt definiert:
1
F = mit Normierungsvariable k = 200. (5.6)
1 + k E
(PD
Die Selektion wird mit der Regel Roulette wheel simuliert. Folgende Werte wer¬den gewählt: Populationsgröße = 16, Crossoverwahrscheinlichkeit = 0.8, Mutation-swahrscheinlichkeit = 0.1.
Der Ablauf ist im Algorithmus 5.1 beschrieben. Schritt 2 von 5.1 ist im Algorithmus 5.2 beschrieben.
Algorithmus 5.1 : Ablauf des GA
Data : Anlagen, Nachfrage ... Result : Output
1 Initialisiere Startpopulation
2 Evaluiere jedes Chromosom
3 Ordne alle Chromosomen nach ihrer Fitness
4 Wähle beste Eltern für die Reproduktion
5 Wende Crossover und Mutation an
6 Evaluiere die neue Generation und ersetze die schwächsten durch die stärksten
7 if Konvergenzkriterium (5.5) erfüllt then
8 Gebe die Population aus
9 else
10 Gehe zu 4
5.3. Zahlenbeispiel
Algorithmus 5.2 : Auswertung der Individuen
1 Bestimme )nm mit (5.3)
2 Bestimme )act mit (5.4)
3 i = 1
4 Berechne Pi mit (*)
5 if Pi in einer verbotenen Zone then
6 if Phase der steigenden Auslastung then
7 Pi=Ppz+
8 else
9 Pi=Ppz−
10 if i die letzte Anlage then
11 Berechne Ploss mit (5.2)
12 else
13 i= i+ 1 und gehe zu 4
14 Berechne die Fitness mit (5.5) und (5.6)
5.3. Zahlenbeispiel
Die Autoren von [27] konstruieren ein Stromnetz aus 3 Kraftwerken mit Eigenschaften, die in Tabellen 5.1 und 5.2 zusammengefasst sind.
/ /
Anlage P min P max ai in $ MW 2 bi in $ MW 2 ci in $
i i
1 50 250 0.00525 8.663 328.13
2 5 150 0.00609 10.04 136.91
3 15 100 0.00592 9.76 59.16
Tabelle 5.1.: Koeffizienten und Leistungsgrenzen
Anlage P 0
i URi MW h DRi MW h
( / ) ( / ) Prohibited zones (MW)
1 215 55 95 [105,117],[165,177]
2 72 55 78 [50,60],[92,102]
3 98 45 64 [25,32],[60,67]
Tabelle 5.2.: Ramp rates und prohibited zones
5.3. Zahlenbeispiel
Die Matrix Bi3 gibt die Verlustkoeffizienten an:
Bi,3 = ⎛ ⎞
0.000136 0.0000175 0.000184
⎜ ⎟
⎜ ⎟
⎜ 0.0000175 0.000154 0.000283 ⎟
⎝ ⎠
0.000184 0.000283 0.00161 .
Die Nachfrage beträgt 300 MW, als Konvergenzkriterium wird < 0.3MW(0.001 * Nachfrage) gewählt.
Die Berechnungen wurden am Rechner mit dem 486-33 Prozessor durchgeführt. In der zweiten Generation ist das Konvergenzkriterium erfüllt und die optimale Lösung ist gefunden (nur die ersten drei Individuen werden aufgelistet):
Rank Chromosom )nm ε(MW) Fitness
1 1000111010 0.556641 3.4625 0.8125
2 1001010101 0.583008 6.5155 0.6972
3 1000101011 0.541992 9.7654 0.6057
Tabelle 5.3.: Initialisation
Rank Chromosom )nm ε(MW) Fitness
1 1001000000 0.562500 1.2371 0.9241
2 1001001000 0.570313 1.7308 0.8966
3 1000111010 0.556641 3.4625 0.8125
Tabelle 5.4.: 1. Iteration
Rank Chromosom )nm ε(MW) Fitness
1 1001000011 0.565430 0.1192 0.9921
2 1001000010 0.564453 0.4898 0.9684
3 1001000000 0.562500 1.2317 0.9241
Tabelle 5.5.: 2. Iteration
Hier sind die Ergebnisse2:
2Rundungsfehler möglich
5.4. Anwendung auf das westalgerische Stromnetz
· P1 = 194.265 MW, P2 = 49.999 MW, P3 = 79.627 MW
· Ploss = 24.011 MW
· λact = 10.7028 $/MWh
· Zeit = 0.016s.
5.4. Anwendung auf das westalgerische Stromnetz
Die Autoren von [25] wenden diesen Algorithmus auf das westalgerische Stromnetz mit zwei Kraftwerken:
Abbildung 5.3.: Stromnetz von Westalgerien Die Fitnessfunktion ist gegeben durch
1
F=
1
1 +
>i fi
Die Parameter werden wie folgt gewählt: Populationsgöße = 10, Crossoverwahrschein-lichkeit = 0.88, Mutationswahrscheinlichkeit = 0.1. Die Kostenfunktionen lauten:
f1 (P1) = 0.85P12+ 150P1 + 2000,
f2 (P2) = 1.7P22+ 250P2 + 3000,
5.4. Anwendung auf das westalgerische Stromnetz
mit
30 < P1 < 510 MW,
10 < P2 < 70 MW,
PD = 505 MW,
Fall 1: Ploss = 15.94,
Fall 2: Ploss = 0.0189P1 + 0.0924P2.
Zur Vereinfachung wird hier unterstellt, dass es keine verbotenen Zonen gibt. Im ersten Fall werden die Übertragungsverluste als konstant angenommen mit Ploss = 15.94MW. Die Ergebnisse sind in der Tabelle 5.6 zusammengefasst. Im zweiten Fall werden die Ver-luste als lineare Funktion in Pi, i = 1, 2 angenommen, mit Ploss = 0.0189P1 + 0.0924P2. Die Ergebnisse sind in der Tabelle 5.7 zu sehen. Die Berechnungen werden mit konvexen Lösungsmethoden und den tatsächlichen Entscheidungen des Netzbetreibers verglichen.
GA Fletcher-Reevesa Fletcherb Sonelgazc
P1 450.95 466.64 469.93 465.94
P2 69 54.25 49.98 55
Ploss 15.94 15.94 15.94 15.94
Kosten (Nm3/h) 241764 278649 279940 278319
Berechnungszeit 0.05 0.01 0.01 /
Generation 304 8 3 /
Tabelle 5.6.: Fall 1
Im ersten Fall mit konstanten Verlusten liefert der genetische Algorithmus im Vergleich zu konvexen Methoden und zu Sonelgaz bei vergleichbarer Laufzeit deutlich niedrigere Kosten mit Jahresersparnis von über 2%. Im zweiten Fall können dank der linearen Form die Verluste deutlich gesenkt werden, obwohl sie höher sind als bei den konvexen Methoden. Die Kosten sind dagegen im GA niedriger und die Jahresersparnis beläuft sich auf ca. 3%. Somit hat der genetische Algorithmus den entscheidenden Vorteil gegenüber den konvexen Optimierungsmethoden.
Als Alternative zu üblichen innerbetrieblichen Lösungen des Problems Economic Dis-patch kann die Verwendung der auf genetischem Algorithmus basierten Methode bessere
aNichtlineare konjugierte Gradienten Methode bBFGS-Verfahren bzw. DFP-Formel
cDaten des algerischen Netzbetreibers Sonelgaz
5.4. Anwendung auf das westalgerische Stromnetz
GA Fletcher-Reeves Fletcher Sonelgaz
P1 449.98 465.374 468.92 465.94
P2 70 53.36 49.5 55
Ploss 14.97 13.72 13.43 15.94
Kosten (Nm3/h) 270444 277067 278779 278319
Berechnungszeit 0.1 0.05 0.01 /
Generation 76 25 3 /
Tabelle 5.7.: Fall 2
Ergebnisse liefern. P. Chen, H-C. Chang [27] zeigen, dass die vorgestellte GA-basierte Methode insbesondere für große Stromnetze bessere Ergebnisse liefert und schneller ist, als die weit verbreitete Lambda-Iterationsmethode. Sie behaupten außerdem, dass die letzte solche Besonderheiten wie ramp-rate limits, verbotene Zonen und Übertra-gungsverluste nicht berücksichtigen kann.
Darüber hinaus ist die Prognose der Energienachfrage mithilfe genetischer Algorith-men möglich (siehe z.B. H. Ceylan und H. K. Ozturk [5]). Es ist damit möglich, beide Ansätze zur Optimierung der Energieproduktion zu verwenden.
6. Weitere Ansätze und
Anwendungen auf die numerische
Lösung von
Differentialgleichungen
6.1. Grammatical Evolution
Für die folgende Theorie zur numerischen Lösung von Differentialgleichungen wird nun eine weitere Verfeinerung der Evolutionary Computing–Familie eingeführt: den Ansatz der sogenannten Grammatical Evolution im Folgenden als GE bezeichnet.
Diese weitere Verfeinerung liegt in der formulierten Zielsetzung der Lösung von Differ-entialgleichungen bzw. dynamischen Systemen begründet: Zwar ist es bisher möglich, Iterationslösungen von GAs sinnvoll mit Hilfe eines Alphabets zu kodieren, für die Darstellung von Differentialgleichungen ist es aber darüber hinaus erforderlich, Funk¬tionen und ihre Ableitungen darzustellen. Eine hilfreiche Lösung dieses Problems bi¬etet das Instrument der GE, erstmals vorgestellt von Koza [22]); das im Folgenden vorgestellte Konzept basiert auf einer Arbeit der Autoren Tsoulos und Lagaris [32]. Der Grundbaustein von GE ist eine kontextfreie Grammatik in Backus–Naur–Form, d.h. ein Quadrupel {N, T, P, S}, wobei N eine Menge von Nichtterminalsymbolen, T eine Menge von Terminalsymbolen, P eine Menge von Produktions– bzw. Ableitungsregeln sowie S E N ein Startsymbol bezeichne. Die Individuen sind nun wie bei GAs als Vektoren von Integers kodiert, wobei jedes Integer für eine Ableitungsregel stehe. Die Auswahl der Ableitungsregel erfolgt in zwei Schritten: zunächst wird von links nach rechts des ungeordneten Chromosomenvektors ein Integer nach dem anderen eingelesen und ihm
6.1. Grammatical Evolution
sein tatsächlicher Wert V zugewiesen. Im Anschluss erfolgt die Auswahl nach folgender Rechnung:
Rule = V mod NR, (6.1)
wobei NR die Anzahl der möglichen Ableitungsregeln je Nichtterminale bezeichne. Dieses technische Vorgehen soll nun anhand eines Beispiels verdeutlicht werden.
Beispiel 6.1
Sei der Vektor j = [16, 3, 7, 4, 10, 28, 24, 1, 2, 4] gegeben.
Behauptung 6.2
j kodiert die Funktion Mi(x) = log(x2).
Beweis
Anwenden der Regel (6.1) zeigt diese Behauptung: die Sequenz, die j erzeugt, ist in der Tabelle 6.2 zusammengefasst, wobei die Tabelle 6.1 die im Folgenden verwendete Grammatik visualisiert.
Zum leichteren Verständnis werden nun die ersten Folgenglieder näher erläutert: Die Nichtterminalsymbole dieser speziellen Grammatik lauten
6.1. Grammatical Evolution
S::=
(
x (4)
y (5)
|z (6)
|- (1)
|× (2)
|/ (3)
sin (0)
|cos (1)
|exp (2)
|log (3)
|1 (1)
|2 (2)
|3 (3)
|4 (4)
|5 (5)
|6 (6)
|7 (7)
|8 (8)
|9 (9)
Tabelle 6.1.: Die Grammatik aus Beispiel 6.1
7 mod 7 = 0, d.h.
willkürlicher Vektor ist. ii
Es ist leicht einzusehen, dass die Anwendung genetischer Operatoren auf dieses Ver¬fahren leicht möglich ist. Gemäß der Definition der Operatoren werden beim Crossover zwei Integervektoren k und l zufällig ausgewählt und Teile beider Chromosomen aus¬getauscht, wohingegen beim Mutationsoperator nur ein Vektor stochastisch modifiziert
6.1. Grammatical Evolution
wird. Im Folgenden werden die durch die Grammatik erstellten analytischen Funktionen Modelle genannt und es gelte die Schreibweise Mi(x).
String Chromosom
log(
log(
log(x
log(x*
log(x*
log(x*
log(x*1
log(x*1*
Tabelle 6.2.: Die Dekodierungssequenz aus Beispiel 6.1
Nun mag die Frage aufkommen, was passiert, wenn ein ungültiges Modell erstellt wurde. Ein ungültiges Modell wäre bspw. eine Funktion, die abhängt von mehr als einer Variablen und die Kodierung einer Differentialgleichung sein soll, die nur von einer Variablen abhängt (dies ist durchaus möglich, schließlich enthält obige Grammatik die Ausdrücke x, y und z). Für diesen und weitere Fälle wird das Chromosom verworfen und ein neues generiert, solange ein wohldefiniertes Modell erstellt wurde. Des Weiteren können sogenannte wrapping events implementiert werden. Wrapping events kehren zurück zum ersten Integer des Chromosoms und haben zwei sinnvolle Anwendungen: Einerseits falls das Ende des Chromosoms erreicht ist, aber noch nicht alle Nichtter-minalsymbole auf ein Terminalsymbol abgebildet wurden; und andererseits zollen sie dem fundamentalen Grundgedanken evolutionärer Programmierung Rechnung: der Er¬höhung der Fitness durch ein Ausschöpfen eines möglichst großen Bereichs der Defi¬nitionsmenge, also der Menge aller realisierbaren Modelle gegeben diese Grammatik. Durch wrapping events wird also die Wahrscheinlichkeit erhöht, eine möglichst große Anzahl verschiedener Modelle zu erzeugen, um damit eine gute Fitness der Lösung des gesamten dynamischen Systems zu erzielen.
Nach diesem einführenden Beispiel in das Konzept der GE soll nun gezeigt werden, wie dieses Konzept dazu verwendet werden kann, Differentialgleichungen zu lösen.
6.2. Anwendung: Lösen von Differentialgleichungen
6.2. Anwendung: Lösen von Differentialgleichungen
Im Folgenden soll ein Verfahren vorgestellt werden, das nach Tsoulos und Lagaris [32] eine erfolgreiche Methode zur Lösung von sowohl gewöhnlichen (ODEs) als auch par¬tiellen Differentialgleichungen (PDEs) sein soll. Die einfachste Form des Verfahrens dient der Lösung von linearen ODEs. Um eine gegebene ODE zu lösen, arbeitet die Methode wie ein GA: zunächst werden im Rahmen der Initialisierung die Grammatik, d.h. ins¬besondere das zur Verfügung stehenden Repertoire an Funktionen bestimmt sowie die Startpopulation generiert; im Anschluss kommt es zur Auswahl von Eltern–Vektoren, zur Anwendung von genetischen Operatoren und zur Fitnessauswertung; diese Schleifen werden ausgeführt bis ein Abbruchkriterium erfüllt ist und die ODE im besten Fall gelöst ist. Da in den vorangegangenen Kapiteln bisher ausschließlich Probleme betrachtet wur¬den, bei denen Funktionen an bestimmten Iterationslösungen ausgewertet werden, jetzt aber Funktionen selbst Iterationslösungen darstellen, liegt besonderes Augenmerk auf der Fitnessauswertung. Eine eindeutig lösbare ODE besteht aus zwei Teilen: der Dif¬ferentialgleichung selbst sowie den Anfangs– bzw. Randbedingungen. Dementsprechend besteht die Auswertung des Fitnesswertes ebenso aus zwei Teilen: Es sei eine ODE in folgender allgemeiner Form
x, y, y', y'', . . . , y(n_1), y(n)
f = 0, x E [a,b]
mit Anfangsbedingungen
gi x, y, y', y'', . . . , y (n_1)x=ti = 0, i = 1,...,n, ti E {a,b}
gegeben.
Die Fitnessauswertung arbeitet mit den Schritten wie im Algorithmus 6.1.
Eine genauere Analyse dieses Algorithmus zeigt: offensichtlich besteht der letztendliche Fitnesswert vi aus zwei Teilen, E und P. Dieser Aspekt rührt daher, dass die gegebene ODE ebenso aus zwei Komponenten besteht. Der Ausdruck
steht für die gegebene ODE selbst, nur dass die vorher erstellten Modelle Mi und ihre Ableitungen die Ableitungen der gesuchten Lösung in der ODE ersetzen und die Modelle an den Stützstellen ausgewertet und im Anschluss über alle N Stützstellen aufsummiert werden. Des Weiteren erfolgt im Anschluss eine Quadrierung dieser mit
6.2. Anwendung: Lösen von Differentialgleichungen
Algorithmus 6.1 : Algorithmus zur Fitnessauswertung
1 Wähle N äquidistante Stützstellen
2 for alle Chromosomen i do
3 Konstruiere die jeweiligen Mi(x)
4 Berechne E(Mi) = PN−1 (f (xj, Mi(xj), Mi(n)(xj)))2j=0
5 Berechne die sogenannte zugehörige Zielfunktion:
P(Mi) = A g2 (x, Mi(n−1)) =ti, , A > 0
x
6 Berechne den Fitnesswert von i wie folgt: vi = E(Mi) + P(Mi)
7 end
Hilfe der Modelle modifizierten ODE. Hier ist es wichtig darauf hinzuweisen, dass die gegebene ODE f in homogener Form gegeben ist. Implizit misst der Ausdruck E(Mi) damit den quadratischen Abstand der Iterations- zur tatsächlichen Lösung. Analog dazu misst die Zielfunktion
wie gut die Iterationslösung die Anfangsbedingungen erfüllt: Auch hier werden diese ersetzt durch die erstellen Modelle sowie ihre Ableitungen, dann quadriert und über die Anfangsbedingungen aufsummiert. Des Weiteren wird dieser Wert mit einem positiven Parameter λ gewichtet. Dieser liegt im Ermessen des Programmierers und spiegelt die Relevanz wider, die den Anfangsbedingungen bei der Lösung des Problems beigemessen wird. So stelle man sich eine einfache lineare Funktion durch eine einzige gegebene An-fangsbedingung vor; diese wäre eine perfekte Approximation an jene Anfangsbedingung mit quadratischem Abstand von Null. Daher erscheint es sinnvoll, dem Term E(Mi) eine größere Relevanz zuzuordnen; Koza bspw. schlägt einen Wert λ = 0,25 vor, d.h. er misst dem Term E eine Relevanz von 75 Prozent und dem Term P eine Relevanz von 25 Prozent bei der Lösung einer ODE bei [22]. Der letztendliche Fitnesswert ergibt sich dann als Summe beider Komponenten. Durch die Messung der quadratischen Abstände ist ersichtlich: eine Funktion mit geringerem Wert v ist im Sinne des Algorithmus fitter als eine Funktion mit großem v.
Im Folgenden soll ein Beispiel zur Anwendung der Methode betrachtet werden. Es sei folgende ODE mit dazugehörigen Anfangswerten gegeben:
y00 + 100y = 0, x E [0, 1], y(0) = 0, y0(0) = 10.
6.2. Anwendung: Lösen von Differentialgleichungen
Das Verfahren soll mit 10 Stützstellen im Intervall [0, 1] arbeiten und es sei der
Beispielvektor c = [7, 2, 10, 4, 4, 2, 12, 20, 30, 5] gegeben. Dieser Vektor ist die Kodierung
für folgende Modellfunktion
Mc(x) = exp(x) + sin(x)
mit den ersten beiden Ableitungen
M0c(x) = exp(x) + cos(x), M0'
c (x) = exp(x) − sin(x).
Mit diesen Informationen lässt sich der Fitnesswert des Vektors c leicht berechnen:
9
E(Mc) = X (101 exp(xj) + 99 sin(xj))2
j=0
P(Mc) = ). ((Mc(0) − y(0))2 + ) (M0c(0) − y0(0))2
( (exp(0) + sin(0) − 0)2 + (exp(0) + cos(0) − 10)2)
= λ = 65)
9
Þ vc = X (101 exp(xj) + 99 sin(xj))2 + 65X j=0
Dies ist eine sehr große Zahl und damit eine schlechte Iterationslösung, allerdings dient diese Rechnung auch nur als Beispiel und ist natürlich kein vollständiger Durchlauf des Algorithmus. Vielmehr wurde hier exemplarisch eine Auswertung des Fitnesswertes vorgenommen, aber keinerlei genetische Operatoren angewendet. Tsoulos und Lagaris haben dieses Verfahren mehrmalig auf verschiedene ODEs angewendet und experimentell untersucht. Dabei wurde das Verfahren auf neun verschiedene ODEs erster und zweiter Ordnung jeweils 30 Mal mit der folgenden Initialisierung angewendet: Mutationsrate: 5 Prozent, Replikationsrate: 10 Prozent, Population: 1000, Chromosomenlänge: 50, maximale Anzahl Generationen: 2000, Abbruchkriterium: Zielfitness von 10−7. Eine dieser neun ODEs wird nun genauer betrachtet: sei y00 = −100y, y(0) = 0, y'(0) = 10, x E [0, 1] gegeben. Die analytische Lösung dieses Problems lautet y(x) = sin(10x). Die beiden Autoren haben eine Iterationslösung
y22(x) = sin((sin(−log(4)x((−cos(cos(exp(7)))exp(cos(6))) − 5))))
mit einem Fitnesswert von v = 4200, 5 in Generation 22 dokumentiert. In Generation 27 hat sich der Fitnesswert bereits auf v = 517, 7 erheblich verbessert; und in Generation 59 hat der Algorithmus die exakte Lösung gefunden. Ihre experimentellen Resultate für ODEs haben die beiden Autoren in Tabelle 6.3 zusammengefasst.
6.3. Die Methode zur Lösung von PDEs
ODE min max avg
ODE1 8 1453 653
ODE2 52 1816 742
ODE3 23 1598 705
ODE4 14 1158 714
ODE5 89 1189 441
ODE6 37 1806 451
ODE7 42 1242 444
ODE8 3 702 66
ODE9 59 1050 411
Tabelle 6.3.: Experimentelle Resultate für ODEs
Hierbei geben die Spalten die minimale, die maximale sowie die durchschnittliche An¬zahl an Generationen, d.h. an Durchläufen des Algorithmus an. Es ist bemerkenswert, dass die maximale Anzahl in jedem Fall die Zahl 2000 unterschreitet, d.h., dass das Verfahren die korrekte analytische Lösung in allen Fällen gefunden hat.
6.3. Die Methode zur Lösung von PDEs
Wie bereits erwähnt, lässt sich dieses Verfahren mit kleinen Modifikationen leicht auf die Lösung von PDEs übertragen. Dieses modifizierte Verfahren soll nun im Folgenden kurz vorgestellt werden. Da das Verfahren für PDEs große Ähnlichkeiten zu dem Ansatz für lineare ODEs aufweist, sei hier nur auf den grundätzlichen Algorithmus 6.2 verwiesen; dieser arbeitet wie nachfolgend beschrieben.
Laut Tsoulos und Lagaris [32] ist die Übertragung auf PDEs höherer Ordnung mit anderen Randbedingungen leicht möglich. Neben der Anwendung von GE zur Lösung von linearen ODEs sowie PDEs, lässt sich das Konzept ebenso leicht auf nichtlineare ODEs sowie auf Systeme von Differentialgleichungen übertragen. Die Experimente der
Autoren zeigen, dass falls eine Lösung in geschlossener Form existiert und falls die
verfügbare Grammatik über ein ausreichendes Repertoire an Funktionen verfügt die Methode diese Lösung findet. Bemerkenswert ist weiterhin, dass für den Fall, dass keine Lösung analytisch anzugeben ist, diese Methode eine Approximationslösung in geschlossener Form angibt. Nicht weniger beachtlich ist, dass das Verfahren bereits er¬folgreich in verschiedenen Anwendungsgebieten wie in der Symbolischen Regression, bei der Vorhersage von Finanzdaten () sowie in der Robotik eingesetzt wurde. Darüber hinaus forschen die Autoren an weiteren PDEs mit modifizierten Grammatiken.
6.4. Genetische Programmierung
Algorithmus 6.2 : Lösung von PDEs
x E [x0, x1], y E [y0, y1] und Dirichlet-Randbedingungen
Ψ(x0, y) = f0(y), Ψ(x1, y) = f1(y), Ψ(x, y0) = g0(x), Ψ(x, y1) = g1(x)
1 Wähle N2 äquidistante Stützstellen aus dem Rechteck [x0, x1] x [y0, y1], davon Nx äquidistante Punkte auf dem Rand in x = x0 und in x = x1 sowie Ny Punkte auf dem Rand in y = y0 und in y = y1.
2 for alle Chromosomen i do
3 Erstelle ein Modell Mi(x, y)
4 Berechne die Größe:
E(Mi) = N2
X fxj,yj,Mi(xj,yj),∂
∂xMi(xj,yj),
j=1(6.2)∂ ∂2 02
OyMi(x Ox2 Mi(xj, yj),ay2 Mi(xj, yj) I2
5 Berechne die zugehörigen Zielfunktionen:
P1(Mi) = PNx
j=1 (Mi(x0,yj) − f0(yj))2,P2(Mi) = PNx
j=1(Mi(x1,yj) − f1(yj))2 ,
P3(Mi) = iNy
j=1 (Mi(xj, y0) − g0(xj))2 , P4(Mi) =PNy
j=1 (Mi(xj, y1) − g1(xj))2
6 Berechne den Fitnesswert von i wie folgt: vi = E(Mi) + λ E4l=1Pl(Mi).
7 end
6.4. Genetische Programmierung
Wie bereits im ersten Abschnitt erwähnt, sind Genetische Algorithmen als ein Teil dessen zu sehen, was man unter der Terminologie Evolutionary Computing subsumieren kann; neben den bereits bekannten GAs ist ein weiteres Mitglied dieser Familie das sogenan¬nte Genetic Programming, im Folgenden GP genannt. In den folgenden Abschnitten soll gezeigt werden, dass es sinnvoll ist, den Fokus weg von GAs hin auf diesen neuen Ansatz zu lenken, da mit letzterem die numerische Lösung von dynamischen Systemen bedeutend einfacher sein kann. Zwischen den beiden Konzepten GA und GP existieren zwei fundamentale Unterschiede, wobei sich der hohe Verwandtschaftsgrad zwischen GAs und GP im ersten technischen Unterschied wiederspiegelt: Genetische Programmierung ist aus der Implementierungsperspektive betrachtet prinzipiell ein GA, mit dem Unter¬schied, dass die Kodierung nicht mit Hilfe von Bit Strings, sondern mit Hilfe von Parse Trees oder auch Syntaxbäumen erfolgt. Ein Beispiel im folgenden Abschnitt wird diese
6.4. Genetische Programmierung
technische Feinheit genauer beleuchten. Der zweite fundamentale Unterschied zwischen GP und GA betrifft weniger die technischen Attribute als die Art der Probleme, zu deren Lösung beide Konzepte sinnvoll eingesetzt werden: Wie bereits in der einleit¬enden Motivation dieser Arbeit dargestellt, werden GAs generell dazu eingesetzt, eine gegebene Funktion numerisch zu optimieren. Jene Funktion kann durchaus beliebig komplex in Bezug auf ihre Optimierbarkeit sein; d.h. sie muss weder differenzierbar, noch stetig oder in analytischer Form gegeben sein. Sie kann sogar einen diskreten Def¬initionsbereich besitzen, wie es beim Travelling Salesman Problem der Fall ist. In all diesen Fällen jedoch ist eine Funktion, d.h. ein gegebenes Problem vorhanden, das mit¬tels eines GAs numerisch optimiert werden soll. GP verwendet dagegen einen anderen Ansatz: Generell wird dieses Konzept dazu eingesetzt, ein bestimmtes abstraktes Mod¬ell, z.B. eine Funktion oder eine Formel zu finden, die eine maximale Fitness aufweisen. Vor dem Hintergrund, dass GP dazu eingesetzt werden soll, dynamische Systeme zu lösen, erscheint dieser neue Ansatz damit durchaus vielversprechend; schließlich ist eine Differentialgleichung eine Gleichung, deren Lösungen unbekannte Funktionen sind. Ein Beispiel soll diese durchaus abstrakte Grundüberlegung konkretisieren.
Beispiel 6.3 (Bonitätsproblem)
Um Kreditausfallrisiken zu minimieren, müssen Banken die Zahlungszuverlässigkeit ihrer Kunden, d.h. der Kreditnehmer untersuchen und pflegen. Diese Information soll dazu genutzt werden, um ein Modell erstellen, das sogenannte gute von schlechten Kunden ab-grenzt. In einem weiteren Schritt kann das Modell dazu genutzt werden, die zukünftige Zahlungsmoral von Bestands– und Neukunden vorauszusagen. Das generierte Modell kann verschiedene, zur Vorhersage von Zahlungsausfällen nützliche Inputdaten verwen¬den, bspw. das jährliche Gehalt, die Anzahl Kinder des Kunden, den Familienstand etc. Das simple Modell soll einen binären Output bestimmen, wobei 0 eine schlechte und 1 eine gute Zahlungsmoral kodiert. Selbstverständlich arbeiten Kreditinstitute und Ratingagenturen in der Praxis mit weitaus mehr Inputvariablen und einer größeren Klas¬sifizierungsmenge vergleiche z.B. die Bandbreite an Bonitätsnoten von D bis AAA der größten Ratingagentur Standard and Poors das einfache Beispiel hier erfüllt dennoch seinen Zweck, den Begriff Modell im Kontext von GP zu erklären. Tabelle 6.4 fasst dieses konstruierte Beispiel zusammen.
6.4. Genetische Programmierung
Kundennr Janresgehalt brutto Anzahl Kinder Familinstand Bonität
00001 45000 2 verheiratet 0
00002 30000 0 ledig 1
00003 40000 1 verheiratet 1
00004 60000 2 geschieden 1
... ... ... ... ...
10000 50000 2 verheiratet 1
Tabelle 6.4.: Beispieldaten für das Bonitätsproblem
Ein mögliches Modell, das das Problem der Bonitätsvorhersage löst, ist z.B. folgendes: IF(Anzahl Kinder = 2)AND(Jahresgehalt > 80.000)THEN gut ELSE schlecht.
Abstrakter kann man das allgemeine Modell zur Lösung des Problems wie folgt beschreiben:
IF Formel THEN gut ELSE schlecht.
Die Variable Formel ist die einzige Unbekannte in diesem Modell. Das Ziel in diesem Problem ist es also, die optimale Formel zu finden; diese wiederum bestimmt das op¬timale Modell und damit die Lösung des Bonitätsproblems. Technisch ist das Bonität-sproblem als Suchproblem über den Raum aller möglichen Formeln definiert. Eine genaue Definition dieses Raumes soll hier außer Acht gelassen werden; vielmehr steht hier die motivierende Beschreibung des Problems im Mittelpunkt.
Ausführungen zur Kodierung folgen im nächsten Abschnitt. Jedoch kann man die Fit-ness dieses Modells beschreiben: Die Qualität einer Formel wird definiert als der Anteil der durch das obige Modell korrekt klassifizierten Kunden. Dieses Beispiel legt nahe, dass es auch vor dem Hintergrund von GP die bekannten Metaphern aus der Darwin-schen Evolutionsbiologie gibt, schließlich liegen GP und GAs als Teile der Evolutionary Computing–Familie nicht allzu sehr auseinander. Tabelle 6.5 fasst diesen Aspekt kurz zusammen.
Kodierung Syntaxbäume
Crossover Austausch von Teilbäumen
Mutation Zufallstausch im Syntaxbaum
Phänotyp Formel
Chromosom, Gen Syntaxbaum, Knoten
Tabelle 6.5.: Zusammenfassung der wichtigsten Attribute der GP
6.4. Genetische Programmierung
+
* —
/
y +
5 1
Abbildung 6.1.: Syntaxbaumkodierung eines arithmetischen Ausdrucks
Kodierung:
Wie bereits erwähnt erfolgt die Kodierung bei GP methodisch mit Hilfe von Syntaxbäu¬men. Prinzipiell ist die Kodierung von vielen verschiedenen Ausdrücken denkbar je nach Problemstellung und Syntax. Hervorzuheben hierbei sind drei zentralen Klassen von Ausdrücken: arithmetische Formeln, logische Formeln sowie Programmcode; diese Aus-drücke kommen besonders häufig zum Einsatz, weshalb deren Kodierung im Folgenden (teilweise) im Fokus stehen soll.
Betrachte folgenden einfachen arithmetischen Ausdruck:
( )
(+ (* 2 ) − (+ x 3) (/ y (+ 5 1))) 2 + (x + 3) − y
5+1
Die erste Teil vor dem Implikationspfeil gibt den Ausdruck in sogenannter symbolic expression an; wie man anhand von Abbildung 6.1 erkennen kann, ermöglicht diese Notation einen leichten Zugang zur äquivalenten Kodierung als Syntaxbaum: Das + ist der globale Operator, daher thront er als Wurzel des Baumes ganz oben. Links und rechts verzweigen die beiden Summanden 2ir bzw. (x + 3) − (y/(5 + 1)) als Äste ab, wobei der erste Summand mit den beiden Blättern 2 und ir beendet ist, wohingegen der rechte Ast am Knoten weiter verzweigt wird.
Die Kodierung von logischen Formeln sowie von Programmcode ist sehr ähnlich; für eine nähere Erklärung siehe (Eiben und Smith [6]).
Im Folgenden soll noch kurz auf die genetischen Operatoren Mutation sowie Crossover eingegangen werden. Analog zu GAs erfolgt der Mutationsvorgang über den zufälligen Austausch einer Teilmenge von Genen bei einem Exemplar der bestehenden Popula¬tion. Der einzige Unterschied zum GA ist, dass die Gene hier nicht aus Strings, sondern aus Teilbäumen bestehen: Es wird zufällig ein Knoten des Syntaxbaumes, der die El-terngeneration repräsentiert, ausgewählt und durch einen wiederum zufällig generierten
6.4. Genetische Programmierung
neuen Teilbaum ersetzt. Abbildung 6.2 visualisiert diesen Prozess anhand des obigen arithmetischen Ausdrucks.
+
Abbildung 6.2.: Syntaxbaumkodierung und Mutation
Der Mutationsvorgang bei GP hat also zwei Parameter: einerseits die Mutationsrate, d.h. die Wahrscheinlichkeit, mit der ein Gen verändert wird sowie die Wahrschein¬lichkeit, mit der ein Knoten des Syntaxbaumes als Wurzel des neuen Teilbaums aus¬gewählt wird. Darüber hinaus ist es erwähnenswert, dass damit die neue Teilbaumlänge die Länge des Elternbaumes übersteigen kann.
Im Folgenden soll noch ein bestimmter Rekombinationsvorgang betrachtet werden, das sogenannte Subtree Crossover. Auch hier läuft der Vorgang prinzipiell analog zum Crossover bei GAs ab: Es werden im Rahmen der Vorauswahl zwei Elternbäume aus-gewählt und Teile beider Individuen werden ab einem zufällig ausgewählten Knoten getauscht, wie Abbildung 6.3 unterstreicht.
( )
2ir + (x + 3) − y (P arent1) ⇒ 2ir + ((x + 3) − 3a) (Child1)
5+1 ( y )
3a(3 + (y + 12)) (Parent2) ⇒ (3 + (y + 12)) (Child2)
5+1
Auch der Rekombinationsvorgang beim GP besteht damit aus zwei Parametern: der Replikationsrate, wobei diese die Wahrscheinlichkeit angibt, mit der Crossover auf ein Chromosomen angewendet wird; und der Wahrscheinlichkeit, mit der ein innerer Knoten jedes Elternbaumes als Crossover-Punkt ausgewählt wird. Es gilt auch hier, dass die neue Teilbaumlänge die der Elternbäume übersteigen kann. Es wurde also gezeigt, dass GP und GA tatsächlich zwei sehr ähnliche Konzepte darstellen. Sieht man von der unterschiedlichen Datenstruktur beider Verfahren ab, so ist der einzige technische
fundamentale Unterschied, dass bei GP die Chromosomen in ihrer Länge gemessen
anhand der Anzahl an Knoten - variieren können, wohingegen die Stringlänge bei GAs in der Regel fix ist (Eiben und Smith [6]).
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
Abbildung 6.3.: Syntaxbaumkodierung und Subtree Crossover
Diese intendiert kurze Einführung in die GP und ihren Instrumenten bildet die Grund¬lage zum besseren Verständnis für die folgenden Anwendungen auf numerische Lö-sungsmethoden von Konvektions-Diffusions-Gleichungen.
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
Genetic Programming (GP) soll nun dazu verwendet werden, die Konvektion–Diffusion– Gleichung (KDG) in einer Variablen zu lösen. Der Zweck der Methoden, die im weiteren Verlauf vorgestellt werden, soll es nicht sein, eine eindimensionale KDG zu lösen, da deren Lösung einfach zu bestimmen ist. Die Methoden sollen stattdessen komplexere (partielle) Differentialgleichungen lösen, werden aber zunächst an der KDG getestet. Die KDG, die hier gelöst wird, lautet
x E [0, 1]
T(0) = 1
T(1) = 0
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
mit der Peclet–Zahl Pe. Die analytische Lösung T ist gegeben durch
Die folgende Abbildung 6.4 zeigt die Lösungen der KDG für verschiedene Peclet–Zahlen.
Abbildung 6.4.: Analytische Lösung der KDG
In den nächsten drei Abschnitten werden drei Methoden beschrieben, wie man eine (zumindest näherungsweise) Lösung der obigen KDG mittels GP berechnen kann. Die drei Methoden basieren auf einer Approximation der analytischen Lösung durch math¬ematische Funktionen wie z.B. Polynome, Sinus–, oder Exponentialfunktionen. Die er¬ste Methode ist von Howard und Roberts [16] und legt bereits zu Beginn die Art der approximierenden Funktion als Polynom mit variabler Länge fest, das die Randwertbe-dingungen a priori erfüllt. Aufgabe der GP–Methode ist es dann, die Koeffizienten und somit den Grad eines Polynoms so zu bestimmen, dass das Polynom eine möglichst gute Approximation darstellt. Die zweite Methode, die von Koza [22] stammt, legt die ap-proximierende Funktion nicht so starr fest. Diese GP–Methode findet eine Komposition von typischen mathematischen Funktionen (+, -, exp, sin, ...) mit der Variablen x als approximierende Lösung. Die Randwertbedingungen sind in dieser Variante nicht a priori erfüllt. In der dritten Methode, die aus Howard und Kolibal [15] ist, wird die GP– Methode mit stochastischer Bernstein Interpolation durchgeführt. Im weiteren Verlauf
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
werden die drei Methoden genauer vorgestellt.
6.5.1. Methode nach Howard und Roberts
Howard und Roberts [16] legen die approximierende Funktion Tˆ a priori durch
Tˆ(x) = x(1 − x)p(x) + (1 − x) (6.5)
mit
p(x) = a0 + a1x + a2x2 + a3x3 ...
fest, wobei p ein Polynom variabler Länge ist. Die Wahl dieser Definition hat zwei Vorteile. Die Randwertbedingungen sind durch diese Definition a priori erfüllt, da
Tˆ(0) = 0(1 − 0)p(0) + (1 − 0) = 1
Tˆ(1) = 1(1 − 1)p(1) + (1 − 1) = 0
ist. Da (6.5) ein Polynom ist, sind die ersten beiden Ableitungen ebenfalls Polynome:
0Tˆ
= x(1 − x)
∂x
+ (1 − 2x)p − 1
∂x
02 T;
ax2 = x(1 − x) ∂2p ∂x2
+ 2(1 − 2x) ∂p ∂x 2p
Diese Tatsache wird sich bei der Bestimmung der Fitness noch als nützlich erweisen. Die Aufgabe der GP-Methode ist es nun, Koeffizienten
[a0, a1, a2, . . . , an]
und somit ein Polynom p vom Grad n E N mit
p(x) = a0 + a1x + a2x2 + · · · + anxn
zu finden, so dass die entsprechende Funktion Tˆ eine gute approximative Lösung der KDG ist.
Die Methode von Howard und Roberts nutzt GP, um die KDG zu lösen. Man könnte das Prinzip dieser Methode jedoch auch nutzen, um die KDG mit GA zu lösen. Da
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
das Polynom von variabler Länge ist und die Methode selbst einen Grad festlegen soll, muss man eine GP–Methode anwenden. Wenn aber der Grad und somit die Anzahl der Koeffizienten vorab festgelegt wird, kann man auch GA anwenden, um die KDG zu lösen.
Im weiteren Verlauf wird die genaue Durchführung der GP–Methode beschrieben. Die approximierenden Funktionen werden, wie es in GP üblich ist, durch Bäume dargestellt. Jeder Baum repräsentiert also jeweils einen Lösungsvorschlag Tˆ , so dass es das Ziel der GP–Methode ist, den Baum mit der besten Fitness zu finden. Für die Methode von Howard und Roberts bedeutet dies, dass jeder Baum einen Koeffizientenvektor repräsen¬tiert. Die GP–Methode durchläuft folgende standardmäßige Schritte 6.3:
Algorithmus 6.3 : Ablauf der GP
1 Setzte k = 1 und bestimme max. Anzahl von Durchläufen kmax
2 Erstelle eine Startpopulation von Bäumen
3 Berechne Fitness der Bäume der Startpopulation und speichere den fittesten in Tˆ
4 while Tˆ keine ausreichende Approximation und k < kmax do
5 GP–Operatoren erzeugen neue Generation von Bäumen
6 Berechne Fitness jedes Baumes der Generation und speichere den fittesten
inTˆ
7 Setzte k=k+1
8 end
Die einzelnen Schritte werden nun von uns näher erläutert.
Schritt 1 und 4: Abbruchkriterium
Die obige GP Methode berechnet Approximationen, bis eine von ihnen gut genug oder eine Höchstzahl an Durchläufen erreicht worden ist. Howard und Roberts geben in ihrer Arbeit jedoch nicht an, welche Höchstzahl sie verwendet haben.
Schritt 2: Startpopulation
Bevor die Startpopulation erstellt wird, müssen einige Parameter der Methode bestimmt werden. Die folgende Tabelle 6.6 zeigt die Wahl von Howard und Roberts.
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
Parameter Wahl
Population
Blätter (terminals) Funktionen
max. Baumgröße { 0,
+, -, ADD,
Wm1, Wm2, 8000
1
255, }
2255,...,1
*, / (x/0:=1) WRITE, BACK
Rm1, Rm2, m1, m2
1000
Tabelle 6.6.: Parameter
Howard und Roberts wählen eine Größe der Startpopulation von 8000 Bäumen, die jeweils eine approximative Lösung T ˆ darstellen. Die Blätter (terminals) ihrer Bäume
{ }
sind numerische Konstanten aus der Menge 0, 1
255, 2 255, . . . , 1 . Die restlichen Knoten stellen Funktionen dar, die die Konstanten verändern (mathematische Grundfunktionen +, -, *, / mit x/0:=1) und sie in den Vektor des jeweiligen Baumes schreiben (übrige Funktionen). Zwei Zeiger C und L, die auf Stellen des Vektors zeigen und durch die Funk¬tionen verschoben werden können, bestimmen, an welche Stellen die Werte geschrieben werden. Ein festzulegender Wert Lmax begrenzt die Länge des Koeffizientenvektors. Als nächstes werden die Funktionen erklärt, die keine mathematischen Grundfunktio¬nen sind. Die folgende Graphik zeigt ein Beispiel eines Baumes, mit dessen Hilfe zwei Funktionen anschaulich erklärt werden, wobei die Zeiger zu Beginn L = C = 0 sind und Lmax > 1 beliebig ist.
+
ADD WRITE
0.6 0.8 1 0.4
Abbildung 6.5.: Beispiel eines Baumes
Die Funktion ADD hat zwei Argumente, übergibt ihr zweites und schreibt ihr erstes Argument an die Stelle des Vektors, auf den der Zeiger L zeigt. Anschließend wird L erhöht, sofern L < Lmax gilt, und die Zeiger gleichgesetzt (C = L). Dies bedeutet für das Beispiel, dass ADD 0.8 an die Funktion + übergibt, den Wert 0.6 an die erste Stelle des noch leeren Koeffizientenvektors schreibt (L = 0) und C = L = 1 setzt. Die Funktion WRITE verhält sich ähnlich. Sie hat zwei Argumente, übergibt ihr zweites und schreibt ihr erstes Argument an die Stelle C des Vektors. Im Fall C < Lmax wird
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
C und im Fall C > L zusätzlich noch der Zeiger L um eine Einheit erhöht. WRITE übergibt im Beispiel also 0.4 an die Funktion +, schreibt den Wert 1 an die zweite Stelle des Koeffizientenvektors (C = 1) und erhöht die Zeiger auf C = L = 2. Der Koeffzientenvektor des Baumes ist also [0.6, 1] und die zugehörige Approximation somit Tˆ(x) = x(1 − x)(0.6 + 1x) + (1 − x) mit p(x) = 0.6 + 1x.
Die Funktion BACK hat ebenfalls zwei Argumente. Im Gegensatz zu den vorherigen beiden Funktionen schreibt BACK jedoch keines ihrer Argumente in den Koeffizienten-vektor, sondern übergibt lediglich ihr erstes Argument und verringert den Zeiger C um eine Einheit (C > 0 vorausgesetzt). Die übrigen Funktionen haben einen anderen Zweck. Die Funktion Wm1 bzw. Wm2 hat zwei Argumente, übergibt ihr erstes und speichert ihr zweites Argument in m1 bzw. m2. Der gespeicherte Wert wird von der Funktion Rm1 bzw. Rm2 abgerufen, die ihrerseits ihre Argumente ignoriert.
Als nächstes wird noch erklärt, wie ein Baum genau entsteht. Howard und Roberts legen eine maximale Baumgröße von 1000 Knoten fest. In der Literatur findet man häufig zwei Methoden, um einen Baum aufzubauen (Koza [22]). Eine Möglichkeit ist es, so lange gleichmäßig Knoten an einen Baum anzufügen, bis dieser die maximale Größe erreicht hat. In dieser Full Method genannten Variante haben alle Bäume die maximale Größe. Die zweite Möglichkeit wird Grow Method genannt und lässt Bäume unterschiedlicher Größe zu. Es wird zunächst eine Funktion als Wurzel gewählt. Der nächste Knoten ist dann mit einer zu wählenden Wahrscheinlichkeit P ein Knoten, der kein Blatt ist (also eine Funktion), und mit einer Wahrscheinlichkeit 1-P ein Blatt (also eine Konstante). Dieses Vorgehen wird durchgeführt, bis nur noch Blätter am Ende des Baumes vorzufinden sind. Da die durchschnittliche Baumgröße, wie man unten bei der Darstellung der Ergebnisse sehen kann, bei den verschiedenen Durchläufen von Howard und Roberts variieren, verwenden beide nicht die Full Method, um ihre Bäume aufzustellen.
Schritt 3 und 6: Fitnessfunktion
Die GP–Methode soll einen Baum bzw. eine Approximation finden, die die KDG möglichst gut löst. Die Güte einer Approximation wird durch die Fitnessfunktion bes¬timmt. Howard und Roberts wählen dafür folgende Funktion:
0 1 2
Z 1 @∂2 Tˆ
F = − ∂x2 − P e∂T A dx (6.6)
0 ∂x
Der Integrand entspricht dem quadratischen Fehler (im Vergleich zur rechten Seite der KDG), den die Approximation T ˆ beim Einsetzen in die KDG verursacht. Ein höherer
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
Fitnesswert bedeutet in diesem Fall eine bessere Approximation. Die analytische Lösung der KDG ist die Maximalstelle der Fitnessfunktion und führt zu F = 0. Eine Betra¬chtung der Definition der Fitnessfunktion macht deutlich, warum Howard und Roberts die obige Darstellung der Approximation Tˆ gewählt haben. Da T ˆ ein Polynom ist, sind die Ableitungen ebenfalls Polynome und die Fitnessfunktion ist somit analytisch berechenbar.
Schritt 5: GP–Operatoren
Wenn innerhalb einer Generation keine ausreichend gute Approxmation gefunden worden ist, wird mittels GP–Operatoren eine neue Generation erstellt. Die Details zu den GP– Operatoren sind in der folgenden Tabelle dargestellt.
Parameter Wahl
Operatoren 80% Crossover, 20% Clone
Kill Tournament Größe: 2
Breed Tournament Größe: 4
Tabelle 6.7.: Parameter zu GP-Operatoren
Howard und Roberts verwenden den Crossover– und den Clone Operator. Bei Ver-wendung des Crossover Operators wird in zwei Bäumen jeweils ein Knoten zufällig gewählt und die darunterliegenden Subbäume getauscht. Beim Clone Operator wird ein gewählter Baum direkt in die neue Generation übernommen. Die Auswahl der Bäume für diese Operatoren erfolgt durch ein Breed Tournament der Größe vier. Es werden in diesem Fall vier Bäume zufällig aus der aktuellen Generation gezogen, wobei der fitteste von ihnen für den Operator gewählt wird. In der klassischen Variante des GP kommen die gewählten Bäume in eine zunächst leere Generation, bis diese so viele Bäume wie die vorherige Generation enthält. In der Steady–State Variante, die Howard und Roberts verwenden, werden die Bäume wieder in die alte Generation zurückgesetzt, wobei jeder zurückeingesetzte Baum einen Baum aus der alten Generation verdrängt. In einem Kill Tournament werden zwei Bäume aus der alten Generation gezogen und der schwäch¬ste Baum wird ersetzt. Eine neue Generation ist dann entstanden, wenn die Methode genau so viele Bäume eingesetzt hat, wie in der alten Generation gewesen sind (siehe z.B. Kinnear [20]).
Ergebnisse
Howard und Roberts haben die Methode mit den obigen Einstellungen verwendet, um eine Lösung der KDG für verschiedene Peclet–Zahlen Pe zu finden. Die Ergebnisse der
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
erfolgreichen Durchläufe sind in den drei folgenden Tabellen 6.8 6.10 dargestellt.
Pe pop gens best F avg tree mins
5 8000 42 -0.000175 55 24.23
20 8000 80 -0.003145 348 72.11
50 8000 80 -0.042962 830 139.61
100 2000 300 -0.258476 958 228.86
Tabelle 6.8.: Informationen
I aI+0 aI+1 aI+2 aI+3
a0 0.964706 0.882353 0.713725 0.717647
a4 0.170588 0.137255 0.447377
Tabelle 6.9.: Pe=5, 7 Koeffizienten
I aI+0 aI+1 aI+2 aI+3
a0 0.996078 0.972549 1.24314 0.615686
a4 0.752941 1.77528 0.980392 0.588235
a8 0.745098 0.611765 2.24359 0.931927
a12 0.586275 0.858824 0.462745 0.92549
a16 0.74902 0.374761 0.374761 0.733333
a20 0.374761 0.410748 0.0313725 0.588235
a24 0.0313725 0.0313725
Tabelle 6.10.: Pe=20, 26 Koeffizienten
Die erste Tabelle zeigt die gewählte Populationsgröße, die Anzahl der benötigten Gen-erationen, den höchsten erzielten Fitnesswert, die durchschnittliche Baumgröße und die Dauer eines Durchlaufes für die verschiedenen Peclet–Zahlen. Eine Approximation der analytischen Lösung durch Polynome wird, wie man an der obigen graphischen Darstel¬lung sieht, mit steigender Peclet–Zahl immer schwieriger, da die Steigung in der Nähe von x = 1 vom Betrag her immer größer wird und im Bereich vorher gegen Null geht. Dieses Problem spiegelt sich in den Ergebnissen der GP–Methode wider. Wenn die Peclet–Zahl steigt, sinkt die Fitness der besten Approximation eines Durchlaufes. Die
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
analytische Lösung kann also immer schlechter durch die Polynome der GP–Methode approximiert werden. Neben der Abweichung von der analytischen Lösung steigt auch der Aufwand der Methode mit zunehmender Peclet–Zahl. Um eine ausreichende Ap-proximation zu erlangen, werden mehr Generationen und Zeit benötigt. Der Grund für diesen Anstieg ist, dass mit zunehmendem Pe Polynome höheren Grades benötigt werden, um eine ausreichende approximative Lösung zu erhalten.
Die beiden anderen Tabellen zeigen jeweils die Koeffizienten für einen Durchlauf mit Pe = 5 und Pe = 20. Werden für Pe = 5 nur sieben Koeffizienten benötigt, sind es bei Pe = 20 schon 26. Obwohl die Anzahl der Koeffizienten steigt, sinkt die Fitness der besten Approximation, da die höhere Peclet–Zahl die Approximation erschwert. Wie erreicht die Methode, dass der Grad der Polynome mit zunehmendem Pe steigt? Ein Blick in die erste Tabelle zeigt, dass die durchschnittliche Größe eines Baumes mit zunehmendem Pe steigt. Der Grund dafür ist, dass größere Bäume während des Durch-laufes eine höhere Fitness aufweisen bei größerem Pe. Deshalb werden bei der Wahl für die GP–Operatoren größere Bäume tendenziell bevorzugt mit der Konsequenz, dass die durchschnittliche Größe in einem Durchlauf steigt. Die folgende Abbildung 6.6 vergle¬icht die beste gefundene Approximation zu Pe = 5 bzw. Pe = 20 mit der analytischen Lösung. Obwohl die Güte der Approximation mit zunehmendem Pe sinkt, liefert die GP–Methode nach Howard und Roberts trotzdem eine relativ gute Approximation der analytischen Lösung für Pe 20, wie man in der Abbildung zu Pe = 5 und an der Skalierung der Achsen in der Abbildung für Pe = 20 sehen kann.
Abbildung 6.6.: Approximationen
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
Modifikationen
Da die oben dargestellte Basismethode nicht geeignet ist, um eine vernünftige Approx-imation für relativ große Peclet–Zahlen (z.B. Pe > 50) zu erhalten, erwähnen Howard und Roberts einige Möglichkeiten, um dieses Problem zu beheben. Eine Möglichkeit basiert auf der folgenden Idee. Zunächst wendet man die GP–Methode auf die KDG mit einer kleineren Peclet–Zahl als eigentlich vorgesehen an. Die letzte Population dieses Durchlaufes wird dann als Startpopulation für die KDG mit der eigentlich hohen Peclet– Zahl verwendet. Howard und Roberts überprüfen außerdem, ob die Wahl anderer nu¬merischer Konstanten bessere Approximationen liefert. In der Basismethode stammen die Konstanten aus dem Intervall [0, 1]. Durchführungen der GP–Methoden mit den Intervallen [0, 0.001] und [−1, 1] führen zu einer Veränderung der Größe der einzelnen Koeffizienten, ohne jedoch eine Verbesserung der Güte oder des Aufwandes zu erreichen. Die Performance der Methode zur Lösung der KDG könnte eventuell dadurch verbessert werden, dass man anstatt der obigen Polynome zum Beispiel Chebyshev oder Legendre Polynome verwendet, die ebenfalls leicht zu differenzieren und intergrieren sind.
Die obige GP–Methode ist verwendet worden, um eine eindimensionale KDG zu lösen. Wünschenswerter wäre jedoch eine Methode, die partielle Differentialgleichungen lösen kann. Howard und Roberts machen einen Vorschlag, wie man mit ihrer Methode folgende KDG in zwei Variablen lösen könnte:
∂2 Tˆ ∂2 Tˆ ∂T ∂T
P ⎛ ⎞ = 0
∂x2 + e ⎝+
∂y2 ∂x ∂y
⎠
T(x = 0) = 1
T(x = 1) = 0
T(y = 0) = 0
T(y = 1) = 0
Als Approximation kann man in diesem Fall dann
Tˆ = xy(1 − x)(1 − y)p + exp(−Fx2)
wählen, wobei p nun ein Polynom in xiy ist und F relativ groß gewählt werden muss (z.B. F > 104), damit die Randbedingungen Tˆ(x = 1) = Tˆ(y = 0) = Tˆ(y = 1) = 0 zumindest näherungsweise erfüllt sind.
Insgesamt ist die GP–Methode von Howard und Roberts gut geeignet, um eine eindi-mensionale KDG zu lösen. Um jedoch kompliziertere (partielle) Differentialgleichungen
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
zu lösen, benötigt man eine Weiterentwicklung dieser vorgestellten Methode oder einen anderen Ansatz.
6.5.2. Methode nach Koza
In diesem Abschnitt wird die Methode aus Koza [22] vorgestellt, mit der man mit Hilfe von GP die KDG lösen kann. Das grundlegende Prinzip ist identisch zur Methode von Howard und Roberts. Zunächst erstellt man eine Startpopulation von Bäumen, die jeweils eine approximative Lösung der KDG darstellen. Mit Hilfe von GP–Operatoren werden dann so lange neue Generationen erstellt, bis einer der Bäume eine ausreichende Fitness hat. Der Unterschied zwischen beiden Methoden liegt in der Darstellung der Bäume und in der Berechnung der Fitness.
Zunächst wird erklärt, wie die Bäume bei Koza aussehen. Koza legt die approximative Lösung nicht a priori (z.B. als Polynom) fest, sondern lässt die GP–Methode aus den vorgegebenen mathematischen Funktionen selber Approximationen erstellen. Ein Baum repräsentiert nun nicht mehr einen Koeffizientenvektor für ein Polynom sondern direkt eine mathematische Funktion. Die folgende Tabelle 6.11 zeigt, welche Funktionen in der Methode von Koza verwendet werden könnten.
Parameter Wahl
Blätter (terminals) [0,1], x
Funktionen +, -, *, sin, cos, exp
Tabelle 6.11.: Parameter, Methode nach Koza
Damit ein Baum insgesamt eine Funktion darstellen kann, wird die Menge, aus der die Blätter bestimmt werden, um die Variable x erweitert. Das folgende Beispiel liefert einen Baum, der in der Methode nach Koza gemäß der obigen Tabelle vorkommen kann.
+
sin exp
1 x
Abbildung 6.7.: Beispiel eines Baumes
Dieser Baum liefert als approximierende FunktionTˆ mit Tˆ(x) = sin(1) + exp(x). Der zweite Unterschied besteht in der Berechnung der Fitness. In dieser Methode wird
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
nun nicht die erste und zweite Ableitung der Funktion analytisch bestimmt, um dann die Fitnessfunktion (6.6) anzuwenden. Koza bestimmt stattdessen in dem Intervall [0,1] aus der obigen Tabelle 6.11 200 Stellen xi mit i ∈ 0,...,199 und xi < xi+1.
(und xi,
). Damit werden dann die 200 Punkte xi3
− Pe ~ A∂
berechnet. Wenn Tˆ die exakte Lösung der KDG ist, sind die zweiten Koordinaten der Punkte Null. Falls Tˆ es nicht ist, sind sie ungleich Null. Koza wählt deshalb folgende Fitnessfunktion:
1 99
X− +0125*200lT(0)—t0)1+0.125*2001T(1)—i1(1)1
j=0
(6.7) Die letzten beiden Summanden bestrafen die Fehler bzgl. der Randwertbedingungen, die bei Koza nicht a priori gesichert sind. Der erste Summand bestraft die Fehler beim Lösen der KDG in den einzelnen Punkten. Es werden folglich die betragsmäßigen Abwe¬ichungen von Null in den 200 Punkten aufaddiert. In diesem Fall bedeutet ein kleinerer Fitnesswert eine bessere Approximation, wobei das Minimum bei F=0 liegt. Das Prinzip der Größe der Startpopulation, der maximalen Größe eines Baumes, der GP–Operatoren und der Tournaments ändert sich nicht im Vergleich zur Methode nach Howard und Roberts.
Es bleibt nur noch die Frage, wie das numerische Differenzieren genau durchge¬führt wird. Zunächst werden die entsprechenden Funktionswerte xi, 1'1(xi) bestimmt.
schnitt der Sekantensteigungen von xi_1, t(xi_1) nach xi, 1'7(4 und xi, t(xi) nach
02 T;
xi+1, t(xi+1). Den Wert
∂x2 erhält man analog, indem man die Steigungen der
x=xi
(Sekanten durch die Punkte xi_1,
) und (xi ) ax bzw. (xi, ax
x
(und xi+1, ∂ Tˆ ∂x x=xi+1 ) berechnet. Für die Randwerte x0 und x199 nimmt man einfach
die Steigungen zu x1 bzw. x198. Damit ist die grobe Struktur der Methode von Koza
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
erklärt.
Durchführungen dieser Methode zur Lösung von Differentialgleichungen konvergieren jedoch oftmals nur langsam, da die Fitnessberechnung einen relativ hohen Aufwand erfordert. Die Methode von Koza muss also ebenfalls weiterentwickelt werden, um kom-plexere (partielle) Differentialgleichungen lösen zu können.
6.5.3. Methode nach Howard und Kolibal
In diesem Abschnitt wird eine dritte Methode zur Lösung von DGL und PDGL vorgestellt. Das GP nach Howard und Roberts (Abschnitt 6.5.1) wird mit der stochastischen Bern¬stein Interpolation (SBI), die von Kolibal entwickelt und patentiert wurde, kombiniert. SBI ist eine neue noch nicht komplett erforschte Methode, die aber einen großen An¬wendungsbereich hat. Man kann sie für Daten-Approximation oder -Interpolation, bei Wiederherstellung der Daten, bei der Filtration der Datenfehler oder bei Daten/-Bild Entfaltungen benutzen. Das ist aufgrund einiger Eigenschaften von SBI möglich. Diese Technik ist robust, man bekommt immer eine numerische Konvergenz, hierar¬chisch aufgebaut, d.h. Daten können zu beliebigem Zeitpunkt hinzugefügt werden, um besseres Resultat zu bekommen, und sie ist nicht harmonisch, sondern spektral, d.h. es ist möglich, eine Folge von Komponenten zu konstruieren, deren Summe zu dem gle¬ichen Resultat konvergiert. Diese neue Methode nach Howard und Kolibal [15] wird auf
eindimensionale KDG (Gleichung (6.3)) angewandt dadurch steigt das Konvergen-
zverhalten von O(h4) auf O(h6) im Inneren des Intervalls und zum Rande hin sinkt sie auf O(h1). Die Methode nach Howard und Kolibal lässt sich auch auf mehrdimensionale Probleme übertragen.
Vorbereitungen für den Algorithmus
Für den Algorithmus werden die Bernstein–Funktionen benötigt. Diese Funktionen sind eine Erweiterung der Bernsteinpolynome, wobei die Binomialverteilung zu der Stan-dardnormalverteilung zusammen mit der entsprechenden Änderung der Variablen wird. Dadurch teilen diese Funktionen die Approximationseigenschaften der Bernsteinpoly-nome.
Sei die Funktion f(x) gegeben für xk E [0,1], so dass gilt f(xk) = yk für k = 1,... ,n.
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
Bernstein Funktionen Kn(x) sind gegeben durch:
wobei f stückweise konstant auf (zk−1, zk) mit dem Wert yk und z0 = zk =
(xk+1 + xk)/2 für k = 1, . . . , n − 1, und zn = ∞ ist. Q(x) ist der freie Parameter oder Eigenglättungsparameter, der das Verhalten der Funktion zwischen den Stützstellen bestimmt. Im Bezug auf KDG bestimmt o die Geschwindigkeit, mit der die Lösung entwickelt wird. Große Diffusionswerte verursachen große Änderung, und kleine Diffu-sionswerte — kleine Änderung. erf(x) ist die Gaußsche Fehlerfunktion:
Weiterhin wird eine quadratische (nxn)–Matrix Ann = (ajk), wobei für die Koeffizienten gilt:
[erf (zk+1 xj) + erf (xj zk)1 (6.8)
Q(x) .Q(x)
Aufgrund der Fehlerfunktion in der Definition sind die Interpolanten über das ganze Intervall unendlich oft differenzierbar und glatter, als z.B. die trigonometrischen Funk¬tionen.
Algorithmus 6.4 : Algorithmus zur Lösung der KDG nach Howard und Kolibal
Data : GP entwickelt die Inputdaten {(xi, yi)}, i = 0, . . . , n − 1
1 Konvertiere Daten auf das Intervall [0, 1]
2 Konstruiere Matrix Ann (hängt von xi und freiem Parameter ab) durch die Bernsteinfunktion, gemäß ( Gleichung (6.8))
3 Berechne A−1
nn
4 Konstruiere augmentierende Matrix
Gleichung (6.8))
5 Berechne ˜AmnA−1
nny, um Outputdaten lyil, i = 0, . . . , 7n − 1 zu erhalten
6 Konvertiere Outputdaten lyil zurück (Lösung der KDG)
7 Berechne Ableitungen
8 Berechne Fitnessfunktion
Die Methode arbeitet mit verschiedenen Längen der Input- und Outputdaten und ist zeitaufwendiger, als die Methode nach Howard und Roberts.
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
Lösung der KDG
Im ersten Schritt wird GP nach Howard und Roberts angewendet mit identischen Pa-rametern aus Abschnitt 6.5.1 und liefert einen Datensatz {(xj, yj)}. GP verteilt gleich-mäßig die berechneten Punkte xj mit Werten yj auf dem Intervall (0, 1) und fügt die beiden Randbedienungen (0, 1) und (1,0) hinzu. Weiterhin wird beim Experimentieren festgestellt, dass GP nach Howard und Roberts die Freiheit für Position der Punkte xi schätzt und SBI Verfahren von den Pufferzonen außerhalb des Lösungsintervalls profi-tiert. Deswegen wird folgende Modifikation implementiert. Intervall der Lösung wurde von (0, 1) auf (0.2, 0.8) umgebaut und der GP-Methode wurde erlaubt, die Position der Punkte überall (unter bestimmten Einschränkungen) im Intervall (0, 1) zu wählen. Der Algorithmus für diese Informationen wird folgend implementiert:
Algorithmus 6.5 : Modifikation des Algorithmus 6.4
1 Auswertung des Individuums
2 if Länge des resultierenden Vektors ungerade then
3 verwerfe seinen letzten Term
4 x0 = 0 ist fix und y0 wird auf das erste Element des resultierenden Vektors gesetzt
5 xn = 1 ist fix und yn wird auf das zweite Element des resultierenden Vektors gesetzt
6 Das Tupel des Vektorsvivi+1 nimmt zugehörige xj und yj im Intervall (0, 1)
7 Das Intervall (0, 1) wird in 97 Abschnitte unterteilt, damit die Differenz zwischen Stützstellen xj − xi+1 nicht mehr als 0.01 beträgt
8 Die Position dieser Abschnitte berücksichtigt die neuen Randbedingungen (0.2, 1.0) und (0.8, 0.0)
KDG wird mit 0.6 skaliert und somit lautet die Gleichung (6.3) wie folgt:
∂2T ∂T
0.6 ∂x2 − Pe ∂x = 0 x E [0, 1]
T(0.2) = 1 T(0.8) = 0
Die Fitness wird über die Outputdaten {yj}, i = 0, . . . , m − 1 berechnet. Eine Menge von 241 Punkten wird durch proportionale Verteilung der Punkte im Intervall (0.2, 0.8) variiert. Es wird mit zwei Fitnessfunktionen experimentiert:
1. Identisch mit der Fitnessfunktion aus GP–Methode nach Howard und Roberts
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
(siehe 6.6)
⎛ ⎞2
0.8 ⎝0.6∂2T
F = − ∂x2 − P e∂T ⎠dx
0.2 ∂x
2. Unterteile das Intervall (0.2, 0.8) in 60 Abschnitte. Jeder Abschnitt hat 4 gleich-
⎛mäßige Lücken und beinhaltet 5 Ausgangspunkte. Das Quadrat von ⎝0.6∂2T
∂x2 − P e∂T ∂x
kann mit Hilfe der Fünf–Punkte Boole's rule integriert werden
Ableitung der augmetierenden Matrix ˜Amn bestimmen, z.B. durch Gradientenversion einer Matrix. Eine alternative Methode ist gegeben, durch berechnen der Ableitungen mit Hilfe Vorwärts–, Rückwärts oder zentraler Differenzenquotienten.
Ergebnisse
Beim Experimentieren stellen sich zwei typische Ergebnisse heraus. Die Abbildung
Abbildung 6.8.: GP–SBI–Methode, Pe=1
Quelle: D. Howard, J. Kolibal [15]
6.8 zeigt Ergebnis der Outputdaten yi durch numerische Approximation der GP–SBI– Methode für eine sehr niedrige Peclet–Zahl (Pe = 1), wobei das Verfahren noch nicht modifiziert wurde. Die magentafarbige Linie ist die analytische Lösung und die blaue Linie ist SBI–Approximation. Der eingefügte Graph zeigt die Position und den Wert
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
der Interpolationspunkte. Es wird zu den Randbedingungen zwei interne Punkte durch GP–Methode entwickelt. Die Outputdaten werden bei 200 zufälligen Orten erzielt. o wird konstant auf 1 gehalten.
In den Abbildungen 6.9 und 6.10 ist Pe = 4, o = 1 und das Intervall ist modifiziert. Die
Abbildung 6.9 veranschaulicht auch die von GP–Methode entwickelte Stützstellen bei x = 0 und x = 1 und drei weitere entwickelte Punkte innerhalb des Intervalls (0.2, 0.8), die modifizierten Randbedingungen sind konstante Werte, die von Evolution nicht be¬troffen sind. Der eingefügte Graph zeigt die Abweichung der Werte vom absolutem Wert der KDG an jedem Outputpunkt. Dabei werden andere Randbedingungen, als 0 bzw. 0.2 und 1 bzw. 0.8 gewählt. Das zeigt, dass beliebige Randbedingungen sich mit diesem Verfahren durchführen lassen werden. Die Abbildung 6.10 zeigt den Graph der Ableitung durch SBI–Approximation (blau) im Vergleich zu der analytischen Lösung der Ableitungen (magenta). Diagramm auf der linken Seite ist die erste und auf der rechten
Seite die zweite Ableitung. Die durch GP–Methode formulierte Approximation der
ersten Ableitung ist näher zu der analytischen Lösung, als die zweite Ableitung. Der
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
Fehler in den Endbereichen ist mehr ausgeprägt.
In den Abbildungen zuvor wird o konstant gehalten, da sonst die Punkte nicht inter¬poliert werden und die Randbedingungen nicht erfüllt werden können. Aber ab Pe>10 ist es notwendig, den Wert von o zu verkleinern. Ein Wert von o = 1 führt zu keinen guten Konvergenzeigenschaften bei höheren Peclet–Zahlen. Die Verringerung des Wertes von o bringt signifikante Ergebnisse hervor, die nahe an analytischen Lösungen bei ho¬hem Wert von Pe sind. Die Abbildung 6.11 zeigt die Ergebnisse eines erfolgreichen Testlaufes bei einer hohen Peclet–Zahl (Pe = 20). Der blaue Graph ist die GP–SBI–
6.5. GP zur Lösung der Konvektion–Diffusion–Gl.
Methode und magentafarbiger Graph ist die analytische Lösung. Im Uhrzeigersinn von oben links: Lösung; erste Ableitung; zweite Ableitung; absoluter Fehler; entwickelte und vorgeschriebene Stützstellen. Die Ergebnisse werden durch eine geeignete Koordi-natentransformation vom Intervall (0.2, 0.8) auf (0, 1) zurück skaliert. Diese Abbildung zeigt zwei wichtige Beobachtungen. Im Graph mit dem absoluten Fehler sieht man, wie das GP–Methode den Pufferraum (0.8, 1) verwendet um die Funktion nach oben zu drehen. GP–Methode entwickelt mehr Punkte als erforderlich und einfach häufte sie in die Region nahe bei 1.0 an, um so effektiv diese in einen einzigen Punkt zu drehen. In Testläufen bei Pe = 100 und GP–Population von 10000 ergibt sich ein sehr niedriger Wert für o (o = 0.01). Obwohl der Fehler recht hoch ist, zeigen die resultierenden Werte richtige Form sowohl für die Lösung, als auch für die Ableitungen.
Die durch SBI berechnete Ableitungen sind in allen Testläufen glatter, als die mit finiten Differenzen gewonnenen Ableitungen. Ableitungen, die durch die finite Differen¬zen berechnet werden, weisen eine gewisse Schwingung auf, die als Quelle der Differenzen zur Fitnessfunktion übersetzt werden können. Der Grund liegt in der Tatsache, dass durch SBI berechnete Ableitungen in fast allen Experimenten verwendet werden, und somit glatter sind, als die lokalisierten Ableitungen der finiten Differenzen. Die GP– Methode ist auch mit aus den finiten Differenzen gewonnenen Ableitungen erfolgreich.
7. Symbolische Regression in der
Zeitreihenanalyse
Bei der Regressionsanalyse besteht die Aufgabe darin, eine abhängige Variable durch eine oder mehrere unabhängige Variablen zu beschreiben. Bei der linearen Regression ist die Methode der kleinsten Quadrate das Standardverfahren. In den Wirtschaftswis¬senschaften können viele Prozesse nur schlecht durch lineare Zusammenhänge erklärt werden. Oft liegen vielen Tatsachen komplexere Beziehungen zugrunde. Die linearen Modelle sind in diesen Fällen hilfslos. Ein Ansatz, der nichtlineare Zusammenhänge konstruiert, ist die symbolische Regression. Sie ermöglicht es, die abhängige Variable durch beliebige Funktionen zu erklären, z.B. y =18x2− sin3(x) o.a. Ein wesentlicher Vorteil der symbolischen Regression ist, dass man keine Theorien für die Erklärung der Zusammenhänge zwischen den Variablen unterstellen muss. Sie findet die Struktur der Zusammenhänge rein numerisch.
Nach I. Zelinka, L. Nolle, Z. Oplatkova [17] existieren heutzutage drei Methoden, die für die symbolische Regression verwendet werden können: GP, grammatikalische Evolution ud analytische Programmierung. Die beiden letzten Methoden basieren auf GP, wobei GP für die symbolische Regression am häufigsten verwendet wird. John R. Koza [23] hat als erster vorgeschlagen, die Aufgaben der symbolischen Regression mithilfe genetischer Programmierung zu lösen. Im Prinzip, ist die genetische Programmierung genau die symbolische Regression.
Sehr oft wird die GP für die symbolische Regression für Anwendungen in der Fi-nanzwirtschaft benutzt. M. Santini und A. Tettamazi [29], 2000, benutzen in Rahmen des Wettbewerbs „Dow Jones Prediction“ diese Kombination für die Vorhersage des In¬dex Dow Jones Industrial Average und belegen Platz 1. Die Aufgabe bestand in der Prognose des Index am Tagesschluss für die Periode vom 19. Juni 2000 bis einschließlich 30. Juni 2000 (10 Tage). Die Bewertung der Ergebnisse erfolgte mit der Formel
1 0
X
t=1 11 − t
10 |xt − ˆxt| (7.1)
7.1. Der Algorithmus
wobei xt den Schlussstand am Tag t bezeichnet und ˆxt die Prognose.
7.1. Der Algorithmus
Die Idee ist, eine Funktion mit der symbolischen Regression aus 116 beobachteten Re-alisierungen des Index (vom 3. Januar 2000 bis 17. Juni 2000) zu konstruieren und mit ihr die 10 Prognosewerte zu bestimmen. Die Autoren konstruieren eine Daten¬menge D = x(i)0 ≤i ≤a, wobei a > 0, x(0) = (x1,..., xT) als die beobachtete
Reihe und x(i) = (xn , 4) für 1 i a als eine Hilfsreihe. Das Ziel ist ein
Vektor e = (e1,..., eh), wobei h die Horizontlänge beschreibt (hier h = 10), und ej : D 1,...,T→IR, so dass ej(D, T) Prognose für xT+j mit 1 ≤j ≤h ist.
Es werden
· Terminalsymbole = [−c,c]⊆
· Function Set = +, −, , /, pow, sqrt, log, exp∪U0iadata[i], diff[i], ave[i]
definiert, wobei für α ∈IR und 0 ≤t ≤T, ∀i gilt
data[i](α,t)=x(i)
t−α,
diff[i](α, t) =x(i)
t−x(i)
t−α,
Et
1 ave[i](α, t) = xk .
(i)
1 +αk=t−α
Damit wird sicher gestellt, dass die Prognosen auf den Daten aus der Vergangen¬heit basieren. Jede Population besteht aus N Individuen, wobei jedes Individuum ein Vektor e = (e1,..., eh) ist. Jedes ej, j = 1...h, ist ein String in der umgekehrten pol¬nischen Notation. Jedem Terminalsymbol und jedem Element des Function Sets wird eine Auswahlwahrscheinlichkeit zugewiesen, und jedes ej wird aus diesen Elementen zusammengebaut.
Es wird eine Distanzfunktion δ: 1182 →IR+ mit δ(x, y) = x −y/xdefiniert. Weiter wird für δ der mittlere Fehler für e zur Zeit 0 ≤t ≤T
h
e(e,t) =E δ (xt+j, ej(D, t)) h j=1
7.2. Ergebnisse und Folgerungen
und der mittlere Fehler für e über die ganze Zeit als
1 T−h 1 T− h h
fs(e) = T − h Ee(e, t) =
E E
h(T - h) 6 (xt+j, ej(D, t))
t=1 t=1 j=1
definiert. Ist δ > 0, folgt fs(e) > 0 Ve. fs(e) wird als die Fitness verwendet, fa(e) = (1 + fs(e))−1 — als angepasste Fitness und f(e) = fs(e)/ P ê fs(e) — als normierte Fitness.
Für die Selektion wird die Methode Truncation Selection verwendet, wobei LNiosj– viele (ρs = 0.1) Individuen für die Erzeugung nächster Generation ausgewählt werden. Für Crossover werden zwei Individuen e = (e1,..., eh) und e' = (e' 1,. . . , e'h) ausgesucht, und in jedem Paar ej, e'j für 1 < j < h werden ej, e'j mit Wahrscheinlichkeit ρc = 0.6 ausgetauscht. Die Mutation eines Individuums e = (e1,..., eh) wird als Austausch zweier Werte ej, ek, 1 < j, k < h, j =6 k durchgeführt, d.h. eine Art Crossover innerhalb eines Individuums. Die Werte werden mit Wahrscheinlichkeit 1/|e| ausgewählt und mit Wahrscheinlichkeit ρm = 0.2 ausgetauscht. Da man nicht annehmen kann, dass die Zeitreihe stationäre ist, wird eine kurze Stichprobe der Werte vom 3. Januar 2000 bis 17. Juni 2000 ausgewählt. In der Tabelle 7.1 werden alle Parameter zusammengefasst:
Parameter Bedeutung Wert
T Größe der Stichprobe 116
h Prognosehorizont 10
N Individuenanzahl 500
G Generationenanzahl 20
m Tiefe des Ausdrucks 5
δ Abstandsfunktion |x − y|/|x|
ρs Selektionsverhältnis 0.1
ρc Crossoverwahscheinlichkeit 0.6
ρm Mutationswahrscheinlichkeit 0.1
Tabelle 7.1.: Parameterwerte für Dow Jones IA Prognose
7.2. Ergebnisse und Folgerungen
In der Tabelle 7.2 sind die Ergebnisse dargestellt. Die Lösung ist charakterisiert durch relativ niedrige Abweichung der Prognose vom tatsächlichen Wert. Die letzte Spalte wird mit der vorletzten nach der Gleichung 7.1 berechnet.
'grob gesprochen bedeutet Stationarität die Abwesenheit eines Trends in der Zeitreihe
7.2. Ergebnisse und Folgerungen
10449.300 10557.80 -108.5 108.5 108.5
10449.300 10435.20 14.1 12.69 121.19
10449.300 10497.70 -48.4 38.72 159.91
10417.374 10376.10 41.274 28.8918 188.8018
10476.400 10404.80 71.6 42.96 231.7618
10535.523 10542.99 -7.467 3.7335 235.4953
10435.246 10504.46 -69.214 27.6856 263.1809
10461.648 10527.79 -66.142 19.8426 283.0235
10418.607 10398.04 20.567 4.1134 287.1369
10421.874 10447.89 -26.016 2.6016 289.7385
Tabelle 7.2.: Ergebnisse des Verfahrens von M. Santini und A. Tettamazi
Leider werden die expliziten Lösungen der Autoren und ihrer Konkurrenten nicht veröffentlicht und man kann deswegen keinen direkten Vergleich unterschiedlicher Mod-elle machen. Allerdings hat diese Methode gegenüber anderen Verfahren den 1. Platz im Wettbewerb belegt. M. Santini und A. Tettamazi behaupten, dass weitere Ver¬feinerungen des Verfahrens mithilfe der Modellierung der Markov-Prozesse, ARCH und GARCH-Modellen, die Ergebnisse nicht verbessern konnten.
Miroslav Kfúčik, Jana Juriová und Marian Kfúčik [24] haben das Bündel symbolische Regression - genetische Programmierung mit einem verbreiteten makroökonomischen ARIMA-Modell2 verglichen. Sie behaupten, dass die Weltwirtschaft strukturellen Än¬derungen bevorsteht und die konventionellen Modelle, wie z.B. das ARIMA-Modell, eventuell, somit nicht mehr in Lage sein werden, die Zusammenhänge zwischen volk-swirtschaftlichen Größen zu bestimmen oder zu erklären. Sie sind in ihren Analysen zu dem Schluss gekommen, dass die symbolische Regression mit der genetischen Pro-grammierung mindestens nicht schlechtere Ergebnisse, als ARIMA liefert, und sogar manchmal bessere. Dabei ist nicht zu vergessen, dass keine strukturellen Zusammen¬hänge zwischen den Variablen angenommen werden müssen. Das heißt, im Hinblick auf die möglichen Änderungen der Wirtschaftssysteme, kann die Hälfte der Arbeit erspart werden: man kann die beobachteten Daten sofort mit symbolischer Regression erklären, ohne eine Theorie unterstellen zu müssen Alles in allem lässt sich sagen, dass die sym-bolische Regression mithilfe der genetischen Programmierung ein sehr gutes Verfahren in den Wirtschaftswissenschaften sein kann, dass es vielversprechend ist und dass es viel an Bedeutung gewinnen kann.
t.
8. Zusammenfassung und
Schlussfolgerungen
Diese Arbeit behandelt Theorie und Anwendungsmöglichkeiten genetischer Algorith¬men, welche ein universelles und robustes Optimierungsverfahren darstellen, wodurch eine Anwendung auf eine Vielzahl an Problemstellungen möglich ist. Des Weiteren werden die Möglichkeiten der genetischen Programmierung, beispielsweise zur Lösung von Systemen partieller Differentialgleichungen, aufgezeigt. Genetische Algorithmen abstrahieren die Terminologie der biologischen Evolution nach Charles Robert Darwin basierend auf der Erkenntnis, dass die Anpassung von Lebewesen an die vorherrschen¬den, oder auch sich ändernden Umweltbedingungen als Optimierungsprozess aufgefasst werden kann. Die Prozesse der natürlichen Evolution werden in Form von Operatoren auf Werte des Definitionsbereichs mit dem Ziel einer optimalen Anpassung an das zu-grundeliegende Problem angewendet. Im Gegensatz zu vielen anderen Optimierungsver¬fahren behelfen sich genetische Algorithmen der populationsbasierten Suche, bei der eine festgelegte Anzahl an Startpunkten generiert wird, welche sich während des Durchlaufs dem Optimum annähern sollen. Die Population ist der Garant für die Robustheit des Verfahrens; außerdem resultiert daraus der Vorteil der Überwindung lokaler Extrema während eines Suchdurchlaufs. Zudem weisen genetische Algorithmen einen geringen Implementierungsaufwand auf.
Kapitel 2 gab eine Einführung in genetische Algorithmen, deren Hintergrund und den grundsätzlichen Aufbau. Außerdem wurden die zur Verfügung stehenden Instrumente näher erläutert. Kapitel 3 stellte die zugrundeliegende Theorie durch Herleitung des Theorems von John H. Holland vor. Damit wurde im Ansatz ersichtlich, wie ein genetis¬cher Algorithmus den Suchraum durchläuft. Es wurde klar, dass das Theorem zum völligen Verständnis der Funktionsweise leider nicht ausreichend ist, da es in der Real¬ität an der Komplexität, der dynamischen Natur und modularen Struktur genetischer Algorithmen scheitert. Des Weiteren gibt es für die Wahl der einzelnen verfügbaren Instrumente keine geschlossene Theorie, weshalb man entweder auf Erfahrungswerte aus der Literatur oder eigene Tests angewiesen ist.
Dieser Problemstellung widmete sich das 4. Kapitel, welches mithilfe numerischer Tests einige Erkenntnisse und Gesetzmäßigkeiten für die Anwendung erzielen sollte. Die Tests zu ein- und mehrdimensionalen Funktionen zeigten deutlich, dass die Wahl der einzelnen Methoden und Parameter von dem zugrundeliegenden Problem abhängig ist — was nicht weiter verwunderlich ist. Sie zeigten auch, dass die Wahl der Selektions-methode und der Populationsgröße eine übergeordnete Rolle spielt. Unter einer großen Population und der Tournament Selection konnten in den Tests recht schnell akzeptable Ergebnisse erzielt werden die Wahl der weiteren Parameter war dann eher nebensäch¬lich. Entgegen der Theorie aus Kapitel 3 trug eine höhere Mutationswahrscheinlichkeit in den meisten Fällen zu einer Verbesserung der Lösungsquote bei, was wiederum dur¬chaus verblüffend war. Unter der Kombination des genetischen Algorithmus mit dem Hill-Climbing Verfahren konnte zudem eine Zeitersparnis erreicht werden, was jedoch eine sorgsame Wahl der Parameter und der Schrittweite voraussetzt. Für weiterführende Aussagen und zur Festigung der erzielten Resultate sind weitere Tests für eine Vielzahl an Funktionen und Problemstellungen erforderlich.
Kapitel 5 zeigt eine Möglichkeit der Verwendung genetischer Algorithmen in der Wirtschaft mittels eines expliziten Beispiels auf. Ein genetischer Algorithmus wird für die Lösung des Problems Economic Dispatch entwickelt. Hier ist das ein nichtkon¬vexes Optimierungsproblem. Es müssen die Kosten der Stromproduktion an mehreren Kraftwerken minimiert werden. Zusätzliche Restriktionen (ramp rate limits, prohibited zones und Übertragungsverluste) machen konvexe Optimierungsverfahren und andere konventionellen Verfahren, wie z.B. Lambda-Iterations-Methode, unbrauchbar. Es stellt sich hieraus, dass auch unter Vereinfachungen des Problems durch Weglassen dieser Re¬striktionen der GA bessere Ergebnisse liefert, als andere Methoden. Die Kodierung des Problems macht das Verfahren vor allem für große Stromnetze attraktiv.
Im Anschluss wurde im Rahmen des Kapitel 6 eine Anwendungsmöglichkeit auf die numerische Lösung von Differentialgleichungen aufgezeigt. Hier lag zunächst der Fokus auf der Analyse zweier mit genetischen Algorithmen verwandter Konzepte: Grammatical Evolution sowie Genetische Programmierung. Das hier vorgestellte Verfahren ist eine ef-fektive Möglichkeit zur Lösung verschiedener dynamischer Systeme, insbesondere zur Lö-sung von ODES in linearer und nichtlinearer Form sowie von PDEs. Der Grundgedanke des Verfahrens ist es, mit Hilfe einer gegebenen kontextfreien Grammatik Funktionen in geschlossener Form zu generieren, diese mit Hilfe genetischer Operatoren zu modifizieren und so ein gegebenes dynamisches System zu lösen. Auch wenn dieses Verfahren struk-turelle Unterschiede zu einem genetischen Algorithmus aufweist, so zeigt sich die nahe Verwandtschaft beispielsweise anhand von Kodierung, der Durchführung von genetis-
chen Operatoren sowie dem Ziel der Fitnessmaximierung. GP arbeitet strukturell wie ein genetischer Algorithmus, nur dass die Kodierung mit Hilfe von parse trees realisiert wird. Der große Vorteil von GP ist dabei die Möglichkeit, sehr abstrakte Problemstel¬lungen zu lösen, wie anhand des Bonitätsproblems 6.3 gezeigt wurde. Auf Basis dieser theoretischen Einführung wurde im Folgenden GP auf die numerischen Lösungspoten¬ziale der Konvektion–Diffusion–Gleichung (KDG) angewendet. Hier wurden zunächst drei Methoden vorgestellt, mit denen man allgemeine PDEs lösen kann. Diese Metho¬den lösen die KDG relativ gut; um jedoch komplexere PDEs zu lösen, bedarf es weiterer Forschung und einer Weiterentwicklung dieser Ansätze.
In Kapitel 7 wird eine Anwendung der genetischen Programmierung und der symbol¬ischen Regression in der Finanzwirtschaft vorgestellt. Es wird ein Verfahren entwickelt, das den zukünftigen Stand (10 Tage) des Aktienindex Dow Jones Industrial Average prognostiziert. Die Methode hat sich gegen viele andere durchgesetzt und erreichte 2000 den 1. Platz im Wettbewerb „Dow Jones Prediction'.
Alles in allem ist das Feld der Einsatzmöglichkeiten weitläufig. Dank ihrer Struktur und Robustheit können genetische Algorithmen auf eine Vielzahl von Problemstellungen erfolgreich angewendet werden. Anwendungsgebiete der GA sind z.B.: Finanzwirtschaft [26] und [3], Course Timetable [2], Spieltheorie [18], Behavioral Economics [4], Job scheduling [13], Ressourcenwirtschaft [8] und viele weitere. Andere Anwendungsgebi¬ete der GP und der symbolischen Regression: Economics [31] und [24], Finanzwirtschaft [29]. Das Feld kann zudem durch Modifikationen am Algorithmus auf neue Probleme ausgeweitet werden, was ebenso für die genetische Programmierung gilt. Nachteilig sind natürlich der höhere Zeitaufwand und die benötigte Speicherkapazität, schließlich können die genetischen Operatoren nur auf eine multiple Anzahl von Definitionspunk¬ten angewendet werden. Je größer die Populationsgröße, umso schneller terminiert der Algorithmus. Der große Vorteil genetischer Algorithmen und der genetischen Program¬mierung ist gleichzeitig auch der größte Nachteil: robuste und vielfältige Verfahren wer¬den meist solange eingesetzt, bis sie von einem spezialisierten Verfahren abgelöst werden oder ihren Kostenvorteil verlieren. Existiert kein spezialisiertes Verfahren, wäre die Entwicklung mit Kosten verbunden, was die Nutzung eines bestehenden Verfahrens na¬helegt. Existieren aber bereits spezialisierte Verfahren, haben genetische Algorithmen gegenüber diesen wegen ihrer geringeren Effizienz oftmals das Nachsehen, womöglich auch in der Kostenfrage (hierbei müssen natürlich wesentlich mehr Faktoren berück¬sichtigt werden). Durch die fortschreitende Entwicklung werden genetische Algorithmen mit der Zeit wohl nach und nach von spezialisierten Verfahren verdrängt werden, denn Speicherkapazität und Zeit sind wichtige Kostenfaktoren. Dem kann nur gegengesteuert
werden, indem neue Einsatzgebiete erschlossen werden oder der bestehende Algorith¬mus für gewisse Problemstellungen perfektioniert wird. Trotz ihrer Nachteile stellen genetische Algorithmen einen interessanten und alternativen Ansatz zur Lösung unter¬schiedlichster Problemstellungen dar und mit neuen Problemen entstehen auch neue Chancen.
Keine Kommentare:
Kommentar veröffentlichen
Hinweis: Nur ein Mitglied dieses Blogs kann Kommentare posten.