Simulated annealing

Simulated annealing (SA) is een generiek, probabilistisch heuristiek optimalisatiealgoritme gebruikt om een benadering van het globale optimum van een gegeven functie in een grote zoekruimte te vinden. Het is onafhankelijk van elkaar uitgevonden door S. Kirkpatrick, C. D. Gelatt en M. P. Vecchi in 1983 en door V. Cerny in 1985.

De naam en inspiratie komen van een Engelse term, 'annealing' (uitgloeien), binnen de metaalbewerking. Het betreft een techniek waarbij metaal verhit wordt en daarna gecontroleerd afgekoeld om de grootte van de kristallen binnen het materiaal te vergroten en daarmee het aantal defecten te verkleinen.

Achtergronden

Simulated annealing is een optimalisatiealgoritme uit de combinatoriek. Het is een verbetering van de familie van lokale zoekalgoritmen, gebaseerd op het annealing-proces dat gebruikt wordt binnen de metallurgie.

Combinatoriek en Local Search

De combinatoriek is een deelveld van de wiskunde dat gaat over de oplossing van bepaalde klassen van problemen, de zogeheten combinatorische optimalisatieproblemen.

Een combinatorisch optimaliseringsprobleem is een probleem dat gaat over het vinden van de beste combinatie van een aantal variabelen of variabele-waarden. Een dergelijk probleem wordt gedefinieerd op een zogeheten toestandsruimte -- een voorbeeld hiervan is het handelsreizigersprobleem, waarbij de toestandsruimte bestaat uit alle mogelijke combinaties van steden dusdanig dat die combinaties een pad vormen.

Bij ieder element van de toestandsruimte hoort een zogeheten waardering. Deze waarderingen geven aan wat de oplossing "waard is" of wat hij "kost"; de waarderingen maken het mogelijk om oplossingen kwalitatief met elkaar te vergelijken. Een combinatorisch optimaliseringsprobleem komt er in het kort op neer om uit de toestandsruimte de oplossing te vinden met de hoogste waarde (dit noemen we een maximalisatieprobleem) hetzij de laagste kosten (minimalisatieprobleem).

Een combinatorisch optimaliseringsprobleem zoals hierboven beschreven wordt in de wiskunde wel genoteerd als het tupel

(S,f)

Met $S$ de toestandsruimte (state space) en $f:S\rightarrow \mathbb {R}$ een functie die voor iedere oplossing de waardering geeft. Combinatorische optimaliseringsproblemen zoeken voor een gegeven tupel dus naar $i_{opt}\in S$ met $\forall _{i\in S}:f(i_{opt})\leq f(i)$ (in het geval van minimalisatie; voor maximalisatie is het net andersom). Deze $i_{opt}$ heet overigens een (het, als er maar een is) globaal optimum van het probleem.

Er zijn twee soorten algoritmen die een dergelijk probleem oplossen:

Exacte algoritmen: Deze algoritmen vinden gegarandeerd een globaal optimum.
Benaderingsalgoritmen of heuristische algoritmen: Deze algoritmen doorzoeken op een intelligente manier de toestandsruimte vanuit een bepaalde begintoestand. Ze vinden altijd wel een oplossing, maar niet noodzakelijk het globale optimum.

Het nut van de heuristische aanpak van dergelijke problemen bestaat erin dat niet voor ieder probleem een exact algoritme bekend is dat het probleem binnen redelijke tijd oplost. Voor ieder probleem is het mogelijk een exacte oplossing te bepalen, maar soms duurt dat onmogelijk lang (zie ook complexiteitstheorie). Benaderingsalgoritmen bieden in deze gevallen vaak een werkbaar alternatief, met het risico dat de beste oplossing niet gevonden wordt.

Een zeer bekende familie van heuristische algoritmes is de familie van Local Search algoritmes. Deze algoritmes zijn populair omdat ze redelijke resultaten opleveren en algemeen toepasbaar zijn -- ze zijn niet toegespitst op een specifiek probleem.

