Il problema del routing dei veicoli soluzioni esatte ed euristiche

Routing dei veicoli soluzioni esatte ed euristiche

Comprendere come risolvere problemi di routing complessi con Python

Il problema del routing dei veicoli (VRP) mira a determinare il miglior insieme di percorsi da effettuare da una flotta di veicoli per servire un determinato insieme di clienti. A causa delle sue numerose applicazioni e degli aspetti combinatori complessi, è uno dei problemi più studiati nella Ricerca Operativa e nell’Ottimizzazione Numerica.

Il problema del routing dei veicoli con capacità limitata (CVRP) è una delle varianti più comuni, poiché introduce veicoli con capacità di carico limitata e possibilmente vincoli di durata/distanza. Altre varianti comuni introducono anche depositi multipli, flotte eterogenee, prelievi e consegne e vincoli di finestra temporale.

Gli aspetti combinatori di questi problemi sono tali che considerando un semplice insieme di 15 punti, ci sono 6 × 10¹¹ percorsi possibili che li collegano (Dantzig & Ramser, 1959). Pertanto, alcune applicazioni del mondo reale rimarrebbero impraticabili fino a quando la ricerca computazionale e algoritmica non avanzerà negli ultimi decenni. Gli algoritmi di Branch-Cut-and-Price sono stati in grado di dimostrare l’ottimalità delle istanze CVRP con qualche centinaio di clienti (Fukasawa et al., 2006; Pecin et al., 2017), e le metaeuristiche all’avanguardia combinate con tecniche di ricerca locale possono fornire soluzioni di buona qualità (a volte ottimali) per queste istanze in pochi secondi (Vidal et al., 2012; Vidal, 2022).

In questo articolo, introdurremo il problema del routing dei veicoli con capacità (e durata) vincolate e lo risolveremo utilizzando la Programmazione Intera Mista (MIP) e algoritmi (meta)euristici specializzati. Nella prima parte, useremo il modulo Python AML Pyomo, con il risolutore HiGHS, mentre nella seconda parte utilizzeremo il pacchetto Google OR Tools.

• Cosa è l’Intelligenza Aziendale?
• Incontra LEVER un semplice approccio di intelligenza artificiale per migliorare la generazione di codice da linguaggio imparando a verificare i programmi generati con i loro risultati di esecuzione.
• Perché dovresti interessarti all’EU AI Act oggi

Il lettore interessato alle applicazioni del mondo reale più che agli aspetti teorici del problema potrebbe scorrere rapidamente la sezione MIP e prestare maggiore attenzione alle sezioni Specialized (Meta)Heuristics e Useful Extensions.

Coloro interessati a comprendere in dettaglio la formulazione MIP ma non ancora familiari con l’ottimizzazione numerica potrebbero trovare utile dare un’occhiata alle mie storie precedenti sulla Programmazione Lineare e sul metodo Branch & Bound prima di procedere con questa.

Come al solito, è possibile trovare il codice completo in questo repository git.

Ora, immergiamoci!

Programmazione Intera Mista

La formulazione matematica presentata in questa sezione utilizzerà le stesse equazioni presentate da Toth & Vigo (2002) nel modello lì definito come “formulazione del flusso di veicoli a tre indici”.

Consideriamo un insieme V di nodi (richieste e deposito) e un insieme K di veicoli. Utilizzeremo le lettere minuscole i e j per indicare gli indici dei nodi e la lettera minuscola k per indicare gli indici dei veicoli. Poiché questo modello è valido per il caso asimmetrico, assumiamo che i nodi facciano parte di un grafo diretto completo G(V, A) con archi A. In questo problema, c’è un singolo nodo deposito indicizzato da 0 e tutti i veicoli hanno la stessa capacità Q. Consideriamo due gruppi di variabili decisionali:

