{ "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 }