operations_research/OR_HTML_04/programme/Skalierung_Kondition.py

#!/usr/bin/env python3

# Skalierung_Kondition.py
"""
Kapitel Fundament: Was die Konditionszahl kappa(A) praktisch bedeutet - und wie man
ein schlecht skaliertes Modell wieder gesund rechnet.

Drei Experimente:

  1. Fehlerverstaerkung: Wie stark schlaegt eine winzige Datenunsicherheit auf
     die Loesung durch? kappa(A) ist genau die Obergrenze dieses Faktors.
  2. Ruiz-Equilibrierung: Ein Modell, in dem Euro-Betraege (1e7) und
     Tonnen-Angaben (1e-3) in derselben Matrix stehen, wird durch Zeilen- und
     Spaltenskalierung um Groessenordnungen besser konditioniert.
  3. Toleranzen: Warum eine Binaervariable mit dem Wert 0.99999998 niemals
     mit int() gerundet werden darf.

Benoetigt: numpy, scipy
"""

from __future__ import annotations

import numpy as np
from scipy.optimize import linprog

RNG = np.random.default_rng(42)


# --- Experiment 1: kappa(A) als Fehlerverstaerker ---------------------------

def fehlerverstaerkung(A: np.ndarray, versuche: int = 200,
                       stoerung: float = 1e-10) -> tuple[float, float]:
    """Stoert die rechte Seite b relativ um 'stoerung' in zufaellige Richtungen
    und misst, um welchen Faktor sich der Fehler in der Loesung x vergroessert.

    Liefert (Median, Maximum) der Verstaerkung. Die Theorie sagt: Das Maximum
    kann bis kappa(A) betragen - und nur bis dahin.
    """
    x_wahr = np.ones(A.shape[1])
    b = A @ x_wahr
    faktoren = []
    for _ in range(versuche):
        richtung = RNG.normal(size=b.size)
        richtung /= np.linalg.norm(richtung)
        b_gestoert = b + stoerung * np.linalg.norm(b) * richtung
        x_gestoert = np.linalg.solve(A, b_gestoert)
        rel_x = np.linalg.norm(x_gestoert - x_wahr) / np.linalg.norm(x_wahr)
        faktoren.append(rel_x / stoerung)
    return float(np.median(faktoren)), float(np.max(faktoren))


def zeige_experiment_1() -> None:
    print("=" * 74)
    print("  1. KONDITIONSZAHL ALS FEHLERVERSTAERKER")
    print("=" * 74)
    print("Zwei Gleichungssysteme, beide exakt loesbar mit x = (1, 1).")
    print("Die rechte Seite wird um relativ 1e-10 gestoert - so viel Unsicherheit")
    print("steckt in JEDER gemessenen Betriebszahl allemal.\n")

    modelle = {
        "gut konditioniert":     np.array([[2.0, 1.0], [1.0, 3.0]]),
        "fast parallele Zeilen": np.array([[1.0, 1.0], [1.0, 1.0 + 1e-8]]),
    }
    print(f"{'Matrix':<24} {'kappa(A)':>12} {'Verst. median':>14} {'Verst. max':>12}")
    print("-" * 74)
    for name, A in modelle.items():
        median, maximum = fehlerverstaerkung(A)
        print(f"{name:<24} {np.linalg.cond(A):>12.2e} {median:>14.2e} {maximum:>12.2e}")

    print("\nLesart: Bei der zweiten Matrix wird aus einem Datenfehler in der")
    print("10. Nachkommastelle ein Loesungsfehler in der 2. Nachkommastelle.")
    print("Das Modell ist mathematisch korrekt - und praktisch wertlos.")


# --- Experiment 2: Ruiz-Equilibrierung --------------------------------------

