Kapitel 17: Supply-Chain und Energieeinsatz unter Unsicherheit

#!/usr/bin/env python3

# Kraftwerkseinsatz.py
"""
Kapitel Supply-Chain: Welche Bloecke laufen morgen? Und was, wenn kein Wind weht?

Die Kraftwerkseinsatzplanung (englisch Unit Commitment) ist die Aufgabe, an der
sich in der Energiewirtschaft alles entscheidet: Fuer jede Stunde des naechsten
Tages muss feststehen, welche Bloecke am Netz sind. Ein Kernblock braucht acht
Stunden Mindestlaufzeit und 40.000 Euro Anfahrkosten - wer ihn abschaltet, hat
ihn fuer den Rest des Tages verloren.

Das Besondere liegt in der Zeitstruktur, und es geht ueber das zweistufige
Modell aus dem Kapitel Unsicherheit hinaus:

    STUFE 1   Das AN/AUS je Block und Stunde. Binaer, und es steht am Vorabend
              fest - bevor irgendjemand weiss, wie viel Wind morgen weht.
    STUFE 2   Die Fahrweise: wie viel jeder laufende Block liefert. Das darf
              sich stundenweise an die Wirklichkeit anpassen.

Die erste Stufe ist also GANZZAHLIG und szenariouebergreifend gleich, die
zweite kontinuierlich und je Szenario verschieden. Genau diese Kombination
macht das Problem interessant - und sie ist der Grund, warum ein Plan, der auf
den Wind-Erwartungswert gerechnet wurde, in der Wirklichkeit teuer wird.

Das Programm zeigt drei Plaene auf denselben Daten:

  1. Deterministisch: gerechnet mit dem Wind-ERWARTUNGSWERT, danach gegen 40
     Szenarien ausgewertet.
  2. Zweistufig: der Commitment-Plan sieht alle 40 Szenarien.
  3. Mit Versorgungssicherheit: die Nichtdeckung in den schlechtesten Faellen
     wird begrenzt - und der Preis dafuer in Euro je vermiedener MWh
     ausgewiesen.

Benoetigt: numpy, ortools
"""

from __future__ import annotations

import time

import numpy as np
from ortools.linear_solver import pywraplp

# (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 = 24
SZENARIEN = 40
FLAUTENANTEIL = 0.15
ALPHA = 0.90                  # die schlechtesten 10 % der Szenarien
LASTABWURF = 3_000.0          # EUR je nicht gedeckter MWh ("value of lost load")
NIEDRIGER_ABWURFPREIS = 300.0  # zum Vergleich: ein Marktpreisdeckel
SAAT = 4

_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 erzeuge_windszenarien(anzahl: int = SZENARIEN, saat: int = SAAT):
    """Windeinspeisung je Szenario und Stunde.

    Zwei Zutaten, die zusammen den Unterschied machen: eine breite
    lognormale Streuung des Tagesniveaus - und in 15 % der Faelle eine
    DUNKELFLAUTE, in der praktisch gar kein Wind weht. Solche seltenen,
    extremen Faelle sind der Grund, warum der Erwartungswert als
    Planungsgrundlage nicht genuegt.
    """
    rng = np.random.default_rng(saat)
    tagesniveau = rng.lognormal(0, 0.40, (anzahl, 1))
    wind = np.clip(WIND_ERWARTUNG * tagesniveau
                   + rng.normal(0, 25, (anzahl, STUNDEN)), 0, 420)
    ist_flaute = rng.random(anzahl) < FLAUTENANTEIL
    wind[ist_flaute] *= 0.05
    return wind, ist_flaute


