Phase 7.2: Wasserfalldiagramm - wofuer der Erwartungswert-Plan bezahlt

Setzt den umsetzbaren Teil von Paket 3 aus Verbesserungen_03.md um
(Vorher-Nachher-Wasserfall, Business Impact). Neuer Generator
bilder_04/erzeuge_wirkung.py; das Supply-Chain-Kapitel bekommt damit sein
erstes Diagramm ueberhaupt.

Der Anlass: Die Tabelle nannte nur die Summen - 874.870 EUR gegen
351.356 EUR, 149 % mehr. Wofuer genau, stand nirgends.

Der Generator rechnet die Instanz erneut, wie es die Konvention verlangt.
Die Probe fiel gut aus: Der unabhaengige Nachbau reproduziert beide
abgedruckten Summen exakt, dazu 178,7 MWh Fehlmenge und 28/40 Szenarien mit
Abwurf. Kein Auseinanderlaufen - der erste Generator dieser Reihe, der
nichts findet.

Was die Zerlegung zeigt und die Tabelle nicht: Die Anfahrkosten sind in
beiden Plaenen gleich (58.000 EUR). Der Erwartungswert-Plan spart also nicht
dort, wo man es vermuten wuerde. Er spart 12.725 EUR Brennstoff und bezahlt
dafuer 536.239 EUR Lastabwurf - das 42-fache. "149 % mehr" ist eine Zahl,
ueber die man streiten kann; "wir sparen 12.725 EUR und riskieren
536.239 EUR" ist eine, ueber die man entscheidet.

Die Trennung von Brennstoff- und Abwurfkosten ist der einzige Zusatz
gegenueber bewerte() im Buchprogramm; dort wird sie nicht gebraucht.

Ein eigener Fehler unterwegs: Ein .replace(",", ".") fuer deutsche
Tausenderpunkte lag auf dem ganzen Titelsatz und machte aus "Brennstoff,
Preis" ein "Brennstoff. Preis". Jetzt gibt es euro(), das nur auf Zahlen
angewandt wird.

SVG byteidentisch ueber zwei Laeufe. Diagramme jetzt 33, davon 19 mit
Generator (16 Skripte). PDF unveraendert 746 Seiten.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

bilder_04/erzeuge_wirkung.py

