Laden...

[Artikel] Simulated Annealing

Erstellt von gfoidl vor 14 Jahren Letzter Beitrag vor 14 Jahren 35.354 Views
gfoidl Themenstarter:in
6.911 Beiträge seit 2009
vor 14 Jahren
[Artikel] Simulated Annealing


(Abbildung1: Beispielanwendung zur Minimierung einer rellen Funktion)

Simulated Annealing

Simualted Annealing (SA) ist ein meta-heuristisches Optimierungsverfahren welches auf kontinuierliche und auf diskrete Probleme angewandt werden kann. Es besitzt den großen Vorteil dass lokale Minima verlassen werden können.

Dieser Artikel ist in zwei Teilen verfasst. Ausführliche Grundlagen sind in einem separaten Teil als PDF verfasst worden. Dieser Teil beinhaltet eine Kurzeinführung und zwei typsiche Beispiele für SA:*Minimierung einer reellen Funktion als Vertreter für ein kontinuierliches Problem *Problem des Handelsreisenden als Vertreter für ein diskretes Problem

Exkurs Wärmebehandlung (hilfreich für das Verständnis)
In der Metallurgie wird ein Metall erhitzt und somit in einen hohen Energiezustand versetzt. Durch langsames Abkühlen (Stunden, Tage, je nach Größe des Metallstücks) haben die Kristalle ausreichend Zeit sich im Gefüge neu anzuordnen und ein stabiles Gefüge zu bilden. Nach der Abkühlung liegt ein Gefügezustand mit minimaler Energie vor, der Nahe am Optimum ist. Aufgrund des niedrigen Energiezustandes wird dies auch als Weichglühen bezeichnet was sich dadurch äußert dass das Metall weich ist 😉

Einführung
Simulated Annealing wurde inspiriert von der Wärmebehandlung von Metallen - dem sogenannten Weichglühen. Daher kommt auch die englische Bezeichnung dieses Verfahrens. Während andere Verfahren zum großen Teil in lokale Minima hängen bleiben können, ist es eine besondere Stärke dieses Algorithmus aus diesen wieder herauszufinden. Die am Ende gefunden Lösung stellt bei unendlicher Abkühlungszeit mit Sicherheit das globale Minimum (oder Maximum) dar, während bei endlicher Abkühlzeit dies nicht immer der Fall sein muss.

Beim langsamen Abkühlen eines Metalls ordnen sich die Kristalle so an dass sich im Gefüge ein Zustand mit möglichst niedriger Energie einstellt. Tritt allerdings eine zu schnelle Abkühlung auf, so haben die Kristalle nicht die Zeit den minimalen Energiezustand zu finden und das System bleibt in einem lokalen Minimum hängen. Ein niedrigerer Energieszustand ist in diesem Fall gleichbedeutend mit einem stabilen Endzustand, also einem stabilen Gefüge im Metall.

Nachfolgend eine grobe Skizzierung des Algorithmus (für ein genauere Betrachtung sei auf den separaten Teil "Einführung in den Simulated Annealing Algorithmus" verwiesen), wo bewusst das Analogon zur Metallurgie aufrecht erhalten wurde, um sich die Schritte bildlicher vorstellen zu können.1.Nimm eine zufällige Anfangslösung x aus dem Definitionsbereich des Problems. 1.Erwärme auf Glühtemperatur und somit die Energie des Systems. 1.Durch den hohen Energiezustand können sich die Kristalle um ihre Anfangsposition bewegen, d.h. es wird eine zufällige Verschiebung dx vorgenommen. 1.Ist der Energiezustand der neuen Position niedriger als der alte Energiezustand so wird dieser verwendet.
Das wäre bereits ein Algorithmus zur Suche nach einem Minimum, allerdings mit dem gravierenden Nachteil, dass der Algorithmus sich aus einem lokalen Minimum nicht mehr befreien kann, da er stets zu niedrigeren Energiewerten geht - bis er in einem stabilen Endzustand (= lokales Minimum) hängen bleibt.
Um dieses Problem zu lösen erlaubt SA auch einen Sprung zu höheren Energiezuständen. Diese werden allerdings nur mit einer bestimmten Wahrscheinlichkeit durchgeführt. Diese Wahrscheinlichkeit ist abhängig von der aktuellen Temperatur des Systems sowie der Energiedifferenz und folgt der Boltzmann-Gibbs-Verteilung.
Sinkt die Temperatur gegen 0 so wird die Wahrscheinlichkeit für einen Sprung auf ein höheres (Energie-) Niveau immer geringer. Für kleine Energieunterschiede ist die Wahrscheinlichkeit wiederum nahezu bei 100% und somit kann sich x in seiner näheren Umgebung noch frei bewegen, während für große Energieunterschiede ein Sprung zwar möglich, aber sehr unwahrscheinlich ist. 1.Senke die Temperatur und der Zyklus mit der zufälligen Verschiebung beginnt erneut. 1.Diese Schritte werden so lange durchgeführt bis die Endtemperatur erreicht ist. Dann beginnt der gesamte Zyklus erneut und wird wiederholt bis die Lösungen konvergieren oder sonst ein Abbruchkriterium erreicht wurde.