def plane(wind: np.ndarray, gewichte, cvar_grenze: float | None = None,
          abwurfpreis: float = None, zeitlimit: float = 300.0):
    """Commitment (Stufe 1) und Fahrweise je Szenario (Stufe 2) in einem Modell.

    'wind' hat die Form (Szenarien, Stunden). Mit einer einzigen Zeile
    Windprognose wird daraus die deterministische Planung.

    'cvar_grenze' begrenzt den CVaR der Nichtdeckung ueber die Szenarien -
    dieselbe Rockafellar-Uryasev-Konstruktion wie im Kapitel CVaR, nur dass
    hier keine Verluste in Euro, sondern Megawattstunden begrenzt werden.
    """
    abwurfpreis = LASTABWURF if abwurfpreis is None else abwurfpreis
    anzahl_szenarien = wind.shape[0]
    anzahl_bloecke = len(KRAFTWERKE)
    solver = pywraplp.Solver.CreateSolver("SCIP")
    solver.SetTimeLimit(int(zeitlimit * 1000))

    # --- Stufe 1: binaer, szenariouebergreifend gleich --------------------
    laeuft = [[solver.BoolVar(f"laeuft_{k}_{t}") for t in range(STUNDEN)]
              for k in range(anzahl_bloecke)]
    faehrt_an = [[solver.BoolVar(f"start_{k}_{t}") for t in range(STUNDEN)]
                 for k in range(anzahl_bloecke)]

    # --- Stufe 2: kontinuierlich, je Szenario -----------------------------
    leistung = [[[solver.NumVar(0, KRAFTWERKE[k][2], f"p_{j}_{k}_{t}")
                  for t in range(STUNDEN)] for k in range(anzahl_bloecke)]
                for j in range(anzahl_szenarien)]
    nichtdeckung = [[solver.NumVar(0, solver.infinity(), f"y_{j}_{t}")
                     for t in range(STUNDEN)] for j in range(anzahl_szenarien)]

    for k, (_, pmin, pmax, _, _, mindestlaufzeit) in enumerate(KRAFTWERKE):
        for t in range(STUNDEN):
            # Anfahren erkennen: aus im Vortakt, an im aktuellen
            vorher = laeuft[k][t - 1] if t > 0 else 0
            solver.Add(faehrt_an[k][t] >= laeuft[k][t] - vorher)
            # Mindestlaufzeit: wer anfaehrt, laeuft die naechsten Stunden weiter
            for spaeter in range(t, min(STUNDEN, t + mindestlaufzeit)):
                solver.Add(laeuft[k][spaeter] >= faehrt_an[k][t])
            # Ein laufender Block liefert zwischen Mindest- und Nennleistung,
            # ein stehender gar nichts. Das koppelt beide Stufen.
            for j in range(anzahl_szenarien):
                solver.Add(leistung[j][k][t] >= pmin * laeuft[k][t])
                solver.Add(leistung[j][k][t] <= pmax * laeuft[k][t])

    for j in range(anzahl_szenarien):
        for t in range(STUNDEN):
            solver.Add(sum(leistung[j][k][t] for k in range(anzahl_bloecke))
                       + float(wind[j, t]) + nichtdeckung[j][t] >= LAST[t])

    # --- Versorgungssicherheit als CVaR-Schranke --------------------------
    if cvar_grenze is not None:
        schwelle = solver.NumVar(-solver.infinity(), solver.infinity(), "schwelle")
        ueberschuss = [solver.NumVar(0, solver.infinity(), f"u_{j}")
                       for j in range(anzahl_szenarien)]
        for j in range(anzahl_szenarien):
            solver.Add(ueberschuss[j] >= sum(nichtdeckung[j][t]
                                             for t in range(STUNDEN)) - schwelle)
        solver.Add(schwelle + (1.0 / (anzahl_szenarien * (1 - ALPHA)))
                   * sum(ueberschuss) <= cvar_grenze)

    anfahrkosten = sum(KRAFTWERKE[k][4] * faehrt_an[k][t]
                       for k in range(anzahl_bloecke) for t in range(STUNDEN))
    betriebskosten = sum(
        gewichte[j] * sum(KRAFTWERKE[k][3] * leistung[j][k][t]
                          for k in range(anzahl_bloecke) for t in range(STUNDEN))
        for j in range(anzahl_szenarien))
    abwurfkosten = sum(gewichte[j] * abwurfpreis
                       * sum(nichtdeckung[j][t] for t in range(STUNDEN))
                       for j in range(anzahl_szenarien))
    solver.Minimize(anfahrkosten + betriebskosten + abwurfkosten)

    beginn = time.perf_counter()
    status = solver.Solve()
    dauer = time.perf_counter() - beginn
    if status not in (pywraplp.Solver.OPTIMAL, pywraplp.Solver.FEASIBLE):
        return None, dauer, False
    plan = np.array([[laeuft[k][t].solution_value() for t in range(STUNDEN)]
                     for k in range(anzahl_bloecke)]).round()
    return plan, dauer, status == pywraplp.Solver.OPTIMAL


