Kapitel 7: Constraint Programming mit CP-SAT — Logik, Scheduling und Zuweisung

#!/usr/bin/env python3

# Propagation_Demo.py
"""
Kapitel CP-SAT: Propagation sichtbar machen.
Vergleicht die Zahl der geprueften Kombinationen bei roher Aufzaehlung
mit der Zahl der Verzweigungen, die CP-SAT tatsaechlich braucht.
"""

import itertools

from ortools.sat.python import cp_model


def rohe_aufzaehlung(n: int, obergrenze: int) -> tuple[int, int]:
    """Zaehlt ALLE Kombinationen und prueft jede einzeln (Brute Force)."""
    geprueft = 0
    loesungen = 0
    for kombination in itertools.product(range(1, obergrenze + 1), repeat=n):
        geprueft += 1
        if len(set(kombination)) == n and sum(kombination) == 3 * n:
            loesungen += 1
    return geprueft, loesungen


def mit_cp_sat(n: int, obergrenze: int) -> tuple[int, int, float]:
    """Dasselbe Problem deklarativ: alle verschieden, Summe = 3n."""
    modell = cp_model.CpModel()
    x = [modell.NewIntVar(1, obergrenze, f"x{i}") for i in range(n)]
    modell.AddAllDifferent(x)                 # globales Constraint
    modell.Add(sum(x) == 3 * n)

    loeser = cp_model.CpSolver()
    loeser.parameters.enumerate_all_solutions = True

    class Zaehler(cp_model.CpSolverSolutionCallback):
        def __init__(self):
            super().__init__()
            self.anzahl = 0

        def on_solution_callback(self):
            self.anzahl += 1

    zaehler = Zaehler()
    loeser.Solve(modell, zaehler)
    return zaehler.anzahl, loeser.NumBranches(), loeser.WallTime()


if __name__ == "__main__":
    print("=" * 84)
    print("  PROPAGATION: WIE VIEL ARBEIT SPART SICH DER SOLVER?")
    print("  Aufgabe: n Zahlen aus 1..G, alle verschieden, Summe = 3n")
    print("=" * 84)
    print(f"{'n':>3} {'G':>4} | {'Brute Force':>14} | {'CP-SAT':>10} | "
          f"{'Ersparnis':>11} | {'Loesungen':>10}")
    print("-" * 84)

    for n, grenze in [(4, 8), (5, 10), (6, 12), (7, 14)]:
        kombis, treffer_bf = rohe_aufzaehlung(n, grenze)
        treffer_cp, verzweigungen, dauer = mit_cp_sat(n, grenze)
        assert treffer_bf == treffer_cp, "Beide Verfahren muessen gleich viele finden!"
        ersparnis = kombis / max(verzweigungen, 1)
        print(f"{n:>3} {grenze:>4} | {kombis:>14,} | {verzweigungen:>10,} | "
              f"{ersparnis:>10.0f}x | {treffer_cp:>10,}")

    print("-" * 84)
    print("'Brute Force' = alle Kombinationen. 'CP-SAT' = tatsaechliche Verzweigungen.")
    print("Der Rest wurde durch Propagation ausgeschlossen, ohne ihn anzusehen.")
    print("=" * 84)

#!/usr/bin/env python3

# CP_SAT_Vertretungssystem.py
"""
Kapitel CP-SAT: Dynamisches Vertretungs- und Einsatzplanungssystem mit CP-SAT.

Eigenschaften:
  * Solver-Status wird ueber SolverStatus aus or_kern.py ausgewertet - dieselbe
    Fallunterscheidung wie bei jedem anderen Solver im Buch
  * Strafkosten-Zerlegung wird ausgewiesen (Erklaerbarkeit, Kapitel Praxisfallen)
  * Fairness-Kriterium ergaenzt (Spannweite der Belastung minimieren)
  * Vorabdiagnose auf offensichtliche Unloesbarkeit
  * Abnahmepruefung des fertigen Plans OHNE den Solver zu fragen

Benoetigt: ortools, pandas, pydantic (ueber or_kern)
"""

from __future__ import annotations

from ortools.sat.python import cp_model
import pandas as pd

from or_kern import Loesung, SolverStatus, status_von_cpsat

# --- 1. Datenbasis --------------------------------------------------------
SLOTS = [1, 2, 3, 4]
FAECHER = ["Mathematik", "Physik", "Mathematik", "Informatik"]
KLASSEN = ["Klasse 8a", "Klasse 10b", "Klasse 7a", "Klasse 9c"]

PERSONAL = ["Frau_Mueller", "Herr_Schmidt", "Frau_Albrecht", "Herr_Bauer", "Frau_Koch"]

QUALIFIKATION = {
    "Frau_Mueller":  {"Mathematik", "Physik"},
    "Herr_Schmidt":  {"Informatik", "Mathematik"},
    "Frau_Albrecht": {"Physik", "Informatik"},
    "Herr_Bauer":    {"Mathematik"},
    "Frau_Koch":     {"Mathematik", "Physik", "Informatik"},
}