Verbildlicht kann sich das auch so vorgestellt werden:
Angenommen, man sucht in einer zweidimensionalen Landschaft den (global) tiefsten Punkt. Die Landschaft selbst besteht aus vielen unterschiedlich tiefen Dellen. Die einfache Suchstrategie (suche den nächsten tiefsten Punkt) entspricht dem Verhalten einer Kugel, welche in dieser Landschaft ausgesetzt wird. Sie rollt zum nächsten lokalen Minimum und bleibt dort. Bei der simulierten Auskühlung wird der Kugel immer wieder ein Stoß versetzt. Dieser ist stark genug, um die Kugel aus einer flachen Delle (lokales Minimum) zu entfernen, reicht aber nicht aus, um aus dem globalen Minimum zu fliehen.

Dieses Verfahren ist somit eine Heuristik. D.h. es ist ein Verfahren das oft sehr gut funktioniert, aber keine Garantie für eine optimale oder auch nur gute Lösung gibt. Simulated Annealing ist eine Entwicklung die aus dem Metropolisalgorithmus hervorgegangen ist.

Basisimplementierung
Der zugrunde liegende Algorithmus ist als abstrakte Basisklasse implementiert.

(Abbildung2: Basisklasse)


/* Copyright © Günther M. FOIDL 2009
 *
 * Lizenz: Common Development and Distribution License (CDDL)
 * [URL]http://www.opensource.org/licenses/cddl1.php[/URL]
 */
using System;

