operations_research/bilder_04/erzeuge_pareto_front.py

#!/usr/bin/env python3

# erzeuge_pareto_front.py
"""
Erzeugt das Pareto-Front-Diagramm zum Mehrziel-Kapitel in zwei Fassungen:

    bilder_04/kap_mehrziel_pareto.svg          statisch, fuer das PDF
    bilder_04/plotly/kap_mehrziel_pareto.html  interaktiv, fuer die Website

Die Instanz ist identisch mit Mehrziel_Pareto.py, und gerechnet wird hier
erneut - so koennen Diagramm und Buchtext nicht auseinanderlaufen.

Das Diagramm zeigt genau die Aussage des Kapitels: Die untere konvexe Huelle
(gestrichelt) ist alles, was eine gewichtete Summe je erreichen kann. Die
Punkte darueber sind Pareto-optimal und trotzdem fuer kein Gewicht zu haben.

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

Benoetigt: numpy, scipy, matplotlib, plotly
"""

from __future__ import annotations

import os

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import plotly.graph_objects as go
from scipy.optimize import linprog

# 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")
PLOTLY_VERZ = os.path.join(BILDER, "plotly")

INDIGO, CYAN, AMBER = "#4338ca", "#0891b2", "#b45309"
BAHN_TRASSEN = 5


# --- Instanz und Modell: wortgleich zu Mehrziel_Pareto.py -------------------

def erzeuge_sendungen(anzahl: int = 12, saat: int = 5):
    rng = np.random.default_rng(saat)
    kosten = np.zeros((anzahl, 3))
    co2 = np.zeros((anzahl, 3))
    for i in range(anzahl):
        grund_kosten = rng.integers(600, 2400)
        grund_co2 = rng.integers(400, 1800)
        kosten[i] = [grund_kosten,
                     grund_kosten * rng.uniform(0.55, 0.85),
                     grund_kosten * rng.uniform(0.70, 1.00)]
        co2[i] = [grund_co2,
                  grund_co2 * rng.uniform(0.15, 0.35),
                  grund_co2 * rng.uniform(0.40, 0.70)]
    return np.round(kosten).astype(int), np.round(co2).astype(int)


KOSTEN, CO2 = erzeuge_sendungen()
N = len(KOSTEN)


def plane(ziel, co2_grenze=None):
    gleichungen = np.zeros((N, N * 3))
    for i in range(N):
        gleichungen[i, i * 3:(i + 1) * 3] = 1.0
    ungleichungen = [[1.0 if j == 1 else 0.0 for _ in range(N) for j in range(3)]]
    grenzen = [float(BAHN_TRASSEN)]
    if co2_grenze is not None:
        ungleichungen.append(CO2.reshape(-1).astype(float))
        grenzen.append(float(co2_grenze))
    ergebnis = linprog(ziel.reshape(-1).astype(float),
                       A_ub=ungleichungen, b_ub=grenzen,
                       A_eq=gleichungen, b_eq=np.ones(N),
                       bounds=(0, 1), integrality=1, method="highs")
    if not ergebnis.success:
        return None
    plan = np.round(ergebnis.x).astype(int)
    return int(KOSTEN.reshape(-1) @ plan), int(CO2.reshape(-1) @ plan)


def pareto_front():
    front = []
    grenze = plane(KOSTEN)[1]
    while True:
        ergebnis = plane(KOSTEN, co2_grenze=grenze)
        if ergebnis is None:
            return front
        front.append(ergebnis)
        grenze = ergebnis[1] - 1


def untere_huelle(punkte):
    huelle = []
    for punkt in sorted(punkte):
        while len(huelle) >= 2:
            (x1, y1), (x2, y2) = huelle[-2], huelle[-1]
            if (x2 - x1) * (punkt[1] - y1) - (y2 - y1) * (punkt[0] - x1) <= 0:
                huelle.pop()
            else:
                break
        huelle.append(punkt)
    return huelle


def tiefste_delle(front, huelle):
    """Der unerreichbare Punkt mit dem groessten Abstand zur Huelle.

    Liefert (Punkt, linke Huellenecke, rechte Huellenecke, Abstand). Genau an
    diesem Punkt ist die Luecke am besten zu sehen, und genau dort wird sie
    beschriftet - gesucht wird sie, nicht ausgewaehlt.

    Gemessen wird SENKRECHT zur Huellenkante, nicht senkrecht nach oben. Der
    senkrechte Abstand ist der Abstand eines Punktes zu einer Geraden; der
    Abstand in y-Richtung haengt dagegen davon ab, welche Zielgroesse man auf
    die Hochachse legt - bei vertauschten Achsen faende er eine andere Delle."""
    unerreichbar = [p for p in front if p not in huelle]
    if not unerreichbar:
        return None
    bestes = None
    for punkt in unerreichbar:
        for links, rechts in zip(huelle, huelle[1:]):
            if not links[0] <= punkt[0] <= rechts[0]:
                continue
            if rechts[0] == links[0]:
                continue
            steigung = (rechts[1] - links[1]) / (rechts[0] - links[0])
            auf_huelle = links[1] + steigung * (punkt[0] - links[0])
            abstand = (punkt[1] - auf_huelle) / np.hypot(1.0, steigung)
            if bestes is None or abstand > bestes[3]:
                bestes = (punkt, links, rechts, abstand)
    return bestes


