#!/usr/bin/env python3 # Chance_Constraints.py """ Kapitel Unsicherheit: Wahrscheinlichkeitsbeschraenkungen. Das Management fragt selten "was ist im Mittel am besten?", sondern "mit welcher Sicherheit haelt der Plan?". Genau das formuliert eine Chance Constraint: P(Versorgung >= Bedarf) >= 1 - alpha. Gezeigt werden beide Wege dorthin, an derselben Instanz: (a) analytisch - unter Normalverteilungsannahme wird daraus eine Second-Order-Cone-Bedingung, loesbar mit CVXPY (b) szenariobasiert - Big-M mit Binaervariablen, fuer beliebige empirische Verteilungen Und die beiden Messungen, auf die es ankommt: Was kostet ein Prozentpunkt Versorgungssicherheit - und haelt die Zusage auch dann, wenn die Verteilung nicht normal ist? Solver: CLARABEL (Kegel) und SciPy/HiGHS (gemischt-ganzzahlig). Weder ortools noch ein direkter highspy-Import, damit alles in einem Prozess laeuft. """ import cvxpy as cp import numpy as np from scipy.stats import norm # --- Der Kraftwerkspark ---------------------------------------------------- # Die Verfuegbarkeit ist der Anteil, den eine installierte MW im Mittel # wirklich liefert: bei Wind und Sonne klein und stark schwankend, bei Gas und # Biomasse gross und stabil. Die Ausbaugrenze ist der Standort - Flaeche, # Genehmigung, Brennstoffversorgung. TECHNIK = ["Gaskraftwerk", "Windpark", "Solarpark", "Biomasse"] KOSTEN = np.array([65_000.0, 24_000.0, 12_000.0, 60_000.0]) # EUR je MW und Jahr VERFUEGBAR = np.array([0.92, 0.35, 0.18, 0.85]) STREUUNG = np.array([0.05, 0.16, 0.10, 0.04]) GRENZE = np.array([400.0, 250.0, 600.0, 350.0]) # MW # Wind und Sonne sind leicht gegenlaeufig: Ein Tiefdruckgebiet bringt Wind und # Wolken zugleich. Genau diese Korrelation macht die Bedingung zu einem Kegel # und nicht zu einer Summe unabhaengiger Einzelzuschlaege. KORRELATION = np.array([ [1.00, 0.00, 0.00, 0.05], [0.00, 1.00, -0.25, 0.00], [0.00, -0.25, 1.00, 0.00], [0.05, 0.00, 0.00, 1.00], ]) SIGMA = np.diag(STREUUNG) @ KORRELATION @ np.diag(STREUUNG) WURZEL = np.linalg.cholesky(SIGMA) # L mit L @ L.T == SIGMA BEDARF = 500.0 # MW, die gesichert bereitstehen muessen N = len(TECHNIK) # Die Kaeltewelle mit Dunkelflaute: selten, aber sie trifft alles zugleich. # Wind und Sonne brechen fast vollstaendig weg - und, das ist der Punkt, das # Gaskraftwerk liefert ebenfalls weniger, weil bei Frost der Netzdruck faellt. # Eine Kovarianzmatrix mit Korrelationen um null kann das nicht ausdruecken. P_KAELTEWELLE = 0.08 EINBRUCH = np.array([0.72, 0.10, 0.15, 0.88]) def mittelwertplan(): """Plant mit den Erwartungswerten - ignoriert die Streuung vollstaendig.""" x = cp.Variable(N, nonneg=True) problem = cp.Problem(cp.Minimize(KOSTEN @ x), [VERFUEGBAR @ x >= BEDARF, x <= GRENZE]) problem.solve(solver=cp.CLARABEL) return problem.value, x.value def chance_constraint_analytisch(alpha): """ P(a^T x >= BEDARF) >= 1 - alpha unter a ~ N(VERFUEGBAR, SIGMA). Aequivalent zu: VERFUEGBAR^T x - z * ||L^T x||_2 >= BEDARF mit z = Phi^-1(1 - alpha). Das ist eine Second-Order-Cone-Bedingung: Der Sicherheitszuschlag ist eine Norm ueber x, kein fester Aufschlag je Anlage. Deshalb belohnt sie Mischung - zwei gegenlaeufige Quellen schwanken gemeinsam weniger als jede fuer sich. """ x = cp.Variable(N, nonneg=True) z = norm.ppf(1 - alpha) bedingung = VERFUEGBAR @ x - z * cp.norm(WURZEL.T @ x, 2) >= BEDARF problem = cp.Problem(cp.Minimize(KOSTEN @ x), [bedingung, x <= GRENZE]) problem.solve(solver=cp.CLARABEL) return problem.value, x.value def chance_constraint_szenarien(a, alpha): """ Dieselbe Zusage ohne Verteilungsannahme: Fuer jedes Szenario s sagt eine Binaervariable z_s, ob es verletzt werden darf. a_s^T x >= BEDARF - M * z_s fuer alle s sum_s z_s <= alpha * S M = BEDARF ist die kleinstmoegliche gueltige Schranke, denn a_s^T x >= 0: Groesser kann eine Verletzung gar nicht ausfallen. Ein unnoetig grosses M wuerde die LP-Relaxierung aufweichen und die Suche verlangsamen - die Big-M-Falle aus dem Kapitel Gemischt-ganzzahlige Optimierung. """ S = len(a) x = cp.Variable(N, nonneg=True) z = cp.Variable(S, boolean=True) problem = cp.Problem( cp.Minimize(KOSTEN @ x), [a @ x >= BEDARF - BEDARF * z, cp.sum(z) <= alpha * S, x <= GRENZE]) problem.solve(solver=cp.SCIPY) return problem.value, x.value def ziehe_wetter(anzahl, seed): """Mischverteilung: normales Wetter, mit P_KAELTEWELLE ein Einbruch.""" rng = np.random.default_rng(seed) a = rng.multivariate_normal(VERFUEGBAR, SIGMA, size=anzahl) getroffen = rng.random(anzahl) < P_KAELTEWELLE a[getroffen] *= EINBRUCH return np.clip(a, 0.0, 1.0) def sicherheit(x, a): """Anteil der Szenarien, in denen der Plan den Bedarf deckt.""" return float(np.mean(a @ x >= BEDARF)) KOPF = (f"{'Plan':<20} {'Kosten (EUR)':>13} {'Gas':>4} {'Wind':>4} " f"{'Sol':>4} {'Bio':>4} {'Normalwelt':>11} {'echte Welt':>11}") def zeile(name, kosten, x, in_normal, in_echt): mix = " ".join(f"{w:4.0f}" for w in x) return (f"{name:<20} {kosten:>13,.0f} {mix} " f"{in_normal * 100:>9.2f} % {in_echt * 100:>9.2f} %") if __name__ == "__main__": print("=" * 78) print(" WAHRSCHEINLICHKEITSBESCHRAENKUNGEN IM KRAFTWERKSPARK") print("=" * 78) print(f"Gesichert bereitzustellen: {BEDARF:.0f} MW\n") print(f"{'Technologie':<14} {'EUR/MW':>8} {'verfuegbar':>11} {'Streuung':>9} " f"{'Grenze':>8} {'EUR je erw. MW':>15}") print("-" * 78) for i, name in enumerate(TECHNIK): print(f"{name:<14} {KOSTEN[i]:>8,.0f} {VERFUEGBAR[i]:>10.2f} " f"{STREUUNG[i]:>8.2f} {GRENZE[i]:>7.0f} " f"{KOSTEN[i] / VERFUEGBAR[i]:>15,.0f}") print(f"\nVollausbau liefert im Mittel {VERFUEGBAR @ GRENZE:.0f} MW, " f"in der Kaeltewelle {(VERFUEGBAR * EINBRUCH) @ GRENZE:.0f} MW.") # Zwei grosse, unabhaengige Testmengen. An ihnen wird JEDER Plan gemessen - # einmal in der unterstellten Normalwelt, einmal in der echten Verteilung. echt = ziehe_wetter(200_000, seed=771) normalwelt = np.random.default_rng(4711).multivariate_normal( VERFUEGBAR, SIGMA, size=200_000) print("\n" + "=" * 78) print(" (1) Mittelwertplan und Chance Constraints im Vergleich") print("=" * 78) print(KOPF) print("-" * 78) k_mittel, x_mittel = mittelwertplan() stufen = [("Mittelwert", k_mittel, sicherheit(x_mittel, normalwelt))] print(zeile("Mittelwert", k_mittel, x_mittel, sicherheit(x_mittel, normalwelt), sicherheit(x_mittel, echt))) plaene = {} for alpha in (0.20, 0.10, 0.05, 0.01): kosten, x = chance_constraint_analytisch(alpha) plaene[alpha] = (kosten, x) quote_normal = sicherheit(x, normalwelt) stufen.append((f"Zusage {(1 - alpha) * 100:.0f} %", kosten, quote_normal)) print(zeile(f"Zusage {(1 - alpha) * 100:.0f} %", kosten, x, quote_normal, sicherheit(x, echt))) print("\n In der Normalwelt trifft jede Zusage ihren Wert - das Verfahren") print(" rechnet richtig. Die Spalte 'echte Welt' kommt in Teil (3).") print("\n" + "=" * 78) print(" (2) Was kostet ein Prozentpunkt Versorgungssicherheit?") print("=" * 78) print(f"{'von -> nach':<26} {'Prozentpunkte':>13} {'Mehrkosten':>13} " f"{'EUR je Punkt':>13}") print("-" * 78) for (n1, k1, q1), (n2, k2, q2) in zip(stufen, stufen[1:]): d_punkte = (q2 - q1) * 100 d_kosten = k2 - k1 print(f"{n1 + ' -> ' + n2:<26} {d_punkte:>13.2f} {d_kosten:>13,.0f} " f"{d_kosten / d_punkte:>13,.0f}") erst = (stufen[1][1] - stufen[0][1]) / ((stufen[1][2] - stufen[0][2]) * 100) letzt = (stufen[-1][1] - stufen[-2][1]) / ((stufen[-1][2] - stufen[-2][2]) * 100) print(f"\n Der letzte Prozentpunkt kostet das {letzt / erst:.1f}-fache des ersten.") print(" Sicherheit ist konvex bepreist - genau deshalb muss jemand") print(" entscheiden, wie viel davon das Unternehmen kaufen will.") print("\n" + "=" * 78) print(" (3) Und wenn die Verteilung nicht normal ist?") print("=" * 78) k_soc, x_soc = plaene[0.05] print(f"Die echte Welt kennt die Kaeltewelle mit Dunkelflaute: in " f"{P_KAELTEWELLE * 100:.0f} % der Faelle liefern") print(f"Wind {EINBRUCH[1] * 100:.0f} %, Sonne {EINBRUCH[2] * 100:.0f} % " f"und - entscheidend - auch das Gaskraftwerk nur " f"{EINBRUCH[0] * 100:.0f} %") print("des Ueblichen. Alles bricht gleichzeitig ein.\n") a_bau = ziehe_wetter(400, seed=20260908) k_sz, x_sz = chance_constraint_szenarien(a_bau, 0.05) print(KOPF) print("-" * 78) print(zeile("SOC, Zusage 95 %", k_soc, x_soc, sicherheit(x_soc, normalwelt), sicherheit(x_soc, echt))) print(zeile("Szenarien, 95 %", k_sz, x_sz, sicherheit(x_sz, normalwelt), sicherheit(x_sz, echt))) print(f"\n Der SOC-Plan verspricht 95 % und haelt " f"{sicherheit(x_soc, echt) * 100:.1f} %. Nicht das Verfahren ist") print(f" falsch, sondern die Annahme: In der Normalwelt liefert derselbe") print(f" Plan {sicherheit(x_soc, normalwelt) * 100:.1f} %.") print(f"\n Der Szenarioplan ({len(a_bau)} Szenarien, {len(a_bau)} " f"Binaervariablen) kommt auf") print(f" {sicherheit(x_sz, echt) * 100:.1f} % - und kostet dafuer " f"{(k_sz / k_soc - 1) * 100:+.1f} %. Er kauft keine Windkraft mehr:") print(" Was im Ernstfall ausfaellt, hilft der Zusage nicht.") print("=" * 78)