namespace gfoidl.ComputationalIntelligence.SimulatedAnnealing
{
	/// <summary>
	/// Stellt die Basisimplementierung des Algorithmus "Simulated Annealing"
	/// dar.
	/// </summary>
	/// <typeparam name="T">Der Typ des Problems.</typeparam>
	public abstract class SimulatedAnnealing<T>
	{
		#region Felder
		/// <summary>
		/// Zufallsgenerator.
		/// </summary>
		protected Random _rnd = new Random();
		#endregion
		//---------------------------------------------------------------------
		#region Eigenschaften
		/// <summary>
		/// Die Starttemperatur. Typsich: 500
		/// </summary>
		public double StartTemperature { get; set; }
		//---------------------------------------------------------------------
		private double _stopTemperatur = 0;
		/// <summary>
		/// Die Endtemperatur. Der Standardwert ist 0.
		/// </summary>
		public double StopTemperature
		{
			get { return _stopTemperatur; }
			set { _stopTemperatur = value; }
		}
		//---------------------------------------------------------------------
		/// <summary>
		/// Die Anzahl der Zyklen die pro Trainings-Iteration verwendet
		/// werden. Typsich: 100
		/// </summary>
		public int Cycles { get; set; }
		//---------------------------------------------------------------------
		/// <summary>
		/// Die aktuelle Energie (der Fehler).
		/// </summary>
		public double Energy { get; set; }
		//---------------------------------------------------------------------
		/// <summary>
		/// Die aktuelle Temperatur.
		/// </summary>
		/// <remarks>
		/// Die Temperatur entspricht der Wahrscheinlichkeit, mit der sich
		/// ein Zwischenergebnis der Optimierung auch verschlechtern darf.
		/// </remarks>
		public double Temperature { get; protected set; }
		//---------------------------------------------------------------------
		/// <summary>
		/// Der Lösungsvektor (beste Lösung).
		/// </summary>
		public abstract T[] Array { get; set; }
		#endregion
		//---------------------------------------------------------------------
		#region Methoden
		/// <summary>
		/// Diese Methode ermittelt den Fehler der aktuellen Lösung.
		/// </summary>
		/// <returns></returns>
		public abstract double DetermineEnergy();
		//---------------------------------------------------------------------
		/// <summary>
		/// Führt einen Zyklus Erwärmung-Abhkühlung durch.
		/// </summary>
		/// <remarks>
		/// Für eine ausführliche Beschreibung siehe
		/// <see cref="SimulatedAnnealing&lt;T>"/>
		/// </remarks>
		public void Anneal()
		{
			// Aktuellen Zustand ermitteln (Anfangszustand):
			double currentEnergy = this.DetermineEnergy();

			// Bisher bester Zustand speichern:
			T[] bestArray = this.GetArrayCopy();
			double bestEnergy = currentEnergy;

			// Erwärmen auf die Glühtemperatur:
			this.Temperature = this.StartTemperature;

			// (Langsame) Abkühlung durchführen:
			for (int i = 0; i < this.Cycles; i++)
			{
				// Aktuellen Zustand merken (für Rücksetzung wenn dieser
				// nicht übernommen wird):
				T[] currentArray = this.GetArrayCopy();

				// Kristalle zufällig anordnen - zu zufälligen Nachbarpunkt
				// gehen. D.h. der aktuelle Zustand wird verändert.
				this.Randomize();

				// Energie des neuen aktuellen Zustand ermitteln:
				double newEnergy = this.DetermineEnergy();

				// Hat sich ein neuer bester Zustand eingestellt?
				// Wenn ja diesen übernehmen:
				if (newEnergy < bestEnergy)
				{
					bestEnergy = newEnergy;
					bestArray = this.GetArrayCopy();
				}

				// Energieänderung im Vergleich zum vorigen Zustand:
				double deltaEnergy = newEnergy - currentEnergy;

				// Wahrscheinlichkeit der Übernahme des neuen Zustands.
				// Je höher die Temperatur desto wahrscheinlicher und
				// je kleiner die Energiedifferenz desto wahrscheinlicher.
				double probality = Math.Exp(-deltaEnergy / this.Temperature);

				// Zustand übernehmen:
				if (probality > _rnd.NextDouble())
				{
					currentEnergy = newEnergy;
					// Zustand ist bereits übernommen, denn die Verschiebung
					// wurde auf den aktuellen Zustand durchgeführt.
				}
				// Wenn der Zustand nicht übernommen wurde auf den vorigen
				// aktuellen Zustand zurücksetzen. Es wird nicht auf den
				// besten Zustand zurückgesetzt damit eine bessere
				// Exploration des Suchgebiets stattfindet.
				else
				{
					this.Array = currentArray;
				}

				// Temperatur senken:
				this.Temperature = ReduceTemperature(i);
			}

			// Beste Lösung verwenden. Die Originalversion verwendet die
			// aktuelle Lösung. Dies ist aber eine "klassische" Optimierung
			// die für jedes meta-heuristische Verfahren angewandt werden
			// kann. Entspricht zwar nicht ganz dem Analogon zur
			// Wärmebehandlung denn das physische System kann sich
			// keinen Zustand merken.
			this.Energy = bestEnergy;
			this.Array = bestArray;
		}
		//---------------------------------------------------------------------
		/// <summary>
		/// Senkt die Temperatur gemäß Abkühlungskurve.
		/// </summary>
		/// <param name="cycle">Die Zahl der vergangenen Zyklen.</param>
		/// <returns>Die abgesenkte Temperatur.</returns>
		protected virtual double ReduceTemperature(int cycle)
		{
			return this.StartTemperature * (1d - (double)(++cycle) / this.Cycles);
		}
		//---------------------------------------------------------------------
		/// <summary>
		/// Gibt eine Kopie des Vektors zurück.
		/// </summary>
		/// <returns>Kopie des Vektors.</returns>
		protected abstract T[] GetArrayCopy();
		//---------------------------------------------------------------------
		/// <summary>
		/// Ermittelt einen zufällige Position in der Nachbarschaft der
		/// aktuellen Position (Nachbar).
		/// </summary>
		/// <remarks>
		/// Wie im realen Vorbild kann diese Verschiebung von der Temperatur
		/// abhängig sein. Muss es aber nicht.
		/// <para>
		/// Für die Wahl des Nachbarn muss berücksichtigt werden dass nach ein
		/// paar Iterationen der aktuelle Zustand bereits eine niedrige Energie
		/// hat. Deshalb kann als grobe Regel angegeben werden, dass der
		/// Generator so arbeiten soll, dass ein Nachbar mit ähnlicher Energie
		/// wie der aktuelle Zustand gewählt wird.
		/// </para>
		/// </remarks>
		protected abstract void Randomize();
		#endregion
	}
}