• x_{i, j, k}: È una variabile binaria che indica un arco attivo dal nodo i al nodo j eseguito dal veicolo k.
• y_{i, k}: È una variabile binaria che indica che la richiesta del nodo i viene soddisfatta dal veicolo k.

Consideriamo come obiettivo minimizzare il valore di costo associato agli archi attivi. La durata totale o la distanza sono esempi comuni. Supponiamo che il costo di attraversamento dell’arco i, j sia cᵢⱼ. La funzione obiettivo può essere formulata come segue.

Dobbiamo anche includere vincoli che garantiscono:

• Ogni cliente i viene visitato una volta, quindi ha un arco attivo che parte da esso e uno che arriva ad esso.
• Se una qualsiasi variabile dell’arco indicizzata dal veicolo k entra in un nodo i o ne esce, la richiesta q di questo nodo viene assegnata al veicolo k.
• La richiesta totale assegnata a un veicolo non deve superare la sua capacità Q.
• Esattamente |K| nodi partono dal deposito e arrivano al deposito.
• Non ci sono sottotour… Tuttavia, il numero di sottotour è potenzialmente troppo grande per essere enumerato dall’inizio. Approfondiremo ulteriormente come farlo.

Come di consueto per i tutorial Python, iniziamo la parte pratica importando le librerie utilizzate in questa sezione:

import timefrom itertools import cycleimport numpy as npfrom scipy.spatial.distance import pdist, squareformimport matplotlib.pyplot as pltimport matplotlib as mplimport networkx as nximport pyomo.environ as pyofrom pyomo.contrib.appsi.solvers.highs import Highs

E ora istanziamo un problema casuale con N nodi di domanda. In questo esempio, si suppone che il nodo del deposito sia il primo nodo (indice 0), quindi assicuriamo che anche la sua domanda corrispondente sia zero.

np.random.seed(42)  # I risultati dovrebbero essere sempre gli stessiN = 10domande = np.random.randint(1, 10, size=N)domande[0] = 0capacità = 15n_veicoli = 4coordinate = np.random.rand(N, 2)distanze = squareform(pdist(coordinate, metric="euclidean"))distanze = np.round(distanze, decimals=4)  # evitare errori numerici

Il numero di veicoli necessari può essere calcolato utilizzando un problema di imballaggio binario. Un esempio di come eseguirlo è incluso anche nel codice sorgente completo.

Ci sono due approcci per modellare un problema in pyomo: modelli astratti e concreti. Nel primo approccio, le espressioni algebriche del problema vengono definite prima di fornire alcuni valori dei dati, mentre nel secondo, l’istanza del modello viene creata immediatamente mentre i suoi elementi vengono definiti. Puoi trovare ulteriori informazioni su questi approcci nella documentazione della libreria o nel libro di Bynum et al. (2021). In questo articolo, adotteremo la formulazione del modello concreto.

modello = pyo.ConcreteModel()

Istanziare i set di nodi di domanda V, archi A e veicoli K. Nota che il nodo del deposito è incluso nell’insieme di nodi V come nella formulazione matematica originale.

modello.V = pyo.Set(initialize=range(len(domande)))modello.A = pyo.Set(initialize=[(i, j) for i in modello.V for j in modello.V if i != j])modello.K = pyo.Set(initialize=range(n_veicoli))

Ora i nostri parametri per capacità, carico di domanda e costi degli archi.

modello.Q = pyo.Param(initialize=capacità)modello.q = pyo.Param(modello.V, initialize={i: d for (i, d) in enumerate(domande)})modello.c = pyo.Param(modello.A, initialize={(i, j): distanze[i, j] for (i, j) in modello.A})

Dobbiamo anche includere i vincoli precedentemente elencati. Iniziamo implementandoli utilizzando la firma Pyomo usuale: funzione(modello, *dominio).

