Nous soussignés Anthony LÉONARD et Luc LIBRALESSO déclarons sur l’honneur que ce rapport est le fruit d’un travail personnel, en binôme, que nous n’avons ni contrefait, ni falsifié, ni copié tout ou partie de l’œuvre d’autrui afin de la faire passer pour nôtre. Toutes les sources d’information utilisées et les citations d’auteur ont été mentionnées confor- mément aux usages en vigueur. Nous sommes conscient(e)s que le fait de ne pas citer une source ou de ne pas la citer claire- ment et complètement est constitutif de plagiat, que le plagiat est considéré comme une faute grave au sein de l’Université, pouvant être sévèrement sanctionnée par la loi.

Explications

Introduction

Le recuit simulé est une méthode d’optimisation originaire de la métallurgie qui consiste à chauffer un morceau de métal et le refroidir doucement pour que les atomes se placent de façon régulière. Le même principe est utilisé en optimisation. C’est cet algorithme que nous implémenterons dans ce TP.

Étude de l’algorithme

A chaque iteration, \(\Delta\) sera négatif si l’action a amélioré l’ordonnancement, et positif si l’action à déterioré l’ordonnancement.

Pour chaque nouvelle perturbation de notre solution, nous avons la condition d’acceptation qui est la suivante :

\(Random < exp(\frac{-\Delta}{T})\)

\(Random\) allant de 0 à 1. Si \(\frac{-\Delta}{T}\) est positif, la perturbation sera gardée à tous les coups. Cependant, si \(\frac{-\Delta}{T}\) est négatif, le changement pourra ou pas être accepté selon la valeur de \(Random\)

Influence des différents paramètres

Commençons par nous intéresser à \(\Delta\) :

Quand \(\Delta\) est négatif : \(-\frac{\Delta}{T}\) sera positif et quand \(\Delta\) est positif, \(-\frac{\Delta}{T}\) sera négatif. Donc \(\exp{-\frac{\Delta}{T}}\) sera compris entre 0 et 1 quand \(\Delta\) est positif et plus grand que 1 quand \(\Delta\) est négatif.

Par conséquent, du moment que la solution est améliorée, on garde le changement, mais on s’autorise aussi à accepter un changement qui déteriore la solution dans certains cas. A savoir que plus la perturbation déteriore la solution, moins il a de chances d’être accepté.

Interessons nous maintenant à \(T\) :

Quand \(T\) est grand, \(\frac{-\Delta}{T}\) se rapproche de \(0\), peu importe le signe de \(\Delta\). Donc \(exp(\frac{-\Delta}{T})\) tendra vers 1.

Par conséquent, quand \(T\) est grand, la probabilité d’accepter un changement, même si il déteriore la solution est proche de 1.

Cependant, quand \(T\) est petit, \(\Delta\) agit d’avantage sur la probabilité d’accepter un changement. Et donc, quand \(T\) est proche de 0, la probabilité d’accepter un changement alors qu’il déteriore la solution est quasiment nulle.

Lien avec la métallurgie :

Quand la température est élevée, les atomes bougent rapidement et de manière chaotique. Quand la température est plus faible, les atomes suivent d’avantage de règles et s’organisent régulièrement par rapport à leurs voisins.

Intuition sur pourquoi cet algorithme marche bien :

Quand on a une température basse, on se rapproche d’un minimum local. Quand on a une température élevée, on peut sortir d’un minimum local plus facilement.

Intuitivement, le fait de faire diminuer la température doucement nous permet de sortir des minimums locaux qui ne sont pas très bas et donc d’avoir une solution meilleure qu’en faisant descendre rapidement la température (on cherche un minimum local sans se donner la possibilité d’en sortir).

Si \(T\) décroit trop vite : La solution générée sera moins bonne, mais l’algorithme plus rapide. Si \(T\) décroit trop lentement : Nous resterons dans un état où nous ne chercherons pas à améliorer la solution trop longtemps ce qui n’est pas intéressant.

On note qu’un juste milieu est à trouver dans la façon dont \(T\) décroit.

Algorithme de génération

