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