operations_research/OR_HTML_04/programme/Betriebsueberwachung.py

#!/usr/bin/env python3

# Betriebsueberwachung.py
"""
Kapitel Praxisfallen: Das Modell, das still aufgehoert hat zu optimieren.

Ein Optimierungsmodell geht in Betrieb und laeuft jede Nacht. Es liefert
zuverlaessig eine Zahl, die Job-Dauer bleibt konstant, es gibt keine
Fehlermeldung. Alles sieht gut aus.

Waehrenddessen waechst das Unternehmen. Aus 15 Lagern werden 40, dann 80. Das
Zeitlimit von drei Sekunden - vor zwei Jahren grosszuegig bemessen - reicht
irgendwann nicht mehr. Der Solver bricht ab und liefert die beste Loesung, die
er bis dahin gefunden hat.

Das ist nicht falsch. Es ist sogar richtig so. Aber es ist etwas anderes als
das, was das Modell einmal geleistet hat - und niemand merkt es, weil im
Protokoll nur die Kosten stehen.

Dieses Programm zeigt die Verschlechterung und die Ueberwachung, die sie
sichtbar macht.

WICHTIG: ortools wird hier nicht importiert (Konflikt mit highspy, siehe
Kapitel Oekosystem).

Benoetigt: numpy, highspy
"""

from __future__ import annotations

import time
from dataclasses import dataclass

import numpy as np
import highspy

ZEITLIMIT = 3.0                  # das Budget, das der Nachtjob hat
GAP_WARNSCHWELLE = 0.02          # ab 2 % Gap wollen wir es wissen
ZEITAUSLASTUNG_WARNSCHWELLE = 0.9


@dataclass
class Laufprotokoll:
    """Was nach JEDEM Produktivlauf protokolliert gehoert.

    Die meisten Systeme schreiben nur 'kosten'. Genau deshalb faellt eine
    schleichende Verschlechterung ueber Monate nicht auf.
    """
    zeitpunkt: str
    lager: int
    kunden: int
    status: str
    kosten: float
    gap: float
    laufzeit: float
    knoten: int

    @property
    def warnungen(self) -> list[str]:
        meldungen = []
        if self.status != "Optimal":
            meldungen.append(f"nicht beweisbar optimal ({self.status})")
        if self.gap > GAP_WARNSCHWELLE:
            meldungen.append(f"Gap {self.gap:.1%} ueber Schwelle "
                             f"{GAP_WARNSCHWELLE:.0%}")
        if self.laufzeit > ZEITAUSLASTUNG_WARNSCHWELLE * ZEITLIMIT:
            meldungen.append(f"Zeitbudget zu {self.laufzeit / ZEITLIMIT:.0%} "
                             f"ausgeschoepft")
        return meldungen


def plane_netzwerk(n_lager: int, n_kunden: int, zeitlimit: float,
                   saat: int = 7) -> tuple[str, float, float, float, int]:
    """Standortplanung wie in Kapitel MILP, nur mit wachsender Groesse."""
    rng = np.random.default_rng(saat)
    fixkosten = rng.uniform(3000, 9000, n_lager)
    transport = rng.uniform(5, 60, (n_lager, n_kunden))
    bedarf = rng.uniform(10, 60, n_kunden)
    kapazitaet = np.full(n_lager, bedarf.sum() * 0.22)

    modell = highspy.Highs()
    modell.setOptionValue("output_flag", False)
    modell.setOptionValue("time_limit", zeitlimit)

    anzahl_x = n_lager * n_kunden
    unendlich = highspy.kHighsInf
    modell.addVars(anzahl_x, np.zeros(anzahl_x), np.full(anzahl_x, unendlich))
    modell.addVars(n_lager, np.zeros(n_lager), np.ones(n_lager))
    for i in range(n_lager):
        modell.changeColIntegrality(anzahl_x + i, highspy.HighsVarType.kInteger)
        modell.changeColCost(anzahl_x + i, fixkosten[i])
        for j in range(n_kunden):
            modell.changeColCost(i * n_kunden + j, transport[i, j])

    for j in range(n_kunden):
        index = np.array([i * n_kunden + j for i in range(n_lager)], dtype=np.int32)
        modell.addRow(bedarf[j], bedarf[j], len(index), index, np.ones(len(index)))
    for i in range(n_lager):
        index = np.array([i * n_kunden + j for j in range(n_kunden)] + [anzahl_x + i],
                         dtype=np.int32)
        werte = np.concatenate([np.ones(n_kunden), [-kapazitaet[i]]])
        modell.addRow(-unendlich, 0.0, len(index), index, werte)

    t0 = time.perf_counter()
    modell.run()
    laufzeit = time.perf_counter() - t0
    info = modell.getInfo()
    return (modell.modelStatusToString(modell.getModelStatus()),
            info.objective_function_value, info.mip_gap, laufzeit,
            info.mip_node_count)