Anwendung für die Minimierung einer reellen Funktion
Abbildung1 zeigt die Beispielanwendung.

Von der Basisimplementierung muss abgeleitet werden um eine Klasse für das konkrete Problem - hier die Minimierung einer reellen Funktion - zu erhalten. Es ist lediglich nötig die abstrakten Member der Basisklasse zu implementieren. Der restliche Algorithmus ist in der Basisklasse generisch implementiert.

Da die Problemdimension niedrig ist (2 dimensional) wird für die Wahl eines Nachbarn einfach eine neue Zufallsposition aus dem Definitionsbereich der Funkionen gewählt. Im Beispiel besitzen alle Funktion einen Definitionsbereich von [-3,3]x[-3,3].


/* Copyright © Günther M. FOIDL 2009
 *
 * Lizenz: Common Development and Distribution License (CDDL)
 * [URL]http://www.opensource.org/licenses/cddl1.php[/URL]
 */
using gfoidl.ComputationalIntelligence.SimulatedAnnealing;
using SimulatedAnnealing.FunctionMinimizing.Function;

namespace SimulatedAnnealing.FunctionMinimizing
{
	public sealed class FunctionMinimizingSimulatedAnnealing :
		SimulatedAnnealing<double>
	{
		private double[] _vector;
		private FunctionBase _function;
		//---------------------------------------------------------------------
		public override double[] Array
		{
			get { return _vector; }
			set { _vector = value; }
		}
		//---------------------------------------------------------------------
		public FunctionMinimizingSimulatedAnnealing(FunctionBase function)
		{
			_function = function;
			this.StartTemperature = 500;
			this.Cycles = 100;

			// [-3,3]...Definitionsbereich:
			double x = _rnd.NextDouble() * 6 - 3;
			double y = _rnd.NextDouble() * 6 - 3;

			_vector = new double[] { x, y };
		}
		//---------------------------------------------------------------------
		public override double DetermineEnergy()
		{
			return _function.Function(_vector);
		}
		//---------------------------------------------------------------------
		protected override double[] GetArrayCopy()
		{
			double[] tmp = new double[_vector.Length];
			_vector.CopyTo(tmp, 0);

			return tmp;
		}
		//---------------------------------------------------------------------
		protected override void Randomize()
		{
			// [-3,3]...Definitionsbereich:
			double x = _rnd.NextDouble() * 6 - 3;
			double y = _rnd.NextDouble() * 6 - 3;

			_vector = new double[] { x, y };
		}
	}
}

Im angehängten Beispielprojekt stehen ein paar Funktionen zur Auswahl bereit und die Anfangslösung - der Startpunkt - kann durch einen Mausklick bestimmt werden.

Das Beispiel meines Artikels "Particle swarm optimization for function optimization" ist mit dem angehängten Beispiel vergleichbar. Dabei ist auch ersichtlich dass SA für diese Art der Probleme wesentlich schneller ist, die Partikelschwärme jedoch interessanter zu beobachten sind 😉

Anwendung für das Problem des Handelsreisenden


(Abbildung3: Beispielanwendung für das Problem des Handelsreisenden. Die Lösung dauerte für die 100 Städte wenige Sekunden.)

Das Problem des Handelsreisenden - auch Traveling Salesman Problem (TSP) genannt - ist ein Vertreter für diskrete kombinatorische Probleme. Für nähere Informationen sei auf die Literatur verwiesen.

Während ein Brute-Force-Ansatz für hunderte Städte auf heutigen CPUs mehrere Jahre benötigen würde und viele andere Verfahren zur Lösung des TSP praktisch nur bis 50...100 Städte verwendbar sind schafft SA eine Tour mit 500 Städten in wenigen Sekunden.

Bevor der Code für die grundlegende Lösung erstellt wird müssen die Rahmenbedingunen zur Lösung aufgestellt werden:*Städte werden durch Vektoren repräsentiert *als Distanzmaß wird die euklidische Distanz verwendet *eine Route wird durch den Indexvektor der Städte dargestellt *SA wird auf diesen Indexvektor angewandt *für die Wahl werden zwei aufeinander folgende Städte vertauscht und die Route angepasst

Für eine etwas mathematischere Behandlung sei wiederum auf den separaten Teil verwiesen.

Für die Darstellung einer Stadt wird folgende Klasse verwendet:


using System;