def bewerte(plan: np.ndarray, wind: np.ndarray, abwurfpreis: float = None):
    """Commitment steht fest - nur noch die Fahrweise je Szenario optimieren.

    Das ist die ehrliche Bewertung eines Plans: Er wird der Wirklichkeit
    ausgesetzt und darf nur noch das anpassen, was sich am Tag selbst
    anpassen laesst.
    """
    abwurfpreis = LASTABWURF if abwurfpreis is None else abwurfpreis
    anzahl_bloecke = len(KRAFTWERKE)
    kosten, fehlmengen = [], []
    for j in range(wind.shape[0]):
        solver = pywraplp.Solver.CreateSolver("GLOP")
        leistung = [[solver.NumVar(0, KRAFTWERKE[k][2], f"p_{k}_{t}")
                     for t in range(STUNDEN)] for k in range(anzahl_bloecke)]
        fehlt = [solver.NumVar(0, solver.infinity(), f"y_{t}")
                 for t in range(STUNDEN)]
        for k, (_, pmin, pmax, _, _, _) in enumerate(KRAFTWERKE):
            for t in range(STUNDEN):
                solver.Add(leistung[k][t] >= pmin * plan[k, t])
                solver.Add(leistung[k][t] <= pmax * plan[k, t])
        for t in range(STUNDEN):
            solver.Add(sum(leistung[k][t] for k in range(anzahl_bloecke))
                       + float(wind[j, t]) + fehlt[t] >= LAST[t])
        solver.Minimize(
            sum(KRAFTWERKE[k][3] * leistung[k][t]
                for k in range(anzahl_bloecke) for t in range(STUNDEN))
            + sum(abwurfpreis * fehlt[t] for t in range(STUNDEN)))
        solver.Solve()
        kosten.append(solver.Objective().Value())
        fehlmengen.append(sum(f.solution_value() for f in fehlt))

    anfahrten = sum(KRAFTWERKE[k][4]
                    * max(0.0, plan[k, t] - (plan[k, t - 1] if t > 0 else 0.0))
                    for k in range(anzahl_bloecke) for t in range(STUNDEN))
    return np.array(kosten) + anfahrten, np.array(fehlmengen)


def cvar(werte: np.ndarray, alpha: float = ALPHA) -> float:
    """Mittelwert der schlechtesten (1-alpha) Faelle."""
    schwelle = np.quantile(werte, alpha)
    schlechteste = werte[werte >= schwelle]
    return float(schlechteste.mean()) if len(schlechteste) else 0.0


def zeige(name: str, plan, wind, abwurfpreis: float = None) -> dict:
    kosten, fehl = bewerte(plan, wind, abwurfpreis)
    kennzahlen = {"mittel": kosten.mean(), "max": kosten.max(),
                  "fehl_mittel": fehl.mean(), "fehl_cvar": cvar(fehl),
                  "betroffen": int((fehl > 0.01).sum()),
                  "blockstunden": int(plan.sum())}
    print(f"  {name:<32} {kennzahlen['mittel']:>10,.0f} {kennzahlen['max']:>11,.0f} "
          f"{kennzahlen['fehl_mittel']:>9.1f} {kennzahlen['fehl_cvar']:>10.1f} "
          f"{kennzahlen['betroffen']:>6}/{len(fehl)}")
    return kennzahlen


