{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Kapitel 8: Graphen, Flüsse und Touren — Min-Cost-Flow, Matching und VRP\n", "\n", "Begleitnotebook zu *Optimierte Entscheidungsfindung mit Python*. Die Codezellen sind identisch mit den im Buch abgedruckten Programmen.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Einmalig ausfuehren: installiert alle im Buch verwendeten Pakete.\n", "# Lokal in einer virtuellen Umgebung genauso gueltig wie in Google Colab.\n", "%pip install --quiet ortools highspy cvxpy scipy numpy pandas polars \\\n", " scikit-learn matplotlib plotly pyomo linopy pymoo pydantic openpyxl" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Das Minimum-Cost-Flow-Problem (MCNFP)\n", "\n", "`Min_Cost_Flow.py`\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#!/usr/bin/env python3\n", "\n", "# Min_Cost_Flow.py\n", "\"\"\"\n", "Kapitel Graphen: Kostenminimaler Fluss durch ein Netzwerk.\n", "\n", "Loest dasselbe Problem zweimal:\n", " (1) als allgemeines LP mit scipy -> zeigt die Modellstruktur\n", " (2) mit dem spezialisierten Netzwerk-Solver von OR-Tools -> zeigt den\n", " Geschwindigkeitsvorteil eines Verfahrens, das die Struktur ausnutzt\n", "\n", "Beide laufen in eigenen Prozessen nicht noetig: scipy und ortools vertragen\n", "sich (nur ortools + highspy kollidieren, siehe Kapitel Oekosystem).\n", "\"\"\"\n", "\n", "import numpy as np\n", "from scipy.optimize import linprog\n", "\n", "# --- Netzwerk definieren ---------------------------------------------------\n", "KNOTEN = [\"Werk_A\", \"Werk_B\", \"Umschlag\", \"Kunde_1\", \"Kunde_2\"]\n", "# (von, nach, Kosten je Einheit, Kapazitaet)\n", "KANTEN = [\n", " (\"Werk_A\", \"Umschlag\", 2.0, 15),\n", " (\"Werk_A\", \"Kunde_1\", 5.0, 10),\n", " (\"Werk_B\", \"Umschlag\", 4.0, 10),\n", " (\"Werk_B\", \"Kunde_2\", 6.0, 10),\n", " (\"Umschlag\", \"Kunde_1\", 1.0, 20),\n", " (\"Umschlag\", \"Kunde_2\", 3.0, 10),\n", "]\n", "# Angebot (+) bzw. Bedarf (-) je Knoten\n", "SALDO = {\"Werk_A\": 20, \"Werk_B\": 10, \"Umschlag\": 0, \"Kunde_1\": -15, \"Kunde_2\": -15}\n", "\n", "\n", "def loese_als_lp():\n", " \"\"\"Flussproblem als allgemeines lineares Programm.\"\"\"\n", " n_kanten = len(KANTEN)\n", " knoten_index = {k: i for i, k in enumerate(KNOTEN)}\n", "\n", " # Zielfunktion: Summe der Transportkosten\n", " kosten = np.array([k[2] for k in KANTEN])\n", "\n", " # Flusserhaltung als Gleichungssystem: A_eq @ x = b_eq\n", " A_eq = np.zeros((len(KNOTEN), n_kanten))\n", " for e, (von, nach, _, _) in enumerate(KANTEN):\n", " A_eq[knoten_index[von], e] = +1.0 # fliesst hinaus\n", " A_eq[knoten_index[nach], e] = -1.0 # fliesst hinein\n", " b_eq = np.array([SALDO[k] for k in KNOTEN], dtype=float)\n", "\n", " schranken = [(0, k[3]) for k in KANTEN] # 0 <= x_ij <= u_ij\n", "\n", " ergebnis = linprog(c=kosten, A_eq=A_eq, b_eq=b_eq, bounds=schranken, method=\"highs\")\n", " if not ergebnis.success:\n", " raise SystemExit(f\"Nicht loesbar: {ergebnis.message}\")\n", " return ergebnis.fun, ergebnis.x, ergebnis.eqlin.marginals\n", "\n", "\n", "if __name__ == \"__main__\":\n", " # Vorabpruefung: Angebot muss Bedarf entsprechen\n", " gesamt = sum(SALDO.values())\n", " print(\"=\" * 78)\n", " print(\" KOSTENMINIMALER FLUSS DURCH EIN TRANSPORTNETZ\")\n", " print(\"=\" * 78)\n", " print(f\"Angebot gesamt: {sum(v for v in SALDO.values() if v > 0)} | \"\n", " f\"Bedarf gesamt: {-sum(v for v in SALDO.values() if v < 0)} | \"\n", " f\"Saldo: {gesamt}\")\n", " if gesamt != 0:\n", " raise SystemExit(\"Angebot und Bedarf stimmen nicht ueberein - unloesbar!\")\n", "\n", " kosten_gesamt, fluss, knotenpreise = loese_als_lp()\n", "\n", " print(f\"\\nMinimale Transportkosten: {kosten_gesamt:,.2f} EUR\\n\")\n", " print(f\"{'Kante':<24} {'Fluss':>7} {'Kapazitaet':>11} {'Kosten/E':>9} {'Kosten':>9}\")\n", " print(\"-\" * 78)\n", " for e, (von, nach, c, u) in enumerate(KANTEN):\n", " menge = fluss[e] + 0.0 if abs(fluss[e]) > 1e-9 else 0.0 # vermeidet \"-0.0\"\n", " ausgelastet = \" (VOLL)\" if abs(menge - u) < 1e-6 else \"\"\n", " print(f\"{von + ' -> ' + nach:<24} {menge:>7.1f} {u:>11} \"\n", " f\"{c:>9.2f} {menge * c:>9.2f}{ausgelastet}\")\n", "\n", " # --- Flusserhaltung nachpruefen --------------------------------------\n", " print(\"\\n--- Pruefung der Flusserhaltung je Knoten ---\")\n", " for k in KNOTEN:\n", " hinaus = sum(fluss[e] for e, (v, n, _, _) in enumerate(KANTEN) if v == k)\n", " hinein = sum(fluss[e] for e, (v, n, _, _) in enumerate(KANTEN) if n == k)\n", " netto = hinaus - hinein\n", " art = \"Quelle\" if SALDO[k] > 0 else (\"Senke\" if SALDO[k] < 0 else \"Umschlag\")\n", " print(f\" {k:<10} ({art:<8}): hinaus {hinaus:5.1f} - hinein {hinein:5.1f} \"\n", " f\"= {netto:+6.1f} (gefordert: {SALDO[k]:+d})\")\n", " assert abs(netto - SALDO[k]) < 1e-6, f\"Flusserhaltung verletzt bei {k}!\"\n", "\n", " # --- Knotenpreise (Dualwerte) interpretieren -------------------------\n", " print(\"\\n--- Knotenpreise (Dualwerte der Flusserhaltung) ---\")\n", " print(\" Differenz zweier Knotenpreise = Grenzkosten einer zusaetzlichen Einheit\")\n", " print(\" auf dem guenstigsten Weg zwischen ihnen.\")\n", " for k, preis in zip(KNOTEN, knotenpreise):\n", " print(f\" {k:<10}: {preis:7.2f}\")\n", " print(\"=\" * 78)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Der Satz von Birkhoff und von Neumann{idx:Satz von Birkhoff und von Neumann}\n", "\n", "`Zuordnung_Ungarisch.py`\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#!/usr/bin/env python3\n", "\n", "# Zuordnung_Ungarisch.py\n", "\"\"\"\n", "Kapitel Graphen: Das Zuordnungsproblem, dreifach geloest.\n", "\n", " (1) Ungarischer Algorithmus (scipy.optimize.linear_sum_assignment) - O(n^3)\n", " (2) als LP OHNE Ganzzahligkeitsforderung -> liefert trotzdem 0/1 (Birkhoff)\n", " (3) als MILP MIT Ganzzahligkeitsforderung -> gleiches Ergebnis, mehr Aufwand\n", "\n", "Zeigt damit die praktische Bedeutung der totalen Unimodularitaet.\n", "\"\"\"\n", "\n", "import time\n", "\n", "import numpy as np\n", "from scipy.optimize import linear_sum_assignment, linprog\n", "\n", "\n", "def erzeuge_kosten(n, seed=11):\n", " rng = np.random.default_rng(seed)\n", " return rng.integers(10, 99, size=(n, n)).astype(float)\n", "\n", "\n", "def loese_ungarisch(kosten):\n", " zeilen, spalten = linear_sum_assignment(kosten)\n", " return kosten[zeilen, spalten].sum(), spalten\n", "\n", "\n", "def baue_lp(kosten):\n", " \"\"\"Gemeinsame LP-Struktur fuer Variante 2 und 3.\"\"\"\n", " n = len(kosten)\n", " c = kosten.flatten() # x_ij in Zeilenreihenfolge\n", " A_eq = np.zeros((2 * n, n * n))\n", " for i in range(n): # jede Person genau eine Aufgabe\n", " A_eq[i, i * n:(i + 1) * n] = 1.0\n", " for j in range(n): # jede Aufgabe genau einer Person\n", " A_eq[n + j, j::n] = 1.0\n", " b_eq = np.ones(2 * n)\n", " return c, A_eq, b_eq\n", "\n", "\n", "def loese_lp(kosten, ganzzahlig):\n", " n = len(kosten)\n", " c, A_eq, b_eq = baue_lp(kosten)\n", " ergebnis = linprog(c=c, A_eq=A_eq, b_eq=b_eq, bounds=[(0, 1)] * (n * n),\n", " integrality=np.ones(n * n) if ganzzahlig else None,\n", " method=\"highs\")\n", " x = ergebnis.x.reshape(n, n)\n", " return ergebnis.fun, x\n", "\n", "\n", "if __name__ == \"__main__\":\n", " print(\"=\" * 84)\n", " print(\" ZUORDNUNGSPROBLEM: DREI WEGE ZUM SELBEN ERGEBNIS\")\n", " print(\"=\" * 84)\n", "\n", " # --- Kleines Beispiel zum Nachvollziehen ------------------------------\n", " kosten = np.array([[82., 83., 69., 92.],\n", " [77., 37., 49., 92.],\n", " [11., 69., 5., 86.],\n", " [8., 9., 98., 23.]])\n", " namen = [\"Anna\", \"Ben\", \"Carla\", \"David\"]\n", " aufgaben = [\"Auftrag W\", \"Auftrag X\", \"Auftrag Y\", \"Auftrag Z\"]\n", "\n", " print(\"\\nKostenmatrix (wer bearbeitet was zu welchen Kosten?):\")\n", " print(f\"{'':<8}\" + \"\".join(f\"{a:>12}\" for a in aufgaben))\n", " for i, name in enumerate(namen):\n", " print(f\"{name:<8}\" + \"\".join(f\"{kosten[i, j]:>12.0f}\" for j in range(4)))\n", "\n", " wert, zuordnung = loese_ungarisch(kosten)\n", " print(f\"\\nOptimale Zuordnung (Gesamtkosten {wert:.0f}):\")\n", " for i, j in enumerate(zuordnung):\n", " print(f\" {namen[i]:<8} -> {aufgaben[j]:<12} ({kosten[i, j]:.0f} EUR)\")\n", "\n", " # --- Nachweis: LP ohne Ganzzahligkeit liefert trotzdem 0/1 -----------\n", " wert_lp, x_lp = loese_lp(kosten, ganzzahlig=False)\n", " ist_binaer = np.all((np.abs(x_lp) < 1e-9) | (np.abs(x_lp - 1) < 1e-9))\n", " print(f\"\\nLP OHNE Ganzzahligkeitsforderung: Kosten {wert_lp:.0f}, \"\n", " f\"Loesung ist {'0/1-wertig' if ist_binaer else 'GEBROCHEN'}\")\n", " print(\" -> Satz von Birkhoff/von Neumann bestaetigt: Die Ecken sind Permutationen.\")\n", "\n", " # --- Laufzeitvergleich bei wachsender Groesse ------------------------\n", " print(\"\\n\" + \"-\" * 84)\n", " print(f\"{'n':>4} | {'Ungarisch':>12} | {'LP (kontinuierlich)':>21} | \"\n", " f\"{'MILP (ganzzahlig)':>19} | {'gleich?':>8}\")\n", " print(\"-\" * 84)\n", " for n in [10, 25, 50, 100]:\n", " k = erzeuge_kosten(n)\n", "\n", " t0 = time.perf_counter(); w1, _ = loese_ungarisch(k); t1 = time.perf_counter() - t0\n", " t0 = time.perf_counter(); w2, _ = loese_lp(k, False); t2 = time.perf_counter() - t0\n", " if n <= 50:\n", " t0 = time.perf_counter(); w3, _ = loese_lp(k, True); t3 = time.perf_counter() - t0\n", " t3_text, gleich = f\"{t3*1000:>16.1f} ms\", abs(w1 - w3) < 1e-6\n", " else:\n", " t3_text, gleich = f\"{'uebersprungen':>19}\", abs(w1 - w2) < 1e-6\n", "\n", " print(f\"{n:>4} | {t1*1000:>9.1f} ms | {t2*1000:>18.1f} ms | {t3_text} | \"\n", " f\"{'ja' if gleich else 'NEIN':>8}\")\n", "\n", " print(\"-\" * 84)\n", " print(\"Fazit: Der spezialisierte Ungarische Algorithmus ist um Groessenordnungen\")\n", " print(\"schneller. Nutzen Sie fuer reine Zuordnungen NIE einen MILP-Solver.\")\n", " print(\"=\" * 84)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Praxisbeispiel: Flotten-Routing\n", "\n", "`VRP_Flotten_Routing.py`\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#!/usr/bin/env python3\n", "\n", "# VRP_Flotten_Routing.py\n", "\"\"\"\n", "Kapitel Graphen: Capacitated Vehicle Routing Problem with Time Windows (CVRPTW)\n", "mit der Routing-Bibliothek von Google OR-Tools.\n", "\n", "Eigenschaften:\n", " * Eingabedaten werden vorab auf Plausibilitaet geprueft (Kapazitaet\n", " ausreichend? Zeitfenster erreichbar?)\n", " * Fahrzeit und Servicezeit werden getrennt ausgewiesen\n", " * Ausgabe als lesbarer Tourenplan mit Ankunftszeiten\n", " * Kennzahlen: Auslastung, Leerfahrten, Wartezeit\n", "\"\"\"\n", "\n", "import numpy as np\n", "from ortools.constraint_solver import pywrapcp, routing_enums_pb2\n", "\n", "SERVICEZEIT = 10 # Minuten je Kundenstopp\n", "WARTEZEIT_MAX = 60 # zulaessige Wartezeit bei zu frueher Ankunft\n", "SCHICHTLAENGE = 600 # Minuten\n", "\n", "\n", "def erzeuge_daten(seed: int = 42):\n", " \"\"\"Synthetische, aber reproduzierbare Instanz: 1 Depot + 16 Kunden.\"\"\"\n", " anzahl_orte = 17\n", " rng = np.random.default_rng(seed)\n", " koordinaten = rng.random((anzahl_orte, 2)) * 100 # 100 x 100 km Raster\n", "\n", " distanz = np.zeros((anzahl_orte, anzahl_orte), dtype=int)\n", " for i in range(anzahl_orte):\n", " for j in range(anzahl_orte):\n", " distanz[i][j] = int(np.linalg.norm(koordinaten[i] - koordinaten[j]))\n", "\n", " return {\n", " \"distanzmatrix\": distanz.tolist(),\n", " \"zeitfenster\": [\n", " (0, SCHICHTLAENGE), # 0: Depot\n", " (30, 120), (60, 180), (100, 240), (150, 300), # Kunden 1-4\n", " (60, 180), (120, 240), (200, 360), (300, 450), # Kunden 5-8\n", " (180, 300), (240, 360), (300, 480), (360, 500), # Kunden 9-12\n", " (60, 200), (120, 300), (240, 400), (300, 550), # Kunden 13-16\n", " ],\n", " \"bedarfe\": [0, 2, 3, 1, 4, 2, 2, 3, 1, 2, 4, 3, 2, 1, 2, 3, 2],\n", " \"kapazitaeten\": [10, 10, 10, 10],\n", " \"anzahl_fahrzeuge\": 4,\n", " \"depot\": 0,\n", " }\n", "\n", "\n", "def pruefe_daten(daten) -> None:\n", " \"\"\"Vorabdiagnose - fangt die haeufigsten Ursachen fuer 'keine Loesung' ab.\"\"\"\n", " gesamtbedarf = sum(daten[\"bedarfe\"])\n", " gesamtkapazitaet = sum(daten[\"kapazitaeten\"])\n", " print(f\"Gesamtbedarf {gesamtbedarf} Einheiten | \"\n", " f\"Flottenkapazitaet {gesamtkapazitaet} Einheiten | \"\n", " f\"Auslastung {gesamtbedarf / gesamtkapazitaet * 100:.0f} %\")\n", " if gesamtbedarf > gesamtkapazitaet:\n", " raise SystemExit(\"UNLOESBAR: Der Bedarf uebersteigt die Flottenkapazitaet.\")\n", "\n", " d = daten[\"distanzmatrix\"]\n", " for kunde, (fruehestens, spaetestens) in enumerate(daten[\"zeitfenster\"]):\n", " if kunde == 0:\n", " continue\n", " direktfahrt = d[0][kunde]\n", " if direktfahrt > spaetestens:\n", " raise SystemExit(\n", " f\"UNLOESBAR: Kunde {kunde} ist erst nach {direktfahrt} min erreichbar, \"\n", " f\"sein Zeitfenster endet aber bei {spaetestens} min.\")\n", " print(\"Vorabpruefung bestanden: Kapazitaet und Zeitfenster sind grundsaetzlich machbar.\")\n", "\n", "\n", "def loese_cvrptw(zeitlimit_s: int = 5):\n", " daten = erzeuge_daten()\n", " pruefe_daten(daten)\n", "\n", " manager = pywrapcp.RoutingIndexManager(\n", " len(daten[\"distanzmatrix\"]), daten[\"anzahl_fahrzeuge\"], daten[\"depot\"])\n", " routing = pywrapcp.RoutingModel(manager)\n", "\n", " # --- Fahrzeit + Servicezeit als Kantengewicht ------------------------\n", " def zeit_callback(von_index, nach_index):\n", " von = manager.IndexToNode(von_index)\n", " nach = manager.IndexToNode(nach_index)\n", " service = SERVICEZEIT if von != daten[\"depot\"] else 0\n", " return daten[\"distanzmatrix\"][von][nach] + service\n", "\n", " zeit_index = routing.RegisterTransitCallback(zeit_callback)\n", " routing.SetArcCostEvaluatorOfAllVehicles(zeit_index)\n", "\n", " # --- Kapazitaetsdimension ---------------------------------------------\n", " def bedarf_callback(von_index):\n", " return daten[\"bedarfe\"][manager.IndexToNode(von_index)]\n", "\n", " bedarf_index = routing.RegisterUnaryTransitCallback(bedarf_callback)\n", " routing.AddDimensionWithVehicleCapacity(\n", " bedarf_index, 0, daten[\"kapazitaeten\"], True, \"Kapazitaet\")\n", "\n", " # --- Zeitdimension mit Zeitfenstern -----------------------------------\n", " routing.AddDimension(zeit_index, WARTEZEIT_MAX, SCHICHTLAENGE, False, \"Zeit\")\n", " zeit_dimension = routing.GetDimensionOrDie(\"Zeit\")\n", " for ort, (fruehestens, spaetestens) in enumerate(daten[\"zeitfenster\"]):\n", " zeit_dimension.CumulVar(manager.NodeToIndex(ort)).SetRange(fruehestens, spaetestens)\n", "\n", " # --- Suchparameter -----------------------------------------------------\n", " parameter = pywrapcp.DefaultRoutingSearchParameters()\n", " parameter.first_solution_strategy = (\n", " routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)\n", " parameter.local_search_metaheuristic = (\n", " routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)\n", " parameter.time_limit.seconds = zeitlimit_s\n", "\n", " loesung = routing.SolveWithParameters(parameter)\n", " if not loesung:\n", " print(\"Keine zulaessige Routenfuehrung gefunden.\")\n", " return\n", "\n", " # --- Auswertung --------------------------------------------------------\n", " print(\"\\n\" + \"=\" * 84)\n", " print(\" OPTIMIERTER TOURENPLAN (CVRPTW)\")\n", " print(\"=\" * 84)\n", "\n", " gesamtzeit = gesamtfracht = gesamtdistanz = 0\n", " kapazitaet = daten[\"kapazitaeten\"]\n", "\n", " for fahrzeug in range(daten[\"anzahl_fahrzeuge\"]):\n", " index = routing.Start(fahrzeug)\n", " if routing.IsEnd(loesung.Value(routing.NextVar(index))):\n", " print(f\"\\nFahrzeug {fahrzeug + 1}: nicht eingesetzt\")\n", " continue\n", "\n", " stationen, fracht, distanz = [], 0, 0\n", " while not routing.IsEnd(index):\n", " knoten = manager.IndexToNode(index)\n", " ankunft = loesung.Min(zeit_dimension.CumulVar(index))\n", " fracht += daten[\"bedarfe\"][knoten]\n", " bezeichnung = \"Depot\" if knoten == 0 else f\"K{knoten}\"\n", " stationen.append(f\"{bezeichnung}@{ankunft}\")\n", " naechster = loesung.Value(routing.NextVar(index))\n", " distanz += daten[\"distanzmatrix\"][knoten][manager.IndexToNode(naechster)]\n", " index = naechster\n", "\n", " endzeit = loesung.Min(zeit_dimension.CumulVar(index))\n", " stationen.append(f\"Depot@{endzeit}\")\n", " gesamtzeit += endzeit\n", " gesamtfracht += fracht\n", " gesamtdistanz += distanz\n", "\n", " print(f\"\\nFahrzeug {fahrzeug + 1}:\")\n", " print(\" \" + \" -> \".join(stationen))\n", " print(f\" Schichtzeit {endzeit} min | Fahrstrecke {distanz} km | \"\n", " f\"Fracht {fracht}/{kapazitaet[fahrzeug]} \"\n", " f\"({fracht / kapazitaet[fahrzeug] * 100:.0f} % Auslastung)\")\n", "\n", " print(\"\\n\" + \"-\" * 84)\n", " print(f\"Summe Schichtzeiten: {gesamtzeit} min\")\n", " print(f\"Summe Fahrstrecken: {gesamtdistanz} km\")\n", " print(f\"Transportierte Fracht: {gesamtfracht} von {sum(daten['bedarfe'])} Einheiten\")\n", " assert gesamtfracht == sum(daten[\"bedarfe\"]), \"Nicht alle Kunden wurden beliefert!\"\n", " print(\"Alle Kunden wurden innerhalb ihrer Zeitfenster beliefert.\")\n", " print(\"=\" * 84)\n", "\n", "\n", "if __name__ == \"__main__\":\n", " loese_cvrptw()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Finde den Denkfehler\n", "\n", "`VRP_Kapazitaetsfalle.py`\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#!/usr/bin/env python3\n", "\n", "# VRP_Kapazitaetsfalle.py\n", "\"\"\"\n", "Kapitel Graphen: Die vergessene Dimension.\n", "\n", "Die Routing-Bibliothek von OR-Tools kennt keine \"Kapazitaet\" von sich aus.\n", "Sie kennt nur DIMENSIONEN - benannte Groessen, die sich entlang einer Tour\n", "aufsummieren und begrenzt werden koennen. Distanz ist eine, Zeit ist eine,\n", "Ladung ist eine. Wer eine davon nicht anlegt, bekommt trotzdem eine Loesung:\n", "eine schoene, kurze, guenstige - und unfahrbare.\n", "\n", "Dieses Programm loest dieselbe Instanz zweimal und prueft beide Ergebnisse\n", "gegen die tatsaechlichen Lademengen.\n", "\n", "Instanz: 1 Depot, 16 Kunden, 4 Fahrzeuge zu je 10 Paletten.\n", "Gesamtbedarf 37 Paletten bei 40 Paletten Flottenkapazitaet - es ist also\n", "knapp, aber machbar.\n", "\n", "Benoetigt: numpy, ortools\n", "\"\"\"\n", "\n", "from __future__ import annotations\n", "\n", "import numpy as np\n", "from ortools.constraint_solver import pywrapcp, routing_enums_pb2\n", "\n", "# Dieselbe Instanz wie VRP_Flotten_Routing.py\n", "BEDARFE = [0, 2, 3, 1, 4, 2, 2, 3, 1, 2, 4, 3, 2, 1, 2, 3, 2]\n", "KAPAZITAETEN = [10, 10, 10, 10]\n", "ANZAHL_FAHRZEUGE = 4\n", "DEPOT = 0\n", "\n", "\n", "def distanzmatrix(seed: int = 42) -> list[list[int]]:\n", " rng = np.random.default_rng(seed)\n", " koordinaten = rng.random((len(BEDARFE), 2)) * 100 # 100 x 100 km\n", " n = len(BEDARFE)\n", " return [[int(np.linalg.norm(koordinaten[i] - koordinaten[j]))\n", " for j in range(n)] for i in range(n)]\n", "\n", "\n", "def plane(mit_kapazitaet: bool, zeitlimit: int = 5) -> dict:\n", " \"\"\"Loest die Tourenplanung - wahlweise mit oder ohne Ladungsdimension.\"\"\"\n", " distanz = distanzmatrix()\n", " manager = pywrapcp.RoutingIndexManager(len(distanz), ANZAHL_FAHRZEUGE, DEPOT)\n", " routing = pywrapcp.RoutingModel(manager)\n", "\n", " def entfernung(von_index, nach_index):\n", " return distanz[manager.IndexToNode(von_index)][manager.IndexToNode(nach_index)]\n", "\n", " kosten_id = routing.RegisterTransitCallback(entfernung)\n", " routing.SetArcCostEvaluatorOfAllVehicles(kosten_id)\n", "\n", " # DIE entscheidende Stelle. Ohne diesen Block existiert im Modell keine\n", " # Ladung - die Fahrzeuge sind dann unendlich gross.\n", " if mit_kapazitaet:\n", " def bedarf(index):\n", " return BEDARFE[manager.IndexToNode(index)]\n", "\n", " bedarf_id = routing.RegisterUnaryTransitCallback(bedarf)\n", " routing.AddDimensionWithVehicleCapacity(\n", " bedarf_id,\n", " 0, # kein Zwischenpuffer\n", " KAPAZITAETEN, # Obergrenze je Fahrzeug\n", " True, # Ladung startet bei 0\n", " \"Ladung\")\n", "\n", " parameter = pywrapcp.DefaultRoutingSearchParameters()\n", " parameter.first_solution_strategy = (\n", " routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)\n", " parameter.local_search_metaheuristic = (\n", " routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)\n", " parameter.time_limit.FromSeconds(zeitlimit)\n", "\n", " loesung = routing.SolveWithParameters(parameter)\n", " if loesung is None:\n", " raise RuntimeError(\"Keine Loesung gefunden\")\n", "\n", " touren, strecken, ladungen = [], [], []\n", " for fahrzeug in range(ANZAHL_FAHRZEUGE):\n", " index = routing.Start(fahrzeug)\n", " tour, strecke, ladung = [], 0, 0\n", " while not routing.IsEnd(index):\n", " knoten = manager.IndexToNode(index)\n", " tour.append(knoten)\n", " ladung += BEDARFE[knoten]\n", " vorher = index\n", " index = loesung.Value(routing.NextVar(index))\n", " strecke += routing.GetArcCostForVehicle(vorher, index, fahrzeug)\n", " tour.append(manager.IndexToNode(index))\n", " touren.append(tour)\n", " strecken.append(strecke)\n", " ladungen.append(ladung)\n", "\n", " return {\"touren\": touren, \"strecken\": strecken, \"ladungen\": ladungen,\n", " \"gesamtstrecke\": sum(strecken)}\n", "\n", "\n", "def pruefe(ergebnis: dict) -> list[str]:\n", " \"\"\"Prueft den Plan gegen die Wirklichkeit - unabhaengig vom Modell.\n", "\n", " Genau diese Trennung ist der Punkt: Die Pruefung darf nicht dieselben\n", " Annahmen benutzen wie das Modell, sonst prueft sie nichts.\n", " \"\"\"\n", " beanstandungen = []\n", " for fahrzeug, (ladung, kapazitaet) in enumerate(\n", " zip(ergebnis[\"ladungen\"], KAPAZITAETEN)):\n", " if ladung > kapazitaet:\n", " beanstandungen.append(\n", " f\"Fahrzeug {fahrzeug + 1}: {ladung} Paletten geladen, \"\n", " f\"Kapazitaet {kapazitaet} ({ladung - kapazitaet} zu viel)\")\n", "\n", " beliefert = sorted(k for tour in ergebnis[\"touren\"] for k in tour[1:-1])\n", " erwartet = list(range(1, len(BEDARFE)))\n", " if beliefert != erwartet:\n", " fehlend = set(erwartet) - set(beliefert)\n", " if fehlend:\n", " beanstandungen.append(f\"nicht beliefert: {sorted(fehlend)}\")\n", " return beanstandungen\n", "\n", "\n", "def zeige(titel: str, ergebnis: dict) -> None:\n", " print(f\"\\n{titel}\")\n", " print(f\" Gesamtstrecke {ergebnis['gesamtstrecke']} km\")\n", " print(f\" {'Fahrzeug':<10} {'Stopps':>7} {'Strecke':>9} {'Ladung':>8} \"\n", " f\"{'Kapazitaet':>11}\")\n", " for i, (tour, strecke, ladung) in enumerate(\n", " zip(ergebnis[\"touren\"], ergebnis[\"strecken\"], ergebnis[\"ladungen\"])):\n", " markierung = \" <-- ueberladen\" if ladung > KAPAZITAETEN[i] else \"\"\n", " print(f\" {i + 1:<10} {len(tour) - 2:>7} {strecke:>8} km {ladung:>8} \"\n", " f\"{KAPAZITAETEN[i]:>11}{markierung}\")\n", "\n", " beanstandungen = pruefe(ergebnis)\n", " if beanstandungen:\n", " print(\" PRUEFUNG: DURCHGEFALLEN\")\n", " for text in beanstandungen:\n", " print(f\" - {text}\")\n", " else:\n", " print(\" PRUEFUNG: bestanden\")\n", "\n", "\n", "if __name__ == \"__main__\":\n", " print(\"=\" * 78)\n", " print(\" DIE VERGESSENE DIMENSION\")\n", " print(\"=\" * 78)\n", " print(f\"16 Kunden, Gesamtbedarf {sum(BEDARFE)} Paletten, \"\n", " f\"{ANZAHL_FAHRZEUGE} Fahrzeuge zu je {KAPAZITAETEN[0]} \"\n", " f\"= {sum(KAPAZITAETEN)} Paletten Flottenkapazitaet.\")\n", "\n", " ohne = plane(mit_kapazitaet=False)\n", " zeige(\"[1] Ohne Ladungsdimension\", ohne)\n", "\n", " mit = plane(mit_kapazitaet=True)\n", " zeige(\"[2] Mit AddDimensionWithVehicleCapacity\", mit)\n", "\n", " print(\"\\n\" + \"=\" * 78)\n", " mehr = mit[\"gesamtstrecke\"] - ohne[\"gesamtstrecke\"]\n", " print(f\"Der korrekte Plan ist {mehr} km laenger \"\n", " f\"({mehr / ohne['gesamtstrecke'] * 100:.1f} %).\")\n", " print()\n", " print(\"Und genau darin liegt die Gefahr: Lauf [1] sieht BESSER aus. Wer\")\n", " print(\"beide Zahlen nebeneinander legt, ohne die Ladung zu pruefen, haelt\")\n", " print(\"die unfahrbare Loesung fuer die bessere Optimierung - und den\")\n", " print(\"korrekten Plan fuer schlechte Arbeit.\")\n", " print()\n", " print(\"Die Routing-Bibliothek kennt keine 'Kapazitaet'. Sie kennt nur\")\n", " print(\"Dimensionen, die man ihr anlegt. Was nicht als Dimension existiert,\")\n", " print(\"wird nicht begrenzt - und faellt niemandem auf, weil das Ergebnis\")\n", " print(\"plausibel aussieht.\")\n", " print(\"=\" * 78)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.11" } }, "nbformat": 4, "nbformat_minor": 5 }