Continuons la discussion de ce matin, avec l’exemple des capitales américaines. On peut trouver une telle liste liste en ligne sur wikipedia,

> capital=read.table(
+ "http://freakonometrics.free.fr/US-capital.csv",
+ header=TRUE,sep=";")
> capital=capital[-c(2,11),]

On va enlever Hawaï et l’Alaska, qui forcent à faire de trop grands détours. Dans TSP, sous R, il y a des listes de villes américaines, avec les coordonnées. En bricolant un peu, on va récupérer les coordonnées de la plupart des capitales,

> library(maps)
> library(maptools)
> library(sp)
> library(TSP)
> Lcap=paste(as.character(capital[,1]),
+ as.character(capital[,2]),sep=", ")
> M=as.matrix(USCA312)
> L=labels(M)[[1]]
> Lcap=Lcap[which(Lcap%in%L)]
> k=which(L%in%Lcap)
> M=M[k,]
> M=M[,k]
> listeUSA=TSP(M,L[k])

Et cette fois on est prêt,

> tour=solve_TSP(listeUSA, method = "nn")

on a lancé l’algorithme, et on peut faire un dessin,

> COORD=data.frame(USCA312_coords[k,])
> COORD=COORD[as.numeric(tour),] 
> map('state',fill=TRUE,col="yellow")
> lines(COORD$coords.x1,COORD$coords.x2,
+ lwd=3,col="blue")
> points(COORD$coords.x1,COORD$coords.x2,
+ pch = 19, cex = 1.2, col = "red")

pas mal, n’est pas ? Bon, mais comme d’habitude, rien ne garantie que l’on ait trouvé le trajet optimal… En fait, le soucis est que si on relance l’algorithme, on retombe sur un autre trajet,

voire encore un autre si on lancer une nouvelle fois l’algorithme,

C’est pénible, n’est ce pas ? En fait, si on lance 1,000 fois l’algorithme on obtient la distribution (des distances optimales) suivante

autrement dit, on a encore beaucoup de variabilité, avec une distance optimale qui peut varier d’environ 10% entre deux boucles d’algorithmes. Cela dit, cela reste faible par rapport à un trajet au hasard,

Et d’ailleurs, si on compare la distance obtenue en tirant les villes au hasard, on arrive à un taux de l’ordre de 27% (qui peut être comparé à ce que nous avions obtenu en tirant au hasard dans le carré unité)

>  HAZARD=OPTIMAL=RATIO=rep(NA,500)
>  for(s in 1:500){
+   HAZARD[s]=tsp.longueur(M,sample(1:nrow(M)))
+   tour=solve_TSP(listeUSA, method = "nn")
+   OPTIMAL[s]=tsp.longueur(M,as.numeric(tour))
+   RATIO[s]=HAZARD[s]/
+            OPTIMAL[s]
+ }

Bref, compte tenu de la (grande) dispersion des résultats optimaux obtenu, on pourrait dire que cet algorithme est assez mauvais. Même s’il est très rapide. Un autre avantage est que l’on peut spécifier la ville de départ, par exemple. Ici, si on souhaite partir de Floride, on peut obtenir le trajet suivant

> tour=solve_TSP(listeUSA, method = "nn",
+ control=list(start=41))

Passionnant n’est-ce pas ? Et encore, ça ne va pas m’aider à trouver le parcours optimal pour les enfants, car pour les enfants, il faut spécifier le point de départ (voire le point d’arrivée), il faut en plus rester sur la route (voire traverser aux passages pour piétons). On peut aussi vouloir faire une pause pipi au milieu de la promenade…. Bref, on a pas mal de contraintes ! Mais tout est expliqué dans le livre fabuleux de William Cook, In Pursuit of the Traveling Salesman, qui reprend toute l’histoire des algorithmes les plus fins pour résoudre ce joli problème de maths…

Bon, ben cette fois ça y est, c’est l’Halloween… et comme tous les ans, mon rôle devrait se limiter à creuser les citrouilles, décorer la maison… et attendre à la porte de la maison que les enfants des voisins viennent me dévaliser ma réserve de bonbons ! Je pourrais à la rigueur revendiquer un petit quelque chose dans le déguisement de mon fils (on a passé pas mal de temps à feuilleter les premiers tomes de Walking Deads, histoire de mieux saisir l’essence des personnages de morts-vivants, en plus de la Zombie Walk la fin de semaine passée, en centre-ville).

Mais cette année, pas de billet sur les zombies comme l’an passé, mais un vrai problème pratique: comment optimiser son trajet, quand on sait où sont les bonnes maisons, où on a de fortes chances d’avoir plein de bonbons (en particulier chez les parents des amis) ? C’est assez proche de l’idée de passer par un maximum de maisons, et ce à une heure raisonnable (histoire qu’il reste des bonbons à récupérer).

