L'algorithme du recuit simulé

Le phénomène physique

Une première approche

Lorsque nous avons abordé l'étude du pendule simple, nous avons constaté que celui-ci avait tendance à toujours revenir dans la position où son énergie potentielle était minimum, c'est à dire à son altitude minimum. C'est une règle fondamentale en physique, un système tend vers la configuration où son énergie totale est la plus faible possible.
Dans le cas du pendule, c'est assez simple : il n'y a qu'un point de sa trajectoire où son énergie est minimum : le point bas de sa trajectoire, avec une vitesse nulle évidemment...
Mais imaginons maintenant un cas un peu plus compliqué, toujours dans le plan. Considérons une balle posée au sommet d'une colline dont le profil ressemble à celui du schéma ci-dessous. Son équilibre est instable et au moindre coup de vent, la balle va glisser sur le long des flancs pour rejoindre la position d'énergie potentielle minimum, le pied de la colline...

Notre colline présente quelques faux-plats ! La balle peut donc se trouver bloquée dans une cuvette et ne pas atteindre le bas de la colline. Pour la faire sortir de sa cuvette et lui permettre de continuer sa descente, il faut lui fournir un peu d'énergie.

Cette situation est assez courante. Il n'est pas toujours facile d'amener un système vers son état d'énergie minimum sans le coincer dans des états locaux dans lesquels son énergie n'est pas minimum.

Les opérations de recuit en métallurgie

Considérons un morceau de métal, pur pour nous simplifier la vie, que l'on chauffe. Les atomes de métal dans ce morceau, sont organisés en un réseau dont le motif varie en fonction du métal. Chauffer le métal revient à faire vibrer plus ou moins fort les atomes autour de leur position moyenne. Le réseau se déforme et l'ensemble des atomes acquièrent de l'énergie. Lorsque le morceau de métal refroidit, on s'attend à ce que les atomes reviennent sagement à leur position d'origine et que tout revienne dans l'ordre, c'est à dire que le morceau de métal revienne à son énergie initiale.

Oui, mais cela ne se passe pas exactement comme ça... Si le refroidissement est brutal, le réseau métallique peut se bloquer dans des configurations qui ne sont pas la configuration initiale, que l'on posera comme celle d'énergie minimum. Des tensions se crééent, qui ont généralement un effet néfaste sur les caractéristiques mécaniques du métal, sauf cas particuliers où l'on cherche cet effet, comme une trempe par exemple. On constate aussi le phénomène dans les verres. Je me souviens de quelques déboires dus à ces tensions lorsque je tentais dans ma jeunesse de souder des tubes de verre...

Que faire pour éviter ce phénomène ? Il faut refroidir très lentement le métal et si des tensions sont apparues, on le réchauffe puis le refroidi très lentement : c'est l'opération de recuit. Pourquoi refroidir très lentement ? Le but est d'obtenir qu'à une température donnée lors du refroidissement, le morceau de métal soit à l'équilibre thermique correspondant à cette température. Parce que cet équilibre thermique correspond à l'état d'énergie minimum à la température en question. On laisse donc le métal évoluer d'un état d'équilibre vers un autre état d'équilibre, avec une énergie toujours plus faible. Ainsi, on obtient un réseau métallique pratiquement exempt de défauts.

Généralisation

Si l'on généralise, l'opération de recuit consiste à mettre en oeuvre tous les moyens nécessaires pour atteindre un minimum qui ne soit pas un minimum local, l'équivalent de nos cuvettes sur les flancs de notre colline. En mathématique, on appelle ce genre d'opération une optimisation.
La recherche du minimum global d'une fonction dans un domaine donné n'est pas chose simple, surtout dans le cas de fonctions à plusieurs variables. Le plus grand risque est de confondre un minimum local sur la domaine avec le minimum global du domaine, en quelque sorte rester bloqué au fond d'une cuvette sur le flanc de la colline.

L'algorithme du recuit simulé

Des méthodes pour rechercher le minimum global d'une fonction

