operations_research/Notebooks_04/qp-nlp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kapitel 11: Quadratische und nichtlineare Optimierung — KKT, Lagrange, Konvexität\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": [
    "## Konvexität: die präzise Aussage\n",
    "\n",
    "`QP_Grundlagen.py`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "\n",
    "# QP_Grundlagen.py\n",
    "\"\"\"\n",
    "Kapitel QP/NLP: Die drei Faelle aus der Konvexitaets-Tabelle (Abschnitt\n",
    "'Das quadratische Programm') an\n",
    "einem Mini-QP demonstriert: P positiv definit, P (singulaer) semidefinit,\n",
    "P mit negativem Eigenwert.\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "import cvxpy as cp\n",
    "\n",
    "\n",
    "def loese_qp(P, q, name):\n",
    "    n = len(q)\n",
    "    w = cp.Variable(n)\n",
    "    ziel = cp.Minimize(0.5 * cp.quad_form(w, P) + q @ w)\n",
    "    bedingungen = [cp.sum(w) == 1, w >= 0]\n",
    "    problem = cp.Problem(ziel, bedingungen)\n",
    "\n",
    "    print(f\"\\n--- {name} ---\")\n",
    "    eigenwerte = np.linalg.eigvalsh(P)\n",
    "    print(f\"Eigenwerte von P: {np.round(eigenwerte, 4)}\")\n",
    "    print(f\"DCP-konvex (CVXPY-Pruefung)? {problem.is_dcp()}\")\n",
    "\n",
    "    if not problem.is_dcp():\n",
    "        print(\"-> CVXPY lehnt das Problem ab, BEVOR ueberhaupt ein Solver laeuft.\")\n",
    "        return\n",
    "\n",
    "    problem.solve()\n",
    "    print(f\"Status: {problem.status}\")\n",
    "    print(f\"w* = {np.round(w.value, 4)}\")\n",
    "    print(f\"Zielwert = {problem.value:.6f}\")\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    q = np.zeros(2)\n",
    "\n",
    "    # Fall 1: P positiv definit -> eindeutiges Minimum\n",
    "    P_definit = np.array([[2.0, 0.5], [0.5, 1.0]])\n",
    "    loese_qp(P_definit, q, \"P positiv definit\")\n",
    "\n",
    "    # Fall 2: P singulaer/semidefinit (zwei \"identische\" Assets) -> unendlich viele Minima\n",
    "    P_semidefinit = np.array([[1.0, 1.0], [1.0, 1.0]])\n",
    "    loese_qp(P_semidefinit, q, \"P positiv semidefinit (singulaer)\")\n",
    "\n",
    "    # Fall 3: P mit negativem Eigenwert -> nicht konvex\n",
    "    P_indefinit = np.array([[1.0, 2.0], [2.0, 1.0]])\n",
    "    loese_qp(P_indefinit, q, \"P indefinit (negativer Eigenwert)\")\n",
    "\n",
    "    print(\"\\n--- Nachweis: 'unendlich viele Minima' im semidefiniten Fall ---\")\n",
    "    for punkt in [np.array([1.0, 0.0]), np.array([0.0, 1.0]), np.array([0.3, 0.7])]:\n",
    "        wert = 0.5 * punkt @ P_semidefinit @ punkt\n",
    "        print(f\"  w = {punkt} -> Zielwert = {wert:.4f}  (identisch, obwohl w verschieden)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Die vier KKT-Bedingungen\n",
    "\n",
    "`KKT_Nachweis.py`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "\n",
    "# KKT_Nachweis.py\n",
    "\"\"\"\n",
    "Kapitel QP/NLP: Die KKT-Bedingungen numerisch nachpruefen.\n",
    "\n",
    "Loest ein QP mit CVXPY, liest die Dualwerte aus und prueft alle vier\n",
    "KKT-Bedingungen einzeln nach. Das ist zugleich eine Vorlage fuer die\n",
    "Qualitaetssicherung eigener Modelle.\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "import cvxpy as cp\n",
    "\n",
    "\n",
    "def baue_gueltige_kovarianz(vola, korrelationen, seed=0):\n",
    "    \"\"\"\n",
    "    Baut Sigma = D * C * D aus Volatilitaeten und einer Korrelationsmatrix.\n",
    "    Dieses Vorgehen ist konstruktionsbedingt positiv semidefinit - im\n",
    "    Gegensatz zum nachtraeglichen Ueberschreiben der Diagonalen, das die\n",
    "    positive Semidefinitheit zerstoeren kann.\n",
    "    \"\"\"\n",
    "    C = np.array(korrelationen, dtype=float)\n",
    "    assert np.allclose(C, C.T), \"Korrelationsmatrix muss symmetrisch sein.\"\n",
    "    assert np.allclose(np.diag(C), 1.0), \"Diagonale der Korrelationsmatrix muss 1 sein.\"\n",
    "    eigen = np.linalg.eigvalsh(C)\n",
    "    assert eigen.min() > -1e-10, (\n",
    "        f\"Korrelationsmatrix ist nicht positiv semidefinit \"\n",
    "        f\"(kleinster Eigenwert {eigen.min():.4f}). Solche Korrelationen sind unmoeglich.\")\n",
    "    D = np.diag(vola)\n",
    "    return D @ C @ D\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    # --- Ein kleines Portfolio-QP ----------------------------------------\n",
    "    vola = np.array([0.20, 0.14, 0.30])                 # Volatilitaeten\n",
    "    korr = [[1.00, 0.30, 0.10],\n",
    "            [0.30, 1.00, 0.25],\n",
    "            [0.10, 0.25, 1.00]]\n",
    "    Sigma = baue_gueltige_kovarianz(vola, korr)\n",
    "    mu = np.array([0.09, 0.05, 0.13])                   # erwartete Renditen\n",
    "    lam = 3.0                                           # Risikoaversion\n",
    "\n",
    "    eigenwerte = np.linalg.eigvalsh(Sigma)\n",
    "    print(\"=\" * 74)\n",
    "    print(\"  KKT-BEDINGUNGEN AM PORTFOLIO-QP\")\n",
    "    print(\"=\" * 74)\n",
    "    print(f\"Eigenwerte von Sigma: {np.round(eigenwerte, 6)}\")\n",
    "    print(f\"  -> positiv definit: {bool(eigenwerte.min() > 0)}  \"\n",
    "          f\"(Problem ist streng konvex, Loesung eindeutig)\")\n",
    "    print(f\"  -> Konditionszahl:  {eigenwerte.max() / eigenwerte.min():.2f}\")\n",
    "\n",
    "    # --- Modell: min  lam/2 * w'Sigma w - mu'w   u.d.N. sum(w)=1, w>=0 ----\n",
    "    n = len(mu)\n",
    "    w = cp.Variable(n)\n",
    "    ziel = cp.Minimize(0.5 * lam * cp.quad_form(w, Sigma) - mu @ w)\n",
    "    budget = cp.sum(w) == 1\n",
    "    nichtnegativ = w >= 0\n",
    "    problem = cp.Problem(ziel, [budget, nichtnegativ])\n",
    "    problem.solve()\n",
    "\n",
    "    w_opt = w.value\n",
    "    nu = budget.dual_value                    # Multiplikator der Gleichung\n",
    "    lam_i = nichtnegativ.dual_value           # Multiplikatoren der Ungleichungen\n",
    "\n",
    "    print(f\"\\nStatus: {problem.status}\")\n",
    "    print(f\"Optimale Gewichte: {np.round(w_opt, 6)}\")\n",
    "    print(f\"Zielwert:          {problem.value:.6f}\")\n",
    "    print(f\"Multiplikator der Budgetgleichung (nu): {nu:.6f}\")\n",
    "    print(f\"Multiplikatoren der w>=0-Bedingungen:   {np.round(lam_i, 6)}\")\n",
    "\n",
    "    # --- KKT-Bedingungen einzeln pruefen ---------------------------------\n",
    "    print(\"\\n--- Pruefung der vier KKT-Bedingungen ---\")\n",
    "\n",
    "    # 1. Stationaritaet:  grad f - lambda + nu*1 = 0\n",
    "    #    f(w) = lam/2 w'Sigma w - mu'w   ->   grad f = lam*Sigma w - mu\n",
    "    #    g_i(w) = -w_i <= 0              ->   grad g_i = -e_i\n",
    "    #    h(w)   = sum(w) - 1 = 0         ->   grad h   = 1\n",
    "    grad_f = lam * (Sigma @ w_opt) - mu\n",
    "    stationaritaet = grad_f - lam_i + nu * np.ones(n)\n",
    "    print(f\"1. Stationaritaet   : max|Residuum| = {np.abs(stationaritaet).max():.2e}\")\n",
    "\n",
    "    # 2. Primale Zulaessigkeit\n",
    "    print(f\"2. Primal zulaessig : sum(w)-1 = {w_opt.sum()-1:.2e}, \"\n",
    "          f\"min(w) = {w_opt.min():.2e}\")\n",
    "\n",
    "    # 3. Duale Zulaessigkeit\n",
    "    print(f\"3. Dual zulaessig   : min(lambda) = {lam_i.min():.2e}  (muss >= 0 sein)\")\n",
    "\n",
    "    # 4. Komplementaerer Schlupf: lambda_i * w_i = 0\n",
    "    print(f\"4. Kompl. Schlupf   : max|lambda_i * w_i| = \"\n",
    "          f\"{np.abs(lam_i * w_opt).max():.2e}\")\n",
    "\n",
    "    alle_ok = (np.abs(stationaritaet).max() < 1e-6\n",
    "               and abs(w_opt.sum() - 1) < 1e-8\n",
    "               and w_opt.min() > -1e-8\n",
    "               and lam_i.min() > -1e-8\n",
    "               and np.abs(lam_i * w_opt).max() < 1e-6)\n",
    "    print(f\"\\nAlle vier KKT-Bedingungen erfuellt: {alle_ok}\")\n",
    "\n",
    "    # --- Interpretation von nu -------------------------------------------\n",
    "    print(\"\\n--- Was bedeutet nu? ---\")\n",
    "    print(\"nu ist der Schattenpreis des Budgets: Um so viel aendert sich der\")\n",
    "    print(\"Zielwert, wenn man statt 100 % nur 99 % investieren duerfte.\")\n",
    "    problem2 = cp.Problem(cp.Minimize(0.5 * lam * cp.quad_form(w, Sigma) - mu @ w),\n",
    "                          [cp.sum(w) == 1.01, w >= 0])\n",
    "    problem2.solve()\n",
    "    print(f\"  Vorhergesagt (nu * 0.01): {nu * 0.01:+.6f}\")\n",
    "    print(f\"  Tatsaechlich gemessen:    {problem2.value - problem.value:+.6f}\")\n",
    "    print(\"=\" * 74)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Praxisfall: Entropie-maximierte Kapitalallokation\n",
    "\n",
    "`Entropie_Maximierte_Allokation.py`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "\n",
    "# Entropie_Maximierte_Allokation.py\n",
    "\"\"\"\n",
    "Kapitel QP/NLP: Nichtlineare Optimierung (NLP) mit scipy SLSQP.\n",
    "Modell: Mean-Variance-Portfolio mit Shannon-Entropie-Diversifikation.\n",
    "\n",
    "Achtung, haeufiger Fehler: Sigma per A@A.T zu erzeugen und anschliessend die\n",
    "Diagonale zu ueberschreiben zerstoert die positive Semidefinitheit - die\n",
    "Matrix kann dadurch einen negativen Eigenwert bekommen und unmoegliche\n",
    "Korrelationen (> 1) implizieren. Dieses Programm baut Sigma stattdessen als\n",
    "D * C * D aus Volatilitaeten und einer echten Korrelationsmatrix -\n",
    "konstruktionsbedingt immer PSD.\n",
    "\n",
    "Zusaetzlich: Multistart, weil SLSQP nur lokale Optima findet.\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "ASSETS = [\"Tech-Aktien\", \"Rohstoffe\", \"US-Treasuries\", \"Krypto\"]\n",
    "MU = np.array([0.14, 0.07, 0.03, 0.22])          # erwartete Jahresrenditen\n",
    "VOLA = np.array([0.20, 0.14, 0.07, 0.40])        # Volatilitaeten p.a.\n",
    "KORRELATION = np.array([\n",
    "    [1.00, 0.25, -0.10, 0.55],\n",
    "    [0.25, 1.00, 0.05, 0.20],\n",
    "    [-0.10, 0.05, 1.00, -0.15],\n",
    "    [0.55, 0.20, -0.15, 1.00],\n",
    "])\n",
    "\n",
    "ALPHA = 1.0        # Gewicht des Ertrags\n",
    "BETA = 0.015       # Gewicht der Entropie-Diversifikation\n",
    "UNTERGRENZE = 0.001\n",
    "OBERGRENZE = 0.80\n",
    "\n",
    "\n",
    "def baue_kovarianz() -> np.ndarray:\n",
    "    \"\"\"Sigma = D * C * D. Immer PSD, wenn C eine gueltige Korrelationsmatrix ist.\"\"\"\n",
    "    eig_c = np.linalg.eigvalsh(KORRELATION)\n",
    "    if eig_c.min() < -1e-10:\n",
    "        raise ValueError(f\"Korrelationsmatrix unmoeglich (Eigenwert {eig_c.min():.4f}).\")\n",
    "    D = np.diag(VOLA)\n",
    "    Sigma = D @ KORRELATION @ D\n",
    "    eig_s = np.linalg.eigvalsh(Sigma)\n",
    "    assert eig_s.min() > 0, \"Sigma nicht positiv definit!\"\n",
    "    return Sigma\n",
    "\n",
    "\n",
    "SIGMA = baue_kovarianz()\n",
    "\n",
    "\n",
    "def zielfunktion(w, alpha=ALPHA, beta=BETA):\n",
    "    \"\"\"Risiko - Ertrag - Entropiepraemie (wird minimiert).\"\"\"\n",
    "    varianz = 0.5 * w @ SIGMA @ w\n",
    "    ertrag = MU @ w\n",
    "    entropie = -np.sum(w * np.log(w + 1e-12))\n",
    "    return varianz - alpha * ertrag - beta * entropie\n",
    "\n",
    "\n",
    "def gradient(w, alpha=ALPHA, beta=BETA):\n",
    "    \"\"\"Exakter analytischer Gradient - beschleunigt und stabilisiert SLSQP.\"\"\"\n",
    "    grad_varianz = SIGMA @ w                       # d/dw von 0.5 w'Sigma w\n",
    "    grad_ertrag = -alpha * MU\n",
    "    grad_entropie = beta * (np.log(w + 1e-12) + 1.0)\n",
    "    return grad_varianz + grad_ertrag + grad_entropie\n",
    "\n",
    "\n",
    "def optimiere(startpunkt, alpha=ALPHA, beta=BETA):\n",
    "    \"\"\"alpha und beta werden durchgereicht - so bleibt die Funktion seiteneffektfrei.\"\"\"\n",
    "    return minimize(\n",
    "        zielfunktion, startpunkt, args=(alpha, beta), jac=gradient,\n",
    "        method=\"SLSQP\", bounds=[(UNTERGRENZE, OBERGRENZE)] * len(MU),\n",
    "        constraints=({\"type\": \"eq\",\n",
    "                      \"fun\": lambda w: np.sum(w) - 1.0,\n",
    "                      \"jac\": lambda w: np.ones(len(MU))}),\n",
    "        options={\"ftol\": 1e-12, \"maxiter\": 300})\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    n = len(MU)\n",
    "    eig = np.linalg.eigvalsh(SIGMA)\n",
    "\n",
    "    print(\"=\" * 74)\n",
    "    print(\"   NICHTLINEARE ENTROPIE-OPTIMIERTE ASSET-ALLOKATION\")\n",
    "    print(\"=\" * 74)\n",
    "    print(\"Pruefung der Kovarianzmatrix:\")\n",
    "    print(f\"  Eigenwerte:          {np.round(eig, 6)}\")\n",
    "    print(f\"  Positiv definit:     {bool(eig.min() > 0)}\")\n",
    "    print(f\"  Groesste Korrelation ausserhalb der Diagonale: \"\n",
    "          f\"{np.abs(KORRELATION - np.eye(n)).max():.2f}  (muss <= 1 sein)\")\n",
    "\n",
    "    # --- Multistart: SLSQP findet nur lokale Optima -----------------------\n",
    "    rng = np.random.default_rng(7)\n",
    "    startpunkte = [np.ones(n) / n]                                  # Gleichgewichtung\n",
    "    for _ in range(9):\n",
    "        z = rng.random(n) + 0.05\n",
    "        startpunkte.append(z / z.sum())\n",
    "\n",
    "    ergebnisse = [optimiere(s) for s in startpunkte]\n",
    "    erfolgreich = [r for r in ergebnisse if r.success]\n",
    "    bestes = min(erfolgreich, key=lambda r: r.fun)\n",
    "    zielwerte = np.array([r.fun for r in erfolgreich])\n",
    "\n",
    "    print(f\"\\nMultistart mit {len(startpunkte)} Startpunkten:\")\n",
    "    print(f\"  Erfolgreich konvergiert: {len(erfolgreich)}\")\n",
    "    print(f\"  Spannweite der Zielwerte: {zielwerte.max() - zielwerte.min():.2e}\")\n",
    "    print(\"  -> \" + (\"alle Startpunkte fuehren zum selben Optimum (Indiz fuer Konvexitaet)\"\n",
    "                     if zielwerte.max() - zielwerte.min() < 1e-6\n",
    "                     else \"ACHTUNG: verschiedene lokale Optima gefunden!\"))\n",
    "\n",
    "    w_opt = bestes.x\n",
    "    rendite = MU @ w_opt\n",
    "    volatilitaet = np.sqrt(w_opt @ SIGMA @ w_opt)\n",
    "    entropie = -np.sum(w_opt * np.log(w_opt))\n",
    "\n",
    "    print(f\"\\nKonvergenz:                {bestes.message}\")\n",
    "    print(f\"Iterationen:               {bestes.nit}\")\n",
    "    print(f\"Erwartete Jahresrendite:   {rendite * 100:6.2f} %\")\n",
    "    print(f\"Erwartete Volatilitaet:    {volatilitaet * 100:6.2f} %\")\n",
    "    print(f\"Diversifikations-Entropie: {entropie:.4f} \"\n",
    "          f\"(Maximum bei Gleichgewichtung: {np.log(n):.4f})\")\n",
    "\n",
    "    print(\"\\n\" + pd.DataFrame({\n",
    "        \"Asset\": ASSETS,\n",
    "        \"Erw. Rendite\": [f\"{r*100:.1f} %\" for r in MU],\n",
    "        \"Volatilitaet\": [f\"{v*100:.1f} %\" for v in VOLA],\n",
    "        \"Gewicht\": [f\"{w*100:6.2f} %\" for w in w_opt],\n",
    "    }).to_string(index=False))\n",
    "\n",
    "    # --- Vergleich: was passiert ohne Entropieterm? ----------------------\n",
    "    ohne = optimiere(np.ones(n) / n, beta=0.0)\n",
    "    print(\"\\n--- Wirkung des Entropieterms ---\")\n",
    "    print(f\"  Mit Entropie  (beta={BETA}): Gewichte {np.round(w_opt * 100, 1)}\")\n",
    "    print(f\"  Ohne Entropie (beta=0):     Gewichte {np.round(ohne.x * 100, 1)}\")\n",
    "    print(\"  Der Entropieterm zieht Kapital aus der Spitzenposition heraus,\")\n",
    "    print(\"  ohne dass eine harte Obergrenze noetig waere.\")\n",
    "    print(\"=\" * 74)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Was ein lokales Optimum praktisch bedeutet\n",
    "\n",
    "`Lokale_Optima_Multistart.py`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "\n",
    "# Lokale_Optima_Multistart.py\n",
    "\"\"\"\n",
    "Kapitel QP/NLP: Was passiert, wenn die Konvexitaet fehlt.\n",
    "\n",
    "CVXPY verweigert nicht-konvexe Probleme - das ist sein Schutzmechanismus.\n",
    "scipy.optimize.minimize verweigert nichts. Es rechnet, meldet 'success: True'\n",
    "und liefert ein Ergebnis. Nur ist das Ergebnis dann kein Optimum, sondern\n",
    "irgendein lokales Minimum, das vom Startpunkt abhaengt.\n",
    "\n",
    "Beispiel aus dem Einkauf: 700 Tonnen Rohstoff werden auf vier Lieferanten\n",
    "verteilt. Jeder gewaehrt einen MENGENRABATT - der Stueckpreis faellt, je mehr\n",
    "man bei ihm bestellt:\n",
    "\n",
    "    preis_i(q) = basis_i * (1 - rabatt_i * (1 - exp(-q / skala_i)))\n",
    "\n",
    "Genau das macht die Zielfunktion nicht-konvex: Grosse Bestellungen lohnen\n",
    "sich ueberproportional, es gibt also mehrere sinnvolle \"Cluster\"-Loesungen -\n",
    "und dazwischen schlechtere Taeler.\n",
    "\n",
    "Das Programm zeigt drei Dinge:\n",
    "  1. Ein einzelner Lauf liefert ein plausibles Ergebnis - ohne jede Warnung.\n",
    "  2. 200 Startpunkte foerdern mehrere verschiedene lokale Optima zutage.\n",
    "  3. Der Unterschied zwischen bestem und schlechtestem betraegt hier 8 %.\n",
    "\n",
    "Benoetigt: numpy, scipy\n",
    "\"\"\"\n",
    "\n",
    "from __future__ import annotations\n",
    "\n",
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "# Vier Lieferanten\n",
    "NAMEN = [\"Nord AG\", \"Ost GmbH\", \"Sued KG\", \"West SE\"]\n",
    "BASISPREIS = np.array([50.0, 47.0, 53.0, 45.0])     # EUR je Tonne ohne Rabatt\n",
    "MAX_RABATT = np.array([0.30, 0.18, 0.35, 0.12])     # hoechstmoeglicher Rabatt\n",
    "RABATT_SKALA = np.array([120.0, 260.0, 90.0, 300.0])  # wie schnell er greift\n",
    "KAPAZITAET = np.array([400.0, 400.0, 400.0, 400.0])\n",
    "BEDARF = 700.0\n",
    "\n",
    "RNG = np.random.default_rng(0)\n",
    "\n",
    "\n",
    "def gesamtkosten(menge: np.ndarray) -> float:\n",
    "    \"\"\"Einkaufskosten bei mengenabhaengigem Stueckpreis.\n",
    "\n",
    "    Der Rabatt waechst mit der Bestellmenge und laeuft gegen MAX_RABATT.\n",
    "    Dadurch ist der Stueckpreis fallend - und die Gesamtkostenfunktion\n",
    "    nicht mehr konvex.\n",
    "    \"\"\"\n",
    "    stueckpreis = BASISPREIS * (1 - MAX_RABATT * (1 - np.exp(-menge / RABATT_SKALA)))\n",
    "    return float(stueckpreis @ menge)\n",
    "\n",
    "\n",
    "NEBENBEDINGUNGEN = [{\"type\": \"eq\", \"fun\": lambda q: q.sum() - BEDARF}]\n",
    "GRENZEN = [(0.0, k) for k in KAPAZITAET]\n",
    "\n",
    "\n",
    "def optimiere_von(startpunkt: np.ndarray):\n",
    "    \"\"\"Ein Lauf von einem gegebenen Startpunkt aus.\"\"\"\n",
    "    return minimize(gesamtkosten, startpunkt, method=\"SLSQP\",\n",
    "                    bounds=GRENZEN, constraints=NEBENBEDINGUNGEN)\n",
    "\n",
    "\n",
    "def zufaelliger_start() -> np.ndarray:\n",
    "    \"\"\"Zufaellige Aufteilung, die den Bedarf bereits erfuellt.\"\"\"\n",
    "    anteil = RNG.random(len(NAMEN))\n",
    "    return anteil / anteil.sum() * BEDARF\n",
    "\n",
    "\n",
    "def zeige_plan(titel: str, menge: np.ndarray, kosten: float) -> None:\n",
    "    print(f\"\\n{titel}\")\n",
    "    for name, m in zip(NAMEN, menge):\n",
    "        anteil = m / BEDARF * 100\n",
    "        print(f\"    {name:<10} {m:7.1f} t  ({anteil:4.1f} %)\")\n",
    "    print(f\"    {'Gesamtkosten':<10} {kosten:9,.2f} EUR\")\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    print(\"=\" * 78)\n",
    "    print(\"  NICHT-KONVEX: DERSELBE CODE, VERSCHIEDENE ERGEBNISSE\")\n",
    "    print(\"=\" * 78)\n",
    "    print(f\"{BEDARF:.0f} t Rohstoff auf {len(NAMEN)} Lieferanten mit Mengenrabatt.\")\n",
    "\n",
    "    # --- 1. Ein einziger Lauf, so wie man es zuerst schreibt --------------\n",
    "    erster = optimiere_von(np.full(len(NAMEN), BEDARF / len(NAMEN)))\n",
    "    print(f\"\\n[1] EIN Lauf, Startpunkt 'gleichmaessig verteilt'\")\n",
    "    print(f\"    scipy meldet: success={erster.success}, \"\n",
    "          f\"'{erster.message}'\")\n",
    "    zeige_plan(\"    Ergebnis:\", erster.x, erster.fun)\n",
    "    print(\"\\n    Nichts an dieser Ausgabe deutet darauf hin, dass etwas fehlt.\")\n",
    "\n",
    "    # --- 1b. Der kaufmaennisch naheliegende Startpunkt --------------------\n",
    "    # \"Kaufe bei den beiden Lieferanten mit dem guenstigsten Basispreis\" -\n",
    "    # West SE (45) und Ost GmbH (47), jeweils bis zur Kapazitaetsgrenze.\n",
    "    guenstigste = np.argsort(BASISPREIS)[:2]\n",
    "    kaufmaennisch = np.zeros(len(NAMEN))\n",
    "    rest = BEDARF\n",
    "    for i in guenstigste:\n",
    "        kaufmaennisch[i] = min(KAPAZITAET[i], rest)\n",
    "        rest -= kaufmaennisch[i]\n",
    "    zweiter = optimiere_von(kaufmaennisch)\n",
    "    print(f\"\\n[1b] EIN Lauf, Startpunkt 'die zwei mit dem guenstigsten Basispreis'\")\n",
    "    print(f\"    scipy meldet: success={zweiter.success}\")\n",
    "    zeige_plan(\"    Ergebnis:\", zweiter.x, zweiter.fun)\n",
    "    print(f\"\\n    Dasselbe Programm, derselbe Aufruf, ein anderer Startpunkt -\")\n",
    "    print(f\"    und {zweiter.fun - erster.fun:,.2f} EUR Unterschied \"\n",
    "          f\"({(zweiter.fun / erster.fun - 1) * 100:.1f} %).\")\n",
    "\n",
    "    # --- 2. Multistart: dasselbe Problem, viele Startpunkte ---------------\n",
    "    laeufe = []\n",
    "    for _ in range(200):\n",
    "        ergebnis = optimiere_von(zufaelliger_start())\n",
    "        if ergebnis.success:\n",
    "            laeufe.append((float(ergebnis.fun), ergebnis.x))\n",
    "\n",
    "    # Ergebnisse, die sich um weniger als 1 Cent unterscheiden, sind dasselbe\n",
    "    # lokale Optimum - zusammenfassen, sonst zaehlt man Rundungsrauschen.\n",
    "    optima: list[tuple[float, np.ndarray]] = []\n",
    "    for wert, plan in sorted(laeufe, key=lambda t: t[0]):\n",
    "        if not optima or abs(wert - optima[-1][0]) > 0.01:\n",
    "            optima.append((wert, plan))\n",
    "\n",
    "    print(\"\\n\" + \"=\" * 78)\n",
    "    print(f\"[2] 200 zufaellige Startpunkte -> {len(laeufe)} erfolgreiche Laeufe\")\n",
    "    print(f\"    darunter {len(optima)} VERSCHIEDENE lokale Optima:\")\n",
    "    print()\n",
    "    print(f\"    {'Rang':>5} {'Kosten':>13} {'Abstand zum besten':>20}   Aufteilung (t)\")\n",
    "    print(\"    \" + \"-\" * 70)\n",
    "    bester = optima[0][0]\n",
    "    for rang, (wert, plan) in enumerate(optima, start=1):\n",
    "        abstand = (wert / bester - 1) * 100\n",
    "        aufteilung = \" \".join(f\"{m:5.0f}\" for m in plan)\n",
    "        print(f\"    {rang:>5} {wert:>13,.2f} {abstand:>19.2f} %   {aufteilung}\")\n",
    "\n",
    "    # --- 3. Was das kostet ------------------------------------------------\n",
    "    schlechtester = optima[-1]\n",
    "    print(\"\\n\" + \"=\" * 78)\n",
    "    print(\"  WAS AUF DEM SPIEL STEHT\")\n",
    "    print(\"=\" * 78)\n",
    "    zeige_plan(\"Bester gefundener Plan:\", optima[0][1], optima[0][0])\n",
    "    zeige_plan(\"Schlechtestes lokales Optimum:\", schlechtester[1], schlechtester[0])\n",
    "    unterschied = schlechtester[0] - optima[0][0]\n",
    "    print(f\"\\n  Unterschied: {unterschied:,.2f} EUR \"\n",
    "          f\"({unterschied / optima[0][0] * 100:.1f} %)\")\n",
    "    print(f\"\\n  Lauf [1]  (gleichmaessiger Start):     {erster.fun:>10,.2f} EUR\")\n",
    "    print(f\"  Lauf [1b] (kaufmaennischer Start):    {zweiter.fun:>10,.2f} EUR\"\n",
    "          f\"   <- {(zweiter.fun / optima[0][0] - 1) * 100:.1f} % ueber dem besten\")\n",
    "    print()\n",
    "    print(\"  Bemerkenswert: Der kaufmaennisch NAHELIEGENDE Startpunkt fuehrt in\")\n",
    "    print(\"  das schlechteste Ergebnis von allen - schlechter als jedes der 200\")\n",
    "    print(\"  zufaellig gefundenen lokalen Optima. Wer beim guenstigsten\")\n",
    "    print(\"  Basispreis anfaengt, uebersieht, dass hier der Mengenrabatt\")\n",
    "    print(\"  entscheidet und nicht der Listenpreis.\")\n",
    "\n",
    "    print()\n",
    "    print(\"Drei Konsequenzen fuer die Praxis:\")\n",
    "    print(\"  1. 'success: True' heisst bei nicht-konvexen Problemen NICHT 'optimal'.\")\n",
    "    print(\"     Es heisst nur: 'Ich bin an einer Stelle angekommen, an der es in\")\n",
    "    print(\"     keine Richtung mehr bergab geht.'\")\n",
    "    print(\"  2. Ein einzelner Lauf ist wertlos. Nehmen Sie viele Startpunkte und\")\n",
    "    print(\"     berichten Sie die STREUUNG mit - sie ist Ihre einzige Auskunft\")\n",
    "    print(\"     darueber, wie zerklueftet die Landschaft ist.\")\n",
    "    print(\"  3. Auch Multistart liefert KEINE Garantie. Dass hier nichts unter\")\n",
    "    print(f\"     {bester:,.2f} EUR gefunden wurde, beweist nicht, dass es nichts gibt.\")\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
}