906 lines
44 KiB
Text
906 lines
44 KiB
Text
|
{
|
||
|
|
"cells": [
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"# Kapitel 9: Metaheuristiken — wenn der exakte Solver aussteigt\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": [
|
||
|
|
"## Simulated Annealing\n",
|
||
|
|
"\n",
|
||
|
|
"`Simulated_Annealing.py`\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"#!/usr/bin/env python3\n",
|
||
|
|
"\n",
|
||
|
|
"# Simulated_Annealing.py\n",
|
||
|
|
"\"\"\"\n",
|
||
|
|
"Kapitel Metaheuristiken: Simulated Annealing an der Lackieranlage.\n",
|
||
|
|
"\n",
|
||
|
|
"Aufgabe: n Auftraege in eine Reihenfolge bringen. Zwischen zwei Auftraegen\n",
|
||
|
|
"faellt eine Ruestzeit an - die Anlage muss gereinigt werden. Wie lange das\n",
|
||
|
|
"dauert, haengt von BEIDEN Auftraegen ab: Ein Wechsel innerhalb derselben\n",
|
||
|
|
"Produktfamilie kostet fast nichts, ein Wechsel von Dunkel nach Hell ist teuer,\n",
|
||
|
|
"und zwischen manchen Familien ist der Reinigungsaufwand schlicht hoch. Gesucht\n",
|
||
|
|
"ist die Reihenfolge mit der kleinsten Summe der Ruestzeiten.\n",
|
||
|
|
"\n",
|
||
|
|
"Das Programm zeigt drei Dinge:\n",
|
||
|
|
"\n",
|
||
|
|
" 1. Warum die Faustregel \"immer der billigste naechste Auftrag\" in eine\n",
|
||
|
|
" Sackgasse laeuft - und wie weit man sie mit lokaler Suche verbessert.\n",
|
||
|
|
" 2. Wie gross der Beitrag der Verschlechterungen WIRKLICH ist. Die Antwort\n",
|
||
|
|
" faellt bescheidener aus als das Lehrbuch verspricht, und genau das ist\n",
|
||
|
|
" der Grund, den Vergleich immer mitzurechnen.\n",
|
||
|
|
" 3. Dass die Starttemperatur kein Beiwerk ist: zu heiss macht die Suche\n",
|
||
|
|
" nutzlos, und man sieht es an einer einzigen Kennzahl kommen.\n",
|
||
|
|
"\n",
|
||
|
|
"Zum Budget: Der Abkuehlplan laeuft ueber eine feste ZUGZAHL, nicht ueber eine\n",
|
||
|
|
"Sekundenzahl. Im Betrieb hat man zwar ein Zeitbudget - fuer einen VERGLEICH ist\n",
|
||
|
|
"die Zugzahl aber das ehrlichere Mass: Sie ist auf jeder Maschine dieselbe, und\n",
|
||
|
|
"die abgedruckten Zahlen unten lassen sich damit nachrechnen. Die gemessene\n",
|
||
|
|
"Laufzeit steht trotzdem dabei.\n",
|
||
|
|
"\n",
|
||
|
|
"Benoetigt: numpy\n",
|
||
|
|
"\"\"\"\n",
|
||
|
|
"\n",
|
||
|
|
"from __future__ import annotations\n",
|
||
|
|
"\n",
|
||
|
|
"import math\n",
|
||
|
|
"import time\n",
|
||
|
|
"\n",
|
||
|
|
"import numpy as np\n",
|
||
|
|
"\n",
|
||
|
|
"ZUEGE = 400_000 # Zugbudget je Lauf (statt Sekunden: reproduzierbar)\n",
|
||
|
|
"SAAT = 11 # feste Saat -> reproduzierbare Instanz\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"# --- 1. Die Instanz ---------------------------------------------------------\n",
|
||
|
|
"\n",
|
||
|
|
"def erzeuge_ruestmatrix(n: int, saat: int = SAAT) -> np.ndarray:\n",
|
||
|
|
" \"\"\"Ruestzeiten in Minuten zwischen je zwei Auftraegen.\n",
|
||
|
|
"\n",
|
||
|
|
" Die Zahl der Produktfamilien waechst mit der Auftragszahl (eine Familie je\n",
|
||
|
|
" zehn Auftraege). Sonst wuerde das Problem mit wachsendem n LEICHTER: Bei\n",
|
||
|
|
" fester Familienzahl haette der Plan irgendwann nur noch ein paar grosse\n",
|
||
|
|
" Bloecke, und jede Faustregel faende ihn.\n",
|
||
|
|
" \"\"\"\n",
|
||
|
|
" rng = np.random.default_rng(saat)\n",
|
||
|
|
" familien = max(6, n // 10)\n",
|
||
|
|
" # Reinigungsaufwand zwischen den Familien - unregelmaessig und asymmetrisch,\n",
|
||
|
|
" # so wie im Betrieb: Von Klarlack auf Rot ist etwas anderes als umgekehrt.\n",
|
||
|
|
" zwischen = rng.integers(8, 60, (familien, familien))\n",
|
||
|
|
" np.fill_diagonal(zwischen, 2)\n",
|
||
|
|
"\n",
|
||
|
|
" familie = rng.integers(0, familien, n)\n",
|
||
|
|
" farbe = rng.integers(0, 10, n) # 0 = weiss ... 9 = schwarz\n",
|
||
|
|
"\n",
|
||
|
|
" matrix = np.zeros((n, n), dtype=np.int64)\n",
|
||
|
|
" for i in range(n):\n",
|
||
|
|
" for j in range(n):\n",
|
||
|
|
" if i != j:\n",
|
||
|
|
" # Dunkel -> hell kostet zusaetzlich: die Anlage muss heller\n",
|
||
|
|
" # werden, das braucht mehr Spuelgaenge.\n",
|
||
|
|
" matrix[i, j] = (zwischen[familie[i], familie[j]]\n",
|
||
|
|
" + 3 * max(0, farbe[i] - farbe[j]))\n",
|
||
|
|
" return matrix\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def gesamtruestzeit(reihe: list[int], matrix: np.ndarray) -> int:\n",
|
||
|
|
" \"\"\"Summe der Ruestzeiten einer Reihenfolge (offene Kette, keine Rundreise).\"\"\"\n",
|
||
|
|
" return int(sum(matrix[reihe[k], reihe[k + 1]] for k in range(len(reihe) - 1)))\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def faustregel(matrix: np.ndarray) -> list[int]:\n",
|
||
|
|
" \"\"\"Immer der billigste noch offene Auftrag - 'naechster Nachbar'.\n",
|
||
|
|
"\n",
|
||
|
|
" So plant ein Meister von Hand, und es ist keine schlechte Regel. Ihr Fehler\n",
|
||
|
|
" ist die Kurzsichtigkeit: Sie spart am Anfang und laesst die teuren Wechsel\n",
|
||
|
|
" fuer das Ende uebrig, wo keine Wahl mehr bleibt.\n",
|
||
|
|
" \"\"\"\n",
|
||
|
|
" n = len(matrix)\n",
|
||
|
|
" offen = set(range(1, n))\n",
|
||
|
|
" reihe = [0]\n",
|
||
|
|
" while offen:\n",
|
||
|
|
" naechster = min(offen, key=lambda j: matrix[reihe[-1], j])\n",
|
||
|
|
" reihe.append(naechster)\n",
|
||
|
|
" offen.discard(naechster)\n",
|
||
|
|
" return reihe\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"# --- 2. Der Zug und seine Kostenaenderung -----------------------------------\n",
|
||
|
|
"\n",
|
||
|
|
"def delta_verschieben(reihe: list[int], matrix: np.ndarray,\n",
|
||
|
|
" von: int, nach: int):\n",
|
||
|
|
" \"\"\"Auftrag von Position 'von' an Position 'nach' versetzen.\n",
|
||
|
|
"\n",
|
||
|
|
" Liefert (Kostenaenderung, Restliste, Einfuegeposition) oder None, wenn der\n",
|
||
|
|
" Zug nichts aendert.\n",
|
||
|
|
"\n",
|
||
|
|
" DAS IST DIE WICHTIGSTE FUNKTION DES PROGRAMMS. Sie berechnet die\n",
|
||
|
|
" Kostenaenderung aus HOECHSTENS SECHS Matrixeintraegen, statt die\n",
|
||
|
|
" Reihenfolge neu durchzusummieren. Der Unterschied ist nicht kosmetisch:\n",
|
||
|
|
" Die volle Summe kostet O(n) je Zug, diese Rechnung O(1). Bei 500 Auftraegen\n",
|
||
|
|
" sind das rund 500-mal mehr geprueste Zuege im selben Zeitbudget - und eine\n",
|
||
|
|
" Metaheuristik lebt von der Zahl der Zuege.\n",
|
||
|
|
"\n",
|
||
|
|
" Herausgerissen wird der Auftrag aus zwei Kanten, seine Nachbarn ruecken\n",
|
||
|
|
" zusammen (eine neue Kante). Eingefuegt wird er zwischen zwei andere\n",
|
||
|
|
" Nachbarn, deren bisherige Kante dadurch verschwindet.\n",
|
||
|
|
" \"\"\"\n",
|
||
|
|
" n = len(reihe)\n",
|
||
|
|
" if von == nach or nach == von + 1:\n",
|
||
|
|
" return None\n",
|
||
|
|
"\n",
|
||
|
|
" heraus = 0\n",
|
||
|
|
" if von > 0:\n",
|
||
|
|
" heraus += matrix[reihe[von - 1], reihe[von]]\n",
|
||
|
|
" if von < n - 1:\n",
|
||
|
|
" heraus += matrix[reihe[von], reihe[von + 1]]\n",
|
||
|
|
" if 0 < von < n - 1: # die Nachbarn ruecken zusammen\n",
|
||
|
|
" heraus -= matrix[reihe[von - 1], reihe[von + 1]]\n",
|
||
|
|
"\n",
|
||
|
|
" rest = reihe[:von] + reihe[von + 1:]\n",
|
||
|
|
" stelle = nach if nach < von else nach - 1\n",
|
||
|
|
"\n",
|
||
|
|
" hinein = 0\n",
|
||
|
|
" if stelle > 0:\n",
|
||
|
|
" hinein += matrix[rest[stelle - 1], reihe[von]]\n",
|
||
|
|
" if stelle < len(rest):\n",
|
||
|
|
" hinein += matrix[reihe[von], rest[stelle]]\n",
|
||
|
|
" if 0 < stelle < len(rest): # die aufgetrennte Kante faellt weg\n",
|
||
|
|
" hinein -= matrix[rest[stelle - 1], rest[stelle]]\n",
|
||
|
|
"\n",
|
||
|
|
" return int(hinein - heraus), rest, stelle\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"# --- 3. Die Suche -----------------------------------------------------------\n",
|
||
|
|
"\n",
|
||
|
|
"class Verlauf:\n",
|
||
|
|
" \"\"\"Was ein Lauf ausser dem Ergebnis noch verraet.\"\"\"\n",
|
||
|
|
"\n",
|
||
|
|
" def __init__(self) -> None:\n",
|
||
|
|
" self.schritte = 0\n",
|
||
|
|
" self.verschlechterungen_geprueft = 0\n",
|
||
|
|
" self.verschlechterungen_genommen = 0\n",
|
||
|
|
" self.schlechtester_zustand = 0 # wie weit die Suche abgedriftet ist\n",
|
||
|
|
" self.sekunden = 0.0 # nur zur Information, nicht zum Vergleich\n",
|
||
|
|
"\n",
|
||
|
|
" @property\n",
|
||
|
|
" def annahmequote(self) -> float:\n",
|
||
|
|
" if not self.verschlechterungen_geprueft:\n",
|
||
|
|
" return 0.0\n",
|
||
|
|
" return self.verschlechterungen_genommen / self.verschlechterungen_geprueft\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def suche(matrix: np.ndarray, zuege: int, start_temperatur: float,\n",
|
||
|
|
" end_temperatur: float = 0.05, saat: int = 1,\n",
|
||
|
|
" bergsteigen: bool = False) -> tuple[list[int], int, Verlauf]:\n",
|
||
|
|
" \"\"\"Simulated Annealing (oder reines Bergsteigen, wenn bergsteigen=True).\n",
|
||
|
|
"\n",
|
||
|
|
" Der einzige Unterschied zwischen beiden steht in der Annahmezeile: Das\n",
|
||
|
|
" Bergsteigen nimmt nur Verbesserungen, SA nimmt Verschlechterungen mit einer\n",
|
||
|
|
" Wahrscheinlichkeit an, die mit der Zeit gegen null geht.\n",
|
||
|
|
" \"\"\"\n",
|
||
|
|
" rng = np.random.default_rng(saat)\n",
|
||
|
|
" n = len(matrix)\n",
|
||
|
|
" reihe = faustregel(matrix)\n",
|
||
|
|
" kosten = gesamtruestzeit(reihe, matrix)\n",
|
||
|
|
" beste, beste_kosten = reihe[:], kosten\n",
|
||
|
|
" verlauf = Verlauf()\n",
|
||
|
|
" verlauf.schlechtester_zustand = kosten\n",
|
||
|
|
"\n",
|
||
|
|
" start = time.perf_counter()\n",
|
||
|
|
" for zug in range(zuege):\n",
|
||
|
|
" verlauf.schritte += 1\n",
|
||
|
|
" # Geometrischer Abkuehlplan ueber das Zugbudget.\n",
|
||
|
|
" temperatur = start_temperatur * (end_temperatur / start_temperatur) \\\n",
|
||
|
|
" ** (zug / zuege)\n",
|
||
|
|
"\n",
|
||
|
|
" von = int(rng.integers(0, n))\n",
|
||
|
|
" nach = int(rng.integers(0, n + 1))\n",
|
||
|
|
" ergebnis = delta_verschieben(reihe, matrix, von, nach)\n",
|
||
|
|
" if ergebnis is None:\n",
|
||
|
|
" continue\n",
|
||
|
|
" aenderung, rest, stelle = ergebnis\n",
|
||
|
|
"\n",
|
||
|
|
" if aenderung <= 0:\n",
|
||
|
|
" annehmen = True\n",
|
||
|
|
" else:\n",
|
||
|
|
" verlauf.verschlechterungen_geprueft += 1\n",
|
||
|
|
" # Metropolis-Kriterium: je groesser die Verschlechterung und je\n",
|
||
|
|
" # kaelter es ist, desto unwahrscheinlicher.\n",
|
||
|
|
" annehmen = (not bergsteigen\n",
|
||
|
|
" and rng.random() < math.exp(-aenderung / temperatur))\n",
|
||
|
|
" verlauf.verschlechterungen_genommen += annehmen\n",
|
||
|
|
"\n",
|
||
|
|
" if annehmen:\n",
|
||
|
|
" reihe = rest[:stelle] + [reihe[von]] + rest[stelle:]\n",
|
||
|
|
" kosten += aenderung\n",
|
||
|
|
" verlauf.schlechtester_zustand = max(verlauf.schlechtester_zustand,\n",
|
||
|
|
" kosten)\n",
|
||
|
|
" if kosten < beste_kosten:\n",
|
||
|
|
" beste, beste_kosten = reihe[:], kosten\n",
|
||
|
|
"\n",
|
||
|
|
" verlauf.sekunden = time.perf_counter() - start\n",
|
||
|
|
" return beste, beste_kosten, verlauf\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"if __name__ == \"__main__\":\n",
|
||
|
|
" N = 200\n",
|
||
|
|
" matrix = erzeuge_ruestmatrix(N)\n",
|
||
|
|
" start_reihe = faustregel(matrix)\n",
|
||
|
|
" start_kosten = gesamtruestzeit(start_reihe, matrix)\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"=\" * 78)\n",
|
||
|
|
" print(\" SIMULATED ANNEALING AN DER LACKIERANLAGE\")\n",
|
||
|
|
" print(\"=\" * 78)\n",
|
||
|
|
" print(f\"{N} Auftraege, {max(6, N // 10)} Produktfamilien, \"\n",
|
||
|
|
" f\"Zugbudget {ZUEGE:,} je Lauf.\\n\")\n",
|
||
|
|
" print(f\"Faustregel 'billigster naechster Auftrag': \"\n",
|
||
|
|
" f\"{start_kosten:,} Minuten Ruestzeit\")\n",
|
||
|
|
"\n",
|
||
|
|
" # --- Wie gross sind die Zuege ueberhaupt? ---------------------------\n",
|
||
|
|
" # Ohne diese Zahl kann man die Temperatur nicht waehlen: Das\n",
|
||
|
|
" # Metropolis-Kriterium vergleicht die Verschlechterung MIT der Temperatur.\n",
|
||
|
|
" rng = np.random.default_rng(0)\n",
|
||
|
|
" stichprobe = []\n",
|
||
|
|
" for _ in range(5000):\n",
|
||
|
|
" e = delta_verschieben(start_reihe, matrix,\n",
|
||
|
|
" int(rng.integers(0, N)), int(rng.integers(0, N + 1)))\n",
|
||
|
|
" if e:\n",
|
||
|
|
" stichprobe.append(abs(e[0]))\n",
|
||
|
|
" print(f\"Typische Zuggroesse |Delta|: Median {np.median(stichprobe):.0f}, \"\n",
|
||
|
|
" f\"90 %-Quantil {np.quantile(stichprobe, 0.9):.0f} Minuten\")\n",
|
||
|
|
"\n",
|
||
|
|
" # --- Bergsteigen gegen Annealing -------------------------------------\n",
|
||
|
|
" print(\"\\n\" + \"-\" * 78)\n",
|
||
|
|
" print(\"Nur Verbesserungen annehmen (Bergsteigen) gegen Annealing:\\n\")\n",
|
||
|
|
" _, berg, berg_verlauf = suche(matrix, ZUEGE, 1.0, bergsteigen=True)\n",
|
||
|
|
" _, gluehen, gluehen_verlauf = suche(matrix, ZUEGE, 1.0)\n",
|
||
|
|
" print(f\" {'Faustregel (Start)':<34} {start_kosten:>7,} Minuten\")\n",
|
||
|
|
" print(f\" {'Bergsteigen':<34} {berg:>7,} Minuten \"\n",
|
||
|
|
" f\"({(start_kosten - berg) / start_kosten * 100:5.1f} % besser)\")\n",
|
||
|
|
" print(f\" {'Simulated Annealing':<34} {gluehen:>7,} Minuten \"\n",
|
||
|
|
" f\"({(start_kosten - gluehen) / start_kosten * 100:5.1f} % besser)\")\n",
|
||
|
|
" print(f\"\\n Beide pruefen genau {berg_verlauf.schritte:,} Zuege \"\n",
|
||
|
|
" f\"({berg_verlauf.sekunden:.1f} bzw. {gluehen_verlauf.sekunden:.1f} Sekunden).\")\n",
|
||
|
|
" print(f\" Annealing nimmt davon {gluehen_verlauf.annahmequote * 100:.2f} % der \"\n",
|
||
|
|
" f\"VERSCHLECHTERUNGEN an, Bergsteigen keine.\")\n",
|
||
|
|
" print(\"\\n Bemerkenswert ist, wie klein der Abstand ist: Das simple\")\n",
|
||
|
|
" print(f\" Bergsteigen holt hier bereits {(start_kosten - berg) / (start_kosten - gluehen) * 100:.0f} % dessen, was Annealing\")\n",
|
||
|
|
" print(\" schafft. Wer eine Metaheuristik einsetzt, sollte diesen Vergleich\")\n",
|
||
|
|
" print(\" IMMER mitrechnen - sonst schreibt man einer aufwendigen Methode\")\n",
|
||
|
|
" print(\" gut, was schon die einfache geliefert haette.\")\n",
|
||
|
|
"\n",
|
||
|
|
" # --- Die Temperatur ist kein Beiwerk ---------------------------------\n",
|
||
|
|
" print(\"\\n\" + \"-\" * 78)\n",
|
||
|
|
" print(\"Und warum die Starttemperatur ausgemessen gehoert:\\n\")\n",
|
||
|
|
" print(f\" {'Start-T':>8} {'Ruestzeit':>11} {'gg. Faustregel':>15} \"\n",
|
||
|
|
" f\"{'angenommene':>12} {'schlechtester':>14}\")\n",
|
||
|
|
" print(f\" {'':>8} {'':>11} {'':>15} {'Verschlecht.':>12} {'Zwischenwert':>14}\")\n",
|
||
|
|
" print(\" \" + \"-\" * 64)\n",
|
||
|
|
" ergebnisse, verlaeufe = {}, {}\n",
|
||
|
|
" for temperatur in (0.5, 1.0, 2.0, 4.0, 8.0):\n",
|
||
|
|
" _, wert, verlauf = suche(matrix, ZUEGE, temperatur)\n",
|
||
|
|
" ergebnisse[temperatur] = wert\n",
|
||
|
|
" verlaeufe[temperatur] = verlauf\n",
|
||
|
|
" gewinn = (start_kosten - wert) / start_kosten * 100\n",
|
||
|
|
" print(f\" {temperatur:>8.1f} {wert:>11,} {gewinn:>13.1f} % \"\n",
|
||
|
|
" f\"{verlauf.annahmequote * 100:>11.2f} % \"\n",
|
||
|
|
" f\"{verlauf.schlechtester_zustand:>14,}\")\n",
|
||
|
|
"\n",
|
||
|
|
" beste_temperatur = min(ergebnisse, key=ergebnisse.get)\n",
|
||
|
|
" heiss = verlaeufe[8.0]\n",
|
||
|
|
" print(f\"\\n Bester Wert bei T0 = {beste_temperatur}.\")\n",
|
||
|
|
" print(\"\\n Die letzten beiden Spalten erklaeren, warum es nach oben kippt.\")\n",
|
||
|
|
" print(f\" Bei T0 = 8 gehen nur {heiss.annahmequote * 100:.2f} % der Verschlechterungen\")\n",
|
||
|
|
" print(\" durch - wenig genug, dass man es fuer harmlos halten koennte. Es\")\n",
|
||
|
|
" print(f\" genuegt aber: Die Suche driftet bis auf {heiss.schlechtester_zustand:,} Minuten ab, das\")\n",
|
||
|
|
" print(f\" {heiss.schlechtester_zustand / start_kosten:.1f}-fache der Startloesung, und findet im Zugbudget nicht\")\n",
|
||
|
|
" print(\" mehr zurueck. Vom Ergebnis bleibt fast nichts uebrig: Die Suche hat\")\n",
|
||
|
|
" print(f\" in {ZUEGE:,} Zuegen nie etwas GESEHEN, das mehr als\")\n",
|
||
|
|
" print(f\" {(start_kosten - ergebnisse[8.0]) / start_kosten * 100:.1f} % unter der Startloesung lag.\")\n",
|
||
|
|
" print(\"\\n Bemerkenswert ist die andere Richtung: Die besten Werte entstehen\")\n",
|
||
|
|
" print(\" bei Annahmequoten nahe null. Die Lehrbuchregel 'anfangs sollen\")\n",
|
||
|
|
" print(\" 20 bis 50 Prozent der Verschlechterungen durchgehen' fuehrt hier in\")\n",
|
||
|
|
" print(\" die Irre - sie stammt aus Problemen mit sehr kleinen Zuggroessen.\")\n",
|
||
|
|
" print(\" Bei Ruestzeiten in Minuten ist ein einziger schlechter Zug teuer.\")\n",
|
||
|
|
" print(\"\\n Was stattdessen traegt: Den SCHLECHTESTEN Zwischenwert protokollieren\")\n",
|
||
|
|
" print(\" und T0 so waehlen, dass er die Startloesung nicht wesentlich\")\n",
|
||
|
|
" print(\" ueberschreitet. Diese eine Zahl haette hier auf Anhieb gezeigt, dass\")\n",
|
||
|
|
" print(\" T0 = 8 unbrauchbar ist - ohne dass man das Endergebnis abwarten muss.\")\n",
|
||
|
|
" print(\"=\" * 78)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## Ab wann lohnt es sich?\n",
|
||
|
|
"\n",
|
||
|
|
"`Metaheuristik_vs_Exakt.py`\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"#!/usr/bin/env python3\n",
|
||
|
|
"\n",
|
||
|
|
"# Metaheuristik_vs_Exakt.py\n",
|
||
|
|
"\"\"\"\n",
|
||
|
|
"Kapitel Metaheuristiken: Ab wann lohnt sich eine Metaheuristik? Die Antwort ist eine Zahl.\n",
|
||
|
|
"\n",
|
||
|
|
"Die verbreitete Regel lautet: \"Kleine Probleme exakt, grosse mit\n",
|
||
|
|
"Metaheuristiken.\" Sie ist richtig - aber voellig nutzlos, solange niemand sagt,\n",
|
||
|
|
"wo 'gross' anfaengt. Dieses Programm misst es fuer die Ruestzeitenaufgabe aus\n",
|
||
|
|
"Simulated_Annealing.py nach, und das Ergebnis ueberrascht in beide Richtungen:\n",
|
||
|
|
"\n",
|
||
|
|
" * Bei 60 Auftraegen loest CP-SAT BEWEISBAR optimal, in wenigen Sekunden.\n",
|
||
|
|
" Die Metaheuristik ist hier schlicht die schlechtere Wahl - sie liefert ein\n",
|
||
|
|
" schwaecheres Ergebnis und weiss nicht einmal, wie schwach.\n",
|
||
|
|
" * Bei 200 Auftraegen liegt CP-SAT ohne Optimalitaetsbeweis vorne.\n",
|
||
|
|
" * Erst bei 500 Auftraegen kippt es: Dort gibt CP-SAT im Zeitlimit eine\n",
|
||
|
|
" Loesung aus, die SCHLECHTER ist als die Faustregel eines Meisters.\n",
|
||
|
|
"\n",
|
||
|
|
"Der Punkt des Kapitels steht in der letzten Spalte: Die Schranke. Nur der\n",
|
||
|
|
"exakte Solver sagt, wie gut die Loesung sein KANN. Eine Metaheuristik allein\n",
|
||
|
|
"liefert eine Zahl ohne Massstab - man weiss nie, ob man 2 % oder 40 % neben dem\n",
|
||
|
|
"Optimum liegt.\n",
|
||
|
|
"\n",
|
||
|
|
"ZUR REPRODUZIERBARKEIT: Die Spalte 'Annealing' ist auf jeder Maschine gleich -\n",
|
||
|
|
"sie laeuft ueber ein festes ZUGbudget. Die CP-SAT-Spalten sind es nicht: Ein\n",
|
||
|
|
"Zeitlimit ist ein Wanduhr-Limit, und wie weit der Solver in 30 Sekunden kommt,\n",
|
||
|
|
"haengt von der Maschine ab. CP-SAT kennt zwar ein deterministisches Zeitmass\n",
|
||
|
|
"(max_deterministic_time), das aber in einer Groessenordnung rechnet, die fuer\n",
|
||
|
|
"ein Buchbeispiel unbrauchbar langsam ist. Die Zahlen unten sind also\n",
|
||
|
|
"hardwareabhaengig - die AUSSAGE der Tabelle ist es nicht, sie ist auf\n",
|
||
|
|
"verschiedenen Groessen und in mehreren Laeufen stabil.\n",
|
||
|
|
"\n",
|
||
|
|
"WICHTIG: Hier wird nur ortools importiert, nicht highspy (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",
|
||
|
|
"ZUEGE = 400_000 # Zugbudget der Metaheuristik (wie Simulated_Annealing.py)\n",
|
||
|
|
"ZEITLIMIT = 30.0 # Zeitbudget des exakten Solvers\n",
|
||
|
|
"GROESSEN = (60, 200, 500)\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"# --- Instanz, Faustregel und Suche: identisch zu Simulated_Annealing.py -----\n",
|
||
|
|
"\n",
|
||
|
|
"def erzeuge_ruestmatrix(n: int, saat: int = 11) -> np.ndarray:\n",
|
||
|
|
" rng = np.random.default_rng(saat)\n",
|
||
|
|
" familien = max(6, n // 10)\n",
|
||
|
|
" zwischen = rng.integers(8, 60, (familien, familien))\n",
|
||
|
|
" np.fill_diagonal(zwischen, 2)\n",
|
||
|
|
" familie = rng.integers(0, familien, n)\n",
|
||
|
|
" farbe = rng.integers(0, 10, n)\n",
|
||
|
|
" matrix = np.zeros((n, n), dtype=np.int64)\n",
|
||
|
|
" for i in range(n):\n",
|
||
|
|
" for j in range(n):\n",
|
||
|
|
" if i != j:\n",
|
||
|
|
" matrix[i, j] = (zwischen[familie[i], familie[j]]\n",
|
||
|
|
" + 3 * max(0, farbe[i] - farbe[j]))\n",
|
||
|
|
" return matrix\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def gesamtruestzeit(reihe: list[int], matrix: np.ndarray) -> int:\n",
|
||
|
|
" return int(sum(matrix[reihe[k], reihe[k + 1]] for k in range(len(reihe) - 1)))\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def faustregel(matrix: np.ndarray) -> list[int]:\n",
|
||
|
|
" n = len(matrix)\n",
|
||
|
|
" offen = set(range(1, n))\n",
|
||
|
|
" reihe = [0]\n",
|
||
|
|
" while offen:\n",
|
||
|
|
" naechster = min(offen, key=lambda j: matrix[reihe[-1], j])\n",
|
||
|
|
" reihe.append(naechster)\n",
|
||
|
|
" offen.discard(naechster)\n",
|
||
|
|
" return reihe\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def delta_verschieben(reihe, matrix, von, nach):\n",
|
||
|
|
" n = len(reihe)\n",
|
||
|
|
" if von == nach or nach == von + 1:\n",
|
||
|
|
" return None\n",
|
||
|
|
" heraus = 0\n",
|
||
|
|
" if von > 0:\n",
|
||
|
|
" heraus += matrix[reihe[von - 1], reihe[von]]\n",
|
||
|
|
" if von < n - 1:\n",
|
||
|
|
" heraus += matrix[reihe[von], reihe[von + 1]]\n",
|
||
|
|
" if 0 < von < n - 1:\n",
|
||
|
|
" heraus -= matrix[reihe[von - 1], reihe[von + 1]]\n",
|
||
|
|
" rest = reihe[:von] + reihe[von + 1:]\n",
|
||
|
|
" stelle = nach if nach < von else nach - 1\n",
|
||
|
|
" hinein = 0\n",
|
||
|
|
" if stelle > 0:\n",
|
||
|
|
" hinein += matrix[rest[stelle - 1], reihe[von]]\n",
|
||
|
|
" if stelle < len(rest):\n",
|
||
|
|
" hinein += matrix[reihe[von], rest[stelle]]\n",
|
||
|
|
" if 0 < stelle < len(rest):\n",
|
||
|
|
" hinein -= matrix[rest[stelle - 1], rest[stelle]]\n",
|
||
|
|
" return int(hinein - heraus), rest, stelle\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def annealing(matrix: np.ndarray, zuege: int = ZUEGE,\n",
|
||
|
|
" start_temperatur: float = 1.0, saat: int = 1) -> tuple[int, float]:\n",
|
||
|
|
" \"\"\"Dieselbe Suche wie in Simulated_Annealing.py, hier nur als Vergleichswert.\"\"\"\n",
|
||
|
|
" import math\n",
|
||
|
|
"\n",
|
||
|
|
" rng = np.random.default_rng(saat)\n",
|
||
|
|
" n = len(matrix)\n",
|
||
|
|
" reihe = faustregel(matrix)\n",
|
||
|
|
" kosten = gesamtruestzeit(reihe, matrix)\n",
|
||
|
|
" beste_kosten = kosten\n",
|
||
|
|
"\n",
|
||
|
|
" t0 = time.perf_counter()\n",
|
||
|
|
" for zug in range(zuege):\n",
|
||
|
|
" temperatur = start_temperatur * (0.05 / start_temperatur) ** (zug / zuege)\n",
|
||
|
|
" von = int(rng.integers(0, n))\n",
|
||
|
|
" nach = int(rng.integers(0, n + 1))\n",
|
||
|
|
" ergebnis = delta_verschieben(reihe, matrix, von, nach)\n",
|
||
|
|
" if ergebnis is None:\n",
|
||
|
|
" continue\n",
|
||
|
|
" aenderung, rest, stelle = ergebnis\n",
|
||
|
|
" if aenderung <= 0 or rng.random() < math.exp(-aenderung / temperatur):\n",
|
||
|
|
" reihe = rest[:stelle] + [reihe[von]] + rest[stelle:]\n",
|
||
|
|
" kosten += aenderung\n",
|
||
|
|
" beste_kosten = min(beste_kosten, kosten)\n",
|
||
|
|
" return beste_kosten, time.perf_counter() - t0\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"# --- Der exakte Weg: CP-SAT mit AddCircuit ---------------------------------\n",
|
||
|
|
"\n",
|
||
|
|
"def exakt(matrix: np.ndarray, zeitlimit: float = ZEITLIMIT):\n",
|
||
|
|
" \"\"\"Reihenfolge als Rundreise mit einem kostenlosen Hilfsknoten.\n",
|
||
|
|
"\n",
|
||
|
|
" AddCircuit verlangt einen geschlossenen Kreis. Unsere Aufgabe ist aber eine\n",
|
||
|
|
" offene Kette: Der erste Auftrag hat keinen Vorgaenger, der letzte keinen\n",
|
||
|
|
" Nachfolger. Der uebliche Kniff ist ein zusaetzlicher Hilfsknoten n, dessen\n",
|
||
|
|
" Kanten nichts kosten. Der Kreis laeuft dann Hilfsknoten -> Kette ->\n",
|
||
|
|
" Hilfsknoten, und die Kosten des Kreises sind genau die der Kette.\n",
|
||
|
|
"\n",
|
||
|
|
" AddCircuit ist der Grund, warum CP-SAT hier so lange mithaelt: Die\n",
|
||
|
|
" Nebenbedingung sorgt selbst dafuer, dass keine Teilkreise entstehen. Wer\n",
|
||
|
|
" dasselbe als MILP formuliert, braucht dafuer entweder exponentiell viele\n",
|
||
|
|
" Ungleichungen oder die schwache MTZ-Formulierung (Kapitel Graphen).\n",
|
||
|
|
" \"\"\"\n",
|
||
|
|
" n = len(matrix)\n",
|
||
|
|
" modell = cp_model.CpModel()\n",
|
||
|
|
" kanten, ziel = [], []\n",
|
||
|
|
" for i in range(n + 1):\n",
|
||
|
|
" for j in range(n + 1):\n",
|
||
|
|
" if i == j:\n",
|
||
|
|
" continue\n",
|
||
|
|
" aktiv = modell.NewBoolVar(f\"kante_{i}_{j}\")\n",
|
||
|
|
" kanten.append((i, j, aktiv))\n",
|
||
|
|
" if i < n and j < n: # Hilfsknoten kostet nichts\n",
|
||
|
|
" ziel.append(int(matrix[i, j]) * aktiv)\n",
|
||
|
|
" modell.AddCircuit(kanten)\n",
|
||
|
|
" modell.Minimize(sum(ziel))\n",
|
||
|
|
"\n",
|
||
|
|
" loeser = cp_model.CpSolver()\n",
|
||
|
|
" loeser.parameters.max_time_in_seconds = zeitlimit\n",
|
||
|
|
" loeser.parameters.num_workers = 8\n",
|
||
|
|
" loeser.parameters.random_seed = 1\n",
|
||
|
|
" t0 = time.perf_counter()\n",
|
||
|
|
" status = loeser.Solve(modell)\n",
|
||
|
|
" dauer = time.perf_counter() - t0\n",
|
||
|
|
" return (loeser.StatusName(status), int(loeser.ObjectiveValue()),\n",
|
||
|
|
" int(loeser.BestObjectiveBound()), dauer)\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"if __name__ == \"__main__\":\n",
|
||
|
|
" print(\"=\" * 88)\n",
|
||
|
|
" print(\" AB WANN LOHNT SICH DIE METAHEURISTIK?\")\n",
|
||
|
|
" print(\"=\" * 88)\n",
|
||
|
|
" print(f\"Dieselbe Aufgabe in drei Groessen. Metaheuristik: {ZUEGE:,} Zuege. \"\n",
|
||
|
|
" f\"Exakt: {ZEITLIMIT:.0f} s Zeitlimit.\\n\")\n",
|
||
|
|
"\n",
|
||
|
|
" print(f\"{'Auftraege':>10} {'Faustregel':>11} {'Annealing':>10} \"\n",
|
||
|
|
" f\"{'CP-SAT':>9} {'Status':>10} {'Schranke':>9} {'CP-Zeit':>9}\")\n",
|
||
|
|
" print(\"-\" * 88)\n",
|
||
|
|
"\n",
|
||
|
|
" zeilen = []\n",
|
||
|
|
" for n in GROESSEN:\n",
|
||
|
|
" matrix = erzeuge_ruestmatrix(n)\n",
|
||
|
|
" start = gesamtruestzeit(faustregel(matrix), matrix)\n",
|
||
|
|
" heuristisch, _ = annealing(matrix)\n",
|
||
|
|
" status, wert, schranke, dauer = exakt(matrix)\n",
|
||
|
|
" zeilen.append((n, start, heuristisch, wert, status, schranke, dauer))\n",
|
||
|
|
" print(f\"{n:>10} {start:>11,} {heuristisch:>10,} {wert:>9,} \"\n",
|
||
|
|
" f\"{status:>10} {schranke:>9,} {dauer:>8.1f}s\")\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"-\" * 88)\n",
|
||
|
|
" print(\"\\nWas in jeder Zeile steht:\\n\")\n",
|
||
|
|
" for n, start, heur, wert, status, schranke, _ in zeilen:\n",
|
||
|
|
" if status == \"OPTIMAL\":\n",
|
||
|
|
" urteil = (f\"CP-SAT beweist das Optimum. Annealing liegt \"\n",
|
||
|
|
" f\"{(heur - wert) / wert * 100:.1f} % daneben - und haette\")\n",
|
||
|
|
" zweite = \" das ohne den exakten Lauf nicht erfahren.\"\n",
|
||
|
|
" elif wert <= heur:\n",
|
||
|
|
" urteil = (f\"CP-SAT liegt {(heur - wert) / heur * 100:.1f} % vor Annealing, \"\n",
|
||
|
|
" f\"ohne Beweis (Gap {(wert - schranke) / wert * 100:.1f} %).\")\n",
|
||
|
|
" zweite = \" Auch ohne Optimalitaetsbeweis ist es die bessere Wahl.\"\n",
|
||
|
|
" else:\n",
|
||
|
|
" urteil = (f\"UMSCHLAGPUNKT: CP-SAT ist {(wert - heur) / heur * 100:.1f} % \"\n",
|
||
|
|
" f\"SCHLECHTER als Annealing\")\n",
|
||
|
|
" zweite = (f\" und sogar {(wert - start) / start * 100:.1f} % schlechter \"\n",
|
||
|
|
" f\"als die Faustregel des Meisters.\")\n",
|
||
|
|
" print(f\" {n:>4} Auftraege: {urteil}\")\n",
|
||
|
|
" print(zweite)\n",
|
||
|
|
"\n",
|
||
|
|
" n, start, heur, wert, status, schranke, _ = zeilen[-1]\n",
|
||
|
|
" print(\"\\n\" + \"=\" * 88)\n",
|
||
|
|
" print(\" DIE SPALTE, DIE MAN NICHT WEGLASSEN DARF\")\n",
|
||
|
|
" print(\"=\" * 88)\n",
|
||
|
|
" print(\"Die Schranke ist der eigentliche Ertrag des exakten Solvers - auch dann,\")\n",
|
||
|
|
" print(\"wenn seine LOESUNG unbrauchbar ist.\")\n",
|
||
|
|
" print()\n",
|
||
|
|
" print(f\"Bei {n} Auftraegen liefert Annealing {heur:,} Minuten. Klingt gut: \"\n",
|
||
|
|
" f\"{(start - heur) / start * 100:.1f} %\")\n",
|
||
|
|
" print(\"besser als die Faustregel. Die Schranke aus demselben CP-SAT-Lauf, dessen\")\n",
|
||
|
|
" print(f\"Loesung wir gerade verworfen haben, sagt aber: Unter {schranke:,} Minuten geht\")\n",
|
||
|
|
" print(f\"es nicht. Zwischen {schranke:,} und {heur:,} liegen {(heur - schranke) / heur * 100:.0f} % - so viel Luft\")\n",
|
||
|
|
" print(\"kann noch nach oben sein.\")\n",
|
||
|
|
" print()\n",
|
||
|
|
" print(\"Das ist der Unterschied zwischen 'besser als vorher' und 'gut'. Eine\")\n",
|
||
|
|
" print(\"Metaheuristik allein kann nur das Erste. Deshalb gehoert auch dann ein\")\n",
|
||
|
|
" print(\"exakter Lauf dazu, wenn man am Ende die heuristische Loesung einsetzt:\")\n",
|
||
|
|
" print(\"nicht wegen seiner Loesung, sondern wegen seiner Schranke.\")\n",
|
||
|
|
" print(\"=\" * 88)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## Large Neighborhood Search\n",
|
||
|
|
"\n",
|
||
|
|
"`Large_Neighborhood_Search.py`\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"#!/usr/bin/env python3\n",
|
||
|
|
"\n",
|
||
|
|
"# Large_Neighborhood_Search.py\n",
|
||
|
|
"\"\"\"\n",
|
||
|
|
"Kapitel Metaheuristiken: Zerstoeren und exakt reparieren - der Solver im Dienst der Heuristik.\n",
|
||
|
|
"\n",
|
||
|
|
"Metaheuristik_vs_Exakt.py endet mit einem Widerspruch. Bei 500 Auftraegen gilt:\n",
|
||
|
|
"\n",
|
||
|
|
" * Der exakte Solver scheitert an der GROESSE - seine Loesung nach 30 Sekunden\n",
|
||
|
|
" ist schlechter als die Faustregel eines Meisters.\n",
|
||
|
|
" * Die Metaheuristik kommt weiter, laesst aber laut Schranke noch rund ein\n",
|
||
|
|
" Viertel liegen. Ihre Zuege sind zu kleinteilig: Sie verschiebt jeweils\n",
|
||
|
|
" EINEN Auftrag und kann eine ganze Passage nicht auf einmal umbauen.\n",
|
||
|
|
"\n",
|
||
|
|
"Large Neighborhood Search loest den Widerspruch, indem sie beide einsetzt -\n",
|
||
|
|
"jeden fuer das, was er kann:\n",
|
||
|
|
"\n",
|
||
|
|
" ZERSTOEREN Ein Stueck des Plans herausbrechen (hier: ein zusammen-\n",
|
||
|
|
" haengendes Fenster von 20 Auftraegen).\n",
|
||
|
|
" REPARIEREN Genau dieses Stueck EXAKT neu optimieren. 20 Auftraege sind\n",
|
||
|
|
" fuer CP-SAT eine Kleinigkeit, 500 sind es nicht.\n",
|
||
|
|
" UEBERNEHMEN Nur behalten, wenn der Gesamtplan besser wurde.\n",
|
||
|
|
"\n",
|
||
|
|
"Der Solver bekommt also nicht mehr das ganze Problem, sondern immer wieder ein\n",
|
||
|
|
"kleines. Das ist der ganze Trick, und er ist in der Praxis der wichtigste\n",
|
||
|
|
"Baustein dieses Kapitels.\n",
|
||
|
|
"\n",
|
||
|
|
"ZUR REPRODUZIERBARKEIT: Alle Laeufe haben eine feste RUNDENzahl, nicht ein\n",
|
||
|
|
"Zeitbudget - die Ergebniswerte sind damit auf jeder Maschine gleich. Die\n",
|
||
|
|
"Laufzeiten daneben sind hardwareabhaengig.\n",
|
||
|
|
"\n",
|
||
|
|
"Benoetigt: numpy, ortools\n",
|
||
|
|
"\"\"\"\n",
|
||
|
|
"\n",
|
||
|
|
"from __future__ import annotations\n",
|
||
|
|
"\n",
|
||
|
|
"import math\n",
|
||
|
|
"import time\n",
|
||
|
|
"\n",
|
||
|
|
"import numpy as np\n",
|
||
|
|
"from ortools.sat.python import cp_model\n",
|
||
|
|
"\n",
|
||
|
|
"FENSTER = 20 # so viele Auftraege werden je Runde herausgebrochen\n",
|
||
|
|
"RUNDEN = 120\n",
|
||
|
|
"ZUEGE = 400_000 # Zugbudget des Annealings (wie Simulated_Annealing.py)\n",
|
||
|
|
"N = 500\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"# --- Instanz und Faustregel (wie in Simulated_Annealing.py) ----------------\n",
|
||
|
|
"\n",
|
||
|
|
"def erzeuge_ruestmatrix(n: int, saat: int = 11) -> np.ndarray:\n",
|
||
|
|
" rng = np.random.default_rng(saat)\n",
|
||
|
|
" familien = max(6, n // 10)\n",
|
||
|
|
" zwischen = rng.integers(8, 60, (familien, familien))\n",
|
||
|
|
" np.fill_diagonal(zwischen, 2)\n",
|
||
|
|
" familie = rng.integers(0, familien, n)\n",
|
||
|
|
" farbe = rng.integers(0, 10, n)\n",
|
||
|
|
" matrix = np.zeros((n, n), dtype=np.int64)\n",
|
||
|
|
" for i in range(n):\n",
|
||
|
|
" for j in range(n):\n",
|
||
|
|
" if i != j:\n",
|
||
|
|
" matrix[i, j] = (zwischen[familie[i], familie[j]]\n",
|
||
|
|
" + 3 * max(0, farbe[i] - farbe[j]))\n",
|
||
|
|
" return matrix\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def gesamtruestzeit(reihe: list[int], matrix: np.ndarray) -> int:\n",
|
||
|
|
" return int(sum(matrix[reihe[k], reihe[k + 1]] for k in range(len(reihe) - 1)))\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def faustregel(matrix: np.ndarray) -> list[int]:\n",
|
||
|
|
" n = len(matrix)\n",
|
||
|
|
" offen = set(range(1, n))\n",
|
||
|
|
" reihe = [0]\n",
|
||
|
|
" while offen:\n",
|
||
|
|
" naechster = min(offen, key=lambda j: matrix[reihe[-1], j])\n",
|
||
|
|
" reihe.append(naechster)\n",
|
||
|
|
" offen.discard(naechster)\n",
|
||
|
|
" return reihe\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def delta_verschieben(reihe, matrix, von, nach):\n",
|
||
|
|
" n = len(reihe)\n",
|
||
|
|
" if von == nach or nach == von + 1:\n",
|
||
|
|
" return None\n",
|
||
|
|
" heraus = 0\n",
|
||
|
|
" if von > 0:\n",
|
||
|
|
" heraus += matrix[reihe[von - 1], reihe[von]]\n",
|
||
|
|
" if von < n - 1:\n",
|
||
|
|
" heraus += matrix[reihe[von], reihe[von + 1]]\n",
|
||
|
|
" if 0 < von < n - 1:\n",
|
||
|
|
" heraus -= matrix[reihe[von - 1], reihe[von + 1]]\n",
|
||
|
|
" rest = reihe[:von] + reihe[von + 1:]\n",
|
||
|
|
" stelle = nach if nach < von else nach - 1\n",
|
||
|
|
" hinein = 0\n",
|
||
|
|
" if stelle > 0:\n",
|
||
|
|
" hinein += matrix[rest[stelle - 1], reihe[von]]\n",
|
||
|
|
" if stelle < len(rest):\n",
|
||
|
|
" hinein += matrix[reihe[von], rest[stelle]]\n",
|
||
|
|
" if 0 < stelle < len(rest):\n",
|
||
|
|
" hinein -= matrix[rest[stelle - 1], rest[stelle]]\n",
|
||
|
|
" return int(hinein - heraus), rest, stelle\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def annealing(matrix: np.ndarray, zuege: int = ZUEGE,\n",
|
||
|
|
" start_temperatur: float = 1.0, saat: int = 1) -> list[int]:\n",
|
||
|
|
" rng = np.random.default_rng(saat)\n",
|
||
|
|
" n = len(matrix)\n",
|
||
|
|
" reihe = faustregel(matrix)\n",
|
||
|
|
" kosten = gesamtruestzeit(reihe, matrix)\n",
|
||
|
|
" beste, beste_kosten = reihe[:], kosten\n",
|
||
|
|
" for zug in range(zuege):\n",
|
||
|
|
" temperatur = start_temperatur * (0.05 / start_temperatur) ** (zug / zuege)\n",
|
||
|
|
" von, nach = int(rng.integers(0, n)), int(rng.integers(0, n + 1))\n",
|
||
|
|
" ergebnis = delta_verschieben(reihe, matrix, von, nach)\n",
|
||
|
|
" if ergebnis is None:\n",
|
||
|
|
" continue\n",
|
||
|
|
" aenderung, rest, stelle = ergebnis\n",
|
||
|
|
" if aenderung <= 0 or rng.random() < math.exp(-aenderung / temperatur):\n",
|
||
|
|
" reihe = rest[:stelle] + [reihe[von]] + rest[stelle:]\n",
|
||
|
|
" kosten += aenderung\n",
|
||
|
|
" if kosten < beste_kosten:\n",
|
||
|
|
" beste, beste_kosten = reihe[:], kosten\n",
|
||
|
|
" return beste\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"# --- Der Reparaturschritt: ein kleines Problem, exakt geloest --------------\n",
|
||
|
|
"\n",
|
||
|
|
"def repariere_fenster(reihe: list[int], matrix: np.ndarray,\n",
|
||
|
|
" a: int, b: int) -> list[int] | None:\n",
|
||
|
|
" \"\"\"Ordnet reihe[a:b] optimal neu; die beiden Raender bleiben, wo sie sind.\n",
|
||
|
|
"\n",
|
||
|
|
" Die Raender festzuhalten ist entscheidend. Ohne sie waere das Teilproblem\n",
|
||
|
|
" ein anderes als der Ausschnitt aus dem Gesamtplan: Der Uebergang vom\n",
|
||
|
|
" Vorgaenger in das Fenster und aus dem Fenster in den Nachfolger gehoert\n",
|
||
|
|
" zu den Kosten dazu. Genau dafuer steht der Hilfsknoten k - er vertritt\n",
|
||
|
|
" beide Raender in einem.\n",
|
||
|
|
" \"\"\"\n",
|
||
|
|
" innen = reihe[a:b]\n",
|
||
|
|
" k = len(innen)\n",
|
||
|
|
" vorgaenger = reihe[a - 1] if a > 0 else None\n",
|
||
|
|
" nachfolger_rand = reihe[b] if b < len(reihe) else None\n",
|
||
|
|
"\n",
|
||
|
|
" modell = cp_model.CpModel()\n",
|
||
|
|
" kanten, ziel = [], []\n",
|
||
|
|
" for i in range(k + 1):\n",
|
||
|
|
" for j in range(k + 1):\n",
|
||
|
|
" if i == j:\n",
|
||
|
|
" continue\n",
|
||
|
|
" aktiv = modell.NewBoolVar(f\"kante_{i}_{j}\")\n",
|
||
|
|
" kanten.append((i, j, aktiv))\n",
|
||
|
|
" if i < k and j < k:\n",
|
||
|
|
" ziel.append(int(matrix[innen[i], innen[j]]) * aktiv)\n",
|
||
|
|
" elif i == k and j < k and vorgaenger is not None:\n",
|
||
|
|
" ziel.append(int(matrix[vorgaenger, innen[j]]) * aktiv)\n",
|
||
|
|
" elif j == k and i < k and nachfolger_rand is not None:\n",
|
||
|
|
" ziel.append(int(matrix[innen[i], nachfolger_rand]) * aktiv)\n",
|
||
|
|
" modell.AddCircuit(kanten)\n",
|
||
|
|
" modell.Minimize(sum(ziel))\n",
|
||
|
|
"\n",
|
||
|
|
" loeser = cp_model.CpSolver()\n",
|
||
|
|
" loeser.parameters.num_workers = 1\n",
|
||
|
|
" loeser.parameters.random_seed = 1\n",
|
||
|
|
" loeser.parameters.max_time_in_seconds = 10.0\n",
|
||
|
|
" status = loeser.Solve(modell)\n",
|
||
|
|
" if status not in (cp_model.OPTIMAL, cp_model.FEASIBLE):\n",
|
||
|
|
" # Auch das gehoert dazu: Wenn das Teilproblem nicht loest, bleibt der\n",
|
||
|
|
" # Plan, wie er war. Eine LNS-Runde darf scheitern.\n",
|
||
|
|
" return None\n",
|
||
|
|
"\n",
|
||
|
|
" nachfolger = {i: j for i, j, aktiv in kanten if loeser.Value(aktiv)}\n",
|
||
|
|
" neue_folge, aktuell = [], nachfolger[k]\n",
|
||
|
|
" while aktuell != k:\n",
|
||
|
|
" neue_folge.append(innen[aktuell])\n",
|
||
|
|
" aktuell = nachfolger[aktuell]\n",
|
||
|
|
" return reihe[:a] + neue_folge + reihe[b:]\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def lns(reihe: list[int], matrix: np.ndarray, runden: int = RUNDEN,\n",
|
||
|
|
" fenster: int = FENSTER, saat: int = 3):\n",
|
||
|
|
" \"\"\"Zerstoeren, exakt reparieren, uebernehmen - so lange das Budget reicht.\"\"\"\n",
|
||
|
|
" rng = np.random.default_rng(saat)\n",
|
||
|
|
" kosten = gesamtruestzeit(reihe, matrix)\n",
|
||
|
|
" verbesserungen = 0\n",
|
||
|
|
" start = time.perf_counter()\n",
|
||
|
|
" for _ in range(runden):\n",
|
||
|
|
" a = int(rng.integers(0, len(reihe) - fenster))\n",
|
||
|
|
" neu = repariere_fenster(reihe, matrix, a, a + fenster)\n",
|
||
|
|
" if neu is None:\n",
|
||
|
|
" continue\n",
|
||
|
|
" neue_kosten = gesamtruestzeit(neu, matrix)\n",
|
||
|
|
" if neue_kosten < kosten:\n",
|
||
|
|
" reihe, kosten = neu, neue_kosten\n",
|
||
|
|
" verbesserungen += 1\n",
|
||
|
|
" return reihe, kosten, verbesserungen, time.perf_counter() - start\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"if __name__ == \"__main__\":\n",
|
||
|
|
" matrix = erzeuge_ruestmatrix(N)\n",
|
||
|
|
" start_reihe = faustregel(matrix)\n",
|
||
|
|
" start_kosten = gesamtruestzeit(start_reihe, matrix)\n",
|
||
|
|
" SCHRANKE = 1768 # aus dem CP-SAT-Lauf in Metaheuristik_vs_Exakt.py\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"=\" * 84)\n",
|
||
|
|
" print(\" ZERSTOEREN UND EXAKT REPARIEREN\")\n",
|
||
|
|
" print(\"=\" * 84)\n",
|
||
|
|
" print(f\"{N} Auftraege. Je Runde werden {FENSTER} aufeinanderfolgende Auftraege\")\n",
|
||
|
|
" print(f\"herausgebrochen und exakt neu geordnet, {RUNDEN} Runden lang.\\n\")\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"Erst die Metaheuristik allein, dann LNS auf ihrem Ergebnis:\\n\")\n",
|
||
|
|
" sa_reihe = annealing(matrix)\n",
|
||
|
|
" sa_kosten = gesamtruestzeit(sa_reihe, matrix)\n",
|
||
|
|
"\n",
|
||
|
|
" lns_ab_faustregel = lns(start_reihe[:], matrix)\n",
|
||
|
|
" lns_ab_annealing = lns(sa_reihe[:], matrix)\n",
|
||
|
|
"\n",
|
||
|
|
" print(f\" {'Verfahren':<34} {'Ruestzeit':>10} {'ueber Schranke':>15}\")\n",
|
||
|
|
" print(\" \" + \"-\" * 62)\n",
|
||
|
|
" for name, wert in [\n",
|
||
|
|
" (\"Faustregel (ohne Solver)\", start_kosten),\n",
|
||
|
|
" (\"CP-SAT allein, 30 s\", 2628),\n",
|
||
|
|
" (\"Simulated Annealing\", sa_kosten),\n",
|
||
|
|
" (\"LNS ab Faustregel\", lns_ab_faustregel[1]),\n",
|
||
|
|
" (\"Annealing, dann LNS\", lns_ab_annealing[1])]:\n",
|
||
|
|
" print(f\" {name:<34} {wert:>10,} {(wert - SCHRANKE) / SCHRANKE * 100:>13.1f} %\")\n",
|
||
|
|
" print(f\" {'untere Schranke (CP-SAT)':<34} {SCHRANKE:>10,} {0.0:>13.1f} %\")\n",
|
||
|
|
"\n",
|
||
|
|
" print(f\"\\n LNS ab Faustregel: {lns_ab_faustregel[2]} von {RUNDEN} Runden brachten \"\n",
|
||
|
|
" f\"eine Verbesserung ({lns_ab_faustregel[3]:.0f} s).\")\n",
|
||
|
|
" print(f\" LNS ab Annealing: {lns_ab_annealing[2]} von {RUNDEN} Runden \"\n",
|
||
|
|
" f\"({lns_ab_annealing[3]:.0f} s).\")\n",
|
||
|
|
"\n",
|
||
|
|
" # --- Die Fenstergroesse ----------------------------------------------\n",
|
||
|
|
" # Verglichen wird bei ANNAEHERND GLEICHER ZEIT, nicht bei gleicher\n",
|
||
|
|
" # Rundenzahl: Eine Runde mit Fenster 30 kostet ein Vielfaches einer Runde\n",
|
||
|
|
" # mit Fenster 10. Wer nach Runden vergleicht, misst nur, dass groessere\n",
|
||
|
|
" # Teilprobleme mehr finden - und uebersieht, was sie kosten. Die\n",
|
||
|
|
" # Rundenzahlen unten sind so gewaehlt, dass alle drei Laeufe in derselben\n",
|
||
|
|
" # Groessenordnung liegen.\n",
|
||
|
|
" print(\"\\n\" + \"-\" * 84)\n",
|
||
|
|
" print(\"Die Fenstergroesse ist die eine Stellschraube - und sie hat ein\")\n",
|
||
|
|
" print(\"Optimum in der Mitte. Gleiche Zeit, verschiedene Fenster:\\n\")\n",
|
||
|
|
" print(f\" {'Fenster':>8} {'Runden':>8} {'Ruestzeit':>11} {'Verbesser-':>12} {'Zeit':>8}\")\n",
|
||
|
|
" print(f\" {'':>8} {'':>8} {'':>11} {'ungen':>12} {'':>8}\")\n",
|
||
|
|
" print(\" \" + \"-\" * 52)\n",
|
||
|
|
" fenster_ergebnis = {}\n",
|
||
|
|
" for fenster, runden in ((10, 4000), (20, 150), (30, 20)):\n",
|
||
|
|
" _, wert, verbesserungen, dauer = lns(sa_reihe[:], matrix,\n",
|
||
|
|
" runden=runden, fenster=fenster)\n",
|
||
|
|
" fenster_ergebnis[fenster] = wert\n",
|
||
|
|
" print(f\" {fenster:>8} {runden:>8,} {wert:>11,} {verbesserungen:>12} \"\n",
|
||
|
|
" f\"{dauer:>7.0f}s\")\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"\\n Nach unten ist die Grenze nicht die Zeit, sondern die SAETTIGUNG:\")\n",
|
||
|
|
" print(f\" Mit Fenster 10 bleiben von {4000:,} Runden nur eine Handvoll\")\n",
|
||
|
|
" print(\" Verbesserungen uebrig. Zwanzig Auftraege lassen sich in einem Zug\")\n",
|
||
|
|
" print(\" umbauen, zehn nicht - und was ein Fenster nicht umbauen kann, findet\")\n",
|
||
|
|
" print(\" es auch in beliebig vielen Runden nicht.\")\n",
|
||
|
|
" print(\"\\n Nach oben ist die Grenze die Rechenzeit: Fenster 30 findet pro Runde\")\n",
|
||
|
|
" print(\" mehr, kommt in derselben Zeit aber nur auf einen Bruchteil der Runden\")\n",
|
||
|
|
" print(f\" und landet bei {fenster_ergebnis[30]:,} statt {fenster_ergebnis[20]:,}. Eine Reparatur mit 30 Auftraegen\")\n",
|
||
|
|
" print(\" ist eben genau das Problem, an dem der exakte Solver im Grossen\")\n",
|
||
|
|
" print(\" scheitert - nur eine Nummer kleiner.\")\n",
|
||
|
|
" print(\"\\n Diese Tabelle gehoert an den Anfang jedes LNS-Projekts. Sie zu raten\")\n",
|
||
|
|
" print(\" statt zu messen ist der haeufigste Fehler beim Einsatz des Verfahrens.\")\n",
|
||
|
|
"\n",
|
||
|
|
" print(\"\\n\" + \"=\" * 84)\n",
|
||
|
|
" print(\" WAS DAS HEISST\")\n",
|
||
|
|
" print(\"=\" * 84)\n",
|
||
|
|
" bester = lns_ab_annealing[1]\n",
|
||
|
|
" print(f\"Der beste Plan kommt aus der Kombination: {bester:,} Minuten,\")\n",
|
||
|
|
" print(f\"{(start_kosten - bester) / start_kosten * 100:.1f} % unter der Faustregel und \"\n",
|
||
|
|
" f\"{(sa_kosten - bester) / sa_kosten * 100:.1f} % unter dem reinen Annealing.\")\n",
|
||
|
|
" print()\n",
|
||
|
|
" print(\"Weder der exakte Solver noch die Metaheuristik allein kommen dorthin.\")\n",
|
||
|
|
" print(\"Der Solver scheitert an der Groesse, die Metaheuristik an der Kleinheit\")\n",
|
||
|
|
" print(\"ihrer Zuege. LNS gibt dem Solver Teilprobleme in einer Groesse, die er\")\n",
|
||
|
|
" print(\"beherrscht, und der Metaheuristik die Umbauten, die sie nicht kann.\")\n",
|
||
|
|
" print()\n",
|
||
|
|
" print(f\"Und die Ehrlichkeit zum Schluss: Auch {bester:,} liegt noch \"\n",
|
||
|
|
" f\"{(bester - SCHRANKE) / SCHRANKE * 100:.0f} % ueber der\")\n",
|
||
|
|
" print(\"Schranke. Ob dort wirklich noch so viel Luft ist oder ob die Schranke\")\n",
|
||
|
|
" print(\"nur schwach ist, weiss man nicht - das ist die Lage, in der man mit\")\n",
|
||
|
|
" print(\"einer Metaheuristik arbeitet. Man kennt seine Loesung, nicht ihren\")\n",
|
||
|
|
" print(\"Abstand zum Optimum.\")\n",
|
||
|
|
" print(\"=\" * 84)"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"metadata": {
|
||
|
|
"kernelspec": {
|
||
|
|
"display_name": "Python 3",
|
||
|
|
"language": "python",
|
||
|
|
"name": "python3"
|
||
|
|
},
|
||
|
|
"language_info": {
|
||
|
|
"name": "python",
|
||
|
|
"version": "3.11"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"nbformat": 4,
|
||
|
|
"nbformat_minor": 5
|
||
|
|
}
|