def ruiz_equilibrierung(A: np.ndarray, durchlaeufe: int = 20
                        ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Skaliert A iterativ so, dass alle Zeilen- und Spaltenmaxima nahe 1
    liegen (Ruiz 2001).

    In jedem Durchlauf wird jede Zeile durch die Wurzel ihres Betragsmaximums
    geteilt, danach jede Spalte. Das Verfahren konvergiert schnell und braucht
    keinerlei Wissen ueber die Bedeutung der Zahlen - genau deshalb steckt es
    in jedem ernsthaften Solver als Vorverarbeitung.

    Liefert (A_skaliert, zeilenfaktor r, spaltenfaktor c) mit
        A_skaliert = diag(r) @ A @ diag(c)
    """
    m, n = A.shape
    r = np.ones(m)
    c = np.ones(n)
    A_s = A.astype(float).copy()

    for _ in range(durchlaeufe):
        zeilen_max = np.abs(A_s).max(axis=1)
        zeilen_max[zeilen_max == 0] = 1.0
        d_r = 1.0 / np.sqrt(zeilen_max)
        A_s = d_r[:, None] * A_s
        r *= d_r

        spalten_max = np.abs(A_s).max(axis=0)
        spalten_max[spalten_max == 0] = 1.0
        d_c = 1.0 / np.sqrt(spalten_max)
        A_s = A_s * d_c[None, :]
        c *= d_c

    return A_s, r, c


def zeige_experiment_2() -> None:
    print("\n" + "=" * 74)
    print("  2. RUIZ-EQUILIBRIERUNG: EINHEITEN GERADERUECKEN")
    print("=" * 74)
    print("Ein Produktionsmodell, in dem vier Ressourcen in voellig")
    print("verschiedenen Einheiten gemessen werden:")
    print("  Zeile 1: Kapitalbindung in Euro          (Groessenordnung 1e7)")
    print("  Zeile 2: Katalysatorverbrauch in Tonnen  (Groessenordnung 1e-3)")
    print("  Zeile 3: Energie in Wattsekunden         (Groessenordnung 1e5)")
    print("  Zeile 4: Ausschussquote als Anteil       (Groessenordnung 1e-2)\n")

    # Quadratisch gewaehlt, damit die Loesung eindeutig ist und die
    # Ruecktransformation unten wirklich etwas beweist.
    n = 4
    grundmatrix = RNG.uniform(0.5, 2.0, size=(n, n))
    einheiten = np.array([1e7, 1e-3, 1e5, 1e-2])
    A = grundmatrix * einheiten[:, None]

    A_s, r, c = ruiz_equilibrierung(A)

    print(f"{'':<28} {'kappa(A)':>12} {'groesster Eintrag':>18} "
          f"{'kleinster Eintrag':>18}")
    print("-" * 74)
    for name, matrix in [("vor der Skalierung", A), ("nach Ruiz-Equilibrierung", A_s)]:
        betraege = np.abs(matrix)
        print(f"{name:<28} {np.linalg.cond(matrix):>12.2e} "
              f"{betraege.max():>18.2e} {betraege.min():>18.2e}")

    # Gegenprobe: Das skalierte Modell beschreibt dasselbe Problem. Wer x_s
    # loest, erhaelt die urspruengliche Loesung durch x = c * x_s.
    x_wahr = RNG.uniform(1.0, 5.0, size=n)
    b = A @ x_wahr
    b_s = r * b
    x_s = np.linalg.solve(A_s, b_s)
    x_zurueck = c * x_s
    print(f"\nRuecktransformation x = c * x_s: groesste Abweichung zur wahren "
          f"Loesung {np.abs(x_zurueck - x_wahr).max():.2e}")
    print("Die Skalierung ist also verlustfrei - sie aendert nur die Zahlen,")
    print("nicht das Problem.")


# --- Experiment 3: Toleranzen und der int()-Fehler --------------------------

def zeige_experiment_3() -> None:
    print("\n" + "=" * 74)
    print("  3. TOLERANZEN: WARUM int() DIE FALSCHE RUNDUNG IST")
    print("=" * 74)

    # Ein LP, dessen Optimum bei x = 1 liegt, aber vom Solver nur bis auf
    # seine Toleranz getroffen wird.
    ergebnis = linprog(c=[-1.0], A_ub=[[1.0]], b_ub=[1.0],
                       bounds=[(0, None)], method="highs")
    wert = float(ergebnis.x[0])
    print(f"Solver liefert x = {wert!r}")

    # Typische Werte, wie sie aus MILP-Solvern zurueckkommen.
    beispiele = [0.99999998, 1.00000002, 0.49999999, 2.9999999]
    print(f"\n{'Solverwert':>14} {'int()':>8} {'round()':>9} {'Kommentar'}")
    print("-" * 74)
    kommentare = {
        0.99999998: "int() macht aus einer JA- eine NEIN-Entscheidung",
        1.00000002: "hier ginge int() zufaellig gut - Verlass ist keiner",
        0.49999999: "echt unentschieden: Modell oder Toleranz pruefen!",
        2.9999999:  "3 Maschinen werden zu 2 - der Plan geht nicht auf",
    }
    for wert_b in beispiele:
        print(f"{wert_b:>14.8f} {int(wert_b):>8} {round(wert_b):>9} "
              f"{kommentare[wert_b]}")

    print("\nRichtige Vorgehensweise: gegen die Solver-Toleranz pruefen,")
    print("dann erst runden - und den Zweifelsfall melden statt still zu raten.")

    def sichere_ganzzahl(wert: float, toleranz: float = 1e-6) -> int:
        naechste = round(wert)
        if abs(wert - naechste) > toleranz:
            raise ValueError(
                f"{wert} ist {abs(wert - naechste):.2e} von der naechsten ganzen "
                f"Zahl entfernt - das ist mehr als die Toleranz {toleranz}. "
                "Ganzzahligkeit im Modell pruefen.")
        return naechste

    for wert_b in beispiele:
        try:
            print(f"  sichere_ganzzahl({wert_b}) = {sichere_ganzzahl(wert_b)}")
        except ValueError as fehler:
            print(f"  sichere_ganzzahl({wert_b}) -> ValueError: {fehler}")


if __name__ == "__main__":
    zeige_experiment_1()
    zeige_experiment_2()
    zeige_experiment_3()
    print("\n" + "=" * 74)
    print("Merksatz: Skalieren Sie Ihre Daten, BEVOR der Solver sie sieht -")
    print("und runden Sie Solver-Ergebnisse NIE ohne Toleranzpruefung.")
    print("=" * 74)