def archi_in(modello, i):    if i == modello.V.first():        return sum(modello.x[:, i, :]) == len(modello.K)    else:        return sum(modello.x[:, i, :]) == 1.0def archi_out(modello, i):    if i == modello.V.first():        return sum(modello.x[i, :, :]) == len(modello.K)    else:        return sum(modello.x[i, :, :]) == 1.0def assegnazione_veicolo(modello, i, k):    return sum(modello.x[:, i, k]) == modello.y[i, k]def assegnazione_veicolo_complementare(modello, i, k):    return sum(modello.x[i, :, k]) == modello.y[i, k]def vincolo_capacità(modello, k):    return sum(modello.y[i, k] * modello.q[i] for i in modello.V) <= modello.Q

E poi li incorporiamo nel nostro modello.

modello.archi_in = pyo.Constraint(modello.V, rule=archi_in)modello.archi_out = pyo.Constraint(modello.V, rule=archi_out)modello.assegnazione_veicolo = pyo.Constraint(modello.V, modello.K, rule=assegnazione_veicolo)modello.assegnazione_veicolo_complementare = pyo.Constraint(modello.V, modello.K, rule=assegnazione_veicolo_complementare)modello.vincolo_capacità = pyo.Constraint(modello.K, rule=vincolo_capacità)

Noto che non ho ancora incluso i vincoli di eliminazione dei sottocicli. Dovremmo considerare tutte le possibili permutazioni dei N nodi prendendo k alla volta e queste potrebbero essere proibitivamente grandi da enumerare anche per istanze di dimensioni moderate. In alternativa, nella nostra procedura di soluzione, includeremo in modo ricorsivo i vincoli di eliminazione dei sottocicli ogni volta che viene trovata una nuova soluzione se verifichiamo che questa soluzione produca sottocicli. In alcuni solver commerciali, questi sono chiamati “vincoli pigri” e possono essere incorporati direttamente nel solver tramite callback.

Per prima cosa creiamo una funzione che, data una sottotour, tutti i nodi rimanenti, un nodo dalla sottotour e un veicolo, restituisce un’espressione Pyomo corrispondente alla formulazione matematica precedentemente indicata. Inoltre, includiamo un ConstraintList a cui aggiungeremo nuovi elementi man mano che procediamo con la soluzione.

def subtour_elimination(model, S, Sout, h, k):    nodes_out = sum(model.x[i, j, k] for i in S for j in Sout)    return model.y[h, k] <= nodes_outmodel.subtour_elimination = pyo.ConstraintList()

Dobbiamo creare alcune funzioni che, dati una soluzione, restituiscano sottotour creati (se esistono). Per farlo, creeremo prima una lista di archi attivi nel modello utilizzando la funzione find_arcs. Questa lista verrà utilizzata per creare un grafo diretto incompleto utilizzando la classe DiGraph di Networkx. La funzione find_subtours dovrebbe restituire una lista di insiemi di componenti collegate.

def find_arcs(model):    arcs = []    for i, j in model.A:        for k in model.K:            if np.isclose(model.x[i, j, k].value, 1):                arcs.append((i, j))    return arcsdef find_subtours(arcs):    G = nx.DiGraph(arcs)    subtours = list(nx.strongly_connected_components(G))    return subtours

Il nostro obiettivo è eliminare gruppi di componenti collegate che non includono il nodo deposito. Quindi nel prossimo passaggio, creeremo funzioni che iterano sulla lista di insiemi e includono nuovi vincoli se l’insieme di componenti non include il nodo deposito. Questo utilizzerà il metodo add della classe ConstraintList.

def eliminate_subtours(model, subtours):    proceed = False    for S in subtours:        if 0 not in S:            proceed = True            Sout = {i for i in model.V if i not in S}            for h in S:                for k in model.K:                    model.subtour_elimination.add(                      subtour_elimination(model, S, Sout, h, k)                    )    return proceed

E ora abbiamo tutto pronto per proporre una procedura di soluzione che risolve iterativamente l’MIP, verifica se la soluzione corrente ha sottotour e, in caso affermativo, include nuovi vincoli per eliminarli. In caso contrario, la soluzione trovata è ottimale.

