Table des matières PDF

Recherche d'un minimum à une dimension

1. Introduction

Soit une fonction f d'une variable réelle et à valeurs réelles, admettant un minimum absolu sur l'intervalle ]a,b[, par exemple :

from matplotlib.pyplot import *
import numpy
def f(x):
    return x**4+2*x**3-12*x**2-2*x+77.37
a=-5
b=4
x=numpy.linspace(a,b,1000)
y=f(x)
figure()
plot(x,y)
xlabel("x")
ylabel("f(x)")
grid()
            
figAfigA.pdf

Il s'agit de déterminer le minimum absolu de la fonction. Si l'on part d'un point quelconque de l'intervalle ]a,b[, la recherche d'un minimum vers la pente descendante a autant de chances de conduire au minimum relatif qu'au minimum absolu. Il faut donc un algorithme pour repérer le minimum absolu avant de faire une recherche plus précise. Nous utiliserons pour cela l'algorithme du recuit simulé.

2. Algorithme du recuit simulé

L'algorithme du recuit simulé est inspiré des méthodes de Monte-Carlo appliquées à la statistique de Boltzmann, et en particulier l'algorithme de Metropolis. Supposons que la fonction f définie ci-dessus représente l'énergie d'un système physique en fonction d'un paramètre x, définissant par exemple la structure cristalline d'un solide. D'après la loi de Boltzmann, la probabilité d'un état d'énergie E pour un système à la température T est :

p0 est une constante et kb la constante de Boltzmann.

Lorsque la température est suffisamment basse, le système se place au voisinage d'un état de minimum de l'énergie, soit le minimum absolu (état stable), plus rarement le minimum relatif (état métastable). Lorsque la température est très élevée, toutes les valeurs de x de l'intervalle sont probables. Si le corps est refroidi lentement, il converge vers son état stable. Un refroidissement rapide peut le conduire dans un état métastable, c'est-à-dire un minimum relatif de l'énergie. Ce phénomène est mis à profit en métallurgie pour obtenir des structures cristallines métastables aux propriétés mécaniques intéressantes. En recuisant un cristal se trouvant dans un état métastable puis en le laissant refroidir lentement on peut l'amener dans son état stable.

L'algorithme du recuit simulé (simulated annealing) consiste en fait à simuler le processus de refroidissement très lent, en partant d'une température assez élevée. On devrait donc plutôt l'appeler algorithme du refroidissement simulé. L'algorithme permet de placer la variable x au voisinage du minimum absolu, mais pas de manière certaine. Avec un bon réglage de la vitesse de décroissance de la température, elle permet d'obtenir le minimum absolu avec une forte probabilité.

Pour notre problème, la loi de Boltzmann s'écrit :

Pour une température T donnée, l'algorithme de Metropolis consiste à générer une suite d'états (ici de valeurs de x) qui converge vers un ensemble d'états obéissant à la loi de probabilité de Boltzmann (ou une autre loi). On commence par une valeur x choisie aléatoirement (avec une densité de probabilité uniforme) dans l'intervalle ]a,b[. La suite d'états est générée par l'algorithme suivant :

Dans une simulation de Metropolis d'un système physique, on fait des calculs statistiques avec la suite des états ainsi générée. Dans l'algorithme du recuit simulé, seule la dernière valeur de x est utilisée. L'itération est répétée un grand nombre de fois avec la même température, puis on recommence avec une température un peu plus faible, en partant de la dernière valeur de x obtenue avec la température précédente. Avec une température finale assez faible, l'état final obtenu est très probablement très proche du minimum absolu.

Cet algorithme nécessite des réglages qui dépendent du problème à résoudre : valeur de Δx, nombre d'itérations à chaque température, températures initiale et finale, loi de décroissance de la température.

3. Calcul du minimum par encadrement itératif

L'algorithme du recuit simulé nous donne une valeur x1 proche du minimum absolu, du moins avec une forte probabilité. Pour rechercher avec précision le minimum, on peut à présent utiliser un algorithme déterministe.

Il faut tout d'abord encadrer le minimum en cherchant x1 et x2 tels que f(x2)<f(x1) et f(x2)<f(x3), ce qui garantit la présence du minimum dans l'intervalle [x1,x3]. On considère tout d'abord x2=x1+δx puis on échange éventuellement x1 et x2 pour que f(x2)<f(x1). En se déplaçant plusieurs fois de δx dans la même direction, on recherche ensuite une valeur x3 telle que f(x3)>f(x2).

On dispose à présent d'un encadrement du minimum. On utilise une méthode itérative se répétant tant que |x3-x1|<ε, où ε est une tolérance égale au moins à la racine carrée de la précision des flottants. À chaque itération, une nouvelle valeur x est prise alternativement au milieu de l'intervalle [x1,x2] ou bien au milieu de l'intervalle [x2,x3]. Cette valeur vient remplacer une des trois valeurs (x1,x2,x3) de telle sorte que la valeur x2 vérifie toujours f(x2)<f(x1) et f(x2)<f(x3) afin de maintenir l'encadrement du minimum.

4. Implémentation

On utilise le module random et la fonction random.random(). Pour ne pas générer toujours la même suite de nombres aléatoires, on initialisera le générateur avec random.seed().

Écrire une fonction recuit_simule(f,a,b) qui effectue le recuit simulé pour la fonction f sur l'intervalle [a,b]. La rempérature pourra être variée de T=100 à T=1, avec une cinquantaine d'étapes comportant chacune un millier d'itérations. On pourra adopter le déplacement Δx=0.1. La fonction renvoie x1, la dernière valeur obtenue. Pour chaque température, faire un affichage de cette température et de la valeur de x finale.

Tester plusieurs fois la fonction pour vérifier qu'elle donne le plus souvent une valeur proche du minimum absolu. Observer comment la valeur de x s'approche progressivement du minimu absolu.

Écrire une fonction encadrement(f,x1,delta_x) qui effectue le premier encadrement du minimum et renvoie (x1,x2,x3).

Écrire une fonction recherche_minimum(f,x1,x2,x3,tol) qui effectue la recherche itérative par encadrement et renvoie le minimum. Tester avec une tolérance 10-4.

Mettre en place un test visant à déterminer la probabilité de trouver le minimum absolu et non pas le minimum relatif. Tester différentes vitesses de refroidissement.

import random
def recuit_simule(f,a,b):
    random.seed()
    x=random.random()*(b-a)+a
    y=f(x)
    delta_x=0.1
    Temp=numpy.logspace(2,0,100)
    liste_x=[]
    liste_T=[]
    for T in Temp:
        for i in range(1000):
            u=random.random()
            if u<0.5:
                x1=x+delta_x 
            else:
                x1=x-delta_x
            y1=f(x1)
            if y1<y:
                x=x1
                y=y1
            else:
                u=random.random()
                if u<numpy.exp(-(y1-y)/T):
                    x=x1
                    y=y1 
        liste_x.append(x)
        liste_T.append(T)
    return (x,liste_x,liste_T)
    
def encadrement(f,x1,delta_x):
    y1=f(x1)
    direction=1
    x2=x1+direction*delta_x
    y2=f(x2)
    if y2>y1:
        (x,y)=(x1,y1)
        (x1,y1)=(x2,y2)
        (x2,y2)=(x,y)
        direction=-1
    x3=x2+direction*delta_x
    y3=f(x3)
    while y3<y2:
        x3 += direction*delta_x
        y3 = f(x3)
    return (x1,x2,x3)

def recherche_minimum(f,a,b,c,tol):
    fa=f(a)
    fb=f(b)
    fc=f(c)
    while abs(a-c)>tol:
        x=(a+b)/2
        fx=f(x)
        if fx>fb:
            (a,b,c)=(x,b,c)
        else:
            (a,b,c)=(a,x,b)
        fa=f(a)
        fb=f(b)
        fc=f(c)
        x=(b+c)/2
        fx=f(x)
        if fx>fb:
            (a,b,c)=(a,b,x)
        else:
            (a,b,c)=(b,x,c)
        fa=f(a)
        fb=f(b)
        fc=f(c)
    return (a+c)/2
              
a=-5
b=4
(x1,x,T)=recuit_simule(f,a,b)
figure()
plot(T,x,"ko")
grid()
xlabel("T")
ylabel("x")
xscale('log')
              
figBfigB.pdf
(x1,x2,x3)=encadrement(f,x1,0.01)
xmin=recherche_minimum(f,x1,x2,x3,1e-4)
              
print(xmin)
--> -3.281830596496343

Calcul du taux de réussite :

n=0
N=500
for k in range(N):
    (x1,x,T)=recuit_simule(f,a,b)
    if x1<-1:
        n += 1
p=n/N
              
print(p)
--> 1.0
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.