Recherche d'un minimum à une dimension

1. Introduction

Soit une fonction f d'une variable réelle et à valeurs réelles, admettant un minimum dans 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()
            

figA.pdf

Il s'agit de déterminer le minimum de la fonction, appelé aussi minimum absolu, situé à $x\approx -3,3$. Cette fonction possède aussi un minimum relatif à $x\approx 2$.

La recherche d'un minimum consiste à partir d'un point quelconque et à suivre la pente descendante jusqu'à parvenir au minimum. En l'absence d'informations plus précises, le point initial de la recherche est choisi aléatoirement (avec une densité uniforme) sur l'intervalle ]a,b[. Cependant, la recherche d'un minimum en suivant la pente descendante a autant de chances de conduire au minimum relatif qu'au minimum absolu. Il faut donc un algorithme pour isoler un intervalle dans lequel se situe le minimum absolu, avant de faire une recherche plus précise de ce minimum. 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 :

$$p=p_{0}exp(-\frac{E}{k_{b}T})\text{(1)}$$

où p₀ est une constante et k_b 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 équiprobables. 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 le nommer 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 :

$$p\left(x\right)=p_{0}exp(-\frac{f\left(x\right)}{T})\text{(2)}$$

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 :

• Choix d'un nouvel état $x_1=x+\Delta x$ ou $x_1=x-\Delta x$. La direction est aléatoire mais le déplacement $\Delta x$ est fixé au départ, assez petit par rapport à la largeur de l'intervalle.
• Si le nouvel état est plus probable que le précédent, c'est-à-dire si f(x₁)<f(x), alors ce nouvel état est accepté.
• Si le nouvel état est moins probable que le précédent, il est accepté avec une probabilité égale au rapport des probabilités des deux états, c'est-à-dire $exp(-(f\left(x_1\right)-f\left(x\right))/T)$. Pour ce faire, on tire élatoirement un nombre réel avec une densité de probabilité uniforme sur l'intervalle [0,1] et on le compare à cette probabilité.

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. Recherche du minimum par encadrement itératif

L'algorithme du recuit simulé nous donne une valeur x₁ proche du minimum absolu, du moins avec une forte probabilité. Pour rechercher avec précision le minimum, on utilise un algorithme déterministe.

Il faut tout d'abord encadrer le minimum en cherchant x₁, x₂ et x₃ tels que f(x₂)<f(x₁) et f(x₂)<f(x₃), ce qui garantit la présence du minimum dans l'intervalle [x₁,x₃].

Soit $\delta x$ un déplacement assez petit (initialement positif). On considère tout d'abord x₂=x₁+δx puis on échange éventuellement x₁ et x₂ pour que f(x₂)<f(x₁) (auquel cas on change le signe de $\delta x$). En se déplaçant plusieurs fois de δx, on recherche ensuite une valeur x₃ telle que f(x₃)>f(x₂). L'algorithme est illustré sur la figure suivante pour trois cas différents.

Figure pleine page

L'algorithme doit fonctionner quelle que soit la position de x₁ par rapport au minimum et quelle que soit la valeur de $\delta x$, à condition que celle-ci soit assez petite (en valeur absolue) pour qu'un déplacement de $\delta x$ ne fasse pas sortir de la cuvette où se situe le minimum recherché.

On dispose à présent d'un encadrement du minimum par x₁ et x₃. On utilise une méthode itérative se répétant tant que |x₃-x₁|>ε, 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 [x₁,x₂] ou bien au milieu de l'intervalle [x₂,x₃]. Cette valeur vient remplacer une des trois valeurs (x₁,x₂,x₃) de telle sorte que l'on ait toujours f(x₂)<f(x₁) et f(x₂)<f(x₃) afin de maintenir l'encadrement du minimum. Cet algorithme est illustré sur la figure suivante :

Figure pleine page

4. Implémentation

Pour le tirage des nombres aléatoires avec une densité de probabilité uniforme sur l'intervalle [0,1[, on utilise la fonction random.random().

[1] Écrire une fonction recuit_simule(f,a,b) qui effectue le recuit simulé pour la fonction f sur l'intervalle [a,b]. On pourra faire varier la température 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 x₁, la dernière valeur obtenue. Pour chaque température, faire un affichage de cette température et de la valeur de x finale.

[2] 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 minimum absolu.

[3] Écrire une fonction encadrement(f,x1,delta_x) qui effectue le premier encadrement du minimum et renvoie (x₁,x₂,x₃).

[4] É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.

[5] 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')
              

figB.pdf

(x1,x2,x3)=encadrement(f,x1,0.01)
xmin=recherche_minimum(f,x1,x2,x3,1e-4)
              

print(xmin)
--> -3.281834395306403

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)
--> 0.998