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>
1049 lines
47 KiB
Text
1049 lines
47 KiB
Text
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Kapitel 6: Gemischt-ganzzahlige Optimierung — Diskrete Entscheidungen und Branch-and-Bound\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": [
|
|
"## Warum Runden fundamental scheitert\n",
|
|
"\n",
|
|
"`Runden_Gegenbeispiel.py`\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"#!/usr/bin/env python3\n",
|
|
"\n",
|
|
"# Runden_Gegenbeispiel.py\n",
|
|
"\"\"\"\n",
|
|
"Kapitel MILP: Wie schlecht ist Runden wirklich?\n",
|
|
"\n",
|
|
"Ersetzt die unbelegte Behauptung \"20-50 % Verlust\" durch eine Messung ueber\n",
|
|
"viele Zufallsinstanzen.\n",
|
|
"\"\"\"\n",
|
|
"\n",
|
|
"import numpy as np\n",
|
|
"from scipy.optimize import linprog\n",
|
|
"\n",
|
|
"\n",
|
|
"def erzeuge_instanz(n, m, rng):\n",
|
|
" \"\"\"Zufaelliges Rucksack-aehnliches MILP mit kleinen Zahlen.\"\"\"\n",
|
|
" c = rng.integers(3, 20, size=n).astype(float) # Ertraege\n",
|
|
" A = rng.integers(1, 9, size=(m, n)).astype(float) # Verbraeuche\n",
|
|
" b = (A.sum(axis=1) * rng.uniform(0.25, 0.45)).round() # knappe Kapazitaeten\n",
|
|
" return c, A, b\n",
|
|
"\n",
|
|
"\n",
|
|
"def loese_lp(c, A, b, ganzzahlig=False):\n",
|
|
" \"\"\"LP-Relaxation oder exaktes MILP ueber HiGHS.\"\"\"\n",
|
|
" n = len(c)\n",
|
|
" ergebnis = linprog(\n",
|
|
" c=-c, A_ub=A, b_ub=b, bounds=[(0, None)] * n,\n",
|
|
" integrality=np.ones(n) if ganzzahlig else None,\n",
|
|
" method=\"highs\")\n",
|
|
" return (-ergebnis.fun, ergebnis.x) if ergebnis.success else (None, None)\n",
|
|
"\n",
|
|
"\n",
|
|
"def abrunden_und_reparieren(x_lp, c, A, b):\n",
|
|
" \"\"\"Naive Strategie: abrunden, dann gierig auffuellen, solange zulaessig.\"\"\"\n",
|
|
" x = np.floor(x_lp + 1e-9)\n",
|
|
" verbessert = True\n",
|
|
" while verbessert: # gierig auffuellen\n",
|
|
" verbessert = False\n",
|
|
" for j in np.argsort(-c): # bester Ertrag zuerst\n",
|
|
" kandidat = x.copy()\n",
|
|
" kandidat[j] += 1\n",
|
|
" if np.all(A @ kandidat <= b + 1e-9):\n",
|
|
" x = kandidat\n",
|
|
" verbessert = True\n",
|
|
" break\n",
|
|
" return c @ x, x\n",
|
|
"\n",
|
|
"\n",
|
|
"if __name__ == \"__main__\":\n",
|
|
" rng = np.random.default_rng(2026)\n",
|
|
" print(\"=\" * 82)\n",
|
|
" print(\" WIE TEUER IST RUNDEN? (200 Zufallsinstanzen je Groesse)\")\n",
|
|
" print(\"=\" * 82)\n",
|
|
" print(f\"{'n x m':>8} | {'Aufrunden unzul.':>17} | {'Abrunden: mittl.':>17} | \"\n",
|
|
" f\"{'schlimmster':>12} | {'gierig':>8}\")\n",
|
|
" print(f\"{'':>8} | {'':>17} | {'Verlust':>17} | {'Fall':>12} | {'Verlust':>8}\")\n",
|
|
" print(\"-\" * 82)\n",
|
|
"\n",
|
|
" for n, m in [(5, 2), (10, 3), (20, 5), (40, 8)]:\n",
|
|
" unzulaessig = 0\n",
|
|
" verluste_ab, verluste_gierig = [], []\n",
|
|
"\n",
|
|
" for _ in range(200):\n",
|
|
" c, A, b = erzeuge_instanz(n, m, rng)\n",
|
|
" z_lp, x_lp = loese_lp(c, A, b, ganzzahlig=False)\n",
|
|
" z_ip, _ = loese_lp(c, A, b, ganzzahlig=True)\n",
|
|
" if z_lp is None or z_ip is None or z_ip <= 0:\n",
|
|
" continue\n",
|
|
"\n",
|
|
" # Variante 1: aufrunden\n",
|
|
" x_auf = np.ceil(x_lp - 1e-9)\n",
|
|
" if np.any(A @ x_auf > b + 1e-9):\n",
|
|
" unzulaessig += 1\n",
|
|
"\n",
|
|
" # Variante 2: abrunden\n",
|
|
" x_ab = np.floor(x_lp + 1e-9)\n",
|
|
" verluste_ab.append(1.0 - (c @ x_ab) / z_ip)\n",
|
|
"\n",
|
|
" # Variante 3: abrunden + gierig auffuellen\n",
|
|
" z_gierig, _ = abrunden_und_reparieren(x_lp, c, A, b)\n",
|
|
" verluste_gierig.append(1.0 - z_gierig / z_ip)\n",
|
|
"\n",
|
|
" print(f\"{n:>3} x {m:<3} | {unzulaessig/2:>15.1f} % | \"\n",
|
|
" f\"{np.mean(verluste_ab)*100:>15.1f} % | \"\n",
|
|
" f\"{np.max(verluste_ab)*100:>10.1f} % | \"\n",
|
|
" f\"{np.mean(verluste_gierig)*100:>6.1f} %\")\n",
|
|
"\n",
|
|
" print(\"-\" * 82)\n",
|
|
" print(\"Lesart: 'Aufrunden unzul.' = Anteil der Faelle, in denen die aufgerundete\")\n",
|
|
" print(\"Loesung eine Nebenbedingung verletzt. 'Verlust' = Abstand zum exakten Optimum.\")\n",
|
|
" print(\"=\" * 82)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Beispiel: Das Rucksackproblem\n",
|
|
"\n",
|
|
"`Rucksack.py`\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"#!/usr/bin/env python3\n",
|
|
"\n",
|
|
"# Rucksack.py\n",
|
|
"\"\"\"\n",
|
|
"Kapitel MILP: Das Rucksackproblem (Knapsack).\n",
|
|
"Zeigt LP-Relaxation, Branch-and-Bound-Ergebnis und den Preis der Ganzzahligkeit\n",
|
|
"an einem Beispiel, das man vollstaendig im Kopf nachvollziehen kann.\n",
|
|
"\"\"\"\n",
|
|
"\n",
|
|
"import numpy as np\n",
|
|
"from scipy.optimize import linprog\n",
|
|
"\n",
|
|
"GEGENSTAENDE = [\"Zelt\", \"Schlafsack\", \"Kocher\", \"Kamera\", \"Buch\", \"Wasserfilter\", \"Seil\"]\n",
|
|
"NUTZEN = np.array([40.0, 35.0, 20.0, 30.0, 8.0, 25.0, 12.0])\n",
|
|
"GEWICHT = np.array([ 6.0, 4.0, 3.0, 2.0, 1.0, 2.0, 3.0])\n",
|
|
"KAPAZITAET = 11.0 # kg\n",
|
|
"\n",
|
|
"\n",
|
|
"def loese(ganzzahlig: bool):\n",
|
|
" n = len(NUTZEN)\n",
|
|
" ergebnis = linprog(\n",
|
|
" c=-NUTZEN, A_ub=[GEWICHT], b_ub=[KAPAZITAET],\n",
|
|
" bounds=[(0, 1)] * n, # jedes Teil hoechstens einmal\n",
|
|
" integrality=np.ones(n) if ganzzahlig else None,\n",
|
|
" method=\"highs\")\n",
|
|
" return -ergebnis.fun, ergebnis.x\n",
|
|
"\n",
|
|
"\n",
|
|
"if __name__ == \"__main__\":\n",
|
|
" z_lp, x_lp = loese(ganzzahlig=False)\n",
|
|
" z_ip, x_ip = loese(ganzzahlig=True)\n",
|
|
"\n",
|
|
" print(\"=\" * 72)\n",
|
|
" print(f\" RUCKSACKPROBLEM (Kapazitaet {KAPAZITAET:.0f} kg)\")\n",
|
|
" print(\"=\" * 72)\n",
|
|
" print(f\"{'Gegenstand':<14} {'Nutzen':>7} {'kg':>5} {'Nutzen/kg':>10} \"\n",
|
|
" f\"{'LP':>7} {'MILP':>6}\")\n",
|
|
" print(\"-\" * 72)\n",
|
|
" for i, name in enumerate(GEGENSTAENDE):\n",
|
|
" # int(round(...)) statt Format \"%.0f\": vermeidet die Ausgabe \"-0\"\n",
|
|
" print(f\"{name:<14} {NUTZEN[i]:>7.0f} {GEWICHT[i]:>5.0f} \"\n",
|
|
" f\"{NUTZEN[i]/GEWICHT[i]:>10.2f} {x_lp[i]:>7.2f} {int(round(x_ip[i])):>6d}\")\n",
|
|
" print(\"-\" * 72)\n",
|
|
" print(f\"{'Gesamtnutzen':<14} {'':<7} {'':<5} {'':<10} {z_lp:>7.2f} {z_ip:>6.0f}\")\n",
|
|
" print(f\"{'Gesamtgewicht':<14} {'':<7} {'':<5} {'':<10} \"\n",
|
|
" f\"{GEWICHT @ x_lp:>7.2f} {GEWICHT @ x_ip:>6.0f}\")\n",
|
|
" print(\"-\" * 72)\n",
|
|
" print(f\"Obere Schranke aus der LP-Relaxation: {z_lp:.2f}\")\n",
|
|
" print(f\"Bestes ganzzahliges Ergebnis: {z_ip:.0f}\")\n",
|
|
" print(f\"Preis der Ganzzahligkeit: {z_lp - z_ip:.2f} \"\n",
|
|
" f\"({(1 - z_ip/z_lp)*100:.1f} %)\")\n",
|
|
"\n",
|
|
" gebrochen = [GEGENSTAENDE[i] for i in range(len(NUTZEN)) if 1e-6 < x_lp[i] < 1 - 1e-6]\n",
|
|
" print(f\"\\nIn der LP-Loesung nur teilweise eingepackt: {gebrochen}\")\n",
|
|
" print(\"Genau hier wuerde Branch-and-Bound verzweigen:\")\n",
|
|
" print(f\" Ast 1: {gebrochen[0]} bleibt ganz zuhause (x=0)\")\n",
|
|
" print(f\" Ast 2: {gebrochen[0]} kommt ganz mit (x=1)\")\n",
|
|
" print(\"=\" * 72)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Praxisfall: Portfolio mit Ordergebühren und Kardinalitätsgrenze\n",
|
|
"\n",
|
|
"`MILP_Portfolio_Fixgebuehren.py`\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"#!/usr/bin/env python3\n",
|
|
"\n",
|
|
"# MILP_Portfolio_Fixgebuehren.py\n",
|
|
"\"\"\"\n",
|
|
"Kapitel MILP: MILP-Portfolio-Selektion mit Fixkosten und Kardinalitaet.\n",
|
|
"\n",
|
|
"CP-SAT-freie, gut lesbare Formulierung ueber scipy/HiGHS (kein manueller\n",
|
|
"CSR-Matrixaufbau, deutlich leichter nachvollziehbar), mit Vergleich gegen die\n",
|
|
"Loesung OHNE Restriktionen.\n",
|
|
"\n",
|
|
"Die Auswertung laeuft ueber SolverStatus und das Loesung-Objekt aus\n",
|
|
"or_kern.py: Alle Statusfaelle werden behandelt, und der MIP-Gap steht im\n",
|
|
"Bericht - statt eines blossen 'success', das nicht verraet, ob der Solver\n",
|
|
"fertig geworden ist oder nur aufgegeben hat.\n",
|
|
"\n",
|
|
"Benoetigt: numpy, pandas, scipy, pydantic (ueber or_kern)\n",
|
|
"\"\"\"\n",
|
|
"\n",
|
|
"from __future__ import annotations\n",
|
|
"\n",
|
|
"import time\n",
|
|
"\n",
|
|
"import numpy as np\n",
|
|
"import pandas as pd\n",
|
|
"from scipy.optimize import linprog\n",
|
|
"\n",
|
|
"from or_kern import Loesung, SolverStatus, status_von_scipy\n",
|
|
"\n",
|
|
"ANLAGEN = [\"US-Aktien\", \"EU-Aktien\", \"Emerging-Markets\",\n",
|
|
" \"Staatsanleihen\", \"Unternehmensanleihen\", \"Rohstoffe\"]\n",
|
|
"RENDITE = np.array([0.11, 0.08, 0.13, 0.03, 0.05, 0.07]) # erwartet, p.a.\n",
|
|
"\n",
|
|
"BUDGET = 100_000.0\n",
|
|
"MIN_POSITION = 10_000.0 # L\n",
|
|
"MAX_POSITION = 40_000.0 # U (dient zugleich als Big-M!)\n",
|
|
"GEBUEHR = 50.0 # F, je aktivierter Position\n",
|
|
"MAX_POSITIONEN = 3 # K\n",
|
|
"\n",
|
|
"N = len(ANLAGEN)\n",
|
|
"\n",
|
|
"\n",
|
|
"def baue_und_loese(zeitlimit: float | None = None) -> Loesung:\n",
|
|
" \"\"\"\n",
|
|
" Variablenreihenfolge: [x_0..x_{N-1}, y_0..y_{N-1}]\n",
|
|
" Zielfunktion (Maximierung -> fuer linprog negiert):\n",
|
|
" max sum(rendite_i * x_i) - GEBUEHR * sum(y_i)\n",
|
|
" \"\"\"\n",
|
|
" c = np.concatenate([-RENDITE, np.full(N, GEBUEHR)]) # negiert = Minimierung\n",
|
|
"\n",
|
|
" # --- Gleichungsnebenbedingung: Budget vollstaendig investiert -----------\n",
|
|
" A_eq = np.zeros((1, 2 * N))\n",
|
|
" A_eq[0, :N] = 1.0\n",
|
|
" b_eq = np.array([BUDGET])\n",
|
|
"\n",
|
|
" zeilen, grenzen = [], []\n",
|
|
"\n",
|
|
" # --- Kardinalitaet: sum(y_i) <= K --------------------------------------\n",
|
|
" zeile = np.zeros(2 * N)\n",
|
|
" zeile[N:] = 1.0\n",
|
|
" zeilen.append(zeile)\n",
|
|
" grenzen.append(MAX_POSITIONEN)\n",
|
|
"\n",
|
|
" # --- Obergrenze (Big-M): x_i - U*y_i <= 0 -----------------------------\n",
|
|
" for i in range(N):\n",
|
|
" zeile = np.zeros(2 * N)\n",
|
|
" zeile[i] = 1.0\n",
|
|
" zeile[N + i] = -MAX_POSITION\n",
|
|
" zeilen.append(zeile)\n",
|
|
" grenzen.append(0.0)\n",
|
|
"\n",
|
|
" # --- Untergrenze: L*y_i - x_i <= 0 (entspricht x_i >= L*y_i) --------\n",
|
|
" for i in range(N):\n",
|
|
" zeile = np.zeros(2 * N)\n",
|
|
" zeile[i] = -1.0\n",
|
|
" zeile[N + i] = MIN_POSITION\n",
|
|
" zeilen.append(zeile)\n",
|
|
" grenzen.append(0.0)\n",
|
|
"\n",
|
|
" A_ub = np.array(zeilen)\n",
|
|
" b_ub = np.array(grenzen)\n",
|
|
"\n",
|
|
" schranken = [(0.0, MAX_POSITION)] * N + [(0.0, 1.0)] * N\n",
|
|
" ganzzahligkeit = np.concatenate([np.zeros(N), np.ones(N)]) # y binaer\n",
|
|
"\n",
|
|
" t0 = time.perf_counter()\n",
|
|
" ergebnis = linprog(c=c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,\n",
|
|
" bounds=schranken, integrality=ganzzahligkeit,\n",
|
|
" method=\"highs\",\n",
|
|
" options={\"time_limit\": zeitlimit} if zeitlimit else None)\n",
|
|
" laufzeit = time.perf_counter() - t0\n",
|
|
"\n",
|
|
" status = status_von_scipy(ergebnis)\n",
|
|
" if not status.brauchbar:\n",
|
|
" return Loesung(status=status, laufzeit=laufzeit)\n",
|
|
"\n",
|
|
" # Zurueck in die Maximierungswelt: Zielwert UND Schranke negieren.\n",
|
|
" # Aus beiden zusammen rechnet das Loesung-Objekt den MIP-Gap aus.\n",
|
|
" x, y = ergebnis.x[:N], np.round(ergebnis.x[N:])\n",
|
|
" return Loesung(\n",
|
|
" status=status,\n",
|
|
" werte={**{name: float(w) for name, w in zip(ANLAGEN, x)},\n",
|
|
" **{f\"aktiv:{name}\": float(w) for name, w in zip(ANLAGEN, y)}},\n",
|
|
" zielwert=float(-ergebnis.fun),\n",
|
|
" schranke=float(-ergebnis.mip_dual_bound),\n",
|
|
" laufzeit=laufzeit)\n",
|
|
"\n",
|
|
"\n",
|
|
"def ohne_restriktionen() -> tuple[float, np.ndarray]:\n",
|
|
" \"\"\"Vergleichsfall: nur Budget, keine Gebuehren/Kardinalitaet/Mindestgroesse.\"\"\"\n",
|
|
" ergebnis = linprog(c=-RENDITE, A_eq=[np.ones(N)], b_eq=[BUDGET],\n",
|
|
" bounds=[(0, None)] * N, method=\"highs\")\n",
|
|
" status = status_von_scipy(ergebnis)\n",
|
|
" if not status.brauchbar:\n",
|
|
" raise SystemExit(f\"Vergleichsfall nicht loesbar: {status.value}\")\n",
|
|
" return -ergebnis.fun, ergebnis.x\n",
|
|
"\n",
|
|
"\n",
|
|
"def pruefe_portfolio(loesung: Loesung, toleranz: float = 1e-6) -> list[str]:\n",
|
|
" \"\"\"Prueft die Loesung gegen die Anforderungen - ohne den Solver zu fragen.\n",
|
|
"\n",
|
|
" Bewusst kein assert: Eine Beanstandungsliste laesst sich protokollieren,\n",
|
|
" weiterreichen und testen. Ein assert verschwindet ausserdem, sobald\n",
|
|
" jemand Python mit -O startet.\n",
|
|
" \"\"\"\n",
|
|
" x = np.array([loesung.werte[name] for name in ANLAGEN])\n",
|
|
" y = np.array([loesung.werte[f\"aktiv:{name}\"] for name in ANLAGEN])\n",
|
|
" beanstandungen: list[str] = []\n",
|
|
"\n",
|
|
" if abs(x.sum() - BUDGET) > 1e-4:\n",
|
|
" beanstandungen.append(f\"Budget nicht exakt investiert: {x.sum():,.2f}\")\n",
|
|
" if y.sum() > MAX_POSITIONEN + toleranz:\n",
|
|
" beanstandungen.append(f\"{y.sum():.0f} Positionen statt hoechstens \"\n",
|
|
" f\"{MAX_POSITIONEN}\")\n",
|
|
" for i, name in enumerate(ANLAGEN):\n",
|
|
" if abs(y[i] - round(y[i])) > toleranz:\n",
|
|
" beanstandungen.append(f\"{name}: y = {y[i]!r} ist nicht ganzzahlig\")\n",
|
|
" elif y[i] and not (MIN_POSITION - toleranz <= x[i]\n",
|
|
" <= MAX_POSITION + toleranz):\n",
|
|
" beanstandungen.append(f\"{name}: {x[i]:,.2f} EUR verletzt die \"\n",
|
|
" f\"Groessengrenzen\")\n",
|
|
" elif not y[i] and x[i] > toleranz:\n",
|
|
" beanstandungen.append(f\"{name}: inaktiv, aber {x[i]:,.2f} EUR \"\n",
|
|
" f\"investiert\")\n",
|
|
" return beanstandungen\n",
|
|
"\n",
|
|
"\n",
|
|
"if __name__ == \"__main__\":\n",
|
|
" loesung = baue_und_loese()\n",
|
|
"\n",
|
|
" print(\"=\" * 88)\n",
|
|
" print(\" OPTIMALE MILP-PORTFOLIO-ALLOKATION MIT FIXGEBUEHREN\")\n",
|
|
" print(\"=\" * 88)\n",
|
|
" print(f\"Budget: {BUDGET:,.0f} EUR | max. {MAX_POSITIONEN} Positionen | \"\n",
|
|
" f\"je {MIN_POSITION:,.0f}-{MAX_POSITION:,.0f} EUR | Gebuehr {GEBUEHR:.0f} EUR\\n\")\n",
|
|
"\n",
|
|
" # Zuerst der Status - erst danach interessieren die Zahlen.\n",
|
|
" if loesung.status.modellfehler:\n",
|
|
" raise SystemExit(f\"Das Modell ist nicht loesbar ({loesung.status.value}). \"\n",
|
|
" f\"Naechster Schritt: Anhang Fehlerdiagnose.\")\n",
|
|
" if not loesung.status.brauchbar:\n",
|
|
" raise SystemExit(f\"Keine Loesung erhalten ({loesung.status.value}). \"\n",
|
|
" f\"Zeitlimit erhoehen oder Modell vereinfachen.\")\n",
|
|
" if loesung.status is SolverStatus.ZULAESSIG:\n",
|
|
" print(f\"ACHTUNG: nicht beweisbar optimal - Gap {loesung.gap:.2%}\\n\")\n",
|
|
"\n",
|
|
" x = np.array([loesung.werte[name] for name in ANLAGEN])\n",
|
|
" y = np.array([loesung.werte[f\"aktiv:{name}\"] for name in ANLAGEN], dtype=int)\n",
|
|
"\n",
|
|
" tabelle = pd.DataFrame({\n",
|
|
" \"Anlage\": ANLAGEN,\n",
|
|
" \"Aktiv\": [\"JA\" if y[i] else \"-\" for i in range(N)],\n",
|
|
" \"Investition (EUR)\": [f\"{x[i]:,.0f}\" for i in range(N)],\n",
|
|
" \"Anteil\": [f\"{x[i]/BUDGET*100:5.1f} %\" for i in range(N)],\n",
|
|
" \"Erw. Rendite\": [f\"{RENDITE[i]*100:4.1f} %\" for i in range(N)],\n",
|
|
" \"Erw. Ertrag (EUR)\": [f\"{x[i]*RENDITE[i]:,.0f}\" for i in range(N)],\n",
|
|
" })\n",
|
|
" print(tabelle.to_string(index=False))\n",
|
|
"\n",
|
|
" brutto = float(RENDITE @ x)\n",
|
|
" gebuehren = float(GEBUEHR * y.sum())\n",
|
|
" print(\"-\" * 88)\n",
|
|
" print(f\"Erwarteter Bruttoertrag: {brutto:>12,.2f} EUR\")\n",
|
|
" print(f\"Ordergebuehren: {-gebuehren:>12,.2f} EUR ({y.sum()} Positionen)\")\n",
|
|
" print(f\"Netto-Erwartungswert: {loesung.zielwert:>12,.2f} EUR\")\n",
|
|
" print(f\"\\n{loesung.als_bericht()}\")\n",
|
|
"\n",
|
|
" # --- Vergleich mit dem unbeschraenkten Fall ---------------------------\n",
|
|
" z_frei, x_frei = ohne_restriktionen()\n",
|
|
" print(\"-\" * 88)\n",
|
|
" print(f\"Zum Vergleich ohne jede Restriktion (alles in den Bestwert): \"\n",
|
|
" f\"{z_frei:,.2f} EUR\")\n",
|
|
" print(f\"Kosten der Realitaet (Gebuehren, Streuung, Mindestgroessen): \"\n",
|
|
" f\"{z_frei - loesung.zielwert:,.2f} EUR \"\n",
|
|
" f\"({(1 - loesung.zielwert/z_frei)*100:.2f} %)\")\n",
|
|
"\n",
|
|
" # --- Alle Nebenbedingungen nachpruefen --------------------------------\n",
|
|
" beanstandungen = pruefe_portfolio(loesung)\n",
|
|
" print(\"-\" * 88)\n",
|
|
" if beanstandungen:\n",
|
|
" raise SystemExit(\"Abnahmepruefung fehlgeschlagen:\\n - \"\n",
|
|
" + \"\\n - \".join(beanstandungen))\n",
|
|
" print(\"Abnahmepruefung: alle Nebenbedingungen geprueft und eingehalten.\")\n",
|
|
" print(\"=\" * 88)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Alle Statusfälle behandeln\n",
|
|
"\n",
|
|
"`Solverstatus_und_Gap.py`\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"#!/usr/bin/env python3\n",
|
|
"\n",
|
|
"# Solverstatus_und_Gap.py\n",
|
|
"\"\"\"\n",
|
|
"Kapitel MILP: Was tun, wenn der Solver nicht fertig wird?\n",
|
|
"\n",
|
|
"Bei einem LP kommt entweder eine optimale Loesung oder eine klare Absage.\n",
|
|
"Bei einem MILP ist der haeufigste Ausgang im Betrieb ein dritter: \"Ich habe\n",
|
|
"eine Loesung, ich weiss aber nicht, ob sie die beste ist - und die Zeit ist\n",
|
|
"um.\" Dieses Programm zeigt, wie man mit diesem Fall umgeht.\n",
|
|
"\n",
|
|
" 1. Alle Statusfaelle explizit behandeln, statt OPTIMAL vorauszusetzen.\n",
|
|
" 2. Den MIP-Gap lesen: Wie weit kann ich hoechstens danebenliegen?\n",
|
|
" 3. Messen, was zusaetzliche Rechenzeit ueberhaupt noch bringt.\n",
|
|
" 4. Warm-Start ausprobieren - und ehrlich messen, ob er etwas bringt.\n",
|
|
"\n",
|
|
"Beispiel: Standortplanung, 45 moegliche Lager, 120 Kunden (5445 Variablen,\n",
|
|
"davon 45 binaer).\n",
|
|
"\n",
|
|
"Benoetigt: numpy, highspy\n",
|
|
"\"\"\"\n",
|
|
"\n",
|
|
"from __future__ import annotations\n",
|
|
"\n",
|
|
"import time\n",
|
|
"from dataclasses import dataclass\n",
|
|
"\n",
|
|
"import numpy as np\n",
|
|
"import highspy\n",
|
|
"\n",
|
|
"RNG = np.random.default_rng(7)\n",
|
|
"\n",
|
|
"N_LAGER, N_KUNDE = 45, 120\n",
|
|
"FIXKOSTEN = RNG.uniform(3000, 9000, N_LAGER)\n",
|
|
"TRANSPORT = RNG.uniform(5, 60, (N_LAGER, N_KUNDE))\n",
|
|
"BEDARF = RNG.uniform(10, 60, N_KUNDE)\n",
|
|
"KAPAZITAET = np.full(N_LAGER, BEDARF.sum() * 0.22)\n",
|
|
"\n",
|
|
"\n",
|
|
"@dataclass\n",
|
|
"class Ergebnis:\n",
|
|
" \"\"\"Alles, was nach einem Solverlauf ausgewertet werden muss - nicht nur\n",
|
|
" der Zielwert.\"\"\"\n",
|
|
" status: str\n",
|
|
" brauchbar: bool # Gibt es ueberhaupt eine zulaessige Loesung?\n",
|
|
" beweisbar_optimal: bool\n",
|
|
" zielwert: float # bester gefundener Wert (Incumbent)\n",
|
|
" schranke: float # beste bewiesene Schranke (Dual Bound)\n",
|
|
" gap: float # relativer Abstand zwischen beiden\n",
|
|
" knoten: int\n",
|
|
" dauer: float\n",
|
|
" loesung: np.ndarray\n",
|
|
"\n",
|
|
"\n",
|
|
"def loese(zeitlimit: float, startloesung: np.ndarray | None = None) -> Ergebnis:\n",
|
|
" \"\"\"Loest das Standortmodell mit Zeitlimit und wertet ALLE Statusfaelle aus.\"\"\"\n",
|
|
" hochschule = highspy.Highs()\n",
|
|
" hochschule.setOptionValue(\"output_flag\", False)\n",
|
|
" hochschule.setOptionValue(\"time_limit\", zeitlimit)\n",
|
|
"\n",
|
|
" anzahl_x = N_LAGER * N_KUNDE\n",
|
|
" unendlich = highspy.kHighsInf\n",
|
|
"\n",
|
|
" hochschule.addVars(anzahl_x, np.zeros(anzahl_x), np.full(anzahl_x, unendlich))\n",
|
|
" hochschule.addVars(N_LAGER, np.zeros(N_LAGER), np.ones(N_LAGER))\n",
|
|
" for i in range(N_LAGER):\n",
|
|
" hochschule.changeColIntegrality(anzahl_x + i, highspy.HighsVarType.kInteger)\n",
|
|
" hochschule.changeColCost(anzahl_x + i, FIXKOSTEN[i])\n",
|
|
" for j in range(N_KUNDE):\n",
|
|
" hochschule.changeColCost(i * N_KUNDE + j, TRANSPORT[i, j])\n",
|
|
"\n",
|
|
" for j in range(N_KUNDE):\n",
|
|
" index = np.array([i * N_KUNDE + j for i in range(N_LAGER)], dtype=np.int32)\n",
|
|
" hochschule.addRow(BEDARF[j], BEDARF[j], len(index), index, np.ones(len(index)))\n",
|
|
"\n",
|
|
" for i in range(N_LAGER):\n",
|
|
" index = np.array([i * N_KUNDE + j for j in range(N_KUNDE)] + [anzahl_x + i],\n",
|
|
" dtype=np.int32)\n",
|
|
" werte = np.concatenate([np.ones(N_KUNDE), [-KAPAZITAET[i]]])\n",
|
|
" hochschule.addRow(-unendlich, 0.0, len(index), index, werte)\n",
|
|
"\n",
|
|
" if startloesung is not None:\n",
|
|
" hochschule.setSolution(len(startloesung),\n",
|
|
" np.arange(len(startloesung), dtype=np.int32),\n",
|
|
" startloesung)\n",
|
|
"\n",
|
|
" t0 = time.perf_counter()\n",
|
|
" hochschule.run()\n",
|
|
" dauer = time.perf_counter() - t0\n",
|
|
"\n",
|
|
" status = hochschule.modelStatusToString(hochschule.getModelStatus())\n",
|
|
" info = hochschule.getInfo()\n",
|
|
"\n",
|
|
" # Der Kern der Sache: Aus dem Status folgt, WAS man mit dem Ergebnis\n",
|
|
" # ueberhaupt anfangen darf.\n",
|
|
" brauchbar = status in (\"Optimal\", \"Time limit reached\", \"Solution limit reached\")\n",
|
|
" beweisbar_optimal = status == \"Optimal\"\n",
|
|
" if status in (\"Infeasible\", \"Unbounded\", \"Primal infeasible or unbounded\"):\n",
|
|
" brauchbar = False\n",
|
|
"\n",
|
|
" return Ergebnis(\n",
|
|
" status=status,\n",
|
|
" brauchbar=brauchbar and info.objective_function_value < unendlich,\n",
|
|
" beweisbar_optimal=beweisbar_optimal,\n",
|
|
" zielwert=info.objective_function_value,\n",
|
|
" schranke=info.mip_dual_bound,\n",
|
|
" gap=info.mip_gap,\n",
|
|
" knoten=info.mip_node_count,\n",
|
|
" dauer=dauer,\n",
|
|
" loesung=np.array(hochschule.getSolution().col_value),\n",
|
|
" )\n",
|
|
"\n",
|
|
"\n",
|
|
"def gieriger_startplan() -> tuple[np.ndarray, float]:\n",
|
|
" \"\"\"Eine Faustregel-Loesung, wie sie ein Disponent von Hand erstellen wuerde:\n",
|
|
" die guenstigsten Lager oeffnen (Fixkosten je Kapazitaetseinheit), dann\n",
|
|
" jeden Kunden dem naechstgelegenen offenen Lager mit Restkapazitaet\n",
|
|
" zuordnen. Kein Solver noetig - und in Sekunden fertig.\n",
|
|
" \"\"\"\n",
|
|
" reihenfolge = np.argsort(FIXKOSTEN / KAPAZITAET)\n",
|
|
" offen: list[int] = []\n",
|
|
" for i in reihenfolge:\n",
|
|
" offen.append(int(i))\n",
|
|
" if KAPAZITAET[offen].sum() >= BEDARF.sum() * 1.05:\n",
|
|
" break\n",
|
|
"\n",
|
|
" rest = KAPAZITAET.copy()\n",
|
|
" x = np.zeros((N_LAGER, N_KUNDE))\n",
|
|
" for j in np.argsort(-BEDARF): # groesste Kunden zuerst\n",
|
|
" for i in sorted(offen, key=lambda i: TRANSPORT[i, j]):\n",
|
|
" menge = min(rest[i], BEDARF[j] - x[:, j].sum())\n",
|
|
" if menge > 1e-9:\n",
|
|
" x[i, j] += menge\n",
|
|
" rest[i] -= menge\n",
|
|
" if abs(x[:, j].sum() - BEDARF[j]) < 1e-9:\n",
|
|
" break\n",
|
|
"\n",
|
|
" y = np.zeros(N_LAGER)\n",
|
|
" y[offen] = 1.0\n",
|
|
" kosten = float((FIXKOSTEN * y).sum() + (TRANSPORT * x).sum())\n",
|
|
" return np.concatenate([x.ravel(), y]), kosten\n",
|
|
"\n",
|
|
"\n",
|
|
"def zeige(titel: str, e: Ergebnis) -> None:\n",
|
|
" print(f\"\\n{titel}\")\n",
|
|
" print(f\" Status {e.status}\")\n",
|
|
" if not e.brauchbar:\n",
|
|
" print(\" -> KEINE verwertbare Loesung. Nicht weiterrechnen!\")\n",
|
|
" return\n",
|
|
" print(f\" bester Plan (Incumbent) {e.zielwert:>12,.2f} EUR\")\n",
|
|
" print(f\" bewiesene Schranke {e.schranke:>12,.2f} EUR\")\n",
|
|
" print(f\" MIP-Gap {e.gap * 100:>12.3f} %\")\n",
|
|
" print(f\" Knoten / Zeit {e.knoten:>12,} / {e.dauer:.2f} s\")\n",
|
|
" if e.beweisbar_optimal:\n",
|
|
" print(\" -> beweisbar optimal\")\n",
|
|
" else:\n",
|
|
" print(f\" -> zulaessig, aber nicht bewiesen optimal. Der wahre Bestwert\")\n",
|
|
" print(f\" liegt zwischen {e.schranke:,.2f} und {e.zielwert:,.2f} EUR.\")\n",
|
|
"\n",
|
|
"\n",
|
|
"if __name__ == \"__main__\":\n",
|
|
" print(\"=\" * 78)\n",
|
|
" print(\" MIP-GAP UND ZEITLIMIT: STANDORTPLANUNG, 45 LAGER, 120 KUNDEN\")\n",
|
|
" print(\"=\" * 78)\n",
|
|
" print(f\"{N_LAGER * N_KUNDE + N_LAGER:,} Variablen, davon {N_LAGER} binaer.\")\n",
|
|
"\n",
|
|
" kurz = loese(2.0)\n",
|
|
" zeige(\"[1] Zeitlimit 2 Sekunden\", kurz)\n",
|
|
"\n",
|
|
" lang = loese(60.0)\n",
|
|
" zeige(\"[2] Zeitlimit 60 Sekunden\", lang)\n",
|
|
"\n",
|
|
" print(\"\\n\" + \"=\" * 78)\n",
|
|
" print(\" WAS BRINGT MEHR RECHENZEIT?\")\n",
|
|
" print(\"=\" * 78)\n",
|
|
" verbesserung = kurz.zielwert - lang.zielwert\n",
|
|
" print(f\"Nach 2 Sekunden: {kurz.zielwert:,.2f} EUR, Gap {kurz.gap * 100:.2f} %\")\n",
|
|
" print(f\"Nach {lang.dauer:.1f} Sekunden: {lang.zielwert:,.2f} EUR, bewiesen optimal\")\n",
|
|
" print(f\"Gewinn durch {lang.dauer - kurz.dauer:.1f} Sekunden mehr Rechenzeit: \"\n",
|
|
" f\"{verbesserung:,.2f} EUR \"\n",
|
|
" f\"({verbesserung / kurz.zielwert * 100:.2f} %)\")\n",
|
|
" print()\n",
|
|
" print(\"Das ist die Frage, die im Betrieb wirklich zaehlt: Der Gap von\")\n",
|
|
" print(f\"{kurz.gap * 100:.1f} % nach 2 Sekunden ist eine GARANTIE - schlechter als\")\n",
|
|
" print(\"dieser Wert kann die Loesung nicht sein. Ob sich die restliche\")\n",
|
|
" print(\"Rechenzeit lohnt, entscheidet nicht der Solver, sondern die Anwendung:\")\n",
|
|
" print(\"Ein naechtlicher Tourenplan darf eine Stunde rechnen, eine Umplanung\")\n",
|
|
" print(\"bei Maschinenausfall hat 30 Sekunden.\")\n",
|
|
"\n",
|
|
" print(\"\\n\" + \"=\" * 78)\n",
|
|
" print(\" BRINGT EIN WARM-START ETWAS?\")\n",
|
|
" print(\"=\" * 78)\n",
|
|
" start, start_kosten = gieriger_startplan()\n",
|
|
" print(f\"Faustregel-Startplan (ohne Solver): {start_kosten:,.2f} EUR\")\n",
|
|
" print(f\"Das sind {(start_kosten / lang.zielwert - 1) * 100:.1f} % ueber dem Optimum.\\n\")\n",
|
|
"\n",
|
|
" warm = loese(60.0, startloesung=start)\n",
|
|
" print(f\"{'':26} {'Zeit':>9} {'Knoten':>9} {'Ziel':>13}\")\n",
|
|
" print(\"-\" * 78)\n",
|
|
" print(f\"{'ohne Warm-Start':<26} {lang.dauer:>8.2f}s {lang.knoten:>9,} \"\n",
|
|
" f\"{lang.zielwert:>13,.2f}\")\n",
|
|
" print(f\"{'mit Warm-Start':<26} {warm.dauer:>8.2f}s {warm.knoten:>9,} \"\n",
|
|
" f\"{warm.zielwert:>13,.2f}\")\n",
|
|
" print(\"-\" * 78)\n",
|
|
" print(\"Ergebnis: praktisch kein Unterschied. Der Grund ist nicht, dass\")\n",
|
|
" print(\"Warm-Starts nichts taugen - sondern dass HiGHS' eigene Heuristiken\")\n",
|
|
" print(\"innerhalb der ersten Sekunde bereits eine BESSERE Loesung finden als\")\n",
|
|
" print(\"unsere Faustregel. Ein Startwert hilft nur, wenn er besser ist als\")\n",
|
|
" print(\"das, was der Solver von allein in derselben Zeit findet.\")\n",
|
|
" print()\n",
|
|
" print(\"Warm-Starts lohnen sich damit vor allem in zwei Faellen:\")\n",
|
|
" print(\" * Sie haben Domaenenwissen, das der Solver nicht hat (siehe\")\n",
|
|
" print(\" Warmstart_Effekt.py - dort halbiert ein Heuristik-Hinweis die Zeit).\")\n",
|
|
" print(\" * Sie planen laufend neu und der gestrige Plan ist fast noch gueltig.\")\n",
|
|
" print(\"In beiden Faellen gilt: MESSEN, nicht glauben.\")\n",
|
|
" print(\"=\" * 78)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Warm-Starts: was sie können und was nicht\n",
|
|
"\n",
|
|
"`Warmstart_Effekt.py`\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"#!/usr/bin/env python3\n",
|
|
"\n",
|
|
"# Warmstart_Effekt.py\n",
|
|
"\"\"\"\n",
|
|
"Kapitel MILP: Wann ein Warm-Start wirklich etwas bringt.\n",
|
|
"\n",
|
|
"Solverstatus_und_Gap.py zeigt einen Fall, in dem ein Startwert NICHTS\n",
|
|
"bringt - HiGHS findet von allein schneller etwas Besseres. Hier der\n",
|
|
"Gegenfall: eine Heuristik, die dem Solver echtes Domaenenwissen liefert.\n",
|
|
"\n",
|
|
"Problem: Lastverteilung. n Auftraege mit bekannter Dauer sind auf m\n",
|
|
"gleichartige Maschinen zu verteilen, sodass die zuletzt fertige Maschine\n",
|
|
"so frueh wie moeglich fertig wird (Makespan-Minimierung).\n",
|
|
"\n",
|
|
"Die Heuristik: LPT (Longest Processing Time first) - laengste Auftraege\n",
|
|
"zuerst, jeder auf die momentan am wenigsten belastete Maschine. Sie ist\n",
|
|
"Jahrzehnte alt, in zwei Zeilen geschrieben und beweisbar nie schlechter\n",
|
|
"als 4/3 des Optimums.\n",
|
|
"\n",
|
|
"Ueber model.AddHint() bekommt CP-SAT diese Loesung als Startpunkt.\n",
|
|
"\n",
|
|
"WICHTIG: highspy wird hier bewusst NICHT importiert - es vertraegt sich\n",
|
|
"nicht mit ortools im selben Prozess (siehe Kapitel Oekosystem).\n",
|
|
"\n",
|
|
"Benoetigt: numpy, ortools\n",
|
|
"\"\"\"\n",
|
|
"\n",
|
|
"from __future__ import annotations\n",
|
|
"\n",
|
|
"import time\n",
|
|
"\n",
|
|
"import numpy as np\n",
|
|
"from ortools.sat.python import cp_model\n",
|
|
"\n",
|
|
"RNG = np.random.default_rng(4)\n",
|
|
"\n",
|
|
"\n",
|
|
"def lpt_heuristik(dauer: np.ndarray, n_maschinen: int) -> tuple[np.ndarray, int]:\n",
|
|
" \"\"\"Longest Processing Time first.\n",
|
|
"\n",
|
|
" Laengste Auftraege zuerst auf die jeweils freieste Maschine legen. Zwei\n",
|
|
" Zeilen, keine Bibliothek, Ergebnis in Mikrosekunden - und erstaunlich\n",
|
|
" nah am Optimum.\n",
|
|
" \"\"\"\n",
|
|
" zuordnung = np.zeros(len(dauer), dtype=int)\n",
|
|
" belegung = np.zeros(n_maschinen)\n",
|
|
" for auftrag in np.argsort(-dauer): # laengster zuerst\n",
|
|
" maschine = int(np.argmin(belegung)) # freieste Maschine\n",
|
|
" zuordnung[auftrag] = maschine\n",
|
|
" belegung[maschine] += dauer[auftrag]\n",
|
|
" return zuordnung, int(belegung.max())\n",
|
|
"\n",
|
|
"\n",
|
|
"def loese(dauer: np.ndarray, n_maschinen: int, zeitlimit: float,\n",
|
|
" hinweis: np.ndarray | None = None) -> tuple[str, int, float]:\n",
|
|
" \"\"\"Exaktes Modell mit CP-SAT, optional mit Startloesung als Hinweis.\"\"\"\n",
|
|
" n_auftraege = len(dauer)\n",
|
|
" obergrenze = int(dauer.sum())\n",
|
|
"\n",
|
|
" modell = cp_model.CpModel()\n",
|
|
" # x[i][k] = 1 <=> Auftrag i laeuft auf Maschine k\n",
|
|
" x = [[modell.NewBoolVar(f\"x_{i}_{k}\") for k in range(n_maschinen)]\n",
|
|
" for i in range(n_auftraege)]\n",
|
|
" for i in range(n_auftraege):\n",
|
|
" modell.AddExactlyOne(x[i]) # jeder Auftrag genau einmal\n",
|
|
"\n",
|
|
" belegung = [modell.NewIntVar(0, obergrenze, f\"last_{k}\")\n",
|
|
" for k in range(n_maschinen)]\n",
|
|
" for k in range(n_maschinen):\n",
|
|
" modell.Add(belegung[k] == sum(int(dauer[i]) * x[i][k]\n",
|
|
" for i in range(n_auftraege)))\n",
|
|
"\n",
|
|
" makespan = modell.NewIntVar(0, obergrenze, \"makespan\")\n",
|
|
" modell.AddMaxEquality(makespan, belegung) # das Maximum ueber alle Maschinen\n",
|
|
" modell.Minimize(makespan)\n",
|
|
"\n",
|
|
" # Der Warm-Start: ein Hinweis pro Variable. CP-SAT muss ihn nicht\n",
|
|
" # befolgen - er nutzt ihn als erste Loesung, wenn er zulaessig ist.\n",
|
|
" if hinweis is not None:\n",
|
|
" for i in range(n_auftraege):\n",
|
|
" for k in range(n_maschinen):\n",
|
|
" modell.AddHint(x[i][k], 1 if hinweis[i] == k else 0)\n",
|
|
"\n",
|
|
" loeser = cp_model.CpSolver()\n",
|
|
" loeser.parameters.max_time_in_seconds = zeitlimit\n",
|
|
" # Ein Arbeiter und fester Startwert, damit die Messung reproduzierbar ist.\n",
|
|
" # Im Produktivbetrieb laesst man beides auf den Standardwerten.\n",
|
|
" loeser.parameters.num_workers = 1\n",
|
|
" loeser.parameters.random_seed = 1\n",
|
|
"\n",
|
|
" t0 = time.perf_counter()\n",
|
|
" status = loeser.Solve(modell)\n",
|
|
" dauer_s = time.perf_counter() - t0\n",
|
|
"\n",
|
|
" if status not in (cp_model.OPTIMAL, cp_model.FEASIBLE):\n",
|
|
" raise RuntimeError(f\"Kein Plan gefunden: {loeser.StatusName(status)}\")\n",
|
|
" return loeser.StatusName(status), int(loeser.ObjectiveValue()), dauer_s\n",
|
|
"\n",
|
|
"\n",
|
|
"if __name__ == \"__main__\":\n",
|
|
" print(\"=\" * 78)\n",
|
|
" print(\" WARM-START: WENN DIE HEURISTIK MEHR WEISS ALS DER SOLVER\")\n",
|
|
" print(\"=\" * 78)\n",
|
|
" print(\"Lastverteilung: Auftraege auf gleichartige Maschinen verteilen,\")\n",
|
|
" print(\"sodass die letzte Maschine so frueh wie moeglich fertig wird.\\n\")\n",
|
|
"\n",
|
|
" print(f\"{'Instanz':<22} {'Variante':<22} {'Makespan':>9} {'Zeit':>9} \"\n",
|
|
" f\"{'Faktor':>8}\")\n",
|
|
" print(\"-\" * 78)\n",
|
|
"\n",
|
|
" for n_auftraege, n_maschinen in [(60, 7), (80, 9)]:\n",
|
|
" dauer = RNG.integers(10, 90, n_auftraege)\n",
|
|
" start, lpt_wert = lpt_heuristik(dauer, n_maschinen)\n",
|
|
" untere_schranke = dauer.sum() / n_maschinen\n",
|
|
"\n",
|
|
" instanz = f\"{n_auftraege} Auftr., {n_maschinen} Masch.\"\n",
|
|
" print(f\"{instanz:<22} {'LPT-Heuristik':<22} {lpt_wert:>9} \"\n",
|
|
" f\"{'< 0.001s':>9} {'':>8}\")\n",
|
|
"\n",
|
|
" _, ziel_kalt, zeit_kalt = loese(dauer, n_maschinen, 60.0)\n",
|
|
" print(f\"{'':<22} {'CP-SAT kalt':<22} {ziel_kalt:>9} \"\n",
|
|
" f\"{zeit_kalt:>8.3f}s {'1,0x':>8}\")\n",
|
|
"\n",
|
|
" _, ziel_warm, zeit_warm = loese(dauer, n_maschinen, 60.0, hinweis=start)\n",
|
|
" print(f\"{'':<22} {'CP-SAT + LPT-Hinweis':<22} {ziel_warm:>9} \"\n",
|
|
" f\"{zeit_warm:>8.3f}s {zeit_kalt / zeit_warm:>7.1f}x\")\n",
|
|
"\n",
|
|
" assert ziel_kalt == ziel_warm, \\\n",
|
|
" \"Der Hinweis darf das Optimum nicht veraendern - nur den Weg dorthin!\"\n",
|
|
" print(f\"{'':<22} {'untere Schranke':<22} {untere_schranke:>9.1f}\")\n",
|
|
" print(\"-\" * 78)\n",
|
|
"\n",
|
|
" print(\"\\nDrei Beobachtungen:\")\n",
|
|
" print(\"1. Der Hinweis aendert das ERGEBNIS nicht - beide Laeufe finden\")\n",
|
|
" print(\" dasselbe Optimum. Er aendert nur, wie lange der Beweis dauert.\")\n",
|
|
" print(\" Genau deshalb ist ein Warm-Start ungefaehrlich: Ein schlechter\")\n",
|
|
" print(\" Hinweis kostet Zeit, er verfaelscht aber nie die Loesung.\")\n",
|
|
" print(\"2. Die LPT-Heuristik liegt schon sehr nah am Optimum. Ihr Wert fuer\")\n",
|
|
" print(\" den Solver liegt weniger in der Qualitaet als darin, dass sie\")\n",
|
|
" print(\" SOFORT da ist - der Solver kann von Beginn an alles verwerfen,\")\n",
|
|
" print(\" was schlechter ist.\")\n",
|
|
" print(\"3. Der Faktor schwankt von Instanz zu Instanz - oben 2,2x und 1,4x -\")\n",
|
|
" print(\" und laesst sich NICHT aus der Problemgroesse ableiten. Er haengt\")\n",
|
|
" print(\" davon ab, wie schnell der Solver von allein eine vergleichbar gute\")\n",
|
|
" print(\" Loesung findet. Das ist die eigentliche Lehre: Ein Warm-Start ist\")\n",
|
|
" print(\" eine Messung wert, keine Glaubensfrage.\")\n",
|
|
" print(\"=\" * 78)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Finde den Denkfehler\n",
|
|
"\n",
|
|
"`Big_M_Falle.py`\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"#!/usr/bin/env python3\n",
|
|
"\n",
|
|
"# Big_M_Falle.py\n",
|
|
"\"\"\"\n",
|
|
"Kapitel MILP: Was ein zu grosses Big-M wirklich anrichtet.\n",
|
|
"\n",
|
|
"Lehrbuecher warnen vor grossem M mit dem Hinweis \"die Relaxation wird\n",
|
|
"schwach, der Solver braucht mehr Knoten\". Das stimmt - ist aber die\n",
|
|
"harmlosere Haelfte der Wahrheit. Die gefaehrlichere: Bei sehr grossem M\n",
|
|
"kann der Solver eine Loesung als ganzzahlig ANNEHMEN, in der die\n",
|
|
"Binaervariablen bei 1e-8 stehen. Dann liefern zugeschaltete Anlagen Ware\n",
|
|
"aus, waehrend das Modell ihre Fixkosten mit 0 verbucht.\n",
|
|
"\n",
|
|
"Beispiel: Standortplanung, 12 moegliche Lager, 40 Kunden.\n",
|
|
" min sum_i fix_i * y_i + sum_ij kosten_ij * x_ij\n",
|
|
" u.d.N. sum_i x_ij = bedarf_j (jeder Kunde beliefert)\n",
|
|
" sum_j x_ij <= M_i * y_i (Lager offen, wenn es liefert)\n",
|
|
" y_i binaer, x_ij >= 0\n",
|
|
"\n",
|
|
"Zwei Laeufe:\n",
|
|
" 1. M knapp gewaehlt (= tatsaechliche Lagerkapazitaet) -> richtig\n",
|
|
" 2. M = 1e7-fach zu gross, Presolve abgeschaltet -> stilles Desaster\n",
|
|
"Beide werden mit derselben Pruefung kontrolliert, die den Fall auffliegen laesst.\n",
|
|
"\n",
|
|
"Benoetigt: numpy, highspy\n",
|
|
"\"\"\"\n",
|
|
"\n",
|
|
"from __future__ import annotations\n",
|
|
"\n",
|
|
"import time\n",
|
|
"\n",
|
|
"import numpy as np\n",
|
|
"import highspy\n",
|
|
"\n",
|
|
"RNG = np.random.default_rng(11)\n",
|
|
"\n",
|
|
"N_LAGER, N_KUNDE = 12, 40\n",
|
|
"FIXKOSTEN = RNG.uniform(2000, 5000, N_LAGER)\n",
|
|
"TRANSPORT = RNG.uniform(5, 40, (N_LAGER, N_KUNDE))\n",
|
|
"BEDARF = RNG.uniform(10, 60, N_KUNDE)\n",
|
|
"KAPAZITAET = BEDARF.sum() * 0.45 # jedes Lager schafft 45 % des Gesamtbedarfs\n",
|
|
"\n",
|
|
"# Toleranz, ab der ein Solver eine Variable als ganzzahlig durchgehen laesst.\n",
|
|
"# HiGHS und die meisten anderen verwenden 1e-6 als Standard.\n",
|
|
"GANZZAHL_TOLERANZ = 1e-6\n",
|
|
"\n",
|
|
"\n",
|
|
"def loese(big_m: float, presolve: str = \"on\") -> dict:\n",
|
|
" \"\"\"Baut und loest das Standortmodell. Liefert Loesung und Solverkennzahlen.\"\"\"\n",
|
|
" hochschule = highspy.Highs()\n",
|
|
" hochschule.setOptionValue(\"output_flag\", False)\n",
|
|
" hochschule.setOptionValue(\"presolve\", presolve)\n",
|
|
" hochschule.setOptionValue(\"time_limit\", 300.0)\n",
|
|
"\n",
|
|
" anzahl_x = N_LAGER * N_KUNDE\n",
|
|
" unendlich = highspy.kHighsInf\n",
|
|
"\n",
|
|
" # Spalten: erst alle x_ij, dann die y_i\n",
|
|
" hochschule.addVars(anzahl_x, np.zeros(anzahl_x), np.full(anzahl_x, unendlich))\n",
|
|
" hochschule.addVars(N_LAGER, np.zeros(N_LAGER), np.ones(N_LAGER))\n",
|
|
" for i in range(N_LAGER):\n",
|
|
" hochschule.changeColIntegrality(anzahl_x + i,\n",
|
|
" highspy.HighsVarType.kInteger)\n",
|
|
" hochschule.changeColCost(anzahl_x + i, FIXKOSTEN[i])\n",
|
|
" for j in range(N_KUNDE):\n",
|
|
" hochschule.changeColCost(i * N_KUNDE + j, TRANSPORT[i, j])\n",
|
|
"\n",
|
|
" # Jeder Kunde wird genau beliefert\n",
|
|
" for j in range(N_KUNDE):\n",
|
|
" index = np.array([i * N_KUNDE + j for i in range(N_LAGER)], dtype=np.int32)\n",
|
|
" hochschule.addRow(BEDARF[j], BEDARF[j], len(index), index,\n",
|
|
" np.ones(len(index)))\n",
|
|
"\n",
|
|
" # Die Kopplung: sum_j x_ij - M * y_i <= 0\n",
|
|
" for i in range(N_LAGER):\n",
|
|
" index = np.array([i * N_KUNDE + j for j in range(N_KUNDE)] + [anzahl_x + i],\n",
|
|
" dtype=np.int32)\n",
|
|
" werte = np.concatenate([np.ones(N_KUNDE), [-big_m]])\n",
|
|
" hochschule.addRow(-unendlich, 0.0, len(index), index, werte)\n",
|
|
"\n",
|
|
" t0 = time.perf_counter()\n",
|
|
" hochschule.run()\n",
|
|
" dauer = time.perf_counter() - t0\n",
|
|
"\n",
|
|
" loesung = np.array(hochschule.getSolution().col_value)\n",
|
|
" info = hochschule.getInfo()\n",
|
|
" return {\n",
|
|
" \"x\": loesung[:anzahl_x].reshape(N_LAGER, N_KUNDE),\n",
|
|
" \"y\": loesung[anzahl_x:],\n",
|
|
" \"zielwert\": info.objective_function_value,\n",
|
|
" \"knoten\": info.mip_node_count,\n",
|
|
" \"dauer\": dauer,\n",
|
|
" }\n",
|
|
"\n",
|
|
"\n",
|
|
"def pruefe(ergebnis: dict) -> tuple[bool, list[str]]:\n",
|
|
" \"\"\"Die Pruefung, die in jedes MILP-Auswertungsskript gehoert.\n",
|
|
"\n",
|
|
" Sie rechnet die Kosten AUS DER LOESUNG neu aus, statt dem Zielwert des\n",
|
|
" Solvers zu glauben - und vergleicht beide. Genau diese Gegenrechnung\n",
|
|
" entlarvt eine Loesung, in der Binaervariablen bei 1e-8 haengengeblieben\n",
|
|
" sind.\n",
|
|
" \"\"\"\n",
|
|
" beanstandungen = []\n",
|
|
" x, y = ergebnis[\"x\"], ergebnis[\"y\"]\n",
|
|
"\n",
|
|
" # 1. Sind die Binaervariablen wirklich binaer?\n",
|
|
" abstand = np.abs(y - np.round(y))\n",
|
|
" if abstand.max() > GANZZAHL_TOLERANZ:\n",
|
|
" beanstandungen.append(\n",
|
|
" f\"y ist nicht ganzzahlig: groesster Abstand {abstand.max():.2e}\")\n",
|
|
"\n",
|
|
" # 2. Liefert ein Lager, dessen Schalter aus ist?\n",
|
|
" liefert = x.sum(axis=1) > 1e-6\n",
|
|
" geschlossen_aber_aktiv = np.where(liefert & (y < 0.5))[0]\n",
|
|
" if geschlossen_aber_aktiv.size:\n",
|
|
" beanstandungen.append(\n",
|
|
" f\"Lager {geschlossen_aber_aktiv.tolist()} liefern Ware, \"\n",
|
|
" f\"gelten im Modell aber als geschlossen\")\n",
|
|
"\n",
|
|
" # 3. Stimmt der Zielwert mit den echten Kosten ueberein?\n",
|
|
" echte_fixkosten = FIXKOSTEN[liefert].sum()\n",
|
|
" echte_transportkosten = float((TRANSPORT * x).sum())\n",
|
|
" echte_kosten = echte_fixkosten + echte_transportkosten\n",
|
|
" if abs(echte_kosten - ergebnis[\"zielwert\"]) > 1e-4 * max(1.0, echte_kosten):\n",
|
|
" beanstandungen.append(\n",
|
|
" f\"Zielwert {ergebnis['zielwert']:,.2f} weicht von den echten Kosten \"\n",
|
|
" f\"{echte_kosten:,.2f} ab (Differenz {echte_kosten - ergebnis['zielwert']:,.2f})\")\n",
|
|
"\n",
|
|
" return not beanstandungen, beanstandungen\n",
|
|
"\n",
|
|
"\n",
|
|
"def zeige(titel: str, ergebnis: dict) -> None:\n",
|
|
" x, y = ergebnis[\"x\"], ergebnis[\"y\"]\n",
|
|
" liefert = x.sum(axis=1) > 1e-6\n",
|
|
" print(f\"\\n{titel}\")\n",
|
|
" print(f\" Zielwert laut Solver {ergebnis['zielwert']:>14,.2f} EUR\")\n",
|
|
" print(f\" Knoten / Zeit {ergebnis['knoten']:>14,} \"\n",
|
|
" f\"/ {ergebnis['dauer']:.3f} s\")\n",
|
|
" print(f\" Lager mit y = 1 {int((y > 0.5).sum()):>14}\")\n",
|
|
" print(f\" Lager, die tatsaechlich liefern {int(liefert.sum()):>10}\")\n",
|
|
" unter = y[y < 0.5]\n",
|
|
" print(f\" groesster y-Wert unter 0.5 \"\n",
|
|
" f\"{(unter.max() if unter.size else 0.0):>14.3e}\")\n",
|
|
" print(f\" Fixkosten real / verbucht {FIXKOSTEN[liefert].sum():>14,.2f} \"\n",
|
|
" f\"/ {float((FIXKOSTEN * y).sum()):,.2f} EUR\")\n",
|
|
"\n",
|
|
" in_ordnung, beanstandungen = pruefe(ergebnis)\n",
|
|
" if in_ordnung:\n",
|
|
" print(\" PRUEFUNG: bestanden\")\n",
|
|
" else:\n",
|
|
" print(\" PRUEFUNG: DURCHGEFALLEN\")\n",
|
|
" for text in beanstandungen:\n",
|
|
" print(f\" - {text}\")\n",
|
|
"\n",
|
|
"\n",
|
|
"if __name__ == \"__main__\":\n",
|
|
" print(\"=\" * 78)\n",
|
|
" print(\" DIE BIG-M-FALLE: STANDORTPLANUNG, 12 LAGER, 40 KUNDEN\")\n",
|
|
" print(\"=\" * 78)\n",
|
|
" print(f\"Tatsaechliche Lagerkapazitaet: {KAPAZITAET:,.1f} Einheiten.\")\n",
|
|
" print(\"Genau das ist das kleinstmoegliche gueltige M - mehr kann ein Lager\")\n",
|
|
" print(\"ohnehin nicht ausliefern.\")\n",
|
|
"\n",
|
|
" knapp = loese(KAPAZITAET)\n",
|
|
" zeige(\"[1] M = Kapazitaet (richtig gewaehlt)\", knapp)\n",
|
|
"\n",
|
|
" gross = loese(1e7 * KAPAZITAET, presolve=\"off\")\n",
|
|
" zeige(\"[2] M = 10 Millionen mal Kapazitaet, Presolve abgeschaltet\", gross)\n",
|
|
"\n",
|
|
" print(\"\\n\" + \"=\" * 78)\n",
|
|
" print(\" WAS DA PASSIERT IST\")\n",
|
|
" print(\"=\" * 78)\n",
|
|
" fehlbetrag = knapp[\"zielwert\"] - gross[\"zielwert\"]\n",
|
|
" print(f\"Lauf [2] meldet {gross['zielwert']:,.2f} EUR und sieht damit um\")\n",
|
|
" print(f\"{fehlbetrag:,.2f} EUR BESSER aus als die richtige Loesung - ein Ergebnis,\")\n",
|
|
" print(\"ueber das sich jeder Auftraggeber freuen wuerde.\")\n",
|
|
" print()\n",
|
|
" print(\"Der Grund steht in der Zeile 'groesster y-Wert unter 0.5': Die\")\n",
|
|
" print(f\"Schaltervariablen stehen bei rund \"\n",
|
|
" f\"{gross['y'][gross['y'] < 0.5].max():.0e}.\")\n",
|
|
" print(f\"Das ist kleiner als die Ganzzahltoleranz {GANZZAHL_TOLERANZ:.0e}, also gilt\")\n",
|
|
" print(\"y = 0 - 'Lager geschlossen'. Zugleich ist M so gross, dass\")\n",
|
|
" print(\" sum_j x_ij <= M * 1e-8\")\n",
|
|
" print(\"immer noch reichlich Liefermenge erlaubt. Die Lager liefern also,\")\n",
|
|
" print(\"ohne dass ihre Fixkosten je bezahlt werden. Der Fachbegriff dafuer\")\n",
|
|
" print(\"ist 'trickle flow'.\")\n",
|
|
" print()\n",
|
|
" print(\"WICHTIG: Mit eingeschaltetem Presolve (Standard) faellt HiGHS hier\")\n",
|
|
" print(\"nicht darauf herein - es zieht M selbst zurecht. Verlassen Sie sich\")\n",
|
|
" print(\"nicht darauf: Presolve kann das nur, wenn eine implizite Schranke\")\n",
|
|
" print(\"herleitbar ist. Die Pruefung aus pruefe() kostet Millisekunden und\")\n",
|
|
" print(\"funktioniert immer.\")\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
|
|
}
|