320 lines
13 KiB
Python
320 lines
13 KiB
Python
#!/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)
|