Kapitel 6: Gemischt-ganzzahlige Optimierung — Diskrete Entscheidungen und Branch-and-Bound

#!/usr/bin/env python3

# Runden_Gegenbeispiel.py
"""
Kapitel MILP: Wie schlecht ist Runden wirklich?

Ersetzt die unbelegte Behauptung "20-50 % Verlust" durch eine Messung ueber
viele Zufallsinstanzen.
"""

import numpy as np
from scipy.optimize import linprog


def erzeuge_instanz(n, m, rng):
    """Zufaelliges Rucksack-aehnliches MILP mit kleinen Zahlen."""
    c = rng.integers(3, 20, size=n).astype(float)          # Ertraege
    A = rng.integers(1, 9, size=(m, n)).astype(float)      # Verbraeuche
    b = (A.sum(axis=1) * rng.uniform(0.25, 0.45)).round()  # knappe Kapazitaeten
    return c, A, b


def loese_lp(c, A, b, ganzzahlig=False):
    """LP-Relaxation oder exaktes MILP ueber HiGHS."""
    n = len(c)
    ergebnis = linprog(
        c=-c, A_ub=A, b_ub=b, bounds=[(0, None)] * n,
        integrality=np.ones(n) if ganzzahlig else None,
        method="highs")
    return (-ergebnis.fun, ergebnis.x) if ergebnis.success else (None, None)


def abrunden_und_reparieren(x_lp, c, A, b):
    """Naive Strategie: abrunden, dann gierig auffuellen, solange zulaessig."""
    x = np.floor(x_lp + 1e-9)
    verbessert = True
    while verbessert:                       # gierig auffuellen
        verbessert = False
        for j in np.argsort(-c):            # bester Ertrag zuerst
            kandidat = x.copy()
            kandidat[j] += 1
            if np.all(A @ kandidat <= b + 1e-9):
                x = kandidat
                verbessert = True
                break
    return c @ x, x


if __name__ == "__main__":
    rng = np.random.default_rng(2026)
    print("=" * 82)
    print("  WIE TEUER IST RUNDEN? (200 Zufallsinstanzen je Groesse)")
    print("=" * 82)
    print(f"{'n x m':>8} | {'Aufrunden unzul.':>17} | {'Abrunden: mittl.':>17} | "
          f"{'schlimmster':>12} | {'gierig':>8}")
    print(f"{'':>8} | {'':>17} | {'Verlust':>17} | {'Fall':>12} | {'Verlust':>8}")
    print("-" * 82)

    for n, m in [(5, 2), (10, 3), (20, 5), (40, 8)]:
        unzulaessig = 0
        verluste_ab, verluste_gierig = [], []

        for _ in range(200):
            c, A, b = erzeuge_instanz(n, m, rng)
            z_lp, x_lp = loese_lp(c, A, b, ganzzahlig=False)
            z_ip, _ = loese_lp(c, A, b, ganzzahlig=True)
            if z_lp is None or z_ip is None or z_ip <= 0:
                continue

            # Variante 1: aufrunden
            x_auf = np.ceil(x_lp - 1e-9)
            if np.any(A @ x_auf > b + 1e-9):
                unzulaessig += 1

            # Variante 2: abrunden
            x_ab = np.floor(x_lp + 1e-9)
            verluste_ab.append(1.0 - (c @ x_ab) / z_ip)

            # Variante 3: abrunden + gierig auffuellen
            z_gierig, _ = abrunden_und_reparieren(x_lp, c, A, b)
            verluste_gierig.append(1.0 - z_gierig / z_ip)

        print(f"{n:>3} x {m:<3} | {unzulaessig/2:>15.1f} % | "
              f"{np.mean(verluste_ab)*100:>15.1f} % | "
              f"{np.max(verluste_ab)*100:>10.1f} % | "
              f"{np.mean(verluste_gierig)*100:>6.1f} %")

    print("-" * 82)
    print("Lesart: 'Aufrunden unzul.' = Anteil der Faelle, in denen die aufgerundete")
    print("Loesung eine Nebenbedingung verletzt. 'Verlust' = Abstand zum exakten Optimum.")
    print("=" * 82)

#!/usr/bin/env python3

# Rucksack.py
"""
Kapitel MILP: Das Rucksackproblem (Knapsack).
Zeigt LP-Relaxation, Branch-and-Bound-Ergebnis und den Preis der Ganzzahligkeit
an einem Beispiel, das man vollstaendig im Kopf nachvollziehen kann.
"""

import numpy as np
from scipy.optimize import linprog

GEGENSTAENDE = ["Zelt", "Schlafsack", "Kocher", "Kamera", "Buch", "Wasserfilter", "Seil"]
NUTZEN =  np.array([40.0, 35.0, 20.0, 30.0,  8.0, 25.0, 12.0])
GEWICHT = np.array([ 6.0,  4.0,  3.0,  2.0,  1.0,  2.0,  3.0])
KAPAZITAET = 11.0        # kg


