operations_research/OR_HTML_04/programme/Chance_Constraints.py

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