862e92bc7b
Setzt den Isolationsteil von Paket 1 aus Verbesserungen_02.md um. Der Plan
nannte zwei Programme; beim Suchen kam ein drittes dazu, das dasselbe Muster
verwendete.
Ein_System_Vier_Ansaetze.py und Benchmark_Skalierung.py hielten ihre vier
Solvervarianten als Zeichenketten in einem Dictionary und gaben sie an
"python -c" weiter - bei Benchmark_Skalierung.py sogar mit
.format()-Platzhaltern fuer die Instanzgroesse. Aus jeder Variante ist jetzt
eine gewoehnliche Funktion mit lokalem Import geworden.
Solverwechsel_CPSAT_HiGHS.py rief sich selbst ueber sys.argv erneut auf;
auch das entfaellt.
Ausgefuehrt wird ueber einen ProcessPoolExecutor mit zwei Einstellungen, die
beide noetig sind: mp_context "spawn" (frischer Interpreter statt geerbtem
Speicher - unter Linux ist fork der Standard) und max_tasks_per_child=1 (ein
neuer Prozess je Aufgabe; ohne das verwendet der Pool seinen Arbeiter
wieder, und beim zweiten Solver ist der Konflikt zurueck). Nachgemessen:
vier Aufgaben, vier verschiedene PIDs.
Der zweite Punkt hat einen eigenen Warnkasten bekommen, weil der Fehler
leicht zu machen und schwer zu finden ist: Der Absturz kaeme nicht beim
ersten Solver, sondern beim zweiten - und saehe aus wie ein Problem des
zweiten.
Regel 4, dreifach geprueft. Ein_System_Vier_Ansaetze.py: identisch bis auf
die Zeitspalte, einschliesslich der Spannweite 2,41e-08, auf die sich der
Merksatz des Kapitels beruft. Benchmark_Skalierung.py: alle zwoelf
Zielwerte und alle drei Spannweiten bitgleich; Zeiten und Speicher haben
sich verschoben, beide sind im Abdruck seit jeher als hardwareabhaengig
gekennzeichnet. Solverwechsel_CPSAT_HiGHS.py: Ausgabe ohne Zeiten
unveraendert.
Bewusst subprocess bleibt Mutationstest.py: Dort wird pytest auf einer
mutierten Kopie in einem temporaeren Verzeichnis gestartet - ein externes
Werkzeug auf veraenderten Dateien, nicht die Isolation eines Imports.
Neu im Kapitel Oekosystem: ein Abschnitt "Wie die Isolation aussieht, wenn
sie tragen soll" - warum ein Codestring die schlechteste Umsetzung von
"eigener Prozess" ist. Anhang C nennt jetzt ebenfalls ProcessPoolExecutor.
Ein eigener Fehler, gefunden und abgesichert: Ich hatte dem neuen ### ein
{#sec:...}-Label gegeben. ABSCHNITT_RE erkennt nur "## " - das Label waere
nie registriert worden und jeder Verweis darauf ins Leere gelaufen, ohne
Warnung. Label entfernt, --check meldet den Fall jetzt. Gegengetestet.
Stand: 818 Querverweise, 76 Programme, 33 pytest-Tests, PDF 760 Seiten, 69
netzfreie Programme fehlerfrei.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
734 lines
32 KiB
Text
Generated
734 lines
32 KiB
Text
Generated
{
|
|
"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
|
|
}
|