def solve_step(model, solver):    sol = solver.solve(model)    arcs = find_arcs(model)    subtours = find_subtours(arcs)    time.sleep(0.1)    proceed = eliminate_subtours(model, subtours)    return sol, proceed def solve(model, solver):    proceed = True    while proceed:        sol, proceed = solve_step(model, solver)    return sol

Ora istanziamo il solver e risolviamo il nostro modello. Il solver Highs è disponibile in Pyomo (controlla gli import) se il pacchetto highspy è installato. Assicurati quindi di eseguire un pip install highspy.

solver = Highs()solver.highs_options = {    "log_file": "Highs.log",    "mip_heuristic_effort": 0.2,    "mip_detect_symmetry": True,    "mip_rel_gap": 1e-6,}solution = solve(model, solver)

Un’altra funzione per trovare i tour creati e siamo pronti per visualizzare i risultati.

def find_tours(model):    tours = []    for k in model.K:        node = 0        tours.append([0])        while True:            for j in model.V:                if (node, j) in model.A:                    if np.isclose(model.x[node, j, k].value, 1):                        node = j                        tours[-1].append(node)                        break            if node == 0:                break    return tours

Risultati sorprendenti per la piccola istanza con un totale di 10 nodi. Tuttavia, anche per questa piccola istanza, il solver ha impiegato quasi mezzo minuto per ottenere la soluzione e la complessità aumenta significativamente con più punti di domanda. Fortunatamente, esistono algoritmi specializzati disponibili pubblicamente per trovare soluzioni di alta qualità per istanze molto più grandi in breve tempo computazionale. Diamo un’occhiata a loro nella prossima sezione.

Euristiche specializzate (Meta)Heuristics

Nel corso degli anni sono state proposte diverse euristiche specializzate (meta) per varianti del VRP. La maggior parte di esse si basa fortemente su algoritmi di ricerca locale in modo che, dato una soluzione, vengano provate diverse perturbazioni per migliorarne sequenzialmente il costo fino a quando non sia più possibile migliorare ulteriormente nel vicinato dato. Quando si utilizzano Google Or Tools, utilizzeremo anche metodi di ricerca locale associati ad algoritmi costruttivi.

In questa sezione, verrà utilizzata l’istanza 150d di Rochat e Taillard (1995). I suoi dati sono stati ottenuti dal CVRPLIB. Questa istanza ha 150 clienti e un nodo di deposito, quindi sicuramente non saremmo in grado di risolverla utilizzando la strategia MIP precedentemente presentata.

Iniziamo ancora una volta importando le librerie utilizzate.

from itertools import cycleimport numpy as npimport pandas as pdfrom scipy.spatial.distance import pdist, squareformimport matplotlib.pyplot as pltimport matplotlib as mplfrom ortools.constraint_solver import routing_enums_pb2from ortools.constraint_solver import pywrapcp

E carichiamo i dati del problema da un file che include le coordinate e la domanda di ogni nodo.

dataset = pd.read_csv("./data/tai150d.csv", index_col=0)coordinates = dataset.loc[:, ["x", "y"]]demands = dataset.d.valuescapacity = 1874n_vehicles = 15N = coordinates.shape[0]distances = squareform(pdist(coordinates, metric="euclidean"))distances = np.round(distances, decimals=4)

Il primo passo per utilizzare il risolutore OR Tools VRP è istanziare un gestore di routing e un modello.

# Creare il gestore di routing per il numero di nodi, il numero di veicoli, il nodo di depositomanager = pywrapcp.RoutingIndexManager(    N, n_vehicles, 0)# Creare il modello di routingrouting = pywrapcp.RoutingModel(manager)

