operations_research/OR_HTML_04/Notebooks_04/oekosystem.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kapitel 3: Das Python-Ökosystem für OR — Solver, Bindings und Modellierungsschichten\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": [
    "## Die Entscheidung in drei Fragen\n",
    "\n",
    "`Solver_Wahl.py`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "\n",
    "# Solver_Wahl.py\n",
    "\"\"\"\n",
    "Kapitel Oekosystem: Entscheidungshilfe zur Solverwahl.\n",
    "Beantwortet drei Fragen und empfiehlt Bibliothek + Backend.\n",
    "\"\"\"\n",
    "\n",
    "from dataclasses import dataclass\n",
    "\n",
    "\n",
    "@dataclass\n",
    "class Problem:\n",
    "    diskrete_entscheidungen: bool   # Ja/Nein-Variablen, Zuordnungen, Reihenfolgen?\n",
    "    zielfunktion: str               # \"linear\" | \"quadratisch\" | \"konvex\" | \"beliebig\"\n",
    "    routing: bool = False           # Fahrzeuge, Depots, Zeitfenster?\n",
    "    groesse: str = \"klein\"          # \"klein\" (<10^3 Var.) | \"gross\"\n",
    "    wiederholte_laeufe: bool = False\n",
    "\n",
    "\n",
    "def empfehle(p: Problem) -> tuple[str, str, str]:\n",
    "    \"\"\"Gibt (Bibliothek, Backend, Begründung) zurück.\"\"\"\n",
    "    if p.routing:\n",
    "        return (\"ortools.constraint_solver (Routing)\", \"Routing Engine\",\n",
    "                \"Spezialisierte Heuristiken für VRP/TSP schlagen jedes selbstgebaute Modell.\")\n",
    "\n",
    "    if p.diskrete_entscheidungen:\n",
    "        if p.zielfunktion in (\"linear\",) and not p.routing:\n",
    "            if p.groesse == \"gross\":\n",
    "                return (\"highspy\", \"HiGHS Branch-and-Cut\",\n",
    "                        \"MILP mit ökonomischer Struktur (Fixkosten, Kardinalität).\")\n",
    "            return (\"ortools.sat (CP-SAT)\", \"CP-SAT\",\n",
    "                    \"Logische Regeln und Zuweisungen: CP-SAT propagiert sehr effizient.\")\n",
    "        return (\"ortools.sat (CP-SAT)\", \"CP-SAT\",\n",
    "                \"Diskrete Struktur dominiert; CP-SAT verarbeitet auch nichtlineare Logik.\")\n",
    "\n",
    "    if p.zielfunktion == \"linear\":\n",
    "        if p.wiederholte_laeufe or p.groesse == \"gross\":\n",
    "            return (\"highspy\", \"HiGHS Dual Simplex\",\n",
    "                    \"Modell einmal aufbauen, Parameter ändern, wiederholt lösen.\")\n",
    "        return (\"scipy.optimize.linprog\", \"HiGHS\",\n",
    "                \"Kleinstes Setup, keine zusätzliche Abhängigkeit.\")\n",
    "\n",
    "    if p.zielfunktion in (\"quadratisch\", \"konvex\"):\n",
    "        return (\"cvxpy\", \"Clarabel / OSQP\",\n",
    "                \"Konvexität wird automatisch geprüft; Portfolio-Standard.\")\n",
    "\n",
    "    return (\"scipy.optimize.minimize\", \"SLSQP / trust-constr\",\n",
    "            \"Nicht konvex: nur lokales Optimum, Startpunkt variieren und vergleichen!\")\n",
    "\n",
    "\n",
    "BEISPIELE = {\n",
    "    \"Vertretungsplan Schule\":        Problem(True,  \"linear\", groesse=\"klein\"),\n",
    "    \"Produktionsplanung (LP)\":       Problem(False, \"linear\", groesse=\"klein\"),\n",
    "    \"Produktionsplanung, 50k Var.\":  Problem(False, \"linear\", groesse=\"gross\"),\n",
    "    \"Portfolio Markowitz\":           Problem(False, \"quadratisch\"),\n",
    "    \"Portfolio mit max. 5 Titeln\":   Problem(True,  \"quadratisch\"),\n",
    "    \"Liefertouren mit Zeitfenstern\": Problem(True,  \"linear\", routing=True),\n",
    "    \"Entropie-Allokation (NLP)\":     Problem(False, \"beliebig\"),\n",
    "    \"Backtest, 60x neu optimieren\":  Problem(False, \"konvex\", wiederholte_laeufe=True),\n",
    "}\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    print(\"=\" * 92)\n",
    "    print(\"  SOLVER-EMPFEHLUNG\")\n",
    "    print(\"=\" * 92)\n",
    "    for name, p in BEISPIELE.items():\n",
    "        lib, backend, grund = empfehle(p)\n",
    "        print(f\"\\n{name}\")\n",
    "        print(f\"   -> Bibliothek: {lib}\")\n",
    "        print(f\"      Backend:    {backend}\")\n",
    "        print(f\"      Grund:      {grund}\")\n",
    "    print(\"\\n\" + \"=\" * 92)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wie die Isolation aussieht, wenn sie tragen soll\n",
    "\n",
    "`Ein_System_Vier_Ansaetze.py`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "\n",
    "# Ein_System_Vier_Ansaetze.py\n",
    "\"\"\"\n",
    "Kapitel Oekosystem: Dasselbe LP in vier Bibliotheken.\n",
    "   max 10*x1 + 15*x2 + 25*x3\n",
    "   u.d.N. x1 + x2 + 2*x3 <= 40\n",
    "          2*x1 + 3*x2 + x3 <= 50\n",
    "          x >= 0\n",
    "\n",
    "Deckt scipy.optimize, highspy, CVXPY und OR-Tools/GLOP ab, mit Kreuzvergleich\n",
    "am Ende.\n",
    "\n",
    "WICHTIG: Jeder Solver laeuft in einem EIGENEN Prozess, weil sich ortools und\n",
    "highspy auf vielen Systemen nicht gemeinsam importieren lassen (beide bringen\n",
    "eine eigene HiGHS-Kopie mit -> Symbolkonflikt).\n",
    "\n",
    "Die Isolation besorgt ein ProcessPoolExecutor. Drei Einstellungen ergeben\n",
    "zusammen die Garantie:\n",
    "\n",
    "  mp_context \"spawn\"     Der Kindprozess startet mit einem FRISCHEN\n",
    "                         Interpreter, statt den Speicher des Elternprozesses\n",
    "                         zu erben. Was hier schon importiert ist, ist dort\n",
    "                         nicht importiert. Mit dem Standard \"fork\" auf Linux\n",
    "                         waere das nicht so.\n",
    "  max_tasks_per_child=1  Jede Aufgabe bekommt einen NEUEN Prozess. Ohne das\n",
    "                         wuerde der Pool seinen Arbeiter wiederverwenden - und\n",
    "                         beim zweiten Solver waere der Konflikt zurueck.\n",
    "  max_workers=1          Haelt die vier Laeufe nacheinander. Nicht aus\n",
    "                         Vorsicht, sondern damit die gemessenen Zeiten\n",
    "                         vergleichbar bleiben.\n",
    "\n",
    "Jeder Solver steht in einer eigenen Funktion mit LOKALEM Import. Das ist der\n",
    "Unterschied zu einem Codestring, den man an 'python -c' uebergibt: Die\n",
    "Funktion laesst sich einzeln aufrufen, testen und vom Editor pruefen - ein\n",
    "String nicht.\n",
    "\n",
    "Benoetigt: scipy, highspy, cvxpy, ortools\n",
    "\"\"\"\n",
    "\n",
    "import multiprocessing\n",
    "import time\n",
    "from concurrent.futures import ProcessPoolExecutor\n",
    "\n",
    "ERWARTET = 530.0            # Ergebnis der Handrechnung zum Produktionsprogramm\n",
    "\n",
    "# Die Instanz - einmal notiert, von allen vier Funktionen benutzt.\n",
    "ZIEL = [10.0, 15.0, 25.0]\n",
    "MATRIX = [[1, 1, 2], [2, 3, 1]]\n",
    "KAPAZITAET = [40.0, 50.0]\n",
    "\n",
    "\n",
    "def loese_mit_scipy() -> tuple[float, list[float]]:\n",
    "    from scipy.optimize import linprog\n",
    "    ergebnis = linprog(c=[-w for w in ZIEL],          # linprog MINIMIERT -> negieren\n",
    "                       A_ub=MATRIX, b_ub=KAPAZITAET,\n",
    "                       bounds=[(0, None)] * 3, method=\"highs\")\n",
    "    return -ergebnis.fun, list(ergebnis.x)\n",
    "\n",
    "\n",
    "def loese_mit_highspy() -> tuple[float, list[float]]:\n",
    "    import highspy\n",
    "    import numpy as np\n",
    "    h = highspy.Highs()\n",
    "    h.setOptionValue(\"output_flag\", False)\n",
    "    h.addVars(3, np.zeros(3), np.full(3, highspy.kHighsInf))\n",
    "    h.changeObjectiveSense(highspy.ObjSense.kMaximize)\n",
    "    for j, wert in enumerate(ZIEL):\n",
    "        h.changeColCost(j, wert)\n",
    "    # CSR-Format: starts[i] = Beginn von Zeile i in indices/values\n",
    "    h.addRows(2, np.full(2, -highspy.kHighsInf), np.array(KAPAZITAET), 6,\n",
    "              np.array([0, 3], dtype=np.int32),\n",
    "              np.array([0, 1, 2, 0, 1, 2], dtype=np.int32),\n",
    "              np.array([float(w) for zeile in MATRIX for w in zeile]))\n",
    "    h.run()\n",
    "    return (h.getInfo().objective_function_value,\n",
    "            list(h.getSolution().col_value[:3]))\n",
    "\n",
    "\n",
    "def loese_mit_cvxpy() -> tuple[float, list[float]]:\n",
    "    import cvxpy as cp\n",
    "    import numpy as np\n",
    "    x = cp.Variable(3, nonneg=True)\n",
    "    problem = cp.Problem(cp.Maximize(np.array(ZIEL) @ x),\n",
    "                         [np.array(MATRIX) @ x <= np.array(KAPAZITAET)])\n",
    "    problem.solve()\n",
    "    return float(problem.value), [float(v) for v in x.value]\n",
    "\n",
    "\n",
    "def loese_mit_ortools() -> tuple[float, list[float]]:\n",
    "    from ortools.linear_solver import pywraplp\n",
    "    s = pywraplp.Solver.CreateSolver(\"GLOP\")\n",
    "    x = [s.NumVar(0, s.infinity(), f\"x{j+1}\") for j in range(3)]\n",
    "    for i, kapazitaet in enumerate(KAPAZITAET):\n",
    "        s.Add(sum(MATRIX[i][j] * x[j] for j in range(3)) <= kapazitaet)\n",
    "    s.Maximize(sum(ZIEL[j] * x[j] for j in range(3)))\n",
    "    s.Solve()\n",
    "    return s.Objective().Value(), [v.solution_value() for v in x]\n",
    "\n",
    "\n",
    "ANSAETZE = {\n",
    "    \"scipy.optimize.linprog\": loese_mit_scipy,\n",
    "    \"highspy (natives HiGHS)\": loese_mit_highspy,\n",
    "    \"cvxpy\": loese_mit_cvxpy,\n",
    "    \"ortools / GLOP\": loese_mit_ortools,\n",
    "}\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    print(\"=\" * 78)\n",
    "    print(\"  EIN SYSTEM - VIER ANSAETZE (je eigener Prozess)\")\n",
    "    print(\"=\" * 78)\n",
    "    print(f\"{'Bibliothek':<26} {'Z*':>10} {'x1':>7} {'x2':>7} {'x3':>7} {'Zeit':>10}\")\n",
    "    print(\"-\" * 78)\n",
    "\n",
    "    werte = []\n",
    "    # Ein Pool, vier Aufgaben, vier frische Prozesse. Der Kontext muss\n",
    "    # \"spawn\" sein - siehe Modulkommentar.\n",
    "    with ProcessPoolExecutor(\n",
    "            max_workers=1,\n",
    "            mp_context=multiprocessing.get_context(\"spawn\"),\n",
    "            max_tasks_per_child=1) as pool:\n",
    "        for name, funktion in ANSAETZE.items():\n",
    "            beginn = time.perf_counter()\n",
    "            try:\n",
    "                wert, x = pool.submit(funktion).result(timeout=120)\n",
    "            except Exception as fehler:                     # Bibliothek fehlt o. Ae.\n",
    "                print(f\"{name:<26} nicht verfuegbar: {str(fehler)[:40]}\")\n",
    "                continue\n",
    "            dauer = time.perf_counter() - beginn\n",
    "            werte.append(wert)\n",
    "            print(f\"{name:<26} {wert:>10.2f} {x[0]:>7.2f} {x[1]:>7.2f} {x[2]:>7.2f} \"\n",
    "                  f\"{dauer:>8.2f} s\")\n",
    "\n",
    "    print(\"-\" * 78)\n",
    "    spanne = max(werte) - min(werte)\n",
    "    print(f\"Spannweite zwischen den Bibliotheken: {spanne:.2e}\")\n",
    "    print(f\"Abweichung zur Handrechnung ({ERWARTET:.0f}):        \"\n",
    "          f\"{abs(werte[0] - ERWARTET):.2e}\")\n",
    "    assert spanne < 1e-6, \"Die Bibliotheken widersprechen sich!\"\n",
    "    assert abs(werte[0] - ERWARTET) < 1e-6, \"Ergebnis weicht von der Handrechnung ab!\"\n",
    "    print(\"Alle Wege fuehren zum selben, von Hand bestaetigten Optimum.\")\n",
    "    print(\"(Die Zeiten enthalten Prozessstart und Import - sie messen NICHT die\")\n",
    "    print(\" reine Solverleistung. Die Uebungsaufgabe 'Laufzeitvergleich' trennt beides.)\")\n",
    "    print(\"=\" * 78)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Die Schichten im Überblick\n",
    "\n",
    "`Modellierungsschichten.py`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "\n",
    "# Modellierungsschichten.py\n",
    "\"\"\"\n",
    "Kapitel Oekosystem: Pyomo und Linopy - zwei Modellierungsschichten fuer grosse Modelle.\n",
    "\n",
    "Geloest wird dasselbe Produktionsproblem wie im Vierfach-Vergleich:\n",
    "\n",
    "    max 10*x1 + 15*x2 + 25*x3\n",
    "    u.d.N.  x1 +  x2 + 2*x3 <= 40\n",
    "           2*x1 + 3*x2 +  x3 <= 50\n",
    "           x >= 0\n",
    "\n",
    "Handrechnung: Z* = 530 bei x = (0, 12, 14).\n",
    "\n",
    "Der Vergleich zeigt die beiden Denkweisen:\n",
    "  * Pyomo  - algebraisch, indexbasiert, Industriestandard fuer Grossmodelle,\n",
    "             trennt Modellstruktur sauber von den Daten (AbstractModel).\n",
    "  * Linopy - beschriftete Arrays (xarray): eine Zeile Code erzeugt Tausende\n",
    "             Nebenbedingungen auf einmal, ohne Python-Schleife.\n",
    "\n",
    "Beide bringen KEINEN eigenen Solver mit; hier rechnet in beiden Faellen HiGHS.\n",
    "\n",
    "WICHTIG: ortools wird in diesem Programm bewusst NICHT importiert - es\n",
    "vertraegt sich nicht mit der HiGHS-Kopie, die Pyomo und Linopy laden\n",
    "(siehe die Stolperfalle im Abschnitt 'Ein System - vier Programmieransaetze').\n",
    "\n",
    "Benoetigt: pyomo, linopy, xarray, pandas, highspy, numpy\n",
    "\"\"\"\n",
    "\n",
    "from __future__ import annotations\n",
    "\n",
    "import time\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pyomo.environ as pyo\n",
    "import xarray as xr\n",
    "import linopy\n",
    "\n",
    "ERWARTET = 530.0                     # Ergebnis der Handrechnung\n",
    "\n",
    "PRODUKTE = [\"Standard\", \"Komfort\", \"Premium\"]\n",
    "RESSOURCEN = [\"Material\", \"Montage\"]\n",
    "\n",
    "DECKUNGSBEITRAG = np.array([10.0, 15.0, 25.0])\n",
    "VERBRAUCH = np.array([[1.0, 1.0, 2.0],       # Material je Produkt\n",
    "                      [2.0, 3.0, 1.0]])      # Montage je Produkt\n",
    "VORRAT = np.array([40.0, 50.0])\n",
    "\n",
    "\n",
    "# --- Pyomo: algebraisch und indexbasiert ------------------------------------\n",
    "\n",
    "def loese_mit_pyomo() -> tuple[float, list[float], float]:\n",
    "    \"\"\"Pyomo denkt in Mengen und Indizes, wie ein Mathematiker es aufschreibt.\n",
    "\n",
    "    `Constraint(RESSOURCEN, rule=...)` erzeugt fuer JEDES Element der Menge\n",
    "    eine Nebenbedingung - das ist das 'fuer alle i' der Formelsprache,\n",
    "    unmittelbar in Code uebersetzt.\n",
    "    \"\"\"\n",
    "    t0 = time.perf_counter()\n",
    "\n",
    "    modell = pyo.ConcreteModel(name=\"Produktionsprogramm\")\n",
    "    modell.P = pyo.Set(initialize=PRODUKTE)\n",
    "    modell.R = pyo.Set(initialize=RESSOURCEN)\n",
    "\n",
    "    modell.db = pyo.Param(modell.P, initialize=dict(zip(PRODUKTE, DECKUNGSBEITRAG)))\n",
    "    modell.a = pyo.Param(modell.R, modell.P, initialize={\n",
    "        (r, p): VERBRAUCH[i, j]\n",
    "        for i, r in enumerate(RESSOURCEN) for j, p in enumerate(PRODUKTE)})\n",
    "    modell.vorrat = pyo.Param(modell.R, initialize=dict(zip(RESSOURCEN, VORRAT)))\n",
    "\n",
    "    modell.x = pyo.Var(modell.P, domain=pyo.NonNegativeReals)\n",
    "\n",
    "    modell.ziel = pyo.Objective(\n",
    "        expr=sum(modell.db[p] * modell.x[p] for p in modell.P),\n",
    "        sense=pyo.maximize)\n",
    "\n",
    "    def kapazitaet(m, r):\n",
    "        return sum(m.a[r, p] * m.x[p] for p in m.P) <= m.vorrat[r]\n",
    "\n",
    "    modell.kapazitaet = pyo.Constraint(modell.R, rule=kapazitaet)\n",
    "\n",
    "    ergebnis = pyo.SolverFactory(\"appsi_highs\").solve(modell)\n",
    "    dauer = time.perf_counter() - t0\n",
    "\n",
    "    status = ergebnis.solver.termination_condition\n",
    "    if status != pyo.TerminationCondition.optimal:\n",
    "        raise RuntimeError(f\"Pyomo meldet Status: {status}\")\n",
    "\n",
    "    return (float(pyo.value(modell.ziel)),\n",
    "            [float(pyo.value(modell.x[p])) for p in PRODUKTE],\n",
    "            dauer)\n",
    "\n",
    "\n",
    "# --- Linopy: beschriftete Arrays --------------------------------------------\n",
    "\n",
    "def loese_mit_linopy() -> tuple[float, list[float], float]:\n",
    "    \"\"\"Linopy denkt in beschrifteten Arrays (xarray).\n",
    "\n",
    "    Der entscheidende Unterschied: `(verbrauch * x).sum(\"produkt\") <= vorrat`\n",
    "    ist EINE Zeile und erzeugt so viele Nebenbedingungen, wie die Dimension\n",
    "    'ressource' Eintraege hat. Bei 2 Ressourcen faellt das nicht auf, bei\n",
    "    200 000 schon - dort entstehen sie als Matrixoperation statt in einer\n",
    "    Python-Schleife.\n",
    "    \"\"\"\n",
    "    t0 = time.perf_counter()\n",
    "\n",
    "    # Benannte Indizes statt blosser Listen: Dadurch heissen die Achsen\n",
    "    # 'produkt' und 'ressource', und xarray fuehrt sie beim Rechnen von allein\n",
    "    # richtig zusammen. Ohne Namen vergibt linopy 'dim_0', und man muss\n",
    "    # spaeter umbenennen - eine haeufige Stolperstelle.\n",
    "    produkt = pd.Index(PRODUKTE, name=\"produkt\")\n",
    "    ressource = pd.Index(RESSOURCEN, name=\"ressource\")\n",
    "\n",
    "    modell = linopy.Model()\n",
    "    modell.add_variables(lower=0, coords=[produkt], name=\"menge\")\n",
    "    x = modell.variables[\"menge\"]\n",
    "\n",
    "    db = xr.DataArray(DECKUNGSBEITRAG, coords=[produkt])\n",
    "    verbrauch = xr.DataArray(VERBRAUCH, coords=[ressource, produkt])\n",
    "    vorrat = xr.DataArray(VORRAT, coords=[ressource])\n",
    "\n",
    "    # EINE Zeile - sie erzeugt so viele Nebenbedingungen, wie die Achse\n",
    "    # 'ressource' Eintraege hat. Genau das ist der Punkt.\n",
    "    modell.add_constraints((verbrauch * x).sum(\"produkt\") <= vorrat,\n",
    "                           name=\"kapazitaet\")\n",
    "    modell.add_objective((db * x).sum(), sense=\"max\")\n",
    "\n",
    "    modell.solve(solver_name=\"highs\", output_flag=False)\n",
    "    dauer = time.perf_counter() - t0\n",
    "\n",
    "    if modell.termination_condition != \"optimal\":\n",
    "        raise RuntimeError(f\"Linopy meldet Status: {modell.termination_condition}\")\n",
    "\n",
    "    loesung = modell.variables[\"menge\"].solution.to_series()\n",
    "    return (float(modell.objective.value),\n",
    "            [float(loesung[p]) for p in PRODUKTE],\n",
    "            dauer)\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    print(\"=\" * 80)\n",
    "    print(\"  MODELLIERUNGSSCHICHTEN FUER GROSSE MODELLE\")\n",
    "    print(\"=\" * 80)\n",
    "    print(f\"{'Schicht':<14} {'Z*':>10} {'Standard':>10} {'Komfort':>10} \"\n",
    "          f\"{'Premium':>10} {'Zeit':>10}\")\n",
    "    print(\"-\" * 80)\n",
    "\n",
    "    ergebnisse = []\n",
    "    for name, loeser in [(\"Pyomo\", loese_mit_pyomo), (\"Linopy\", loese_mit_linopy)]:\n",
    "        ziel, mengen, dauer = loeser()\n",
    "        ergebnisse.append(ziel)\n",
    "        print(f\"{name:<14} {ziel:>10.2f} {mengen[0]:>10.2f} {mengen[1]:>10.2f} \"\n",
    "              f\"{mengen[2]:>10.2f} {dauer:>8.2f} s\")\n",
    "\n",
    "    print(\"-\" * 80)\n",
    "    for ziel in ergebnisse:\n",
    "        assert abs(ziel - ERWARTET) < 1e-6, \\\n",
    "            f\"Abweichung von der Handrechnung: {ziel} statt {ERWARTET}\"\n",
    "    print(f\"Beide stimmen mit der Handrechnung ueberein (Z* = {ERWARTET:.0f}).\")\n",
    "    print()\n",
    "    print(\"Wann welche Schicht?\")\n",
    "    print(\"  Pyomo  -> wenn Modellstruktur und Daten getrennt bleiben sollen,\")\n",
    "    print(\"            wenn nichtlineare Terme oder MINLP dazukommen koennen,\")\n",
    "    print(\"            wenn spaeter ein kommerzieller Solver angebunden wird.\")\n",
    "    print(\"  Linopy -> wenn die Daten ohnehin als beschriftete Arrays vorliegen\")\n",
    "    print(\"            (Energiesystem-, Netz- und Zeitreihenmodelle) und das\")\n",
    "    print(\"            Modell zehntausende gleichartige Nebenbedingungen hat.\")\n",
    "    print(\"=\" * 80)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vier Stufen an einem Transportproblem\n",
    "\n",
    "`Vektorisierte_Modellgenerierung.py`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "\n",
    "# Vektorisierte_Modellgenerierung.py\n",
    "\"\"\"\n",
    "Kapitel Oekosystem: Warum der Solver oft gar nicht der Engpass ist.\n",
    "\n",
    "In realen Projekten geht ein grosser Teil der Rechenzeit nicht ins Loesen,\n",
    "sondern ins AUFBAUEN des Modells. Dieses Programm misst das an einem\n",
    "Transportproblem wachsender Groesse in vier Stufen:\n",
    "\n",
    "    A  Modellierungsschicht, Nebenbedingung fuer Nebenbedingung (OR-Tools)\n",
    "    B  Matrix direkt, aber mit Python-Schleifen ueber die Eintraege (COO)\n",
    "    C  Matrix vektorisiert ueber Kronecker-Produkte (NumPy/SciPy)\n",
    "    D  Daten kommen als lange Tabelle, aufbereitet mit Polars\n",
    "\n",
    "Alle Varianten loesen dasselbe Problem und muessen denselben Zielwert\n",
    "liefern - das wird am Ende geprueft.\n",
    "\n",
    "Benoetigt: numpy, scipy, ortools; Variante D zusaetzlich polars (optional).\n",
    "\"\"\"\n",
    "\n",
    "from __future__ import annotations\n",
    "\n",
    "import importlib.util\n",
    "import time\n",
    "\n",
    "import numpy as np\n",
    "import scipy.sparse as sp\n",
    "from ortools.linear_solver import pywraplp\n",
    "from scipy.optimize import linprog\n",
    "\n",
    "HAT_POLARS = importlib.util.find_spec(\"polars\") is not None\n",
    "\n",
    "\n",
    "def erzeuge_daten(m: int, n: int, saat: int = 3\n",
    "                  ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:\n",
    "    \"\"\"Transportproblem: m Werke, n Kunden.\n",
    "\n",
    "    Liefert (kosten[m, n], angebot[m], bedarf[n]). Das Gesamtangebot liegt\n",
    "    20 % ueber dem Gesamtbedarf, damit das Modell sicher loesbar ist.\n",
    "    \"\"\"\n",
    "    rng = np.random.default_rng(saat)\n",
    "    kosten = rng.uniform(1.0, 20.0, size=(m, n))\n",
    "    bedarf = rng.uniform(10.0, 50.0, size=n)\n",
    "    angebot = np.full(m, 1.2 * bedarf.sum() / m)\n",
    "    return kosten, angebot, bedarf\n",
    "\n",
    "\n",
    "# --- Variante A: Modellierungsschicht, Bedingung fuer Bedingung -------------\n",
    "\n",
    "def loese_mit_modellierungsschicht(kosten: np.ndarray, angebot: np.ndarray,\n",
    "                                   bedarf: np.ndarray) -> tuple[float, float, float]:\n",
    "    \"\"\"So schreibt man ein Transportproblem zuerst hin - gut lesbar, nah an\n",
    "    der mathematischen Formulierung, jede Nebenbedingung ein eigener Aufruf.\n",
    "\n",
    "    Jedes `s.Add(sum(...))` baut in Python einen Ausdrucksbaum aus m bzw. n\n",
    "    Termen auf und uebergibt ihn einzeln an die C++-Schicht. Das ist der\n",
    "    Preis der Bequemlichkeit - und er waechst linear mit der Modellgroesse.\n",
    "    \"\"\"\n",
    "    m, n = kosten.shape\n",
    "\n",
    "    t0 = time.perf_counter()\n",
    "    s = pywraplp.Solver.CreateSolver(\"GLOP\")\n",
    "    x = [[s.NumVar(0, s.infinity(), f\"x_{i}_{j}\") for j in range(n)]\n",
    "         for i in range(m)]\n",
    "    for i in range(m):\n",
    "        s.Add(sum(x[i][j] for j in range(n)) <= angebot[i])\n",
    "    for j in range(n):\n",
    "        s.Add(sum(x[i][j] for i in range(m)) >= bedarf[j])\n",
    "    s.Minimize(sum(kosten[i][j] * x[i][j] for i in range(m) for j in range(n)))\n",
    "    t_aufbau = time.perf_counter() - t0\n",
    "\n",
    "    t0 = time.perf_counter()\n",
    "    status = s.Solve()\n",
    "    t_loesen = time.perf_counter() - t0\n",
    "    if status != pywraplp.Solver.OPTIMAL:\n",
    "        raise RuntimeError(f\"Solver-Status: {status}\")\n",
    "    return t_aufbau, t_loesen, s.Objective().Value()\n",
    "\n",
    "\n",
    "# --- Varianten B bis D: Matrix selbst bauen, dann SciPy/HiGHS ---------------\n",
    "\n",
    "def baue_mit_schleifen(kosten: np.ndarray, angebot: np.ndarray,\n",
    "                       bedarf: np.ndarray):\n",
    "    \"\"\"Die Nebenbedingungsmatrix als COO-Tripel (Zeile, Spalte, Wert), erzeugt\n",
    "    in verschachtelten Python-Schleifen.\n",
    "\n",
    "    Schon deutlich naeher am Blech als Variante A - es entsteht kein\n",
    "    Ausdrucksbaum mehr. Die Schleife selbst bleibt aber Python.\n",
    "    \"\"\"\n",
    "    m, n = kosten.shape\n",
    "    zeilen: list[int] = []\n",
    "    spalten: list[int] = []\n",
    "    werte: list[float] = []\n",
    "    rechte_seite: list[float] = []\n",
    "\n",
    "    for i in range(m):                       # Angebot je Werk\n",
    "        for j in range(n):\n",
    "            zeilen.append(i)\n",
    "            spalten.append(i * n + j)\n",
    "            werte.append(1.0)\n",
    "        rechte_seite.append(float(angebot[i]))\n",
    "\n",
    "    for j in range(n):                       # Bedarf je Kunde, als -x <= -bedarf\n",
    "        for i in range(m):\n",
    "            zeilen.append(m + j)\n",
    "            spalten.append(i * n + j)\n",
    "            werte.append(-1.0)\n",
    "        rechte_seite.append(-float(bedarf[j]))\n",
    "\n",
    "    A_ub = sp.csr_matrix((werte, (zeilen, spalten)), shape=(m + n, m * n))\n",
    "    return A_ub, np.array(rechte_seite), kosten.ravel()\n",
    "\n",
    "\n",
    "def baue_vektorisiert(kosten: np.ndarray, angebot: np.ndarray,\n",
    "                      bedarf: np.ndarray):\n",
    "    \"\"\"Dieselbe Matrix ohne eine einzige Schleife - ueber Kronecker-Produkte.\n",
    "\n",
    "    Die Angebotsmatrix ist  kron(I_m, 1_n^T):  je Werk eine Zeile mit Einsen\n",
    "    an genau den n Spalten dieses Werks. Die Bedarfsmatrix ist\n",
    "    kron(1_m^T, I_n). Beide entstehen in je einem Aufruf und sind sofort\n",
    "    duennbesetzt.\n",
    "    \"\"\"\n",
    "    m, n = kosten.shape\n",
    "    angebots_matrix = sp.kron(sp.identity(m, format=\"csr\"), np.ones((1, n)))\n",
    "    bedarfs_matrix = sp.kron(np.ones((1, m)), sp.identity(n, format=\"csr\"))\n",
    "\n",
    "    A_ub = sp.vstack([angebots_matrix, -bedarfs_matrix], format=\"csr\")\n",
    "    b_ub = np.concatenate([angebot, -bedarf])\n",
    "    return A_ub, b_ub, kosten.ravel()\n",
    "\n",
    "\n",
    "def baue_mit_polars(kosten: np.ndarray, angebot: np.ndarray, bedarf: np.ndarray):\n",
    "    \"\"\"Der realistische Fall: Die Kosten kommen als LANGE Tabelle\n",
    "    (werk, kunde, kosten) aus Datenbank, Data Lake oder CSV-Datei.\n",
    "\n",
    "    Polars berechnet den Spaltenindex jeder Variablen in einem einzigen\n",
    "    Spaltenausdruck - ohne Python-Schleife ueber die Zeilen. Genau so baut man\n",
    "    Modelle aus Millionen Tabellenzeilen.\n",
    "    \"\"\"\n",
    "    import polars as pl\n",
    "\n",
    "    m, n = kosten.shape\n",
    "    tabelle = pl.DataFrame({\n",
    "        \"werk\": np.repeat(np.arange(m), n),\n",
    "        \"kunde\": np.tile(np.arange(n), m),\n",
    "        \"kosten\": kosten.ravel(),\n",
    "    }).with_columns(\n",
    "        (pl.col(\"werk\") * n + pl.col(\"kunde\")).alias(\"var_index\")\n",
    "    )\n",
    "\n",
    "    var_index = tabelle[\"var_index\"].to_numpy()\n",
    "    werk = tabelle[\"werk\"].to_numpy()\n",
    "    kunde = tabelle[\"kunde\"].to_numpy()\n",
    "    eins = np.ones(var_index.size)\n",
    "\n",
    "    angebots_matrix = sp.csr_matrix((eins, (werk, var_index)), shape=(m, m * n))\n",
    "    bedarfs_matrix = sp.csr_matrix((eins, (kunde, var_index)), shape=(n, m * n))\n",
    "\n",
    "    A_ub = sp.vstack([angebots_matrix, -bedarfs_matrix], format=\"csr\")\n",
    "    b_ub = np.concatenate([angebot, -bedarf])\n",
    "    return A_ub, b_ub, tabelle[\"kosten\"].to_numpy()\n",
    "\n",
    "\n",
    "def messe_matrixvariante(bauer, kosten, angebot, bedarf\n",
    "                         ) -> tuple[float, float, float]:\n",
    "    \"\"\"Liefert (Aufbauzeit, Loesezeit, Zielwert) fuer die Varianten B bis D.\"\"\"\n",
    "    t0 = time.perf_counter()\n",
    "    A_ub, b_ub, c = bauer(kosten, angebot, bedarf)\n",
    "    t_aufbau = time.perf_counter() - t0\n",
    "\n",
    "    t0 = time.perf_counter()\n",
    "    ergebnis = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=(0, None), method=\"highs\")\n",
    "    t_loesen = time.perf_counter() - t0\n",
    "\n",
    "    if not ergebnis.success:\n",
    "        raise RuntimeError(f\"Solver-Status: {ergebnis.message}\")\n",
    "    return t_aufbau, t_loesen, float(ergebnis.fun)\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    print(\"=\" * 88)\n",
    "    print(\"  MODELLAUFBAU: WO DIE ZEIT WIRKLICH HINGEHT\")\n",
    "    print(\"=\" * 88)\n",
    "    print(\"Transportproblem mit m Werken und n Kunden -> m*n Variablen.\\n\")\n",
    "\n",
    "    varianten: list[tuple[str, object]] = [\n",
    "        (\"A: OR-Tools, Add() je NB\", None),        # Sonderfall, eigener Messpfad\n",
    "        (\"B: COO in Schleifen\", baue_mit_schleifen),\n",
    "        (\"C: NumPy vektorisiert\", baue_vektorisiert),\n",
    "    ]\n",
    "    if HAT_POLARS:\n",
    "        varianten.append((\"D: Polars-Tabelle\", baue_mit_polars))\n",
    "    else:\n",
    "        print(\"Hinweis: polars nicht installiert - Variante D wird uebersprungen.\\n\")\n",
    "\n",
    "    # Aufwaermlauf: Der erste Aufruf bezahlt Importe und einmalige\n",
    "    # Initialisierungen. Wer den mitmisst, vergleicht Startkosten statt\n",
    "    # Rechenarbeit - ein klassischer Benchmark-Fehler.\n",
    "    aufwaerm = erzeuge_daten(10, 10)\n",
    "    loese_mit_modellierungsschicht(*aufwaerm)\n",
    "    for _, bauer in varianten[1:]:\n",
    "        bauer(*aufwaerm)\n",
    "\n",
    "    print(f\"{'Groesse':<26} {'Variante':<26} {'Aufbau':>9} {'Loesen':>9} \"\n",
    "          f\"{'Aufbauanteil':>13}\")\n",
    "    print(\"-\" * 88)\n",
    "\n",
    "    for m, n in [(40, 40), (120, 120), (250, 250)]:\n",
    "        kosten, angebot, bedarf = erzeuge_daten(m, n)\n",
    "        zielwerte = []\n",
    "        for nummer, (name, bauer) in enumerate(varianten):\n",
    "            if bauer is None:\n",
    "                t_aufbau, t_loesen, ziel = loese_mit_modellierungsschicht(\n",
    "                    kosten, angebot, bedarf)\n",
    "            else:\n",
    "                t_aufbau, t_loesen, ziel = messe_matrixvariante(\n",
    "                    bauer, kosten, angebot, bedarf)\n",
    "            zielwerte.append(ziel)\n",
    "            anteil = 100.0 * t_aufbau / (t_aufbau + t_loesen)\n",
    "            groesse = f\"{m}x{n} = {m*n:,} Variablen\" if nummer == 0 else \"\"\n",
    "            print(f\"{groesse:<26} {name:<26} {t_aufbau:>8.3f}s {t_loesen:>8.3f}s \"\n",
    "                  f\"{anteil:>12.0f} %\")\n",
    "\n",
    "        # Alle Varianten muessen dasselbe Problem beschreiben.\n",
    "        spanne = max(zielwerte) - min(zielwerte)\n",
    "        assert spanne < 1e-6 * max(abs(z) for z in zielwerte), \\\n",
    "            f\"Varianten widersprechen sich: {zielwerte}\"\n",
    "        print(f\"{'':<26} {'-> Zielwert (alle gleich)':<26} {zielwerte[0]:>9.2f}\"\n",
    "              f\"   Spanne {spanne:.1e}\")\n",
    "        print(\"-\" * 88)\n",
    "\n",
    "    print(\"\\nZwei Lehren aus der Tabelle:\")\n",
    "    print(\"1. Bei Variante A geht mehr Zeit in den AUFBAU als ins Loesen. Wer hier\")\n",
    "    print(\"   einen schnelleren Solver kauft, beschleunigt den kleineren Teil.\")\n",
    "    print(\"2. Zwischen B und C liegt keine andere Mathematik, nur eine andere\")\n",
    "    print(\"   Schreibweise derselben Matrix - Schleife gegen Kronecker-Produkt.\")\n",
    "    print(\"\\nDie Lesbarkeit von Variante A ist trotzdem viel wert: Fangen Sie dort an,\")\n",
    "    print(\"und vektorisieren Sie erst, wenn die Messung es verlangt.\")\n",
    "    print(\"=\" * 88)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}