VORBELASTUNG = {                      # bereits geleistete Stunden heute
    "Frau_Mueller": 1, "Herr_Schmidt": 0, "Frau_Albrecht": 2,
    "Herr_Bauer": 0, "Frau_Koch": 1,
}

MAX_VERTRETUNGEN = 2
STRAFE_VORBELASTUNG = 50              # je Stunde Vorbelastung und Zuweisung
STRAFE_LOCH = 80                      # je Freistundenloch
STRAFE_UNFAIRNESS = 30                # je Einheit Spannweite der Gesamtbelastung


def pruefe_grundsaetzliche_loesbarkeit() -> bool:
    """Vorabdiagnose: Gibt es fuer jeden Slot ueberhaupt qualifiziertes Personal?"""
    ok = True
    for s, fach in enumerate(FAECHER):
        kandidaten = [p for p in PERSONAL if fach in QUALIFIKATION[p]]
        if not kandidaten:
            print(f"  UNLOESBAR: Fuer Slot {SLOTS[s]} ({fach}) gibt es niemanden.")
            ok = False
    kapazitaet = len(PERSONAL) * MAX_VERTRETUNGEN
    if kapazitaet < len(SLOTS):
        print(f"  UNLOESBAR: Kapazitaet {kapazitaet} < {len(SLOTS)} offene Stunden.")
        ok = False
    return ok


def pruefe_plan(plan: dict[int, str]) -> list[str]:
    """Prueft den fertigen Plan gegen die harten Regeln - ohne den Solver.

    Der Solver kann nur pruefen, was ihm gesagt wurde. Diese Funktion prueft
    gegen die ANFORDERUNG und benutzt dafuer bewusst keine Modellvariable
    (siehe Kapitel Praxisfallen). Leere Liste heisst bestanden.
    """
    beanstandungen: list[str] = []
    if sorted(plan) != list(range(len(SLOTS))):
        beanstandungen.append(f"nicht jede Stunde genau einmal besetzt: {sorted(plan)}")
        return beanstandungen
    for s, person in plan.items():
        if FAECHER[s] not in QUALIFIKATION[person]:
            beanstandungen.append(f"{person} ist nicht fuer {FAECHER[s]} qualifiziert")
    for person in PERSONAL:
        anzahl = sum(1 for p in plan.values() if p == person)
        if anzahl > MAX_VERTRETUNGEN:
            beanstandungen.append(f"{person} hat {anzahl} Vertretungen "
                                  f"(hoechstens {MAX_VERTRETUNGEN})")
    return beanstandungen


