operations_research/bilder_04/erzeuge_wirkung.py

#!/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

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",):
        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)