operations_research/OR_HTML_04/programme/Warmstart_Effekt.py

#!/usr/bin/env python3

# Warmstart_Effekt.py
"""
Kapitel MILP: Wann ein Warm-Start wirklich etwas bringt.

Solverstatus_und_Gap.py zeigt einen Fall, in dem ein Startwert NICHTS
bringt - HiGHS findet von allein schneller etwas Besseres. Hier der
Gegenfall: eine Heuristik, die dem Solver echtes Domaenenwissen liefert.

Problem: Lastverteilung. n Auftraege mit bekannter Dauer sind auf m
gleichartige Maschinen zu verteilen, sodass die zuletzt fertige Maschine
so frueh wie moeglich fertig wird (Makespan-Minimierung).

Die Heuristik: LPT (Longest Processing Time first) - laengste Auftraege
zuerst, jeder auf die momentan am wenigsten belastete Maschine. Sie ist
Jahrzehnte alt, in zwei Zeilen geschrieben und beweisbar nie schlechter
als 4/3 des Optimums.

Ueber model.AddHint() bekommt CP-SAT diese Loesung als Startpunkt.

WICHTIG: highspy wird hier bewusst NICHT importiert - es vertraegt sich
nicht mit ortools im selben Prozess (siehe Kapitel Oekosystem).

Benoetigt: numpy, ortools
"""

from __future__ import annotations

import time

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

RNG = np.random.default_rng(4)


def lpt_heuristik(dauer: np.ndarray, n_maschinen: int) -> tuple[np.ndarray, int]:
    """Longest Processing Time first.

    Laengste Auftraege zuerst auf die jeweils freieste Maschine legen. Zwei
    Zeilen, keine Bibliothek, Ergebnis in Mikrosekunden - und erstaunlich
    nah am Optimum.
    """
    zuordnung = np.zeros(len(dauer), dtype=int)
    belegung = np.zeros(n_maschinen)
    for auftrag in np.argsort(-dauer):            # laengster zuerst
        maschine = int(np.argmin(belegung))       # freieste Maschine
        zuordnung[auftrag] = maschine
        belegung[maschine] += dauer[auftrag]
    return zuordnung, int(belegung.max())


def loese(dauer: np.ndarray, n_maschinen: int, zeitlimit: float,
          hinweis: np.ndarray | None = None) -> tuple[str, int, float]:
    """Exaktes Modell mit CP-SAT, optional mit Startloesung als Hinweis."""
    n_auftraege = len(dauer)
    obergrenze = int(dauer.sum())

    modell = cp_model.CpModel()
    # x[i][k] = 1  <=>  Auftrag i laeuft auf Maschine k
    x = [[modell.NewBoolVar(f"x_{i}_{k}") for k in range(n_maschinen)]
         for i in range(n_auftraege)]
    for i in range(n_auftraege):
        modell.AddExactlyOne(x[i])                # jeder Auftrag genau einmal

    belegung = [modell.NewIntVar(0, obergrenze, f"last_{k}")
                for k in range(n_maschinen)]
    for k in range(n_maschinen):
        modell.Add(belegung[k] == sum(int(dauer[i]) * x[i][k]
                                      for i in range(n_auftraege)))

    makespan = modell.NewIntVar(0, obergrenze, "makespan")
    modell.AddMaxEquality(makespan, belegung)     # das Maximum ueber alle Maschinen
    modell.Minimize(makespan)

    # Der Warm-Start: ein Hinweis pro Variable. CP-SAT muss ihn nicht
    # befolgen - er nutzt ihn als erste Loesung, wenn er zulaessig ist.
    if hinweis is not None:
        for i in range(n_auftraege):
            for k in range(n_maschinen):
                modell.AddHint(x[i][k], 1 if hinweis[i] == k else 0)

    loeser = cp_model.CpSolver()
    loeser.parameters.max_time_in_seconds = zeitlimit
    # Ein Arbeiter und fester Startwert, damit die Messung reproduzierbar ist -
    # der Seed allein genuegt dafuer NICHT (Kapitel Constraint Programming).
    # Im Produktivbetrieb laesst man beides auf den Standardwerten.
    loeser.parameters.num_workers = 1
    loeser.parameters.random_seed = 1

    t0 = time.perf_counter()
    status = loeser.Solve(modell)
    dauer_s = time.perf_counter() - t0

    if status not in (cp_model.OPTIMAL, cp_model.FEASIBLE):
        raise RuntimeError(f"Kein Plan gefunden: {loeser.StatusName(status)}")
    return loeser.StatusName(status), int(loeser.ObjectiveValue()), dauer_s