def plane_vertretung() -> list[dict[str, str]] | None:
    if not pruefe_grundsaetzliche_loesbarkeit():
        return None

    modell = cp_model.CpModel()

    # Entscheidungsvariablen: x[person, slot] = 1  <=>  Person uebernimmt Slot
    x = {(p, s): modell.NewBoolVar(f"zuweisung_{p}_slot{s+1}")
         for p in PERSONAL for s in range(len(SLOTS))}

    # --- HARTE NEBENBEDINGUNGEN ------------------------------------------
    # H1: Jede Stunde genau einmal besetzen
    for s in range(len(SLOTS)):
        modell.AddExactlyOne(x[p, s] for p in PERSONAL)

    # H2: Qualifikation - unqualifizierte Zuweisung ausschliessen
    for s, fach in enumerate(FAECHER):
        for p in PERSONAL:
            if fach not in QUALIFIKATION[p]:
                modell.Add(x[p, s] == 0)

    # H3: Hoechstens MAX_VERTRETUNGEN Stunden je Person
    for p in PERSONAL:
        modell.Add(sum(x[p, s] for s in range(len(SLOTS))) <= MAX_VERTRETUNGEN)

    # (H4 "nicht an zwei Orten gleichzeitig" ist durch H1 bereits erfuellt:
    #  jeder Slot hat genau eine Person, und eine Person kann pro Slot nur
    #  eine Variable auf 1 setzen.)

    # --- WEICHE ZIELE ----------------------------------------------------
    strafterme = []

    # W1: Vorbelastete Personen schonen
    strafe_last = []
    for p in PERSONAL:
        zuweisungen = sum(x[p, s] for s in range(len(SLOTS)))
        strafe_last.append(VORBELASTUNG[p] * STRAFE_VORBELASTUNG * zuweisungen)
    strafterme.extend(strafe_last)

    # W2: Freistundenloecher vermeiden (Muster: arbeitet, frei, arbeitet)
    loch_variablen = []
    for p in PERSONAL:
        for s in range(len(SLOTS) - 2):
            loch = modell.NewBoolVar(f"loch_{p}_{s}")
            # Reifizierung: loch == 1  <=>  (x_s AND NOT x_{s+1} AND x_{s+2})
            modell.AddBoolAnd([x[p, s], x[p, s + 1].Not(), x[p, s + 2]]).OnlyEnforceIf(loch)
            modell.AddBoolOr([x[p, s].Not(), x[p, s + 1], x[p, s + 2].Not()]) \
                  .OnlyEnforceIf(loch.Not())
            loch_variablen.append(loch)
            strafterme.append(loch * STRAFE_LOCH)

    # W3: Fairness - Spannweite der Gesamtbelastung minimieren
    gesamtlast = {}
    for p in PERSONAL:
        last = modell.NewIntVar(0, len(SLOTS) + max(VORBELASTUNG.values()), f"last_{p}")
        modell.Add(last == VORBELASTUNG[p] + sum(x[p, s] for s in range(len(SLOTS))))
        gesamtlast[p] = last
    max_last = modell.NewIntVar(0, 10, "max_last")
    min_last = modell.NewIntVar(0, 10, "min_last")
    modell.AddMaxEquality(max_last, list(gesamtlast.values()))
    modell.AddMinEquality(min_last, list(gesamtlast.values()))
    spannweite = modell.NewIntVar(0, 10, "spannweite")
    modell.Add(spannweite == max_last - min_last)
    strafterme.append(spannweite * STRAFE_UNFAIRNESS)

    modell.Minimize(sum(strafterme))

    # --- LOESEN -----------------------------------------------------------
    loeser = cp_model.CpSolver()
    loeser.parameters.max_time_in_seconds = 5.0
    # Reproduzierbarkeit: ein Arbeiter, fester Startwert. Ohne das kann CP-SAT
    # bei mehreren gleich guten Plaenen von Lauf zu Lauf einen anderen liefern
    # (siehe JobShop_Intervalle.py im naechsten Abschnitt).
    loeser.parameters.num_workers = 1
    loeser.parameters.random_seed = 1
    rohstatus = loeser.Solve(modell)

    # Der Rueckgabewert wird in die gemeinsame Sprache uebersetzt. Ab hier
    # sieht die Auswertung genauso aus wie bei HiGHS, GLOP oder CVXPY.
    status = status_von_cpsat(rohstatus)
    loesung = Loesung(
        status=status,
        zielwert=loeser.ObjectiveValue() if status.brauchbar else None,
        schranke=loeser.BestObjectiveBound() if status.brauchbar else None,
        laufzeit=loeser.WallTime())

    print("=" * 78)
    print("     DYNAMISCHER VERTRETUNGSPLAN (CP-SAT OPTIMIERT)")
    print("=" * 78)

    if status.modellfehler:
        print(f"Kein zulaessiger Plan moeglich (Status: {status.value}).")
        print("Die harten Regeln widersprechen sich. Naechster Schritt: eine Regel")
        print("weich machen (Kapitel Praxisfallen) oder den Konflikt eingrenzen")
        print("(Anhang Fehlerdiagnose).")
        return None
    if not status.brauchbar:
        print(f"Solver lieferte kein Ergebnis (Status: {status.value}).")
        print("Zeitlimit erhoehen oder das Modell vereinfachen.")
        return None

    guete = ("beweisbar optimal" if status is SolverStatus.OPTIMAL
             else "zulaessig, aber nicht bewiesen - Gap beachten")
    print(f"Solver-Status: {loeser.StatusName(rohstatus)} -> {status.value} ({guete})")
    print(f"{loesung.als_bericht()}\n")

    # --- Plan ausgeben ----------------------------------------------------
    zuweisung = {s: next(p for p in PERSONAL if loeser.Value(x[p, s]) == 1)
                 for s in range(len(SLOTS))}
    plan = [{
        "Stunde": f"Std. {SLOTS[s]}",
        "Klasse": KLASSEN[s],
        "Fach": FAECHER[s],
        "Vertretung": zuweisung[s],
        "Vorbelastung": f"{VORBELASTUNG[zuweisung[s]]} Std.",
    } for s in range(len(SLOTS))]
    print(pd.DataFrame(plan).to_string(index=False))

    # --- Abnahmepruefung gegen die Anforderung ----------------------------
    beanstandungen = pruefe_plan(zuweisung)
    if beanstandungen:
        print("\nABNAHMEPRUEFUNG FEHLGESCHLAGEN:")
        for beanstandung in beanstandungen:
            print(f"  - {beanstandung}")
        return None
    print("\nAbnahmepruefung: bestanden (Besetzung, Qualifikation, Hoechstzahl).")

    # --- Erklaerbarkeit: Strafkosten aufschluesseln -----------------------
    kosten_last = sum(VORBELASTUNG[p] * STRAFE_VORBELASTUNG
                      * sum(loeser.Value(x[p, s]) for s in range(len(SLOTS)))
                      for p in PERSONAL)
    kosten_loecher = sum(loeser.Value(v) for v in loch_variablen) * STRAFE_LOCH
    kosten_fairness = loeser.Value(spannweite) * STRAFE_UNFAIRNESS

    print("\n--- Woraus bestehen die Strafkosten? ---")
    print(f"  Vorbelastung geschont:  {kosten_last:5.0f} Punkte")
    print(f"  Freistundenloecher:     {kosten_loecher:5.0f} Punkte "
          f"({sum(loeser.Value(v) for v in loch_variablen)} Loch/Loecher)")
    print(f"  Fairness (Spannweite {loeser.Value(spannweite)}): {kosten_fairness:5.0f} Punkte")
    print(f"  {'Summe':<23} {loeser.ObjectiveValue():5.0f} Punkte")

    print("\n--- Auslastung nach Optimierung ---")
    for p in PERSONAL:
        heute = sum(loeser.Value(x[p, s]) for s in range(len(SLOTS)))
        print(f"  * {p:<15}: {VORBELASTUNG[p]} vorher + {heute} Vertretung "
              f"= {VORBELASTUNG[p] + heute} Stunden")

    print(f"\nSuchstatistik: {loeser.NumBranches()} Verzweigungen, "
          f"{loeser.NumConflicts()} Konflikte, {loeser.WallTime():.3f} s")
    print("=" * 78)
    return plan