Local Search algoritmes beginnen bij een gegeven toestand i uit de toestandsruimte $S$ en bekijken dan de "buurtoestanden" van i. De buurtoestanden van i zijn de toestanden die op de een of andere manier "dicht bij i in de buurt liggen" -- in het geval van het handelsreizigersprobleem bijvoorbeeld, gegeven een bepaalde route, de routes die je kunt vinden door in die eerste route net twee paren steden andersom aan elkaar te knopen. Local Search algoritmen bekijken deze buurruimte en kiezen daaruit de oplossing met de beste waardering. Dat wordt volgende toestand. Daarna beginnen ze opnieuw met zoeken naar de volgende, betere toestand enzovoorts, enzovoorts, totdat er in de lokale buurruimte geen betere toestand is. Dit is een lokale optimum en mogelijk ook het globale optimum.

Het probleem met Local Search algoritmen is dat hun werking eigenlijk neerkomt op het geblinddoekt rondlopen door een heuvellandschap, waarbij je onder de rand van je blinddoek net één stap om je heen kunt kijken. Als je begint aan de hoge rand van een diepe vallei en steeds die stap zet die het hardst naar beneden gaat, kom je aan het eind van de rit op het laagste punt uit. Maar in een ingewikkelder landschap kan aan het eind van je wandeling ook blijken dat je in een heel klein kuiltje boven op een hoge berg staat -- lokaal minimum, maar zeker niet het best bereikbare.

Om deze reden zijn er vele soorten lokale zoekalgoritmen. Ieder heeft een eigen strategie om het landschap te doorlopen en wellicht ook een strategie om een goed beginpunt te kiezen. Simulated annealing is een voorbeeld van een dergelijk algoritme, dat probeert het resultaat te verbeteren door een betere strategie om de volgende stap te kiezen. In tegenstelling tot wat hierboven beschreven is, geldt namelijk dat simulated annealing niet altijd de stap zet die omlaag gaat -- soms gaat simulated annealing ook wel weer eens omhoog, om te kijken of er verderop een dieper dal ligt.

Annealing

Simulated annealing is gebaseerd op het annealing-proces uit de metallurgie.

Annealing wordt toegepast in de metallurgie om een metalen object zo te maken dat er geen scheurtjes of fouten in het materiaal zitten -- dat wil zeggen, zo te maken dat de moleculen van het metaal perfect in een raamwerk zitten zonder gaten in het raamwerk.

Om dit te bereiken wordt het object in een verwarmingsbad gedompeld en precies op smelttemperatuur gebracht. Daarna wordt het object heel langzaam, in precies gecontroleerde stappen, weer afgekoeld. Dit heeft het effect dat het energieniveau van alle moleculen in het object ongeveer gelijk blijft gedurende het afkoelen -- thermisch evenwicht. Door dit evenwicht te bewaren tijdens het afkoelen ontstaat er nergens in het object een "energiepiek", een plek met moleculen die hogere energie hebben dan omliggende moleculen en dus ten opzichte van het rooster extra gaan bewegen. Door het object in thermische balans af te laten koelen bereiken alle moleculen tegelijkertijd (zo ongeveer) hun laagste energieniveau en wordt een object verkregen waarin alle moleculen netjes op hun plek in het rooster "tot rust" gekomen zijn.

Rond 1953 werden de eerste computersimulaties van annealing ontwikkeld. Hierbij werd van Monte Carlo technieken gebruikgemaakt om het gedrag van moleculen in het verhitte object te bepalen. De toegepaste analyse van het probleem resulteerde in het Metropolis algoritme voor het simuleren van het annealing proces.

Het Metropolis algoritme beschouwt toestanden die bestaan uit de posities van moleculen ten opzichte van het raster. Iedere toestand wordt gewaardeerd met de totale energie van die toestand. Het algoritme draait totdat het systeem een toestand bereikt heeft waarin geen verbetering mogelijk is, de toestand met de laagste energie.

Daarmee is echter de kous niet af. Het Metropolis algoritme vangt namelijk een bijzondere eigenschap van het annealing-proces: niet iedere toestandsovergang bij het afkoelen verlaagt de totale energie. Soms gaat een molecuul (of meerdere) toch nog "dansen" in het rooster, waardoor de volgende toestand niet een lagere energie heeft dan de vorige. Om dit te vangen, bedachten de makers van Metropolis de volgende regel voor toestandsovergang:

