operations_research/OR_HTML_04/Notebooks_04/unsicherheit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kapitel 12: Optimierung unter Unsicherheit — Monte-Carlo, Stochastik, Robustheit\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": [
    "## Der Fluch des Durchschnitts\n",
    "\n",
    "`Fluch_des_Durchschnitts.py`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "\n",
    "# Fluch_des_Durchschnitts.py\n",
    "\"\"\"\n",
    "Kapitel Unsicherheit: Der Fluch des Durchschnitts (die Handrechnung dazu) als Scan ueber alle\n",
    "moeglichen Kapazitaeten - zeigt, dass das Optimum nicht beim Mittelwert liegt.\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "SZENARIEN = np.array([100, 250, 500])\n",
    "WAHRSCHEINLICHKEITEN = np.array([0.5, 0.3, 0.2])\n",
    "PREIS_VORAB = 40.0\n",
    "PREIS_ZUKAUF = 120.0\n",
    "PREIS_VERWALTUNG = 5.0\n",
    "\n",
    "ERWARTUNGSWERT = float(SZENARIEN @ WAHRSCHEINLICHKEITEN)\n",
    "\n",
    "\n",
    "def erwartete_kosten(x):\n",
    "    kosten = np.where(\n",
    "        SZENARIEN >= x,\n",
    "        PREIS_VORAB * x + PREIS_ZUKAUF * (SZENARIEN - x),\n",
    "        PREIS_VORAB * x + PREIS_VERWALTUNG * (x - SZENARIEN),\n",
    "    )\n",
    "    return float(kosten @ WAHRSCHEINLICHKEITEN)\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    print(\"=\" * 70)\n",
    "    print(\"  DER FLUCH DES DURCHSCHNITTS\")\n",
    "    print(\"=\" * 70)\n",
    "    print(f\"Erwarteter Bedarf (Mittelwert): {ERWARTUNGSWERT:.1f}\")\n",
    "    print(f\"Kosten bei naiver Planung x={ERWARTUNGSWERT:.0f}: \"\n",
    "          f\"{erwartete_kosten(ERWARTUNGSWERT):,.2f} EUR\\n\")\n",
    "\n",
    "    x_werte = np.arange(0, 501, 1)\n",
    "    kosten_werte = np.array([erwartete_kosten(x) for x in x_werte])\n",
    "    x_optimal = x_werte[np.argmin(kosten_werte)]\n",
    "    kosten_optimal = kosten_werte.min()\n",
    "\n",
    "    print(f\"Optimales x (durch Scan gefunden): {x_optimal}\")\n",
    "    print(f\"Kosten beim Optimum:                {kosten_optimal:,.2f} EUR\")\n",
    "    print(f\"Ersparnis gegenueber naiver Planung: \"\n",
    "          f\"{erwartete_kosten(ERWARTUNGSWERT) - kosten_optimal:,.2f} EUR\\n\")\n",
    "\n",
    "    print(\"Kosten fuer ausgewaehlte x zum Vergleich:\")\n",
    "    for x in [100, 225, 250, 300, 500]:\n",
    "        markierung = \"  <- Mittelwert\" if x == 225 else (\"  <- Optimum\" if x == 250 else \"\")\n",
    "        print(f\"  x={x:4d}: {erwartete_kosten(x):>12,.2f} EUR{markierung}\")\n",
    "\n",
    "    print(\"\\nDas Optimum liegt exakt auf einem Szenariowert (250 = 'Volatil'),\")\n",
    "    print(\"nicht beim Mittelwert 225 - typisch fuer asymmetrische Kostenfunktionen.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monte-Carlo-Simulation\n",
    "\n",
    "`Monte_Carlo.py`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "\n",
    "# Monte_Carlo.py\n",
    "\"\"\"\n",
    "Kapitel Unsicherheit: Monte-Carlo-Bewertung von Kapazitaetsplaenen.\n",
    "\n",
    "Monte Carlo OPTIMIERT nicht - es BEWERTET. Der Nutzen liegt darin, dass man\n",
    "beliebige Kennzahlen ablesen kann: Erwartungswert, Quantile, Ausfallwahr-\n",
    "scheinlichkeit, Worst Case. Genau diese Groessen braucht man, um zwischen\n",
    "Plaenen zu entscheiden.\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "KOSTEN_VORAB = 40.0        # EUR je Einheit, im Voraus gekauft\n",
    "KOSTEN_SPOT = 120.0        # EUR je Einheit, kurzfristig zugekauft\n",
    "KOSTEN_LEERLAUF = 5.0      # EUR je ungenutzter Einheit\n",
    "\n",
    "ANZAHL_ZIEHUNGEN = 100_000\n",
    "\n",
    "\n",
    "def ziehe_bedarf(rng, ziehungen):\n",
    "    \"\"\"\n",
    "    Bedarfsmodell: Mischverteilung aus Normalbetrieb und seltenen Lastspitzen.\n",
    "    Realistischer als eine reine Normalverteilung - Krisen sind selten,\n",
    "    aber extrem (fat tail).\n",
    "    \"\"\"\n",
    "    normal = rng.normal(loc=150, scale=40, size=ziehungen)\n",
    "    spitze = rng.normal(loc=450, scale=80, size=ziehungen)\n",
    "    ist_spitze = rng.random(ziehungen) < 0.15            # 15 % Lastspitzen\n",
    "    return np.maximum(np.where(ist_spitze, spitze, normal), 0.0)\n",
    "\n",
    "\n",
    "def kosten_fuer(kapazitaet, bedarf):\n",
    "    \"\"\"Gesamtkosten je Szenario fuer eine gegebene Vorabkapazitaet.\"\"\"\n",
    "    unterdeckung = np.maximum(bedarf - kapazitaet, 0.0)\n",
    "    ueberdeckung = np.maximum(kapazitaet - bedarf, 0.0)\n",
    "    return (KOSTEN_VORAB * kapazitaet\n",
    "            + KOSTEN_SPOT * unterdeckung\n",
    "            + KOSTEN_LEERLAUF * ueberdeckung)\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    rng = np.random.default_rng(2026)\n",
    "    bedarf = ziehe_bedarf(rng, ANZAHL_ZIEHUNGEN)\n",
    "\n",
    "    print(\"=\" * 88)\n",
    "    print(f\"  MONTE-CARLO-BEWERTUNG ({ANZAHL_ZIEHUNGEN:,} Szenarien)\")\n",
    "    print(\"=\" * 88)\n",
    "    print(f\"Bedarfsverteilung: Mittelwert {bedarf.mean():.1f} | \"\n",
    "          f\"Median {np.median(bedarf):.1f} | \"\n",
    "          f\"95%-Quantil {np.percentile(bedarf, 95):.1f} | \"\n",
    "          f\"Maximum {bedarf.max():.1f}\")\n",
    "    print(\"Der Median liegt deutlich unter dem Mittelwert - die Verteilung ist\")\n",
    "    print(\"rechtsschief. Genau hier fuehrt Planung mit dem Mittelwert in die Irre.\\n\")\n",
    "\n",
    "    print(f\"{'Kapazitaet':>10} | {'Erw. Kosten':>12} | {'Median':>10} | \"\n",
    "          f\"{'95%-Quantil':>12} | {'Unterdeckung':>12}\")\n",
    "    print(\"-\" * 88)\n",
    "\n",
    "    kandidaten = [150, 200, 225, 250, 300, 350, 400]\n",
    "    ergebnisse = []\n",
    "    for kapazitaet in kandidaten:\n",
    "        kosten = kosten_fuer(kapazitaet, bedarf)\n",
    "        p_unterdeckung = float(np.mean(bedarf > kapazitaet))\n",
    "        ergebnisse.append((kapazitaet, kosten.mean(), p_unterdeckung))\n",
    "        print(f\"{kapazitaet:>10} | {kosten.mean():>12,.0f} | \"\n",
    "              f\"{np.median(kosten):>10,.0f} | {np.percentile(kosten, 95):>12,.0f} | \"\n",
    "              f\"{p_unterdeckung*100:>11.1f} %\")\n",
    "\n",
    "    beste = min(ergebnisse, key=lambda t: t[1])\n",
    "    print(\"-\" * 88)\n",
    "    print(f\"Bester Kandidat: Kapazitaet {beste[0]} mit erwarteten Kosten \"\n",
    "          f\"{beste[1]:,.0f} EUR\")\n",
    "\n",
    "    # --- Feinsuche ueber ein Raster ---------------------------------------\n",
    "    raster = np.arange(100, 500, 5)\n",
    "    erwartete = np.array([kosten_fuer(k, bedarf).mean() for k in raster])\n",
    "    optimum = raster[int(np.argmin(erwartete))]\n",
    "    print(f\"Feinsuche (Raster 100..500): Optimum bei Kapazitaet {optimum}, \"\n",
    "          f\"Kosten {erwartete.min():,.0f} EUR\")\n",
    "\n",
    "    # --- Vergleich mit der naiven Mittelwertplanung ----------------------\n",
    "    naiv = int(round(bedarf.mean()))\n",
    "    kosten_naiv = kosten_fuer(naiv, bedarf).mean()\n",
    "    kosten_opt = kosten_fuer(optimum, bedarf).mean()\n",
    "    print(\"\\n--- Fluch des Durchschnitts, gemessen ---\")\n",
    "    print(f\"  Planung mit Mittelwert ({naiv}):  {kosten_naiv:,.0f} EUR\")\n",
    "    print(f\"  Monte-Carlo-Optimum    ({optimum}):  {kosten_opt:,.0f} EUR\")\n",
    "    print(f\"  Mehrkosten der naiven Planung:  {kosten_naiv - kosten_opt:,.0f} EUR \"\n",
    "          f\"({(kosten_naiv/kosten_opt - 1)*100:.1f} %)\")\n",
    "    print(\"=\" * 88)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zweistufige stochastische Programmierung\n",
    "\n",
    "`Stochastische_Optimierung.py`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "\n",
    "# Stochastische_Optimierung.py\n",
    "\"\"\"\n",
    "Kapitel Unsicherheit: Two-Stage Stochastic Programming mit CVXPY.\n",
    "\n",
    "Eigenschaften:\n",
    "  * Die Analyse am Ende wird BERECHNET statt fest verdrahtet\n",
    "  * Vergleich gegen drei Alternativen: Mittelwert, Worst Case, perfekte Voraussicht\n",
    "  * Kennzahl EVPI (Wert perfekter Information) wird ausgewiesen\n",
    "\"\"\"\n",
    "\n",
    "import cvxpy as cp\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "SZENARIEN = [\"Ruhig\", \"Volatil\", \"Crash\"]\n",
    "WAHRSCHEINLICHKEIT = np.array([0.50, 0.30, 0.20])\n",
    "BEDARF = np.array([100.0, 250.0, 500.0])\n",
    "\n",
    "KOSTEN_VORAB = 40.0\n",
    "KOSTEN_SPOT = 120.0\n",
    "KOSTEN_LEERLAUF = 5.0\n",
    "\n",
    "S = len(SZENARIEN)\n",
    "\n",
    "\n",
    "def loese_stochastisch():\n",
    "    \"\"\"Zweistufiges Modell: eine Vorabentscheidung, szenarioabhaengige Korrektur.\"\"\"\n",
    "    x = cp.Variable(nonneg=True, name=\"Basiskapazitaet\")        # Stufe 1\n",
    "    y_spot = cp.Variable(S, nonneg=True, name=\"Spot_Zukauf\")    # Stufe 2\n",
    "    y_leer = cp.Variable(S, nonneg=True, name=\"Leerlauf\")       # Stufe 2\n",
    "\n",
    "    # Kopplung: Basis + Zukauf - Leerlauf == Bedarf (je Szenario)\n",
    "    nebenbedingungen = [x + y_spot[s] - y_leer[s] == BEDARF[s] for s in range(S)]\n",
    "\n",
    "    erwartete_korrektur = sum(\n",
    "        WAHRSCHEINLICHKEIT[s] * (KOSTEN_SPOT * y_spot[s] + KOSTEN_LEERLAUF * y_leer[s])\n",
    "        for s in range(S))\n",
    "\n",
    "    problem = cp.Problem(cp.Minimize(KOSTEN_VORAB * x + erwartete_korrektur),\n",
    "                         nebenbedingungen)\n",
    "    problem.solve()\n",
    "    return problem, x, y_spot, y_leer\n",
    "\n",
    "\n",
    "def kosten_bei(kapazitaet):\n",
    "    \"\"\"Erwartete Gesamtkosten fuer eine fest vorgegebene Kapazitaet.\"\"\"\n",
    "    unter = np.maximum(BEDARF - kapazitaet, 0.0)\n",
    "    ueber = np.maximum(kapazitaet - BEDARF, 0.0)\n",
    "    je_szenario = KOSTEN_VORAB * kapazitaet + KOSTEN_SPOT * unter + KOSTEN_LEERLAUF * ueber\n",
    "    return float(WAHRSCHEINLICHKEIT @ je_szenario)\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    problem, x, y_spot, y_leer = loese_stochastisch()\n",
    "    kapazitaet = float(x.value)\n",
    "    mittelwert = float(WAHRSCHEINLICHKEIT @ BEDARF)\n",
    "\n",
    "    print(\"=\" * 82)\n",
    "    print(\"     STOCHASTISCHE TWO-STAGE OPTIMIERUNG (KAPAZITAETSPLANUNG)\")\n",
    "    print(\"=\" * 82)\n",
    "    print(f\"Status: {problem.status}\")\n",
    "    print(f\"Erwarteter Bedarf (Mittelwert):   {mittelwert:.1f} Einheiten\")\n",
    "    print(f\"Optimale Stufe-1-Kapazitaet x*:   {kapazitaet:.1f} Einheiten\")\n",
    "    print(f\"Minimale erwartete Gesamtkosten:  {problem.value:,.2f} EUR\\n\")\n",
    "\n",
    "    tabelle = pd.DataFrame({\n",
    "        \"Szenario\": SZENARIEN,\n",
    "        \"Wahrsch.\": [f\"{p*100:.0f} %\" for p in WAHRSCHEINLICHKEIT],\n",
    "        \"Bedarf\": BEDARF,\n",
    "        \"Basis genutzt\": [min(kapazitaet, b) for b in BEDARF],\n",
    "        \"Spot-Zukauf\": np.round(y_spot.value, 1),\n",
    "        \"Leerlauf\": np.round(y_leer.value, 1),\n",
    "        \"Kosten (EUR)\": [f\"{KOSTEN_VORAB*kapazitaet + KOSTEN_SPOT*y_spot.value[s] + KOSTEN_LEERLAUF*y_leer.value[s]:,.0f}\"\n",
    "                         for s in range(S)],\n",
    "    })\n",
    "    print(tabelle.to_string(index=False))\n",
    "\n",
    "    # --- Vergleich mit Alternativstrategien (berechnet, nicht behauptet) --\n",
    "    print(\"\\n\" + \"-\" * 82)\n",
    "    print(\"Vergleich verschiedener Planungsstrategien:\")\n",
    "    print(f\"{'Strategie':<34} {'Kapazitaet':>11} {'Erw. Kosten':>14} {'Mehrkosten':>13}\")\n",
    "    print(\"-\" * 82)\n",
    "\n",
    "    optimal = problem.value\n",
    "    strategien = [\n",
    "        (\"Stochastisch optimal\", kapazitaet),\n",
    "        (\"Naiv: Mittelwert einsetzen\", mittelwert),\n",
    "        (\"Vorsichtig: Worst Case abdecken\", float(BEDARF.max())),\n",
    "        (\"Optimistisch: Bestfall\", float(BEDARF.min())),\n",
    "    ]\n",
    "    for name, kap in strategien:\n",
    "        kosten = kosten_bei(kap)\n",
    "        print(f\"{name:<34} {kap:>11.1f} {kosten:>14,.0f} \"\n",
    "              f\"{kosten - optimal:>+13,.0f}\")\n",
    "\n",
    "    # --- EVPI: Was waere perfekte Voraussicht wert? ----------------------\n",
    "    # Bei perfekter Information wuerde man je Szenario genau den Bedarf kaufen.\n",
    "    kosten_perfekt = float(WAHRSCHEINLICHKEIT @ (KOSTEN_VORAB * BEDARF))\n",
    "    evpi = optimal - kosten_perfekt\n",
    "    print(\"-\" * 82)\n",
    "    print(f\"Kosten bei perfekter Voraussicht:      {kosten_perfekt:>10,.0f} EUR\")\n",
    "    print(f\"Wert perfekter Information (EVPI):     {evpi:>10,.0f} EUR \"\n",
    "          f\"({evpi/optimal*100:.1f} % der Kosten)\")\n",
    "    print(\"  -> So viel duerfte eine perfekte Bedarfsprognose hoechstens kosten.\")\n",
    "\n",
    "    # --- Automatische Interpretation --------------------------------------\n",
    "    print(\"-\" * 82)\n",
    "    if kapazitaet > mittelwert + 1e-6:\n",
    "        print(f\"Analyse: Der Solver waehlt {kapazitaet:.0f} Einheiten und damit MEHR als\")\n",
    "        print(f\"den Mittelwert ({mittelwert:.0f}), weil Unterdeckung ({KOSTEN_SPOT:.0f} EUR)\")\n",
    "        print(f\"deutlich teurer ist als Leerlauf ({KOSTEN_LEERLAUF:.0f} EUR).\")\n",
    "    elif kapazitaet < mittelwert - 1e-6:\n",
    "        print(f\"Analyse: Der Solver waehlt {kapazitaet:.0f} und damit WENIGER als den\")\n",
    "        print(f\"Mittelwert ({mittelwert:.0f}) - Leerlauf ist hier teurer als Zukauf.\")\n",
    "    else:\n",
    "        print(\"Analyse: Kapazitaet entspricht dem Mittelwert (symmetrische Kosten).\")\n",
    "    print(\"=\" * 82)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Robuste Optimierung: gegen den Worst Case absichern\n",
    "\n",
    "`Robuste_Optimierung.py`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "\n",
    "# Robuste_Optimierung.py\n",
    "\"\"\"\n",
    "Kapitel Unsicherheit: Robuste Portfolio-Optimierung.\n",
    "\n",
    "Vergleicht drei Haltungen zur Unsicherheit:\n",
    "  (a) nominal   - vertraut den Punktschaetzungen blind\n",
    "  (b) robust    - sichert gegen Box-Unsicherheit ab (Worst Case)\n",
    "  (c) stochastisch - optimiert den Erwartungswert ueber Szenarien\n",
    "\n",
    "und misst den \"Preis der Robustheit\": Wie viel Ertrag kostet die Absicherung\n",
    "im Normalfall - und wie viel Verlust erspart sie im Ernstfall?\n",
    "\"\"\"\n",
    "\n",
    "import cvxpy as cp\n",
    "import numpy as np\n",
    "\n",
    "ASSETS = [\"Aktien Welt\", \"Anleihen\", \"Rohstoffe\", \"Immobilien\"]\n",
    "MU_SCHAETZUNG = np.array([0.085, 0.030, 0.055, 0.060])   # Punktschaetzung\n",
    "UNSICHERHEIT = np.array([0.040, 0.008, 0.045, 0.025])    # +/- delta je Titel\n",
    "VOLA = np.array([0.17, 0.05, 0.22, 0.12])\n",
    "KORR = np.array([\n",
    "    [1.00, -0.15, 0.35, 0.55],\n",
    "    [-0.15, 1.00, -0.05, 0.10],\n",
    "    [0.35, -0.05, 1.00, 0.25],\n",
    "    [0.55, 0.10, 0.25, 1.00],\n",
    "])\n",
    "SIGMA = np.diag(VOLA) @ KORR @ np.diag(VOLA)\n",
    "LAMBDA = 4.0        # Risikoaversion\n",
    "N = len(ASSETS)\n",
    "\n",
    "\n",
    "def optimiere(mu_effektiv):\n",
    "    \"\"\"Standard-Mean-Variance mit vorgegebenem Renditevektor.\"\"\"\n",
    "    w = cp.Variable(N, nonneg=True)\n",
    "    ziel = cp.Maximize(mu_effektiv @ w - 0.5 * LAMBDA * cp.quad_form(w, SIGMA))\n",
    "    problem = cp.Problem(ziel, [cp.sum(w) == 1])\n",
    "    problem.solve()\n",
    "    return w.value\n",
    "\n",
    "\n",
    "def kennzahlen(w, mu):\n",
    "    ertrag = float(mu @ w)\n",
    "    risiko = float(np.sqrt(w @ SIGMA @ w))\n",
    "    return ertrag, risiko\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    print(\"=\" * 88)\n",
    "    print(\"  ROBUSTE vs. NOMINALE PORTFOLIO-OPTIMIERUNG\")\n",
    "    print(\"=\" * 88)\n",
    "    print(f\"{'Asset':<14} {'Erw. Rendite':>14} {'Unsicherheit':>14} \"\n",
    "          f\"{'Worst Case':>12} {'Volatilitaet':>13}\")\n",
    "    print(\"-\" * 88)\n",
    "    for i, name in enumerate(ASSETS):\n",
    "        print(f\"{name:<14} {MU_SCHAETZUNG[i]*100:>13.1f} % \"\n",
    "              f\"{'+/- ' + format(UNSICHERHEIT[i]*100, '.1f') + ' %':>14} \"\n",
    "              f\"{(MU_SCHAETZUNG[i]-UNSICHERHEIT[i])*100:>11.1f} % \"\n",
    "              f\"{VOLA[i]*100:>12.1f} %\")\n",
    "\n",
    "    # (a) nominal: vertraut den Schaetzungen\n",
    "    w_nominal = optimiere(MU_SCHAETZUNG)\n",
    "    # (b) robust: rechnet mit dem Worst Case der Box-Unsicherheitsmenge\n",
    "    w_robust = optimiere(MU_SCHAETZUNG - UNSICHERHEIT)\n",
    "\n",
    "    print(\"\\n\" + \"-\" * 88)\n",
    "    print(f\"{'':<14} {'nominal':>22} {'robust':>22}\")\n",
    "    print(\"-\" * 88)\n",
    "    for i, name in enumerate(ASSETS):\n",
    "        print(f\"{name:<14} {w_nominal[i]*100:>21.1f} % {w_robust[i]*100:>21.1f} %\")\n",
    "\n",
    "    # --- Bewertung in beiden Welten --------------------------------------\n",
    "    mu_worst = MU_SCHAETZUNG - UNSICHERHEIT\n",
    "    e_nom_gut, r_nom = kennzahlen(w_nominal, MU_SCHAETZUNG)\n",
    "    e_rob_gut, r_rob = kennzahlen(w_robust, MU_SCHAETZUNG)\n",
    "    e_nom_schlecht, _ = kennzahlen(w_nominal, mu_worst)\n",
    "    e_rob_schlecht, _ = kennzahlen(w_robust, mu_worst)\n",
    "\n",
    "    print(\"\\n\" + \"-\" * 88)\n",
    "    print(f\"{'Bewertung':<34} {'nominales Portfolio':>22} {'robustes Portfolio':>22}\")\n",
    "    print(\"-\" * 88)\n",
    "    print(f\"{'Ertrag, wenn Schaetzung stimmt':<34} {e_nom_gut*100:>21.2f} % \"\n",
    "          f\"{e_rob_gut*100:>21.2f} %\")\n",
    "    print(f\"{'Ertrag im Worst Case':<34} {e_nom_schlecht*100:>21.2f} % \"\n",
    "          f\"{e_rob_schlecht*100:>21.2f} %\")\n",
    "    print(f\"{'Volatilitaet':<34} {r_nom*100:>21.2f} % {r_rob*100:>21.2f} %\")\n",
    "\n",
    "    print(\"\\n\" + \"-\" * 88)\n",
    "    print(f\"Preis der Robustheit (Ertragsverzicht im Normalfall): \"\n",
    "          f\"{(e_nom_gut - e_rob_gut)*100:+.2f} Prozentpunkte\")\n",
    "    print(f\"Nutzen der Robustheit (Vorteil im Worst Case):        \"\n",
    "          f\"{(e_rob_schlecht - e_nom_schlecht)*100:+.2f} Prozentpunkte\")\n",
    "    verhaeltnis = ((e_rob_schlecht - e_nom_schlecht)\n",
    "                   / max(e_nom_gut - e_rob_gut, 1e-9))\n",
    "    print(f\"Verhaeltnis Nutzen/Preis: {verhaeltnis:.2f}\")\n",
    "    print(\"  -> Werte > 1 bedeuten: Die Absicherung bringt im Ernstfall mehr,\")\n",
    "    print(\"     als sie im Normalfall kostet.\")\n",
    "    print(\"=\" * 88)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Der Fall: ein Kraftwerkspark mit 500 MW Zusage\n",
    "\n",
    "`Chance_Constraints.py`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "\n",
    "# Chance_Constraints.py\n",
    "\"\"\"\n",
    "Kapitel Unsicherheit: Wahrscheinlichkeitsbeschraenkungen.\n",
    "\n",
    "Das Management fragt selten \"was ist im Mittel am besten?\", sondern \"mit\n",
    "welcher Sicherheit haelt der Plan?\". Genau das formuliert eine Chance\n",
    "Constraint:  P(Versorgung >= Bedarf) >= 1 - alpha.\n",
    "\n",
    "Gezeigt werden beide Wege dorthin, an derselben Instanz:\n",
    "  (a) analytisch      - unter Normalverteilungsannahme wird daraus eine\n",
    "                        Second-Order-Cone-Bedingung, loesbar mit CVXPY\n",
    "  (b) szenariobasiert - Big-M mit Binaervariablen, fuer beliebige\n",
    "                        empirische Verteilungen\n",
    "\n",
    "Und die beiden Messungen, auf die es ankommt: Was kostet ein Prozentpunkt\n",
    "Versorgungssicherheit - und haelt die Zusage auch dann, wenn die Verteilung\n",
    "nicht normal ist?\n",
    "\n",
    "Solver: CLARABEL (Kegel) und SciPy/HiGHS (gemischt-ganzzahlig). Weder ortools\n",
    "noch ein direkter highspy-Import, damit alles in einem Prozess laeuft.\n",
    "\"\"\"\n",
    "\n",
    "import cvxpy as cp\n",
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "\n",
    "# --- Der Kraftwerkspark ----------------------------------------------------\n",
    "# Die Verfuegbarkeit ist der Anteil, den eine installierte MW im Mittel\n",
    "# wirklich liefert: bei Wind und Sonne klein und stark schwankend, bei Gas und\n",
    "# Biomasse gross und stabil. Die Ausbaugrenze ist der Standort - Flaeche,\n",
    "# Genehmigung, Brennstoffversorgung.\n",
    "TECHNIK = [\"Gaskraftwerk\", \"Windpark\", \"Solarpark\", \"Biomasse\"]\n",
    "KOSTEN = np.array([65_000.0, 24_000.0, 12_000.0, 60_000.0])   # EUR je MW und Jahr\n",
    "VERFUEGBAR = np.array([0.92, 0.35, 0.18, 0.85])\n",
    "STREUUNG = np.array([0.05, 0.16, 0.10, 0.04])\n",
    "GRENZE = np.array([400.0, 250.0, 600.0, 350.0])               # MW\n",
    "\n",
    "# Wind und Sonne sind leicht gegenlaeufig: Ein Tiefdruckgebiet bringt Wind und\n",
    "# Wolken zugleich. Genau diese Korrelation macht die Bedingung zu einem Kegel\n",
    "# und nicht zu einer Summe unabhaengiger Einzelzuschlaege.\n",
    "KORRELATION = np.array([\n",
    "    [1.00,  0.00,  0.00, 0.05],\n",
    "    [0.00,  1.00, -0.25, 0.00],\n",
    "    [0.00, -0.25,  1.00, 0.00],\n",
    "    [0.05,  0.00,  0.00, 1.00],\n",
    "])\n",
    "SIGMA = np.diag(STREUUNG) @ KORRELATION @ np.diag(STREUUNG)\n",
    "WURZEL = np.linalg.cholesky(SIGMA)          # L mit L @ L.T == SIGMA\n",
    "\n",
    "BEDARF = 500.0          # MW, die gesichert bereitstehen muessen\n",
    "N = len(TECHNIK)\n",
    "\n",
    "# Die Kaeltewelle mit Dunkelflaute: selten, aber sie trifft alles zugleich.\n",
    "# Wind und Sonne brechen fast vollstaendig weg - und, das ist der Punkt, das\n",
    "# Gaskraftwerk liefert ebenfalls weniger, weil bei Frost der Netzdruck faellt.\n",
    "# Eine Kovarianzmatrix mit Korrelationen um null kann das nicht ausdruecken.\n",
    "P_KAELTEWELLE = 0.08\n",
    "EINBRUCH = np.array([0.72, 0.10, 0.15, 0.88])\n",
    "\n",
    "\n",
    "def mittelwertplan():\n",
    "    \"\"\"Plant mit den Erwartungswerten - ignoriert die Streuung vollstaendig.\"\"\"\n",
    "    x = cp.Variable(N, nonneg=True)\n",
    "    problem = cp.Problem(cp.Minimize(KOSTEN @ x),\n",
    "                         [VERFUEGBAR @ x >= BEDARF, x <= GRENZE])\n",
    "    problem.solve(solver=cp.CLARABEL)\n",
    "    return problem.value, x.value\n",
    "\n",
    "\n",
    "def chance_constraint_analytisch(alpha):\n",
    "    \"\"\"\n",
    "    P(a^T x >= BEDARF) >= 1 - alpha  unter a ~ N(VERFUEGBAR, SIGMA).\n",
    "\n",
    "    Aequivalent zu:  VERFUEGBAR^T x - z * ||L^T x||_2 >= BEDARF\n",
    "    mit z = Phi^-1(1 - alpha). Das ist eine Second-Order-Cone-Bedingung: Der\n",
    "    Sicherheitszuschlag ist eine Norm ueber x, kein fester Aufschlag je Anlage.\n",
    "    Deshalb belohnt sie Mischung - zwei gegenlaeufige Quellen schwanken\n",
    "    gemeinsam weniger als jede fuer sich.\n",
    "    \"\"\"\n",
    "    x = cp.Variable(N, nonneg=True)\n",
    "    z = norm.ppf(1 - alpha)\n",
    "    bedingung = VERFUEGBAR @ x - z * cp.norm(WURZEL.T @ x, 2) >= BEDARF\n",
    "    problem = cp.Problem(cp.Minimize(KOSTEN @ x), [bedingung, x <= GRENZE])\n",
    "    problem.solve(solver=cp.CLARABEL)\n",
    "    return problem.value, x.value\n",
    "\n",
    "\n",
    "def chance_constraint_szenarien(a, alpha):\n",
    "    \"\"\"\n",
    "    Dieselbe Zusage ohne Verteilungsannahme: Fuer jedes Szenario s sagt eine\n",
    "    Binaervariable z_s, ob es verletzt werden darf.\n",
    "\n",
    "        a_s^T x >= BEDARF - M * z_s        fuer alle s\n",
    "        sum_s z_s <= alpha * S\n",
    "\n",
    "    M = BEDARF ist die kleinstmoegliche gueltige Schranke, denn a_s^T x >= 0:\n",
    "    Groesser kann eine Verletzung gar nicht ausfallen. Ein unnoetig grosses M\n",
    "    wuerde die LP-Relaxierung aufweichen und die Suche verlangsamen - die\n",
    "    Big-M-Falle aus dem Kapitel Gemischt-ganzzahlige Optimierung.\n",
    "    \"\"\"\n",
    "    S = len(a)\n",
    "    x = cp.Variable(N, nonneg=True)\n",
    "    z = cp.Variable(S, boolean=True)\n",
    "    problem = cp.Problem(\n",
    "        cp.Minimize(KOSTEN @ x),\n",
    "        [a @ x >= BEDARF - BEDARF * z, cp.sum(z) <= alpha * S, x <= GRENZE])\n",
    "    problem.solve(solver=cp.SCIPY)\n",
    "    return problem.value, x.value\n",
    "\n",
    "\n",
    "def ziehe_wetter(anzahl, seed):\n",
    "    \"\"\"Mischverteilung: normales Wetter, mit P_KAELTEWELLE ein Einbruch.\"\"\"\n",
    "    rng = np.random.default_rng(seed)\n",
    "    a = rng.multivariate_normal(VERFUEGBAR, SIGMA, size=anzahl)\n",
    "    getroffen = rng.random(anzahl) < P_KAELTEWELLE\n",
    "    a[getroffen] *= EINBRUCH\n",
    "    return np.clip(a, 0.0, 1.0)\n",
    "\n",
    "\n",
    "def sicherheit(x, a):\n",
    "    \"\"\"Anteil der Szenarien, in denen der Plan den Bedarf deckt.\"\"\"\n",
    "    return float(np.mean(a @ x >= BEDARF))\n",
    "\n",
    "\n",
    "KOPF = (f\"{'Plan':<20} {'Kosten (EUR)':>13}  {'Gas':>4} {'Wind':>4} \"\n",
    "        f\"{'Sol':>4} {'Bio':>4}  {'Normalwelt':>11} {'echte Welt':>11}\")\n",
    "\n",
    "\n",
    "def zeile(name, kosten, x, in_normal, in_echt):\n",
    "    mix = \" \".join(f\"{w:4.0f}\" for w in x)\n",
    "    return (f\"{name:<20} {kosten:>13,.0f}  {mix}  \"\n",
    "            f\"{in_normal * 100:>9.2f} % {in_echt * 100:>9.2f} %\")\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    print(\"=\" * 78)\n",
    "    print(\"  WAHRSCHEINLICHKEITSBESCHRAENKUNGEN IM KRAFTWERKSPARK\")\n",
    "    print(\"=\" * 78)\n",
    "    print(f\"Gesichert bereitzustellen: {BEDARF:.0f} MW\\n\")\n",
    "    print(f\"{'Technologie':<14} {'EUR/MW':>8} {'verfuegbar':>11} {'Streuung':>9} \"\n",
    "          f\"{'Grenze':>8} {'EUR je erw. MW':>15}\")\n",
    "    print(\"-\" * 78)\n",
    "    for i, name in enumerate(TECHNIK):\n",
    "        print(f\"{name:<14} {KOSTEN[i]:>8,.0f} {VERFUEGBAR[i]:>10.2f}  \"\n",
    "              f\"{STREUUNG[i]:>8.2f} {GRENZE[i]:>7.0f} \"\n",
    "              f\"{KOSTEN[i] / VERFUEGBAR[i]:>15,.0f}\")\n",
    "    print(f\"\\nVollausbau liefert im Mittel {VERFUEGBAR @ GRENZE:.0f} MW, \"\n",
    "          f\"in der Kaeltewelle {(VERFUEGBAR * EINBRUCH) @ GRENZE:.0f} MW.\")\n",
    "\n",
    "    # Zwei grosse, unabhaengige Testmengen. An ihnen wird JEDER Plan gemessen -\n",
    "    # einmal in der unterstellten Normalwelt, einmal in der echten Verteilung.\n",
    "    echt = ziehe_wetter(200_000, seed=771)\n",
    "    normalwelt = np.random.default_rng(4711).multivariate_normal(\n",
    "        VERFUEGBAR, SIGMA, size=200_000)\n",
    "\n",
    "    print(\"\\n\" + \"=\" * 78)\n",
    "    print(\"  (1) Mittelwertplan und Chance Constraints im Vergleich\")\n",
    "    print(\"=\" * 78)\n",
    "    print(KOPF)\n",
    "    print(\"-\" * 78)\n",
    "\n",
    "    k_mittel, x_mittel = mittelwertplan()\n",
    "    stufen = [(\"Mittelwert\", k_mittel, sicherheit(x_mittel, normalwelt))]\n",
    "    print(zeile(\"Mittelwert\", k_mittel, x_mittel,\n",
    "                sicherheit(x_mittel, normalwelt), sicherheit(x_mittel, echt)))\n",
    "\n",
    "    plaene = {}\n",
    "    for alpha in (0.20, 0.10, 0.05, 0.01):\n",
    "        kosten, x = chance_constraint_analytisch(alpha)\n",
    "        plaene[alpha] = (kosten, x)\n",
    "        quote_normal = sicherheit(x, normalwelt)\n",
    "        stufen.append((f\"Zusage {(1 - alpha) * 100:.0f} %\", kosten, quote_normal))\n",
    "        print(zeile(f\"Zusage {(1 - alpha) * 100:.0f} %\", kosten, x,\n",
    "                    quote_normal, sicherheit(x, echt)))\n",
    "\n",
    "    print(\"\\n  In der Normalwelt trifft jede Zusage ihren Wert - das Verfahren\")\n",
    "    print(\"  rechnet richtig. Die Spalte 'echte Welt' kommt in Teil (3).\")\n",
    "\n",
    "    print(\"\\n\" + \"=\" * 78)\n",
    "    print(\"  (2) Was kostet ein Prozentpunkt Versorgungssicherheit?\")\n",
    "    print(\"=\" * 78)\n",
    "    print(f\"{'von -> nach':<26} {'Prozentpunkte':>13} {'Mehrkosten':>13} \"\n",
    "          f\"{'EUR je Punkt':>13}\")\n",
    "    print(\"-\" * 78)\n",
    "    for (n1, k1, q1), (n2, k2, q2) in zip(stufen, stufen[1:]):\n",
    "        d_punkte = (q2 - q1) * 100\n",
    "        d_kosten = k2 - k1\n",
    "        print(f\"{n1 + ' -> ' + n2:<26} {d_punkte:>13.2f} {d_kosten:>13,.0f} \"\n",
    "              f\"{d_kosten / d_punkte:>13,.0f}\")\n",
    "    erst = (stufen[1][1] - stufen[0][1]) / ((stufen[1][2] - stufen[0][2]) * 100)\n",
    "    letzt = (stufen[-1][1] - stufen[-2][1]) / ((stufen[-1][2] - stufen[-2][2]) * 100)\n",
    "    print(f\"\\n  Der letzte Prozentpunkt kostet das {letzt / erst:.1f}-fache des ersten.\")\n",
    "    print(\"  Sicherheit ist konvex bepreist - genau deshalb muss jemand\")\n",
    "    print(\"  entscheiden, wie viel davon das Unternehmen kaufen will.\")\n",
    "\n",
    "    print(\"\\n\" + \"=\" * 78)\n",
    "    print(\"  (3) Und wenn die Verteilung nicht normal ist?\")\n",
    "    print(\"=\" * 78)\n",
    "    k_soc, x_soc = plaene[0.05]\n",
    "    print(f\"Die echte Welt kennt die Kaeltewelle mit Dunkelflaute: in \"\n",
    "          f\"{P_KAELTEWELLE * 100:.0f} % der Faelle liefern\")\n",
    "    print(f\"Wind {EINBRUCH[1] * 100:.0f} %, Sonne {EINBRUCH[2] * 100:.0f} % \"\n",
    "          f\"und - entscheidend - auch das Gaskraftwerk nur \"\n",
    "          f\"{EINBRUCH[0] * 100:.0f} %\")\n",
    "    print(\"des Ueblichen. Alles bricht gleichzeitig ein.\\n\")\n",
    "\n",
    "    a_bau = ziehe_wetter(400, seed=20260908)\n",
    "    k_sz, x_sz = chance_constraint_szenarien(a_bau, 0.05)\n",
    "\n",
    "    print(KOPF)\n",
    "    print(\"-\" * 78)\n",
    "    print(zeile(\"SOC, Zusage 95 %\", k_soc, x_soc,\n",
    "                sicherheit(x_soc, normalwelt), sicherheit(x_soc, echt)))\n",
    "    print(zeile(\"Szenarien, 95 %\", k_sz, x_sz,\n",
    "                sicherheit(x_sz, normalwelt), sicherheit(x_sz, echt)))\n",
    "\n",
    "    print(f\"\\n  Der SOC-Plan verspricht 95 % und haelt \"\n",
    "          f\"{sicherheit(x_soc, echt) * 100:.1f} %. Nicht das Verfahren ist\")\n",
    "    print(f\"  falsch, sondern die Annahme: In der Normalwelt liefert derselbe\")\n",
    "    print(f\"  Plan {sicherheit(x_soc, normalwelt) * 100:.1f} %.\")\n",
    "    print(f\"\\n  Der Szenarioplan ({len(a_bau)} Szenarien, {len(a_bau)} \"\n",
    "          f\"Binaervariablen) kommt auf\")\n",
    "    print(f\"  {sicherheit(x_sz, echt) * 100:.1f} % - und kostet dafuer \"\n",
    "          f\"{(k_sz / k_soc - 1) * 100:+.1f} %. Er kauft keine Windkraft mehr:\")\n",
    "    print(\"  Was im Ernstfall ausfaellt, hilft der Zusage nicht.\")\n",
    "    print(\"=\" * 78)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}