Einführung in die Welt der Tourenoptimierung – Visualisierung und Lösungsverfahren in einem Python Jupyter Notebook (2/3)

In diesem Artikel möchte ich euch zeigen, wie ihr Probleme der Tourenoptimierung in einem Python Jupyter Notebook lösen und visualisieren könnt. Am Beispiel eines Fahrradkurierdienst zeige ich außerdem, wie das Grundproblem um gängige Nebenbedingungen erweitert werden kann. Dieser Artikel ist der zweite meiner Tourenoptimierungs-Blogreihe und baut auf dem ersten Blogbeitrag  auf, in dem wir die theoretischen Grundlagen besprochen haben. Mit ein paar Grundkenntnissen solltet ihr aber folgen können, auch ohne diesen Artikel gelesen zu haben.

Anforderungsanalyse für die Tourenoptimierung

Relevant für die Praxis ist in den meisten Fällen nicht die Basisvariante des Traveling Salesman Problem, sondern Abwandlungen von dieser. So gibt es eine ganze Reihe von Beschränkungen und Nebenbedingungen die zusätzlich berücksichtigt werden müssen. Bei vielen Transportproblemen stehen mehrere Fahrer oder Fahrzeuge zur Auswahl und es ist nicht von vornherein festgelegt, wer welchen Zielpunkt anfahren soll. Zusätzlich gibt es häufig Kapazitätsbeschränkungen. Ein Transporteur kann nur eine bestimmte Anzahl von Transportgütern transportieren. Typisch ist außerdem die Priorisierung bestimmter Kunden und die Einhaltung von Lieferzeitfenstern. Weiterhin müssen Unternehmen sicherstellen, dass sie innerhalb kürzester Zeit auf wechselnde Bedingungen reagieren können, z. B. neue kurzfristige Aufträge oder dynamische Verkehrsbedingungen. Je nach Anwender kann auch die Berücksichtigung von wechselnden Wetterbedingungen oder unterschiedlichen Ausliefergeschwindigkeiten von verschiedenen Fahrzeugen oder Fahrern wichtig sein. Diese Liste kann noch beliebig erweitert werden. Die meisten der zuvor aufgezählten Punkte sind zugleich super praxisrelevant und gut erforscht. Daher gibt es viele Bibliotheken, in denen diese Fälle bereits implementiert wurden. Zur reinen mathematischen Optimierung bieten sich verschiedene Open-Source-Lösungen an. Eine dieser Bibliotheken sind die Google-OR-Tools . Im nächsten Abschnitt zeige ich euch, wie ihr mit Hilfe dieser Bibliothek einen einfachen Demo Case (Hier gehts zum Repository ) implementieren könnt und wie dieser in einem Jupyter Notebook mit Hilfe der Folium-Bibliothek visualisieren könnt. 

Zu Beginn der Implementierung erläutere ich, wie der Demo Use Case aussieht. Dabei beginne ich mit einem Grundproblem, löse dieses und erweitere diese einfache Basisanwendung dann im nächsten Schritt um realistischere Anforderungen und Nebenbedingungen.

Grundproblem: Ein Fahrer

Für den Demo Case gehe ich von einem Fahrradkurier aus, der seine Kunden ausgehend von einem Zentrallager beliefern soll. Die Zuteilung der Zielpunkte auf die verschiedenen Fahrer ist bereits erfolgt und Kapazitätseinschränkungen und Lieferzeitfenster müssen zunächst nicht berücksichtigt werden. Der Fokus liegt ausschließlich darauf, 10 Zielpunkte in einer möglichst kurzen Route effizient miteinander zu verbinden und das Ergebnis zu visualisieren. Dafür habe ich 10 zufällige Adressen in Münster gewählt. Das grüne Icon symbolisiert dabei den Start- und Endpunkt. Die 10 blauen Icons symbolisieren die Zielpunkte.

Zur Umsetzung müssen ein paar Bibliotheken importiert und installiert werden. Neben den bekannten Bibliotheken Pandas und NumPy zum Data Handling, benötigen wir Folium zum Erstellen von Karten (ähnlich der oben abgebildeten), Geopy zum Bestimmen der Geokoordinaten unserer Zieladressen, Ortools zur Tourenoptimierung und Pyroutelib3 zur Routenoptimierung.

1import pandas as pd
2import numpy as np
3import folium
4 
5from geopy import Nominatim
6from geopy.distance import great_circle as GC
7 
8from ortools.constraint_solver import pywrapcp, routing_enums_pb2
9from pyroutelib3 import Router