@@ -0,0 +1,253 @@
+#!/usr/bin/env python3
+
+# erzeuge_wirkung.py
+"""
+Erzeugt das Wasserfalldiagramm zum Kapitel Supply-Chain und Energieeinsatz:
+
+    bilder_04/kap_supplychain_wirkung.svg   statisch, fuer das PDF
+    bilder_04/kap_supplychain_wirkung.png   Rasterfassung
+
+Die Frage, die es beantwortet, ist die des Managements: WOFUER genau kostet
+der Erwartungswert-Plan 149 % mehr? Die Tabelle im Kapitel nennt nur die
+Summen (874.870 EUR gegen 351.356 EUR). Der Wasserfall zerlegt die
+Differenz - und die Zerlegung ist der eigentliche Befund:
+
+  * Die ANFAHRKOSTEN sind in beiden Plaenen gleich (58.000 EUR). Der
+    Erwartungswert-Plan spart also nicht dort, wo man es vermuten wuerde.
+  * Er spart 12.725 EUR Brennstoff, weil er weniger Blockstunden vorhaelt.
+  * Und zahlt dafuer 536.239 EUR Lastabwurf.
+
+Das ist das 42-fache dessen, was er einspart. Genau diese Gegenueberstellung
+laesst sich in eine Vorlage uebernehmen, eine Summendifferenz nicht.
+
+Verfahren und Instanz sind wortgleich zu Kraftwerkseinsatz.py, und gerechnet
+wird hier ERNEUT - so koennen Diagramm und Buchtext nicht auseinanderlaufen.
+Der Nachbau reproduziert die im Kapitel abgedruckten Summen exakt; das ist
+die Probe darauf, dass beide dasselbe Modell meinen. Zusaetzlich trennt er
+Brennstoff- und Abwurfkosten, was bewerte() im Buchprogramm nicht tut, weil
+es dort nicht gebraucht wird.
+
+Aufruf (aus dem Repository-Wurzelverzeichnis):
+    python3 bilder_04/erzeuge_wirkung.py
+
+Benoetigt: numpy, matplotlib, ortools
+"""
+
+from __future__ import annotations
+
+import os
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+from ortools.linear_solver import pywraplp
+
+# 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")
+
+INDIGO, CYAN, AMBER, GRUEN = "#4338ca", "#0891b2", "#b45309", "#059669"
+GRAU = "#64748b"
+
+# --- Instanz: wortgleich zu Kraftwerkseinsatz.py ----------------------------
+# (Name, Mindestleistung, Nennleistung, Grenzkosten, Anfahrkosten, Mindestlaufzeit)
+KRAFTWERKE = [
+    ("Kernblock",   200, 600,  22, 40_000, 8),
+    ("Braunkohle",  100, 400,  35, 18_000, 6),
+    ("Steinkohle",   80, 300,  48, 12_000, 4),
+    ("Gas GuD",      50, 250,  72,  6_000, 2),
+    ("Gasturbine",   20, 150, 105,   1_500, 1),
+]
+STUNDEN, SZENARIEN, FLAUTENANTEIL, SAAT = 24, 40, 0.15, 4
+LASTABWURF = 3_000.0          # EUR je nicht gedeckter MWh
+
+_t = np.arange(STUNDEN)
+LAST = 620 + 260 * np.sin((_t - 7) / 24 * 2 * np.pi) + 90 * np.sin((_t - 4) / 12 * 2 * np.pi)
+WIND_ERWARTUNG = np.maximum(0.0, 160 + 130 * np.sin((_t - 14) / 24 * 2 * np.pi))
+
+
+def euro(betrag: float, vorzeichen: bool = False) -> str:
+    """Deutsche Tausenderpunkte. Nur auf ZAHLEN anwenden, nie auf ganze Saetze."""
+    text = f"{betrag:+,.0f}" if vorzeichen else f"{betrag:,.0f}"
+    return text.replace(",", ".") + " €"
+
+
+def windszenarien() -> tuple[np.ndarray, np.ndarray]:
+    rng = np.random.default_rng(SAAT)
+    niveau = rng.lognormal(0, 0.40, (SZENARIEN, 1))
+    wind = np.clip(WIND_ERWARTUNG * niveau
+                   + rng.normal(0, 25, (SZENARIEN, STUNDEN)), 0, 420)
+    flaute = rng.random(SZENARIEN) < FLAUTENANTEIL
+    wind[flaute] *= 0.05
+    return wind, flaute
+
+
+def plane(wind: np.ndarray, gewichte: np.ndarray) -> np.ndarray:
+    """Stufe 1 (An/Aus, szenariouebergreifend) und Stufe 2 (Fahrweise je
+    Szenario) in einem Modell - wie plane() in Kraftwerkseinsatz.py."""
+    n_s, n_b = wind.shape[0], len(KRAFTWERKE)
+    s = pywraplp.Solver.CreateSolver("SCIP")
+    s.SetTimeLimit(300_000)
+    laeuft = [[s.BoolVar(f"l{k}_{t}") for t in range(STUNDEN)] for k in range(n_b)]
+    start = [[s.BoolVar(f"s{k}_{t}") for t in range(STUNDEN)] for k in range(n_b)]
+    p = [[[s.NumVar(0, KRAFTWERKE[k][2], f"p{j}_{k}_{t}") for t in range(STUNDEN)]
+          for k in range(n_b)] for j in range(n_s)]
+    y = [[s.NumVar(0, s.infinity(), f"y{j}_{t}") for t in range(STUNDEN)]
+         for j in range(n_s)]
+
+    for k, (_, pmin, pmax, _, _, mindestlaufzeit) in enumerate(KRAFTWERKE):
+        for t in range(STUNDEN):
+            vorher = laeuft[k][t - 1] if t > 0 else 0
+            s.Add(start[k][t] >= laeuft[k][t] - vorher)
+            for spaeter in range(t, min(STUNDEN, t + mindestlaufzeit)):
+                s.Add(laeuft[k][spaeter] >= start[k][t])
+            for j in range(n_s):
+                s.Add(p[j][k][t] >= pmin * laeuft[k][t])
+                s.Add(p[j][k][t] <= pmax * laeuft[k][t])
+    for j in range(n_s):
+        for t in range(STUNDEN):
+            s.Add(sum(p[j][k][t] for k in range(n_b))
+                  + float(wind[j, t]) + y[j][t] >= LAST[t])
+
+    s.Minimize(
+        sum(KRAFTWERKE[k][4] * start[k][t]
+            for k in range(n_b) for t in range(STUNDEN))
+        + sum(gewichte[j] * sum(KRAFTWERKE[k][3] * p[j][k][t]
+                                for k in range(n_b) for t in range(STUNDEN))
+              for j in range(n_s))
+        + sum(gewichte[j] * LASTABWURF * sum(y[j][t] for t in range(STUNDEN))
+              for j in range(n_s)))
+    s.Solve()
+    return np.array([[laeuft[k][t].solution_value() for t in range(STUNDEN)]
+                     for k in range(n_b)]).round()
+
+
+def zerlege(plan: np.ndarray, wind: np.ndarray) -> dict:
+    """Setzt einen fertigen Plan der Wirklichkeit aus und trennt dabei die
+    Kostenarten - das ist der Zusatz gegenueber bewerte() im Buchprogramm."""
+    n_b = len(KRAFTWERKE)
+    brennstoff, abwurf, fehlmenge = [], [], []
+    for j in range(wind.shape[0]):
+        s = pywraplp.Solver.CreateSolver("GLOP")
+        p = [[s.NumVar(0, KRAFTWERKE[k][2], f"p{k}_{t}") for t in range(STUNDEN)]
+             for k in range(n_b)]
+        y = [s.NumVar(0, s.infinity(), f"y{t}") for t in range(STUNDEN)]
+        for k, (_, pmin, pmax, _, _, _) in enumerate(KRAFTWERKE):
+            for t in range(STUNDEN):
+                s.Add(p[k][t] >= pmin * plan[k, t])
+                s.Add(p[k][t] <= pmax * plan[k, t])
+        for t in range(STUNDEN):
+            s.Add(sum(p[k][t] for k in range(n_b))
+                  + float(wind[j, t]) + y[t] >= LAST[t])
+        s.Minimize(sum(KRAFTWERKE[k][3] * p[k][t]
+                       for k in range(n_b) for t in range(STUNDEN))
+                   + sum(LASTABWURF * y[t] for t in range(STUNDEN)))
+        s.Solve()
+        brennstoff.append(sum(KRAFTWERKE[k][3] * p[k][t].solution_value()
+                              for k in range(n_b) for t in range(STUNDEN)))
+        fehlt = sum(v.solution_value() for v in y)
+        fehlmenge.append(fehlt)
+        abwurf.append(LASTABWURF * fehlt)
+
+    anfahrt = sum(KRAFTWERKE[k][4]
+                  * max(0.0, plan[k, t] - (plan[k, t - 1] if t > 0 else 0.0))
+                  for k in range(n_b) for t in range(STUNDEN))
+    fehlmenge = np.array(fehlmenge)
+    return {"anfahrt": float(anfahrt),
+            "brennstoff": float(np.mean(brennstoff)),
+            "abwurf": float(np.mean(abwurf)),
+            "fehlmenge": float(fehlmenge.mean()),
+            "szenarien_mit_abwurf": int((fehlmenge > 1e-6).sum()),
+            "gesamt": float(anfahrt + np.mean(brennstoff) + np.mean(abwurf))}
+
+
+def zeichne(zwei: dict, erw: dict) -> None:
+    """Wasserfall vom guten zum schlechten Plan: Was wird gespart, was kostet es."""
+    d_anfahrt = erw["anfahrt"] - zwei["anfahrt"]
+    d_brennstoff = erw["brennstoff"] - zwei["brennstoff"]
+    d_abwurf = erw["abwurf"] - zwei["abwurf"]
+
+    stufen = [("Zweistufiger\nPlan", zwei["gesamt"], None),
+              ("Brennstoff", d_brennstoff, "delta"),
+              ("Anfahren", d_anfahrt, "delta"),
+              ("Lastabwurf", d_abwurf, "delta"),
+              ("Erwartungswert-\nPlan", erw["gesamt"], None)]
+
+    figur, achse = plt.subplots(figsize=(10.2, 5.6))
+    unten = 0.0
+    for i, (name, wert, art) in enumerate(stufen):
+        if art is None:
+            achse.bar(i, wert, 0.62, color=INDIGO, zorder=3)
+            achse.text(i, wert + 18_000, euro(wert),
+                       ha="center", va="bottom", fontsize=11, fontweight="bold",
+                       color=INDIGO)
+            unten = wert
+            continue
+        farbe = AMBER if wert > 0 else GRUEN
+        start = unten if wert > 0 else unten + wert
+        achse.bar(i, abs(wert), 0.62, bottom=start, color=farbe, zorder=3)
+        # Nullbalken (hier: Anfahren) sichtbar machen, sonst fehlt die Aussage
+        if abs(wert) < 5_000:
+            achse.plot([i - 0.31, i + 0.31], [start, start], color=GRAU,
+                       linewidth=2.2, zorder=4)
+        achse.text(i, max(start + abs(wert), unten) + 18_000,
+                   euro(wert, vorzeichen=True),
+                   ha="center", va="bottom", fontsize=10.5, color=farbe,
+                   fontweight="bold")
+        unten = unten + wert
+        if i < len(stufen) - 1:
+            achse.plot([i - 0.31, i + 0.69], [unten, unten], color=GRAU,
+                       linewidth=0.9, linestyle=":", zorder=2)
+
+    achse.set_xticks(range(len(stufen)))
+    achse.set_xticklabels([n for n, _, _ in stufen], fontsize=10.5)
+    achse.set_ylabel("Erwartete Tageskosten")
+    achse.set_ylim(0, erw["gesamt"] * 1.18)
+    achse.yaxis.set_major_formatter(
+        plt.FuncFormatter(lambda v, _: f"{v/1000:,.0f}k €".replace(",", ".")))
+    achse.grid(axis="y", linestyle=":", alpha=0.45, zorder=0)
+    for rand in ("top", "right"):
+        achse.spines[rand].set_visible(False)
+
+    verhaeltnis = abs(d_abwurf / d_brennstoff) if d_brennstoff else float("inf")
+    # Zahlen einzeln setzen: Ein .replace(",", ".") auf den ganzen Satz wuerde
+    # auch das Satzkomma treffen.
+    achse.set_title(
+        "Was der Erwartungswert-Plan spart — und was er dafür bezahlt\n"
+        f"Ersparnis {euro(abs(d_brennstoff))} Brennstoff, Preis "
+        f"{euro(d_abwurf)} Lastabwurf: das {verhaeltnis:.0f}-fache",
+        fontsize=12.5, pad=14)
+    achse.text(0.5, -0.155,
+               f"Lastabwurf in {erw['szenarien_mit_abwurf']} von {SZENARIEN} "
+               f"Szenarien statt in {zwei['szenarien_mit_abwurf']}; "
+               f"Anfahrkosten in beiden Plänen gleich.",
+               transform=achse.transAxes, ha="center", fontsize=10, color=GRAU)
+
+    figur.tight_layout()
+    os.makedirs(BILDER, exist_ok=True)
+    for endung in ("svg", "png"):
+        pfad = os.path.join(BILDER, f"kap_supplychain_wirkung.{endung}")
+        figur.savefig(pfad, format=endung, dpi=160,
+                      metadata={"Date": None} if endung == "svg" else None)
+        print(f"geschrieben: {pfad}")
+    plt.close(figur)
+
+
+if __name__ == "__main__":
+    wind, flaute = windszenarien()
+    gleich = np.full(SZENARIEN, 1.0 / SZENARIEN)
+    print(f"{SZENARIEN} Windszenarien, davon {int(flaute.sum())} Dunkelflauten")
+
+    erw = zerlege(plane(WIND_ERWARTUNG.reshape(1, -1), np.array([1.0])), wind)
+    zwei = zerlege(plane(wind, gleich), wind)
+
+    for name, w in (("Erwartungswert", erw), ("zweistufig", zwei)):
+        print(f"  {name:<16} Anfahrt {w['anfahrt']:>9,.0f}  "
+              f"Brennstoff {w['brennstoff']:>9,.0f}  "
+              f"Abwurf {w['abwurf']:>9,.0f}  = {w['gesamt']:>9,.0f}")
+    print(f"  Aufschlag des Erwartungswert-Plans: "
+          f"{erw['gesamt'] / zwei['gesamt'] - 1:+.0%}")
+    zeichne(zwei, erw)