Kapitel 9: Metaheuristiken — wenn der exakte Solver aussteigt

#!/usr/bin/env python3

# Simulated_Annealing.py
"""
Kapitel Metaheuristiken: Simulated Annealing an der Lackieranlage.

Aufgabe: n Auftraege in eine Reihenfolge bringen. Zwischen zwei Auftraegen
faellt eine Ruestzeit an - die Anlage muss gereinigt werden. Wie lange das
dauert, haengt von BEIDEN Auftraegen ab: Ein Wechsel innerhalb derselben
Produktfamilie kostet fast nichts, ein Wechsel von Dunkel nach Hell ist teuer,
und zwischen manchen Familien ist der Reinigungsaufwand schlicht hoch. Gesucht
ist die Reihenfolge mit der kleinsten Summe der Ruestzeiten.

Das Programm zeigt drei Dinge:

  1. Warum die Faustregel "immer der billigste naechste Auftrag" in eine
     Sackgasse laeuft - und wie weit man sie mit lokaler Suche verbessert.
  2. Wie gross der Beitrag der Verschlechterungen WIRKLICH ist. Die Antwort
     faellt bescheidener aus als das Lehrbuch verspricht, und genau das ist
     der Grund, den Vergleich immer mitzurechnen.
  3. Dass die Starttemperatur kein Beiwerk ist: zu heiss macht die Suche
     nutzlos, und man sieht es an einer einzigen Kennzahl kommen.

Zum Budget: Der Abkuehlplan laeuft ueber eine feste ZUGZAHL, nicht ueber eine
Sekundenzahl. Im Betrieb hat man zwar ein Zeitbudget - fuer einen VERGLEICH ist
die Zugzahl aber das ehrlichere Mass: Sie ist auf jeder Maschine dieselbe, und
die abgedruckten Zahlen unten lassen sich damit nachrechnen. Die gemessene
Laufzeit steht trotzdem dabei.

Benoetigt: numpy
"""

from __future__ import annotations

import math
import time

import numpy as np

ZUEGE = 400_000           # Zugbudget je Lauf (statt Sekunden: reproduzierbar)
SAAT = 11                 # feste Saat -> reproduzierbare Instanz


# --- 1. Die Instanz ---------------------------------------------------------

def erzeuge_ruestmatrix(n: int, saat: int = SAAT) -> np.ndarray:
    """Ruestzeiten in Minuten zwischen je zwei Auftraegen.

    Die Zahl der Produktfamilien waechst mit der Auftragszahl (eine Familie je
    zehn Auftraege). Sonst wuerde das Problem mit wachsendem n LEICHTER: Bei
    fester Familienzahl haette der Plan irgendwann nur noch ein paar grosse
    Bloecke, und jede Faustregel faende ihn.
    """
    rng = np.random.default_rng(saat)
    familien = max(6, n // 10)
    # Reinigungsaufwand zwischen den Familien - unregelmaessig und asymmetrisch,
    # so wie im Betrieb: Von Klarlack auf Rot ist etwas anderes als umgekehrt.
    zwischen = rng.integers(8, 60, (familien, familien))
    np.fill_diagonal(zwischen, 2)

    familie = rng.integers(0, familien, n)
    farbe = rng.integers(0, 10, n)          # 0 = weiss ... 9 = schwarz

    matrix = np.zeros((n, n), dtype=np.int64)
    for i in range(n):
        for j in range(n):
            if i != j:
                # Dunkel -> hell kostet zusaetzlich: die Anlage muss heller
                # werden, das braucht mehr Spuelgaenge.
                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:
    """Summe der Ruestzeiten einer Reihenfolge (offene Kette, keine Rundreise)."""
    return int(sum(matrix[reihe[k], reihe[k + 1]] for k in range(len(reihe) - 1)))


def faustregel(matrix: np.ndarray) -> list[int]:
    """Immer der billigste noch offene Auftrag - 'naechster Nachbar'.

    So plant ein Meister von Hand, und es ist keine schlechte Regel. Ihr Fehler
    ist die Kurzsichtigkeit: Sie spart am Anfang und laesst die teuren Wechsel
    fuer das Ende uebrig, wo keine Wahl mehr bleibt.
    """
    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


# --- 2. Der Zug und seine Kostenaenderung -----------------------------------

def delta_verschieben(reihe: list[int], matrix: np.ndarray,
                      von: int, nach: int):
    """Auftrag von Position 'von' an Position 'nach' versetzen.

    Liefert (Kostenaenderung, Restliste, Einfuegeposition) oder None, wenn der
    Zug nichts aendert.

    DAS IST DIE WICHTIGSTE FUNKTION DES PROGRAMMS. Sie berechnet die
    Kostenaenderung aus HOECHSTENS SECHS Matrixeintraegen, statt die
    Reihenfolge neu durchzusummieren. Der Unterschied ist nicht kosmetisch:
    Die volle Summe kostet O(n) je Zug, diese Rechnung O(1). Bei 500 Auftraegen
    sind das rund 500-mal mehr geprueste Zuege im selben Zeitbudget - und eine
    Metaheuristik lebt von der Zahl der Zuege.

    Herausgerissen wird der Auftrag aus zwei Kanten, seine Nachbarn ruecken
    zusammen (eine neue Kante). Eingefuegt wird er zwischen zwei andere
    Nachbarn, deren bisherige Kante dadurch verschwindet.
    """
    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:                       # die Nachbarn ruecken zusammen
        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):                # die aufgetrennte Kante faellt weg
        hinein -= matrix[rest[stelle - 1], rest[stelle]]

    return int(hinein - heraus), rest, stelle