Gegeven een begintoestand, kiezen we een toestand uit de mogelijke buurtoestanden.
Is de waardering van de gekozen toestand lager dan de waardering van de huidige toestand, dan accepteren we hem.
Is de waardering van de gekozen toestand hoger dan de waardering van de huidige toestand, dan accepteren we hem met een zekere waarschijnlijkheid, of wijzen hem af met de bijbehorende complementaire waarschijnlijkheid.
Is de gekozen toestand geaccepteerd, ga dan over naar die toestand en begin opnieuw.
Is de gekozen toestand niet geaccepteerd, dan
1. Als we nog niet de hele buurruimte hebben gezien, probeer dan een andere buurtoestand.
2. Hebben we geen buurtoestanden meer over, stop dan.

Het gedeelte over acceptatie en afwijzing van de overgang van toestand i naar toestand j vingen ze in de volgende uitdrukking:

\mathbb {P} _{T}(X_{i}=j)={\begin{cases}E_{j}\leq E_{i}\rightarrow 1\\E_{i}<E_{j}\rightarrow e^{\frac {E_{i}-E_{j}}{k_{B}\cdot T}}\end{cases}}

Hierin is $X_{i}$ een stochast die staat voor "de opvolger van toestand i". Dat wil zeggen, het bovenstaande leest als "de kans dat de opvolgende toestand van i, toestand j is, is ...". Hierin is $E_{n}$ de energie (of waardering) van toestand n, $T$ de temperatuur van het object en $k_{B}$ de constante van Boltzmann.

Deze simulatie demonstreert goed het fenomeen uit de annealing dat thermisch evenwicht bij iedere temperatuur gekarakteriseerd wordt door de Boltzmann-verdeling -- dit is een statistische verdeling die van iedere, mogelijke toestand zijn waarschijnlijkheid aangeeft bij een gegeven temperatuur T. De kans dat een object bij temperatuur T een gegeven toestand aan zal nemen tijdens het annealen is

\mathbb {P} _{T}(X_{i}=j)={\frac {e^{\frac {-E_{i}}{k_{B}\cdot T}}}{\sum _{j}e^{\frac {-E_{j}}{k_{B}\cdot T}}}}

De noemer in deze term heet de partitie functie en het domein van deze functie is de verzameling van alle, mogelijke toestanden in het systeem. Een dergelijke functie is ook van belang voor de analyse van het gedrag van het simulated annealing algoritme.

Simulated annealing

Het algoritme

Uit het bovenstaande is het simulated annealing algoritme bijzonder eenvoudig te herleiden. Net als bij het echte annealen hebben we een bestaande toestand, een aantal mogelijke toestandsovergangen en een waardering voor iedere overgang. Vanuit de gegeven toestand kiezen we een buurtoestand. Heeft deze een betere waardering, dan accepteren we hem. Zo niet, dan accepteren we hem met een bepaalde waarschijnlijkheid. Al met al komt de regel voor acceptatie neer op