def zeichne_luecke(achse, front, huelle) -> None:
    """Warum die gewichtete Summe die Delle nicht erreicht - und was hilft.

    Gezeichnet werden zwei Geraden derselben Steigung: die Tangente an die
    Huelle und die parallele Gerade durch den unerreichbaren Punkt. Eine
    gewichtete Summe minimiert genau diese Steigung; sie waehlt also immer die
    untere der beiden Geraden und damit nie den Punkt darueber. Der
    waagerechte Strich zeigt die Alternative: eine CO2-Schranke, unter der
    dieser Punkt der billigste ist - so ist die Front oben ueberhaupt
    entstanden."""
    delle = tiefste_delle(front, huelle)
    if delle is None:
        return
    punkt, links, rechts, _abstand = delle
    steigung = ((rechts[1] - links[1]) / (rechts[0] - links[0])
                if rechts[0] != links[0] else 0.0)

    x = np.array(achse.get_xlim())
    achse.plot(x, links[1] + steigung * (x - links[0]), "-", color=INDIGO,
               linewidth=1.1, alpha=0.75, zorder=1)
    achse.plot(x, punkt[1] + steigung * (x - punkt[0]), ":", color=AMBER,
               linewidth=1.6, zorder=2)
    achse.annotate("dieselbe Steigung, eine Stufe höher —\n"
                   "jede gewichtete Summe nimmt die untere",
                   xy=(punkt[0], punkt[1]), xytext=(12, 26),
                   textcoords="offset points", fontsize=7.5, color=AMBER,
                   zorder=6)

    achse.axhline(punkt[1], color="#059669", linestyle="--", linewidth=1.2,
                  alpha=0.9, zorder=1)
    achse.annotate(f"ε-Constraint: CO₂ ≤ {punkt[1]:,} kg".replace(",", ".")
                   + "\nfindet genau diesen Punkt",
                   xy=(x[0], punkt[1]), xytext=(6, -22),
                   textcoords="offset points", fontsize=7.5, color="#059669",
                   zorder=6)


# --- Statische Fassung ------------------------------------------------------

def zeichne_statisch(front, huelle) -> None:
    unerreichbar = [p for p in front if p not in huelle]
    figur, achse = plt.subplots(figsize=(7.2, 4.4))

    hx, hy = zip(*huelle)
    achse.plot(hx, hy, linestyle="--", color=INDIGO, linewidth=1.4, zorder=1,
               label="untere konvexe Hülle")
    achse.scatter(hx, hy, s=70, color=INDIGO, zorder=3,
                  label=f"von Gewichten erreichbar ({len(huelle)})")
    if unerreichbar:
        ux, uy = zip(*unerreichbar)
        achse.scatter(ux, uy, s=90, facecolors="none", edgecolors=AMBER,
                      linewidths=2.0, zorder=4,
                      label=f"Pareto-optimal, aber unerreichbar ({len(unerreichbar)})")

    for x, y in front:
        achse.annotate(f"{y:,} kg".replace(",", "."), (x, y),
                       textcoords="offset points", xytext=(6, 6), fontsize=7,
                       color="#475569")

    achse.set_xlabel("Transportkosten (EUR)")
    achse.set_ylabel("CO₂-Ausstoß (kg)")
    achse.set_title("Pareto-Front: Kosten gegen CO₂", fontsize=11)
    # Erst jetzt, damit die Geraden die endgueltigen Achsengrenzen benutzen.
    achse.autoscale(enable=False)
    zeichne_luecke(achse, front, huelle)
    achse.grid(linestyle=":", alpha=0.5)
    for rand in ("top", "right"):
        achse.spines[rand].set_visible(False)
    achse.legend(fontsize=8, loc="upper right", frameon=False)
    figur.tight_layout()

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


# --- Interaktive Fassung ----------------------------------------------------