if __name__ == "__main__":
    plane_vertretung()

#!/usr/bin/env python3

# JobShop_Intervalle.py
"""
Kapitel CP-SAT: Job-Shop-Scheduling mit Intervallvariablen.
Zeigt NewIntervalVar, AddNoOverlap und die Minimierung des Makespan.
"""

import collections

from ortools.sat.python import cp_model

# (Maschine, Dauer) je Arbeitsgang, in der Reihenfolge des Auftrags
AUFTRAEGE = [
    [(0, 3), (1, 2), (2, 2)],      # Auftrag 0
    [(0, 2), (2, 1), (1, 4)],      # Auftrag 1
    [(1, 4), (2, 3)],              # Auftrag 2
]
MASCHINENNAMEN = ["Fraese", "Dreherei", "Lackiererei"]


def loese_jobshop():
    anzahl_maschinen = 1 + max(m for auftrag in AUFTRAEGE for m, _ in auftrag)
    horizont = sum(dauer for auftrag in AUFTRAEGE for _, dauer in auftrag)

    modell = cp_model.CpModel()
    Gang = collections.namedtuple("Gang", "start ende intervall")
    plaene = {}
    maschinen_intervalle = collections.defaultdict(list)

    # --- Intervallvariablen anlegen --------------------------------------
    for a, auftrag in enumerate(AUFTRAEGE):
        for g, (maschine, dauer) in enumerate(auftrag):
            start = modell.NewIntVar(0, horizont, f"start_{a}_{g}")
            ende = modell.NewIntVar(0, horizont, f"ende_{a}_{g}")
            # Ein Intervall koppelt start + dauer == ende automatisch
            intervall = modell.NewIntervalVar(start, dauer, ende, f"intervall_{a}_{g}")
            plaene[a, g] = Gang(start, ende, intervall)
            maschinen_intervalle[maschine].append(intervall)

    # --- H1: eine Maschine bearbeitet nur einen Gang gleichzeitig --------
    for maschine in range(anzahl_maschinen):
        modell.AddNoOverlap(maschinen_intervalle[maschine])

    # --- H2: Reihenfolge innerhalb eines Auftrags einhalten -------------
    for a, auftrag in enumerate(AUFTRAEGE):
        for g in range(len(auftrag) - 1):
            modell.Add(plaene[a, g + 1].start >= plaene[a, g].ende)

    # --- Ziel: Makespan minimieren ---------------------------------------
    makespan = modell.NewIntVar(0, horizont, "makespan")
    modell.AddMaxEquality(makespan, [plaene[a, len(auftrag) - 1].ende
                                     for a, auftrag in enumerate(AUFTRAEGE)])
    modell.Minimize(makespan)

    loeser = cp_model.CpSolver()
    loeser.parameters.max_time_in_seconds = 10.0
    # Es gibt mehrere Plaene mit demselben kuerzesten Makespan. Welchen CP-SAT
    # findet, haengt sonst davon ab, welcher seiner parallelen Suchstraenge
    # zuerst fertig wird - dasselbe Programm liefert dann von Lauf zu Lauf
    # verschiedene (gleich gute) Plaene. Fuer ein reproduzierbares Buchbeispiel
    # fixieren wir beides. Im Produktivbetrieb laesst man die Standardwerte
    # stehen: Mehrere Arbeiter sind dort deutlich schneller.
    loeser.parameters.num_workers = 1
    loeser.parameters.random_seed = 1
    status = loeser.Solve(modell)

    if status not in (cp_model.OPTIMAL, cp_model.FEASIBLE):
        print(f"Kein Plan gefunden: {loeser.StatusName(status)}")
        return

    print("=" * 74)
    print("  JOB-SHOP-SCHEDULING MIT INTERVALLVARIABLEN")
    print("=" * 74)
    print(f"Status: {loeser.StatusName(status)} | "
          f"Kuerzeste Gesamtdauer (Makespan): {loeser.Value(makespan)} Zeiteinheiten\n")

    # --- Gantt-Diagramm als Text ------------------------------------------
    dauer_gesamt = loeser.Value(makespan)
    print(f"{'Maschine':<14}|" + "".join(f"{t:<3}" for t in range(dauer_gesamt)))
    print("-" * (15 + 3 * dauer_gesamt))
    for maschine in range(anzahl_maschinen):
        zeile = [" . "] * dauer_gesamt
        for a, auftrag in enumerate(AUFTRAEGE):
            for g, (m, dauer) in enumerate(auftrag):
                if m == maschine:
                    beginn = loeser.Value(plaene[a, g].start)
                    for t in range(beginn, beginn + dauer):
                        zeile[t] = f" A{a}"
        print(f"{MASCHINENNAMEN[maschine]:<14}|" + "".join(zeile))

    print("\n--- Detailplan ---")
    for a, auftrag in enumerate(AUFTRAEGE):
        teile = []
        for g, (maschine, dauer) in enumerate(auftrag):
            beginn = loeser.Value(plaene[a, g].start)
            teile.append(f"{MASCHINENNAMEN[maschine]} {beginn}-{beginn + dauer}")
        print(f"  Auftrag {a}: " + "  ->  ".join(teile))

    # --- Untere Schranke zum Vergleich ------------------------------------
    laengster_auftrag = max(sum(d for _, d in a) for a in AUFTRAEGE)
    hoechste_maschinenlast = max(
        sum(d for auftrag in AUFTRAEGE for m, d in auftrag if m == maschine)
        for maschine in range(anzahl_maschinen))
    schranke = max(laengster_auftrag, hoechste_maschinenlast)
    print(f"\nUntere Schranke (laengster Auftrag / hoechste Maschinenlast): {schranke}")
    print(f"Erreichter Makespan: {loeser.Value(makespan)}"
          f"{'  -> beweisbar bestmoeglich' if loeser.Value(makespan) == schranke else ''}")
    print("=" * 74)


