operations_research/bilder_04/erzeuge_vrp_touren.py

#!/usr/bin/env python3

# erzeuge_vrp_touren.py
"""
Erzeugt die Tourenkarte zum Graphen-Kapitel:

    bilder_04/kap_graphen_vrp_touren.svg

Das Vorgaengerbild war ein Schema: zwei Touren, vier Kunden, Zeitfenster in
Uhrzeiten. Die Instanz des Kapitels hat 16 Kunden, vier Fahrzeuge und
Zeitfenster in Minuten - das Bild zeigte also die Bauform, nicht den Fall.

Hier wird die Instanz aus VRP_Flotten_Routing.py geloest und die Loesung als
Karte gezeichnet: Depot, Kunden (Punktgroesse = Bedarf), die vier Touren in
vier Farben, und in der Legende je Fahrzeug Strecke und Auslastung.

ZUR REPRODUZIERBARKEIT: Das Buchprogramm gibt der Suche ein ZEITBUDGET von
fuenf Sekunden. Ein Zeitbudget ist auf einer anderen Maschine ein anderes
Rechenbudget, das Ergebnis also nicht garantiert identisch. Fuer das Bild wird
deshalb zusaetzlich eine feste Loesungszahl als Abbruchkriterium gesetzt und
die Gesamtstrecke ausgegeben - weicht sie vom abgedruckten Lauf ab, faellt es
sofort auf, statt still ins Buch zu wandern.

Aufruf (aus dem Repository-Wurzelverzeichnis):
    python3 bilder_04/erzeuge_vrp_touren.py

Benoetigt: numpy, matplotlib, ortools
"""

from __future__ import annotations

import os
import re

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from ortools.constraint_solver import pywrapcp, routing_enums_pb2

# Reproduzierbare SVG-Ausgabe (siehe erzeuge_titelseite.py)
plt.rcParams["svg.hashsalt"] = "or-mit-python-v04"

BASIS = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
BILDER = os.path.join(BASIS, "bilder_04")

TOURFARBEN = ["#4338ca", "#0891b2", "#b45309", "#059669"]

# --- Instanz: wortgleich zu VRP_Flotten_Routing.py --------------------------
SERVICEZEIT = 10
WARTEZEIT_MAX = 60
SCHICHTLAENGE = 600
SAAT = 42
ZEITBUDGET_S = 5
LOESUNGSGRENZE = 300          # zusaetzliches, maschinenunabhaengiges Kriterium

ZEITFENSTER = [
    (0, SCHICHTLAENGE),
    (30, 120),  (60, 180),  (100, 240), (150, 300),
    (60, 180),  (120, 240), (200, 360), (300, 450),
    (180, 300), (240, 360), (300, 480), (360, 500),
    (60, 200),  (120, 300), (240, 400), (300, 550),
]
BEDARFE = [0, 2, 3, 1, 4, 2, 2, 3, 1, 2, 4, 3, 2, 1, 2, 3, 2]
KAPAZITAETEN = [10, 10, 10, 10]


def erzeuge_daten():
    anzahl_orte = len(ZEITFENSTER)
    rng = np.random.default_rng(SAAT)
    koordinaten = rng.random((anzahl_orte, 2)) * 100
    distanz = np.zeros((anzahl_orte, anzahl_orte), dtype=int)
    for i in range(anzahl_orte):
        for j in range(anzahl_orte):
            distanz[i][j] = int(np.linalg.norm(koordinaten[i] - koordinaten[j]))
    return koordinaten, distanz