# --- 3. Die Suche -----------------------------------------------------------

class Verlauf:
    """Was ein Lauf ausser dem Ergebnis noch verraet."""

    def __init__(self) -> None:
        self.schritte = 0
        self.verschlechterungen_geprueft = 0
        self.verschlechterungen_genommen = 0
        self.schlechtester_zustand = 0      # wie weit die Suche abgedriftet ist
        self.sekunden = 0.0                 # nur zur Information, nicht zum Vergleich

    @property
    def annahmequote(self) -> float:
        if not self.verschlechterungen_geprueft:
            return 0.0
        return self.verschlechterungen_genommen / self.verschlechterungen_geprueft


def suche(matrix: np.ndarray, zuege: int, start_temperatur: float,
          end_temperatur: float = 0.05, saat: int = 1,
          bergsteigen: bool = False) -> tuple[list[int], int, Verlauf]:
    """Simulated Annealing (oder reines Bergsteigen, wenn bergsteigen=True).

    Der einzige Unterschied zwischen beiden steht in der Annahmezeile: Das
    Bergsteigen nimmt nur Verbesserungen, SA nimmt Verschlechterungen mit einer
    Wahrscheinlichkeit an, die mit der Zeit gegen null geht.
    """
    rng = np.random.default_rng(saat)
    n = len(matrix)
    reihe = faustregel(matrix)
    kosten = gesamtruestzeit(reihe, matrix)
    beste, beste_kosten = reihe[:], kosten
    verlauf = Verlauf()
    verlauf.schlechtester_zustand = kosten

    start = time.perf_counter()
    for zug in range(zuege):
        verlauf.schritte += 1
        # Geometrischer Abkuehlplan ueber das Zugbudget.
        temperatur = start_temperatur * (end_temperatur / start_temperatur) \
            ** (zug / zuege)

        von = int(rng.integers(0, n))
        nach = 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:
            annehmen = True
        else:
            verlauf.verschlechterungen_geprueft += 1
            # Metropolis-Kriterium: je groesser die Verschlechterung und je
            # kaelter es ist, desto unwahrscheinlicher.
            annehmen = (not bergsteigen
                        and rng.random() < math.exp(-aenderung / temperatur))
            verlauf.verschlechterungen_genommen += annehmen

        if annehmen:
            reihe = rest[:stelle] + [reihe[von]] + rest[stelle:]
            kosten += aenderung
            verlauf.schlechtester_zustand = max(verlauf.schlechtester_zustand,
                                                kosten)
            if kosten < beste_kosten:
                beste, beste_kosten = reihe[:], kosten

    verlauf.sekunden = time.perf_counter() - start
    return beste, beste_kosten, verlauf


