{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Kapitel 17: Supply-Chain und Energieeinsatz unter Unsicherheit\n", "\n", "Begleitnotebook zu *Optimierte Entscheidungsfindung mit Python*. Die Codezellen sind identisch mit den im Buch abgedruckten Programmen.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Einmalig ausfuehren: installiert alle im Buch verwendeten Pakete.\n", "# Lokal in einer virtuellen Umgebung genauso gueltig wie in Google Colab.\n", "%pip install --quiet ortools highspy cvxpy scipy numpy pandas polars \\\n", " scikit-learn matplotlib plotly pyomo linopy pymoo pydantic openpyxl" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Das Programm\n", "\n", "`Kraftwerkseinsatz.py`\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#!/usr/bin/env python3\n", "\n", "# Kraftwerkseinsatz.py\n", "\"\"\"\n", "Kapitel Supply-Chain: Welche Bloecke laufen morgen? Und was, wenn kein Wind weht?\n", "\n", "Die Kraftwerkseinsatzplanung (englisch Unit Commitment) ist die Aufgabe, an der\n", "sich in der Energiewirtschaft alles entscheidet: Fuer jede Stunde des naechsten\n", "Tages muss feststehen, welche Bloecke am Netz sind. Ein Kernblock braucht acht\n", "Stunden Mindestlaufzeit und 40.000 Euro Anfahrkosten - wer ihn abschaltet, hat\n", "ihn fuer den Rest des Tages verloren.\n", "\n", "Das Besondere liegt in der Zeitstruktur, und es geht ueber das zweistufige\n", "Modell aus dem Kapitel Unsicherheit hinaus:\n", "\n", " STUFE 1 Das AN/AUS je Block und Stunde. Binaer, und es steht am Vorabend\n", " fest - bevor irgendjemand weiss, wie viel Wind morgen weht.\n", " STUFE 2 Die Fahrweise: wie viel jeder laufende Block liefert. Das darf\n", " sich stundenweise an die Wirklichkeit anpassen.\n", "\n", "Die erste Stufe ist also GANZZAHLIG und szenariouebergreifend gleich, die\n", "zweite kontinuierlich und je Szenario verschieden. Genau diese Kombination\n", "macht das Problem interessant - und sie ist der Grund, warum ein Plan, der auf\n", "den Wind-Erwartungswert gerechnet wurde, in der Wirklichkeit teuer wird.\n", "\n", "Das Programm zeigt drei Plaene auf denselben Daten:\n", "\n", " 1. Deterministisch: gerechnet mit dem Wind-ERWARTUNGSWERT, danach gegen 40\n", " Szenarien ausgewertet.\n", " 2. Zweistufig: der Commitment-Plan sieht alle 40 Szenarien.\n", " 3. Mit Versorgungssicherheit: die Nichtdeckung in den schlechtesten Faellen\n", " wird begrenzt - und der Preis dafuer in Euro je vermiedener MWh\n", " ausgewiesen.\n", "\n", "Benoetigt: numpy, ortools\n", "\"\"\"\n", "\n", "from __future__ import annotations\n", "\n", "import time\n", "\n", "import numpy as np\n", "from ortools.linear_solver import pywraplp\n", "\n", "# (Name, Mindestleistung, Nennleistung, Grenzkosten, Anfahrkosten, Mindestlaufzeit)\n", "KRAFTWERKE = [\n", " (\"Kernblock\", 200, 600, 22, 40_000, 8),\n", " (\"Braunkohle\", 100, 400, 35, 18_000, 6),\n", " (\"Steinkohle\", 80, 300, 48, 12_000, 4),\n", " (\"Gas GuD\", 50, 250, 72, 6_000, 2),\n", " (\"Gasturbine\", 20, 150, 105, 1_500, 1),\n", "]\n", "STUNDEN = 24\n", "SZENARIEN = 40\n", "FLAUTENANTEIL = 0.15\n", "ALPHA = 0.90 # die schlechtesten 10 % der Szenarien\n", "LASTABWURF = 3_000.0 # EUR je nicht gedeckter MWh (\"value of lost load\")\n", "NIEDRIGER_ABWURFPREIS = 300.0 # zum Vergleich: ein Marktpreisdeckel\n", "SAAT = 4\n", "\n", "_t = np.arange(STUNDEN)\n", "LAST = 620 + 260 * np.sin((_t - 7) / 24 * 2 * np.pi) + 90 * np.sin((_t - 4) / 12 * 2 * np.pi)\n", "WIND_ERWARTUNG = np.maximum(0.0, 160 + 130 * np.sin((_t - 14) / 24 * 2 * np.pi))\n", "\n", "\n", "def erzeuge_windszenarien(anzahl: int = SZENARIEN, saat: int = SAAT):\n", " \"\"\"Windeinspeisung je Szenario und Stunde.\n", "\n", " Zwei Zutaten, die zusammen den Unterschied machen: eine breite\n", " lognormale Streuung des Tagesniveaus - und in 15 % der Faelle eine\n", " DUNKELFLAUTE, in der praktisch gar kein Wind weht. Solche seltenen,\n", " extremen Faelle sind der Grund, warum der Erwartungswert als\n", " Planungsgrundlage nicht genuegt.\n", " \"\"\"\n", " rng = np.random.default_rng(saat)\n", " tagesniveau = rng.lognormal(0, 0.40, (anzahl, 1))\n", " wind = np.clip(WIND_ERWARTUNG * tagesniveau\n", " + rng.normal(0, 25, (anzahl, STUNDEN)), 0, 420)\n", " ist_flaute = rng.random(anzahl) < FLAUTENANTEIL\n", " wind[ist_flaute] *= 0.05\n", " return wind, ist_flaute\n", "\n", "\n", "def plane(wind: np.ndarray, gewichte, cvar_grenze: float | None = None,\n", " abwurfpreis: float = None, zeitlimit: float = 300.0):\n", " \"\"\"Commitment (Stufe 1) und Fahrweise je Szenario (Stufe 2) in einem Modell.\n", "\n", " 'wind' hat die Form (Szenarien, Stunden). Mit einer einzigen Zeile\n", " Windprognose wird daraus die deterministische Planung.\n", "\n", " 'cvar_grenze' begrenzt den CVaR der Nichtdeckung ueber die Szenarien -\n", " dieselbe Rockafellar-Uryasev-Konstruktion wie im Kapitel CVaR, nur dass\n", " hier keine Verluste in Euro, sondern Megawattstunden begrenzt werden.\n", " \"\"\"\n", " abwurfpreis = LASTABWURF if abwurfpreis is None else abwurfpreis\n", " anzahl_szenarien = wind.shape[0]\n", " anzahl_bloecke = len(KRAFTWERKE)\n", " solver = pywraplp.Solver.CreateSolver(\"SCIP\")\n", " solver.SetTimeLimit(int(zeitlimit * 1000))\n", "\n", " # --- Stufe 1: binaer, szenariouebergreifend gleich --------------------\n", " laeuft = [[solver.BoolVar(f\"laeuft_{k}_{t}\") for t in range(STUNDEN)]\n", " for k in range(anzahl_bloecke)]\n", " faehrt_an = [[solver.BoolVar(f\"start_{k}_{t}\") for t in range(STUNDEN)]\n", " for k in range(anzahl_bloecke)]\n", "\n", " # --- Stufe 2: kontinuierlich, je Szenario -----------------------------\n", " leistung = [[[solver.NumVar(0, KRAFTWERKE[k][2], f\"p_{j}_{k}_{t}\")\n", " for t in range(STUNDEN)] for k in range(anzahl_bloecke)]\n", " for j in range(anzahl_szenarien)]\n", " nichtdeckung = [[solver.NumVar(0, solver.infinity(), f\"y_{j}_{t}\")\n", " for t in range(STUNDEN)] for j in range(anzahl_szenarien)]\n", "\n", " for k, (_, pmin, pmax, _, _, mindestlaufzeit) in enumerate(KRAFTWERKE):\n", " for t in range(STUNDEN):\n", " # Anfahren erkennen: aus im Vortakt, an im aktuellen\n", " vorher = laeuft[k][t - 1] if t > 0 else 0\n", " solver.Add(faehrt_an[k][t] >= laeuft[k][t] - vorher)\n", " # Mindestlaufzeit: wer anfaehrt, laeuft die naechsten Stunden weiter\n", " for spaeter in range(t, min(STUNDEN, t + mindestlaufzeit)):\n", " solver.Add(laeuft[k][spaeter] >= faehrt_an[k][t])\n", " # Ein laufender Block liefert zwischen Mindest- und Nennleistung,\n", " # ein stehender gar nichts. Das koppelt beide Stufen.\n", " for j in range(anzahl_szenarien):\n", " solver.Add(leistung[j][k][t] >= pmin * laeuft[k][t])\n", " solver.Add(leistung[j][k][t] <= pmax * laeuft[k][t])\n", "\n", " for j in range(anzahl_szenarien):\n", " for t in range(STUNDEN):\n", " solver.Add(sum(leistung[j][k][t] for k in range(anzahl_bloecke))\n", " + float(wind[j, t]) + nichtdeckung[j][t] >= LAST[t])\n", "\n", " # --- Versorgungssicherheit als CVaR-Schranke --------------------------\n", " if cvar_grenze is not None:\n", " schwelle = solver.NumVar(-solver.infinity(), solver.infinity(), \"schwelle\")\n", " ueberschuss = [solver.NumVar(0, solver.infinity(), f\"u_{j}\")\n", " for j in range(anzahl_szenarien)]\n", " for j in range(anzahl_szenarien):\n", " solver.Add(ueberschuss[j] >= sum(nichtdeckung[j][t]\n", " for t in range(STUNDEN)) - schwelle)\n", " solver.Add(schwelle + (1.0 / (anzahl_szenarien * (1 - ALPHA)))\n", " * sum(ueberschuss) <= cvar_grenze)\n", "\n", " anfahrkosten = sum(KRAFTWERKE[k][4] * faehrt_an[k][t]\n", " for k in range(anzahl_bloecke) for t in range(STUNDEN))\n", " betriebskosten = sum(\n", " gewichte[j] * sum(KRAFTWERKE[k][3] * leistung[j][k][t]\n", " for k in range(anzahl_bloecke) for t in range(STUNDEN))\n", " for j in range(anzahl_szenarien))\n", " abwurfkosten = sum(gewichte[j] * abwurfpreis\n", " * sum(nichtdeckung[j][t] for t in range(STUNDEN))\n", " for j in range(anzahl_szenarien))\n", " solver.Minimize(anfahrkosten + betriebskosten + abwurfkosten)\n", "\n", " beginn = time.perf_counter()\n", " status = solver.Solve()\n", " dauer = time.perf_counter() - beginn\n", " if status not in (pywraplp.Solver.OPTIMAL, pywraplp.Solver.FEASIBLE):\n", " return None, dauer, False\n", " plan = np.array([[laeuft[k][t].solution_value() for t in range(STUNDEN)]\n", " for k in range(anzahl_bloecke)]).round()\n", " return plan, dauer, status == pywraplp.Solver.OPTIMAL\n", "\n", "\n", "def bewerte(plan: np.ndarray, wind: np.ndarray, abwurfpreis: float = None):\n", " \"\"\"Commitment steht fest - nur noch die Fahrweise je Szenario optimieren.\n", "\n", " Das ist die ehrliche Bewertung eines Plans: Er wird der Wirklichkeit\n", " ausgesetzt und darf nur noch das anpassen, was sich am Tag selbst\n", " anpassen laesst.\n", " \"\"\"\n", " abwurfpreis = LASTABWURF if abwurfpreis is None else abwurfpreis\n", " anzahl_bloecke = len(KRAFTWERKE)\n", " kosten, fehlmengen = [], []\n", " for j in range(wind.shape[0]):\n", " solver = pywraplp.Solver.CreateSolver(\"GLOP\")\n", " leistung = [[solver.NumVar(0, KRAFTWERKE[k][2], f\"p_{k}_{t}\")\n", " for t in range(STUNDEN)] for k in range(anzahl_bloecke)]\n", " fehlt = [solver.NumVar(0, solver.infinity(), f\"y_{t}\")\n", " for t in range(STUNDEN)]\n", " for k, (_, pmin, pmax, _, _, _) in enumerate(KRAFTWERKE):\n", " for t in range(STUNDEN):\n", " solver.Add(leistung[k][t] >= pmin * plan[k, t])\n", " solver.Add(leistung[k][t] <= pmax * plan[k, t])\n", " for t in range(STUNDEN):\n", " solver.Add(sum(leistung[k][t] for k in range(anzahl_bloecke))\n", " + float(wind[j, t]) + fehlt[t] >= LAST[t])\n", " solver.Minimize(\n", " sum(KRAFTWERKE[k][3] * leistung[k][t]\n", " for k in range(anzahl_bloecke) for t in range(STUNDEN))\n", " + sum(abwurfpreis * fehlt[t] for t in range(STUNDEN)))\n", " solver.Solve()\n", " kosten.append(solver.Objective().Value())\n", " fehlmengen.append(sum(f.solution_value() for f in fehlt))\n", "\n", " anfahrten = sum(KRAFTWERKE[k][4]\n", " * max(0.0, plan[k, t] - (plan[k, t - 1] if t > 0 else 0.0))\n", " for k in range(anzahl_bloecke) for t in range(STUNDEN))\n", " return np.array(kosten) + anfahrten, np.array(fehlmengen)\n", "\n", "\n", "def cvar(werte: np.ndarray, alpha: float = ALPHA) -> float:\n", " \"\"\"Mittelwert der schlechtesten (1-alpha) Faelle.\"\"\"\n", " schwelle = np.quantile(werte, alpha)\n", " schlechteste = werte[werte >= schwelle]\n", " return float(schlechteste.mean()) if len(schlechteste) else 0.0\n", "\n", "\n", "def zeige(name: str, plan, wind, abwurfpreis: float = None) -> dict:\n", " kosten, fehl = bewerte(plan, wind, abwurfpreis)\n", " kennzahlen = {\"mittel\": kosten.mean(), \"max\": kosten.max(),\n", " \"fehl_mittel\": fehl.mean(), \"fehl_cvar\": cvar(fehl),\n", " \"betroffen\": int((fehl > 0.01).sum()),\n", " \"blockstunden\": int(plan.sum())}\n", " print(f\" {name:<32} {kennzahlen['mittel']:>10,.0f} {kennzahlen['max']:>11,.0f} \"\n", " f\"{kennzahlen['fehl_mittel']:>9.1f} {kennzahlen['fehl_cvar']:>10.1f} \"\n", " f\"{kennzahlen['betroffen']:>6}/{len(fehl)}\")\n", " return kennzahlen\n", "\n", "\n", "if __name__ == \"__main__\":\n", " wind, ist_flaute = erzeuge_windszenarien()\n", "\n", " print(\"=\" * 88)\n", " print(\" KRAFTWERKSEINSATZ: DER PLAN STEHT, BEVOR DER WIND WEHT\")\n", " print(\"=\" * 88)\n", " print(f\"{len(KRAFTWERKE)} Bloecke, {STUNDEN} Stunden, {SZENARIEN} Windszenarien \"\n", " f\"(davon {int(ist_flaute.sum())} Dunkelflauten).\")\n", " print(f\"Last {LAST.min():.0f} bis {LAST.max():.0f} MW, \"\n", " f\"Wind im Erwartungswert {WIND_ERWARTUNG.min():.0f} bis \"\n", " f\"{WIND_ERWARTUNG.max():.0f} MW.\")\n", " print(f\"Nicht gedeckte Last kostet {LASTABWURF:,.0f} EUR je MWh.\\n\")\n", "\n", " print(f\" {'Planungsgrundlage':<32} {'Kosten':>10} {'Kosten':>11} \"\n", " f\"{'Fehlmenge':>9} {'Fehlmenge':>10} {'Szen. mit':>12}\")\n", " print(f\" {'':<32} {'im Mittel':>10} {'schlimmst.':>11} \"\n", " f\"{'Mittel':>9} {'CVaR 90%':>10} {'Abwurf':>12}\")\n", " print(\" \" + \"-\" * 86)\n", "\n", " # --- 1. Deterministisch ----------------------------------------------\n", " plan_det, dauer_det, _ = plane(WIND_ERWARTUNG.reshape(1, -1), [1.0])\n", " det = zeige(\"nur Wind-Erwartungswert\", plan_det, wind)\n", "\n", " # --- 2. Zweistufig, risikoneutral ------------------------------------\n", " gleich = [1.0 / SZENARIEN] * SZENARIEN\n", " plan_sto, dauer_sto, optimal_sto = plane(wind, gleich)\n", " sto = zeige(\"alle 40 Szenarien\", plan_sto, wind)\n", "\n", " # --- 3. Mit Versorgungssicherheit ------------------------------------\n", " plan_sicher, dauer_sicher, optimal_sicher = plane(wind, gleich, cvar_grenze=0.0)\n", " sicher = zeige(\"+ Nichtdeckung CVaR 90 % = 0\", plan_sicher, wind)\n", "\n", " print(f\"\\n Rechenzeiten: {dauer_det:.1f}s / {dauer_sto:.1f}s / {dauer_sicher:.1f}s \"\n", " f\"(alle beweisbar optimal: {optimal_sto and optimal_sicher})\")\n", "\n", " # --- Was die Zeilen bedeuten -----------------------------------------\n", " print(\"\\n\" + \"-\" * 88)\n", " print(\"Was der Erwartungswert-Plan anrichtet\\n\")\n", " print(f\" Er ist auf dem Papier der billigste - gerechnet auf den mittleren\")\n", " print(f\" Wind kostet er weniger als jeder andere. In der Wirklichkeit liegt\")\n", " print(f\" er {det['mittel'] / sto['mittel'] - 1:+.0%} ueber dem zweistufigen Plan, und sein \"\n", " f\"schlimmster Tag\")\n", " print(f\" kostet {det['max'] / sto['max']:.1f}-mal so viel.\")\n", " print(f\"\\n Der Grund steht in der letzten Spalte: In {det['betroffen']} von {SZENARIEN} \"\n", " f\"Szenarien muss\")\n", " print(f\" Last abgeworfen werden. Was fehlt, sind nicht Kraftwerke, sondern\")\n", " print(f\" ANGEFAHRENE Kraftwerke - und ein Block mit acht Stunden\")\n", " print(f\" Mindestlaufzeit laesst sich um 18 Uhr nicht mehr herbeirufen.\")\n", "\n", " print(f\"\\n Der Unterschied im Plan ist klein:\")\n", " for name, plan in [(\"Erwartungswert\", plan_det), (\"zweistufig\", plan_sto),\n", " (\"mit Sicherheit\", plan_sicher)]:\n", " laufend = plan.sum(axis=0).astype(int)\n", " print(f\" {name:<16} {''.join(str(x) for x in laufend)} \"\n", " f\"({int(plan.sum())} Blockstunden)\")\n", " print(f\"\\n Der zweistufige Plan haelt {sto['blockstunden'] - det['blockstunden']} \"\n", " f\"Blockstunden mehr vor - das genuegt,\")\n", " print(f\" um alle {det['betroffen']} Lastabwuerfe zu vermeiden.\")\n", "\n", " # --- Der Preis der Sicherheit ----------------------------------------\n", " print(\"\\n\" + \"-\" * 88)\n", " print(\"Was Versorgungssicherheit kostet\\n\")\n", " print(f\" Bei {LASTABWURF:,.0f} EUR/MWh kostet die CVaR-Schranke NICHTS: Der\")\n", " print(\" risikoneutrale Plan haelt sie schon ein. Das ist kein Zufall - bei\")\n", " print(\" diesem Preis lohnt sich Vorhaltung bereits im Erwartungswert.\")\n", " print(\"\\n Interessant wird die Schranke dort, wo der Schaden ZU NIEDRIG\")\n", " print(\" bepreist ist. Dieselbe Rechnung mit einem Marktpreisdeckel von\")\n", " print(f\" {NIEDRIGER_ABWURFPREIS:,.0f} EUR/MWh statt {LASTABWURF:,.0f}:\\n\")\n", "\n", " print(f\" {'Planungsgrundlage':<32} {'Kosten':>10} {'Kosten':>11} \"\n", " f\"{'Fehlmenge':>9} {'Fehlmenge':>10} {'Szen. mit':>12}\")\n", " print(f\" {'':<32} {'im Mittel':>10} {'schlimmst.':>11} \"\n", " f\"{'Mittel':>9} {'CVaR 90%':>10} {'Abwurf':>12}\")\n", " print(\" \" + \"-\" * 86)\n", " plan_billig, _, _ = plane(wind, gleich, abwurfpreis=NIEDRIGER_ABWURFPREIS)\n", " billig = zeige(\"risikoneutral, billiger Abwurf\", plan_billig, wind,\n", " NIEDRIGER_ABWURFPREIS)\n", " plan_billig_sicher, _, _ = plane(wind, gleich, cvar_grenze=0.0,\n", " abwurfpreis=NIEDRIGER_ABWURFPREIS)\n", " billig_sicher = zeige(\"+ CVaR 90 % der Fehlmenge = 0\", plan_billig_sicher, wind,\n", " NIEDRIGER_ABWURFPREIS)\n", "\n", " aufpreis = billig_sicher[\"mittel\"] - billig[\"mittel\"]\n", " vermieden = billig[\"fehl_cvar\"] - billig_sicher[\"fehl_cvar\"]\n", " print(f\"\\n Jetzt greift die Schranke: Der risikoneutrale Plan nimmt in\")\n", " print(f\" {billig['betroffen']} von {SZENARIEN} Szenarien einen Lastabwurf in Kauf, weil er bei\")\n", " print(f\" {NIEDRIGER_ABWURFPREIS:,.0f} EUR/MWh billiger ist als das Vorhalten eines Blocks.\")\n", " print(f\"\\n Aufpreis fuer die Sicherheit: {aufpreis:,.0f} EUR je Tag\")\n", " if vermieden > 0.01:\n", " print(f\" Vermiedene Fehlmenge (CVaR 90 %): {vermieden:.1f} MWh\")\n", " print(f\" -> {aufpreis / vermieden:,.0f} EUR je vermiedener MWh\")\n", " print(f\"\\n Diese Zahl ist die Entscheidungsgrundlage - genau wie die Spalte\")\n", " print(f\" 'EUR je kg' im Kapitel Mehrziel. Sie sagt der Aufsicht, was ihre\")\n", " print(f\" Versorgungssicherheitsvorgabe tatsaechlich kostet, statt darueber\")\n", " print(f\" zu streiten, wie wichtig sie ist.\")\n", "\n", " print(\"\\n\" + \"=\" * 88)\n", " print(\" WAS MAN DARAUS MITNIMMT\")\n", " print(\"=\" * 88)\n", " print(\"1. Die erste Stufe ist binaer und traege. Was am Vorabend nicht\")\n", " print(\" angefahren wurde, steht am naechsten Tag nicht zur Verfuegung -\")\n", " print(\" egal, wie hoch der Preis dann steigt.\")\n", " print(\"2. Deshalb ist der Erwartungswert als Planungsgrundlage nicht nur\")\n", " print(\" ungenau, sondern SYSTEMATISCH zu knapp: Er plant fuer einen Tag,\")\n", " print(\" der so nie eintritt, und laesst keine Reserve fuer die Haelfte\")\n", " print(\" der Faelle, in denen es schlechter kommt.\")\n", " print(\"3. Der Ausweg ist keine bessere Windprognose, sondern ein Modell,\")\n", " print(\" das die Szenarien SIEHT. Die Rechenzeit dafuer liegt hier im\")\n", " print(\" Sekundenbereich - der Aufwand ist die Modellierung, nicht der\")\n", " print(\" Solver.\")\n", " print(\"4. Wo der Lastabwurf teuer genug bepreist ist, deckt schon die\")\n", " print(\" Minimierung der ERWARTETEN Kosten die Extremfaelle mit ab. Die\")\n", " print(\" Risikoschranke wird genau dann gebraucht, wenn der Schaden ZU\")\n", " print(\" NIEDRIG bepreist ist - was bei Versorgungssicherheit die Regel\")\n", " print(\" ist, weil ein Marktpreisdeckel nicht den volkswirtschaftlichen\")\n", " print(\" Schaden abbildet. Dann liefert sie den Preis der Vorgabe in Euro\")\n", " print(\" je vermiedener MWh, und darueber laesst sich verhandeln.\")\n", " print(\"=\" * 88)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.11" } }, "nbformat": 4, "nbformat_minor": 5 }