Au lycée, quand vous cherchez l'extremum minimum d'une fonction sur un domaine donné (il faut toujours préciser le domaine !), vous calculez la valeur (ou les valeurs) pour laquelle (lesquelles) la dérivée s'annule et change de signe (passe de négatif à positif). Si cette valeur est unique, pas de problème, vous avez un minimum global sur le domaine. Mais s'il existe plusieurs valeurs pour lesquelles cette dérivée s'annule ? Comment savoir quel est le minimum global dans l'ensemble de ces valeurs ? Le calcul...

Le principe reste le même dans le cas des fonctions à plusieurs variables, en utilisant le gradient de la fonction, qui n'est d'autre que sa dérivée selon une direction donnée (voir ici par exemple). Mais le problème reste le même : comment savoir si le minimum calculé est un minimuml local ou global sur le domaine ?

La méthode du gradient a été largement perfectionnée pour éviter de rester coincé dans un minimum local, par exemple avec la méthode du gradient conjugué. Elle a donné naissance à toute une famille d'algorithmes déterministes pour optimiser une fonction, c'est à dire trouver ses minimums (ou ses maximums d'ailleurs...). Mais on arrive vite aux limites de la méthode, en particulier pour les fonctions à grand nombre de variables.

Il existe une autre famille d'algorithmes pour faire le boulot, que l'on qualifie d'algorithmes stochastiques, parce qu'ils font intervenir le hasard. On y trouve les algorithmes génétiques et aussi les algorithmes de recuit simulé (simulated annealing). C'est à ces derniers que nous allons nous intéresser.

L'algorithme du recuit simulé

Les principes

Comme vous pouvez vous en douter, cet algorithme est fortement inspiré du phénomène physique du recuit ! Pour construire cet algorithme, nous allons imaginer que la fonction à optimiser soit assimilée à l'énergie d'un système physique fictif. A l'initialisation, notre système physique possède une énergie \( E_0\) et est porté à une température \( T_0\). L'énergie initiale du système fictif correspond à la valeur \( x_0\) que l'on applique à la fonction pour débuter la recherche du minimum. La température est un paramètre de l'algorithme.

Puis nous appliquons les lois de la physique statistique à notre système physique fictif. Elles nous apprennent que l'état d'énergie de notre système physique fluctue autour d'une énergie moyenne dépendante de la température (non nulle). Vu de notre fonction, cela veut dire que nous faisons fluctuer notre variable x dans un intervalle \( [x - \delta x, x + \delta x] \) avec \( \delta x \) une valeur aléatoire. Nous procédons à cette opération jusqu'à l'équilibre thermique pour la température donnée, équilibre désigné par un paramètre de notre algorithme. Une fois l'équilibre atteint, on diminue un peu la température et l'on recommence la recherche de l'équilibre thermique avec la nouvelle température. On diminue ainsi la température jusqu'à atteindre une température basse limite, qui est un autre paramètre de l'algorithme. On suppose, et cela se démontre, que l'on atteint ainsi l'état d'énergie minimum du système physique fictif et donc le minimum de notre fonction.

Comment déterminer les fluctuations aléatoires de l'énergie de notre système ? En appliquant la méthode de Monte-Carlo à l'aide de l'algorithme de Nicholas Metropolis ! Oui mais encore ... La physique statistique, Boltzmann plus précisément, nous dit que la probabilité de trouver un système dans un état d'énergie s est proportionnelle à \( e^{-\dfrac{E_s}{T}} \). L'algorithme de Metropolis propose une méthode de fluctuation aléatoire basée sur cette loi:

soit \( x_i \) un état d'énergie \( E_i \), dans notre système
je fais varier \( x_i \) en lui ajoutant une grandeur aléatoire comprise entre 0 et 1, j'obtiens \( x_{i+1} \) d'énergie \( E_{i+1} \)
si \( E_{i+1} \le E_i \) alors \( x_{i+1} \) est le nouvel état de mon système, parce que son énergie diminue
sinon, je ne me bloque pas pour ne pas rester dans un minimum local éventuel, je décide que \( x_{i+1} \) devient mon nouvel état avec une probabilité égale à \( e^{-\dfrac{E_{i+1} - E_i}{k_B T}}\).

Il y a beaucoup de paramètres dans cet algorithme, ce qui en fait sa force et sa complexité de mise en oeuvre. Il faut en particulier choisir la "température" initiale du système et sa loi de décroissance. Il y a quelques règles mais à vrai dire, ces choix tiennent beaucoup de l'expérience... En résumé:

la "température" initiale doit être choisie suffisamment élevée pour que l'équilibre thermique soit rapidement atteint. Cela signifie qu'elle doit être grande par rapport aux variations aléatoires d'énergie du système lors de la recherche de l'équilibre.
la loi de décroissance de la température est un compromis entre la recherche du meilleur résultat possible et le temps de calcul. La théorie (le théorème de Hajek) nous dit que cette loi devrait être en log inverse (1/log(n) si n est le nombre de pas de varaition de température). Mais cela génère des temps de calcul épouvantables ! Alors on choisit généralement une loi de décroissance exponentionnelle du genre \( T = T_0 e^{-\dfrac{t}{\tau}}\), avec \( \tau \) suffisamment grand pour que la décroissance soit lente. A régler dans le programme...
le choix de la température finale : là aussi, c'est un compromis entre la qualité du résultat et le temps de calcul. Et cela dépend de la nature du problème d'optimisation ! Question d'expérience et d'essais multiples...

L'algorithme

Voici le pseudocode de l'algorithme du recuit simulé appliqué à l'optimisation d'une fonction f(x1,x2,...,xn) , en utilisant l'algorithme de Métropolis pour le traitement des fluctuations:

Définir la fonction f, la valeur initiale des variables xi et la valeur initiale T0 de la température
Tant que Ti < Tfinale
- Atteindre l'équilibre du système à Ti
- Calculer l'énergie moyenne du système pour Ti.
- Traiter l'écart d'énergie (Ei - Ei-1) selon l'algorithme de Metropolis
- Sinon, on diminue la température selon la formule choisie.

L'algorithme de recuit simulé est très général et l'on peut l'appliquer dans une grande famille de problèmes d'optimisation, par exemple comment implanter au mieux des composants sur une carte électronique ou comment définir au mieux l'enchainement des tâches dans un process de fabrication. Nous allons en voir une application très classique...

Une application très classique : le problème du voyageur de commerce

Le problème

Un voyageur de commerce (terme ancestral pour désigner un commercial...) doit visiter ses clients qui sont répartis dans un certain nombre de villes du pays. Son patron étant regardant sur les notes de frais, et lui étant regardant sur la fatigue des voyages, il cherche à minimiser la distance à parcourir. Il cherche donc comment organiser sa tournée en ne passant qu'une fois dans chaque ville, en revenant à son point de départ et en minimisant la distance parcourue.

Fixons les idées. Soit E l'ensemble de n villes numérotées de 0 à n-1. J'appelle \( \delta(v_{n+1},v_i)\) la distance connue entre deux villes voisines. La distance totale d que va parcourir le voyageur de commerce en faisant son tour sera donc égale à \( d = \sum_{i=0}^{n-1} \delta(v_{n+1},v_i) \). Le problème du voyageur de commerce consiste "simplement" à optimiser d, c'est à dire à chercher sa valeur minimum.

Pourquoi utiliser l'algorithme du recuit simulé

Raisonnons sur un ensemble E avec peu d'élements, par exemple 3 villes, numérotées 1, 2 et 3. Un moyen de résoudre le problème est d'identifier tous les chemins possibles et d'en calculer la distance. La solution est évidemment le chemin dont la distance est la plus faible.

Les chemins possibles, sachant qu'il faut revenir au point de départ, sont : (1,2,3,1), (1,3,2,1), (2,3,1,2), (2,1,3,2), (3,2,1,3) et (3,1,2,3). Il y a 6 combinaisons de chemins possibles. En fait, le nombre de possibilités s'élève à n! (n factorielle). Pour n=3, cela fait 6 possibilités. Pour 10 villes, cela fait 3,6 millions de combinaisons possibles et pour 20 villes 2,4 10^18 combinaisons possibles ! La méthode exhaustive n'est pas envisageable ! Le problème du voyageur de commerce (PVC) est de complexité non polynomiale, un problème que l'on appelle NP-complet.

C'est pour cela que nous allons utiliser un algorithme qui ne nous donnera peut-être pas le résultat exact mais une très bonne estimation, d'autant meilleure que nous choisirons de privilégier la qualité du résultat sur le temps de calcul...

Adapter le recuit simulé au PVC

La fonction à optimiser est toute trouvée : il s'agit de la distance totale à parcourir par le voyageur, qui est la somme des distances entre chacune des paires de villes voisines.

L'état du système est caractérisé par le trajet entre les différentes villes et la distance totale associée qui est représentée par son énergie. Pour modifier l'état du système, il suffit de permuter l'ordre des villes.

L'algorithme de Metropolis sera adapté : l'énergie du système est la distance d à minimiser, la variation d'énergie étant la différence de distance entre deux villes.

La température est un paramètre de contrôle sans signification fonctionnelle particulière. Mais le choix des températures initiale et finale contraindra l'efficacité de l'algorithme.

Le script Python

En fonction du nombre de villes, les calculs seront lourds. Il convient donc dans ce genre de scripts d'optimiser sa programmation. En Python, on utilisera autant que possible les instructions vectorielles plutôt que des boucles. Dans mon script, plusieurs améliorations sont possibles. Par exemple, je calcule l'énergie totale sur chaque trajet en calculant la distance entre chaque couple de villes. Même si j'utilise des instructions vectorielles et toutes les astuces Python, il serait plus rapide de créer une matrice contenant les distances entre villes calculées une bonne fois pour toutes. On peut aussi affiner le calcul de fluctuation...

Je décris ci-dessous les différentes parties significatives du code, sans m'attarder sur les initialisations diverses et les fonctions d'affichage.

Définition du problème

Les coordonnées des villes (x,y) sont tirées aléatoirement selon la loi uniforme, x et y sont compris entre 0 et 1.

x = random.uniform(size=N)

y = random.uniform(size=N)

Le trajet initial est simplement l'indexation du vecteur des villes selon leurs coordonnées aléatoires initialisées ci-dessus.

trajet = arange(N)

Je sauvegarde le trajet initial pour pouvoir l'afficher à l'issue du calcul.

trajet_init = trajet.copy()

Sélection des paramètres de l'algorithme

Je fixe le nombre de villes à 20, pour que le temps de calcul ne soit pas trop long en phase de mise au point.

N = 20 # nombre de villes

Les paramètres de température sont déclarés :

T0 = 10.0

Tmin = 1e-2

tau = 1e4

Le choix des valeurs de température initiale et de température minimum, celle qui fixe la fin du recuit, nécessite plusieurs essais. En règle générale, la valeur de T0 doit être du même ordre de grandeur que la variation d'énergie.

La valeur de \( \tau \) doit être choisie assez grande, pour que la diminution de la température soit suffisamment lente, mais pas trop pour ne pas avoir une masse trop grande de calculs.

Bien sur, je vous invite à modifier ces paramètres pour vérifier leur influence sur les résultats.

La mesure de l'énergie du système

L'énergie totale du système est mesurée par la somme des distances entre les villes dans l'ordre de parcours du trajet. Dans le principe, il suffit donc de calculer la distance entre les villes consécutives sur le trajet et de sommer ces distances. Dans la pratique, j'essaye d'optimiser au mieux ce calcul. Quelques petites précisions pour éclairer ce code, qui peut sembler obscur:

l'opérateur c_ concatène deux matrices, ici les matrices 1D de chaque coordonnée des villes, afin de construire une matrice des coordonnées des villes dans l'ordre du trajet.
la fonction roll, avec le paramètre -1, opére un shift gauche d'un élément (-1), ce qui me code \( (x_{i+1} - x_i)^2 \) de manière assez rapide.

Ce qui nous donne la fonction EnergieTotale() :

def EnergieTotale():

global trajet

energie = 0.0

coord = c_[x[trajet],y[trajet]]

energie = sum(sqrt(sum((coord - roll(coord,-1,axis=0))**2,axis=1)))

return energie

La création du fluctation dans le système

La fonction Fluctuation "swappe" un certain nombre de segments, fixés entre i et j, dans le trajet, ce qui simule une fluctuation d'énergie dans le système. Tel qu'elle est codée, cette fonction provoque le swap inverse si elle est appelée une seconde fois avec les mêmes valeurs (i,j).

def Fluctuation(i,j):

global trajet

Min = min(i,j)

Max = max(i,j)

trajet[Min:Max] = trajet[Min:Max].copy()[::-1]

return

L'implémentation de l'algorithme de Metropolis

C'est la pure transcription de l'algorithme décrit ci-dessus en pseudo-code :

def Metropolis(E1,E2):

global T

if E1 ? E2:

E2 = E1 # énergie du nouvel état = énergie système

else:

dE = E1-E2

if random.uniform() > exp(-dE/T): # la fluctuation est retenue avec

Fluctuation(i,j) # la proba p. sinon retour trajet antérieur

else:

E2 = E1 # la fluctuation est retenue

return E2

La loi de refroidissement

Là aussi, rien de particulier à dire. Vous pouvez toutefois changer la loi de variation, à condition que cette variation soit lente, sinon l'algorithme de recuit simulé ne fonctionnera pas.

T = T0*exp(-t/tau)

la boucle principale de calcul

Dans la boucle principale, on reprend l'algorithme de recuit simulé. Je provoque une fluctuation aléatoire d'énergie dans le système, en d'autres termes, je provoque une petite modification du trajet entre deux villes différentes :

i = random.random_integers(0,N-1)

j = random.random_integers(0,N-1)

if i == j: continue

Fluctuation(i,j)

Puis je calcule l'énergie du système après cette fluctuation :

Ef = EnergieTotale()

et j'applique l'algorithme de Metropolis :

Ec = Metropolis(Ef,Ec)

J'applique ensuite la loi de refroidissement et je boucle.

Les résultats

J'ai lancé le script avec N, le nombre de villes, égal à 20. C'est très peu et le temps de calcul est court (quelques secondes sur mon Mac). Le temps de calcul restera assez court pour N inférieur à la centaine. Après, cela devient un peu pénible, sur un Mac i5. Je n'ai pas essayé sur un i7.

Voici la première figure affichée : à gauche le trajet au lancement du script. Le trajet entre les différentes villes est aléatoire. A droite, s'affiche le trajet optimisé entre les différentes villes. Vous pouvez toujours vous amuser à vérifier à la main...

Algorithme voyageur de commerce - trajet

La figure suivante montre l'évolution de l'énergie du système (à gauche) et de sa température (à droite) en fonction des itérations de l'algorithme de recuit simulé. Les axes des ordonnées sont gradués en semi-log, ce qui explique la droite d'évolution de la température.

Algorithme voyageur de commerce - Energie

Sur la vue d'évolution de l'énergie en fonction du temps (des itérations), vous remarquez les oscillations aléatoires de la variation d'énergie pendant la phase de recuit et la convergence vers des variations de plus en plus faibles puis nulles lorsque l'état d'équilibre thermique est atteint.

Vous pouvez refaire les calculs avec des paramètres de température initiale et de coefficient de variation différents, et même avec des lois de variation de température différentes.

Les scripts Python

Le script RecuitSimule.py étudié dans cette page est disponible dans le package RecuitSimule.zip