operations_research/OR_HTML_04/programme/Simulated_Annealing.py

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