if __name__ == "__main__":
    wind, ist_flaute = erzeuge_windszenarien()

    print("=" * 88)
    print("  KRAFTWERKSEINSATZ: DER PLAN STEHT, BEVOR DER WIND WEHT")
    print("=" * 88)
    print(f"{len(KRAFTWERKE)} Bloecke, {STUNDEN} Stunden, {SZENARIEN} Windszenarien "
          f"(davon {int(ist_flaute.sum())} Dunkelflauten).")
    print(f"Last {LAST.min():.0f} bis {LAST.max():.0f} MW, "
          f"Wind im Erwartungswert {WIND_ERWARTUNG.min():.0f} bis "
          f"{WIND_ERWARTUNG.max():.0f} MW.")
    print(f"Nicht gedeckte Last kostet {LASTABWURF:,.0f} EUR je MWh.\n")

    print(f"  {'Planungsgrundlage':<32} {'Kosten':>10} {'Kosten':>11} "
          f"{'Fehlmenge':>9} {'Fehlmenge':>10} {'Szen. mit':>12}")
    print(f"  {'':<32} {'im Mittel':>10} {'schlimmst.':>11} "
          f"{'Mittel':>9} {'CVaR 90%':>10} {'Abwurf':>12}")
    print("  " + "-" * 86)

    # --- 1. Deterministisch ----------------------------------------------
    plan_det, dauer_det, _ = plane(WIND_ERWARTUNG.reshape(1, -1), [1.0])
    det = zeige("nur Wind-Erwartungswert", plan_det, wind)

    # --- 2. Zweistufig, risikoneutral ------------------------------------
    gleich = [1.0 / SZENARIEN] * SZENARIEN
    plan_sto, dauer_sto, optimal_sto = plane(wind, gleich)
    sto = zeige("alle 40 Szenarien", plan_sto, wind)

    # --- 3. Mit Versorgungssicherheit ------------------------------------
    plan_sicher, dauer_sicher, optimal_sicher = plane(wind, gleich, cvar_grenze=0.0)
    sicher = zeige("+ Nichtdeckung CVaR 90 % = 0", plan_sicher, wind)

    print(f"\n  Rechenzeiten: {dauer_det:.1f}s / {dauer_sto:.1f}s / {dauer_sicher:.1f}s "
          f"(alle beweisbar optimal: {optimal_sto and optimal_sicher})")

    # --- Was die Zeilen bedeuten -----------------------------------------
    print("\n" + "-" * 88)
    print("Was der Erwartungswert-Plan anrichtet\n")
    print(f"  Er ist auf dem Papier der billigste - gerechnet auf den mittleren")
    print(f"  Wind kostet er weniger als jeder andere. In der Wirklichkeit liegt")
    print(f"  er {det['mittel'] / sto['mittel'] - 1:+.0%} ueber dem zweistufigen Plan, und sein "
          f"schlimmster Tag")
    print(f"  kostet {det['max'] / sto['max']:.1f}-mal so viel.")
    print(f"\n  Der Grund steht in der letzten Spalte: In {det['betroffen']} von {SZENARIEN} "
          f"Szenarien muss")
    print(f"  Last abgeworfen werden. Was fehlt, sind nicht Kraftwerke, sondern")
    print(f"  ANGEFAHRENE Kraftwerke - und ein Block mit acht Stunden")
    print(f"  Mindestlaufzeit laesst sich um 18 Uhr nicht mehr herbeirufen.")

    print(f"\n  Der Unterschied im Plan ist klein:")
    for name, plan in [("Erwartungswert", plan_det), ("zweistufig", plan_sto),
                       ("mit Sicherheit", plan_sicher)]:
        laufend = plan.sum(axis=0).astype(int)
        print(f"    {name:<16} {''.join(str(x) for x in laufend)}  "
              f"({int(plan.sum())} Blockstunden)")
    print(f"\n  Der zweistufige Plan haelt {sto['blockstunden'] - det['blockstunden']} "
          f"Blockstunden mehr vor - das genuegt,")
    print(f"  um alle {det['betroffen']} Lastabwuerfe zu vermeiden.")

    # --- Der Preis der Sicherheit ----------------------------------------
    print("\n" + "-" * 88)
    print("Was Versorgungssicherheit kostet\n")
    print(f"  Bei {LASTABWURF:,.0f} EUR/MWh kostet die CVaR-Schranke NICHTS: Der")
    print("  risikoneutrale Plan haelt sie schon ein. Das ist kein Zufall - bei")
    print("  diesem Preis lohnt sich Vorhaltung bereits im Erwartungswert.")
    print("\n  Interessant wird die Schranke dort, wo der Schaden ZU NIEDRIG")
    print("  bepreist ist. Dieselbe Rechnung mit einem Marktpreisdeckel von")
    print(f"  {NIEDRIGER_ABWURFPREIS:,.0f} EUR/MWh statt {LASTABWURF:,.0f}:\n")

    print(f"  {'Planungsgrundlage':<32} {'Kosten':>10} {'Kosten':>11} "
          f"{'Fehlmenge':>9} {'Fehlmenge':>10} {'Szen. mit':>12}")
    print(f"  {'':<32} {'im Mittel':>10} {'schlimmst.':>11} "
          f"{'Mittel':>9} {'CVaR 90%':>10} {'Abwurf':>12}")
    print("  " + "-" * 86)
    plan_billig, _, _ = plane(wind, gleich, abwurfpreis=NIEDRIGER_ABWURFPREIS)
    billig = zeige("risikoneutral, billiger Abwurf", plan_billig, wind,
                   NIEDRIGER_ABWURFPREIS)
    plan_billig_sicher, _, _ = plane(wind, gleich, cvar_grenze=0.0,
                                     abwurfpreis=NIEDRIGER_ABWURFPREIS)
    billig_sicher = zeige("+ CVaR 90 % der Fehlmenge = 0", plan_billig_sicher, wind,
                          NIEDRIGER_ABWURFPREIS)

    aufpreis = billig_sicher["mittel"] - billig["mittel"]
    vermieden = billig["fehl_cvar"] - billig_sicher["fehl_cvar"]
    print(f"\n  Jetzt greift die Schranke: Der risikoneutrale Plan nimmt in")
    print(f"  {billig['betroffen']} von {SZENARIEN} Szenarien einen Lastabwurf in Kauf, weil er bei")
    print(f"  {NIEDRIGER_ABWURFPREIS:,.0f} EUR/MWh billiger ist als das Vorhalten eines Blocks.")
    print(f"\n  Aufpreis fuer die Sicherheit: {aufpreis:,.0f} EUR je Tag")
    if vermieden > 0.01:
        print(f"  Vermiedene Fehlmenge (CVaR 90 %): {vermieden:.1f} MWh")
        print(f"  -> {aufpreis / vermieden:,.0f} EUR je vermiedener MWh")
        print(f"\n  Diese Zahl ist die Entscheidungsgrundlage - genau wie die Spalte")
        print(f"  'EUR je kg' im Kapitel Mehrziel. Sie sagt der Aufsicht, was ihre")
        print(f"  Versorgungssicherheitsvorgabe tatsaechlich kostet, statt darueber")
        print(f"  zu streiten, wie wichtig sie ist.")

    print("\n" + "=" * 88)
    print("  WAS MAN DARAUS MITNIMMT")
    print("=" * 88)
    print("1. Die erste Stufe ist binaer und traege. Was am Vorabend nicht")
    print("   angefahren wurde, steht am naechsten Tag nicht zur Verfuegung -")
    print("   egal, wie hoch der Preis dann steigt.")
    print("2. Deshalb ist der Erwartungswert als Planungsgrundlage nicht nur")
    print("   ungenau, sondern SYSTEMATISCH zu knapp: Er plant fuer einen Tag,")
    print("   der so nie eintritt, und laesst keine Reserve fuer die Haelfte")
    print("   der Faelle, in denen es schlechter kommt.")
    print("3. Der Ausweg ist keine bessere Windprognose, sondern ein Modell,")
    print("   das die Szenarien SIEHT. Die Rechenzeit dafuer liegt hier im")
    print("   Sekundenbereich - der Aufwand ist die Modellierung, nicht der")
    print("   Solver.")
    print("4. Wo der Lastabwurf teuer genug bepreist ist, deckt schon die")
    print("   Minimierung der ERWARTETEN Kosten die Extremfaelle mit ab. Die")
    print("   Risikoschranke wird genau dann gebraucht, wenn der Schaden ZU")
    print("   NIEDRIG bepreist ist - was bei Versorgungssicherheit die Regel")
    print("   ist, weil ein Marktpreisdeckel nicht den volkswirtschaftlichen")
    print("   Schaden abbildet. Dann liefert sie den Preis der Vorgabe in Euro")
    print("   je vermiedener MWh, und darueber laesst sich verhandeln.")
    print("=" * 88)