Réflexions autour du voyage de commerce (TSP)

Le problème du voyageur de commerce consiste à trouver le plus court chemin passant par toutes les villes. On parle aussi de circuit hamiltonien qui consiste à trouver le plus court chemin passant par tous les noeuds d’un graphe. Ce programme explore quelques solutions approchées et intuitives.

Ce problème est NP-complet à savoir qu’il n’existe pas d’algorithme qui permette de trouver la solution avec un coût polynômial. C’est aussi un problème différent du plus court chemin dans un graphe qui consiste à trouver le plus court chemin reliant deux noeuds d’un graphe (mais pas forcément tous les noeuds de ce graphe).

Des villes tirées au hasard

import random
import matplotlib.pyplot as plt

n = 30
x = [random.random() for _ in range(n)]
y = [random.random() for _ in range(n)]

plt.plot(x, y, "o")

[<matplotlib.lines.Line2D object at 0x7fa94763ab90>]

Un parcours aléatoire de tous les noeuds de graphe donnera quelque chose de très éloigné de la solution optimale :

plt.plot(x + [x[0]], y + [y[0]], "o-")

[<matplotlib.lines.Line2D object at 0x7fa96217dde0>]

Croisements

La première constation est que le chemin ne peut pas être optimal car des arcs se croisent. On en déduit qu’une façon d’améliorer ce chemin est de décroiser certaines parties. On peut par exemple choisir deux points au hasard, retourner la partie du chemin au milieu de ces deux points et voir si la longueur du chemin s’en trouve diminuée. On peut également parcourir toutes les paires de noeuds possibles. C’est ce qui est implémenté ci-dessous.

def longueur(x, y, ordre):
    i = ordre[-1]
    x0, y0 = x[i], y[i]
    d = 0
    for o in ordre:
        x1, y1 = x[o], y[o]
        d += (x0 - x1) ** 2 + (y0 - y1) ** 2
        x0, y0 = x1, y1
    return d


ordre = list(range(len(x)))
print("longueur initiale", longueur(x, y, ordre))

longueur initiale 11.292330591173462

Permutations.

def permutation(x, y, ordre):
    d = longueur(x, y, ordre)
    d0 = d + 1
    it = 1
    while d < d0:
        it += 1
        print("iteration", it, "d=", d)
        d0 = d
        for i in range(0, len(ordre) - 1):
            for j in range(i + 2, len(ordre)):
                r = ordre[i:j].copy()
                r.reverse()
                ordre2 = ordre[:i] + r + ordre[j:]
                t = longueur(x, y, ordre2)
                if t < d:
                    d = t
                    ordre = ordre2
    return ordre


ordre = permutation(x, y, list(range(len(x))))
print("longueur min", longueur(x, y, ordre))
xo = [x[o] for o in ordre + [ordre[0]]]
yo = [y[o] for o in ordre + [ordre[0]]]
plt.plot(xo, yo, "o-")

iteration 2 d= 11.292330591173462
iteration 3 d= 1.341397079008917
iteration 4 d= 0.9548761833186405
iteration 5 d= 0.8631425458159007
iteration 6 d= 0.8456272108828561
longueur min 0.8456272108828561

[<matplotlib.lines.Line2D object at 0x7fa94c0d1180>]

Voilà qui est mieux. Maintenant, supposons que nous faisons une erreur lors du calcul de la distance : nous oublions le dernier arc qui boucle le chemin du dernier noeud au premier.

def longueur(x, y, ordre):
    # on change cette fonction
    d = 0
    for i in range(1, len(ordre)):
        n = ordre[i - 1]
        o = ordre[i]
        x0, y0 = x[n], y[n]
        x1, y1 = x[o], y[o]
        d += (x0 - x1) ** 2 + (y0 - y1) ** 2
    return d


ordre = list(range(len(x)))
print("longueur initiale", longueur(x, y, ordre))

longueur initiale 11.137631508007717

Et graphiquement.

ordre = permutation(x, y, list(range(len(x))))
print("longueur min", longueur(x, y, ordre))
xo = [x[o] for o in ordre]
yo = [y[o] for o in ordre]
plt.plot(xo, yo, "o-")

iteration 2 d= 11.137631508007717
iteration 3 d= 1.8052001053865456
iteration 4 d= 0.9739961042920614
iteration 5 d= 0.8147617077948388
longueur min 0.8147617077948388

[<matplotlib.lines.Line2D object at 0x7fa94bfc29b0>]

Noeud de départ constant

Jusque ici, tout concorde. Le chemin est plus court en ce sens qu’il oublie délibérément l’arc de bouclage que l’algorithme a tendance à choisir grand. Pour gagner du temps de calcul, un développeur se dit que le noeud de départ peut être constant. Après tout, le chemin est une boucle, elle passera toujours par le premier noeud. Qu’il soit en première position ne change rien et puis inverser une moitié, c’est équivalent à inverser l’autre moitié. On fait donc juste une modification :

def longueur(x, y, ordre):
    i = ordre[-1]
    x0, y0 = x[i], y[i]
    d = 0
    for o in ordre:
        x1, y1 = x[o], y[o]
        d += (x0 - x1) ** 2 + (y0 - y1) ** 2
        x0, y0 = x1, y1
    return d


ordre = list(range(len(x)))
print("longueur initiale", longueur(x, y, ordre))


def permutation(x, y, ordre):
    d = longueur(x, y, ordre)
    d0 = d + 1
    it = 1
    while d < d0:
        it += 1
        print("iteration", it, "d=", d, "ordre[0]", ordre[0])
        d0 = d
        for i in range(
            1, len(ordre) - 1
        ):  # on part de 1 et plus de 0, on est sûr que le premier noeud ne bouge pas
            for j in range(i + 2, len(ordre)):
                r = ordre[i:j].copy()
                r.reverse()
                ordre2 = ordre[:i] + r + ordre[j:]
                t = longueur(x, y, ordre2)
                if t < d:
                    d = t
                    ordre = ordre2
    return ordre


ordre = permutation(x, y, list(range(len(x))))
print("longueur min", longueur(x, y, ordre))
xo = [x[o] for o in ordre + [ordre[0]]]
yo = [y[o] for o in ordre + [ordre[0]]]
plt.plot(xo, yo, "o-")
plt.text(xo[0], yo[0], "0", color="r", weight="bold", size="x-large")
plt.text(xo[-2], yo[-2], "N-1", color="r", weight="bold", size="x-large")

longueur initiale 11.292330591173462
iteration 2 d= 11.292330591173462 ordre[0] 0
iteration 3 d= 1.740352647987561 ordre[0] 0
iteration 4 d= 1.0081199538507615 ordre[0] 0
iteration 5 d= 0.9966224243230848 ordre[0] 0
longueur min 0.9966224243230848

Text(0.8728691050558798, 0.522924291603524, 'N-1')

Le résultat attendu n’est pas celui qu’on observe. Est-ce une erreur d’implémentation ou une erreur de raisonnement ? J’étais pourtant sûr que mon raisonnement était correct et j’aurais tort d’en douter. C’est une erreur d’implémentation. Lorsqu’on``for j in range(i+2,len(ordre)):`` et r = ordre[i:j].copy(), on écrit que j va de i+2 inclus à len(ordre) exclu. Puis lorsqu’on écrit ordre[i:j], l’indice j est exclu ! Autrement dit, dans cette implémentation, le premier noeud et le dernier noeud ne bougeront jamais ! On s’empresse de corriger cela.

ordre = list(range(len(x)))
print("longueur initiale", longueur(x, y, ordre))


def permutation(x, y, ordre):
    d = longueur(x, y, ordre)
    d0 = d + 1
    it = 1
    while d < d0:
        it += 1
        print("iteration", it, "d=", d, "ordre[0]", ordre[0])
        d0 = d
        for i in range(
            1, len(ordre) - 1
        ):  # on part de 1 et plus de 0, on est sûr que le premier noeud ne bouge pas
            for j in range(i + 2, len(ordre) + 1):  # correction !
                r = ordre[i:j].copy()
                r.reverse()
                ordre2 = ordre[:i] + r + ordre[j:]
                t = longueur(x, y, ordre2)
                if t < d:
                    d = t
                    ordre = ordre2
    return ordre