if __name__ == "__main__":
    N = 200
    matrix = erzeuge_ruestmatrix(N)
    start_reihe = faustregel(matrix)
    start_kosten = gesamtruestzeit(start_reihe, matrix)

    print("=" * 78)
    print("  SIMULATED ANNEALING AN DER LACKIERANLAGE")
    print("=" * 78)
    print(f"{N} Auftraege, {max(6, N // 10)} Produktfamilien, "
          f"Zugbudget {ZUEGE:,} je Lauf.\n")
    print(f"Faustregel 'billigster naechster Auftrag': "
          f"{start_kosten:,} Minuten Ruestzeit")

    # --- Wie gross sind die Zuege ueberhaupt? ---------------------------
    # Ohne diese Zahl kann man die Temperatur nicht waehlen: Das
    # Metropolis-Kriterium vergleicht die Verschlechterung MIT der Temperatur.
    rng = np.random.default_rng(0)
    stichprobe = []
    for _ in range(5000):
        e = delta_verschieben(start_reihe, matrix,
                              int(rng.integers(0, N)), int(rng.integers(0, N + 1)))
        if e:
            stichprobe.append(abs(e[0]))
    print(f"Typische Zuggroesse |Delta|: Median {np.median(stichprobe):.0f}, "
          f"90 %-Quantil {np.quantile(stichprobe, 0.9):.0f} Minuten")

    # --- Bergsteigen gegen Annealing -------------------------------------
    print("\n" + "-" * 78)
    print("Nur Verbesserungen annehmen (Bergsteigen) gegen Annealing:\n")
    _, berg, berg_verlauf = suche(matrix, ZUEGE, 1.0, bergsteigen=True)
    _, gluehen, gluehen_verlauf = suche(matrix, ZUEGE, 1.0)
    print(f"  {'Faustregel (Start)':<34} {start_kosten:>7,} Minuten")
    print(f"  {'Bergsteigen':<34} {berg:>7,} Minuten  "
          f"({(start_kosten - berg) / start_kosten * 100:5.1f} % besser)")
    print(f"  {'Simulated Annealing':<34} {gluehen:>7,} Minuten  "
          f"({(start_kosten - gluehen) / start_kosten * 100:5.1f} % besser)")
    print(f"\n  Beide pruefen genau {berg_verlauf.schritte:,} Zuege "
          f"({berg_verlauf.sekunden:.1f} bzw. {gluehen_verlauf.sekunden:.1f} Sekunden).")
    print(f"  Annealing nimmt davon {gluehen_verlauf.annahmequote * 100:.2f} % der "
          f"VERSCHLECHTERUNGEN an, Bergsteigen keine.")
    print("\n  Bemerkenswert ist, wie klein der Abstand ist: Das simple")
    print(f"  Bergsteigen holt hier bereits {(start_kosten - berg) / (start_kosten - gluehen) * 100:.0f} % dessen, was Annealing")
    print("  schafft. Wer eine Metaheuristik einsetzt, sollte diesen Vergleich")
    print("  IMMER mitrechnen - sonst schreibt man einer aufwendigen Methode")
    print("  gut, was schon die einfache geliefert haette.")

    # --- Die Temperatur ist kein Beiwerk ---------------------------------
    print("\n" + "-" * 78)
    print("Und warum die Starttemperatur ausgemessen gehoert:\n")
    print(f"  {'Start-T':>8} {'Ruestzeit':>11} {'gg. Faustregel':>15} "
          f"{'angenommene':>12} {'schlechtester':>14}")
    print(f"  {'':>8} {'':>11} {'':>15} {'Verschlecht.':>12} {'Zwischenwert':>14}")
    print("  " + "-" * 64)
    ergebnisse, verlaeufe = {}, {}
    for temperatur in (0.5, 1.0, 2.0, 4.0, 8.0):
        _, wert, verlauf = suche(matrix, ZUEGE, temperatur)
        ergebnisse[temperatur] = wert
        verlaeufe[temperatur] = verlauf
        gewinn = (start_kosten - wert) / start_kosten * 100
        print(f"  {temperatur:>8.1f} {wert:>11,} {gewinn:>13.1f} % "
              f"{verlauf.annahmequote * 100:>11.2f} % "
              f"{verlauf.schlechtester_zustand:>14,}")

    beste_temperatur = min(ergebnisse, key=ergebnisse.get)
    heiss = verlaeufe[8.0]
    print(f"\n  Bester Wert bei T0 = {beste_temperatur}.")
    print("\n  Die letzten beiden Spalten erklaeren, warum es nach oben kippt.")
    print(f"  Bei T0 = 8 gehen nur {heiss.annahmequote * 100:.2f} % der Verschlechterungen")
    print("  durch - wenig genug, dass man es fuer harmlos halten koennte. Es")
    print(f"  genuegt aber: Die Suche driftet bis auf {heiss.schlechtester_zustand:,} Minuten ab, das")
    print(f"  {heiss.schlechtester_zustand / start_kosten:.1f}-fache der Startloesung, und findet im Zugbudget nicht")
    print("  mehr zurueck. Vom Ergebnis bleibt fast nichts uebrig: Die Suche hat")
    print(f"  in {ZUEGE:,} Zuegen nie etwas GESEHEN, das mehr als")
    print(f"  {(start_kosten - ergebnisse[8.0]) / start_kosten * 100:.1f} % unter der Startloesung lag.")
    print("\n  Bemerkenswert ist die andere Richtung: Die besten Werte entstehen")
    print("  bei Annahmequoten nahe null. Die Lehrbuchregel 'anfangs sollen")
    print("  20 bis 50 Prozent der Verschlechterungen durchgehen' fuehrt hier in")
    print("  die Irre - sie stammt aus Problemen mit sehr kleinen Zuggroessen.")
    print("  Bei Ruestzeiten in Minuten ist ein einziger schlechter Zug teuer.")
    print("\n  Was stattdessen traegt: Den SCHLECHTESTEN Zwischenwert protokollieren")
    print("  und T0 so waehlen, dass er die Startloesung nicht wesentlich")
    print("  ueberschreitet. Diese eine Zahl haette hier auf Anhieb gezeigt, dass")
    print("  T0 = 8 unbrauchbar ist - ohne dass man das Endergebnis abwarten muss.")
    print("=" * 78)