Successivamente, includeremo callback per quantificare le dimensioni legate ad archi/bordi e nodi. Il metodo RegisterTransitCallback della nostra istanza di routing può essere utilizzato per quantificare qualsiasi dimensione legata ad archi/bordi, mentre il metodo RegisterUnaryTransitCallback può quantificare valori legati ai nodi.

# Stesso valido per qualsiasi callback relativa ad archi/bordidef distance_callback(from_index, to_index):    from_node = manager.IndexToNode(from_index)    to_node = manager.IndexToNode(to_index)    return distances[from_node, to_node]transit_callback_index = routing.RegisterTransitCallback(distance_callback)# Stesso valido per qualsiasi callback relativa ai nodidef demand_callback(from_index):    from_node = manager.IndexToNode(from_index)    return demands[from_node]demand_callback_index = routing.RegisterUnaryTransitCallback(demand_callback)

Ora, includeremo un vincolo di capacità utilizzando l’indice demand_callback_index precedentemente definito. Si noti che i vincoli di durata avrebbero potuto essere definiti utilizzando la stessa sintassi, passando istanze di RegisterTransitCallback come primo argomento. Inoltre, è importante sottolineare che il modello di routing gestisce flotte eterogenee, quindi dobbiamo passare una lista di valori nel terzo argomento.

# Qualsiasi vincolo associato ai veicoli può avere gli stessi argomentirouting.AddDimensionWithVehicleCapacity(    demand_callback_index,    0,  # margine di capacità nullo    [capacity,] * n_vehicles,  # capacità massima dei veicoli (lista per ogni veicolo)    True,  # inizio cumulativo a zero    'Capacità')

Allo stesso modo, la definizione dell’obiettivo richiede anche una callback come argomento principale. In questo esempio, cerchiamo di minimizzare le distanze definite nell’indice transit_callback_index.

routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

Infine, dobbiamo definire i parametri del solutore e risolvere il nostro modello.

# Impostazione delle strategie euristichesearch_parameters = pywrapcp.DefaultRoutingSearchParameters()search_parameters.first_solution_strategy = (    routing_enums_pb2.FirstSolutionStrategy.CHRISTOFIDES)search_parameters.local_search_metaheuristic = (    routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)search_parameters.time_limit.FromSeconds(300)# Risolvere i problemisolution = routing.SolveWithParameters(search_parameters)

Il seguente pezzo di codice può essere utilizzato per estrarre i percorsi utilizzati nella nostra soluzione.

tours = []for vehicle_id in range(n_vehicles):    index = routing.Start(vehicle_id)    tours.append([])    while not routing.IsEnd(index):        node_index = manager.IndexToNode(index)        previous_index = index        index = solution.Value(routing.NextVar(index))        tours[-1].append(node_index)    else:        node_index = manager.IndexToNode(index)        tours[-1].append(node_index)

E si può accedere al valore obiettivo eseguendo la semplice linea solution.ObjectiveValue(). Utilizzando la configurazione presentata, ho ottenuto un obiettivo di 2679, che è abbastanza vicino al valore ottimale comprovato di 2645 (scarto del 1,2%). I percorsi ottenuti sono rappresentati di seguito.

Il codice completo (inclusi i grafici) è disponibile in questo repository git.

Estensioni utili

La libreria OR Tools è fantastica come risolutore generale per problemi di routing in quanto gestisce diverse varianti del VRP, come finestre temporali, flotte eterogenee e depositi multipli. Tuttavia, gli algoritmi adatti al CVRP canonico possono funzionare ancora meglio. Vale sicuramente la pena dare un’occhiata al pacchetto HGS-CVRP (Vidal, 2022) che combina un algoritmo genetico all’avanguardia con diverse mosse di ricerca locale. L’algoritmo trova la soluzione ottimale per l’istanza tai150d in meno di 20 secondi.