def loese(ganzzahlig: bool):
    n = len(NUTZEN)
    ergebnis = linprog(
        c=-NUTZEN, A_ub=[GEWICHT], b_ub=[KAPAZITAET],
        bounds=[(0, 1)] * n,                       # jedes Teil hoechstens einmal
        integrality=np.ones(n) if ganzzahlig else None,
        method="highs")
    return -ergebnis.fun, ergebnis.x


if __name__ == "__main__":
    z_lp, x_lp = loese(ganzzahlig=False)
    z_ip, x_ip = loese(ganzzahlig=True)

    print("=" * 72)
    print(f"  RUCKSACKPROBLEM  (Kapazitaet {KAPAZITAET:.0f} kg)")
    print("=" * 72)
    print(f"{'Gegenstand':<14} {'Nutzen':>7} {'kg':>5} {'Nutzen/kg':>10} "
          f"{'LP':>7} {'MILP':>6}")
    print("-" * 72)
    for i, name in enumerate(GEGENSTAENDE):
        # int(round(...)) statt Format "%.0f": vermeidet die Ausgabe "-0"
        print(f"{name:<14} {NUTZEN[i]:>7.0f} {GEWICHT[i]:>5.0f} "
              f"{NUTZEN[i]/GEWICHT[i]:>10.2f} {x_lp[i]:>7.2f} {int(round(x_ip[i])):>6d}")
    print("-" * 72)
    print(f"{'Gesamtnutzen':<14} {'':<7} {'':<5} {'':<10} {z_lp:>7.2f} {z_ip:>6.0f}")
    print(f"{'Gesamtgewicht':<14} {'':<7} {'':<5} {'':<10} "
          f"{GEWICHT @ x_lp:>7.2f} {GEWICHT @ x_ip:>6.0f}")
    print("-" * 72)
    print(f"Obere Schranke aus der LP-Relaxation: {z_lp:.2f}")
    print(f"Bestes ganzzahliges Ergebnis:         {z_ip:.0f}")
    print(f"Preis der Ganzzahligkeit:             {z_lp - z_ip:.2f} "
          f"({(1 - z_ip/z_lp)*100:.1f} %)")

    gebrochen = [GEGENSTAENDE[i] for i in range(len(NUTZEN)) if 1e-6 < x_lp[i] < 1 - 1e-6]
    print(f"\nIn der LP-Loesung nur teilweise eingepackt: {gebrochen}")
    print("Genau hier wuerde Branch-and-Bound verzweigen:")
    print(f"  Ast 1: {gebrochen[0]} bleibt ganz zuhause (x=0)")
    print(f"  Ast 2: {gebrochen[0]} kommt ganz mit    (x=1)")
    print("=" * 72)

#!/usr/bin/env python3

# MILP_Portfolio_Fixgebuehren.py
"""
Kapitel MILP: MILP-Portfolio-Selektion mit Fixkosten und Kardinalitaet.

CP-SAT-freie, gut lesbare Formulierung ueber scipy/HiGHS (kein manueller
CSR-Matrixaufbau, deutlich leichter nachvollziehbar), mit Vergleich gegen die
Loesung OHNE Restriktionen.

Die Auswertung laeuft ueber SolverStatus und das Loesung-Objekt aus
or_kern.py: Alle Statusfaelle werden behandelt, und der MIP-Gap steht im
Bericht - statt eines blossen 'success', das nicht verraet, ob der Solver
fertig geworden ist oder nur aufgegeben hat.

Benoetigt: numpy, pandas, scipy, pydantic (ueber or_kern)
"""

from __future__ import annotations

import time

import numpy as np
import pandas as pd
from scipy.optimize import linprog

from or_kern import Loesung, SolverStatus, status_von_scipy

ANLAGEN = ["US-Aktien", "EU-Aktien", "Emerging-Markets",
           "Staatsanleihen", "Unternehmensanleihen", "Rohstoffe"]
RENDITE = np.array([0.11, 0.08, 0.13, 0.03, 0.05, 0.07])   # erwartet, p.a.

BUDGET = 100_000.0
MIN_POSITION = 10_000.0        # L
MAX_POSITION = 40_000.0        # U  (dient zugleich als Big-M!)
GEBUEHR = 50.0                 # F, je aktivierter Position
MAX_POSITIONEN = 3             # K

N = len(ANLAGEN)