def loese(distanz):
    """CVRPTW mit denselben Einstellungen wie im Buchprogramm."""
    manager = pywrapcp.RoutingIndexManager(len(distanz), len(KAPAZITAETEN), 0)
    routing = pywrapcp.RoutingModel(manager)

    def zeit_callback(von_index, nach_index):
        von = manager.IndexToNode(von_index)
        nach = manager.IndexToNode(nach_index)
        service = SERVICEZEIT if von != 0 else 0
        return int(distanz[von][nach]) + service

    zeit_index = routing.RegisterTransitCallback(zeit_callback)
    routing.SetArcCostEvaluatorOfAllVehicles(zeit_index)

    def bedarf_callback(von_index):
        return BEDARFE[manager.IndexToNode(von_index)]

    bedarf_index = routing.RegisterUnaryTransitCallback(bedarf_callback)
    routing.AddDimensionWithVehicleCapacity(
        bedarf_index, 0, KAPAZITAETEN, True, "Kapazitaet")

    routing.AddDimension(zeit_index, WARTEZEIT_MAX, SCHICHTLAENGE, False, "Zeit")
    zeit_dimension = routing.GetDimensionOrDie("Zeit")
    for ort, (fruehestens, spaetestens) in enumerate(ZEITFENSTER):
        zeit_dimension.CumulVar(manager.NodeToIndex(ort)).SetRange(
            fruehestens, spaetestens)

    parameter = pywrapcp.DefaultRoutingSearchParameters()
    parameter.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)
    parameter.local_search_metaheuristic = (
        routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)
    parameter.time_limit.seconds = ZEITBUDGET_S
    parameter.solution_limit = LOESUNGSGRENZE

    loesung = routing.SolveWithParameters(parameter)
    if not loesung:
        raise SystemExit("Keine zulaessige Routenfuehrung gefunden.")

    touren = []
    for fahrzeug in range(len(KAPAZITAETEN)):
        index = routing.Start(fahrzeug)
        if routing.IsEnd(loesung.Value(routing.NextVar(index))):
            continue
        stationen, fracht, strecke = [], 0, 0
        while not routing.IsEnd(index):
            knoten = manager.IndexToNode(index)
            stationen.append(knoten)
            fracht += BEDARFE[knoten]
            vorher = index
            index = loesung.Value(routing.NextVar(index))
            strecke += int(distanz[manager.IndexToNode(vorher)][
                manager.IndexToNode(index)])
        stationen.append(0)
        touren.append({"stationen": stationen, "fracht": fracht,
                       "strecke": strecke,
                       "kapazitaet": KAPAZITAETEN[fahrzeug]})
    return touren


def zeichne(koordinaten, touren) -> None:
    figur, achse = plt.subplots(figsize=(7.6, 6.2))
    gesamt = sum(t["strecke"] for t in touren)

    for nummer, (tour, farbe) in enumerate(zip(touren, TOURFARBEN), start=1):
        weg = koordinaten[tour["stationen"]]
        achse.plot(weg[:, 0], weg[:, 1], "-", color=farbe, linewidth=1.8,
                   alpha=0.85, zorder=2,
                   label=f"Fahrzeug {nummer}: {tour['strecke']} km, "
                         f"{tour['fracht']}/{tour['kapazitaet']} Einheiten")

    for ort in range(1, len(koordinaten)):
        achse.scatter(*koordinaten[ort], s=30 + 26 * BEDARFE[ort],
                      color="#334155", zorder=3)
        achse.annotate(f"K{ort}", koordinaten[ort], textcoords="offset points",
                       xytext=(7, 5), fontsize=7.5, color="#475569")

    achse.scatter(*koordinaten[0], marker="s", s=190, color="#be123c", zorder=4)
    achse.annotate("Depot", koordinaten[0], textcoords="offset points",
                   xytext=(9, 7), fontsize=9.5, color="#be123c",
                   fontweight="bold")

    achse.set_xlabel("km")
    achse.set_ylabel("km")
    achse.set_title(f"Vier Touren für 16 Kunden mit Zeitfenstern — "
                    f"{gesamt} km insgesamt\n"
                    "Punktgröße = Liefermenge", fontsize=10.5)
    achse.grid(linestyle=":", alpha=0.4)
    achse.legend(fontsize=8, loc="upper center", bbox_to_anchor=(0.5, -0.11),
                 ncol=2, frameon=False)
    for seite in ("top", "right"):
        achse.spines[seite].set_visible(False)
    figur.tight_layout()

    os.makedirs(BILDER, exist_ok=True)
    for endung in ("svg",):
        pfad = os.path.join(BILDER, f"kap_graphen_vrp_touren.{endung}")
        figur.savefig(pfad, format=endung, dpi=160, bbox_inches="tight",
                      metadata={"Date": None} if endung == "svg" else None)
        print(f"geschrieben: {pfad}")
    plt.close(figur)