Als nächstes erstelle ich eine Liste von zufällig ausgewählten Adressen in Münster.

1# Define list of target points addresses:
2addresses = [
3    'Sternstraße 12 Münster',
4    'Bohlweg 17 Münster',
5    'Scharnhorststraße 85 Münster',
6    'Geiststraße 46 Münster',
7    'Kanalstraße 14 Münster',
8 
9    'Junkerstraße 2 Münster',
10    'Alter Fischmarkt 11 Münster',
11    'Scharnhorststraße 2 Münster',
12    'Königsstraße 3 Münster',
13    'Fürstenbergstraße 30 Münster',
14]

Um die Geokoordinaten der Adressen abzurufen, nutze ich die geopy-Library. Nominatim ist ein Geocoding-Tool, das basierend auf OpenStreetMaps-Daten die Geokoordinaten von Adressen ermittelt. Dazu iteriere ich über jede Adresse in der Liste und erstellen eine neue Liste, die zu jedem Zielpunkt die Adresse und die Geokoordinaten beinhaltet.

1# Get coordinates and create list of target points including address and coordinates retrieved from Open Street Map Server:
2geolocator = Nominatim(user_agent="TSP_Muenster")
3 
4addresses_geo_data = []
5for address in addresses:
6 
7    loc = geolocator.geocode(address)
8 
9    # Handle cases in which address couldn't be retrieved:
10    if loc is None:
11        print(f"Address couldn't be located: {address}")
12        break
13 
14    addresses_geo_data.append({
15        'address': address,
16        'coords': (loc.latitude,loc.longitude)
17    })

Außerdem erstelle ich noch ein Ausgangslager, das am Münsteraner Hafen liegen soll.

1# Create Warehouse:
2warehouse_address = 'Am Mittelhafen 14 Münster'
3loc = geolocator.geocode(warehouse_address)
4 
5warehouse = {
6    'address': warehouse_address,
7    'coords': (loc.latitude,loc.longitude)
8}

Im nächsten Schritt können wir die bereits zu Beginn gezeigte Karte visualisieren. Dazu verwende ich die Bibliothek Folium, die als Python Wrapper um die beliebte Leaflet Javascript Library zum Erstellen interaktiver Karten verwendet wird.
Im ersten Schritt erstelle ich zunächst eine Karte zentriert auf die Koordinaten des Zentrallagers. Das Zoom-Level stelle ich so hoch, dass es einen schönen Überblick über Münster liefert. Als nächstes erstelle ich ein Icon für das Zentrallager. Neben den Koordinaten, dem Inhalt beim Hovern (Tooltip) und dem Inhalt des Popup-Fensters, lässt sich auch die Farbe für das Icon anpassen. Neben einigen built-in Icons können auch Font-Awesome v4 Icons verwendet werden. Ein solches nutze ich für das Zentrallager. Abschließend iteriere ich über die Zielpunkte, erstelle auch hier jeweils ein Icon und gebe dann die resultierende Karte aus.

1# Create map object:
2map_osm = folium.Map(location=warehouse['coords'], zoom_start=14, tiles='Open Street Map')
3 
4# Create icon for warehouse:
5folium.Marker(
6    location=warehouse['coords'],
7    icon=folium.Icon(color='green', icon='industry', prefix='fa'),
8    popup=warehouse['address'],
9    tooltip=warehouse["address"],
10    draggable=False).add_to(map_osm)
11 
12# Create icons for each address:
13for address in addresses_geo_data:
14    folium.Marker(
15        location=address['coords'],
16        icon=folium.Icon(color='darkblue', icon='home'),
17        popup=address['address'],
18        tooltip=address["address"],
19        draggable=False).add_to(map_osm)
20 
21map_osm

Im nächsten Schritt stelle ich die Distanzmatrix zwischen den Zielpunkten und dem Ausgangspunkt auf. Dazu werden zwei Hilfsfunktionen benötigt. Die erste dient dem Zweck, die Distanz zwischen zwei Punkten zu berechnen. Der Einfachheit halber verzichte ich für den Moment darauf, die genaue Fahrdistanz zu verwenden, sondern nutze die Luftlinie zwischen zwei Punkten als Proxy für die tatsächliche Entfernung. Zum Erstellen der Distanzmatrix verwende ich eine Funktion, die als flexibles Rückgabeformat entweder ein Pandas DataFrame oder ein NumPy Array zurückgibt. Während auf das NumPy Array mit Hilfe der Zielpunkt-Indizes einfach zugegriffen werden kann, kann mit dem Pandas DataFrame eine übersichtliche Tabelle erstellt werden, um vor allem beim Entwickeln schnell überprüfen zu können, ob alles wie erwartet funktioniert.