def baue_und_loese(zeitlimit: float | None = None) -> Loesung:
    """
    Variablenreihenfolge:  [x_0..x_{N-1}, y_0..y_{N-1}]
    Zielfunktion (Maximierung -> fuer linprog negiert):
        max  sum(rendite_i * x_i) - GEBUEHR * sum(y_i)
    """
    c = np.concatenate([-RENDITE, np.full(N, GEBUEHR)])     # negiert = Minimierung

    # --- Gleichungsnebenbedingung: Budget vollstaendig investiert -----------
    A_eq = np.zeros((1, 2 * N))
    A_eq[0, :N] = 1.0
    b_eq = np.array([BUDGET])

    zeilen, grenzen = [], []

    # --- Kardinalitaet: sum(y_i) <= K --------------------------------------
    zeile = np.zeros(2 * N)
    zeile[N:] = 1.0
    zeilen.append(zeile)
    grenzen.append(MAX_POSITIONEN)

    # --- Obergrenze (Big-M):  x_i - U*y_i <= 0 -----------------------------
    for i in range(N):
        zeile = np.zeros(2 * N)
        zeile[i] = 1.0
        zeile[N + i] = -MAX_POSITION
        zeilen.append(zeile)
        grenzen.append(0.0)

    # --- Untergrenze:  L*y_i - x_i <= 0   (entspricht x_i >= L*y_i) --------
    for i in range(N):
        zeile = np.zeros(2 * N)
        zeile[i] = -1.0
        zeile[N + i] = MIN_POSITION
        zeilen.append(zeile)
        grenzen.append(0.0)

    A_ub = np.array(zeilen)
    b_ub = np.array(grenzen)

    schranken = [(0.0, MAX_POSITION)] * N + [(0.0, 1.0)] * N
    ganzzahligkeit = np.concatenate([np.zeros(N), np.ones(N)])   # y binaer

    t0 = time.perf_counter()
    ergebnis = linprog(c=c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,
                       bounds=schranken, integrality=ganzzahligkeit,
                       method="highs",
                       options={"time_limit": zeitlimit} if zeitlimit else None)
    laufzeit = time.perf_counter() - t0

    status = status_von_scipy(ergebnis)
    if not status.brauchbar:
        return Loesung(status=status, laufzeit=laufzeit)

    # Zurueck in die Maximierungswelt: Zielwert UND Schranke negieren.
    # Aus beiden zusammen rechnet das Loesung-Objekt den MIP-Gap aus.
    x, y = ergebnis.x[:N], np.round(ergebnis.x[N:])
    return Loesung(
        status=status,
        werte={**{name: float(w) for name, w in zip(ANLAGEN, x)},
               **{f"aktiv:{name}": float(w) for name, w in zip(ANLAGEN, y)}},
        zielwert=float(-ergebnis.fun),
        schranke=float(-ergebnis.mip_dual_bound),
        laufzeit=laufzeit)


def ohne_restriktionen() -> tuple[float, np.ndarray]:
    """Vergleichsfall: nur Budget, keine Gebuehren/Kardinalitaet/Mindestgroesse."""
    ergebnis = linprog(c=-RENDITE, A_eq=[np.ones(N)], b_eq=[BUDGET],
                       bounds=[(0, None)] * N, method="highs")
    status = status_von_scipy(ergebnis)
    if not status.brauchbar:
        raise SystemExit(f"Vergleichsfall nicht loesbar: {status.value}")
    return -ergebnis.fun, ergebnis.x


def pruefe_portfolio(loesung: Loesung, toleranz: float = 1e-6) -> list[str]:
    """Prueft die Loesung gegen die Anforderungen - ohne den Solver zu fragen.

    Bewusst kein assert: Eine Beanstandungsliste laesst sich protokollieren,
    weiterreichen und testen. Ein assert verschwindet ausserdem, sobald
    jemand Python mit -O startet.
    """
    x = np.array([loesung.werte[name] for name in ANLAGEN])
    y = np.array([loesung.werte[f"aktiv:{name}"] for name in ANLAGEN])
    beanstandungen: list[str] = []

    if abs(x.sum() - BUDGET) > 1e-4:
        beanstandungen.append(f"Budget nicht exakt investiert: {x.sum():,.2f}")
    if y.sum() > MAX_POSITIONEN + toleranz:
        beanstandungen.append(f"{y.sum():.0f} Positionen statt hoechstens "
                              f"{MAX_POSITIONEN}")
    for i, name in enumerate(ANLAGEN):
        if abs(y[i] - round(y[i])) > toleranz:
            beanstandungen.append(f"{name}: y = {y[i]!r} ist nicht ganzzahlig")
        elif y[i] and not (MIN_POSITION - toleranz <= x[i]
                           <= MAX_POSITION + toleranz):
            beanstandungen.append(f"{name}: {x[i]:,.2f} EUR verletzt die "
                                  f"Groessengrenzen")
        elif not y[i] and x[i] > toleranz:
            beanstandungen.append(f"{name}: inaktiv, aber {x[i]:,.2f} EUR "
                                  f"investiert")
    return beanstandungen


