operations_research/OR_HTML_04/programme/Stochastische_Optimierung.py

#!/usr/bin/env python3

# Stochastische_Optimierung.py
"""
Kapitel Unsicherheit: Two-Stage Stochastic Programming mit CVXPY.

Eigenschaften:
  * Die Analyse am Ende wird BERECHNET statt fest verdrahtet
  * Vergleich gegen drei Alternativen: Mittelwert, Worst Case, perfekte Voraussicht
  * Kennzahl EVPI (Wert perfekter Information) wird ausgewiesen
"""

import cvxpy as cp
import numpy as np
import pandas as pd

SZENARIEN = ["Ruhig", "Volatil", "Crash"]
WAHRSCHEINLICHKEIT = np.array([0.50, 0.30, 0.20])
BEDARF = np.array([100.0, 250.0, 500.0])

KOSTEN_VORAB = 40.0
KOSTEN_SPOT = 120.0
KOSTEN_LEERLAUF = 5.0

S = len(SZENARIEN)


def loese_stochastisch():
    """Zweistufiges Modell: eine Vorabentscheidung, szenarioabhaengige Korrektur."""
    x = cp.Variable(nonneg=True, name="Basiskapazitaet")        # Stufe 1
    y_spot = cp.Variable(S, nonneg=True, name="Spot_Zukauf")    # Stufe 2
    y_leer = cp.Variable(S, nonneg=True, name="Leerlauf")       # Stufe 2

    # Kopplung: Basis + Zukauf - Leerlauf == Bedarf (je Szenario)
    nebenbedingungen = [x + y_spot[s] - y_leer[s] == BEDARF[s] for s in range(S)]

    erwartete_korrektur = sum(
        WAHRSCHEINLICHKEIT[s] * (KOSTEN_SPOT * y_spot[s] + KOSTEN_LEERLAUF * y_leer[s])
        for s in range(S))

    problem = cp.Problem(cp.Minimize(KOSTEN_VORAB * x + erwartete_korrektur),
                         nebenbedingungen)
    problem.solve()
    return problem, x, y_spot, y_leer


def kosten_bei(kapazitaet):
    """Erwartete Gesamtkosten fuer eine fest vorgegebene Kapazitaet."""
    unter = np.maximum(BEDARF - kapazitaet, 0.0)
    ueber = np.maximum(kapazitaet - BEDARF, 0.0)
    je_szenario = KOSTEN_VORAB * kapazitaet + KOSTEN_SPOT * unter + KOSTEN_LEERLAUF * ueber
    return float(WAHRSCHEINLICHKEIT @ je_szenario)


if __name__ == "__main__":
    problem, x, y_spot, y_leer = loese_stochastisch()
    kapazitaet = float(x.value)
    mittelwert = float(WAHRSCHEINLICHKEIT @ BEDARF)

    print("=" * 82)
    print("     STOCHASTISCHE TWO-STAGE OPTIMIERUNG (KAPAZITAETSPLANUNG)")
    print("=" * 82)
    print(f"Status: {problem.status}")
    print(f"Erwarteter Bedarf (Mittelwert):   {mittelwert:.1f} Einheiten")
    print(f"Optimale Stufe-1-Kapazitaet x*:   {kapazitaet:.1f} Einheiten")
    print(f"Minimale erwartete Gesamtkosten:  {problem.value:,.2f} EUR\n")

    tabelle = pd.DataFrame({
        "Szenario": SZENARIEN,
        "Wahrsch.": [f"{p*100:.0f} %" for p in WAHRSCHEINLICHKEIT],
        "Bedarf": BEDARF,
        "Basis genutzt": [min(kapazitaet, b) for b in BEDARF],
        "Spot-Zukauf": np.round(y_spot.value, 1),
        "Leerlauf": np.round(y_leer.value, 1),
        "Kosten (EUR)": [f"{KOSTEN_VORAB*kapazitaet + KOSTEN_SPOT*y_spot.value[s] + KOSTEN_LEERLAUF*y_leer.value[s]:,.0f}"
                         for s in range(S)],
    })
    print(tabelle.to_string(index=False))

    # --- Vergleich mit Alternativstrategien (berechnet, nicht behauptet) --
    print("\n" + "-" * 82)
    print("Vergleich verschiedener Planungsstrategien:")
    print(f"{'Strategie':<34} {'Kapazitaet':>11} {'Erw. Kosten':>14} {'Mehrkosten':>13}")
    print("-" * 82)

    optimal = problem.value
    strategien = [
        ("Stochastisch optimal", kapazitaet),
        ("Naiv: Mittelwert einsetzen", mittelwert),
        ("Vorsichtig: Worst Case abdecken", float(BEDARF.max())),
        ("Optimistisch: Bestfall", float(BEDARF.min())),
    ]
    for name, kap in strategien:
        kosten = kosten_bei(kap)
        print(f"{name:<34} {kap:>11.1f} {kosten:>14,.0f} "
              f"{kosten - optimal:>+13,.0f}")

    # --- EVPI: Was waere perfekte Voraussicht wert? ----------------------
    # Bei perfekter Information wuerde man je Szenario genau den Bedarf kaufen.
    kosten_perfekt = float(WAHRSCHEINLICHKEIT @ (KOSTEN_VORAB * BEDARF))
    evpi = optimal - kosten_perfekt
    print("-" * 82)
    print(f"Kosten bei perfekter Voraussicht:      {kosten_perfekt:>10,.0f} EUR")
    print(f"Wert perfekter Information (EVPI):     {evpi:>10,.0f} EUR "
          f"({evpi/optimal*100:.1f} % der Kosten)")
    print("  -> So viel duerfte eine perfekte Bedarfsprognose hoechstens kosten.")

    # --- Automatische Interpretation --------------------------------------
    print("-" * 82)
    if kapazitaet > mittelwert + 1e-6:
        print(f"Analyse: Der Solver waehlt {kapazitaet:.0f} Einheiten und damit MEHR als")
        print(f"den Mittelwert ({mittelwert:.0f}), weil Unterdeckung ({KOSTEN_SPOT:.0f} EUR)")
        print(f"deutlich teurer ist als Leerlauf ({KOSTEN_LEERLAUF:.0f} EUR).")
    elif kapazitaet < mittelwert - 1e-6:
        print(f"Analyse: Der Solver waehlt {kapazitaet:.0f} und damit WENIGER als den")
        print(f"Mittelwert ({mittelwert:.0f}) - Leerlauf ist hier teurer als Zukauf.")
    else:
        print("Analyse: Kapazitaet entspricht dem Mittelwert (symmetrische Kosten).")
    print("=" * 82)