if __name__ == "__main__":
    loese_jobshop()

#!/usr/bin/env python3

# CP_SAT_Statusfaelle.py
"""
Kapitel CP-SAT: Alle fuenf Antworten von CP-SAT - und was jede bedeutet.

Der haeufigste Fehler beim Einsatz von CP-SAT im Betrieb ist eine Zeile wie

    if status == cp_model.OPTIMAL:
        ...verwerte die Loesung...

Sie wirft in drei von fuenf Faellen etwas weg, das man haette gebrauchen
koennen, und verschweigt in einem weiteren Fall einen Modellierungsfehler.

Dieses Programm erzeugt jeden Statusfall absichtlich und zeigt die passende
Reaktion:

    OPTIMAL        beweisbar bestmoeglich          -> ausfuehren
    FEASIBLE       zulaessig, Zeit war um          -> Gap pruefen, dann entscheiden
    INFEASIBLE     kein Plan existiert             -> harte Regeln lockern
    MODEL_INVALID  das Modell ist fehlerhaft       -> Programmierfehler beheben
    UNKNOWN        nichts gefunden, Zeit war um    -> mehr Zeit oder Heuristik

Benoetigt: numpy, ortools
"""

from __future__ import annotations

from dataclasses import dataclass

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


@dataclass
class Auswertung:
    """Was nach einem CP-SAT-Lauf tatsaechlich feststeht."""
    status: str
    hat_loesung: bool
    beweisbar_optimal: bool
    zielwert: float | None
    schranke: float | None
    gap: float | None
    empfehlung: str


def werte_aus(loeser: cp_model.CpSolver, status: int,
              mit_zielfunktion: bool) -> Auswertung:
    """Uebersetzt einen CP-SAT-Status in eine Handlungsempfehlung.

    Genau diese Funktion fehlt in den meisten Projekten. Sie ist kurz, aber
    sie ist der Unterschied zwischen einem Prototyp und einem System, das
    nachts ohne Aufsicht laeuft.
    """
    name = loeser.StatusName(status)
    hat_loesung = status in (cp_model.OPTIMAL, cp_model.FEASIBLE)

    zielwert = schranke = gap = None
    if hat_loesung and mit_zielfunktion:
        zielwert = loeser.ObjectiveValue()
        schranke = loeser.BestObjectiveBound()
        nenner = max(abs(zielwert), 1e-9)
        gap = abs(zielwert - schranke) / nenner

    if status == cp_model.OPTIMAL:
        empfehlung = "Ausfuehren. Besser geht es nachweislich nicht."
    elif status == cp_model.FEASIBLE:
        empfehlung = (f"Brauchbar. Der Plan ist hoechstens {gap * 100:.2f} % "
                      f"schlechter als das theoretisch Bestmoegliche."
                      if gap is not None else
                      "Brauchbar, aber nicht als optimal bewiesen.")
    elif status == cp_model.INFEASIBLE:
        empfehlung = ("KEIN Plan existiert. Das ist eine Aussage ueber Ihre "
                      "harten Regeln, nicht ueber die Rechenzeit: Auch "
                      "unendlich viel Zeit wuerde nichts aendern. Regeln "
                      "lockern oder weich machen (Anhang C).")
    elif status == cp_model.MODEL_INVALID:
        empfehlung = ("Das MODELL ist fehlerhaft, nicht das Problem. "
                      "Typisch: eine Variable mit leerem Wertebereich. "
                      "model.Validate() nennt die Ursache.")
    else:                                    # UNKNOWN
        empfehlung = ("Nichts gefunden - aber auch nicht widerlegt. Mehr Zeit "
                      "geben, das Modell vereinfachen oder mit einer "
                      "Heuristik vorbelegen (AddHint).")

    return Auswertung(
        status=name,
        hat_loesung=hat_loesung,
        beweisbar_optimal=status == cp_model.OPTIMAL,
        zielwert=zielwert, schranke=schranke, gap=gap,
        empfehlung=empfehlung,
    )