if __name__ == "__main__":
    loesung = baue_und_loese()

    print("=" * 88)
    print("      OPTIMALE MILP-PORTFOLIO-ALLOKATION MIT FIXGEBUEHREN")
    print("=" * 88)
    print(f"Budget: {BUDGET:,.0f} EUR | max. {MAX_POSITIONEN} Positionen | "
          f"je {MIN_POSITION:,.0f}-{MAX_POSITION:,.0f} EUR | Gebuehr {GEBUEHR:.0f} EUR\n")

    # Zuerst der Status - erst danach interessieren die Zahlen.
    if loesung.status.modellfehler:
        raise SystemExit(f"Das Modell ist nicht loesbar ({loesung.status.value}). "
                         f"Naechster Schritt: Anhang Fehlerdiagnose.")
    if not loesung.status.brauchbar:
        raise SystemExit(f"Keine Loesung erhalten ({loesung.status.value}). "
                         f"Zeitlimit erhoehen oder Modell vereinfachen.")
    if loesung.status is SolverStatus.ZULAESSIG:
        print(f"ACHTUNG: nicht beweisbar optimal - Gap {loesung.gap:.2%}\n")

    x = np.array([loesung.werte[name] for name in ANLAGEN])
    y = np.array([loesung.werte[f"aktiv:{name}"] for name in ANLAGEN], dtype=int)

    tabelle = pd.DataFrame({
        "Anlage": ANLAGEN,
        "Aktiv": ["JA" if y[i] else "-" for i in range(N)],
        "Investition (EUR)": [f"{x[i]:,.0f}" for i in range(N)],
        "Anteil": [f"{x[i]/BUDGET*100:5.1f} %" for i in range(N)],
        "Erw. Rendite": [f"{RENDITE[i]*100:4.1f} %" for i in range(N)],
        "Erw. Ertrag (EUR)": [f"{x[i]*RENDITE[i]:,.0f}" for i in range(N)],
    })
    print(tabelle.to_string(index=False))

    brutto = float(RENDITE @ x)
    gebuehren = float(GEBUEHR * y.sum())
    print("-" * 88)
    print(f"Erwarteter Bruttoertrag: {brutto:>12,.2f} EUR")
    print(f"Ordergebuehren:          {-gebuehren:>12,.2f} EUR ({y.sum()} Positionen)")
    print(f"Netto-Erwartungswert:    {loesung.zielwert:>12,.2f} EUR")
    print(f"\n{loesung.als_bericht()}")

    # --- Vergleich mit dem unbeschraenkten Fall ---------------------------
    z_frei, x_frei = ohne_restriktionen()
    print("-" * 88)
    print(f"Zum Vergleich ohne jede Restriktion (alles in den Bestwert): "
          f"{z_frei:,.2f} EUR")
    print(f"Kosten der Realitaet (Gebuehren, Streuung, Mindestgroessen): "
          f"{z_frei - loesung.zielwert:,.2f} EUR "
          f"({(1 - loesung.zielwert/z_frei)*100:.2f} %)")

    # --- Alle Nebenbedingungen nachpruefen --------------------------------
    beanstandungen = pruefe_portfolio(loesung)
    print("-" * 88)
    if beanstandungen:
        raise SystemExit("Abnahmepruefung fehlgeschlagen:\n  - "
                         + "\n  - ".join(beanstandungen))
    print("Abnahmepruefung: alle Nebenbedingungen geprueft und eingehalten.")
    print("=" * 88)

#!/usr/bin/env python3

# Solverstatus_und_Gap.py
"""
Kapitel MILP: Was tun, wenn der Solver nicht fertig wird?

Bei einem LP kommt entweder eine optimale Loesung oder eine klare Absage.
Bei einem MILP ist der haeufigste Ausgang im Betrieb ein dritter: "Ich habe
eine Loesung, ich weiss aber nicht, ob sie die beste ist - und die Zeit ist
um." Dieses Programm zeigt, wie man mit diesem Fall umgeht.

  1. Alle Statusfaelle explizit behandeln, statt OPTIMAL vorauszusetzen.
  2. Den MIP-Gap lesen: Wie weit kann ich hoechstens danebenliegen?
  3. Messen, was zusaetzliche Rechenzeit ueberhaupt noch bringt.
  4. Warm-Start ausprobieren - und ehrlich messen, ob er etwas bringt.

Beispiel: Standortplanung, 45 moegliche Lager, 120 Kunden (5445 Variablen,
davon 45 binaer).

Benoetigt: numpy, highspy
"""

from __future__ import annotations

import time
from dataclasses import dataclass

import numpy as np
import highspy

RNG = np.random.default_rng(7)

N_LAGER, N_KUNDE = 45, 120
FIXKOSTEN = RNG.uniform(3000, 9000, N_LAGER)
TRANSPORT = RNG.uniform(5, 60, (N_LAGER, N_KUNDE))
BEDARF = RNG.uniform(10, 60, N_KUNDE)
KAPAZITAET = np.full(N_LAGER, BEDARF.sum() * 0.22)