#!/usr/bin/env python3

# Metaheuristik_vs_Exakt.py
"""
Kapitel Metaheuristiken: Ab wann lohnt sich eine Metaheuristik? Die Antwort ist eine Zahl.

Die verbreitete Regel lautet: "Kleine Probleme exakt, grosse mit
Metaheuristiken." Sie ist richtig - aber voellig nutzlos, solange niemand sagt,
wo 'gross' anfaengt. Dieses Programm misst es fuer die Ruestzeitenaufgabe aus
Simulated_Annealing.py nach, und das Ergebnis ueberrascht in beide Richtungen:

  * Bei 60 Auftraegen loest CP-SAT BEWEISBAR optimal, in wenigen Sekunden.
    Die Metaheuristik ist hier schlicht die schlechtere Wahl - sie liefert ein
    schwaecheres Ergebnis und weiss nicht einmal, wie schwach.
  * Bei 200 Auftraegen liegt CP-SAT ohne Optimalitaetsbeweis vorne.
  * Erst bei 500 Auftraegen kippt es: Dort gibt CP-SAT im Zeitlimit eine
    Loesung aus, die SCHLECHTER ist als die Faustregel eines Meisters.

Der Punkt des Kapitels steht in der letzten Spalte: Die Schranke. Nur der
exakte Solver sagt, wie gut die Loesung sein KANN. Eine Metaheuristik allein
liefert eine Zahl ohne Massstab - man weiss nie, ob man 2 % oder 40 % neben dem
Optimum liegt.

ZUR REPRODUZIERBARKEIT: Die Spalte 'Annealing' ist auf jeder Maschine gleich -
sie laeuft ueber ein festes ZUGbudget. Die CP-SAT-Spalten sind es nicht: Ein
Zeitlimit ist ein Wanduhr-Limit, und wie weit der Solver in 30 Sekunden kommt,
haengt von der Maschine ab. CP-SAT kennt zwar ein deterministisches Zeitmass
(max_deterministic_time), das aber in einer Groessenordnung rechnet, die fuer
ein Buchbeispiel unbrauchbar langsam ist. Die Zahlen unten sind also
hardwareabhaengig - die AUSSAGE der Tabelle ist es nicht, sie ist auf
verschiedenen Groessen und in mehreren Laeufen stabil.

WICHTIG: Hier wird nur ortools importiert, nicht highspy (Kapitel Oekosystem).

Benoetigt: numpy, ortools
"""

