#!/usr/bin/env python3 # Large_Neighborhood_Search.py """ Kapitel Metaheuristiken: Zerstoeren und exakt reparieren - der Solver im Dienst der Heuristik. Metaheuristik_vs_Exakt.py endet mit einem Widerspruch. Bei 500 Auftraegen gilt: * Der exakte Solver scheitert an der GROESSE - seine Loesung nach 30 Sekunden ist schlechter als die Faustregel eines Meisters. * Die Metaheuristik kommt weiter, laesst aber laut Schranke noch rund ein Viertel liegen. Ihre Zuege sind zu kleinteilig: Sie verschiebt jeweils EINEN Auftrag und kann eine ganze Passage nicht auf einmal umbauen. Large Neighborhood Search loest den Widerspruch, indem sie beide einsetzt - jeden fuer das, was er kann: ZERSTOEREN Ein Stueck des Plans herausbrechen (hier: ein zusammen- haengendes Fenster von 20 Auftraegen). REPARIEREN Genau dieses Stueck EXAKT neu optimieren. 20 Auftraege sind fuer CP-SAT eine Kleinigkeit, 500 sind es nicht. UEBERNEHMEN Nur behalten, wenn der Gesamtplan besser wurde. Der Solver bekommt also nicht mehr das ganze Problem, sondern immer wieder ein kleines. Das ist der ganze Trick, und er ist in der Praxis der wichtigste Baustein dieses Kapitels. ZUR REPRODUZIERBARKEIT: Alle Laeufe haben eine feste RUNDENzahl, nicht ein Zeitbudget - die Ergebniswerte sind damit auf jeder Maschine gleich. Die Laufzeiten daneben sind hardwareabhaengig. Benoetigt: numpy, ortools """ from __future__ import annotations import math import time import numpy as np from ortools.sat.python import cp_model FENSTER = 20 # so viele Auftraege werden je Runde herausgebrochen RUNDEN = 120 ZUEGE = 400_000 # Zugbudget des Annealings (wie Simulated_Annealing.py) N = 500 # --- Instanz und Faustregel (wie in Simulated_Annealing.py) ---------------- def erzeuge_ruestmatrix(n: int, saat: int = 11) -> np.ndarray: rng = np.random.default_rng(saat) familien = max(6, n // 10) zwischen = rng.integers(8, 60, (familien, familien)) np.fill_diagonal(zwischen, 2) familie = rng.integers(0, familien, n) farbe = rng.integers(0, 10, n) matrix = np.zeros((n, n), dtype=np.int64) for i in range(n): for j in range(n): if i != j: matrix[i, j] = (zwischen[familie[i], familie[j]] + 3 * max(0, farbe[i] - farbe[j])) return matrix def gesamtruestzeit(reihe: list[int], matrix: np.ndarray) -> int: return int(sum(matrix[reihe[k], reihe[k + 1]] for k in range(len(reihe) - 1))) def faustregel(matrix: np.ndarray) -> list[int]: n = len(matrix) offen = set(range(1, n)) reihe = [0] while offen: naechster = min(offen, key=lambda j: matrix[reihe[-1], j]) reihe.append(naechster) offen.discard(naechster) return reihe def delta_verschieben(reihe, matrix, von, nach): n = len(reihe) if von == nach or nach == von + 1: return None heraus = 0 if von > 0: heraus += matrix[reihe[von - 1], reihe[von]] if von < n - 1: heraus += matrix[reihe[von], reihe[von + 1]] if 0 < von < n - 1: heraus -= matrix[reihe[von - 1], reihe[von + 1]] rest = reihe[:von] + reihe[von + 1:] stelle = nach if nach < von else nach - 1 hinein = 0 if stelle > 0: hinein += matrix[rest[stelle - 1], reihe[von]] if stelle < len(rest): hinein += matrix[reihe[von], rest[stelle]] if 0 < stelle < len(rest): hinein -= matrix[rest[stelle - 1], rest[stelle]] return int(hinein - heraus), rest, stelle def annealing(matrix: np.ndarray, zuege: int = ZUEGE, start_temperatur: float = 1.0, saat: int = 1) -> list[int]: rng = np.random.default_rng(saat) n = len(matrix) reihe = faustregel(matrix) kosten = gesamtruestzeit(reihe, matrix) beste, beste_kosten = reihe[:], kosten for zug in range(zuege): temperatur = start_temperatur * (0.05 / start_temperatur) ** (zug / zuege) von, nach = int(rng.integers(0, n)), int(rng.integers(0, n + 1)) ergebnis = delta_verschieben(reihe, matrix, von, nach) if ergebnis is None: continue aenderung, rest, stelle = ergebnis if aenderung <= 0 or rng.random() < math.exp(-aenderung / temperatur): reihe = rest[:stelle] + [reihe[von]] + rest[stelle:] kosten += aenderung if kosten < beste_kosten: beste, beste_kosten = reihe[:], kosten return beste # --- Der Reparaturschritt: ein kleines Problem, exakt geloest -------------- def repariere_fenster(reihe: list[int], matrix: np.ndarray, a: int, b: int) -> list[int] | None: """Ordnet reihe[a:b] optimal neu; die beiden Raender bleiben, wo sie sind. Die Raender festzuhalten ist entscheidend. Ohne sie waere das Teilproblem ein anderes als der Ausschnitt aus dem Gesamtplan: Der Uebergang vom Vorgaenger in das Fenster und aus dem Fenster in den Nachfolger gehoert zu den Kosten dazu. Genau dafuer steht der Hilfsknoten k - er vertritt beide Raender in einem. """ innen = reihe[a:b] k = len(innen) vorgaenger = reihe[a - 1] if a > 0 else None nachfolger_rand = reihe[b] if b < len(reihe) else None modell = cp_model.CpModel() kanten, ziel = [], [] for i in range(k + 1): for j in range(k + 1): if i == j: continue aktiv = modell.NewBoolVar(f"kante_{i}_{j}") kanten.append((i, j, aktiv)) if i < k and j < k: ziel.append(int(matrix[innen[i], innen[j]]) * aktiv) elif i == k and j < k and vorgaenger is not None: ziel.append(int(matrix[vorgaenger, innen[j]]) * aktiv) elif j == k and i < k and nachfolger_rand is not None: ziel.append(int(matrix[innen[i], nachfolger_rand]) * aktiv) modell.AddCircuit(kanten) modell.Minimize(sum(ziel)) loeser = cp_model.CpSolver() loeser.parameters.num_workers = 1 loeser.parameters.random_seed = 1 loeser.parameters.max_time_in_seconds = 10.0 status = loeser.Solve(modell) if status not in (cp_model.OPTIMAL, cp_model.FEASIBLE): # Auch das gehoert dazu: Wenn das Teilproblem nicht loest, bleibt der # Plan, wie er war. Eine LNS-Runde darf scheitern. return None nachfolger = {i: j for i, j, aktiv in kanten if loeser.Value(aktiv)} neue_folge, aktuell = [], nachfolger[k] while aktuell != k: neue_folge.append(innen[aktuell]) aktuell = nachfolger[aktuell] return reihe[:a] + neue_folge + reihe[b:] def lns(reihe: list[int], matrix: np.ndarray, runden: int = RUNDEN, fenster: int = FENSTER, saat: int = 3): """Zerstoeren, exakt reparieren, uebernehmen - so lange das Budget reicht.""" rng = np.random.default_rng(saat) kosten = gesamtruestzeit(reihe, matrix) verbesserungen = 0 start = time.perf_counter() for _ in range(runden): a = int(rng.integers(0, len(reihe) - fenster)) neu = repariere_fenster(reihe, matrix, a, a + fenster) if neu is None: continue neue_kosten = gesamtruestzeit(neu, matrix) if neue_kosten < kosten: reihe, kosten = neu, neue_kosten verbesserungen += 1 return reihe, kosten, verbesserungen, time.perf_counter() - start if __name__ == "__main__": matrix = erzeuge_ruestmatrix(N) start_reihe = faustregel(matrix) start_kosten = gesamtruestzeit(start_reihe, matrix) SCHRANKE = 1768 # aus dem CP-SAT-Lauf in Metaheuristik_vs_Exakt.py print("=" * 84) print(" ZERSTOEREN UND EXAKT REPARIEREN") print("=" * 84) print(f"{N} Auftraege. Je Runde werden {FENSTER} aufeinanderfolgende Auftraege") print(f"herausgebrochen und exakt neu geordnet, {RUNDEN} Runden lang.\n") print("Erst die Metaheuristik allein, dann LNS auf ihrem Ergebnis:\n") sa_reihe = annealing(matrix) sa_kosten = gesamtruestzeit(sa_reihe, matrix) lns_ab_faustregel = lns(start_reihe[:], matrix) lns_ab_annealing = lns(sa_reihe[:], matrix) print(f" {'Verfahren':<34} {'Ruestzeit':>10} {'ueber Schranke':>15}") print(" " + "-" * 62) for name, wert in [ ("Faustregel (ohne Solver)", start_kosten), ("CP-SAT allein, 30 s", 2628), ("Simulated Annealing", sa_kosten), ("LNS ab Faustregel", lns_ab_faustregel[1]), ("Annealing, dann LNS", lns_ab_annealing[1])]: print(f" {name:<34} {wert:>10,} {(wert - SCHRANKE) / SCHRANKE * 100:>13.1f} %") print(f" {'untere Schranke (CP-SAT)':<34} {SCHRANKE:>10,} {0.0:>13.1f} %") print(f"\n LNS ab Faustregel: {lns_ab_faustregel[2]} von {RUNDEN} Runden brachten " f"eine Verbesserung ({lns_ab_faustregel[3]:.0f} s).") print(f" LNS ab Annealing: {lns_ab_annealing[2]} von {RUNDEN} Runden " f"({lns_ab_annealing[3]:.0f} s).") # --- Die Fenstergroesse ---------------------------------------------- # Verglichen wird bei ANNAEHERND GLEICHER ZEIT, nicht bei gleicher # Rundenzahl: Eine Runde mit Fenster 30 kostet ein Vielfaches einer Runde # mit Fenster 10. Wer nach Runden vergleicht, misst nur, dass groessere # Teilprobleme mehr finden - und uebersieht, was sie kosten. Die # Rundenzahlen unten sind so gewaehlt, dass alle drei Laeufe in derselben # Groessenordnung liegen. print("\n" + "-" * 84) print("Die Fenstergroesse ist die eine Stellschraube - und sie hat ein") print("Optimum in der Mitte. Gleiche Zeit, verschiedene Fenster:\n") print(f" {'Fenster':>8} {'Runden':>8} {'Ruestzeit':>11} {'Verbesser-':>12} {'Zeit':>8}") print(f" {'':>8} {'':>8} {'':>11} {'ungen':>12} {'':>8}") print(" " + "-" * 52) fenster_ergebnis = {} for fenster, runden in ((10, 4000), (20, 150), (30, 20)): _, wert, verbesserungen, dauer = lns(sa_reihe[:], matrix, runden=runden, fenster=fenster) fenster_ergebnis[fenster] = wert print(f" {fenster:>8} {runden:>8,} {wert:>11,} {verbesserungen:>12} " f"{dauer:>7.0f}s") print("\n Nach unten ist die Grenze nicht die Zeit, sondern die SAETTIGUNG:") print(f" Mit Fenster 10 bleiben von {4000:,} Runden nur eine Handvoll") print(" Verbesserungen uebrig. Zwanzig Auftraege lassen sich in einem Zug") print(" umbauen, zehn nicht - und was ein Fenster nicht umbauen kann, findet") print(" es auch in beliebig vielen Runden nicht.") print("\n Nach oben ist die Grenze die Rechenzeit: Fenster 30 findet pro Runde") print(" mehr, kommt in derselben Zeit aber nur auf einen Bruchteil der Runden") print(f" und landet bei {fenster_ergebnis[30]:,} statt {fenster_ergebnis[20]:,}. Eine Reparatur mit 30 Auftraegen") print(" ist eben genau das Problem, an dem der exakte Solver im Grossen") print(" scheitert - nur eine Nummer kleiner.") print("\n Diese Tabelle gehoert an den Anfang jedes LNS-Projekts. Sie zu raten") print(" statt zu messen ist der haeufigste Fehler beim Einsatz des Verfahrens.") print("\n" + "=" * 84) print(" WAS DAS HEISST") print("=" * 84) bester = lns_ab_annealing[1] print(f"Der beste Plan kommt aus der Kombination: {bester:,} Minuten,") print(f"{(start_kosten - bester) / start_kosten * 100:.1f} % unter der Faustregel und " f"{(sa_kosten - bester) / sa_kosten * 100:.1f} % unter dem reinen Annealing.") print() print("Weder der exakte Solver noch die Metaheuristik allein kommen dorthin.") print("Der Solver scheitert an der Groesse, die Metaheuristik an der Kleinheit") print("ihrer Zuege. LNS gibt dem Solver Teilprobleme in einer Groesse, die er") print("beherrscht, und der Metaheuristik die Umbauten, die sie nicht kann.") print() print(f"Und die Ehrlichkeit zum Schluss: Auch {bester:,} liegt noch " f"{(bester - SCHRANKE) / SCHRANKE * 100:.0f} % ueber der") print("Schranke. Ob dort wirklich noch so viel Luft ist oder ob die Schranke") print("nur schwach ist, weiss man nicht - das ist die Lage, in der man mit") print("einer Metaheuristik arbeitet. Man kennt seine Loesung, nicht ihren") print("Abstand zum Optimum.") print("=" * 84)