def loese(baue, zeitlimit: float, mit_zielfunktion: bool) -> Auswertung:
    modell = cp_model.CpModel()
    baue(modell)
    loeser = cp_model.CpSolver()
    loeser.parameters.max_time_in_seconds = zeitlimit
    # Fuer reproduzierbare Buchausgaben; im Betrieb weglassen.
    loeser.parameters.num_workers = 1
    loeser.parameters.random_seed = 1
    status = loeser.Solve(modell)
    return werte_aus(loeser, status, mit_zielfunktion)


# --- Die fuenf Faelle --------------------------------------------------------

def dienstplan_loesbar(modell: cp_model.CpModel) -> None:
    """Vier Personen, fuenf Dienste, hoechstens zwei je Person - loesbar."""
    personen, dienste = 4, 5
    x = [[modell.NewBoolVar(f"x_{i}_{j}") for j in range(dienste)]
         for i in range(personen)]
    for j in range(dienste):
        modell.AddExactlyOne(x[i][j] for i in range(personen))
    for i in range(personen):
        modell.Add(sum(x[i]) <= 2)


def dienstplan_unloesbar(modell: cp_model.CpModel) -> None:
    """Dieselben fuenf Dienste, aber hoechstens EIN Dienst je Person.

    Vier Personen koennen zusammen hoechstens vier Dienste uebernehmen - fuer
    fuenf Dienste reicht das nicht. Kein Solver der Welt findet hier etwas.
    """
    personen, dienste = 4, 5
    x = [[modell.NewBoolVar(f"x_{i}_{j}") for j in range(dienste)]
         for i in range(personen)]
    for j in range(dienste):
        modell.AddExactlyOne(x[i][j] for i in range(personen))
    for i in range(personen):
        modell.Add(sum(x[i]) <= 1)


def modell_fehlerhaft(modell: cp_model.CpModel) -> None:
    """Eine Variable mit leerem Wertebereich: untere Schranke > obere.

    Im Alltag entsteht das durch einen Rechenfehler in den Grenzen, etwa
    NewIntVar(mindestbesetzung, kapazitaet, ...) bei falsch sortierten Daten.
    """
    kaputt = modell.NewIntVar(5, 2, "leerer_bereich")
    modell.Add(kaputt >= 0)


def lastverteilung(modell: cp_model.CpModel) -> None:
    """120 Auftraege auf 11 Maschinen - gross genug, dass der Optimalitaets-
    beweis laenger dauert als das Finden einer sehr guten Loesung."""
    rng = np.random.default_rng(3)
    dauer = rng.integers(100, 9000, 120)
    n, maschinen = len(dauer), 11
    obergrenze = int(dauer.sum())

    x = [[modell.NewBoolVar(f"x_{i}_{k}") for k in range(maschinen)]
         for i in range(n)]
    for i in range(n):
        modell.AddExactlyOne(x[i])
    belegung = [modell.NewIntVar(0, obergrenze, f"l_{k}") for k in range(maschinen)]
    for k in range(maschinen):
        modell.Add(belegung[k] == sum(int(dauer[i]) * x[i][k] for i in range(n)))
    makespan = modell.NewIntVar(0, obergrenze, "makespan")
    modell.AddMaxEquality(makespan, belegung)
    modell.Minimize(makespan)