@dataclass
class Ergebnis:
    """Alles, was nach einem Solverlauf ausgewertet werden muss - nicht nur
    der Zielwert."""
    status: str
    brauchbar: bool          # Gibt es ueberhaupt eine zulaessige Loesung?
    beweisbar_optimal: bool
    zielwert: float          # bester gefundener Wert (Incumbent)
    schranke: float          # beste bewiesene Schranke (Dual Bound)
    gap: float               # relativer Abstand zwischen beiden
    knoten: int
    dauer: float
    loesung: np.ndarray


def loese(zeitlimit: float, startloesung: np.ndarray | None = None) -> Ergebnis:
    """Loest das Standortmodell mit Zeitlimit und wertet ALLE Statusfaelle aus."""
    hochschule = highspy.Highs()
    hochschule.setOptionValue("output_flag", False)
    hochschule.setOptionValue("time_limit", zeitlimit)

    anzahl_x = N_LAGER * N_KUNDE
    unendlich = highspy.kHighsInf

    hochschule.addVars(anzahl_x, np.zeros(anzahl_x), np.full(anzahl_x, unendlich))
    hochschule.addVars(N_LAGER, np.zeros(N_LAGER), np.ones(N_LAGER))
    for i in range(N_LAGER):
        hochschule.changeColIntegrality(anzahl_x + i, highspy.HighsVarType.kInteger)
        hochschule.changeColCost(anzahl_x + i, FIXKOSTEN[i])
        for j in range(N_KUNDE):
            hochschule.changeColCost(i * N_KUNDE + j, TRANSPORT[i, j])

    for j in range(N_KUNDE):
        index = np.array([i * N_KUNDE + j for i in range(N_LAGER)], dtype=np.int32)
        hochschule.addRow(BEDARF[j], BEDARF[j], len(index), index, np.ones(len(index)))

    for i in range(N_LAGER):
        index = np.array([i * N_KUNDE + j for j in range(N_KUNDE)] + [anzahl_x + i],
                         dtype=np.int32)
        werte = np.concatenate([np.ones(N_KUNDE), [-KAPAZITAET[i]]])
        hochschule.addRow(-unendlich, 0.0, len(index), index, werte)

    if startloesung is not None:
        hochschule.setSolution(len(startloesung),
                               np.arange(len(startloesung), dtype=np.int32),
                               startloesung)

    t0 = time.perf_counter()
    hochschule.run()
    dauer = time.perf_counter() - t0

    status = hochschule.modelStatusToString(hochschule.getModelStatus())
    info = hochschule.getInfo()

    # Der Kern der Sache: Aus dem Status folgt, WAS man mit dem Ergebnis
    # ueberhaupt anfangen darf.
    brauchbar = status in ("Optimal", "Time limit reached", "Solution limit reached")
    beweisbar_optimal = status == "Optimal"
    if status in ("Infeasible", "Unbounded", "Primal infeasible or unbounded"):
        brauchbar = False

    return Ergebnis(
        status=status,
        brauchbar=brauchbar and info.objective_function_value < unendlich,
        beweisbar_optimal=beweisbar_optimal,
        zielwert=info.objective_function_value,
        schranke=info.mip_dual_bound,
        gap=info.mip_gap,
        knoten=info.mip_node_count,
        dauer=dauer,
        loesung=np.array(hochschule.getSolution().col_value),
    )


def gieriger_startplan() -> tuple[np.ndarray, float]:
    """Eine Faustregel-Loesung, wie sie ein Disponent von Hand erstellen wuerde:
    die guenstigsten Lager oeffnen (Fixkosten je Kapazitaetseinheit), dann
    jeden Kunden dem naechstgelegenen offenen Lager mit Restkapazitaet
    zuordnen. Kein Solver noetig - und in Sekunden fertig.
    """
    reihenfolge = np.argsort(FIXKOSTEN / KAPAZITAET)
    offen: list[int] = []
    for i in reihenfolge:
        offen.append(int(i))
        if KAPAZITAET[offen].sum() >= BEDARF.sum() * 1.05:
            break

    rest = KAPAZITAET.copy()
    x = np.zeros((N_LAGER, N_KUNDE))
    for j in np.argsort(-BEDARF):                 # groesste Kunden zuerst
        for i in sorted(offen, key=lambda i: TRANSPORT[i, j]):
            menge = min(rest[i], BEDARF[j] - x[:, j].sum())
            if menge > 1e-9:
                x[i, j] += menge
                rest[i] -= menge
            if abs(x[:, j].sum() - BEDARF[j]) < 1e-9:
                break

    y = np.zeros(N_LAGER)
    y[offen] = 1.0
    kosten = float((FIXKOSTEN * y).sum() + (TRANSPORT * x).sum())
    return np.concatenate([x.ravel(), y]), kosten


