#!/usr/bin/env python3 # Mehrziel_Pareto.py """ Kapitel Mehrziel: Kosten gegen CO2 - und warum Gewichte nicht genuegen. Eine Spedition vergibt zwoelf Sendungen an drei Verkehrstraeger: LKW (schnell, teuer, schmutzig), Bahn (billig und sauber, aber nur fuenf Trassen frei) und Kombinierten Verkehr (dazwischen). Zwei Ziele stehen gegeneinander: Transportkosten und CO2-Ausstoss. Der uebliche Reflex ist, beide Ziele zu einem zusammenzuruehren: minimiere Kosten + w * CO2 Das ist bequem, liefert zulaessige Loesungen - und ist unvollstaendig. Bei ganzzahligen Entscheidungen gibt es Kompromisse, die auf diesem Weg GRUNDSAETZLICH nicht erreichbar sind, egal welches w man waehlt. Nicht "schwer zu finden", sondern beweisbar unerreichbar. Das Programm zeigt in vier Teilen: 1. Die beiden Extreme - was jedes Ziel allein kostet. 2. Die lineare Skalarisierung ueber ein feines Gewichtsraster. 3. Die vollstaendige Pareto-Front ueber das eps-Constraint-Verfahren. 4. Den Nachweis, dass die Luecke keine Frage des Rasters ist, sondern Geometrie: Die fehlenden Punkte liegen strikt oberhalb der konvexen Huelle und koennen deshalb von keiner Geraden gestuetzt werden. Benoetigt: numpy, scipy """ from __future__ import annotations import numpy as np from scipy.optimize import linprog TRAEGER = ["LKW", "Bahn", "Kombiniert"] BAHN_TRASSEN = 5 # so viele Sendungen passen hoechstens auf die Bahn SAAT = 5 def erzeuge_sendungen(anzahl: int = 12, saat: int = SAAT): """Kosten und CO2 je Sendung und Verkehrstraeger. Die Bahn ist immer billiger und sauberer als der LKW - der Zielkonflikt entsteht nicht zwischen den Traegern, sondern durch die KNAPPHEIT der Trassen. Genau so sieht es in der Praxis aus. """ 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: np.ndarray, co2_grenze: float | None = None, kosten_grenze: float | None = None): """Jede Sendung genau einem Traeger zuordnen; Trassen sind knapp. 'ziel' ist die zu minimierende Matrix (Kosten, CO2 oder eine Mischung). Die beiden Grenzen sind das Werkzeug fuer eps-Constraint und lexikografische Optimierung - sie machen aus einem Ziel eine Schranke. """ gleichungen = np.zeros((N, N * 3)) for i in range(N): gleichungen[i, i * 3:(i + 1) * 3] = 1.0 # genau ein Traeger 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)) if kosten_grenze is not None: ungleichungen.append(KOSTEN.reshape(-1).astype(float)) grenzen.append(float(kosten_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), plan) def pareto_front(): """Vollstaendige Front ueber das eps-Constraint-Verfahren. Der Ablauf ist der eigentliche Inhalt dieser Funktion: Erst das Kostenminimum bestimmen (der eine Rand der Front), dann die CO2-Schranke schrittweise um genau ein Kilogramm unter den zuletzt erreichten Wert druecken. Jeder Lauf liefert den naechsten Punkt - und wenn keiner mehr zulaessig ist, ist die Front vollstaendig. Das braucht so viele Solveraufrufe, wie die Front Punkte hat. Ein Raster ueber alle moeglichen CO2-Werte braeuchte hier fast tausend. """ front = [] start = plane(KOSTEN) grenze = start[1] while True: ergebnis = plane(KOSTEN, co2_grenze=grenze) if ergebnis is None: break front.append((ergebnis[0], ergebnis[1])) grenze = ergebnis[1] - 1 return front def skalarisierung(gewichte): """Was 'Kosten + w * CO2' fuer viele w hergibt - als Menge von Punkten.""" gefunden = {} for w in gewichte: ergebnis = plane(KOSTEN + w * CO2) if ergebnis: gefunden.setdefault((ergebnis[0], ergebnis[1]), []).append(w) return gefunden def untere_huelle(punkte): """Untere linke konvexe Huelle - genau die Punkte, die eine Gerade stuetzt. Der Zusammenhang, um den es geht: 'Kosten + w*CO2 minimieren' heisst geometrisch, eine Gerade der Steigung -1/w von links unten an die Punktwolke zu schieben. Sie beruehrt immer einen Eckpunkt der konvexen Huelle. Punkte, die oberhalb liegen, werden nie beruehrt - fuer kein w. """ 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 if __name__ == "__main__": print("=" * 80) print(" KOSTEN GEGEN CO2 - UND WARUM GEWICHTE NICHT GENUEGEN") print("=" * 80) print(f"{N} Sendungen, {len(TRAEGER)} Verkehrstraeger, " f"{BAHN_TRASSEN} freie Bahntrassen.\n") # --- 1. Die beiden Extreme ------------------------------------------- guenstigst = plane(KOSTEN) saubersten = plane(CO2) print("1. Was jedes Ziel allein ergibt\n") print(f" {'':<22} {'Kosten':>10} {'CO2':>10}") print(" " + "-" * 44) print(f" {'nur Kosten minimal':<22} {guenstigst[0]:>10,} {guenstigst[1]:>9,} kg") print(f" {'nur CO2 minimal':<22} {saubersten[0]:>10,} {saubersten[1]:>9,} kg") print(f"\n Der Zielkonflikt ist echt, aber klein: " f"{saubersten[0] - guenstigst[0]:,} EUR mehr") print(f" ({(saubersten[0] - guenstigst[0]) / guenstigst[0] * 100:.1f} %) sparen " f"{guenstigst[1] - saubersten[1]:,} kg CO2 " f"({(guenstigst[1] - saubersten[1]) / guenstigst[1] * 100:.1f} %).") print(" Genau solche Zahlen will die Geschaeftsfuehrung sehen - nicht ein") print(" Gewicht, das niemand interpretieren kann.") # --- 2. Die lineare Skalarisierung ----------------------------------- gewichte = np.concatenate([np.linspace(0.0, 3.0, 1201), np.geomspace(3.0, 1000.0, 200)]) gefunden = skalarisierung(gewichte) print("\n" + "-" * 80) print(f"2. Lineare Skalarisierung: 'Kosten + w * CO2' fuer " f"{len(gewichte):,} Gewichte\n") print(f" {'Kosten':>10} {'CO2':>10} {'gefunden bei w':>16}") print(" " + "-" * 46) for (k, c), ws in sorted(gefunden.items()): print(f" {k:>10,} {c:>9,} kg {min(ws):>7.3f} bis {max(ws):>7.3f}") print(f"\n {len(gefunden)} verschiedene Plaene - fuer {len(gewichte):,} Gewichte.") print(" Das Gewicht ist also gar keine Feineinstellung: Weite Bereiche") print(" liefern dasselbe Ergebnis, und dazwischen springt es.") # --- 3. Die vollstaendige Front -------------------------------------- front = pareto_front() print("\n" + "-" * 80) print(f"3. Die vollstaendige Pareto-Front ueber eps-Constraint " f"({len(front)} Solverlaeufe)\n") print(f" {'Kosten':>10} {'CO2':>10} {'Aufpreis':>9} {'CO2-Ersparnis':>14} " f"{'EUR je kg':>10}") print(" " + "-" * 60) for k, c in front: auf = k - guenstigst[0] ersparnis = guenstigst[1] - c preis = auf / ersparnis if ersparnis else 0.0 print(f" {k:>10,} {c:>9,} kg {auf:>8,} {ersparnis:>13,} " f"{preis:>10.2f}") # --- 4. Der Nachweis -------------------------------------------------- huelle = untere_huelle(front) unerreichbar = [p for p in front if p not in huelle] print("\n" + "-" * 80) print("4. Was die Skalarisierung nicht findet\n") erreicht = [p for p in front if p in gefunden] print(f" Pareto-Punkte insgesamt: {len(front)}") print(f" davon von der Skalarisierung gefunden: {len(erreicht)}") print(f" nie gefunden: {len(front) - len(erreicht)}") print() print(" Diese Kompromisse sind fuer KEIN Gewicht erreichbar:\n") print(f" {'Kosten':>10} {'CO2':>10} {'Aufpreis':>9} {'CO2-Ersparnis':>14}") print(" " + "-" * 50) for k, c in unerreichbar: print(f" {k:>10,} {c:>9,} kg {k - guenstigst[0]:>8,} " f"{guenstigst[1] - c:>13,}") stimmt = sorted(unerreichbar) == sorted(p for p in front if p not in gefunden) print(f"\n Gegenprobe ueber die Geometrie: {'bestanden' if stimmt else 'ABWEICHUNG'}") print(" Genau die Punkte, die das Gewichtsraster verfehlt, liegen strikt") print(" oberhalb der unteren konvexen Huelle. Das ist kein Rasterproblem -") print(" eine Gerade, die von links unten an die Wolke geschoben wird,") print(" beruehrt immer einen Eckpunkt der Huelle und nie einen Punkt") print(" darueber. Ein feineres Raster aendert daran nichts.") # --- 5. Lexikografisch ------------------------------------------------- print("\n" + "-" * 80) print("5. Lexikografisch: erst Kosten, dann CO2 im Rahmen eines Budgets\n") print(f" {'Kostenbudget':>14} {'Kosten':>10} {'CO2':>10} {'gegenueber Minimum':>20}") print(" " + "-" * 58) for aufschlag in (0.00, 0.01, 0.02, 0.05, 0.10): budget = guenstigst[0] * (1 + aufschlag) ergebnis = plane(CO2, kosten_grenze=budget) if ergebnis is None: print(f" {aufschlag:>13.0%} unzulaessig") continue print(f" {aufschlag:>13.0%} {ergebnis[0]:>10,} {ergebnis[1]:>9,} kg " f"{guenstigst[1] - ergebnis[1]:>15,} kg weniger") print("\n Das ist die Form, die im Betrieb am ehesten trifft: Nicht 'wie") print(" wichtig ist CO2?', sondern 'wir geben zwei Prozent mehr aus - was") print(" bringt das?'. Die Frage kann ein Kaufmann beantworten.") print("\n" + "=" * 80) print(" WAS MAN DARAUS MITNIMMT") print("=" * 80) print("Ein Gewicht zu setzen heisst, die Entscheidung heimlich zu treffen -") print("und dabei einen Teil der Moeglichkeiten gar nicht erst zu sehen.") print() print("Die Pareto-Front ist die ehrlichere Antwort: Sie legt dem Betrieb") print("alle sinnvollen Kompromisse vor und ueberlaesst ihm die Wahl. Die") print("Spalte 'EUR je kg' macht sie entscheidbar - man vergleicht sie mit") print("dem CO2-Preis, den das Unternehmen ohnehin ansetzt.") print("=" * 80)