1def get_distance(point1: tuple, point2: tuple) -> float:
2 
3    # Air distance in meters:
4    dist = GC(point1, point2).km * 1000
5 
6    return round(dist,2)
7 
8 
9def create_distance_matrix(coords: list, address_name:list, format='DataFrame', verbose=False):
10    # Create empty np-array:
11    dist_array = np.empty((len(coords),len(coords)))
12 
13    # Compute distances:
14    for i in range(0,len(coords)):
15        for j in range(i,len(coords)):
16            if i < j:
17                dist = get_distance(coords[i],coords[j])
18 
19                if verbose:
20                    print(f"Distance between: {address_name[i]} and {address_name[j]}: {dist} km")
21 
22                # Assuming symmetric TSP:
23                dist_array[i][j] = dist
24                dist_array[j][i] = dist
25 
26            elif i == j:
27                dist_array[i][i] = 0
28 
29            else:
30                continue
31 
32    if format == 'NumpyArray':
33        return dist_array
34 
35    else:
36        # Create pandas distance matrix:
37        return pd.DataFrame(data=dist_array, index=address_name, columns=address_name)

Nun verbinde ich die Liste aus Zieladressen und Ausgangspunkt und nutze die zuvor erstellten Funktionen, um die Distanzmatrix zu erstellen.

1# Create list of addresses including warehouse and target points to be used when creating the distance matrix dataframe as column names and index names:
2address_list = [warehouse['address']] + [obj['address'] for obj in addresses_geo_data]
3coords_list = [warehouse['coords']] + [obj['coords'] for obj in addresses_geo_data]
4 
5# Call function twice for returning a numpyArray and a DataFrame:
6dist_matrix_np = create_distance_matrix(address_name=address_list, coords=coords_list, verbose=False, format='NumpyArray')
7dist_matrix_df = create_distance_matrix(address_name=address_list, coords=coords_list, verbose=False, format='DataFrame')

Optimierung mit Google OR-Tools

Jetzt kommt der interessante Teil: Das Finden einer effizienten Tour. Natürlich können ich dafür eine der Funktionen verwenden, die ich im ersten Blogartikel eingeführt hatte. Ich möchte diese Demo aber für einen kleinen Einblick in die Google OR-Tools Library für Tourenoptimierung nutzen. Mit dieser mächtigen Bibliothek können auch die gängigsten Beschränkungen und Nebenbedingungen ohne großen Implementierungsaufwand abgebildet werden.
Zunächst erstelle ich einen IndexManager und ein RoutingModel-Objekt. Der IndexManager dient dazu, die internen Indizes des Solvers in die Indizes der Distanzmatrix zu konvertieren. Dazu übergebe ich die Anzahl der Knoten, also alle 11 Punkte, die Anzahl der Fahrzeuge und den Index unseres Ausgangspunkts. Das RoutingModel-Objekt ist das zentrale Objekt, um das sich ab jetzt so ziemlich alles dreht. Sobald alles konfiguriert ist, kann mit Hilfe der Solve-Methode eine Lösung berechnet werden.

1# The index manager handles the conversion of the solver's internal indices to the location indices of our distance matrix:
2index_manager = pywrapcp.RoutingIndexManager(len(dist_matrix_np), 1, 0) ## 11 nodes, 1 vehicle, 0 warehouse_index
3 
4# The routing model is the central object that we can configure to solve our problem:
5routing_model = pywrapcp.RoutingModel(index_manager)

Vorher muss ich aber noch weitere Konfigurationen treffen. Als nächstes wird eine Callback-Funktion benötigt, welche die Distanz zwischen zwei Punkten zurückgibt. Hier werden die Distanzmatrix und das RoutingModel miteinander verknüpft.
Weiterhin müssen noch die Kosten für das Wählen einer Verknüpfung zwischen zwei Punkten festlegen. Moment mal, habe ich das nicht gerade gemacht? Nicht ganz. In diesem Fall will ich als Kosten die Distanz zwischen zwei Punkten verwenden. Theoretisch könnte man aber hier auch noch weitere oder alternative Faktoren einfließen lassen. Beispielsweise könnten wir statt der Distanz die Fahrzeit und somit auch mögliche Geschwindigkeitsunterschiede bei verschiedenen Fahrern oder Fahrzeugen berücksichtigen.
Als nächstes muss ich eine Optimierungsmethode wählen. Ich wähle dafür die Default-Strategie. Um eine initiale Lösung zu erhalten (First Solution Strategy), setze ich auf den Greedy-Ansatz, ähnlich wie im ersten Teil der Artikelreihe beschrieben.
Abschließend ermittelt die Methode SolveWithParameters die Lösung.