def zeige(titel: str, e: Ergebnis) -> None:
    print(f"\n{titel}")
    print(f"  Status                {e.status}")
    if not e.brauchbar:
        print("  -> KEINE verwertbare Loesung. Nicht weiterrechnen!")
        return
    print(f"  bester Plan (Incumbent)   {e.zielwert:>12,.2f} EUR")
    print(f"  bewiesene Schranke        {e.schranke:>12,.2f} EUR")
    print(f"  MIP-Gap                   {e.gap * 100:>12.3f} %")
    print(f"  Knoten / Zeit             {e.knoten:>12,} / {e.dauer:.2f} s")
    if e.beweisbar_optimal:
        print("  -> beweisbar optimal")
    else:
        print(f"  -> zulaessig, aber nicht bewiesen optimal. Der wahre Bestwert")
        print(f"     liegt zwischen {e.schranke:,.2f} und {e.zielwert:,.2f} EUR.")


if __name__ == "__main__":
    print("=" * 78)
    print("  MIP-GAP UND ZEITLIMIT: STANDORTPLANUNG, 45 LAGER, 120 KUNDEN")
    print("=" * 78)
    print(f"{N_LAGER * N_KUNDE + N_LAGER:,} Variablen, davon {N_LAGER} binaer.")

    kurz = loese(2.0)
    zeige("[1] Zeitlimit 2 Sekunden", kurz)

    lang = loese(60.0)
    zeige("[2] Zeitlimit 60 Sekunden", lang)

    print("\n" + "=" * 78)
    print("  WAS BRINGT MEHR RECHENZEIT?")
    print("=" * 78)
    verbesserung = kurz.zielwert - lang.zielwert
    print(f"Nach 2 Sekunden:  {kurz.zielwert:,.2f} EUR, Gap {kurz.gap * 100:.2f} %")
    print(f"Nach {lang.dauer:.1f} Sekunden: {lang.zielwert:,.2f} EUR, bewiesen optimal")
    print(f"Gewinn durch {lang.dauer - kurz.dauer:.1f} Sekunden mehr Rechenzeit: "
          f"{verbesserung:,.2f} EUR "
          f"({verbesserung / kurz.zielwert * 100:.2f} %)")
    print()
    print("Das ist die Frage, die im Betrieb wirklich zaehlt: Der Gap von")
    print(f"{kurz.gap * 100:.1f} % nach 2 Sekunden ist eine GARANTIE - schlechter als")
    print("dieser Wert kann die Loesung nicht sein. Ob sich die restliche")
    print("Rechenzeit lohnt, entscheidet nicht der Solver, sondern die Anwendung:")
    print("Ein naechtlicher Tourenplan darf eine Stunde rechnen, eine Umplanung")
    print("bei Maschinenausfall hat 30 Sekunden.")

    print("\n" + "=" * 78)
    print("  BRINGT EIN WARM-START ETWAS?")
    print("=" * 78)
    start, start_kosten = gieriger_startplan()
    print(f"Faustregel-Startplan (ohne Solver): {start_kosten:,.2f} EUR")
    print(f"Das sind {(start_kosten / lang.zielwert - 1) * 100:.1f} % ueber dem Optimum.\n")

    warm = loese(60.0, startloesung=start)
    print(f"{'':26} {'Zeit':>9} {'Knoten':>9} {'Ziel':>13}")
    print("-" * 78)
    print(f"{'ohne Warm-Start':<26} {lang.dauer:>8.2f}s {lang.knoten:>9,} "
          f"{lang.zielwert:>13,.2f}")
    print(f"{'mit Warm-Start':<26} {warm.dauer:>8.2f}s {warm.knoten:>9,} "
          f"{warm.zielwert:>13,.2f}")
    print("-" * 78)
    print("Ergebnis: praktisch kein Unterschied. Der Grund ist nicht, dass")
    print("Warm-Starts nichts taugen - sondern dass HiGHS' eigene Heuristiken")
    print("innerhalb der ersten Sekunde bereits eine BESSERE Loesung finden als")
    print("unsere Faustregel. Ein Startwert hilft nur, wenn er besser ist als")
    print("das, was der Solver von allein in derselben Zeit findet.")
    print()
    print("Warm-Starts lohnen sich damit vor allem in zwei Faellen:")
    print("  * Sie haben Domaenenwissen, das der Solver nicht hat (siehe")
    print("    Warmstart_Effekt.py - dort halbiert ein Heuristik-Hinweis die Zeit).")
    print("  * Sie planen laufend neu und der gestrige Plan ist fast noch gueltig.")
    print("In beiden Faellen gilt: MESSEN, nicht glauben.")
    print("=" * 78)

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

#!/usr/bin/env python3

