Vehicle Routing Problem mit Quantum-Annealing

Das Vehicle Routing Problem (VRP) ist ein NP‑schweres Optimierungsproblem in der Logistik. In unserer modernen Welt ist es essenziell, eine ausreichend gute Lösung innerhalb kurzer Zeit zu finden. In den letzten Jahren hat die Methode des Quantencomputings zum Lösen des VRPs mehr Aufmerksamkeit bekommen. Wie in einem Reviewpapier von Liu et al.¹ beschrieben, wird hauptsächlich Quantum-Annealing für Optimierungsprobleme in der Logistik genutzt. Wir haben den Quantum-Annealing-Ansatz in einem vorherigen Blogpost behandelt und mit dem gatter-basierten Framework für Quantencomputing verglichen.
In diesem Blogpost wollen wir zeigen, wie man Quantum Annealing zum Lösen eines VRPs nutzen kann. Wir folgen dem von Feld et al.² vorgeschlagenen Ansatz. In einem VRP hat man einen Fuhrpark, eine Liste von Kunden und ein Depot mit deren Positionen. Die Fahrzeuge haben eine vordefinierte Kapazität und die Kunden eine zu erfüllende Nachfrage nach Produkten. Man muss entscheiden, welches Fahrzeug die Produkte zu welchen Kunden liefert. Gleichzeitig muss die Lieferzeit minimiert werden, während die Fahrzeugkapazität nicht überschritten werden darf. Feld et al. teilen das VRP in zwei Phasen: einer Clusterphase und einer Routephase.
Dieser Post ist Teil des QUICS-Projektes, einer Kollaboration von GWDG, LUH und PTB. QUICS bietet KMU einen niedrigschwelligen Zugang zu Quantencomputing, um Anwendungsfälle mit verschiedener Quantenhardware und -software zu testen. Unserer Code ist verfügbar in unseren D-Wave Apptainer Container. Er ist in dem Ordner /opt/notebooks/VRP/ zu finden. Eine Liste unserer derzeitigen Simulatoren ist hier zu finden. Später wird es die Möglichkeit geben Quantensimulatoren in der HPC-Umgebung des QUICS-Projektes zu nutzen.

Die Clusterphase

Die Clusterphase kann als ein Rucksackproblem betrachtet werden, bei dem man unterschiedlich große Objekte in Rucksäcke mit einer bestimmten Kapazität packt, welche hier die Fahrzeuge sind. Daher wird für jedes Cluster ein Fahrzeug genutzt. Feld et al. lösen das Problem klassisch.
Der erste Teil dieser Phase ist die Clustergenerierung, die wie folgt strukturiert ist:

1. Wähle den Kern des Clusters. Feld et al. schlagen zwei unterschiedliche Methoden vor:
   a) wähl den Kunden mit der größten Nachfrage, oder
   b) wähl den Kunden welcher am weitesten entfernt vom Depot ist.
2. Berechne das geometrische Zentrum des Clusters:
   $$CC(m_k) = \sum\limits_{i=0}^n \frac{v_i^x}{n}, \sum\limits_{i=0}^n \frac{v_i^y}{n}\, ,$$ wobei $m_k$ den $k$-ten Cluster beschreibt.
3. Wähle den Kunden mit der kürzesten Reichweite zum momentanen geometrischen Zentrum des Clusters und füge ihn zum Cluster hinzu falls dessen Nachfrage nicht die restliche Kapazität des Fahrzeuges überschreitet.
4. Berechne $CC(m_k)$ neu und wiederhole (3.).
5. Falls das Fahrzeug voll ist, beende die Erstellung des Clusters und wiederhole von (1.) für das nächste Cluster.

Als Zweites müssen die Cluster verbessert werden. Dafür wählen wir zufällig einen Kunden und berechnen dessen Abstand zum geometrischen Zentrum des ursprünglichen Clusters. Dann überprüfen wir, in welche der anderen Cluster die Nachfrage des Kunden passt und ob der Kunde näher an einem dieser zulässigen Cluster ist. Falls dies der Fall ist, wird der Kunde dem nächstgelegenen Cluster neu zugeordnet. Dies wird mehrmals wiederholt.

Die Routephase

Nach der Clusterphase haben wir für jedes Cluster ein Problem des Handlungsreisenden (“Travelling Salesman Problem”, TSP). In einem TSP gibt es nur ein Fahrzeug (Reisender) und das Ziel ist es, die optimale/kürzeste Route zu finden. Diese Optimierungsprobleme können mithilfe von Quantum‑Annealing‑Geräten gelöst werden. Um das Problem für Quantum Annealing zu kodieren, betrachten wir jedes TSP als Graphen, wobei die Kunden und das Depot die Knoten sind und die Verbindungen zwischen allen Kunden und dem Depot die Kanten sind. Die Abbildung unten veranschaulicht ein TSP zusammen mit zwei Beispiellösungen. Beide Touren sind realisierbar, da jeder Kunde genau einmal aufgesucht wird. Jedoch haben die beiden Touren nicht die gleiche Güte. Die erste Tour hat eine Länge von $17,7$ und die zweite (gute) Tour hat eine Länge von $15,8$. Diese Tourlänge³ ist das Minimierungsziel unseres Beispiel-TSPs.

