633 lines
28 KiB
Text
633 lines
28 KiB
Text
|
{
|
||
|
|
"cells": [
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"# Kapitel 5: Lineare Programmierung — Simplex, Dualität und Schattenpreise\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": [
|
||
|
|
"## Das Simplex-Tableau in Python\n",
|
||
|
|
"\n",
|
||
|
|
"`Simplex_Tableau_LP.py`\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"#!/usr/bin/env python3\n",
|
||
|
|
"\n",
|
||
|
|
"# Simplex_Tableau_LP.py\n",
|
||
|
|
"\"\"\"\n",
|
||
|
|
"Kapitel LP: Vollständige Implementierung des Simplex-Algorithmus (Tableau-Methode).\n",
|
||
|
|
"\n",
|
||
|
|
"GRENZEN DIESER IMPLEMENTIERUNG (bewusst, aus didaktischen Gründen):\n",
|
||
|
|
" * nur Maximierung\n",
|
||
|
|
" * nur \"<=\"-Nebenbedingungen\n",
|
||
|
|
" * alle b_i >= 0 (sonst wäre der Ursprung keine zulässige Startecke und man\n",
|
||
|
|
" bräuchte eine Phase-1-Rechnung mit künstlichen Variablen)\n",
|
||
|
|
"Für den produktiven Einsatz nimmt man HiGHS - dieser Code dient dem Verständnis.\n",
|
||
|
|
"\n",
|
||
|
|
"Voraussetzungen werden geprüft statt stillschweigend angenommen;\n",
|
||
|
|
"Iterationsprotokoll und Schattenpreise werden ausgegeben.\n",
|
||
|
|
"\"\"\"\n",
|
||
|
|
"\n",
|
||
|
|
"import numpy as np\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"class SimplexTableauSolver:\n",
|
||
|
|
" \"\"\"Maximierungs-Standardform: max c^T x u.d.N. A x <= b, x >= 0, b >= 0.\"\"\"\n",
|
||
|
|
"\n",
|
||
|
|
" def __init__(self, c, A, b, variablennamen=None, restriktionsnamen=None):\n",
|
||
|
|
" self.c = np.asarray(c, dtype=float)\n",
|
||
|
|
" self.A = np.asarray(A, dtype=float)\n",
|
||
|
|
" self.b = np.asarray(b, dtype=float)\n",
|
||
|
|
" self.n = len(self.c) # Anzahl Originalvariablen\n",
|
||
|
|
" self.m = len(self.b) # Anzahl Nebenbedingungen\n",
|
||
|
|
"\n",
|
||
|
|
" # --- Voraussetzungen pruefen, statt sie stillschweigend anzunehmen ---\n",
|
||
|
|
" if self.A.shape != (self.m, self.n):\n",
|
||
|
|
" raise ValueError(f\"A hat Form {self.A.shape}, erwartet ({self.m}, {self.n}).\")\n",
|
||
|
|
" if np.any(self.b < 0):\n",
|
||
|
|
" raise ValueError(\n",
|
||
|
|
" \"Mindestens ein b_i ist negativ. Dann ist der Ursprung keine zulässige \"\n",
|
||
|
|
" \"Startecke; dieser Solver benötigt eine Phase-1-Rechnung, die hier \"\n",
|
||
|
|
" \"bewusst nicht implementiert ist. Nutzen Sie scipy.optimize.linprog.\"\n",
|
||
|
|
" )\n",
|
||
|
|
"\n",
|
||
|
|
" self.var_namen = variablennamen or [f\"x{j+1}\" for j in range(self.n)]\n",
|
||
|
|
" self.restr_namen = restriktionsnamen or [f\"R{i+1}\" for i in range(self.m)]\n",
|
||
|
|
"\n",
|
||
|
|
" self.tableau = None\n",
|
||
|
|
" self.basis = None # welche Variable ist in welcher Zeile Basis?\n",
|
||
|
|
" self._baue_starttableau()\n",
|
||
|
|
"\n",
|
||
|
|
" def _baue_starttableau(self):\n",
|
||
|
|
" \"\"\"Zeilen: m Nebenbedingungen + Zielfunktionszeile.\n",
|
||
|
|
" Spalten: n Variablen + m Schlupfvariablen + rechte Seite.\"\"\"\n",
|
||
|
|
" m, n = self.m, self.n\n",
|
||
|
|
" self.tableau = np.zeros((m + 1, n + m + 1))\n",
|
||
|
|
" self.tableau[:m, :n] = self.A # Koeffizienten\n",
|
||
|
|
" self.tableau[:m, n:n + m] = np.eye(m) # Schlupfvariablen\n",
|
||
|
|
" self.tableau[:m, -1] = self.b # rechte Seite\n",
|
||
|
|
" self.tableau[-1, :n] = -self.c # Zielzeile: -c (Maximierung)\n",
|
||
|
|
" self.basis = list(range(n, n + m)) # Start: alle Schlupf in der Basis\n",
|
||
|
|
"\n",
|
||
|
|
" def _spaltenname(self, index):\n",
|
||
|
|
" return self.var_namen[index] if index < self.n else f\"s{index - self.n + 1}\"\n",
|
||
|
|
"\n",
|
||
|
|
" def solve(self, max_iterationen=100, protokoll=True):\n",
|
||
|
|
" m, n = self.m, self.n\n",
|
||
|
|
" if protokoll:\n",
|
||
|
|
" print(f\"{'Iter':>4} | {'eintritt':>9} | {'austritt':>9} | \"\n",
|
||
|
|
" f\"{'Pivot':>8} | {'Z':>12}\")\n",
|
||
|
|
" print(\"-\" * 58)\n",
|
||
|
|
"\n",
|
||
|
|
" for iteration in range(1, max_iterationen + 1):\n",
|
||
|
|
" zielzeile = self.tableau[-1, :-1]\n",
|
||
|
|
"\n",
|
||
|
|
" # 1. Optimalitätsprüfung: alle Koeffizienten >= 0 ?\n",
|
||
|
|
" if np.all(zielzeile >= -1e-9):\n",
|
||
|
|
" if protokoll:\n",
|
||
|
|
" print(\"-\" * 58)\n",
|
||
|
|
" print(f\"Optimum nach {iteration - 1} Pivotschritten erreicht.\")\n",
|
||
|
|
" return self._loesung_auslesen()\n",
|
||
|
|
"\n",
|
||
|
|
" # 2. Pivotspalte: negativster Eintrag (Dantzig-Regel)\n",
|
||
|
|
" pivot_spalte = int(np.argmin(zielzeile))\n",
|
||
|
|
"\n",
|
||
|
|
" # 3. Pivotzeile: minimaler Quotient über POSITIVE Spalteneinträge\n",
|
||
|
|
" spalte = self.tableau[:m, pivot_spalte]\n",
|
||
|
|
" rechte_seite = self.tableau[:m, -1]\n",
|
||
|
|
" quotienten = np.where(spalte > 1e-9, rechte_seite / np.where(spalte > 1e-9, spalte, 1),\n",
|
||
|
|
" np.inf)\n",
|
||
|
|
" pivot_zeile = int(np.argmin(quotienten))\n",
|
||
|
|
" if not np.isfinite(quotienten[pivot_zeile]):\n",
|
||
|
|
" raise ValueError(\n",
|
||
|
|
" f\"Problem ist unbeschraenkt: Variable {self._spaltenname(pivot_spalte)} \"\n",
|
||
|
|
" \"kann beliebig wachsen, ohne eine Bedingung zu verletzen. \"\n",
|
||
|
|
" \"Meist fehlt eine Kapazitaetsbeschraenkung.\"\n",
|
||
|
|
" )\n",
|
||
|
|
"\n",
|
||
|
|
" if protokoll:\n",
|
||
|
|
" print(f\"{iteration:>4} | {self._spaltenname(pivot_spalte):>9} | \"\n",
|
||
|
|
" f\"{self._spaltenname(self.basis[pivot_zeile]):>9} | \"\n",
|
||
|
|
" f\"{self.tableau[pivot_zeile, pivot_spalte]:>8.3f} | \"\n",
|
||
|
|
" f\"{self.tableau[-1, -1]:>12,.2f}\")\n",
|
||
|
|
"\n",
|
||
|
|
" # 4. Pivotoperation (Gauß-Jordan)\n",
|
||
|
|
" pivot_wert = self.tableau[pivot_zeile, pivot_spalte]\n",
|
||
|
|
" self.tableau[pivot_zeile, :] /= pivot_wert\n",
|
||
|
|
" for zeile in range(m + 1):\n",
|
||
|
|
" if zeile != pivot_zeile:\n",
|
||
|
|
" faktor = self.tableau[zeile, pivot_spalte]\n",
|
||
|
|
" self.tableau[zeile, :] -= faktor * self.tableau[pivot_zeile, :]\n",
|
||
|
|
"\n",
|
||
|
|
" self.basis[pivot_zeile] = pivot_spalte\n",
|
||
|
|
"\n",
|
||
|
|
" raise RuntimeError(\"Maximale Iterationszahl ueberschritten (moeglicherweise Zyklus).\")\n",
|
||
|
|
"\n",
|
||
|
|
" def _loesung_auslesen(self):\n",
|
||
|
|
" \"\"\"Basisvariablen tragen den RHS-Wert ihrer Zeile, Nichtbasisvariablen sind 0.\"\"\"\n",
|
||
|
|
" x = np.zeros(self.n + self.m)\n",
|
||
|
|
" for zeile, spalte in enumerate(self.basis):\n",
|
||
|
|
" x[spalte] = self.tableau[zeile, -1]\n",
|
||
|
|
" return x[:self.n], x[self.n:], self.tableau[-1, -1]\n",
|
||
|
|
"\n",
|
||
|
|
" def schattenpreise(self):\n",
|
||
|
|
" \"\"\"Die Zielzeile unter den Schlupfspalten enthält direkt die Dualwerte.\"\"\"\n",
|
||
|
|
" return self.tableau[-1, self.n:self.n + self.m].copy()\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"if __name__ == \"__main__\":\n",
|
||
|
|
" # Modell aus der Simplex-Handrechnung (Bot-Beispiel, Kapitel Einfuehrung)\n",
|
||
|
|
" ertraege = [150.0, 250.0]\n",
|
||
|
|
" matrix = [[2.0, 5.0],\n",
|
||
|
|
" [4.0, 6.0],\n",
|
||
|
|
" [1.0, 0.0]]\n",
|
||
|
|
" kapazitaeten = [40.0, 60.0, 8.0]\n",
|
||
|
|
" var_namen = [\"x_A\", \"x_B\"]\n",
|
||
|
|
" restr_namen = [\"vCPU\", \"RAM\", \"Marktlimit\"]\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"=\" * 58)\n",
|
||
|
|
" print(\" SIMPLEX-TABLEAU: ITERATIONSPROTOKOLL\")\n",
|
||
|
|
" print(\"=\" * 58)\n",
|
||
|
|
"\n",
|
||
|
|
" solver = SimplexTableauSolver(ertraege, matrix, kapazitaeten, var_namen, restr_namen)\n",
|
||
|
|
" x_opt, schlupf, z_opt = solver.solve()\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"\\n\" + \"=\" * 58)\n",
|
||
|
|
" print(\" ERGEBNIS\")\n",
|
||
|
|
" print(\"=\" * 58)\n",
|
||
|
|
" for name, wert in zip(var_namen, x_opt):\n",
|
||
|
|
" print(f\" {name:<12} = {wert:8.4f}\")\n",
|
||
|
|
" print(f\" {'Zielwert Z':<12} = {z_opt:8.2f} EUR\")\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"\\n Ressourcenanalyse:\")\n",
|
||
|
|
" print(f\" {'Ressource':<12} {'Schlupf':>9} {'Status':>22} {'Schattenpreis':>15}\")\n",
|
||
|
|
" print(\" \" + \"-\" * 60)\n",
|
||
|
|
" for name, s, y in zip(restr_namen, schlupf, solver.schattenpreise()):\n",
|
||
|
|
" status = \"ENGPASS (bindend)\" if abs(s) < 1e-9 else \"Reserve vorhanden\"\n",
|
||
|
|
" print(f\" {name:<12} {s:>9.3f} {status:>22} {y:>12.2f} EUR\")\n",
|
||
|
|
"\n",
|
||
|
|
" # Selbstkontrolle: komplementaerer Schlupf muss gelten\n",
|
||
|
|
" for s, y in zip(schlupf, solver.schattenpreise()):\n",
|
||
|
|
" assert abs(s * y) < 1e-6, \"Komplementaerer Schlupf verletzt - Rechenfehler!\"\n",
|
||
|
|
" print(\"\\n Pruefung: komplementaerer Schlupf (s_i * y_i = 0) fuer alle i erfuellt.\")\n",
|
||
|
|
" print(\"=\" * 58)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## Der starke Dualitätssatz\n",
|
||
|
|
"\n",
|
||
|
|
"`Dualitaet_Nachweis.py`\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"#!/usr/bin/env python3\n",
|
||
|
|
"\n",
|
||
|
|
"# Dualitaet_Nachweis.py\n",
|
||
|
|
"\"\"\"\n",
|
||
|
|
"Kapitel LP: Primales und duales Problem unabhaengig loesen und den starken\n",
|
||
|
|
"Dualitaetssatz sowie den komplementaeren Schlupf numerisch nachweisen.\n",
|
||
|
|
"\n",
|
||
|
|
"Modell aus der Simplex-Handrechnung (Bot-Allokation, LP-Relaxation):\n",
|
||
|
|
" max 150*x1 + 250*x2 u.d.N. 2*x1+5*x2<=40, 4*x1+6*x2<=60, x1<=8, x>=0\n",
|
||
|
|
"\"\"\"\n",
|
||
|
|
"\n",
|
||
|
|
"import numpy as np\n",
|
||
|
|
"from scipy.optimize import linprog\n",
|
||
|
|
"\n",
|
||
|
|
"# --- Primales Problem --------------------------------------------------------\n",
|
||
|
|
"C_PRIMAL = np.array([150.0, 250.0])\n",
|
||
|
|
"A_PRIMAL = np.array([[2.0, 5.0], [4.0, 6.0], [1.0, 0.0]])\n",
|
||
|
|
"B_PRIMAL = np.array([40.0, 60.0, 8.0])\n",
|
||
|
|
"\n",
|
||
|
|
"# --- Duales Problem: min b^T y u.d.N. A^T y >= c, y >= 0 ------------------\n",
|
||
|
|
"# linprog kennt nur <=, also A^T y >= c <=> -A^T y <= -c\n",
|
||
|
|
"C_DUAL = B_PRIMAL\n",
|
||
|
|
"A_DUAL = -A_PRIMAL.T\n",
|
||
|
|
"B_DUAL = -C_PRIMAL\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def loese():\n",
|
||
|
|
" primal = linprog(c=-C_PRIMAL, A_ub=A_PRIMAL, b_ub=B_PRIMAL,\n",
|
||
|
|
" bounds=[(0, None)] * 2, method=\"highs\")\n",
|
||
|
|
" dual = linprog(c=C_DUAL, A_ub=A_DUAL, b_ub=B_DUAL,\n",
|
||
|
|
" bounds=[(0, None)] * 3, method=\"highs\")\n",
|
||
|
|
" if not primal.success or not dual.success:\n",
|
||
|
|
" raise SystemExit(\"Primal oder Dual nicht loesbar.\")\n",
|
||
|
|
" return primal, dual\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"if __name__ == \"__main__\":\n",
|
||
|
|
" primal, dual = loese()\n",
|
||
|
|
" x = primal.x\n",
|
||
|
|
" y = dual.x\n",
|
||
|
|
" z_primal = -primal.fun\n",
|
||
|
|
" z_dual = dual.fun\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"=\" * 78)\n",
|
||
|
|
" print(\" PRIMALES PROBLEM\")\n",
|
||
|
|
" print(\"=\" * 78)\n",
|
||
|
|
" print(f\"x* = ({x[0]:.4f}, {x[1]:.4f})\")\n",
|
||
|
|
" print(f\"Z* = {z_primal:.4f}\")\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"\\n\" + \"=\" * 78)\n",
|
||
|
|
" print(\" DUALES PROBLEM\")\n",
|
||
|
|
" print(\"=\" * 78)\n",
|
||
|
|
" print(f\"y* = ({y[0]:.4f}, {y[1]:.4f}, {y[2]:.4f})\")\n",
|
||
|
|
" print(f\"W* = {z_dual:.4f}\")\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"\\n\" + \"=\" * 78)\n",
|
||
|
|
" print(\" STARKER DUALITAETSSATZ: c^T x* == b^T y* ?\")\n",
|
||
|
|
" print(\"=\" * 78)\n",
|
||
|
|
" print(f\" Primal Z* = {z_primal:.6f}\")\n",
|
||
|
|
" print(f\" Dual W* = {z_dual:.6f}\")\n",
|
||
|
|
" differenz = abs(z_primal - z_dual)\n",
|
||
|
|
" print(f\" Differenz = {differenz:.2e} -> \"\n",
|
||
|
|
" f\"{'BESTAETIGT' if differenz < 1e-6 else 'VERLETZT!'}\")\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"\\n\" + \"=\" * 78)\n",
|
||
|
|
" print(\" KOMPLEMENTAERER SCHLUPF: s_i * y_i == 0 fuer alle i ?\")\n",
|
||
|
|
" print(\"=\" * 78)\n",
|
||
|
|
" schlupf = B_PRIMAL - A_PRIMAL @ x\n",
|
||
|
|
" ressourcen = [\"vCPU (s1)\", \"RAM (s2)\", \"Marktlimit (s3)\"]\n",
|
||
|
|
" for name, s, yi in zip(ressourcen, schlupf, y):\n",
|
||
|
|
" produkt = s * yi\n",
|
||
|
|
" print(f\" {name:<16} Schlupf s={s:6.4f} Schattenpreis y={yi:6.4f} \"\n",
|
||
|
|
" f\"s*y={produkt:.2e} {'OK' if abs(produkt) < 1e-6 else 'VERLETZT!'}\")\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"\\nFazit: Das dual geloeste y* stimmt exakt mit den Schattenpreisen\")\n",
|
||
|
|
" print(\"überein, die die Simplex-Rechnung von Hand in der Z-Zeile\")\n",
|
||
|
|
" print(\"ablas - unabhaengig voneinander berechnet, identisches Ergebnis.\")\n",
|
||
|
|
" print(\"=\" * 78)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## Praxisfall: Sensitivitätsanalyse mit korrekten Schattenpreisen\n",
|
||
|
|
"\n",
|
||
|
|
"`Sensitivitaetsanalyse.py`\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"#!/usr/bin/env python3\n",
|
||
|
|
"\n",
|
||
|
|
"# Sensitivitaetsanalyse.py\n",
|
||
|
|
"\"\"\"\n",
|
||
|
|
"Kapitel LP: Schattenpreis- und Sensitivitätsanalyse mit SciPy und HiGHS.\n",
|
||
|
|
"\n",
|
||
|
|
"Achtung: Ohne Vorzeichenumkehr der Dualwerte waere die Handlungsempfehlung\n",
|
||
|
|
"strukturell immer \"kein Zukauf noetig\" - selbst bei harten Engpaessen. Dieses\n",
|
||
|
|
"Programm zeigt die korrekte Vorzeichenbehandlung.\n",
|
||
|
|
"\"\"\"\n",
|
||
|
|
"\n",
|
||
|
|
"import numpy as np\n",
|
||
|
|
"from scipy.optimize import linprog\n",
|
||
|
|
"\n",
|
||
|
|
"# --- Modell: Maximiere Deckungsbeitrag aus 3 Produkten ---------------------\n",
|
||
|
|
"# max 40*x1 + 30*x2 + 50*x3\n",
|
||
|
|
"DECKUNGSBEITRAG = np.array([40.0, 30.0, 50.0])\n",
|
||
|
|
"\n",
|
||
|
|
"# Ressourcenverbrauch je Produkt (Zeile = Ressource, Spalte = Produkt)\n",
|
||
|
|
"VERBRAUCH = np.array([\n",
|
||
|
|
" [2.0, 1.0, 3.0], # Montagezeit\n",
|
||
|
|
" [1.0, 2.0, 1.0], # Lackierzeit\n",
|
||
|
|
" [1.0, 0.5, 2.0], # Qualitätsprüfung\n",
|
||
|
|
"])\n",
|
||
|
|
"KAPAZITAET = np.array([120.0, 80.0, 50.0]) # Stunden\n",
|
||
|
|
"RESSOURCEN = [\"Montage\", \"Lackieren\", \"Qualitätsprüfung\"]\n",
|
||
|
|
"PRODUKTE = [\"Produkt 1\", \"Produkt 2\", \"Produkt 3\"]\n",
|
||
|
|
"\n",
|
||
|
|
"ANGEBOTSPREIS_PRUEFSTUNDE = 18.0 # EUR/h - lohnt sich der Zukauf?\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def analysiere():\n",
|
||
|
|
" # linprog MINIMIERT -> Zielfunktion negieren\n",
|
||
|
|
" ergebnis = linprog(c=-DECKUNGSBEITRAG, A_ub=VERBRAUCH, b_ub=KAPAZITAET,\n",
|
||
|
|
" bounds=[(0, None)] * len(DECKUNGSBEITRAG), method=\"highs\")\n",
|
||
|
|
" if not ergebnis.success:\n",
|
||
|
|
" raise SystemExit(f\"Kein Optimum gefunden: {ergebnis.message}\")\n",
|
||
|
|
"\n",
|
||
|
|
" max_gewinn = -ergebnis.fun\n",
|
||
|
|
" mengen = ergebnis.x\n",
|
||
|
|
" schlupf = ergebnis.slack\n",
|
||
|
|
"\n",
|
||
|
|
" # ------------------------------------------------------------------\n",
|
||
|
|
" # DER ENTSCHEIDENDE PUNKT:\n",
|
||
|
|
" # Weil wir zur Maximierung negiert haben, sind die Dualwerte aus\n",
|
||
|
|
" # linprog fuer \"<=\"-Bedingungen <= 0. Zurueckdrehen!\n",
|
||
|
|
" # ------------------------------------------------------------------\n",
|
||
|
|
" schattenpreise = -ergebnis.ineqlin.marginals\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"=\" * 74)\n",
|
||
|
|
" print(\" PRIMALE UND DUALE ERGEBNISANALYSE (SENSITIVITAET)\")\n",
|
||
|
|
" print(\"=\" * 74)\n",
|
||
|
|
" print(f\"Maximaler Deckungsbeitrag: {max_gewinn:,.2f} EUR\\n\")\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"--- Primalloesung: optimale Produktionsmengen ---\")\n",
|
||
|
|
" for name, menge in zip(PRODUKTE, mengen):\n",
|
||
|
|
" print(f\" * {name}: {menge:8.2f} Stueck\")\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"\\n--- Duale Analyse: Schattenpreise und Auslastung ---\")\n",
|
||
|
|
" for i, name in enumerate(RESSOURCEN):\n",
|
||
|
|
" kapazitaet = KAPAZITAET[i]\n",
|
||
|
|
" genutzt = kapazitaet - schlupf[i]\n",
|
||
|
|
" auslastung = genutzt / kapazitaet * 100\n",
|
||
|
|
" preis = schattenpreise[i]\n",
|
||
|
|
" bindend = abs(schlupf[i]) < 1e-9\n",
|
||
|
|
"\n",
|
||
|
|
" print(f\"\\nRessource '{name}':\")\n",
|
||
|
|
" print(f\" Auslastung: {genutzt:6.1f} / {kapazitaet:6.1f} h ({auslastung:5.1f} %)\"\n",
|
||
|
|
" f\" -> {'ENGPASS' if bindend else 'Reserve: %.1f h' % schlupf[i]}\")\n",
|
||
|
|
" print(f\" Schattenpreis: {preis:6.2f} EUR je zusaetzlicher Stunde\")\n",
|
||
|
|
"\n",
|
||
|
|
" if preis > 1e-9:\n",
|
||
|
|
" print(f\" >> Zusaetzliche Stunden lohnen sich bis zu einem Preis von \"\n",
|
||
|
|
" f\"{preis:.2f} EUR/h.\")\n",
|
||
|
|
" else:\n",
|
||
|
|
" print(f\" >> Kein Zukauf noetig - die Kapazitaet ist nicht erschoepft.\")\n",
|
||
|
|
"\n",
|
||
|
|
" # --- Konkrete Kaufentscheidung ---------------------------------------\n",
|
||
|
|
" preis_pruefung = schattenpreise[RESSOURCEN.index(\"Qualitätsprüfung\")]\n",
|
||
|
|
" marge = preis_pruefung - ANGEBOTSPREIS_PRUEFSTUNDE\n",
|
||
|
|
" print(\"\\n\" + \"-\" * 74)\n",
|
||
|
|
" print(f\"ENTSCHEIDUNG: Pruefstunden werden fuer \"\n",
|
||
|
|
" f\"{ANGEBOTSPREIS_PRUEFSTUNDE:.2f} EUR/h angeboten.\")\n",
|
||
|
|
" print(f\" Schattenpreis: {preis_pruefung:6.2f} EUR/h\")\n",
|
||
|
|
" print(f\" Angebotspreis: {ANGEBOTSPREIS_PRUEFSTUNDE:6.2f} EUR/h\")\n",
|
||
|
|
" print(f\" Marge: {marge:+6.2f} EUR je zugekaufter Stunde\")\n",
|
||
|
|
" print(f\" >> {'ZUKAUFEN' if marge > 0 else 'NICHT ZUKAUFEN'}\")\n",
|
||
|
|
"\n",
|
||
|
|
" # --- Numerische Gegenprobe: Kapazitaet wirklich um 1 erhoehen --------\n",
|
||
|
|
" kapazitaet_plus = KAPAZITAET.copy()\n",
|
||
|
|
" kapazitaet_plus[RESSOURCEN.index(\"Qualitätsprüfung\")] += 1.0\n",
|
||
|
|
" gegenprobe = linprog(c=-DECKUNGSBEITRAG, A_ub=VERBRAUCH, b_ub=kapazitaet_plus,\n",
|
||
|
|
" bounds=[(0, None)] * 3, method=\"highs\")\n",
|
||
|
|
" tatsaechlicher_zuwachs = -gegenprobe.fun - max_gewinn\n",
|
||
|
|
" print(\"\\n--- Gegenprobe: Modell mit +1 Pruefstunde neu geloest ---\")\n",
|
||
|
|
" print(f\" Vorhergesagt (Schattenpreis): {preis_pruefung:8.4f} EUR\")\n",
|
||
|
|
" print(f\" Tatsaechlich gemessen: {tatsaechlicher_zuwachs:8.4f} EUR\")\n",
|
||
|
|
" assert abs(tatsaechlicher_zuwachs - preis_pruefung) < 1e-6, \\\n",
|
||
|
|
" \"Schattenpreis stimmt nicht mit der Messung ueberein!\"\n",
|
||
|
|
" print(\" -> Der Schattenpreis ist bestaetigt.\")\n",
|
||
|
|
"\n",
|
||
|
|
" # --- Komplementaerer Schlupf pruefen --------------------------------\n",
|
||
|
|
" for s, y in zip(schlupf, schattenpreise):\n",
|
||
|
|
" assert abs(s * y) < 1e-6, \"Komplementaerer Schlupf verletzt!\"\n",
|
||
|
|
" print(\"\\nPruefung: komplementaerer Schlupf fuer alle Ressourcen erfuellt.\")\n",
|
||
|
|
" print(\"=\" * 74)\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"if __name__ == \"__main__\":\n",
|
||
|
|
" analysiere()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## Die ehrliche Auskunft: Schattenpreis-Spannen\n",
|
||
|
|
"\n",
|
||
|
|
"`Toleranzen_und_Entartung.py`\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"#!/usr/bin/env python3\n",
|
||
|
|
"\n",
|
||
|
|
"# Toleranzen_und_Entartung.py\n",
|
||
|
|
"\"\"\"\n",
|
||
|
|
"Kapitel LP: Zwei Faelle, in denen man Schattenpreisen NICHT trauen darf.\n",
|
||
|
|
"\n",
|
||
|
|
" 1. Entartung (Degeneriertheit): Mehr Nebenbedingungen sind aktiv, als das\n",
|
||
|
|
" Problem Variablen hat. Dann ist der Schattenpreis nicht eindeutig - zwei\n",
|
||
|
|
" korrekte Solver liefern voellig verschiedene Werte fuer dasselbe Optimum.\n",
|
||
|
|
" 2. Toleranzen: \"ausgelastet\" heisst nie 'schlupf == 0', sondern immer\n",
|
||
|
|
" 'schlupf < toleranz'. Wer auf exakte Gleichheit prueft, baut Berichte,\n",
|
||
|
|
" die zufaellig mal stimmen und mal nicht.\n",
|
||
|
|
"\n",
|
||
|
|
"Teil 3 zeigt die professionelle Antwort auf Fall 1: Statt EINEN Schattenpreis\n",
|
||
|
|
"zu melden, berechnet man seine SPANNE ueber alle optimalen Dualloesungen.\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",
|
||
|
|
"# Ein bewusst entartetes Beispiel: drei Geraden schneiden sich in EINEM Punkt.\n",
|
||
|
|
"# max x1 + x2\n",
|
||
|
|
"# u.d.N. x1 + 2*x2 <= 4 (A)\n",
|
||
|
|
"# 2*x1 + x2 <= 4 (B)\n",
|
||
|
|
"# x1 + x2 <= 8/3 (C) - laeuft genau durch die Ecke (4/3, 4/3)\n",
|
||
|
|
"# In zwei Dimensionen legen schon zwei Geraden eine Ecke fest. Hier sind drei\n",
|
||
|
|
"# aktiv - eine zu viel. Genau das ist Entartung.\n",
|
||
|
|
"C_ZIEL = np.array([-1.0, -1.0]) # linprog minimiert -> negiert\n",
|
||
|
|
"A_UB = np.array([[1.0, 2.0],\n",
|
||
|
|
" [2.0, 1.0],\n",
|
||
|
|
" [1.0, 1.0]])\n",
|
||
|
|
"B_UB = np.array([4.0, 4.0, 8.0 / 3.0])\n",
|
||
|
|
"NAMEN = [\"A: Fraeszeit\", \"B: Schleifzeit\", \"C: Pruefzeit\"]\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def zeige_entartung() -> float:\n",
|
||
|
|
" \"\"\"Loest dasselbe LP mit zwei Verfahren und vergleicht die Dualwerte.\"\"\"\n",
|
||
|
|
" print(\"=\" * 78)\n",
|
||
|
|
" print(\" 1. ENTARTUNG: DERSELBE PLAN, GEGENSAETZLICHE SCHATTENPREISE\")\n",
|
||
|
|
" print(\"=\" * 78)\n",
|
||
|
|
"\n",
|
||
|
|
" print(f\"{'Verfahren':<26} {'x1':>7} {'x2':>7} {'Z*':>9} Schattenpreise\")\n",
|
||
|
|
" print(\"-\" * 78)\n",
|
||
|
|
"\n",
|
||
|
|
" dualwerte = {}\n",
|
||
|
|
" for verfahren, beschreibung in [(\"highs-ds\", \"Dual Simplex\"),\n",
|
||
|
|
" (\"highs-ipm\", \"Innere-Punkte-Verfahren\")]:\n",
|
||
|
|
" ergebnis = linprog(C_ZIEL, A_ub=A_UB, b_ub=B_UB, bounds=(0, None),\n",
|
||
|
|
" method=verfahren)\n",
|
||
|
|
" if not ergebnis.success:\n",
|
||
|
|
" raise RuntimeError(f\"{verfahren}: {ergebnis.message}\")\n",
|
||
|
|
" y = -ergebnis.ineqlin.marginals\n",
|
||
|
|
" dualwerte[verfahren] = y\n",
|
||
|
|
" print(f\"{beschreibung:<26} {ergebnis.x[0]:>7.3f} {ergebnis.x[1]:>7.3f} \"\n",
|
||
|
|
" f\"{-ergebnis.fun:>9.4f} {np.round(y, 4)}\")\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"-\" * 78)\n",
|
||
|
|
" print(\"Beide Zeilen sind RICHTIG: gleicher Plan, gleicher Zielwert, und beide\")\n",
|
||
|
|
" print(\"Dualvektoren erfuellen die Optimalitaetsbedingungen. Trotzdem sagen sie\")\n",
|
||
|
|
" print(\"das Gegenteil:\")\n",
|
||
|
|
" print(f\" Dual Simplex : {NAMEN[2]} ist wertlos, A und B sind je 0,33 EUR wert.\")\n",
|
||
|
|
" print(f\" Innere Punkte: {NAMEN[0]} und {NAMEN[1]} sind wertlos, C ist 1,00 EUR wert.\")\n",
|
||
|
|
" print()\n",
|
||
|
|
" print(\"Wer auf dieser Grundlage eine Maschine kauft, hat eine 50:50-Chance -\")\n",
|
||
|
|
" print(\"abhaengig davon, welches Verfahren der Solver zufaellig gewaehlt hat.\")\n",
|
||
|
|
"\n",
|
||
|
|
" return float(-linprog(C_ZIEL, A_ub=A_UB, b_ub=B_UB, bounds=(0, None)).fun)\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def entartung_erkennen() -> None:\n",
|
||
|
|
" \"\"\"Der Test, der in jedes Auswertungsskript gehoert.\"\"\"\n",
|
||
|
|
" print(\"\\n\" + \"=\" * 78)\n",
|
||
|
|
" print(\" 2. ENTARTUNG ERKENNEN - UND WARUM '== 0' DABEI VERSAGT\")\n",
|
||
|
|
" print(\"=\" * 78)\n",
|
||
|
|
"\n",
|
||
|
|
" ergebnis = linprog(C_ZIEL, A_ub=A_UB, b_ub=B_UB, bounds=(0, None))\n",
|
||
|
|
" schlupf = B_UB - A_UB @ ergebnis.x\n",
|
||
|
|
"\n",
|
||
|
|
" print(f\"{'Nebenbedingung':<18} {'Schlupf':>16} {'== 0 ?':>9} \"\n",
|
||
|
|
" f\"{'< 1e-7 ?':>10}\")\n",
|
||
|
|
" print(\"-\" * 78)\n",
|
||
|
|
" for name, s in zip(NAMEN, schlupf):\n",
|
||
|
|
" print(f\"{name:<18} {s:>16.3e} {str(s == 0.0):>9} {str(abs(s) < 1e-7):>10}\")\n",
|
||
|
|
"\n",
|
||
|
|
" aktiv = int((np.abs(schlupf) < 1e-7).sum())\n",
|
||
|
|
" variablen = A_UB.shape[1]\n",
|
||
|
|
" print(\"-\" * 78)\n",
|
||
|
|
" print(f\"Aktive Nebenbedingungen: {aktiv}, Variablen: {variablen}\")\n",
|
||
|
|
" if aktiv > variablen:\n",
|
||
|
|
" print(f\"=> ENTARTET. {aktiv} aktive Restriktionen bei nur {variablen} \"\n",
|
||
|
|
" \"Variablen bedeuten:\")\n",
|
||
|
|
" print(\" Der Schattenpreis ist nicht eindeutig. Melden Sie eine Spanne,\")\n",
|
||
|
|
" print(\" keinen Einzelwert (siehe Teil 3).\")\n",
|
||
|
|
" else:\n",
|
||
|
|
" print(\"=> nicht entartet, die Dualwerte sind eindeutig.\")\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"\\nBeachten Sie die Spalte '== 0': Ein Schlupf von 4.44e-16 ist\")\n",
|
||
|
|
" print(\"rechnerisch null, aber nicht gleich 0.0. Wer mit '==' prueft,\")\n",
|
||
|
|
" print(\"uebersieht genau die Engpaesse, die er sucht.\")\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def schattenpreis_spanne(zielwert: float, toleranz: float = 1e-9\n",
|
||
|
|
" ) -> list[tuple[float, float]]:\n",
|
||
|
|
" \"\"\"Berechnet fuer jede Nebenbedingung die Spanne ihres Schattenpreises\n",
|
||
|
|
" ueber ALLE optimalen Dualloesungen.\n",
|
||
|
|
"\n",
|
||
|
|
" Die Menge der optimalen Dualloesungen ist selbst ein Polyeder:\n",
|
||
|
|
"\n",
|
||
|
|
" A^T y >= c, y >= 0, b^T y = Z*\n",
|
||
|
|
"\n",
|
||
|
|
" (Dualzulaessigkeit plus starker Dualitaetssatz.) Minimiert und maximiert\n",
|
||
|
|
" man darauf y_i, erhaelt man die exakten Grenzen. Das ist die ehrliche\n",
|
||
|
|
" Auskunft an das Management: nicht 'die Stunde ist 0,33 EUR wert', sondern\n",
|
||
|
|
" 'zwischen 0,00 und 0,33 EUR - der Wert ist aus dem Modell nicht bestimmbar'.\n",
|
||
|
|
" \"\"\"\n",
|
||
|
|
" m = A_UB.shape[0]\n",
|
||
|
|
" # A^T y >= c <=> -A^T y <= -c ; Ziel war max c^T x, in linprog-Notation\n",
|
||
|
|
" # steckt c mit negativem Vorzeichen in C_ZIEL.\n",
|
||
|
|
" c_original = -C_ZIEL\n",
|
||
|
|
" A_dual_ub = -A_UB.T\n",
|
||
|
|
" b_dual_ub = -c_original\n",
|
||
|
|
"\n",
|
||
|
|
" spannen = []\n",
|
||
|
|
" for i in range(m):\n",
|
||
|
|
" richtung = np.zeros(m)\n",
|
||
|
|
" richtung[i] = 1.0\n",
|
||
|
|
" grenzen = []\n",
|
||
|
|
" for vorzeichen in (1.0, -1.0): # 1 = minimieren, -1 = maximieren\n",
|
||
|
|
" ergebnis = linprog(\n",
|
||
|
|
" vorzeichen * richtung,\n",
|
||
|
|
" A_ub=A_dual_ub, b_ub=b_dual_ub,\n",
|
||
|
|
" A_eq=B_UB.reshape(1, -1), b_eq=[zielwert],\n",
|
||
|
|
" bounds=(0, None))\n",
|
||
|
|
" if not ergebnis.success:\n",
|
||
|
|
" raise RuntimeError(f\"Spannenberechnung fehlgeschlagen: \"\n",
|
||
|
|
" f\"{ergebnis.message}\")\n",
|
||
|
|
" grenzen.append(float(ergebnis.x[i]))\n",
|
||
|
|
" spannen.append((min(grenzen), max(grenzen)))\n",
|
||
|
|
" return spannen\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def zeige_spanne(zielwert: float) -> None:\n",
|
||
|
|
" print(\"\\n\" + \"=\" * 78)\n",
|
||
|
|
" print(\" 3. DIE EHRLICHE AUSKUNFT: SCHATTENPREIS-SPANNEN\")\n",
|
||
|
|
" print(\"=\" * 78)\n",
|
||
|
|
"\n",
|
||
|
|
" spannen = schattenpreis_spanne(zielwert)\n",
|
||
|
|
" print(f\"{'Nebenbedingung':<18} {'von':>10} {'bis':>10} Aussage\")\n",
|
||
|
|
" print(\"-\" * 78)\n",
|
||
|
|
" for name, (unten, oben) in zip(NAMEN, spannen):\n",
|
||
|
|
" if oben - unten < 1e-7:\n",
|
||
|
|
" aussage = f\"eindeutig {oben:.2f} EUR\"\n",
|
||
|
|
" elif oben < 1e-7:\n",
|
||
|
|
" aussage = \"sicher wertlos (kein Engpass)\"\n",
|
||
|
|
" else:\n",
|
||
|
|
" aussage = \"NICHT bestimmbar - Spanne melden!\"\n",
|
||
|
|
" print(f\"{name:<18} {unten:>10.4f} {oben:>10.4f} {aussage}\")\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"-\" * 78)\n",
|
||
|
|
" print(\"So berichtet man an Entscheider: 'Eine zusaetzliche Fraesstunde ist\")\n",
|
||
|
|
" print(\"zwischen 0,00 und 0,33 EUR wert - das Modell kann es nicht genauer\")\n",
|
||
|
|
" print(\"sagen, weil drei Engpaesse exakt gleichzeitig binden.' Das ist eine\")\n",
|
||
|
|
" print(\"brauchbare Aussage. Ein erfundener Einzelwert ist es nicht.\")\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"if __name__ == \"__main__\":\n",
|
||
|
|
" zielwert = zeige_entartung()\n",
|
||
|
|
" entartung_erkennen()\n",
|
||
|
|
" zeige_spanne(zielwert)\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"\\n\" + \"=\" * 78)\n",
|
||
|
|
" print(\"Merksatz: Pruefen Sie VOR jeder Sensitivitaetsaussage auf Entartung -\")\n",
|
||
|
|
" print(\"und vergleichen Sie Schlupfwerte nie mit '== 0', sondern mit einer\")\n",
|
||
|
|
" print(\"Toleranz.\")\n",
|
||
|
|
" print(\"=\" * 78)"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"metadata": {
|
||
|
|
"kernelspec": {
|
||
|
|
"display_name": "Python 3",
|
||
|
|
"language": "python",
|
||
|
|
"name": "python3"
|
||
|
|
},
|
||
|
|
"language_info": {
|
||
|
|
"name": "python",
|
||
|
|
"version": "3.11"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"nbformat": 4,
|
||
|
|
"nbformat_minor": 5
|
||
|
|
}
|