# Big_M_Falle.py
"""
Kapitel MILP: Was ein zu grosses Big-M wirklich anrichtet.

Lehrbuecher warnen vor grossem M mit dem Hinweis "die Relaxation wird
schwach, der Solver braucht mehr Knoten". Das stimmt - ist aber die
harmlosere Haelfte der Wahrheit. Die gefaehrlichere: Bei sehr grossem M
kann der Solver eine Loesung als ganzzahlig ANNEHMEN, in der die
Binaervariablen bei 1e-8 stehen. Dann liefern zugeschaltete Anlagen Ware
aus, waehrend das Modell ihre Fixkosten mit 0 verbucht.

Beispiel: Standortplanung, 12 moegliche Lager, 40 Kunden.
    min  sum_i fix_i * y_i  +  sum_ij kosten_ij * x_ij
    u.d.N.  sum_i x_ij = bedarf_j                    (jeder Kunde beliefert)
            sum_j x_ij <= M_i * y_i                  (Lager offen, wenn es liefert)
            y_i binaer,  x_ij >= 0

Zwei Laeufe:
  1. M knapp gewaehlt (= tatsaechliche Lagerkapazitaet)  -> richtig
  2. M = 1e7-fach zu gross, Presolve abgeschaltet        -> stilles Desaster
Beide werden mit derselben Pruefung kontrolliert, die den Fall auffliegen laesst.

Benoetigt: numpy, highspy
"""

from __future__ import annotations

import time

import numpy as np
import highspy

RNG = np.random.default_rng(11)

N_LAGER, N_KUNDE = 12, 40
FIXKOSTEN = RNG.uniform(2000, 5000, N_LAGER)
TRANSPORT = RNG.uniform(5, 40, (N_LAGER, N_KUNDE))
BEDARF = RNG.uniform(10, 60, N_KUNDE)
KAPAZITAET = BEDARF.sum() * 0.45          # jedes Lager schafft 45 % des Gesamtbedarfs

# Toleranz, ab der ein Solver eine Variable als ganzzahlig durchgehen laesst.
# HiGHS und die meisten anderen verwenden 1e-6 als Standard.
GANZZAHL_TOLERANZ = 1e-6


def loese(big_m: float, presolve: str = "on") -> dict:
    """Baut und loest das Standortmodell. Liefert Loesung und Solverkennzahlen."""
    hochschule = highspy.Highs()
    hochschule.setOptionValue("output_flag", False)
    hochschule.setOptionValue("presolve", presolve)
    hochschule.setOptionValue("time_limit", 300.0)

    anzahl_x = N_LAGER * N_KUNDE
    unendlich = highspy.kHighsInf

    # Spalten: erst alle x_ij, dann die y_i
    hochschule.addVars(anzahl_x, np.zeros(anzahl_x), np.full(anzahl_x, unendlich))
    hochschule.addVars(N_LAGER, np.zeros(N_LAGER), np.ones(N_LAGER))
    for i in range(N_LAGER):
        hochschule.changeColIntegrality(anzahl_x + i,
                                        highspy.HighsVarType.kInteger)
        hochschule.changeColCost(anzahl_x + i, FIXKOSTEN[i])
        for j in range(N_KUNDE):
            hochschule.changeColCost(i * N_KUNDE + j, TRANSPORT[i, j])

    # Jeder Kunde wird genau beliefert
    for j in range(N_KUNDE):
        index = np.array([i * N_KUNDE + j for i in range(N_LAGER)], dtype=np.int32)
        hochschule.addRow(BEDARF[j], BEDARF[j], len(index), index,
                          np.ones(len(index)))

    # Die Kopplung: sum_j x_ij - M * y_i <= 0
    for i in range(N_LAGER):
        index = np.array([i * N_KUNDE + j for j in range(N_KUNDE)] + [anzahl_x + i],
                         dtype=np.int32)
        werte = np.concatenate([np.ones(N_KUNDE), [-big_m]])
        hochschule.addRow(-unendlich, 0.0, len(index), index, werte)

    t0 = time.perf_counter()
    hochschule.run()
    dauer = time.perf_counter() - t0

    loesung = np.array(hochschule.getSolution().col_value)
    info = hochschule.getInfo()
    return {
        "x": loesung[:anzahl_x].reshape(N_LAGER, N_KUNDE),
        "y": loesung[anzahl_x:],
        "zielwert": info.objective_function_value,
        "knoten": info.mip_node_count,
        "dauer": dauer,
    }