if __name__ == "__main__":
    # So ist das Unternehmen ueber zwei Jahre gewachsen.
    entwicklung = [("2024-Q1", 15, 40), ("2024-Q3", 25, 70),
                   ("2025-Q1", 40, 110), ("2025-Q3", 60, 160),
                   ("2026-Q1", 80, 220)]

    print("=" * 84)
    print("  WAS DAS PROTOKOLL ZEIGT - UND WAS ES ZEIGEN SOLLTE")
    print("=" * 84)
    print(f"Derselbe Nachtjob, unveraendert, mit {ZEITLIMIT:.0f} Sekunden Zeitlimit.\n")

    protokolle = []
    for zeitpunkt, lager, kunden in entwicklung:
        status, kosten, gap, laufzeit, knoten = plane_netzwerk(
            lager, kunden, ZEITLIMIT)
        protokolle.append(Laufprotokoll(zeitpunkt, lager, kunden, status,
                                        kosten, gap, laufzeit, knoten))

    print("So sieht das ueblich gefuehrte Protokoll aus:\n")
    print(f"  {'Zeitpunkt':<10} {'Kosten':>12}")
    print("  " + "-" * 24)
    for p in protokolle:
        print(f"  {p.zeitpunkt:<10} {p.kosten:>12,.0f}")
    print("\n  -> Die Kosten steigen. Das Unternehmen waechst ja auch.")
    print("     Nichts an dieser Tabelle deutet auf ein Problem hin.")

    print("\n" + "-" * 84)
    print("So sieht ein vollstaendiges Protokoll aus:\n")
    print(f"  {'Zeitpunkt':<10} {'Groesse':>10} {'Status':>20} {'Kosten':>11} "
          f"{'Gap':>8} {'Zeit':>7}")
    print("  " + "-" * 80)
    for p in protokolle:
        groesse = f"{p.lager}x{p.kunden}"
        print(f"  {p.zeitpunkt:<10} {groesse:>10} {p.status:>20} "
              f"{p.kosten:>11,.0f} {p.gap * 100:7.2f} % {p.laufzeit:6.2f}s")

    print("\n" + "-" * 84)
    print("Und so sehen die Warnungen aus, die daraus folgen:\n")
    for p in protokolle:
        if p.warnungen:
            print(f"  {p.zeitpunkt}:  " + "; ".join(p.warnungen))
        else:
            print(f"  {p.zeitpunkt}:  in Ordnung")

    letzte = protokolle[-1]
    erste_warnung = next(p for p in protokolle if p.warnungen)

    print("\n" + "=" * 84)
    print("  WAS DA PASSIERT IST")
    print("=" * 84)
    print(f"Seit {erste_warnung.zeitpunkt} erreicht der Job das Zeitlimit und liefert")
    print("nicht mehr das Optimum, sondern die beste bis dahin gefundene Loesung.")
    print(f"Im letzten Lauf betraegt der Abstand zum Bestmoeglichen "
          f"{letzte.gap:.1%}.")
    print()
    print("Der Job ist nicht abgestuerzt. Er hat keine Fehlermeldung erzeugt.")
    print("Die Laufzeit ist sogar bemerkenswert STABIL geblieben - genau deshalb,")
    print("weil das Zeitlimit greift. Ein Ueberwachungssystem, das auf Abstuerze")
    print("und Laufzeitspitzen achtet, sieht hier nichts.")
    print()
    print("Drei Zahlen gehoeren deshalb in jedes Protokoll eines Optimierungsjobs:")
    print("  * der STATUS   - 'Optimal' oder etwas anderes?")
    print("  * der GAP      - wie weit ist die Loesung vom Bestmoeglichen entfernt?")
    print("  * die LAUFZEIT im Verhaeltnis zum Limit - wie nah am Anschlag?")
    print()
    print("Die dritte ist die frueheste Warnung: Sie steigt, lange bevor der Gap")
    print("sichtbar wird, und gibt Zeit zum Handeln.")
    print("=" * 84)