\mathbb {P} _{c}(accept(j))=\left\{{\begin{matrix}f(j)\leq f(i)&1\\f(i)<f(j)&e^{\frac {f(i)-f(j)}{c}}\end{matrix}}\right.

Hierin is c een controleparameter die de rol van de temperatuur overneemt. Het hele algoritme voert de selectie van buren en overgang van toestanden een aantal keren uit. Bij iedere iteratie wordt de waarde van c verlaagd. We zullen zien dat hierdoor het algoritme steeds "lokaler" wordt, wat wil zeggen dat het algoritme steeds minder grote verslechteringen van de waardering accepteert bij het kiezen van een buur en dus steeds meer gaat neigen naar het kiezen van betere waarderingen naarmate de tijd vordert.

Het gehele algoritme voor een probleem $(S,f)$ ziet er als volgt uit:

 proc Simulated_Annealing(S, f) = 
 |[
  var i, j : element van S;
    c : Real;
    k, l, L : Integer;
  |
  i := <element van S>;
  k, c, L := 0, infinity, <redelijke beginwaarde>;
  
  do NOT(stopcriterium) ->
   l := 1;
   do l <= L ->
    j := <kies een element uit de buren van i in S>;
    if f(j) <= f(i) -> i := j
    [] f(i) < f(j) -> if e^((f(i) - f(j)) / c) > random[0, 1) -> i := j fi
    fi;
   od;
   k := k + 1;
   c := <bereken nieuwe c>;
   L := <bereken nieuwe L bij nieuwe c>;
  od
 ]|

Merk de extra parameter L op in de bovenstaande code. Voor iedere waarde van de controleparameter c worden een aantal toestandsovergangen gesimuleerd, net als gebeurt in het echte annealing bij een bepaalde temperatuur. Bij annealing gebeurt dit om te simuleren wat er in een fysiek object gebeurt, namelijk dat de moleculen zich gedurende afkoelen van de ene temperatuur naar de volgende in een aantal stappen (toestanden) naar een lager energieniveau toe werken.

Door dit getrouw na te spelen in simulated annealing, verhogen we de kans om het globale optimum te vinden. Het komt erop neer dat we niet klakkeloos per controleparameter één toestand accepteren en daarmee verdergaan, maar eerst rond kunnen kijken of er een beter optie is. Soms accepteren we zelfs een slechtere optie dan mogelijk is, met het idee dat we uiteindelijk beter uitkomen. Dat we -- naarmate c een lagere waarde krijgt -- het globale optimum waarschijnlijk beter benaderen, zullen we verderop zien.

Voor nu merken we echter in verfijndere mate op wat we al eerder aanstipten: de directe invloed van c op de berekening is dat simulated annealing het weliswaar toestaat om toestandsovergangen te hebben waarbij de waardering slechter wordt, maar wel steeds minder naarmate de tijd vordert. Dat wil zeggen, naarmate c een lagere waarde krijgt (en we gaan er meestal van uit dat c monotoon daalt), wordt de kans op acceptatie van een overgang naar een toestand met een hogere waardering kleiner. Dit valt direct te zien aan de term

e^{\frac {f(i)-f(j)}{c}}

Als c een grote waarde heeft (zoals initieel), is de waarde van deze term ongeveer 1. Dat betekent dat een overgang naar een toestand met een hogere (slechtere) waardering (f(j) is tenslotte hoger dan f(i)) ongeveer 1 is. Naarmate c kleiner wordt, wordt de uitkomst meer en meer bepaald door de term $f(i)-f(j)$ , die altijd ten hoogste nul is (wederom: het gaat om overgang naar een toestand met hogere waardering). Hoe groter het verschil, hoe lager de waarde van de term (en dus de kans op acceptatie van de overgang). De kans op acceptatie is dan dus hoger naarmate het verschil $f(i)-f(j)$ kleiner is.

Waarschijnlijkheid van het vinden van een oplossing

Bij de bespreking van echte annealing en het Metropolis algoritme hebben we opgemerkt dat het object voor iedere temperatuur een zogeheten thermisch evenwicht bereikt. Dit evenwicht werd gekarakteriseerd door een verdeling van waarschijnlijkheden over toestanden -- voor iedere toestand was er een bepaalde waarschijnlijkheid dat het systeem die toestand aan zou nemen.

Simulated annealing gedraagt zich net zo. Bij iedere waarde van de parameter c worden een aantal toestandsovergangen in de buurruimte bekeken -- dit aantal is de parameter L. Het is aantoonbaar dat, als L maar groot genoeg is voor iedere waarde van c, dat zich voor iedere waarde van c een "thermische balans" instelt. Dat wil zeggen, voor iedere waarde van c en toestand $i\in S$ geldt dat

\mathbb {P} _{c}(X=i)\ {\hat {=}}\ {\frac {e^{\frac {-f(i)}{c}}}{\sum _{j\in S}e^{\frac {-f(j)}{c}}}}

De reden dat deze balans van belang is, is dat ze garandeert dat het simulated annealing-algoritme asymptotisch nadert tot de verzameling van globaal optimale toestanden in de toestandsruimte (deze verzameling noemen we $S_{opt}$ ). Er geldt immers

{\begin{matrix}&{\begin{matrix}\lim \\c\downarrow 0\end{matrix}}{\frac {e^{\frac {-f(i)}{c}}}{\sum _{j\in S}e^{\frac {-f(j)}{c}}}}\\=&\{{\mbox{inbrengen van de waardering van het globale optimum}}\}\\&{\begin{matrix}\lim \\c\downarrow 0\end{matrix}}{\frac {e^{\frac {f_{opt}}{c}}}{e^{\frac {f_{opt}}{c}}}}{\frac {e^{\frac {-f(i)}{c}}}{\sum _{j\in S}e^{\frac {-f(j)}{c}}}}\\=&\\&{\begin{matrix}\lim \\c\downarrow 0\end{matrix}}{\frac {e^{\frac {f_{opt}-f(i)}{c}}}{\sum _{j\in S}e^{\frac {f_{opt}-f(j)}{c}}}}\\=&\{{\mbox{Deze kans is samengesteld uit de kans wanneer toestand i een optimum is en wanneer niet;}}\\&{\mbox{ dit geven we aan met de chi-functie}}\}\\&{\begin{matrix}\lim \\c\downarrow 0\end{matrix}}{\frac {e^{\frac {f_{opt}-f(i)}{c}}}{\sum _{j\in S}e^{\frac {f_{opt}-f(j)}{c}}}}\mathrm {X} _{(S_{opt})}(i)+{\frac {e^{\frac {f_{opt}-f(i)}{c}}}{\sum _{j\in S}e^{\frac {f_{opt}-f(j)}{c}}}}\mathrm {X} _{(S\setminus S_{opt})}(i)\\=&\{{\mbox{Als i niet een optimum is, dan is }}f_{opt}-f(i){\mbox{ niet 0;}}\\&{\mbox{gebruik nu dat }}{\begin{matrix}\lim \\x\downarrow 0\end{matrix}}e^{\frac {a}{x}}\\&{\mbox{ gelijk is aan 0 als a kleiner is dan 0 en 1 als a = 0 }}\}\\&{\begin{matrix}\lim \\c\downarrow 0\end{matrix}}{\frac {e^{\frac {f_{opt}-f(i)}{c}}}{\sum _{j\in S}e^{\frac {f_{opt}-f(j)}{c}}}}\mathrm {X} _{(S_{opt})}(i)+0\\=&\{{\mbox{De taylor-ontwikkeling van }}e^{x}{\mbox{ is }}1+x+{\frac {x^{2}}{2!}}+\ldots \\&{\mbox{Alleen de eerste term van deze expansie is niet 0, want }}f_{opt}-f(i)\\&{\mbox{ is 0 als i een optimum is}}\}\\&{\begin{matrix}\lim \\c\downarrow 0\end{matrix}}{\frac {1}{\sum _{j\in S}e^{\frac {f_{opt}-f(j)}{c}}}}\mathrm {X} _{(S_{opt})}(i)\\=&\{{\mbox{Op de noemer in de vorige term kunnen we ook weer de limiet van hiervoor toepassen.}}\\&{\mbox{ De sommatie wordt dan de sommatie van 1 voor optima en 0 voor niet-optima;}}\\&{\mbox{ dat is het aantal optima in S}}\}\\&{\frac {1}{|S_{opt}|}}\mathrm {X} _{(S_{opt})}(i)\end{matrix}}

Dit wil zeggen: naarmate c kleiner wordt, naderen we een verdeling waarbij de kans dat we het systeem in een gegeven toestand aantreffen uniform verdeeld is over de optimale toestanden en 0 is voor niet-optimale toestanden. Naarmate we c verlagen, wordt de benadering van het optimum steeds beter.

In het bovenstaande hebben we de Χ-functie gebruikt om de kans te verdelen over wel- en niet-optimale toestanden uit de toestandsruimte. Deze functie is als volgt gedefinieerd: