Risolvere i problemi del commesso viaggiatore geografico utilizzando Python

Python per risolvere il problema del commesso viaggiatore geografico

Utilizzare pyconcorde per trovare soluzioni ottimali a problemi di routing del mondo reale

Il famoso problema del commesso viaggiatore (TSP) consiste nel trovare un percorso ottimale tra una collezione di nodi (città) e tornare al punto di partenza. Sembra semplice, ma è impossibile risolverlo con la forza bruta per un elevato numero di nodi, poiché il numero di possibili ordinamenti di n città è n!. Ciò significa che anche solo per 30 città, il numero di viaggi da controllare è di 265.252.859.812.191.058.636.308.480.000.000. I problemi TSP di grandi dimensioni sono impraticabili da risolvere con la forza bruta, anche per computer potenti.

Fortunatamente, sono stati sviluppati alcuni algoritmi che riducono drasticamente il calcolo necessario per risolvere grandi TSP. Un software del genere, Concorde, è stato sviluppato qualche decennio fa per l’uso nella comunità accademica. Sebbene sia piuttosto tecnico da usare come software autonomo e sia destinato solo a specialisti, è stato sviluppato pyconcorde come wrapper Python per Concorde. Una spiegazione degli algoritmi utilizzati in Concorde è al di fuori dello scopo di questo articolo. Tuttavia, approfondiremo il codice necessario per riprodurre questi problemi e le relative soluzioni in Python.

TSP geografico del mondo reale

Come si potrebbe risolvere un problema del commesso viaggiatore del mondo reale, geografico? I punti del mondo reale non sono collegati da semplici linee bidimensionali come nell’immagine sopra. Invece, le caratteristiche geografiche sono collegate da varie possibili rotte, e queste rotte cambieranno a seconda che qualcuno stia camminando, andando in bicicletta o guidando.

Perché un data scientist o un ingegnere del software potrebbe voler risolvere un TSP del mondo reale? Ecco alcuni esempi di casi d’uso:

• IA generativa e MLOps una potente combinazione per lo sviluppo efficiente ed efficace dell’IA
• Un nuovo modo di guardare alla privacy dei dati
• Liste Python vs Array NumPy Un’analisi approfondita della disposizione in memoria e dei vantaggi delle prestazioni

1. Un’azienda che impiega corrieri ha bisogno di un modo per calcolare percorsi ottimali attraverso una città, riducendo al minimo il tempo trascorso sulla strada per ciascun autista.
2. Un tour operator deve trovare il percorso più breve che collega un insieme di destinazioni, entro un determinato periodo di tempo.
3. Un’azienda di smaltimento rifiuti o un’autorità locale deve assegnare le risorse per garantire che i ritiri siano ordinati nel modo più efficiente possibile.

Per risolvere un TSP del mondo reale, è possibile utilizzare la libreria routingpy per trovare percorsi, distanze (in metri) e durate (in secondi) tra punti geografici in coppie [longitudine, latitudine]. In questo articolo descriveremo il metodo che può essere utilizzato per un tale problema.

Percorso di codifica

Qui viene descritta una guida per risolvere un TSP geografico utilizzando Python. La struttura generale del processo di risoluzione del problema è la seguente:

1. Ottenere un elenco di n coordinate come coppie [longitudine, latitudine].
2. Utilizzare un servizio di routing per ottenere una matrice (n x n) di durate reali del mondo reale tra ciascuna di queste coordinate, per il profilo appropriato (a piedi, in bicicletta, guidando un’auto, guidando un veicolo pesante, ecc.). Questa matrice sarà asimmetrica (guidare da A a B non è esattamente l’inverso di B a A).
3. Trasformare la matrice (n x n) in una matrice simmetrica (2n x 2n).
4. Inserire questa matrice nel risolutore Concorde per trovare un ordinamento ottimale delle coordinate.
5. Creare il percorso del mondo reale utilizzando il servizio di routing.
6. Visualizzare i risultati su una mappa.
7. Facoltativamente, creare un file GPX del percorso finale.

Ogni di questi passaggi verrà trattato in dettaglio.

Passaggio 1: Ottenere le coordinate

Per il nostro esempio, considereremo il problema di guidare una macchina tra 79 città nel Regno Unito. Nella mappa sottostante sono mostrate le città del Regno Unito in blu. Un data scientist può trovare le coordinate in vari modi. Se necessario, possono essere trovate manualmente utilizzando Google Maps o Google Earth.

La struttura del codice e i dati utilizzati in questo esempio sono disponibili anche in questo repository di GitHub.

Ecco un file CSV che contiene le coordinate delle città (gb_cities.csv nel repository), e di seguito il codice necessario per importarle utilizzando pandas.

Nome Luogo,Latitudine,LongitudineAberdeen,57.149651,-2.099075Ayr,55.458565,-4.629179Basildon,51.572376,0.470009Bath,51.380001,-2.36Bedford,52.136436,-0.460739...

import pandas as pddf = pd.read_csv('gb_cities.csv')coordinates = df[['Longitudine', 'Latitudine']].valuenames = df['Nome Luogo'].values

Passaggio 2: Utilizzare un servizio di routing per ottenere la matrice della durata

Sono disponibili diversi servizi di routing tramite la libreria routingpy. L’API di Graphhopper include un livello gratuito che consente un utilizzo limitato dalla velocità. Altri router disponibili tramite routingpy sono elencati nella documentazione.

import routingpy as rpimport numpy as npapi_key = # ottenere una chiave gratuita su https://www.graphhopper.com/api = rp.Graphhopper(api_key=api_key)matrix = api.matrix(locations=coordinates, profile='car')durations = np.matrix(matrix.durations)print(durations)

Ecco durations, una matrice 79 x 79 del tempo di guida in secondi tra le coordinate:

matrix([[    0, 10902, 30375, ..., 23380, 25233, 19845],        [10901,     0, 23625, ..., 16458, 18312, 13095],        [30329, 23543,     0, ...,  8835,  9441, 12260],        ...,        [23397, 16446,  9007, ...,     0,  2789,  7924],        [25275, 18324,  9654, ...,  2857,     0,  9625],        [19857, 13071, 12340, ...,  8002,  9632,     0]])

Il tempo di guida tra le città può essere determinato come segue:

1. Ogni riga e colonna corrisponde a una città: Aberdeen è la prima riga e colonna, Ayr è la seconda, Basildon è la terza, e così via.
2. Per trovare il tempo tra Aberdeen e Ayr, guarda la 1a riga, 2a colonna: 10.902 secondi. Il tempo inverso (Ayr verso Aberdeen) è di 10.901 secondi.
3. In generale, il tempo dalla città i-esima alla città j-esima è all’intersezione tra la i-esima riga e la j-esima colonna.

Si noti che la matrice, come previsto, ha zeri lungo la diagonale, poiché ciascun punto è collegato a se stesso con una distanza o durata zero. Inoltre, la matrice non è del tutto simmetrica: le durate di guida tra le città sono probabilmente diverse in direzioni opposte, a causa di diverse disposizioni stradali e punti critici del traffico. Tuttavia, sono globalmente simili, come ci si aspetterebbe.

Passaggio 3: Trasformare la matrice asimmetrica in simmetrica

Prima di utilizzare questa matrice per generare un ordinamento ottimale in pyconcorde, è necessario rendere la matrice simmetrica. Un metodo per trasformare un TSP asimmetrico in TSP simmetrico viene descritto da Jonker e Volgenant (1983): Transforming asymmetric into symmetric traveling salesman problems, Operations Research Letters, 2(4), 161–163. La teoria di questa trasformazione è descritta di seguito. Se desiderato, questa sezione può essere saltata (scorri fino alla sezione intitolata Transforming the geographic asymmetric TSP).

Trasformazione asimmetrica-simmetrica di Jonker/Volgenant

Di seguito è visualizzata un’istanza di TSP asimmetrico con 3 nodi e la sua matrice delle distanze.

matrix([[0, 5, 2],        [7, 0, 4],        [3, 4, 0]])

Ecco uno schema del metodo utilizzato per trasformarlo in un TSP simmetrico:

1. Crea nuovi nodi fantasma, A’, B’ e C’. Collega A ad A’, B a B’ e C a C’ con distanza zero.
2. Collega i nodi con i pesi come segue: A a B è ora rappresentato da A’ a B; B a A è ora B’ a A. B a C è ora B’ a C; C a B è ora C’ a B. C a A è ora C’ a A; A a C è A’ a C.
3. Imposta tutti gli altri pesi degli archi come infiniti, in modo che qualsiasi algoritmo non tenti di viaggiare tra di essi. Poiché questo sarà impraticabile in seguito utilizzando pyconcorde, invece imposta tutti gli altri pesi molto più alti del peso più alto che abbiamo. In questo caso, li impostiamo tutti a 99.

Ecco la matrice delle distanze risultante. L’ordinamento dei nodi nella matrice è: A, B, C, A’, B’, C’.

matrix([[ 0, 99, 99,  0,  7,  3],        [99,  0, 99,  5,  0,  4],        [99, 99,  0,  2,  4,  0],        [ 0,  5,  2,  0, 99, 99],        [ 7,  0,  4, 99,  0, 99],        [ 3,  4,  0, 99, 99,  0]])

Si noti nuovamente che la diagonale è zero, come ci si aspetterebbe, e che la matrice è ora simmetrica. La matrice originale si trova nell’angolo in basso a sinistra della nuova matrice, e la sua trasposta si trova nell’angolo in alto a destra. Nel frattempo, le parti in alto a sinistra e in basso a destra contengono pesi molto alti tra i nodi.

A, B e C (in alto a sinistra) non sono più collegati tra loro (più precisamente, sono collegati ma con pesi molto alti invece di infinito, per scopi pratici). Ciò significa che qualsiasi algoritmo non cercherà di trovare un percorso tra questi nodi. Allo stesso modo, A’, B’ e C’ (in basso a destra) non sono collegati tra loro. Invece, la natura direzionale della rete asimmetrica originale è rappresentata qui dai pesi sui nodi originali A, B e C, insieme ai loro nodi fantasma A’, B’ e C’.

Esiste una corrispondenza uno a uno tra le soluzioni del problema asimmetrico originale e il nuovo TSP simmetrico:

• A — B — C — A corrisponde a A — A’ — B — B’ — C — C’ — A
• A — C — B — A corrisponde a A — A’ — C — C’ — B — B’ — A

In entrambi i casi, i nodi fantasma A’, B’ e C’ si alternano con i nodi originali A, B e C, e ogni nodo originale è adiacente al suo nodo fantasma “partner” (A è adiacente ad A’, e così via).

Trasformazione del TSP asimmetrico geografico

Tornando al nostro esempio pratico, possiamo creare una funzione per trasformare una matrice TSP asimmetrica in una simmetrica:

def symmetricize(m, high_int=None):        # se high_int non è fornito, rendilo uguale a 10 volte il valore massimo:    if high_int is None:        high_int = round(10*m.max())            m_bar = m.copy()    np.fill_diagonal(m_bar, 0)    u = np.matrix(np.ones(m.shape) * high_int)    np.fill_diagonal(u, 0)    m_symm_top = np.concatenate((u, np.transpose(m_bar)), axis=1)    m_symm_bottom = np.concatenate((m_bar, u), axis=1)    m_symm = np.concatenate((m_symm_top, m_symm_bottom), axis=0)        return m_symm.astype(int) # Concorde richiede pesi interi

symmetricize(durations) restituisce:

matrix([[     0, 461120, 461120, ...,  23397,  25275,  19857],        [461120,      0, 461120, ...,  16446,  18324,  13071],        [461120, 461120,      0, ...,   9007,   9654,  12340],        ...,        [ 23397,  16446,   9007, ...,      0, 461120, 461120],        [ 25275,  18324,   9654, ..., 461120,      0, 461120],        [ 19857,  13071,  12340, ..., 461120, 461120,      0]])

Questa matrice 158 x 158 contiene una copia di durations nella parte inferiore sinistra e una copia trasposta nella parte superiore destra. Il valore elevato di 461.120 (10 volte il valore massimo in durations) significa che, per scopi pratici, i nodi con questa durata non sono collegati.

Questa matrice può infine essere fornita a pyconcorde per calcolare un percorso ottimale.

Passaggio 4: Utilizzare il risolutore Concorde

Installazione di pyconcorde

Eseguire i seguenti comandi per installare pyconcorde (l’installazione è disponibile in Linux o Mac OS, ma non in Windows al momento):

virtualenv venv                                  # crea ambienti virtualisource venv/bin/activate                         # attivalogit clone https://github.com/jvkersch/pyconcorde # clona il repocd pyconcorde                                    # cambia directorypip install -e .                                 # installa pyconcorde

Risoluzione del TSP in Python

Ora possiamo importare da concorde in uno script Python.

from concorde.problem import Problemfrom concorde.concorde import Concordedef solve_concorde(matrix):    problem = Problem.from_matrix(matrix)    solver = Concorde()    solution = solver.solve(problem)    print(f'Percorso ottimale: {solution.tour}')    return solution

La nostra matrice simmetrica delle durate può essere fornita a solve_concorde().

durations_symm = symmetricize(durations)solution = solve_concorde(durations_symm)

Ecco l’output della stampa:

Percorso ottimale: [0, 79, 22, 101, 25, 104, 48, 127, 68, 147, 23, 102, 58, 137, 7, 86, 39, 118, 73, 152, 78, 157, 36, 115, 42, 121, 62, 141, 16, 95, 20, 99, 51, 130, 40, 119, 19, 98, 59, 138, 50, 129, 54, 133, 27, 106, 10, 89, 4, 83, 66, 145, 33, 112, 14, 93, 2, 81, 45, 124, 32, 111, 11, 90, 29, 108, 34, 113, 24, 103, 8, 87, 17, 96, 56, 135, 64, 143, 61, 140, 75, 154, 52, 131, 71, 150, 18, 97, 3, 82, 9, 88, 74, 153, 55, 134, 72, 151, 28, 107, 12, 91, 70, 149, 65, 144, 35, 114, 31, 110, 77, 156, 63, 142, 41, 120, 69, 148, 6, 85, 76, 155, 67, 146, 15, 94, 44, 123, 47, 126, 60, 139, 57, 136, 38, 117, 13, 92, 5, 84, 43, 122, 49, 128, 46, 125, 21, 100, 1, 80, 30, 109, 53, 132, 37, 116, 26, 105]

Questa soluzione mostra l’ordinamento dei nodi nel tour ottimale. Nota che questa soluzione, come previsto in precedenza, contiene nodi originali (numerati da 0 a 78) alternati con i loro nodi fantasma partner (da 79 a 157):

• 0 è associato a 79,
• 22 a 101,
• 25 a 104, e così via…

Questo suggerisce che la soluzione ha funzionato correttamente.

Passo 5: Creazione del percorso nel mondo reale

Il prossimo passo è selezionare gli elementi alternati della soluzione (i nodi corrispondenti alle 79 città originali), quindi ordinare le coordinate di conseguenza.

# seleziona gli elementi alternativi: questi corrispondono ai nodi originalistour = solution.tour[::2]# ordina le coordinate originali e i nomicoords_ordered = [coordinates[i].tolist() for i in tour]names_ordered = [names[i] for i in tour]

Ecco i primi nomi delle città in names_ordered, (l’ordinamento reale delle città nel tour ottimale):

['Aberdeen', 'Dundee', 'Edinburgh', 'Newcastle Upon Tyne', 'Sunderland', 'Durham', ...]

Ora aggiungiamo nuovamente la prima città per ottenere un tour completo a circuito chiuso, e infine otteniamo il percorso finale utilizzando l’API delle indicazioni di Graphhopper.

# aggiungi nuovamente la prima città per un tour completo a circuito chiusocoords_ordered_return = coords_ordered + [coords_ordered[0]]# ottieni le indicazioni stradali complete per il percorso ordinatoirections = api.directions(locations=coords_ordered_return, profile='car')

Passo 6: Visualizzazione su una mappa

Vedere il percorso finale su una mappa ci permetterà di avere fiducia nel risultato, oltre a consentirci di utilizzare la soluzione in un contesto pratico. Il seguente codice mostrerà una mappa folium che può essere salvata in HTML.

import foliumdef generate_map(coordinates, names, directions):    # folium ha bisogno di lat, long    coordinates = [(y, x) for (x, y) in coordinates]    route_points = [(y, x) for (x, y) in directions.geometry]    lat_centre = np.mean([x for (x, y) in coordinates])    lon_centre = np.mean([y for (x, y) in coordinates])    centre = lat_centre, lon_centre    m = folium.Map(location=centre, zoom_start=1, zoom_control=False)    # plot the route line    folium.PolyLine(route_points, color='red', weight=2).add_to(m)        # plot each point with a hover tooltip      for i, (point, name) in enumerate(zip(coordinates, names)):        folium.CircleMarker(location=point,                      tooltip=f'{i}: {name}',                      radius=2).add_to(m)    custom_tile_layer = folium.TileLayer(        tiles='http://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}.png',        attr='CartoDB Positron',        name='Positron',        overlay=True,        control=True,        opacity=0.7  # Regola l'opacità per controllare il livello di sfocatura    )    custom_tile_layer.add_to(m)    folium.LayerControl().add_to(m)    sw = (np.min([x for (x, y) in coordinates]), np.min([y for (x, y) in coordinates]))    ne = (np.max([x for (x, y) in coordinates]), np.max([y for (x, y) in coordinates]))    m.fit_bounds([sw, ne])    return mgenerate_map(coords_ordered, names_ordered, directions).save('gb_cities.html')

Il risultato è mostrato nella parte superiore di questo articolo. Clicca qui per visualizzarlo come mappa interattiva. È possibile ingrandire la mappa per vedere più dettagli e passare il mouse sopra le singole città che riveleranno il loro numero nella sequenza del tour. Di seguito è mostrata una parte ingrandita della mappa che mostra il percorso che passa attraverso Sheffield (tra Lincoln e Chesterfield nel tour ottimale).

Passo 7: Opzionale: Creazione di un file GPX

Se il percorso calcolato deve essere seguito nella vita reale, ad esempio su un dispositivo con un GPS (come un telefono o un sistema di navigazione per auto), è possibile creare un file GPX. Questo non fa parte del problema di ottimizzazione, ma è un passaggio opzionale aggiuntivo disponibile se si desidera salvare il percorso in un file. Il file GPX viene creato dalla variabile directions:

def generate_gpx_file(directions, filename):    gpx_template = """<?xml version="1.0" encoding="UTF-8"?>    <gpx version="1.1" xmlns="http://www.topografix.com/GPX/1/1"        xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"        xsi:schemaLocation="http://www.topografix.com/GPX/1/1        http://www.topografix.com/GPX/1/1/gpx.xsd">        <trk>            <name>Traccia</name>            <trkseg>{}</trkseg>        </trk>    </gpx>    """    trkseg_template = """        <trkpt lat="{}" lon="{}"/>    """    trkseg_elements = ""    for point in directions.geometry:        trkseg_elements += trkseg_template.format(point[1], point[0])    gpx_data = gpx_template.format(trkseg_elements)    with open(filename, 'w') as file:        file.write(gpx_data)generate_gpx_file(directions, 'gb_cities.gpx')

È possibile scaricare qui il file GPX per questo problema.

Conclusioni

Abbiamo visto come possiamo combinare i seguenti elementi per risolvere problemi di commesso viaggiatore geografici del mondo reale:

1. Matrici di indicazioni e durata dalla libreria routingpy, specificando un appropriato profilo (mezzo di trasporto).
2. Risolutore Concorde efficiente e potente tramite l’incapsulamento pyconcorde, per fornire un percorso ottimale.
3. Visualizzazione utilizzando folium per creare una mappa.

Il tour in auto mostrato sopra è una soluzione convincente al problema del commesso viaggiatore con 79 città e, secondo il risolutore Concorde, è provabilmente ‘ottimale’. Tuttavia, poiché stiamo lavorando con dati del mondo reale, il risultato finale è valido solo quanto l’input. Ci affidiamo al fatto che la matrice delle durate punto-punto ottenuta da routingpy sia rappresentativa del mondo reale. In realtà, il tempo impiegato per camminare, pedalare o guidare tra i punti dipenderà dall’ora del giorno o dal giorno della settimana. Questa è una limitazione del metodo che abbiamo utilizzato. Un modo per avere maggiore fiducia nel risultato finale sarebbe quello di provare lo stesso metodo con un servizio di routing alternativo. Ogni servizio di routing (Graphhopper, ORS, Valhalla, ecc.) ha la propria API che potrebbe essere utilizzata per un problema di commesso viaggiatore come quello descritto qui, e i risultati potrebbero essere confrontati tra diversi servizi.

Nonostante le limitazioni del mondo reale nel risolvere un problema come questo, la metodologia sopra fornisce un buon punto di partenza per un venditore o corriere che ha bisogno di muoversi in città nel modo più efficiente possibile o per un turista che spera di visitare il maggior numero possibile di luoghi durante il suo tour. Visualizzando i risultati su una mappa interattiva e archiviando il percorso come file GPX, la soluzione è utile sia per l’utente finale che per il data scientist che ha implementato il codice.