if __name__ == "__main__":
    print("=" * 78)
    print("  WARM-START: WENN DIE HEURISTIK MEHR WEISS ALS DER SOLVER")
    print("=" * 78)
    print("Lastverteilung: Auftraege auf gleichartige Maschinen verteilen,")
    print("sodass die letzte Maschine so frueh wie moeglich fertig wird.\n")

    print(f"{'Instanz':<22} {'Variante':<22} {'Makespan':>9} {'Zeit':>9} "
          f"{'Faktor':>8}")
    print("-" * 78)

    for n_auftraege, n_maschinen in [(60, 7), (80, 9)]:
        dauer = RNG.integers(10, 90, n_auftraege)
        start, lpt_wert = lpt_heuristik(dauer, n_maschinen)
        untere_schranke = dauer.sum() / n_maschinen

        instanz = f"{n_auftraege} Auftr., {n_maschinen} Masch."
        print(f"{instanz:<22} {'LPT-Heuristik':<22} {lpt_wert:>9} "
              f"{'< 0.001s':>9} {'':>8}")

        _, ziel_kalt, zeit_kalt = loese(dauer, n_maschinen, 60.0)
        print(f"{'':<22} {'CP-SAT kalt':<22} {ziel_kalt:>9} "
              f"{zeit_kalt:>8.3f}s {'1,0x':>8}")

        _, ziel_warm, zeit_warm = loese(dauer, n_maschinen, 60.0, hinweis=start)
        print(f"{'':<22} {'CP-SAT + LPT-Hinweis':<22} {ziel_warm:>9} "
              f"{zeit_warm:>8.3f}s {zeit_kalt / zeit_warm:>7.1f}x")

        assert ziel_kalt == ziel_warm, \
            "Der Hinweis darf das Optimum nicht veraendern - nur den Weg dorthin!"
        print(f"{'':<22} {'untere Schranke':<22} {untere_schranke:>9.1f}")
        print("-" * 78)

    print("\nDrei Beobachtungen:")
    print("1. Der Hinweis aendert das ERGEBNIS nicht - beide Laeufe finden")
    print("   dasselbe Optimum. Er aendert nur, wie lange der Beweis dauert.")
    print("   Genau deshalb ist ein Warm-Start ungefaehrlich: Ein schlechter")
    print("   Hinweis kostet Zeit, er verfaelscht aber nie die Loesung.")
    print("2. Die LPT-Heuristik liegt schon sehr nah am Optimum. Ihr Wert fuer")
    print("   den Solver liegt weniger in der Qualitaet als darin, dass sie")
    print("   SOFORT da ist - der Solver kann von Beginn an alles verwerfen,")
    print("   was schlechter ist.")
    print("3. Der Faktor schwankt von Instanz zu Instanz - oben 2,2x und 1,4x -")
    print("   und laesst sich NICHT aus der Problemgroesse ableiten. Er haengt")
    print("   davon ab, wie schnell der Solver von allein eine vergleichbar gute")
    print("   Loesung findet. Das ist die eigentliche Lehre: Ein Warm-Start ist")
    print("   eine Messung wert, keine Glaubensfrage.")
    print("=" * 78)