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