Pour tous ceux qui ont fait un peu de programmation dynamique, c’est un problème classique, connu sous le nom du “voyageur de commerce“. Ce problème est né en juillet 1954 avec “the shortest route for a traveling salesman” paru dans Newsweek ci-dessous (que l’on retrouve en détails dans le livre 50 Years of Integer Programming)

Enfin, disons qu’il a été popularisé à cette époque, car Euler évoquait déjà un problème similaire quand il travaillait sur le problème des sept ponts de Königsberg, dans “Solutio problematis ad geometriam situs pertinentis“,

(même si ce problème consiste à étudier l’existence d’un trajet sous des contraintes de passages par des nœuds, et pas vraiment de l’optimisation). On retrouve aussi un exemple en 1832, en Allemagne, dans un livre intitulé Der Handlungsreisende – wie er sein soll und was er zu thun hat, um Auftraege zu erhalten und eines gluecklichen Erfolgs in seinen Geschaeften gewiss zu sein – Von einem alten Commis-Voyageur où un trajet passant par 45 villes en Allemagne et en Suisse devait être optimisé

(le schéma ci-dessus montre celle décrite en 1832, pour un total de 1,285km, mais il est possible de faire une boucle de 1,248km).

En fait le problème a explosé en 1962 (à cette époque, il semble que les gens pouvaient encore penser que faire des maths pouvait être le fun), lors que Procter & Gamble ont offert un prix de 10,000$ à celui qui trouverait le chemin le plus court pour que Toody et Muldoom, les conducteurs de la voiture 54 dans une série populaire à l’époque, passent par 33 villes, aux États-Unis.

Cela dit, Carl Menger avait noté dès 1930 a noté que ce problème devrait etre difficile à résoudre. Facile à écrire, mais difficile à résoudre explicitement. Dans son livre, William Cooke rappelle qu’avec 33 villes, il y aurait

131,565,418,466,846,765,083,609,006,080,000,000

chemins possibles, dont il faudrait calculer la longueur. Et qu’avec l’ordinateur le plus puissant en 2009 (effectuant 1.5 million de milliards d’opérations à la seconde) il faudrait 28,000 milliards d’année pour mener à bien les calculs. Même si on essayait de tenir compte des progrès de l’informatique dans les mille prochaines années, ça risque de prendre du temps !..

En 1967, Jack Edmonds notait “I conjecture that there is no good algorithm for the traveling salesman problem” où “good” est une manière élégante pour dire un temps polynomial en , le nombre de villes. Ce qui s’appelle les problèmes de la classe . Il semblerait que le problème du voyageur de commerce soit un problème  et pas de la classe (cf le chapitre du livre de Sanjeev Arora et Boaz Barak, ou le chapitre 11 des notes de cours de Olivier Bournez, que j’ai pu découvrir grâce à @dmonniaux, ou formellement de la classe  et pas , i.e.  complet), voir aussi la discussion sur wikipedia),

Et depuis, l’institut Clay a proposé $1,000,000 à celui qui prouverait que les problèmes   ne sont pas de type  (ou qui prouveraient qu’ils le sont… laconjecture reste ouverte, semble-t-il). Bref, les enjeux financiers derrière ont considérablement augmenté (on pourra consulter Sipser (2007) pour la petite histoire)

Bon, et si l’histoire est intéressante, ce n’est pas une raison pour la lire en se croisant les bras… essayons de programmer un algorithme (simple) afin d’approcher une solution, car je ne pourrais me contenter de dire à mes enfants “c’est un problème  complet, papa a capitulé“.

Commençons par tirer 5 points au hasard dans le carré unité,

>  n = 5
>  (x = matrix(runif(2*n),nr=n))
[,1]       [,2]
[1,] 0.4142621 0.35600725
[2,] 0.9477231 0.66938653
[3,] 0.2518656 0.19900873
[4,] 0.5314231 0.22954670
[5,] 0.5156522 0.09066741

et commençons par calculer la distance entre les points (i.e. une matrice, ou ici une matrice triangulaire inférieure)

>  (d = dist(x))
1         2         3         4
2 0.6186980
3 0.2258786 0.8399243
4 0.1723919 0.6056111 0.2812205
5 0.2840514 0.7222195 0.2851688 0.1397719

On va mettre ce triangle sous forme de matrice (pleine), ça sera plus simple à manipuler par la suite

>  (d = as.matrix(d))
1         2         3         4         5
1 0.0000000 0.6186980 0.2258786 0.1723919 0.2840514
2 0.6186980 0.0000000 0.8399243 0.6056111 0.7222195
3 0.2258786 0.8399243 0.0000000 0.2812205 0.2851688
4 0.1723919 0.6056111 0.2812205 0.0000000 0.1397719
5 0.2840514 0.7222195 0.2851688 0.1397719 0.0000000