ordre = permutation(x, y, list(range(len(x))))
print("longueur min", longueur(x, y, ordre))
xo = [x[o] for o in ordre + [ordre[0]]]
yo = [y[o] for o in ordre + [ordre[0]]]
plt.plot(xo, yo, "o-")
plt.text(xo[0], yo[0], "0", color="r", weight="bold", size="x-large")
plt.text(xo[-2], yo[-2], "N-1", color="r", weight="bold", size="x-large")

longueur initiale 11.292330591173462
iteration 2 d= 11.292330591173462 ordre[0] 0
iteration 3 d= 1.7214335447427036 ordre[0] 0
iteration 4 d= 1.2431331046193033 ordre[0] 0
iteration 5 d= 1.118929138241379 ordre[0] 0
iteration 6 d= 0.9399537486282151 ordre[0] 0
iteration 7 d= 0.848220111125475 ordre[0] 0
iteration 8 d= 0.8297080670435762 ordre[0] 0
longueur min 0.8297080670435762

Text(0.23407436285281968, 0.23141818184855922, 'N-1')

Pas parfait mais conforme à nos attentes (les miennes en tout cas) ! Soit dit en passant, la première version de l’algorithme laissait déjà le dernier noeud inchangé.

Un peu d’aléatoire en plus

La solution n’est pas parfaite en ce sens que visuellement, on voit que certaines parties du chemin pourraient être facilement améliorées. Mais si la solution était parfaite en toute circonstance, nous aurions trouvé un algorithme à temps polynômial ce qui est impossible. Dans notre cas, l’algorithme produit toujours la même solution car il parcourt les noeuds toujours dans le même sens. Un peu d’aléa devrait l’aider à trouver de meilleures solutions après quelques essais.

# In[8]:


ordre = list(range(len(x)))
print("longueur initiale", longueur(x, y, ordre))


def permutation_rnd(x, y, ordre):
    d = longueur(x, y, ordre)
    d0 = d + 1
    it = 1
    while d < d0:
        it += 1
        print("iteration", it, "d=", d, "ordre[0]", ordre[0])
        d0 = d
        for i in range(1, len(ordre) - 1):
            for j in range(i + 2, len(ordre) + 1):
                ik = random.randint(1, len(ordre) - 1)
                il = random.randint(ik + 1, len(ordre))
                r = ordre[ik:il].copy()
                r.reverse()
                ordre2 = ordre[:ik] + r + ordre[il:]
                t = longueur(x, y, ordre2)
                if t < d:
                    d = t
                    ordre = ordre2
    return ordre


ordre = permutation_rnd(x, y, list(range(len(x))))
print("longueur min", longueur(x, y, ordre))
xo = [x[o] for o in ordre + [ordre[0]]]
yo = [y[o] for o in ordre + [ordre[0]]]
plt.plot(xo, yo, "o-")
plt.text(xo[0], yo[0], "0", color="r", weight="bold", size="x-large")
plt.text(xo[-2], yo[-2], "N-1", color="r", weight="bold", size="x-large")

longueur initiale 11.292330591173462
iteration 2 d= 11.292330591173462 ordre[0] 0
iteration 3 d= 1.8476621280516345 ordre[0] 0
iteration 4 d= 1.3199871306554722 ordre[0] 0
iteration 5 d= 0.8810707833443141 ordre[0] 0
iteration 6 d= 0.8575701970529361 ordre[0] 0
iteration 7 d= 0.8571236004366429 ordre[0] 0
longueur min 0.8571236004366429

Text(0.3037264087729815, 0.03717707930068548, 'N-1')

Ca a l’air de marcher un peu mieux mais quelques aberrations car l’aléatoire n’est pas un parcours systématique de toutes les pairs. Par conséquent, il peut rester des croisements :