from __future__ import annotations

import time

import numpy as np
from ortools.sat.python import cp_model

ZUEGE = 400_000           # Zugbudget der Metaheuristik (wie Simulated_Annealing.py)
ZEITLIMIT = 30.0          # Zeitbudget des exakten Solvers
GROESSEN = (60, 200, 500)


# --- Instanz, Faustregel und Suche: identisch zu 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) -> tuple[int, float]:
    """Dieselbe Suche wie in Simulated_Annealing.py, hier nur als Vergleichswert."""
    import math

    rng = np.random.default_rng(saat)
    n = len(matrix)
    reihe = faustregel(matrix)
    kosten = gesamtruestzeit(reihe, matrix)
    beste_kosten = kosten

    t0 = time.perf_counter()
    for zug in range(zuege):
        temperatur = start_temperatur * (0.05 / start_temperatur) ** (zug / zuege)
        von = int(rng.integers(0, n))
        nach = 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
            beste_kosten = min(beste_kosten, kosten)
    return beste_kosten, time.perf_counter() - t0


# --- Der exakte Weg: CP-SAT mit AddCircuit ---------------------------------

def exakt(matrix: np.ndarray, zeitlimit: float = ZEITLIMIT):
    """Reihenfolge als Rundreise mit einem kostenlosen Hilfsknoten.

    AddCircuit verlangt einen geschlossenen Kreis. Unsere Aufgabe ist aber eine
    offene Kette: Der erste Auftrag hat keinen Vorgaenger, der letzte keinen
    Nachfolger. Der uebliche Kniff ist ein zusaetzlicher Hilfsknoten n, dessen
    Kanten nichts kosten. Der Kreis laeuft dann Hilfsknoten -> Kette ->
    Hilfsknoten, und die Kosten des Kreises sind genau die der Kette.

    AddCircuit ist der Grund, warum CP-SAT hier so lange mithaelt: Die
    Nebenbedingung sorgt selbst dafuer, dass keine Teilkreise entstehen. Wer
    dasselbe als MILP formuliert, braucht dafuer entweder exponentiell viele
    Ungleichungen oder die schwache MTZ-Formulierung (Kapitel Graphen).
    """
    n = len(matrix)
    modell = cp_model.CpModel()
    kanten, ziel = [], []
    for i in range(n + 1):
        for j in range(n + 1):
            if i == j:
                continue
            aktiv = modell.NewBoolVar(f"kante_{i}_{j}")
            kanten.append((i, j, aktiv))
            if i < n and j < n:                  # Hilfsknoten kostet nichts
                ziel.append(int(matrix[i, j]) * aktiv)
    modell.AddCircuit(kanten)
    modell.Minimize(sum(ziel))

    loeser = cp_model.CpSolver()
    loeser.parameters.max_time_in_seconds = zeitlimit
    loeser.parameters.num_workers = 8
    loeser.parameters.random_seed = 1
    t0 = time.perf_counter()
    status = loeser.Solve(modell)
    dauer = time.perf_counter() - t0
    return (loeser.StatusName(status), int(loeser.ObjectiveValue()),
            int(loeser.BestObjectiveBound()), dauer)