Le principe, pour initialiser l’algorithme est de tirer au hasard un trajet, c’est à dire une permutation de nos points,

>  (o = sample(1:n))
[1] 2 4 1 5 3

Il faut alors calculer le 4 distances entre ces 5 villes (on ne prévoit pas de retour à la ville de départ, pour l’instant). Le plus simple est de jouer avec cette matrice de la réarranger,

>  d[o[1:(n-1)],][,o[2:n]]
4         1         5         3
2 0.6056111 0.6186980 0.7222195 0.8399243
4 0.0000000 0.1723919 0.1397719 0.2812205
1 0.1723919 0.0000000 0.2840514 0.2258786
5 0.1397719 0.2840514 0.0000000 0.2851688

et de prendre les 4 valeurs sur la diagonale

>  diag(d[o[1:(n-1)],][,o[2:n]])
[1] 0.6056111 0.1723919 0.2840514 0.2851688

Ce sont les 4 distances que l’on cherche, entre les 5 villes. La distance totale parcourue est alors

>  sum(diag(d[o[1:(n-1)],][,o[2:n]]))
[1] 1.347223

L’idée est ensuite de permuter deux des villes

>  (i=sample(1:n,2))
[1] 3 1

(ici la première et la troisième de la liste, mais on pourrait décider de fixer la première si on doit partir de là) et de voir si la nouvelle distance est plus courte, ou pas,

>  os=o
>  os[i]=o[rev(i)]
>  sum(diag(d[os[1:(n-1)],][,os[2:n]]))
[1] 1.785391

La distance étant plus longue, on ne va pas permuter. On pourrait tenter 1000 fois de permuter deux villes. Si on ne parvient pas à réduire la distance totale, c’est qu’on a atteint une valeur intéressante (à défaut de pouvoir prouver qu’elle est optimale). Et si on peut améliorer la distance, on permute. Le plus simple est donc de faire des petites fonctions. La première reprend le calcul de la distance totale, sur la diagonale de la matrice permutée, à matrice de distance donnée et à permutation donnée,

> tsp.longueur=function(matrice,ordres) {
+ n=length(ordres)
+ sum(diag(matrice[ordres[1:(n-1)],][,ordres[2:n]])) }

 La seconde est juste une boucle,

> tsp.optimal=function(d,N=1000) {
+ d=as.matrix(d)
+ n=ncol(d)
+ o=sample(1:n)
+ v=tsp.longueur(d,o)
+ k=0
+ while(k < N) {
+   i<-sample(1:n,2)
+   os=o
+   os[i]=o[rev(i)]
+   w=tsp.longueur(d,os)
+   if(w < v) {
+     v=w
+     o=os
+     k=0} else {k=k+1}}
+ list(ordre=o,longueur=v)}

On peut se lancer maintenant, avec par exemple 30 points simulés, sur le carré unité,

>  set.seed(1)
>  n=30
>  x=matrix(runif(2*n),nr=n)
>  o=sample(1:n)
>  plot(x,xlim=0:1,ylim=0:1,xlab="",ylab="")
>  lines(x[o,],col='blue')
>  r=tsp.optimal(dist(x))
>  os=r$ordre
>  lines(x[os,],col='blue',lwd=1.5)

avec la première trajectoire, à gauche, qui reste en trait fin sur le dessin de droite,

Pas mal, non ? Bon, en fait le soucis est que si on joue plusieurs fois avec notre fonction, on obtient toujours un graphique différent,

Autrement dit, je trouve (certes) des chemins pas trop long, mais je suis loin de trouver le chemin optimal. Et histoire de clotûrer le débat, notons qu’on peut faire un peu de simulation, afin de comparer la distance obtenue en faisant le trajet au hasard, et un trajet “optimal“,

> RATIO=HASARD=OPTIMAL=matrix(NA,500,10)
> for(m in 1:10){
+ n=5*m
+ for(s in 1:500){
+  x=matrix(runif(2*n),nr=n)
+  r=tsp.opt(dist(x))
+  HAZARD[s,m]=tsp.longueur(as.matrix(dist(x)),1:n)
+  OPTIMAL[s,m]=r$longueur
+  RATIO[s,m]=HAZARD[s,m]/
+             OPTIMAL[s,m]
+ }
+}

Si on compare la longueur moyenne du trajet (sur 500 scénarios) avec un tour au hasard (en rouge) et ou tour optimisé (en bleu), on a

ce qui donne un ratio moyen (ou une moyenne des ratios, les deux sont représentés)

Autrement dit, avec une dizaine de points éparpillés (au hasard) dans le carré unité, un tour au hasard ne sera que deux fois plus long… étonnant non ?

Mais on doit pouvoir aller plus loin, parce que sous R, il y a des fonctions (et des packages) pour faire des algorithmes plus compliqués… à suivre avant la tournée des bonbons de demain soir !