Optimisation de recuit simulé pour le problème du voyageur

Bonjour tout le monde,

J'essaie de résoudre le problème du voyageur, autrement dit soit un ensemble de villes, passer par chacune d'entre elles une et une seule fois, le tout en minimisant la distance totale parcourue. Pour l'instant, je fais ça avec un algorithme de recuit simulé basique (sans redémarrages etc).

Je stocke les distances dans une matrice de taille N*N avec N le nombre de villes, en considérant que la distance entre A et B est la même qu'entre B et A (j'aurais pu utiliser juste une matrice triangulaire mais bon, c'est pas le plus important là). Pour savoir la distance entre la ville $i$ et la ville $j$, je regarde l'élément d'abscisse $i$ et d'ordonnée $j$ de la matrice. D'ailleurs, cette matrice est stockée au format CSV. La voilà:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

0,29,82,46,68,52,72,42,51,55,29,74,23,72,46  
29,0,55,46,42,43,43,23,23,31,41,51,11,52,21  
82,55,0,68,46,55,23,43,41,29,79,21,64,31,51  
46,46,68,0,82,15,72,31,62,42,21,51,51,43,64  
68,42,46,82,0,74,23,52,21,46,82,58,46,65,23  
52,43,55,15,74,0,61,23,55,31,33,37,51,29,59  
72,43,23,72,23,61,0,42,23,31,77,37,51,46,33  
42,23,43,31,52,23,42,0,33,15,37,33,33,31,37  
51,23,41,62,21,55,23,33,0,29,62,46,29,51,11  
55,31,29,42,46,31,31,15,29,0,51,21,41,23,37  
29,41,79,21,82,33,77,37,62,51,0,65,42,59,61  
74,51,21,51,58,37,37,33,46,21,65,0,61,11,55  
23,11,64,51,46,51,51,33,29,41,42,61,0,62,23  
72,52,31,43,65,29,46,31,51,23,59,11,62,0,59  
46,21,51,64,23,59,33,37,11,37,61,55,23,59,0  

Analyse des solutions

J'ai fait quelques tests en laissant tourner l'algo en boucle pendant la nuit jusqu'à trouver de bonnes solutions. La meilleure que j'ai obtenue était [9, 1, 12, 0, 10, 3, 5, 7, 13, 11, 2, 6, 4, 8, 14](ce sont les indices des villes, comptés à partir de 0) ce qui donne 294 km. Ca a mis vraiment longtemps avant de trouver cette solution, du coup je me suis demandé quelle était l'efficacité moyenne de cet algo. Pour cela, je fais 100 calculs de solution et je moyenne les distances obtenues. Il en est ressorti que plus la température initiale est élevée (j'avais mis $10^{50}$) ou plus le taux de refroidissement est faible, et meilleures sont les solutions, ce qui n'a rien de surprenant.

Code

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

#!/usr/bin/python2 -u
# The -u flag is very important, as it prevents output buffering

import random, math, csv, sys, copy, time

class TSP:
   def __init__(self, distancematrix):
       """ distancematrix is a matrix of size N*N, giving the distance between two cities
           A matrix allows us to have different distances from A to B than from B to A, but that wouldn't |make sense.
           Hence we'll use an upper triangular matrix, the rest of the values are undefined and we don't care |about them """
       self.distances = distancematrix
       self.size = len(self.distances[0])
   
   def stillaccept(self, rdist,sdist,t):
       if random.random() < math.exp((rdist - sdist)/t):
           return True
       else:
           return False

   def simulated_annealing(self, init_temp, coolrate):
       s = Solution(self.size) # Initial derpy random solution but that guarantees that all cities are |visited
       best = s
       bestdistance = best.totaldistance(self.distances)
       sdist = bestdistance
       t = init_temp
       while t > 1:
           r = copy.deepcopy(s).tweak()
           rdist = r.totaldistance(self.distances)
           if rdist < sdist or self.stillaccept(rdist,sdist,t):
               s = r
               sdist = rdist
           t *= 1-coolrate # Cool system
           if rdist < bestdistance:
               best = r
               bestdistance = rdist
       return best, bestdistance

class Solution:
   def __init__(self,size):
       self.sol = [i for i in range(0,size)]
       random.shuffle(self.sol)
   
   def tweak(self):
       """ Swap two random cities """
       i,j = random.randrange(0,len(self.sol)), random.randrange(0,len(self.sol))
       temp = self.sol[i]
       self.sol[i] = self.sol[j]
       self.sol[j] = temp
       return self
   
   def totaldistance(self,distancematrix):
       distance = 0
       for i in range(0, len(self.sol)-1):
           distance += distancematrix[self.sol[i]][self.sol[i+1]]
       return distance

   def __repr__(self):
       return self.sol.__repr__()