def zahlpartition(modell: cp_model.CpModel) -> None:
    """46 grosse Zahlen exakt in zwei gleich schwere Haelften teilen.

    Ein beruehmt schweres Problem: Es gibt keine Zielfunktion, entweder es
    geht auf oder nicht - und beides zu entscheiden dauert lange.
    """
    rng = np.random.default_rng(1)
    gewichte = rng.integers(10**6, 2 * 10**6, 46)
    x = [modell.NewBoolVar(f"x_{i}") for i in range(len(gewichte))]
    modell.Add(sum(int(gewichte[i]) * x[i] for i in range(len(gewichte)))
               == int(gewichte.sum() // 2))


def zeige(titel: str, a: Auswertung) -> None:
    print(f"\n{titel}")
    print(f"  Status               {a.status}")
    print(f"  Loesung vorhanden    {'ja' if a.hat_loesung else 'nein'}")
    if a.zielwert is not None:
        print(f"  Zielwert / Schranke  {a.zielwert:,.0f} / {a.schranke:,.0f}"
              f"   (Gap {a.gap * 100:.3f} %)")
    print(f"  -> {a.empfehlung}")


if __name__ == "__main__":
    print("=" * 78)
    print("  DIE FUENF ANTWORTEN VON CP-SAT")
    print("=" * 78)

    zeige("[1] Dienstplan, hoechstens 2 Dienste je Person",
          loese(dienstplan_loesbar, 10.0, mit_zielfunktion=False))

    zeige("[2] Lastverteilung 120 Auftraege / 11 Maschinen, 0,5 s Limit",
          loese(lastverteilung, 0.5, mit_zielfunktion=True))

    zeige("[3] Derselbe Fall mit 10 s Limit",
          loese(lastverteilung, 10.0, mit_zielfunktion=True))

    zeige("[4] Dienstplan, hoechstens 1 Dienst je Person",
          loese(dienstplan_unloesbar, 10.0, mit_zielfunktion=False))

    zeige("[5] Variable mit leerem Wertebereich",
          loese(modell_fehlerhaft, 10.0, mit_zielfunktion=False))

    zeige("[6] Zahlpartition mit 46 grossen Zahlen, 2 s Limit",
          loese(zahlpartition, 2.0, mit_zielfunktion=False))

    print("\n" + "=" * 78)
    print("  DIE DREI, DIE MAN NICHT VERWECHSELN DARF")
    print("=" * 78)
    print("INFEASIBLE  Es gibt keine Loesung. Eine Aussage ueber Ihr MODELL.")
    print("            Mehr Rechenzeit aendert daran nichts.")
    print("UNKNOWN     Es wurde keine gefunden. Eine Aussage ueber die ZEIT.")
    print("            Ob es eine gibt, ist offen.")
    print("MODEL_INVALID  Das Modell ist gar kein gueltiges Modell. Eine")
    print("            Aussage ueber Ihren CODE - immer ein Programmierfehler.")
    print()
    print("Wer diese drei in einem 'else: return None' zusammenfasst, verliert")
    print("genau die Information, die zur Fehlersuche noetig waere.")
    print("=" * 78)

#!/usr/bin/env python3

# Strafgewichte.py
"""
Kapitel CP-SAT: Warum ein Strafgewicht keine Wichtigkeit ausdrueckt, sondern
einen Wechselkurs.

Weiche Regeln kommen als Strafterm in die Zielfunktion. Die Frage, die dabei
regelmaessig falsch beantwortet wird, lautet: "Wie gross muss die Strafe sein?"
Die uebliche Antwort - "je wichtiger, desto groesser" - klingt vernuenftig und
fuehrt zuverlaessig zu Plaenen, die niemand unterschreiben will.

Der Grund: Zwei Strafgewichte legen nicht Wichtigkeiten fest, sondern einen
UMRECHNUNGSKURS. Steht die Serienregel bei 1 und der Wunschfrei-Tag bei 10,
dann hat man dem Solver woertlich gesagt: "Zehn Verstoesse gegen die
Serienregel sind so schlimm wie ein abgelehnter Wunsch." Er nimmt das
ernst - und handelt danach.

Szenario: Schichtplan, 6 Personen, 14 Tage, 4 Besetzungen pro Tag.
  HART   jeder Tag genau 4 Personen; hoechstens 11 Dienste je Person
  WEICH  keine drei Dienste in Folge (Betriebsvereinbarung)
  WEICH  Wunschfrei-Tage einhalten (26 Antraege)

Benoetigt: numpy, ortools
"""

from __future__ import annotations

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

PERSONEN, TAGE, PRO_TAG = 6, 14, 4
MAX_DIENSTE = 11
STRAFE_WUNSCH = 10                    # bleibt fest - variiert wird die Serienstrafe

RNG = np.random.default_rng(2)
WUNSCHFREI = sorted({(p, t) for p in range(PERSONEN) for t in range(TAGE)
                     if RNG.random() < 0.35})


def plane(strafe_serie: int, zeitlimit: float = 25.0) -> dict:
    """Baut den Schichtplan mit dem angegebenen Gewicht fuer die Serienregel."""
    modell = cp_model.CpModel()
    x = {(p, t): modell.NewBoolVar(f"x_{p}_{t}")
         for p in range(PERSONEN) for t in range(TAGE)}

    # --- harte Regeln --------------------------------------------------
    for t in range(TAGE):
        modell.Add(sum(x[p, t] for p in range(PERSONEN)) == PRO_TAG)
    for p in range(PERSONEN):
        modell.Add(sum(x[p, t] for t in range(TAGE)) <= MAX_DIENSTE)

    # --- weiche Regel 1: keine drei Dienste in Folge --------------------
    # serie[p, t] = 1  <=>  Person p arbeitet an t, t+1 UND t+2.
    # Die Ungleichung erzwingt das nur in eine Richtung; das genuegt, weil
    # der Solver serie minimiert - er wuerde die Variable nie freiwillig
    # auf 1 setzen.
    strafterme = []
    serien = []
    for p in range(PERSONEN):
        for t in range(TAGE - 2):
            serie = modell.NewBoolVar(f"serie_{p}_{t}")
            modell.Add(x[p, t] + x[p, t + 1] + x[p, t + 2] <= 2 + serie)
            serien.append(serie)
            strafterme.append(serie * strafe_serie)

    # --- weiche Regel 2: Wunschfrei ------------------------------------
    wunschverstoesse = []
    for p, t in WUNSCHFREI:
        wunschverstoesse.append(x[p, t])
        strafterme.append(x[p, t] * STRAFE_WUNSCH)

    modell.Minimize(sum(strafterme))

    loeser = cp_model.CpSolver()
    loeser.parameters.max_time_in_seconds = zeitlimit
    loeser.parameters.num_workers = 1          # reproduzierbare Buchausgabe
    loeser.parameters.random_seed = 1
    status = loeser.Solve(modell)
    if status not in (cp_model.OPTIMAL, cp_model.FEASIBLE):
        raise RuntimeError(f"Kein Plan: {loeser.StatusName(status)}")

    return {
        "status": loeser.StatusName(status),
        "zielwert": loeser.ObjectiveValue(),
        "serien": sum(loeser.Value(v) for v in serien),
        "wunsch_verletzt": sum(loeser.Value(v) for v in wunschverstoesse),
    }


if __name__ == "__main__":
    print("=" * 78)
    print("  STRAFGEWICHTE SIND WECHSELKURSE, KEINE WICHTIGKEITEN")
    print("=" * 78)
    print(f"{PERSONEN} Personen, {TAGE} Tage, {PRO_TAG} Besetzungen pro Tag "
          f"= {TAGE * PRO_TAG} Dienste")
    print(f"Das sind {TAGE * PRO_TAG / PERSONEN:.1f} Dienste je Person - der Plan ist eng.")
    print(f"{len(WUNSCHFREI)} Wunschfrei-Antraege, Strafe je abgelehntem Wunsch: "
          f"{STRAFE_WUNSCH}\n")

    print(f"{'Strafe je 3er-Serie':>20} {'3er-Serien':>12} {'Wunsch verletzt':>17} "
          f"{'Zielwert':>10}")
    print("-" * 78)
    ergebnisse = {}
    for strafe in (1, 2, 5, 10, 30, 100):
        e = plane(strafe)
        ergebnisse[strafe] = e
        print(f"{strafe:>20} {e['serien']:>12} {e['wunsch_verletzt']:>17} "
              f"{e['zielwert']:>10.0f}")
    print("-" * 78)

    niedrig, hoch = ergebnisse[1], ergebnisse[30]
    print(f"\nBei Strafe 1 baut der Solver {niedrig['serien']} Drei-Tage-Serien ein.")
    print("Die Betriebsvereinbarung steht im Modell - und wird trotzdem")
    print("systematisch verletzt. Das ist kein Fehler des Solvers: Bei einem")
    print(f"Gewicht von 1 gegen {STRAFE_WUNSCH} lohnt sich jede Serie, solange sie")
    print("auch nur einen Zehntel-Wunsch rettet.")
    print(f"\nBei Strafe 30 sind es {hoch['serien']} Serien - dafuer werden")
    print(f"{hoch['wunsch_verletzt']} statt {niedrig['wunsch_verletzt']} Wuensche abgelehnt.")
    print("Beide Plaene sind optimal. Sie beantworten nur verschiedene Fragen.")
    print()
    print("Die Frage lautet also nie 'wie wichtig ist mir diese Regel?', sondern:")
    print("  'Wie viele abgelehnte Wuensche bin ich bereit zu akzeptieren,")
    print("   um eine Drei-Tage-Serie zu vermeiden?'")
    print("Wer darauf keine Zahl nennen kann, hat das Problem noch nicht")
    print("verstanden - und sollte diese Tabelle dem Auftraggeber vorlegen,")
    print("statt ein Gewicht zu raten.")
    print()
    print("Und die wichtigste Konsequenz: Eine Regel, die NIE gebrochen werden")
    print("darf, gehoert nicht in die Zielfunktion, sondern unter die harten")
    print("Nebenbedingungen. Alles, was einen Preis hat, wird irgendwann gekauft.")
    print("=" * 78)