operations_research/OR_HTML_04/Notebooks_04/graphen.ipynb

{
 "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
}