Phase 6.1: Chance Constraints - die Zusage "mit 95 % Sicherheit"
Neuer Abschnitt im Kapitel Unsicherheit plus Chance_Constraints.py (74. Programm). Das Kapitel hatte Monte-Carlo, Zweistufigkeit und Worst-Case-Robustheit; die Wahrscheinlichkeitszusage war die fehlende vierte Antwort - und die, nach der das Management tatsaechlich fragt. Setzt Paket 4 aus Verbesserungen_02.md um. Beide Wege an DERSELBEN Instanz (Kraftwerkspark, 500 MW gesicherte Zusage): analytisch als Second-Order-Cone-Bedingung (CLARABEL) und szenariobasiert als Big-M-MILP (SciPy/HiGHS) - beides in einem Prozess, ohne ortools- oder highspy-Import. Drei gemessene Befunde: * Der Mittelwertplan haelt 50,08 %. Kein Fehler, sondern die Definition des Erwartungswerts. * Sicherheit ist konvex bepreist: 56.289 EUR je Prozentpunkt auf dem Weg zu 80 %, 253.848 EUR zwischen 95 und 99 % - das 4,5-fache. Das Programm rechnet die Tabelle selbst aus, statt sie zu behaupten. * Die Zusage gilt nur fuer die unterstellte Verteilung: Der 95-%-Plan haelt gemessen 87,44 %, sobald die Testverteilung eine Kaeltewelle mit Dunkelflaute enthaelt, in der alles zugleich einbricht - auch das Gaskraftwerk. In einer reinen Normalwelt liefert derselbe Plan 94,97 %. Zwei eigene Fehlgriffe, beide durch Messen aufgefallen und korrigiert: Die erste Kostentarierung ergab eine entartete Loesung (alles ins Gaskraftwerk), womit die Kovarianzmatrix wirkungslos war - und gerade sie begruendet die Kegelform. Und der erste Kaelteeinbruch traf nur Wind und Sonne; der SOC-Plan hatte die ohnehin herausgehalten und war zufaellig robust, das Argument trug nicht. Ehrlich berichtet statt geglaettet: Die Szenariomethode ueberanpasst. Ueber zwoelf Laeufe (S = 200, 400, 800) lag die tatsaechliche Quote zwischen 93,14 % und 96,71 %, und die Spanne wurde von S = 400 auf 800 wieder breiter. Beide Solver bestaetigen optimal bei identischen Kosten - echte Ueberanpassung, kein Solverartefakt. Steht als Warnkasten im Abschnitt. Die Laufzeitmessung ist aus der Ausgabe entfernt: Eine Wanduhrzeit ist nie byteidentisch reproduzierbar (2,2 s / 2,3 s zwischen zwei Laeufen) und haette Regel 4 dauerhaft gebrochen. Danach drei Laeufe zeichengleich, und die abgedruckte Ausgabe stimmt mit dem Lauf des extrahierten Programms ueberein. Mitgezogen: Kapitelkopf, Lernziele, Uebersichtstabelle (drei -> vier Ansaetze), Selbsttest, Zusammenfassung, Vorwort-Programmverzeichnis, eine neue Uebungsaufgabe und ihre Loesung in Anhang A. pyproject.toml korrigiert: cvxpy lag in [finance], wird aber von 13 Programmen gebraucht, darunter dreien in diesem Kapitel - jetzt in der Grundausstattung. Es zieht kein highspy nach, der Solverkonflikt bleibt auf [large-scale] beschraenkt. Stand: 293 Abschnitte, 712 Querverweise, 327 Indexmarken, 74 Programme (0 Fehler), 33 pytest-Tests, PDF 725 Seiten. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This commit is contained in:
dschlueter
parent
62952ba066
commit
6c67e126dd
23 changed files with 9089 additions and 6419 deletions
228
OR_HTML_04/programme/Chance_Constraints.py
Normal file
228
OR_HTML_04/programme/Chance_Constraints.py
Normal file
|
|
@ -0,0 +1,228 @@
|
|||
#!/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)
|
||||
Loading…
Reference in a new issue