Dans cette partie, nous nous intéressons à la génération d’un ensemble de \(n\) tâches de même loi.

On commence par initialiser la seed :

set.seed(42)

Nous pouvons maintenant implémenter différents algorithmes de génération de \(n\) tâches :

gen1 <- function(n){
  return(runif(n))
}

\(gen1\) renvoie un ensmble de \(n\) tâches de durée uniforme sur \([0,1[\).

gen2 <- function(n){
  res <- rnorm(n)
  return(abs(res))
}

\(gen2\) renvoie un ensemble de \(n\) tâches de durée suivant une loi normale pour lesquelles on a pris la valeur absolue.

Comparons maintenant la proportion de tâches d’une certaine taille générés nos algorithmes pour un échantillon de 1000 tâches.

hist(gen1(1000), breaks=30, xlab="durée de la tâche", ylab="nombre d'occurences", main="Fréquence d'apparition des tâches en fonction de leur durée")

hist(gen2(1000), breaks=30, xlab="durée de la tâche", ylab="nombre d'occurences", main="Fréquence d'apparition des tâches en fonction de leur durée")

Étude expérimentale

Nous pouvons maintenant nous intéresser aux politiques de décroissance de température pour une instance du problème avec \(n=20\) et \(K=5\) où \(n\) est le nombre de tâches et \(K\) le nombre de machines. Nous fixons aussi une température max \(Tmax\) à 100.

n <- 20
K <- 5
Tmax <- 100

Schémas de décroissance

Chacune des fonctions de décroissance a deux paramètres :

  • \(Tindex\) qui correspond à l’index de la température demandée

  • \(nbIndex\) correspond au nombre de températures demandées (fixé à 500 par défaut)

Pour les fonctions suivantes, on définit à l’interieur une fonction \(f\) qui représente l’allure de la décroissance (linéaire, sinusoïdale, quadratique, …) qui prend un paramètre \(i\) compris entre 0 et 1 et renvoie un nombre entre 0 et 1. Le reste de la fonction applique une mise à l’échelle de cette valeur. Cela permet de modifier plus facilement l’allure des courbes au besoin.

dec1 <- function(Tindex, nbIndex=500){
  f <- function(i){
    return(i)
  }
  fi <- 1-f(Tindex/(nbIndex-1))
  return( Tmax*( fi ) )
}
dec2 <- function(Tindex, nbIndex=500){
  f <- function(i){
    return(sin(i*pi/2))
  }
  fi <- 1-f(Tindex/(nbIndex-1))
  return( Tmax*( fi ) )
}
dec3 <- function(Tindex, nbIndex=500){
  f <- function(i){
    return(i*i)
  }
  fi <- 1-f(Tindex/(nbIndex-1))
  return( Tmax*( fi ) )
}
dec4 <- function(Tindex, nbIndex=500){
  f <- function(i){
    return(i^3)
  }
  fi <- 1-f(Tindex/(nbIndex-1))
  return( Tmax*( fi ) )
}

Réalisons maintenant un tracé graphique de ces fonctions :

nbIndex = 500
x <- 1:nbIndex
y1 <- dec1(x)
y2 <- dec2(x)
y3 <- dec3(x)
y4 <- dec4(x)
plot(x,y1, type="l", col="blue")
lines(x,y2, col="red")
lines(x,y3, col="green")
lines(x,y4, col="black")

legend("topright", 
       legend=c(
         "dec1 (linéaire)",
         "dec2 (sinusoïdale)",
         "dec3 (quadratique)",
         "dec4 (cubique)"
       ),
       lwd=1, col=c("blue","red","green", "black")
)

Implémentation de l’algorithme du recuit simulé

On commence par générer un ensemble de tâches en utilisant un des générateurs.

taches <- gen1(n)

On définit ensuite un ensemble de fonctions annexes à commencer par l’allocation initiale.

On attribut une tâche à une machine puis on passe à la tâche et la machine suivante. On retourne ensuite une matrice où les lignes correspondent aux machines et les colonnes aux tâches. À chaque intersection, on trouve un 1 si la tâche correspondant à cette colonne est attribuée à la machine représentée par la ligne, sinon on trouve 0.

init_alloc <- function(){
  mat <- matrix(0, nrow = K, ncol = n)
  for(i in 1:n) {
    mat[(i%%(K))+1,i] = 1
  }
  mat
}

Viennent ensuite les fonctions permettant de tirer au hasard une machine et une tâche.

random_task <- function() {
  floor((runif(1)*n)+1);
}

random_machine <- function() {
  floor((runif(1)*K)+1);
}

Une fonction permettant de déplacer une tâche vers une machine donnée en retournant une nouvelle matrice.

move <- function(mat, i, j){
  nmat <- mat;
  for(k in 1:K){
    nmat[k,i] = 0
  }
  nmat[j,i] = 1
  nmat
}

Une fonction calculant le temps global d’exécution d’une configuration donnée. On calcule pour chaque machine la durée d’exécution de ses tâches et l’on renvoit la durée maximale.

makespan <- function(mat) {
  max = 0
  x = 0
  for(i in 1:K){
    for(j in 1:n){
      x = x + (taches[j]*mat[i,j])
    }
    if(x > max) max = x
    x = 0
  }
  max
}

Finalement, on simule le recuit et l’on renvoit la configuration d’ordonnancement trouvée. On lui donne en paramètre la fonction de décroissance mise au point précédemment.

recuit <- function(dec) {
  c <- init_alloc()
  t <- 0
  temp <- dec(t) # Ici on choisit un des schémas de décroissance
  while (temp > 0) { 
    i <- random_task()
    j <- random_machine()
    cp <- move(c, i, j)
    delta <- makespan(cp) - makespan(c)
    if(runif(1) < exp(-(delta/temp)))
        c <- cp
    t <- t+1
    temp <- dec(t)
  }
  c
}

Mesures de performance par la simulation

Prenons une instance de notre problème :

taches <- gen1(n)
taches
##  [1] 0.59086189 0.51631583 0.93469628 0.26051241 0.71070244 0.16814803
##  [7] 0.81987237 0.75010702 0.28999593 0.17981304 0.61729193 0.61230199
## [13] 0.86571886 0.63971645 0.09189931 0.86737715 0.84175433 0.49549081
## [19] 0.07937267 0.02165050

Calculons la durée avec une allocation initiale :

makespan(init_alloc())
## [1] 3.046013

Faisons de même avec le résultat du recuit simulé avec les différentes politiques de décroissance de température :

makespan(recuit(dec1))
## [1] 4.175875
makespan(recuit(dec2))
## [1] 2.614516
makespan(recuit(dec3))
## [1] 2.839879
makespan(recuit(dec4))
## [1] 3.680587

A la vue des résultats : on constate que le le makespan d’une fonction sinusoïdale et quadratique est plus efficace que l’allocation initiale.

Ces résultats ne sont que le fruit d’une seule simulation, il faudrait donc en faire d’avantage. Mais sont quand même surprenants pour deux raisons :

  • Nous nous attendions à ce que la fonction cubique soit plus efficace que les autres et il se trouve que c’est celle qui a donné le deuxième moins bon résultat derrière la fonction linéaire.

  • Nous nous attendions également à avoir un ordre dans la qualité des résultats (par exemple, la fonction sinusoïdale donne un moins bon résultat que la fonction linéaire qui donne un moins bon résultat que la fonction quadratique et ainsi de suite). Ici, ce n’est pas le cas et il se trouve que la pire fonction est la fonction linéaire qui est entre les deux meilleures.

Réalisons d’avantage de tests :

Pour cela, nous créons une nouvelle instance \(nbTests\) fois. puis pour chacune de ces instances, lançons l’algorithme du recuit simulé avec une fonction de décroissance précise, puis réalisons la moyenne du makespan.

nbTests <- 10
moyInit <- 0
moy1 <- 0
moy2 <- 0
moy3 <- 0
moy4 <- 0
for (i in 1:nbTests) {
  taches <- gen1(n)
  moyInit <- moyInit + makespan(init_alloc())
  moy1 <- moy1 + makespan(recuit(dec1))
  moy2 <- moy2 + makespan(recuit(dec2))
  moy3 <- moy3 + makespan(recuit(dec3))
  moy4 <- moy4 + makespan(recuit(dec4))
}
moyInit <- moyInit / nbTests
moy1 <- moy1 / nbTests
moy2 <- moy2 / nbTests
moy3 <- moy3 / nbTests
moy4 <- moy4 / nbTests

Voyons maintenant ces résultats :

moyInit
## [1] 2.675388
moy1
## [1] 3.137522
moy2
## [1] 2.48678
moy3
## [1] 3.339614
moy4
## [1] 3.413452

Ok, on constate que la fonction sinusoïdale reste un bon choix car est la seule qui donne un meilleur résultat que le placement initial.

On note cependant par rapport à la fois précédente qu’on observe bien un ordre : la sinusoïde donne un meilleur résultat que la linéaire qui donne un meilleur résultat que la quadratique et ainsi de suite.

Essayons maintenant de comprendre pourquoi ces résultats :

Le placement initial que nous effectuons est souvent plutôt bon, car les tâches ont des durées plutôt similaires. Or le fait de passer vite à une température faible, nous permet d’en dévier que très légèrement et de ne garder que les placements qui améliorent la solution. C’est pour cela que la sinusoïde marche bien.

Cependant, pour les autres, étant donné qu’ils peuvent accepter des placements qui déteriorent la solution, il y a des chances de tomber dans un autre minimum local qui n’est pas aussi bon que celui de départ.

Nous pouvons maintenant nous intéresser aux résultats du recuit simulé sur l’autre manière de générer des tâches qui crée des tâches plus différentes que le premier. Est-ce que la sinusoïde reste efficace ? Est-ce que le placement initial reste efficace ?

Résultats avec le deuxième algorithme de génération des tâches

Nous reprenons le code précédemment écrit sur le tests des fonctions de décroissance et l’effectuons avec des tâches générées depuis le deuxième algorithme de génération de tâches.

nbTests <- 10
moyInit <- 0
moy1 <- 0
moy2 <- 0
moy3 <- 0
moy4 <- 0
for (i in 1:nbTests) {
  taches <- gen2(n)
  moyInit <- moyInit + makespan(init_alloc())
  moy1 <- moy1 + makespan(recuit(dec1))
  moy2 <- moy2 + makespan(recuit(dec2))
  moy3 <- moy3 + makespan(recuit(dec3))
  moy4 <- moy4 + makespan(recuit(dec4))
}
moyInit <- moyInit / nbTests
moy1 <- moy1 / nbTests
moy2 <- moy2 / nbTests
moy3 <- moy3 / nbTests
moy4 <- moy4 / nbTests

Voyons maintenant ces résultats :

moyInit
## [1] 4.473421
moy1
## [1] 4.647606
moy2
## [1] 3.818959
moy3
## [1] 4.664191
moy4
## [1] 5.579796

Finalement, l’ordre ne change pas, la sinusoïde reste largement devant les autres.

Conclusion

Notre implémentation du recuit simulé a pour avantage d’être exécutable en un temps indépendant du nombre de tâches à ordonnancer (le paramètre d’arrêt étant lié à la température, elle même indépendante du reste). La méthode la plus efficace pour la décroissance de T semble être une décroissance rapide au début et qui tend à devenir de plus en plus lente au fur et à mesure.

Cependant cet algorithme présente quelques inconvénients. Le premier est qu’il ne garanti par d’atteindre un minimum global à la fin de son exécution. De plus de nombreux paramètres sont à configurer (température initiale, vitesse de décroissance, etc…) et influent sur l’efficacité de l’algorithme. Ces paramètres eux-même dépendent de la taille de l’échantillon (décroissance plus lente si beaucoup de tâches par exemple). Il est donc nécessaire de faire des tests et d’ajuster les paramètres au fur et à mesure.

Cependant, on note que notre algorithme, une fois configuré, donne un résultat souvent superieur à l’heuristique qui fait le placement initial. Ce qui montre que le recuit simulé peut être une bonne approche pour donner une solution à une instance de problème difficile tout en maîtrisant le temps d’exécution.

Pour aller plus loin :