b7af2f1d9a
Erster Commit des Strangs "Optimierte Entscheidungsfindung mit Python" (Version 04). Die Historie der 71 Commits bis zur Trennung bleibt im uebergeordneten Repository OR_mit_Python liegen, das ab jetzt nur noch Version_03 (eingefroren) verwaltet und Version_04/ ignoriert. Bewusst kein "git subtree split": Der Pfad Version_04/ existiert erst seit der Verzeichnistrennung, ein Split braechte daher nur 7 der 41 einschlaegigen Commits - eine Teilhistorie, die vollstaendig aussieht und es nicht ist. Stand: 5 Teile, 23 Kapitel, 5 Anhaenge, 292 Abschnitte, 703 Querverweise, 325 Indexmarken, 73 Beispielprogramme, 32 SVGs, 4 Plotly-Figuren, 25 Notebooks, PDF mit 715 Seiten. Zusaetzlich in diesem Commit: * pyproject.toml mit Abhaengigkeitsgruppen finance, large-scale, api, figures, dev, empfehlungen. Die abgedruckte requirements.txt bleibt unveraendert daneben bestehen. ortools steht in der Grundausstattung, highspy erst in [large-scale] - so kann der HiGHS-Symbolkonflikt bei der schlanken Installation gar nicht erst auftreten. * Dabei zwei Funde: graphviz wird von erzeuge_architektur_diagramme.py importiert, fehlt aber in requirements.txt (jetzt in [figures]); pymoo steht in requirements.txt, wird aber von keinem Programm importiert, sondern nur im Kapitel Metaheuristiken empfohlen (jetzt in [empfehlungen]). * NEUER_TITEL.md nach Kritik_und_Verbesserungsvorschlaege/ verschoben - es ist die Vorlage des Titelblatts, kein Bestandteil des Werks. Die beiden Fundstellen in PROGRESS.md und erzeuge_titelseite.py nachgezogen. * PROGRESS.md nannte noch den Untertitel der ersten Fassung; auf den tatsaechlichen aus erzeuge_titelseite.py korrigiert. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
579 lines
25 KiB
Text
Generated
579 lines
25 KiB
Text
Generated
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Kapitel 2: Das mathematische Fundament — Vektoren, Matrizen, 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": [
|
|
"## Dasselbe in Python\n",
|
|
"\n",
|
|
"`Matrixform.py`\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"#!/usr/bin/env python3\n",
|
|
"\n",
|
|
"# Matrixform.py\n",
|
|
"\"\"\"\n",
|
|
"Kapitel Fundament: Von der ausgeschriebenen Form zur Matrixform - und zurück.\n",
|
|
"Zeigt, dass beide Schreibweisen dasselbe Modell beschreiben.\n",
|
|
"\"\"\"\n",
|
|
"\n",
|
|
"import numpy as np\n",
|
|
"from scipy.optimize import linprog\n",
|
|
"\n",
|
|
"# --- Modell in Matrixform -------------------------------------------------\n",
|
|
"# max 3*x1 + 5*x2 u.d.N. x1 <= 4, 2*x2 <= 12, 3*x1 + 2*x2 <= 18, x >= 0\n",
|
|
"c = np.array([3.0, 5.0]) # Ertragsvektor (Maximierung)\n",
|
|
"A = np.array([[1.0, 0.0], # Zeile 1: nur x1 kommt vor\n",
|
|
" [0.0, 2.0], # Zeile 2: nur x2 kommt vor\n",
|
|
" [3.0, 2.0]]) # Zeile 3: beide\n",
|
|
"b = np.array([4.0, 12.0, 18.0])\n",
|
|
"namen = [\"Rohstoff A\", \"Rohstoff B\", \"Maschinenzeit\"]\n",
|
|
"\n",
|
|
"# --- Ausgeschriebene Form maschinell erzeugen ------------------------------\n",
|
|
"def zeige_ausgeschrieben(c, A, b, namen):\n",
|
|
" \"\"\"Druckt die Matrixform als lesbares Ungleichungssystem.\"\"\"\n",
|
|
" terme = \" + \".join(f\"{c[j]:g}*x{j+1}\" for j in range(len(c)))\n",
|
|
" print(f\"max {terme}\")\n",
|
|
" print(\"u.d.N.\")\n",
|
|
" for i in range(A.shape[0]):\n",
|
|
" summanden = \" + \".join(f\"{A[i, j]:g}*x{j+1}\"\n",
|
|
" for j in range(A.shape[1]) if A[i, j] != 0)\n",
|
|
" print(f\" {summanden:<24} <= {b[i]:>5g} ({namen[i]})\")\n",
|
|
" print(f\" x1, ..., x{len(c)} >= 0\")\n",
|
|
"\n",
|
|
"zeige_ausgeschrieben(c, A, b, namen)\n",
|
|
"\n",
|
|
"# --- Zulässigkeit eines Punktes prüfen ------------------------------------\n",
|
|
"def ist_zulaessig(x, A, b, toleranz=1e-9):\n",
|
|
" \"\"\"Prüft A x <= b und x >= 0 komponentenweise.\"\"\"\n",
|
|
" verbrauch = A @ x # Matrix-Vektor-Produkt: alle Zeilen auf einmal\n",
|
|
" return bool(np.all(verbrauch <= b + toleranz) and np.all(x >= -toleranz))\n",
|
|
"\n",
|
|
"for kandidat in [np.array([2.0, 6.0]), np.array([4.0, 3.0]), np.array([4.0, 6.0])]:\n",
|
|
" zulaessig = ist_zulaessig(kandidat, A, b)\n",
|
|
" zielwert = c @ kandidat\n",
|
|
" verbrauch = A @ kandidat\n",
|
|
" print(f\"\\nx = {kandidat} -> A x = {verbrauch} \"\n",
|
|
" f\"{'zulaessig' if zulaessig else 'UNZULAESSIG'}, Z = {zielwert:g}\")\n",
|
|
"\n",
|
|
"# --- Lösen: linprog minimiert, also c negieren -----------------------------\n",
|
|
"ergebnis = linprog(c=-c, A_ub=A, b_ub=b, bounds=[(0, None)] * len(c), method=\"highs\")\n",
|
|
"print(\"\\n\" + \"-\" * 60)\n",
|
|
"print(f\"Optimale Loesung: x* = {np.round(ergebnis.x, 4)}\")\n",
|
|
"print(f\"Optimaler Wert: Z* = {-ergebnis.fun:g}\")\n",
|
|
"print(\"Hinweis: linprog minimiert, deshalb wurde c negiert und das\")\n",
|
|
"print(\" Ergebnis am Ende wieder mit -1 multipliziert.\")"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Konvexität sichtbar machen\n",
|
|
"\n",
|
|
"`Konvexitaet_Demo.py`\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"#!/usr/bin/env python3\n",
|
|
"\n",
|
|
"# Konvexitaet_Demo.py\n",
|
|
"\"\"\"\n",
|
|
"Kapitel Fundament: Konvexität praktisch erfahrbar machen.\n",
|
|
"Teil 1: Die Sehnen-Bedingung numerisch nachprüfen.\n",
|
|
"Teil 2: Zeigen, warum lokale Suche bei nicht-konvexen Funktionen scheitert.\n",
|
|
"\"\"\"\n",
|
|
"\n",
|
|
"import numpy as np\n",
|
|
"from scipy.optimize import minimize\n",
|
|
"\n",
|
|
"# --- Teil 1: Sehnen-Test ---------------------------------------------------\n",
|
|
"def ist_konvex_numerisch(f, unten, oben, proben=2000, seed=0):\n",
|
|
" \"\"\"\n",
|
|
" Prüft die Konvexitätsdefinition an zufälligen Punktepaaren:\n",
|
|
" f(theta*x + (1-theta)*y) <= theta*f(x) + (1-theta)*f(y) ?\n",
|
|
" Findet Gegenbeispiele - beweist aber KEINE Konvexität.\n",
|
|
" \"\"\"\n",
|
|
" rng = np.random.default_rng(seed)\n",
|
|
" schlimmste_verletzung = 0.0\n",
|
|
" for _ in range(proben):\n",
|
|
" x, y = rng.uniform(unten, oben, size=2)\n",
|
|
" theta = rng.uniform(0.0, 1.0)\n",
|
|
" links = f(theta * x + (1 - theta) * y)\n",
|
|
" rechts = theta * f(x) + (1 - theta) * f(y)\n",
|
|
" schlimmste_verletzung = max(schlimmste_verletzung, links - rechts)\n",
|
|
" return schlimmste_verletzung\n",
|
|
"\n",
|
|
"\n",
|
|
"funktionen = {\n",
|
|
" \"f(x) = x^2 (konvex)\": lambda x: x ** 2,\n",
|
|
" \"f(x) = |x| (konvex)\": lambda x: abs(x),\n",
|
|
" \"f(x) = e^x (konvex)\": lambda x: np.exp(x),\n",
|
|
" \"f(x) = x^3 (NICHT konvex)\": lambda x: x ** 3,\n",
|
|
" \"f(x) = x^2+3sin(3x) (NICHT konvex)\": lambda x: x ** 2 + 3 * np.sin(3 * x),\n",
|
|
"}\n",
|
|
"\n",
|
|
"print(\"=\" * 68)\n",
|
|
"print(\" TEIL 1: SEHNEN-TEST (positive Zahl = Konvexität verletzt)\")\n",
|
|
"print(\"=\" * 68)\n",
|
|
"for name, f in funktionen.items():\n",
|
|
" verletzung = ist_konvex_numerisch(f, -3.0, 3.0)\n",
|
|
" urteil = \"konvex (keine Verletzung gefunden)\" if verletzung < 1e-9 \\\n",
|
|
" else f\"NICHT konvex (Verletzung bis {verletzung:.3f})\"\n",
|
|
" print(f\"{name:<34} -> {urteil}\")\n",
|
|
"\n",
|
|
"# --- Teil 2: Lokale Suche von verschiedenen Startpunkten -------------------\n",
|
|
"def wellige_funktion(x):\n",
|
|
" \"\"\"Nicht-konvex: eine Parabel mit aufmodulierter Welle -> viele lokale Minima.\"\"\"\n",
|
|
" return x[0] ** 2 + 3.0 * np.sin(3.0 * x[0])\n",
|
|
"\n",
|
|
"print(\"\\n\" + \"=\" * 68)\n",
|
|
"print(\" TEIL 2: LOKALE SUCHE BEI NICHT-KONVEXER FUNKTION\")\n",
|
|
"print(\"=\" * 68)\n",
|
|
"print(f\"{'Startpunkt':>12} | {'gefundenes Minimum':>20} | {'Funktionswert':>14}\")\n",
|
|
"print(\"-\" * 68)\n",
|
|
"\n",
|
|
"ergebnisse = []\n",
|
|
"for start in [-3.0, -1.5, 0.0, 1.5, 3.0]:\n",
|
|
" res = minimize(wellige_funktion, x0=[start], method=\"BFGS\")\n",
|
|
" ergebnisse.append((start, res.x[0], res.fun))\n",
|
|
" print(f\"{start:>12.1f} | {res.x[0]:>20.4f} | {res.fun:>14.4f}\")\n",
|
|
"\n",
|
|
"bester = min(ergebnisse, key=lambda t: t[2])\n",
|
|
"print(\"-\" * 68)\n",
|
|
"print(f\"Je nach Startpunkt landet derselbe Algorithmus in \"\n",
|
|
" f\"{len({round(e[1], 3) for e in ergebnisse})} verschiedenen Minima.\")\n",
|
|
"print(f\"Das beste gefundene: x = {bester[1]:.4f} mit f = {bester[2]:.4f} \"\n",
|
|
" f\"(Start bei {bester[0]:.1f})\")\n",
|
|
"print(\"Bei einer KONVEXEN Funktion waeren alle Zeilen identisch --\")\n",
|
|
"print(\"der Startpunkt waere voellig gleichgueltig.\")\n",
|
|
"print(\"=\" * 68)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Geometrische Visualisierung des Lösungsraums\n",
|
|
"\n",
|
|
"`Visualisierung_Loesungsraum.py`\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"#!/usr/bin/env python3\n",
|
|
"\n",
|
|
"# Visualisierung_Loesungsraum.py\n",
|
|
"\"\"\"\n",
|
|
"Kapitel Fundament: Geometrische Visualisierung eines 2D-Optimierungsraums.\n",
|
|
"\n",
|
|
"Ecken werden berechnet, auf Zulässigkeit geprüft, bewertet und eingezeichnet.\n",
|
|
"Zusätzlich wird das ganzzahlige Optimum systematisch bestimmt statt behauptet.\n",
|
|
"\"\"\"\n",
|
|
"\n",
|
|
"import itertools\n",
|
|
"import os\n",
|
|
"\n",
|
|
"import numpy as np\n",
|
|
"import matplotlib\n",
|
|
"matplotlib.use(\"Agg\") # kein Bildschirm nötig\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"\n",
|
|
"OUTPUT_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), \"output\")\n",
|
|
"os.makedirs(OUTPUT_DIR, exist_ok=True)\n",
|
|
"\n",
|
|
"# --- Modell (identisch zum Kapitel Einfuehrung) ----------------------------\n",
|
|
"# max 150*xA + 250*xB\n",
|
|
"# u.d.N. 2*xA + 5*xB <= 40 (vCPU)\n",
|
|
"# 4*xA + 6*xB <= 60 (RAM)\n",
|
|
"# 1*xA + 0*xB <= 8 (Marktliquidität)\n",
|
|
"c = np.array([150.0, 250.0])\n",
|
|
"A = np.array([[2.0, 5.0], [4.0, 6.0], [1.0, 0.0]])\n",
|
|
"b = np.array([40.0, 60.0, 8.0])\n",
|
|
"restriktionsnamen = [\"vCPU\", \"RAM\", \"Marktlimit\"]\n",
|
|
"\n",
|
|
"\n",
|
|
"def ist_zulaessig(punkt, tol=1e-7):\n",
|
|
" return np.all(A @ punkt <= b + tol) and np.all(punkt >= -tol)\n",
|
|
"\n",
|
|
"\n",
|
|
"def berechne_ecken():\n",
|
|
" \"\"\"\n",
|
|
" Ecken = Schnittpunkte je zweier Begrenzungsgeraden, die zulässig sind.\n",
|
|
" Begrenzungen sind die 3 Restriktionen plus die beiden Achsen xA=0, xB=0.\n",
|
|
" \"\"\"\n",
|
|
" geraden = [(A[i], b[i]) for i in range(len(b))]\n",
|
|
" geraden.append((np.array([1.0, 0.0]), 0.0)) # xA = 0\n",
|
|
" geraden.append((np.array([0.0, 1.0]), 0.0)) # xB = 0\n",
|
|
"\n",
|
|
" ecken = []\n",
|
|
" for (n1, d1), (n2, d2) in itertools.combinations(geraden, 2):\n",
|
|
" M = np.array([n1, n2])\n",
|
|
" if abs(np.linalg.det(M)) < 1e-9: # parallel -> kein Schnittpunkt\n",
|
|
" continue\n",
|
|
" p = np.linalg.solve(M, np.array([d1, d2]))\n",
|
|
" if ist_zulaessig(p) and not any(np.allclose(p, e) for e in ecken):\n",
|
|
" ecken.append(p)\n",
|
|
" return np.array(ecken)\n",
|
|
"\n",
|
|
"\n",
|
|
"def bestes_ganzzahliges():\n",
|
|
" \"\"\"Vollständige Suche über das kleine Gitter - hier zulässig, weil winzig.\"\"\"\n",
|
|
" bester_wert, bester_punkt = -np.inf, None\n",
|
|
" for xa in range(0, 21):\n",
|
|
" for xb in range(0, 21):\n",
|
|
" p = np.array([float(xa), float(xb)])\n",
|
|
" if ist_zulaessig(p) and c @ p > bester_wert:\n",
|
|
" bester_wert, bester_punkt = c @ p, p\n",
|
|
" return bester_punkt, bester_wert\n",
|
|
"\n",
|
|
"\n",
|
|
"# --- Analyse ---------------------------------------------------------------\n",
|
|
"ecken = berechne_ecken()\n",
|
|
"werte = ecken @ c\n",
|
|
"reihenfolge = np.argsort(-werte)\n",
|
|
"\n",
|
|
"print(\"=\" * 62)\n",
|
|
"print(\" ECKEN DES ZULÄSSIGEN POLYEDERS (nach Zielwert sortiert)\")\n",
|
|
"print(\"=\" * 62)\n",
|
|
"print(f\"{'x_A':>8} {'x_B':>8} {'Z = 150 xA + 250 xB':>24}\")\n",
|
|
"print(\"-\" * 62)\n",
|
|
"for i in reihenfolge:\n",
|
|
" print(f\"{ecken[i, 0]:>8.2f} {ecken[i, 1]:>8.2f} {werte[i]:>24,.2f} EUR\")\n",
|
|
"\n",
|
|
"lp_punkt, lp_wert = ecken[reihenfolge[0]], werte[reihenfolge[0]]\n",
|
|
"ip_punkt, ip_wert = bestes_ganzzahliges()\n",
|
|
"\n",
|
|
"print(\"-\" * 62)\n",
|
|
"print(f\"Kontinuierliches Optimum (LP): x = ({lp_punkt[0]:.2f}, {lp_punkt[1]:.2f}), \"\n",
|
|
" f\"Z = {lp_wert:,.2f} EUR\")\n",
|
|
"print(f\"Ganzzahliges Optimum (IP): x = ({ip_punkt[0]:.0f}, {ip_punkt[1]:.0f}), \"\n",
|
|
" f\"Z = {ip_wert:,.2f} EUR\")\n",
|
|
"print(f\"Preis der Ganzzahligkeit: {lp_wert - ip_wert:,.2f} EUR \"\n",
|
|
" f\"({(1 - ip_wert / lp_wert) * 100:.2f} %)\")\n",
|
|
"print(\"=\" * 62)\n",
|
|
"\n",
|
|
"# --- Zeichnung -------------------------------------------------------------\n",
|
|
"gitter = np.linspace(0, 15, 400)\n",
|
|
"xa_gitter, xb_gitter = np.meshgrid(gitter, gitter)\n",
|
|
"\n",
|
|
"plt.figure(figsize=(10, 8))\n",
|
|
"\n",
|
|
"# Restriktionsgeraden\n",
|
|
"plt.plot(gitter, (40 - 2 * gitter) / 5, color=\"tab:blue\", lw=2,\n",
|
|
" label=r\"$2x_A + 5x_B \\leq 40$ (vCPU)\")\n",
|
|
"plt.plot(gitter, (60 - 4 * gitter) / 6, color=\"tab:green\", lw=2,\n",
|
|
" label=r\"$4x_A + 6x_B \\leq 60$ (RAM)\")\n",
|
|
"plt.axvline(x=8, color=\"tab:orange\", lw=2, label=r\"$x_A \\leq 8$ (Marktlimit)\")\n",
|
|
"\n",
|
|
"# Zulässiger Bereich\n",
|
|
"maske = ((2 * xa_gitter + 5 * xb_gitter <= 40) & (4 * xa_gitter + 6 * xb_gitter <= 60)\n",
|
|
" & (xa_gitter <= 8) & (xa_gitter >= 0) & (xb_gitter >= 0))\n",
|
|
"plt.imshow(maske.astype(int), extent=(0, 15, 0, 15), origin=\"lower\",\n",
|
|
" cmap=\"Greys\", alpha=0.25, aspect=\"auto\")\n",
|
|
"\n",
|
|
"# Höhenlinien der Zielfunktion\n",
|
|
"Z = 150 * xa_gitter + 250 * xb_gitter\n",
|
|
"hoehen = plt.contour(xa_gitter, xb_gitter, Z, levels=[500, 1000, 1500, 2000, 2375],\n",
|
|
" colors=\"purple\", linestyles=\"--\", alpha=0.7)\n",
|
|
"plt.clabel(hoehen, inline=True, fontsize=9, fmt=\"Z = %1.0f EUR\")\n",
|
|
"\n",
|
|
"# Ecken einzeichnen - jetzt werden sie tatsächlich benutzt\n",
|
|
"plt.scatter(ecken[:, 0], ecken[:, 1], s=70, facecolors=\"white\",\n",
|
|
" edgecolors=\"black\", zorder=4, label=\"Ecken des Polyeders\")\n",
|
|
"for e, w in zip(ecken, werte):\n",
|
|
" plt.annotate(f\"({e[0]:.1f}, {e[1]:.1f})\\nZ={w:,.0f}\", (e[0], e[1]),\n",
|
|
" textcoords=\"offset points\", xytext=(6, 6), fontsize=8)\n",
|
|
"\n",
|
|
"plt.scatter([lp_punkt[0]], [lp_punkt[1]], color=\"purple\", marker=\"D\", s=110, zorder=5,\n",
|
|
" label=f\"LP-Optimum ({lp_punkt[0]:.1f}, {lp_punkt[1]:.1f})\")\n",
|
|
"plt.scatter([ip_punkt[0]], [ip_punkt[1]], color=\"red\", s=170, zorder=6,\n",
|
|
" label=f\"Ganzzahliges Optimum ({ip_punkt[0]:.0f}, {ip_punkt[1]:.0f})\")\n",
|
|
"\n",
|
|
"plt.xlim(0, 12)\n",
|
|
"plt.ylim(0, 10)\n",
|
|
"plt.xlabel(\"Anzahl Arbitrage-Bots ($x_A$)\", fontsize=11)\n",
|
|
"plt.ylabel(\"Anzahl Trendfolge-Bots ($x_B$)\", fontsize=11)\n",
|
|
"plt.title(\"Polyeder des zulässigen Bereichs mit Niveaulinien der Zielfunktion\", fontsize=13)\n",
|
|
"plt.grid(True, linestyle=\":\", alpha=0.6)\n",
|
|
"plt.legend(loc=\"upper right\", framealpha=0.9)\n",
|
|
"plt.tight_layout()\n",
|
|
"ziel = os.path.join(OUTPUT_DIR, \"feasible_region_2d.png\")\n",
|
|
"plt.savefig(ziel, dpi=150)\n",
|
|
"print(f\"Visualisierung gespeichert unter '{ziel}'\")"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Das Experiment\n",
|
|
"\n",
|
|
"`Skalierung_Kondition.py`\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"#!/usr/bin/env python3\n",
|
|
"\n",
|
|
"# Skalierung_Kondition.py\n",
|
|
"\"\"\"\n",
|
|
"Kapitel Fundament: Was die Konditionszahl kappa(A) praktisch bedeutet - und wie man\n",
|
|
"ein schlecht skaliertes Modell wieder gesund rechnet.\n",
|
|
"\n",
|
|
"Drei Experimente:\n",
|
|
"\n",
|
|
" 1. Fehlerverstaerkung: Wie stark schlaegt eine winzige Datenunsicherheit auf\n",
|
|
" die Loesung durch? kappa(A) ist genau die Obergrenze dieses Faktors.\n",
|
|
" 2. Ruiz-Equilibrierung: Ein Modell, in dem Euro-Betraege (1e7) und\n",
|
|
" Tonnen-Angaben (1e-3) in derselben Matrix stehen, wird durch Zeilen- und\n",
|
|
" Spaltenskalierung um Groessenordnungen besser konditioniert.\n",
|
|
" 3. Toleranzen: Warum eine Binaervariable mit dem Wert 0.99999998 niemals\n",
|
|
" mit int() gerundet werden darf.\n",
|
|
"\n",
|
|
"Benoetigt: numpy, scipy\n",
|
|
"\"\"\"\n",
|
|
"\n",
|
|
"from __future__ import annotations\n",
|
|
"\n",
|
|
"import numpy as np\n",
|
|
"from scipy.optimize import linprog\n",
|
|
"\n",
|
|
"RNG = np.random.default_rng(42)\n",
|
|
"\n",
|
|
"\n",
|
|
"# --- Experiment 1: kappa(A) als Fehlerverstaerker ---------------------------\n",
|
|
"\n",
|
|
"def fehlerverstaerkung(A: np.ndarray, versuche: int = 200,\n",
|
|
" stoerung: float = 1e-10) -> tuple[float, float]:\n",
|
|
" \"\"\"Stoert die rechte Seite b relativ um 'stoerung' in zufaellige Richtungen\n",
|
|
" und misst, um welchen Faktor sich der Fehler in der Loesung x vergroessert.\n",
|
|
"\n",
|
|
" Liefert (Median, Maximum) der Verstaerkung. Die Theorie sagt: Das Maximum\n",
|
|
" kann bis kappa(A) betragen - und nur bis dahin.\n",
|
|
" \"\"\"\n",
|
|
" x_wahr = np.ones(A.shape[1])\n",
|
|
" b = A @ x_wahr\n",
|
|
" faktoren = []\n",
|
|
" for _ in range(versuche):\n",
|
|
" richtung = RNG.normal(size=b.size)\n",
|
|
" richtung /= np.linalg.norm(richtung)\n",
|
|
" b_gestoert = b + stoerung * np.linalg.norm(b) * richtung\n",
|
|
" x_gestoert = np.linalg.solve(A, b_gestoert)\n",
|
|
" rel_x = np.linalg.norm(x_gestoert - x_wahr) / np.linalg.norm(x_wahr)\n",
|
|
" faktoren.append(rel_x / stoerung)\n",
|
|
" return float(np.median(faktoren)), float(np.max(faktoren))\n",
|
|
"\n",
|
|
"\n",
|
|
"def zeige_experiment_1() -> None:\n",
|
|
" print(\"=\" * 74)\n",
|
|
" print(\" 1. KONDITIONSZAHL ALS FEHLERVERSTAERKER\")\n",
|
|
" print(\"=\" * 74)\n",
|
|
" print(\"Zwei Gleichungssysteme, beide exakt loesbar mit x = (1, 1).\")\n",
|
|
" print(\"Die rechte Seite wird um relativ 1e-10 gestoert - so viel Unsicherheit\")\n",
|
|
" print(\"steckt in JEDER gemessenen Betriebszahl allemal.\\n\")\n",
|
|
"\n",
|
|
" modelle = {\n",
|
|
" \"gut konditioniert\": np.array([[2.0, 1.0], [1.0, 3.0]]),\n",
|
|
" \"fast parallele Zeilen\": np.array([[1.0, 1.0], [1.0, 1.0 + 1e-8]]),\n",
|
|
" }\n",
|
|
" print(f\"{'Matrix':<24} {'kappa(A)':>12} {'Verst. median':>14} {'Verst. max':>12}\")\n",
|
|
" print(\"-\" * 74)\n",
|
|
" for name, A in modelle.items():\n",
|
|
" median, maximum = fehlerverstaerkung(A)\n",
|
|
" print(f\"{name:<24} {np.linalg.cond(A):>12.2e} {median:>14.2e} {maximum:>12.2e}\")\n",
|
|
"\n",
|
|
" print(\"\\nLesart: Bei der zweiten Matrix wird aus einem Datenfehler in der\")\n",
|
|
" print(\"10. Nachkommastelle ein Loesungsfehler in der 2. Nachkommastelle.\")\n",
|
|
" print(\"Das Modell ist mathematisch korrekt - und praktisch wertlos.\")\n",
|
|
"\n",
|
|
"\n",
|
|
"# --- Experiment 2: Ruiz-Equilibrierung --------------------------------------\n",
|
|
"\n",
|
|
"def ruiz_equilibrierung(A: np.ndarray, durchlaeufe: int = 20\n",
|
|
" ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:\n",
|
|
" \"\"\"Skaliert A iterativ so, dass alle Zeilen- und Spaltenmaxima nahe 1\n",
|
|
" liegen (Ruiz 2001).\n",
|
|
"\n",
|
|
" In jedem Durchlauf wird jede Zeile durch die Wurzel ihres Betragsmaximums\n",
|
|
" geteilt, danach jede Spalte. Das Verfahren konvergiert schnell und braucht\n",
|
|
" keinerlei Wissen ueber die Bedeutung der Zahlen - genau deshalb steckt es\n",
|
|
" in jedem ernsthaften Solver als Vorverarbeitung.\n",
|
|
"\n",
|
|
" Liefert (A_skaliert, zeilenfaktor r, spaltenfaktor c) mit\n",
|
|
" A_skaliert = diag(r) @ A @ diag(c)\n",
|
|
" \"\"\"\n",
|
|
" m, n = A.shape\n",
|
|
" r = np.ones(m)\n",
|
|
" c = np.ones(n)\n",
|
|
" A_s = A.astype(float).copy()\n",
|
|
"\n",
|
|
" for _ in range(durchlaeufe):\n",
|
|
" zeilen_max = np.abs(A_s).max(axis=1)\n",
|
|
" zeilen_max[zeilen_max == 0] = 1.0\n",
|
|
" d_r = 1.0 / np.sqrt(zeilen_max)\n",
|
|
" A_s = d_r[:, None] * A_s\n",
|
|
" r *= d_r\n",
|
|
"\n",
|
|
" spalten_max = np.abs(A_s).max(axis=0)\n",
|
|
" spalten_max[spalten_max == 0] = 1.0\n",
|
|
" d_c = 1.0 / np.sqrt(spalten_max)\n",
|
|
" A_s = A_s * d_c[None, :]\n",
|
|
" c *= d_c\n",
|
|
"\n",
|
|
" return A_s, r, c\n",
|
|
"\n",
|
|
"\n",
|
|
"def zeige_experiment_2() -> None:\n",
|
|
" print(\"\\n\" + \"=\" * 74)\n",
|
|
" print(\" 2. RUIZ-EQUILIBRIERUNG: EINHEITEN GERADERUECKEN\")\n",
|
|
" print(\"=\" * 74)\n",
|
|
" print(\"Ein Produktionsmodell, in dem vier Ressourcen in voellig\")\n",
|
|
" print(\"verschiedenen Einheiten gemessen werden:\")\n",
|
|
" print(\" Zeile 1: Kapitalbindung in Euro (Groessenordnung 1e7)\")\n",
|
|
" print(\" Zeile 2: Katalysatorverbrauch in Tonnen (Groessenordnung 1e-3)\")\n",
|
|
" print(\" Zeile 3: Energie in Wattsekunden (Groessenordnung 1e5)\")\n",
|
|
" print(\" Zeile 4: Ausschussquote als Anteil (Groessenordnung 1e-2)\\n\")\n",
|
|
"\n",
|
|
" # Quadratisch gewaehlt, damit die Loesung eindeutig ist und die\n",
|
|
" # Ruecktransformation unten wirklich etwas beweist.\n",
|
|
" n = 4\n",
|
|
" grundmatrix = RNG.uniform(0.5, 2.0, size=(n, n))\n",
|
|
" einheiten = np.array([1e7, 1e-3, 1e5, 1e-2])\n",
|
|
" A = grundmatrix * einheiten[:, None]\n",
|
|
"\n",
|
|
" A_s, r, c = ruiz_equilibrierung(A)\n",
|
|
"\n",
|
|
" print(f\"{'':<28} {'kappa(A)':>12} {'groesster Eintrag':>18} \"\n",
|
|
" f\"{'kleinster Eintrag':>18}\")\n",
|
|
" print(\"-\" * 74)\n",
|
|
" for name, matrix in [(\"vor der Skalierung\", A), (\"nach Ruiz-Equilibrierung\", A_s)]:\n",
|
|
" betraege = np.abs(matrix)\n",
|
|
" print(f\"{name:<28} {np.linalg.cond(matrix):>12.2e} \"\n",
|
|
" f\"{betraege.max():>18.2e} {betraege.min():>18.2e}\")\n",
|
|
"\n",
|
|
" # Gegenprobe: Das skalierte Modell beschreibt dasselbe Problem. Wer x_s\n",
|
|
" # loest, erhaelt die urspruengliche Loesung durch x = c * x_s.\n",
|
|
" x_wahr = RNG.uniform(1.0, 5.0, size=n)\n",
|
|
" b = A @ x_wahr\n",
|
|
" b_s = r * b\n",
|
|
" x_s = np.linalg.solve(A_s, b_s)\n",
|
|
" x_zurueck = c * x_s\n",
|
|
" print(f\"\\nRuecktransformation x = c * x_s: groesste Abweichung zur wahren \"\n",
|
|
" f\"Loesung {np.abs(x_zurueck - x_wahr).max():.2e}\")\n",
|
|
" print(\"Die Skalierung ist also verlustfrei - sie aendert nur die Zahlen,\")\n",
|
|
" print(\"nicht das Problem.\")\n",
|
|
"\n",
|
|
"\n",
|
|
"# --- Experiment 3: Toleranzen und der int()-Fehler --------------------------\n",
|
|
"\n",
|
|
"def zeige_experiment_3() -> None:\n",
|
|
" print(\"\\n\" + \"=\" * 74)\n",
|
|
" print(\" 3. TOLERANZEN: WARUM int() DIE FALSCHE RUNDUNG IST\")\n",
|
|
" print(\"=\" * 74)\n",
|
|
"\n",
|
|
" # Ein LP, dessen Optimum bei x = 1 liegt, aber vom Solver nur bis auf\n",
|
|
" # seine Toleranz getroffen wird.\n",
|
|
" ergebnis = linprog(c=[-1.0], A_ub=[[1.0]], b_ub=[1.0],\n",
|
|
" bounds=[(0, None)], method=\"highs\")\n",
|
|
" wert = float(ergebnis.x[0])\n",
|
|
" print(f\"Solver liefert x = {wert!r}\")\n",
|
|
"\n",
|
|
" # Typische Werte, wie sie aus MILP-Solvern zurueckkommen.\n",
|
|
" beispiele = [0.99999998, 1.00000002, 0.49999999, 2.9999999]\n",
|
|
" print(f\"\\n{'Solverwert':>14} {'int()':>8} {'round()':>9} {'Kommentar'}\")\n",
|
|
" print(\"-\" * 74)\n",
|
|
" kommentare = {\n",
|
|
" 0.99999998: \"int() macht aus einer JA- eine NEIN-Entscheidung\",\n",
|
|
" 1.00000002: \"hier ginge int() zufaellig gut - Verlass ist keiner\",\n",
|
|
" 0.49999999: \"echt unentschieden: Modell oder Toleranz pruefen!\",\n",
|
|
" 2.9999999: \"3 Maschinen werden zu 2 - der Plan geht nicht auf\",\n",
|
|
" }\n",
|
|
" for wert_b in beispiele:\n",
|
|
" print(f\"{wert_b:>14.8f} {int(wert_b):>8} {round(wert_b):>9} \"\n",
|
|
" f\"{kommentare[wert_b]}\")\n",
|
|
"\n",
|
|
" print(\"\\nRichtige Vorgehensweise: gegen die Solver-Toleranz pruefen,\")\n",
|
|
" print(\"dann erst runden - und den Zweifelsfall melden statt still zu raten.\")\n",
|
|
"\n",
|
|
" def sichere_ganzzahl(wert: float, toleranz: float = 1e-6) -> int:\n",
|
|
" naechste = round(wert)\n",
|
|
" if abs(wert - naechste) > toleranz:\n",
|
|
" raise ValueError(\n",
|
|
" f\"{wert} ist {abs(wert - naechste):.2e} von der naechsten ganzen \"\n",
|
|
" f\"Zahl entfernt - das ist mehr als die Toleranz {toleranz}. \"\n",
|
|
" \"Ganzzahligkeit im Modell pruefen.\")\n",
|
|
" return naechste\n",
|
|
"\n",
|
|
" for wert_b in beispiele:\n",
|
|
" try:\n",
|
|
" print(f\" sichere_ganzzahl({wert_b}) = {sichere_ganzzahl(wert_b)}\")\n",
|
|
" except ValueError as fehler:\n",
|
|
" print(f\" sichere_ganzzahl({wert_b}) -> ValueError: {fehler}\")\n",
|
|
"\n",
|
|
"\n",
|
|
"if __name__ == \"__main__\":\n",
|
|
" zeige_experiment_1()\n",
|
|
" zeige_experiment_2()\n",
|
|
" zeige_experiment_3()\n",
|
|
" print(\"\\n\" + \"=\" * 74)\n",
|
|
" print(\"Merksatz: Skalieren Sie Ihre Daten, BEVOR der Solver sie sieht -\")\n",
|
|
" print(\"und runden Sie Solver-Ergebnisse NIE ohne Toleranzpruefung.\")\n",
|
|
" print(\"=\" * 74)"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Python 3",
|
|
"language": "python",
|
|
"name": "python3"
|
|
},
|
|
"language_info": {
|
|
"name": "python",
|
|
"version": "3.11"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 5
|
|
}
|