Zuerst muss der Hamiltonian des Problems durch die Kombination der Nebenbedingungen und der Zielfunktion definiert werden. Dafür definieren wir die binären Variablen $x_{uj}$, welche $1$ sind, falls der Kunde $u$ an Position $j$ der Tour ist. Die Fahrzeiten oder Entfernungen zwischen allen Kunden können in einer Matrix $D$ aufgelistet werden. Dann kann die zu minimierende Zielfunktion geschrieben werden als

wobei $E$ das Set aller Kanten/Verbindungen zwischen den Kunden ist und $n$ ist die Anzahl an Kunden. Die erste Summe läuft über all möglichen Verbindungen des Problemgraphen. In der zweiten Summe sind der nullte Halt und der $n+1$-te Halt das Depot. Dies stellt sicher, dass kürzere Touren bevorzugt werden, indem der Term $l\cdot B$ proportional zur Gesamtlänge der Tour $l$ anwächst.
Es gibt drei Nebenbedingungen in einem TSP: Die Kundenhäufigkeitsbedingung, die Tourpositionsbedingung, und die Reihenfolge-der-Kunden-Bedingung. Die Kundenhäufigkeitsbedingung stellt sicher, dass jeder Kunde genau einmal in der Tour vorkommt:

Die Tourpositionsbedingung stellt sicher, dass jede Position in der Tour genau einmal einem Kunden zugeordnet wird:

Die Reihenfolge-der-Kunden-Bedingung stellt sicher, dass die Reihenfolge der Kunden möglich ist:

Diese Bedingung muss nur genutzt werden, wenn nicht jeder Kunde mit jedem anderen Kunden verbunden ist. Der komplette Hamiltonian des Problems lautet $H = H_O + H_{C_1} + H_{C_2} + H_{C_3}$. Die Parameter $B$ und $A$ müssen sorgfältig gewählt werden, sodass die Nebenbedingungen eingehalten werden, z.B. $0< B\cdot \max(D_{uv}) <A$.
Die Beschreibung mit einem Hamiltonian ist äquivalent zu einem QUBO‑Problem. In unserem Code benutzen wir die QUBO‑Formulierung als Input für das Quantum‑Annealing‑Gerät. Nach dem Anwenden der binomischen Formeln und dem Ausmultiplizieren der Summen erhält man die folgende Implementierung:

def objective_function(Q, B, x_var_index, travel_time_dict, customers, num_customers):
    """
    Q: QUBO matrix
    B: Penalty of the objective function
    x_var_index: Dictionary of the indices of the binary x variable
    travel_time_dict: contains all the travel times of all edges. Should be a dictionary with keys "uv"
    customers: The customer names / numbers of the TSP
    num_customers: Number of customers (not including the depot)
    """
    
    # Depot -> first customer:
    for u in customers:
        travel_time = travel_time_dict[("0", str(u))]
        idx = x_var_index[(u, 1)]
        Q[idx, idx] += B * travel_time

    # Customer -> customer transitions:
    for u in customers:
        for v in customers:
            if u != v:
                travel_time = travel_time_dict[(u, v)]
                for j in range(1, num_customers):
                    idx_1 = x_var_index[(u, j)]
                    idx_2 = x_var_index[(v, j + 1)]
                    Q[idx_1, idx_2] += B * travel_time

    # Last customer -> depot:
    for u in customers:
        travel_time = travel_time_dict[(str(u), "0")]
        idx = x_var_index[(u, num_customers)]
        Q[idx, idx] += B * travel_time
    return Q

def customer_occurance_constraint(Q, A, x_var_index, customers, num_customers):
    """
    Q: QUBO matrix
    A: Penalty of the constraint
    x_var_index: Dictionary of the indices of the binary x variable
    customers: The customer names / numbers of the TSP
    num_customers: Number of customers (not including the depot)
    """
    
    for u in customers:
        for j in range (1, num_customers+1):
            idx_1 = x_var_index[(u, j)]
            Q[idx_1, idx_1] -= A 

            for j_prime in range(1, num_customers+1):
                if j < j_prime:
                    idx_2 = x_var_index[(u, j_prime)]
                    Q[idx_1, idx_2] += 2*A 
    return Q