1# Function that returns the distance between two points. In this case the easiest solution is to use the distance matrix:
2def distance_callback(from_index, to_index):
3    # Convert from routing variable Index to distance matrix NodeIndex.
4    from_node = index_manager.IndexToNode(from_index)
5    to_node = index_manager.IndexToNode(to_index)
6    return dist_matrix_np[from_node][to_node]
7 
8# Create Callback that is needed to connect our routing model object with the distance matrix:
9transit_callback_index = routing_model.RegisterTransitCallback(distance_callback)
10 
11# The Arc Cost Evaluator tells the model how to compute the costs of choosing one arc and thus driving from one point A to another point B:
12routing_model.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)
13 
14# Set the search strategy to the default strategy:
15search_parameters = pywrapcp.DefaultRoutingSearchParameters()
16 
17# Use the 'greedy approach' to create an initial solution that is then improved:
18search_parameters.first_solution_strategy = (
19    routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)

Zur Ausgabe des berechneten Tourenplans benötige ich eine weitere Hilfsfunktion. Hier iteriere ich über die zuvor berechnete Lösung und erstelle eine Liste, die die geplante Reihenfolge enthält. Da ich die Funktion später auch für Touren mit mehreren Fahrzeugen nutze, berücksichtige ich das gleich mit und erstelle für jedes Fahrzeug eine weitere eigene Liste, die die geplante Tour enthält.

1def get_tours(solution, routing, manager):
2    # Get vehicle routes and store them in a two dimensional array whose
3    # i,j entry is the jth location visited by vehicle i along its route.
4    all_tours = []
5    for vehicle in range(routing.vehicles()):
6        index = routing.Start(vehicle)
7        curr_tour = [manager.IndexToNode(index)]
8        while not routing.IsEnd(index):
9            index = solution.Value(routing.NextVar(index))
10            curr_tour.append(manager.IndexToNode(index))
11        all_tours.append(curr_tour)
12    return all_tours
13 
14tour_plan = get_tours(solution, routing_model, index_manager)

Jetzt fehlt noch die Visualisierung. Dazu erweitere ich den bestehenden Code um eine Schleife, die die einzelnen Zielpunkte der Tour miteinander verbindet. Auch hier abstrahiere ich von den tatsächlichen Routen zwischen zwei Zielpunkten und zeichne stattdessen eine gerade Linie zwischen zwei Punkten.

1## Create map object:
2map_osm = folium.Map(location=warehouse['coords'], zoom_start=14, tiles='Open Street Map')
3 
4## Create icon for warehouse:
5folium.Marker(
6    location=warehouse['coords'],
7    icon=folium.Icon(color='green', icon='industry', prefix='fa'),
8    popup=warehouse['address'],
9    tooltip=warehouse["address"],
10    draggable=False).add_to(map_osm)
11 
12## Create icons for each address:
13for address in addresses_geo_data:
14    folium.Marker(
15        location=address['coords'],
16        icon=folium.Icon(color='darkblue', icon='home'),
17        popup=address['address'],
18        tooltip=address["address"],
19        draggable=False).add_to(map_osm)
20 
21## Create connections between target points:
22# Iterate over tour plan:
23for tour in tour_plan:
24    # Iterate from first to penultimate element:
25    for i in range(0,len(tour) - 1):
26        coords_point_a = coords_list[tour[i]]
27        coords_point_b = coords_list[tour[i+1]]
28 
29        folium.PolyLine([coords_point_a, coords_point_b],
30                        color="black",
31                        weight=3).add_to(map_osm)
32 
33map_osm

Das sieht doch schon ganz schick aus, oder?

Im nächsten Artikel der Reihe erweitere ich den Demo Case um die exakte Routen zwischen zwei Zielpunkten, statt weiter die Luftlinien als Annäherung für die tatsächliche Entfernung zu nutzen. Außerdem zeige ich, wie wir das Grundproblem um realistische Anforderungen wie mehrere Fahrzeuge und Kapazitätsbeschränkungen erweitern können.