# --- Karte ------------------------------------------------------------------
#
# Dieselbe Loesung ein zweites Mal, diesmal auf einem echten Kartenhintergrund.
# Das war lange als "blockiert" vermerkt, weil die Instanz ein synthetisches
# 100 x 100-km-Raster ohne Geokoordinaten ist: Eine Karte schien entweder eine
# neue Instanz zu verlangen (und damit zwei Dutzend neue abgedruckte Zahlen)
# oder eine willkuerliche Verankerung.
#
# Das war eine falsche Alternative. Die azimutal-aequidistante Projektion um
# den Schwerpunkt des Rasters bildet die Punkte SO auf Laenge und Breite ab,
# dass die Grosskreisdistanzen mit den Rasterdistanzen uebereinstimmen - ueber
# alle 136 Punktpaare bis auf 0,001 %. Das Modell bleibt unberuehrt, kein
# abgedruckter Wert aendert sich, und die Karte zeigt trotzdem echte Lagen.
# pruefe_projektion() unten weist das bei jedem Lauf nach.
ERDRADIUS_KM = 6371.0088
ANKER_LAT, ANKER_LON = 51.00, 9.20      # Raum Kassel; zwei Zahlen genuegen,
                                        # um die Instanz woanders hinzulegen
KARTEN = os.path.join(BILDER, "karten")
KARTENNAME = "kap_graphen_vrp_touren"


def projiziere(koordinaten):
    """Rasterkoordinaten (km) -> (Breite, Laenge), azimutal-aequidistant um den
    Schwerpunkt. Formel: Zielpunkt aus Startpunkt, Azimut und Bogendistanz."""
    mitte = koordinaten.mean(axis=0)
    dx, dy = (koordinaten - mitte).T
    bogen = np.hypot(dx, dy) / ERDRADIUS_KM        # Winkeldistanz im Bogenmass
    azimut = np.arctan2(dx, dy)                    # 0 = Nord, im Uhrzeigersinn
    phi0 = np.radians(ANKER_LAT)
    breite = np.arcsin(np.sin(phi0) * np.cos(bogen)
                       + np.cos(phi0) * np.sin(bogen) * np.cos(azimut))
    laenge = np.radians(ANKER_LON) + np.arctan2(
        np.sin(azimut) * np.sin(bogen) * np.cos(phi0),
        np.cos(bogen) - np.sin(phi0) * np.sin(breite))
    return np.degrees(breite), np.degrees(laenge)


def pruefe_projektion(koordinaten, breite, laenge) -> float:
    """Vergleicht jede Grosskreisdistanz mit der Rasterdistanz und bricht ab,
    wenn die Karte die Instanz verzerrt darstellen wuerde. Die Zahl steht im
    Buch - sie wird hier bewiesen, nicht behauptet."""
    def grosskreis(i, j):
        p1, p2 = np.radians(breite[i]), np.radians(breite[j])
        dl = np.radians(laenge[j] - laenge[i])
        h = (np.sin((p2 - p1) / 2) ** 2
             + np.cos(p1) * np.cos(p2) * np.sin(dl / 2) ** 2)
        return 2 * ERDRADIUS_KM * np.arcsin(np.sqrt(h))

    schlimmster = 0.0
    for i in range(len(koordinaten)):
        for j in range(i + 1, len(koordinaten)):
            raster = float(np.linalg.norm(koordinaten[i] - koordinaten[j]))
            if raster > 0:
                schlimmster = max(schlimmster,
                                  abs(grosskreis(i, j) - raster) / raster * 100)
    if schlimmster > 0.01:
        raise SystemExit(
            f"Karte: Die Projektion verzerrt die Instanz um bis zu "
            f"{schlimmster:.3f} % - erwartet sind hoechstens 0,01 %. "
            f"Anker oder Projektion pruefen.")
    return schlimmster


# Von Foliums zehn CDN-Einbindungen braucht diese Karte genau zwei: leaflet.js
# und leaflet.css. jQuery, Bootstrap, Glyphicons, FontAwesome und
# awesome-markers gehoeren zu Foliums Standardmarkern - hier kommen nur
# CircleMarker, PolyLine und Tooltip vor, die reines Leaflet sind.
LEAFLET_JS = "https://cdn.jsdelivr.net/npm/leaflet@1.9.3/dist/leaflet.js"
LEAFLET_CSS = "https://cdn.jsdelivr.net/npm/leaflet@1.9.3/dist/leaflet.css"
EINBINDUNG_RE = re.compile(
    r'[ \t]*<(?:script src|link rel="stylesheet" href)="(https?://[^"]+)"'
    r'[^>]*>(?:</script>)?\n?')
KENNUNG_RE = re.compile(r"\b([a-z_]+)_([0-9a-f]{32})\b")