def tour_position_constraint(Q, A, x_var_index, customers, num_customers):
    """
    Q: QUBO matrix
    A: Penalty of the constraint
    x_var_index: Dictionary of the indices of the binary x variable
    customers: The customer names / numbers of the TSP
    num_customers: Number of customers (not including the depot)
    """
    
    for j in range(1, num_customers+1):
        for u in customers:
            idx_1 = x_var_index[(u, j)]
            Q[idx_1, idx_1] -= A 

            for u_prime in customers:
                if u < u_prime:
                    idx_2 = x_var_index[(u_prime, j)]
                    Q[idx_1, idx_2] += 2*A 
    return Q

Umsetzung mit Quantum-Annealing

Um diesen Ansatz zu testen, nutzen wir das CMT1-Datensatz von Christofides et al. als Beispiel. Zuerst kreieren wir die CLusters:

from create_clusters import cluster_generation, cluster_improvement

num_cluster_improve_iterations = 10000

initial_clusters = cluster_generation(customer_positions, customer_demand, depot_position, vehicle_capacity, use_demand_for_core_customer=False)
final_clusters = cluster_improvement(initial_clusters, num_cluster_improve_iterations, customer_positions, customer_demand, vehicle_capacity, depot_position)

Diese Python‑Funktionen können in unserem Code eingesehen werden. Dann wird für jedes Cluster ein TSP gelöst. Wir erzeugen die QUBO‑Matrix wie oben beschrieben und kreieren dann ein “binary quadratic model” mit der dimod-Bibliothek als Input für den Simulator. Danach benutzen wir simuliertes Annealing mit der Funktion neal.sampler.SimulatedAnnealingSampler, um das Problem zu lösen.

import dimod
from neal.sampler import SimulatedAnnealingSampler

B = 1
for cluster_idx in list(final_clusters.keys()):
    customers = final_clusters[cluster_idx]
    num_customers = len(customers)
    
    ...
    
    A = 1000 * num_customers * B * max(travel_time_dict.values())
    
    # Get the QUBO matrix:
    Q = construct_qubo_matrix_for_TSP(A, B, customers, num_customers, travel_time_dict, set_of_edges)
    qubo = qubo_from_numpy(Q)

    # Create the BQM:
    bqm = dimod.BQM.from_qubo(qubo)
    
    # Use simulated annealing:
    sampler = SimulatedAnnealingSampler()
    response = sampler.sample(bqm, num_reads=3000, num_sweeps=6000)
    
    best_sample = response.first.sample 

Wenn wir diesen Code ausführen, erhalten wir eine Gesamtlänge der Tour zwischen $1200$ und $1400$. Die beste bekannte Lösung für dieses Beispiel ist $524,61$. Wie man sehen kann, kommen wir mit unserem Experiment nicht in die Nähe der besten Lösung. Das Problem unserer Lösung ist, dass die Lösungen nicht unbedingt zulässig sind. Manche Kunden kommen mehrmals in der Tour vor. Um hochqualitative Lösungen für dieses mittelgroße Problem zu erhalten, muss man die Anzahl an “reads” (num_reads) und die Anzahl an “sweeps” (num_sweeps) des Simulators erhöhen. Jedoch erhöht dies die Berechnungszeit signifikant. Für unsere Wahl der Parameter braucht die Berechnung ungefähr $315,\text{s}$. Wenn wir unser Zeitergebnis mit dem von Feld et al. vergleichen, sehen wir, dass unsere Berechnung deutlich länger dauert. Dessen Berechnung auf dem D-Wave 2000Q Modell brauchte mit demselben Datensatz nur $15.792,\text{s}$. Folglich ist die Berechnung auf einem echten Quantum‑Annealer deutlich schneller. Des Weiteren ist deren Quantum‑Annealing‑Ergebnis sehr nahe an der besten bekannten Lösung.
In Feld et al. vergleichen sie ihre Berechnung mit einer klassischen Methode, die $1,474,\text{s}$ benötigte. Der Quantum‑Annealer braucht deutlich länger als die klassische Lösung, was auf den erhöhten Rechenaufwand zurückzuführen ist. Dieser Aufwand ist eine Folge der Einbettung des Problems auf den Hardware‑Graphen. Daher ist dieses Einbetten ein wesentlicher Teil der Quantum‑Annealing‑Forschung. Betrachtet man nur die tatsächliche Berechnungszeit auf dem QPU, braucht die Berechnung von Feld et al. nur $0,110,\text{s}$.
Zusammenfassend lässt sich sagen, dass echte Quantum‑Annealer zu ausreichend guten Lösungen führen können, auch wenn deren Rechenleistung noch nicht so gut ist wie die klassischer Methoden. Jedoch haben Quantum‑Annealer ein großes Potenzial, falls der Rechenaufwand reduziert werden kann. Außerdem führen Annealing‑Simulatoren in angemessener Zeit nicht zu zulässigen Lösungen. Solche Simulatoren können aber nützlich zum Testen von Algorithmen und Code sein, bevor ein echter Quantum‑Annealer genutzt wird.