namespace SimulatedAnnealing.TravelingSalesMan
{
	public class City
	{
		public int XPos { get; private set; }
		public int YPos { get; private set; }
		//---------------------------------------------------------------------
		public City(int x, int y)
		{
			this.XPos = x;
			this.YPos = y;
		}
		//---------------------------------------------------------------------
		public int Distance(City other)
		{
			return Distance(other.XPos, other.YPos);
		}
		//---------------------------------------------------------------------
		public int Distance(int x, int y)
		{
			int deltaX = this.XPos - x;
			int deltaY = this.YPos - y;

			return (int)Math.Sqrt(deltaX * deltaX + deltaY * deltaY);
		}
	}
}

Die konkrete Implementierung des SA für das TS-Problem:


/* Copyright © Günther M. FOIDL 2009
 *
 * Lizenz: Common Development and Distribution License (CDDL)
 * [URL]http://www.opensource.org/licenses/cddl1.php[/URL]
 */
using gfoidl.ComputationalIntelligence.SimulatedAnnealing;

namespace SimulatedAnnealing.TravelingSalesMan
{
	public sealed class TravelingSalesmanSimulatedAnnealing : SimulatedAnnealing<int>
	{
		private City[] _cities;
		//---------------------------------------------------------------------
		private int[] _path;
		public override int[] Array
		{
			get { return _path; }
			set { _path = value; }
		}
		//---------------------------------------------------------------------
		public TravelingSalesmanSimulatedAnnealing(
			City[] cities,
			double startTemp,
			int cycles)
		{
			this.Temperature = startTemp;
			this.StartTemperature = startTemp;
			this.Cycles = cycles;

			_cities = cities;
			_path = new int[_cities.Length];
		}
		//---------------------------------------------------------------------
		public override double DetermineEnergy()
		{
			double cost = 0;
			for (int i = 1; i < _cities.Length; i++)
			{
				double dist = _cities[_path[i - 1]].Distance(_cities[_path[i]]);
				cost += dist;
			}

			return cost;
		}
		//---------------------------------------------------------------------
		protected override int[] GetArrayCopy()
		{
			int[] result = new int[_path.Length];
			_path.CopyTo(result, 0);

			return result;
		}
		//---------------------------------------------------------------------
		/// <summary>
		/// Vertauscht Städte der Tour entsprechend der aktuellen Temperatur.
		/// </summary>
		protected override void Randomize()
		{
			for (int i = 0; i < this.Temperature; i++)
			{
				int index1 = _rnd.Next(_path.Length);
				int index2 = _rnd.Next(_path.Length);
				double d =
					Distance(index1, index1 + 1) +
					Distance(index2, index2 + 1) -
					Distance(index1, index2) -
					Distance(index1 + 1, index2 + 1);

				if (d > 0)
				{
					if (index2 < index1)
					{
						int tmp = index1;
						index1 = index2;
						index2 = tmp;
					}

					for (; index2 > index1; index2--)
					{
						int tmp = _path[index1 + 1];
						_path[index1 + 1] = _path[index2];
						_path[index2] = tmp;
						index1++;
					}
				}
			}
		}
		//---------------------------------------------------------------------
		/// <summary>
		/// Ermittelt die Distanz zwischen 2 Städten.
		/// </summary>
		/// <param name="i">Index der 1. Stadt.</param>
		/// <param name="j">Index der 2. Stadt.</param>
		/// <returns>Distanz zwischen beiden Städten.</returns>
		private double Distance(int i, int j)
		{
			int c1 = _path[i % _path.Length];
			int c2 = _path[j % _path.Length];

			return _cities[c1].Distance(_cities[c2]);
		}
	}
}

Die Beispielanwendung hierfür ist ebenfalls angehängt.

Verwandte Verfahren
Andere Verfahren der Heuristik um ähnliche Problem zu lösen wären:*genetische Algorithmen *Partikelschwarmoptimierung *Ameisenalgorithmus *Bienenalgorithmus *und Abwandlungen der genannten

Siehe auch meine Blogbeiträge über Computation Intelligence und natürlich den separaten Teil 😉

Viel Spass!

mfG Gü

Schlagwörter: Optimierung, Minimum, Minimumsuche, Heuristik, Simulated Annealing, Traveling Salesman, Problem des Handelsreisenden

Stellt fachliche Fragen bitte im Forum, damit von den Antworten alle profitieren. Daher beantworte ich solche Fragen nicht per PM.

"Alle sagten, das geht nicht! Dann kam einer, der wusste das nicht - und hat's gemacht!"