Riguardo ad alcune applicazioni reali, è probabile che sia necessario fare affidamento sulle distanze stradali (o durate) anziché sulle distanze euclidee per collegare le posizioni. Google fornisce un’interfaccia a pagamento molto utile per farlo, che puoi controllare in questo tutorial. Tuttavia, se stai cercando alternative open source, vale la pena dare un’occhiata all’API di OpenStreetMap. Alcune richieste utili sono:

• https://router.project-osrm.org/table/v1/driving/<LOCATIONS>
• http://router.project-osrm.org/route/v1/car/<LOCATIONS>?overview=false&steps=true

In cui <LOCATIONS> dovrebbe essere un elenco di coppie di longitudine e latitudine separate da virgole all’interno delle coppie e da punti e virgola tra diverse coppie. Puoi anche specificare sorgenti e destinazioni nella richiesta della tabella, il che è utile nel caso in cui la tabella completa sia troppo grande da gestire in una singola richiesta.

Oltre a effettuare calcoli di routing precisi, la visualizzazione può essere uno strumento importante. La libreria Python folium può essere molto utile per farlo. Vale sicuramente la pena dare un’occhiata.

Ulteriori letture

In precedenza in questo articolo abbiamo implementato un modello MIP esatto per il CVRP, che non è adatto per istanze di dimensioni moderate. Tuttavia, gli algoritmi che combinano la generazione di colonne con il Branch and Cut sono stati efficaci nella risoluzione di istanze con fino a qualche centinaio di clienti. Vale la pena dare un’occhiata ai documenti di ricerca di Fukasawa et al. (2006) e Pecin et al. (2017).

Coloro interessati a un’introduzione precedente alla generazione di colonne potrebbero trovarla nel mio precedente articolo VoAGI.

Riguardo alle meta-euristiche, i documenti di Vidal et al. (2012) e Vidal (2022) sono fantastici. Entrambi sono disponibili anche sotto forma di rapporti tecnici, con link disponibili nel repository HGS-CVRP.

Conclusioni

In questo articolo sono state presentate due approcci per risolvere il problema di routing dei veicoli capacitated (CVRP): la programmazione lineare intera e le (meta)euristiche. La prima alternativa è stata utilizzata per risolvere un’istanza di piccole dimensioni in cui è stata efficace, anche se non è in grado di gestire istanze di dimensioni moderate o grandi. Il secondo approccio è stato utilizzato per risolvere un problema impegnativo della letteratura con 150 clienti per il quale il risolutore ha trovato una soluzione di buona qualità con uno scarto del 1,2% rispetto all’ottimo conosciuto entro 300 secondi.

Riferimenti

Bynum, M. L. et al., 2021. Pyomo-optimization modeling in python. Springer.

Dantzig, G. B., & Ramser, J. H., 1959. The truck dispatching problem. Management science, 6(1), 80–91.

Fukasawa, R., Longo, H., Lysgaard, J., Aragão, M. P. D., Reis, M., Uchoa, E., & Werneck, R. F., 2006. Robust branch-and-cut-and-price for the capacitated vehicle routing problem. Mathematical programming, 106, 491–511.

Pecin, D., Pessoa, A., Poggi, M., & Uchoa, E., 2017. Miglioramento del branch-cut-and-price per il routing dei veicoli con capacità. Mathematical Programming Computation, 9, 61–100.

Rochat, Y., & Taillard, É. D., 1995. Diversificazione e intensificazione probabilistica nella ricerca locale per il routing dei veicoli. Journal of heuristics, 1, 147–167.

Toth, P., & Vigo, D., 2002. Panoramica dei problemi di routing dei veicoli. The vehicle routing problem, 1–26.

Vidal, T., 2022. Ricerca genetica ibrida per il CVRP: implementazione open-source e quartiere SWAP*. Computers & Operations Research, 140, 105643.

Vidal, T., Crainic, T. G., Gendreau, M., Lahrichi, N., & Rei, W., 2012. Un algoritmo genetico ibrido per problemi di routing multidepot e periodici dei veicoli. Operations Research, 60(3), 611–624.