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