def zeichne_interaktiv(front, huelle) -> None:
    unerreichbar = [p for p in front if p not in huelle]
    guenstigst = front[0]
    figur = go.Figure()

    hx, hy = zip(*huelle)
    figur.add_trace(go.Scatter(
        x=hx, y=hy, mode="lines", line=dict(color=INDIGO, dash="dash", width=1.5),
        name="untere konvexe Hülle", hoverinfo="skip"))

    def beschriftung(punkt):
        aufpreis = punkt[0] - guenstigst[0]
        ersparnis = guenstigst[1] - punkt[1]
        preis = f"{aufpreis / ersparnis:.2f} EUR je kg" if ersparnis else "Kostenminimum"
        return (f"<b>{punkt[0]:,} EUR</b><br>{punkt[1]:,} kg CO₂<br>"
                f"Aufpreis {aufpreis:,} EUR<br>"
                f"Ersparnis {ersparnis:,} kg<br>{preis}<extra></extra>")

    figur.add_trace(go.Scatter(
        x=hx, y=hy, mode="markers", marker=dict(size=11, color=INDIGO),
        name=f"von Gewichten erreichbar ({len(huelle)})",
        hovertemplate=[beschriftung(p) for p in huelle]))

    if unerreichbar:
        ux, uy = zip(*unerreichbar)
        figur.add_trace(go.Scatter(
            x=ux, y=uy, mode="markers",
            marker=dict(size=13, color="rgba(0,0,0,0)",
                        line=dict(color=AMBER, width=2.5)),
            name=f"Pareto-optimal, aber unerreichbar ({len(unerreichbar)})",
            hovertemplate=[beschriftung(p) for p in unerreichbar]))

    figur.update_layout(
        title="Pareto-Front: Kosten gegen CO₂",
        xaxis=dict(title="Transportkosten (EUR)"),
        yaxis=dict(title="CO₂-Ausstoß (kg)"),
        height=430, margin=dict(l=70, r=20, t=55, b=50),
        paper_bgcolor="rgba(0,0,0,0)", plot_bgcolor="rgba(0,0,0,0)",
        legend=dict(orientation="h", y=-0.2),
        hovermode="closest")

    os.makedirs(PLOTLY_VERZ, exist_ok=True)
    pfad = os.path.join(PLOTLY_VERZ, "kap_mehrziel_pareto.html")
    with open(pfad, "w", encoding="utf-8") as datei:
        datei.write(figur.to_html(full_html=False, include_plotlyjs=False,
                                  div_id="kap-mehrziel-pareto"))
    print(f"geschrieben: {pfad}")


if __name__ == "__main__":
    front = pareto_front()
    huelle = untere_huelle(front)
    print(f"Pareto-Punkte: {len(front)}, davon gestuetzt: {len(huelle)}, "
          f"unerreichbar: {len(front) - len(huelle)}")

    # Die beiden Behauptungen des Luecken-Bildes nachrechnen, statt sie
    # hinzuschreiben.
    delle = tiefste_delle(front, huelle)
    if delle is not None:
        punkt, links, rechts, abstand = delle
        steigung = (rechts[1] - links[1]) / (rechts[0] - links[0])

        def auf_tangente(x):
            return links[1] + steigung * (x - links[0])

        # Eine Stuetzgerade der UNTEREN Huelle liegt unter allen Punkten und
        # beruehrt die beiden Ecken der Kante. Beides zusammen macht sie zur
        # Tangente; nur eines davon reicht nicht.
        for ecke in huelle:
            if ecke[1] < auf_tangente(ecke[0]) - 1e-6:
                raise SystemExit(
                    f"Pareto-Luecke: Die Huellenecke {ecke} liegt UNTER der "
                    f"gezeichneten Geraden - dann stuetzt sie die Huelle "
                    f"nicht.")
        for ecke in (links, rechts):
            if abs(ecke[1] - auf_tangente(ecke[0])) > 1e-6:
                raise SystemExit(
                    f"Pareto-Luecke: Die Gerade beruehrt die Huellenecke "
                    f"{ecke} nicht - dann ist sie keine Tangente an diese "
                    f"Kante.")
        if punkt[1] <= auf_tangente(punkt[0]) + 1e-6:
            raise SystemExit(
                f"Pareto-Luecke: {punkt} liegt nicht ueber der Tangente - "
                f"dann waere der Punkt von einer gewichteten Summe erreichbar.")
        gefunden = plane(KOSTEN, co2_grenze=punkt[1])
        if gefunden != punkt:
            raise SystemExit(
                f"Pareto-Luecke: Die Schranke CO2 <= {punkt[1]} liefert "
                f"{gefunden}, nicht {punkt}. Die Bildunterschrift zum "
                f"ε-Constraint waere falsch.")
        # Die Zahl getrennt formatieren: Ein .replace(",", ".") auf den ganzen
        # f-String erwischt auch das Komma im Tupel.
        abstand_text = f"{abstand:.1f}".replace(".", ",")
        print(f"  tiefste Delle: {punkt}, Abstand zur Hülle "
              f"{abstand_text} (senkrecht)")
        print(f"  Tangente stützt die Hülle, der Punkt liegt darüber")
        print(f"  ε-Constraint CO₂ ≤ {punkt[1]} liefert genau diesen Punkt")

    zeichne_statisch(front, huelle)
    zeichne_interaktiv(front, huelle)