def pruefe(ergebnis: dict) -> tuple[bool, list[str]]:
    """Die Pruefung, die in jedes MILP-Auswertungsskript gehoert.

    Sie rechnet die Kosten AUS DER LOESUNG neu aus, statt dem Zielwert des
    Solvers zu glauben - und vergleicht beide. Genau diese Gegenrechnung
    entlarvt eine Loesung, in der Binaervariablen bei 1e-8 haengengeblieben
    sind.
    """
    beanstandungen = []
    x, y = ergebnis["x"], ergebnis["y"]

    # 1. Sind die Binaervariablen wirklich binaer?
    abstand = np.abs(y - np.round(y))
    if abstand.max() > GANZZAHL_TOLERANZ:
        beanstandungen.append(
            f"y ist nicht ganzzahlig: groesster Abstand {abstand.max():.2e}")

    # 2. Liefert ein Lager, dessen Schalter aus ist?
    liefert = x.sum(axis=1) > 1e-6
    geschlossen_aber_aktiv = np.where(liefert & (y < 0.5))[0]
    if geschlossen_aber_aktiv.size:
        beanstandungen.append(
            f"Lager {geschlossen_aber_aktiv.tolist()} liefern Ware, "
            f"gelten im Modell aber als geschlossen")

    # 3. Stimmt der Zielwert mit den echten Kosten ueberein?
    echte_fixkosten = FIXKOSTEN[liefert].sum()
    echte_transportkosten = float((TRANSPORT * x).sum())
    echte_kosten = echte_fixkosten + echte_transportkosten
    if abs(echte_kosten - ergebnis["zielwert"]) > 1e-4 * max(1.0, echte_kosten):
        beanstandungen.append(
            f"Zielwert {ergebnis['zielwert']:,.2f} weicht von den echten Kosten "
            f"{echte_kosten:,.2f} ab (Differenz {echte_kosten - ergebnis['zielwert']:,.2f})")

    return not beanstandungen, beanstandungen


def zeige(titel: str, ergebnis: dict) -> None:
    x, y = ergebnis["x"], ergebnis["y"]
    liefert = x.sum(axis=1) > 1e-6
    print(f"\n{titel}")
    print(f"  Zielwert laut Solver        {ergebnis['zielwert']:>14,.2f} EUR")
    print(f"  Knoten / Zeit               {ergebnis['knoten']:>14,}  "
          f"/ {ergebnis['dauer']:.3f} s")
    print(f"  Lager mit y = 1             {int((y > 0.5).sum()):>14}")
    print(f"  Lager, die tatsaechlich liefern {int(liefert.sum()):>10}")
    unter = y[y < 0.5]
    print(f"  groesster y-Wert unter 0.5  "
          f"{(unter.max() if unter.size else 0.0):>14.3e}")
    print(f"  Fixkosten real / verbucht   {FIXKOSTEN[liefert].sum():>14,.2f} "
          f"/ {float((FIXKOSTEN * y).sum()):,.2f} EUR")

    in_ordnung, beanstandungen = pruefe(ergebnis)
    if in_ordnung:
        print("  PRUEFUNG: bestanden")
    else:
        print("  PRUEFUNG: DURCHGEFALLEN")
        for text in beanstandungen:
            print(f"    - {text}")


if __name__ == "__main__":
    print("=" * 78)
    print("  DIE BIG-M-FALLE: STANDORTPLANUNG, 12 LAGER, 40 KUNDEN")
    print("=" * 78)
    print(f"Tatsaechliche Lagerkapazitaet: {KAPAZITAET:,.1f} Einheiten.")
    print("Genau das ist das kleinstmoegliche gueltige M - mehr kann ein Lager")
    print("ohnehin nicht ausliefern.")

    knapp = loese(KAPAZITAET)
    zeige("[1] M = Kapazitaet (richtig gewaehlt)", knapp)

    gross = loese(1e7 * KAPAZITAET, presolve="off")
    zeige("[2] M = 10 Millionen mal Kapazitaet, Presolve abgeschaltet", gross)

    print("\n" + "=" * 78)
    print("  WAS DA PASSIERT IST")
    print("=" * 78)
    fehlbetrag = knapp["zielwert"] - gross["zielwert"]
    print(f"Lauf [2] meldet {gross['zielwert']:,.2f} EUR und sieht damit um")
    print(f"{fehlbetrag:,.2f} EUR BESSER aus als die richtige Loesung - ein Ergebnis,")
    print("ueber das sich jeder Auftraggeber freuen wuerde.")
    print()
    print("Der Grund steht in der Zeile 'groesster y-Wert unter 0.5': Die")
    print(f"Schaltervariablen stehen bei rund "
          f"{gross['y'][gross['y'] < 0.5].max():.0e}.")
    print(f"Das ist kleiner als die Ganzzahltoleranz {GANZZAHL_TOLERANZ:.0e}, also gilt")
    print("y = 0 - 'Lager geschlossen'. Zugleich ist M so gross, dass")
    print("  sum_j x_ij <= M * 1e-8")
    print("immer noch reichlich Liefermenge erlaubt. Die Lager liefern also,")
    print("ohne dass ihre Fixkosten je bezahlt werden. Der Fachbegriff dafuer")
    print("ist 'trickle flow'.")
    print()
    print("WICHTIG: Mit eingeschaltetem Presolve (Standard) faellt HiGHS hier")
    print("nicht darauf herein - es zieht M selbst zurecht. Verlassen Sie sich")
    print("nicht darauf: Presolve kann das nur, wenn eine implizite Schranke")
    print("herleitbar ist. Die Pruefung aus pruefe() kostet Millisekunden und")
    print("funktioniert immer.")
    print("=" * 78)