if __name__ == "__main__":
    print("=" * 88)
    print("  AB WANN LOHNT SICH DIE METAHEURISTIK?")
    print("=" * 88)
    print(f"Dieselbe Aufgabe in drei Groessen. Metaheuristik: {ZUEGE:,} Zuege. "
          f"Exakt: {ZEITLIMIT:.0f} s Zeitlimit.\n")

    print(f"{'Auftraege':>10} {'Faustregel':>11} {'Annealing':>10} "
          f"{'CP-SAT':>9} {'Status':>10} {'Schranke':>9} {'CP-Zeit':>9}")
    print("-" * 88)

    zeilen = []
    for n in GROESSEN:
        matrix = erzeuge_ruestmatrix(n)
        start = gesamtruestzeit(faustregel(matrix), matrix)
        heuristisch, _ = annealing(matrix)
        status, wert, schranke, dauer = exakt(matrix)
        zeilen.append((n, start, heuristisch, wert, status, schranke, dauer))
        print(f"{n:>10} {start:>11,} {heuristisch:>10,} {wert:>9,} "
              f"{status:>10} {schranke:>9,} {dauer:>8.1f}s")

    print("-" * 88)
    print("\nWas in jeder Zeile steht:\n")
    for n, start, heur, wert, status, schranke, _ in zeilen:
        if status == "OPTIMAL":
            urteil = (f"CP-SAT beweist das Optimum. Annealing liegt "
                      f"{(heur - wert) / wert * 100:.1f} % daneben - und haette")
            zweite = "     das ohne den exakten Lauf nicht erfahren."
        elif wert <= heur:
            urteil = (f"CP-SAT liegt {(heur - wert) / heur * 100:.1f} % vor Annealing, "
                      f"ohne Beweis (Gap {(wert - schranke) / wert * 100:.1f} %).")
            zweite = "     Auch ohne Optimalitaetsbeweis ist es die bessere Wahl."
        else:
            urteil = (f"UMSCHLAGPUNKT: CP-SAT ist {(wert - heur) / heur * 100:.1f} % "
                      f"SCHLECHTER als Annealing")
            zweite = (f"     und sogar {(wert - start) / start * 100:.1f} % schlechter "
                      f"als die Faustregel des Meisters.")
        print(f"  {n:>4} Auftraege: {urteil}")
        print(zweite)

    n, start, heur, wert, status, schranke, _ = zeilen[-1]
    print("\n" + "=" * 88)
    print("  DIE SPALTE, DIE MAN NICHT WEGLASSEN DARF")
    print("=" * 88)
    print("Die Schranke ist der eigentliche Ertrag des exakten Solvers - auch dann,")
    print("wenn seine LOESUNG unbrauchbar ist.")
    print()
    print(f"Bei {n} Auftraegen liefert Annealing {heur:,} Minuten. Klingt gut: "
          f"{(start - heur) / start * 100:.1f} %")
    print("besser als die Faustregel. Die Schranke aus demselben CP-SAT-Lauf, dessen")
    print(f"Loesung wir gerade verworfen haben, sagt aber: Unter {schranke:,} Minuten geht")
    print(f"es nicht. Zwischen {schranke:,} und {heur:,} liegen {(heur - schranke) / heur * 100:.0f} % - so viel Luft")
    print("kann noch nach oben sein.")
    print()
    print("Das ist der Unterschied zwischen 'besser als vorher' und 'gut'. Eine")
    print("Metaheuristik allein kann nur das Erste. Deshalb gehoert auch dann ein")
    print("exakter Lauf dazu, wenn man am Ende die heuristische Loesung einsetzt:")
    print("nicht wegen seiner Loesung, sondern wegen seiner Schranke.")
    print("=" * 88)

#!/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)