ordre = permutation_rnd(x, y, list(range(len(x))))
print("longueur min", longueur(x, y, ordre))
xo = [x[o] for o in ordre + [ordre[0]]]
yo = [y[o] for o in ordre + [ordre[0]]]
plt.plot(xo, yo, "o-")
plt.text(xo[0], yo[0], "0", color="r", weight="bold", size="x-large")
plt.text(xo[-2], yo[-2], "N-1", color="r", weight="bold", size="x-large")

iteration 2 d= 11.292330591173462 ordre[0] 0
iteration 3 d= 1.913964651299246 ordre[0] 0
iteration 4 d= 1.294046515504576 ordre[0] 0
iteration 5 d= 1.2023764900880596 ordre[0] 0
iteration 6 d= 1.1982004152982153 ordre[0] 0
iteration 7 d= 1.187839655906766 ordre[0] 0
iteration 8 d= 1.1639354407467406 ordre[0] 0
longueur min 1.1639354407467406

Text(0.5484190781132359, 0.363880751818936, 'N-1')

Pour éviter cela, on peut imposer un nombre d’itérations minimum et recommencer plusieurs à partir d’ordre initiaux aléatoires :

def permutation_rnd(x, y, ordre, miniter):
    d = longueur(x, y, ordre)
    d0 = d + 1
    it = 1
    while d < d0 or it < miniter:
        it += 1
        d0 = d
        for i in range(1, len(ordre) - 1):
            for j in range(i + 2, len(ordre) + 1):
                ik = random.randint(1, len(ordre) - 1)
                il = random.randint(ik + 1, len(ordre))
                r = ordre[ik:il].copy()
                r.reverse()
                ordre2 = ordre[:ik] + r + ordre[il:]
                t = longueur(x, y, ordre2)
                if t < d:
                    d = t
                    ordre = ordre2
    return ordre


def n_permutation(x, y, miniter):
    ordre = list(range(len(x)))
    bordre = ordre.copy()
    d0 = longueur(x, y, ordre)
    for i in range(0, 20):
        print("iteration", i, "d=", d0)
        random.shuffle(ordre)
        ordre = permutation_rnd(x, y, ordre, 20)
        d = longueur(x, y, ordre)
        if d < d0:
            d0 = d
            bordre = ordre.copy()
    return bordre

La distance initiale.

ordre = list(range(len(x)))
print("longueur initiale", longueur(x, y, ordre))

longueur initiale 11.292330591173462

La longueur obtenue.

ordre = n_permutation(x, y, 20)
print("longueur min", longueur(x, y, ordre))
xo = [x[o] for o in ordre + [ordre[0]]]
yo = [y[o] for o in ordre + [ordre[0]]]
plt.plot(xo, yo, "o-")
plt.text(xo[0], yo[0], "0", color="r", weight="bold", size="x-large")
plt.text(xo[-2], yo[-2], "N-1", color="r", weight="bold", size="x-large")


# C'est mieux.

iteration 0 d= 11.292330591173462
iteration 1 d= 0.8669800334772021
iteration 2 d= 0.8669800334772021
iteration 3 d= 0.8669800334772021
iteration 4 d= 0.8669800334772021
iteration 5 d= 0.8669800334772021
iteration 6 d= 0.8669800334772021
iteration 7 d= 0.8669800334772021
iteration 8 d= 0.8449800420196413
iteration 9 d= 0.8449800420196413
iteration 10 d= 0.8449800420196413
iteration 11 d= 0.8449800420196413
iteration 12 d= 0.8449800420196413
iteration 13 d= 0.8449800420196413
iteration 14 d= 0.8237571704989521
iteration 15 d= 0.822760461350098
iteration 16 d= 0.822760461350098
iteration 17 d= 0.822760461350098
iteration 18 d= 0.822760461350098
iteration 19 d= 0.822760461350098
longueur min 0.822760461350098

Text(0.9551541397429326, 0.4454818330167861, 'N-1')