strmat = list(csv.reader(open(sys.argv[1], "r"), delimiter=','))
matrix = [ map(float, strmat[i]) for i in range (0, len(strmat)) ]
tsp = TSP(matrix)
mean = 0.0
for i in range(0,100):
   s, d = tsp.simulated_annealing(10**float(sys.argv[2]), float(sys.argv[3]))
   mean += d / 100.0
print("Average: {}".format(mean))

#bestsol, bestdist = tsp.simulated_annealing(10**float(sys.argv[2]), float(sys.argv[3]))
#while True:
#    newsol, newdist = tsp.simulated_annealing(float(sys.argv[2]), float(sys.argv[3]))
#    if newdist < bestdist:
#        bestsol = newsol
#        bestdist = newdist
#        print("[{}] Best distance so far: {} for {}".format(time.strftime("%H:%M"),bestdist, bestsol))

Pour faire tourner ce code, il faut faire ./tsp.py <fichier csv> <exposant de la température> <taux de refroidissement> (j'ai mis respectivement 50 et 0.1 pour les tests je crois). Si vous voulez que j'explique plus comment j'ai procédé, dites-le vu que je commente très peu.

Maintenant, je voudrais améliorer encore les performances, autrement qu'en modifiant température et taux de refroidissement, car la valeur moyenne de 433 km est franchement éloignée des 300 km presque optimaux. La page Wikipédia suggère de faire des redémarrages, mais qu'est-ce que je peux faire d'autre pour améliorer l'algo ou le code ?

Merci d'avance

Je vais me pencher sur la fonction tweak() du coup. Ton papier m'intéresse, ça m'aiderait sans doute à visualiser les possibilités de mutations.

En fait le 1 c'est un nombre parfaitement pifométrique pour décider quand arrêter. Comme j'ai remarqué que les solutions sont rarement très bonnes, j'ai une étape qui correspondrait peut-être à ce que tu appelles recherche locale:

1
2
3
4
5
6
7

bestsol, bestdist = tsp.simulated_annealing(10**float(sys.argv[2]), float(sys.argv[3]))
while True:
    newsol, newdist = tsp.simulated_annealing(float(sys.argv[2]), float(sys.argv[3]))
    if newdist < bestdist:
        bestsol = newsol
        bestdist = newdist
        print("[{}] Best distance so far: {} for {}".format(time.strftime("%H:%M"),bestdist, bestsol))

Ca tourne en boucle infinie pour l'instant, mais je pourrais faire en sorte de l'arrêter qund je suis suffisamment proche des 300km. Est-ce que c'était ça ce que tu voulais dire ?

Pour faire tourner ce code, il faut faire ./tsp.py <fichier csv> <exposant de la température> <taux de refroidissement> (j'ai mis respectivement 0.50 et 0.1 pour les tests je crois). Si vous voulez que j'explique plus comment j'ai procédé, dites-le vu que je commente très peu.

J'imagine que tu veux dire 50.0 et 0.1, non ?

A ta place, je prendrai 0.01. Avec ce réglage, j'obtiens ces stats (sur 100 essais) :

• min : 341
• max : 429
• avg : 394.5

L'essentiel de l'optimisation de ce genre d'algo repose justement sur le choix des paramètres, tu va avoir du mal à faire autrement. Sinon, ton code est certainement améliorable, mais le seul impact sera sur la vitesse d’exécution (et encore je ne sais pas combien on peut gagner), pas sur la qualité des solutions.

Pour les autres astuces possibles, je n'en connais pas en particulier, mais je ne suis absolument pas expert. Tu pourrais par exemple essayer une décroissance linéaire de la température.

Pour la corrélation avec les algos génétiques, attention quand même: en SA il faut un solution voisine, ce qui n'est pas forcément toujours le cas avec les mutation génétiques (où le but est de combiner des génomes, et pas de trouver une solution proche).

Ouais je me suis planté, c'était 50 et 0.1. J'ai fait un nouveau benchmark avec 1000 essais et pour une température de $10^{50}$:

C'est marrant, je trouve un maximum et une moyenne bien plus élevés que toi pour 0.01. Par contre 0.001 c'est bien trop lent même si ça donne de très bon résultats.

C'est très curieux. Je viens de tester avec 0.001 sur 100 essais, et j'obtiens :

• min : 328
• max : 396
• avg : 367.5
• std : 14

L'écart me semble beaucoup trop grand pour être explicable simplement. La seule différence notable est que j'utilise python3, mais je ne vois pas pourquoi ça aurait une influence (peut-être un autre type de random ?).

Si tu veux, voici ton code source que j'ai modifié à ma sauce (mais sans toucher à l'algo, en principe).

EDIT: Ah, en fait il y a une petite différence au niveau du tweak: dans ton cas, il y a une légère probabilité qu'il renvoie la même solution.

J'ai mis random.sample dans ma fonction tweak mais j'ai toujours cette différence. c'est vraiment bizarre, il faut que je regarde ça de plus près. Sinon, je dois te dire que je préfère de loin ton code

Je vais me pencher sur la fonction tweak() du coup. Ton papier m'intéresse, ça m'aiderait sans doute à visualiser les possibilités de mutations.

Voici l'article d'origine.

Ce n'est pas beaucoup plus détaillé que ma citation pour les parties qui t'intéressent. Je n'arrive plus à retrouver les sources qui détaillent la dernière méthode.

Ca tourne en boucle infinie pour l'instant, mais je pourrais faire en sorte de l'arrêter quand je suis suffisamment proche des 300km. Est-ce que c'était ça ce que tu voulais dire ?

J'ai du mal à comprendre ce que fait exactement ton code, désolé. Pour la condition d’arrêt je pensais juste à expliciter la liste des voisins à chaque itérations (il y a $n(n-1)$ voisins différents possibles), et à sélectionner le meilleur. C'est un algorithme glouton, qui ressemble à un hill climbing. Ton algorithme s'arrête uniquement lorsqu'il n'existe aucun voisin immédiat avec un meilleur trajet. Ainsi tu es sûr d'arriver assez vite dans un optimum local (pas global).

Dans un premier temps, je te conseille de travailler sur ton cooling schedule. Ici tu as un cooling arithmétique qui est rarement, en pratique, le plus efficace, par rapport à un cooling géométrique de type $T_n = \beta T_{n-1}$ pour un $\beta \in ]0,1[$.
Tu peux également tester un cooling logarithmique de type $T_n = \frac{T_0}{(n-1)^2}$ assez utilisé aussi.

En règle général, le gros problème est toujours le parameter tuning sur ce genre de méthodes approchées. Du coup, pour vraiment calibrer et tirer le meilleur de ton algorithme sur ce problème particulier (tu ne réussiras pas à avoir un bon algorithme sur un gros ensemble de problèmes ou d'instances : No Free Lunch Theorem), tu n'as pas le choix :

• Plan d'expérience
• Analyse de variance (ANOVA)
• Optimisation offline automatique via un framework comme ParamILS

Quelque chose de plus intérêt pour toi et pour l'exercice serait de te tourner vers de l'optimisation online de paramètres, qui lorsque les stratégies sont bien conçues, outperforms toujours des stratégies hors ligne ou statique. Par ailleurs, ces stratégies sont bien plus dans l'esprit dynamique de ces méthodes stochastiques.

Pourtant quand je fais le refroidissement, je fais t *= 1-coolingrate, pour moi c'est géométrique.

D'ailleurs, j'ai retapé mon code en m'inspirant fortement des modifications de yoch, et je trouve la même chose que lui maintenant. Je vais tenter d'ajouter la recherche locale à la fin vu que ça a l'air de bien améliorer les résultats.

Je n'avais pas vu l'etoile devant le signe egal. Ce qui de toute maniere n'aurait pas eu de sens. J'ai lu trop vite. Ceci dit, la remarque reste valide.

J'essaierai le refroidissement logarithmique quand j'aurai fini la recherche gloutonne. En tout cas, j'apprécie beaucoup votre aide !

EDIT:

J'ai implémenté la recherche gloutonne. C'est assez efficace et parfois ça donne de très très bons résultats:

1
2
3
4
5

λ ./travelingsalesman.py test.txt 50 0.001
Best: 319.0 for [9, 7, 11, 13, 2, 6, 14, 8, 4, 1, 12, 0, 10, 3, 5]
Worst: 395.0 for [11, 2, 6, 4, 1, 8, 7, 14, 12, 0, 10, 9, 13, 5, 3]
Average: 367.27
Greedy local search gives 268.0 for [7, 9, 13, 11, 2, 6, 4, 8, 14, 1, 12, 0, 10, 3, 5]

Voilà le code modifié:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119

#!/usr/bin/python2 -u
# The -u flag is very important, as it prevents output buffering

import random, math, csv, sys, itertools, time

class TSP:
   def __init__(self, distancematrix):
       """ distancematrix is a matrix of size N*N, giving the distance between two cities
           A matrix allows us to have different distances from A to B than from B to A, but that wouldn't make sense.
           Hence we'll use an upper triangular matrix, the rest of the values are undefined and we don't care about them """
       self.distances = distancematrix
       self.size = len(self.distances[0]) # Number of cities
   
   def stillaccept(self, rdist,sdist,t):
       return random.random() < math.exp((rdist - sdist)/t)
   
   def tweak(self, solution):
       """ Swap two random cities """
       solution = solution[::] # Copy solution
       i,j = random.sample(range(len(solution)), 2)
       solution[i], solution[j] = solution[j], solution[i]
       return solution

   def totaldistance(self,solution):
       return sum(self.distances[i][j] for i, j in zip(solution, solution[1:]))

   def seed_solution(self):
       s = list(range(len(self.distances[0])))
       random.shuffle(s)
       return s


class SimulatedAnnealing:
   def __init__(self, problem, coolingrate, inittemp):
      self.pb = problem
      self.cRate = coolingrate
      self.initT = inittemp

   def run(self):
       pb = self.pb
       s = pb.seed_solution()
       best = s
       bestdistance = pb.totaldistance(s)
       sdist = bestdistance
       t = self.initT
       while t > 1:
           r = pb.tweak(s)
           rdist = pb.totaldistance(r)
           if rdist < sdist or pb.stillaccept(rdist,sdist,t):
               s = r
               sdist = rdist
           t *= 1-self.cRate # Cool system
           if rdist < bestdistance:
               best = r
               bestdistance = rdist
       return best, bestdistance

class GreedyLocalSearch:
   def __init__(self, problem):
       self.pb = problem

   def swap(self, s, swap):
       i = swap[0]
       j = swap[1]
       s = s[::]
       s[i], s[j] = s[j], s[i]
       return s
   
   def bestNeighbor(self, s, sdist):
       # Yield better neighbor
       size = len(s)
       possibleSwaps = itertools.product(range(size), range(size))
       best = s
       bdist = sdist
       for x in possibleSwaps:
           new = self.swap(best, x)
           newdist = self.pb.totaldistance(new)
           if newdist < bdist:
               best = new; bdist = newdist
       return best, bdist

   def run(self, seed):
       pb = self.pb
       size = len(seed)
       best = seed
       bdist = pb.totaldistance(seed)
       new, newdist = self.bestNeighbor(best, bdist)
       while new != best and newdist != bdist: # Loop while there is a better neighbor
           best, bdist = new, newdist
           new, newdist = self.bestNeighbor(best, bdist) 
       return best, bdist

strmat = list(csv.reader(open(sys.argv[1], "r"), delimiter=','))
matrix = [ map(float, strmat[i]) for i in range (0, len(strmat)) ]

tsp = TSP(matrix)
sa = SimulatedAnnealing(tsp, float(sys.argv[3]), 10**float(sys.argv[2]))

bestsol, bestdist = sa.run()
worstsol, worstdist = bestsol, bestdist
mean = bestdist/100
for i in range(0,99):
   newsol, newdist = sa.run()
   mean += newdist/100
   if newdist < bestdist:
       bestsol = newsol
       bestdist = newdist
   if newdist > worstdist:
       worstsol = newsol
       worstdist = newdist
print("Best: {} for {}".format(bestdist, bestsol))
print("Worst: {} for {}".format(worstdist, worstsol))
print("Average: {}".format(mean))

# Now that we have a good solution, let's apply greedy local search to it

greedy = GreedyLocalSearch(tsp)
best, bdist = greedy.run(bestsol)
print("Greedy local search gives {} for {}".format(bdist, best))

Juste une petite mise à jour pour vous montrer un graphique que j'ai fait sur l'amélioration par algorithme glouton. SA c'est le recuit simulé et SA+GA c'est recuit simulé avec recherche gloutonne après. C'est un violin plot, ça permet de voir la densité de probabilité en plus du reste des informations d'une simple boîte à moustaches. J'obtiens ceci:

C'est assez parlant je trouve. Il faudrait que je fasse un autre graphique avec les delta.