def normalisiere(html: str) -> str:
    """Macht Foliums Ausgabe byteweise reproduzierbar und netzunabhaengig.

    Zwei Eingriffe:
    1. Folium vergibt jedem Element eine uuid4 - jeder Lauf erzeugte sonst eine
       andere Datei und einen unlesbaren git-diff. Dieselbe Sorgfalt wie
       svg.hashsalt bei den matplotlib-Bildern.
    2. Die beiden gebrauchten Leaflet-Dateien werden lokal eingebunden (sie
       liegen in web_04/assets/ und landen beim Bau in OR_HTML_04/assets/),
       die uebrigen acht Einbindungen entfallen. Ohne Netz fehlen dann nur die
       Kartenkacheln - Marker und Touren bleiben sichtbar."""
    reihenfolge: dict[str, str] = {}

    def kennung(m: re.Match) -> str:
        hexwert = m.group(2)
        if hexwert not in reihenfolge:
            reihenfolge[hexwert] = f"{len(reihenfolge) + 1:02d}"
        return f"{m.group(1)}_{reihenfolge[hexwert]}"

    html = KENNUNG_RE.sub(kennung, html)

    def einbindung(m: re.Match) -> str:
        url = m.group(1)
        if url == LEAFLET_JS:
            return '    <script src="../../assets/leaflet.js"></script>\n'
        if url == LEAFLET_CSS:
            return ('    <link rel="stylesheet" '
                    'href="../../assets/leaflet.css"/>\n')
        return ""

    return EINBINDUNG_RE.sub(einbindung, html)


def zeichne_karte(koordinaten, touren, abweichung: float) -> None:
    import folium

    breite, laenge = projiziere(koordinaten)
    karte = folium.Map(location=[float(breite.mean()), float(laenge.mean())],
                       zoom_start=9, tiles="OpenStreetMap",
                       control_scale=True)

    for nummer, (tour, farbe) in enumerate(zip(touren, TOURFARBEN), start=1):
        weg = [[float(breite[s]), float(laenge[s])] for s in tour["stationen"]]
        folium.PolyLine(
            weg, color=farbe, weight=3, opacity=0.85,
            tooltip=(f"Fahrzeug {nummer}: {tour['strecke']} km, "
                     f"{tour['fracht']}/{tour['kapazitaet']} Einheiten")
        ).add_to(karte)

    # Besuchsreihenfolge je Kunde - dieselbe Information, die im SVG als
    # Linienzug steckt, hier als Zahl am Punkt.
    reihenfolge = {}
    for nummer, tour in enumerate(touren, start=1):
        for stelle, ort in enumerate(tour["stationen"][1:-1], start=1):
            reihenfolge[ort] = (nummer, stelle)

    for ort in range(1, len(koordinaten)):
        fahrzeug, stelle = reihenfolge.get(ort, (0, 0))
        fenster = ZEITFENSTER[ort]
        folium.CircleMarker(
            [float(breite[ort]), float(laenge[ort])],
            radius=4 + 1.4 * BEDARFE[ort],
            color=TOURFARBEN[fahrzeug - 1] if fahrzeug else "#334155",
            fill=True, fill_opacity=0.75, weight=2,
            tooltip=(f"K{ort} — Bedarf {BEDARFE[ort]}, "
                     f"Zeitfenster {fenster[0]}–{fenster[1]} min<br>"
                     f"Fahrzeug {fahrzeug}, Stopp {stelle}"),
        ).add_to(karte)

    folium.CircleMarker(
        [float(breite[0]), float(laenge[0])], radius=9, color="#be123c",
        fill=True, fill_opacity=0.9, weight=3,
        tooltip=f"Depot — Schicht 0–{SCHICHTLAENGE} min",
    ).add_to(karte)

    os.makedirs(KARTEN, exist_ok=True)
    pfad = os.path.join(KARTEN, f"{KARTENNAME}.html")
    with open(pfad, "w", encoding="utf-8") as datei:
        datei.write(normalisiere(karte.get_root().render()))
    print(f"geschrieben: {pfad} "
          f"(Projektionsabweichung hoechstens {abweichung:.4f} %)")


if __name__ == "__main__":
    koordinaten, distanz = erzeuge_daten()
    touren = loese(distanz)
    gesamt = sum(t["strecke"] for t in touren)
    for nummer, tour in enumerate(touren, start=1):
        stationen = " -> ".join("Depot" if s == 0 else f"K{s}"
                                for s in tour["stationen"])
        print(f"Fahrzeug {nummer}: {tour['strecke']:3} km, "
              f"{tour['fracht']}/{tour['kapazitaet']} | {stationen}")
    print(f"Gesamtstrecke: {gesamt} km")
    zeichne(koordinaten, touren)
    breite, laenge = projiziere(koordinaten)
    zeichne_karte(koordinaten, touren,
                  pruefe_projektion(koordinaten, breite, laenge))