Kapitel 2: Das mathematische Fundament — Vektoren, Matrizen, Konvexität

2.1 In 5 Minuten gelöst

Dieses Kapitel handelt von Geometrie. Bevor wir Vektoren und Matrizen einführen, sehen Sie die zentrale Einsicht in Aktion — ganz ohne Solver, mit acht Zeilen NumPy.

Das Optimum lautet (0, 8) mit Z = 40 — und Sie haben es gefunden, ohne einen einzigen Solver zu starten. Vier Kandidaten statt unendlich vieler Punkte: Das ist der Fundamentalsatz der linearen Optimierung, und er ist der Grund, warum sich lineare Programme überhaupt zuverlässig lösen lassen.

Warum funktioniert das? Weil der zulässige Bereich ein konvexes Polyeder ist — ein Vielflächner ohne Diagonalen nach innen — und weil eine lineare Zielfunktion darauf ihren größten Wert immer am Rand annimmt, genauer: in einer Ecke. Beides begründet dieses Kapitel. Der Simplex-Algorithmus in Kapitel 5 tut im Kern nichts anderes als das Programm oben, nur klüger: Er probiert nicht alle Ecken durch, sondern läuft gezielt von Ecke zu Ecke bergauf.

⚠️ Und wo ist der Haken? Bei zwei Variablen gibt es 6 Eckenkandidaten, bei 50 Variablen und 50 Bedingungen sind es \binom{100}{50} \approx 10^{29}. Der Fundamentalsatz sagt uns wo wir suchen müssen — nicht, dass die Suche leicht wird.

2.2 Lernziele

Nach diesem Kapitel können Sie …

1. … ein Optimierungsmodell von der „Zeile-für-Zeile“-Schreibweise in die kompakte Matrixform \max \mathbf{c}^\top\mathbf{x} u. d. N. \mathbf{A}\mathbf{x} \le \mathbf{b} übersetzen — und zurück.
2. … erklären, warum der zulässige Bereich eines linearen Programms ein Polyeder ist.
3. … den Fundamentalsatz der linearen Optimierung anwenden, um Kandidaten für das Optimum zu finden, ohne alles durchzuprobieren.
4. … prüfen, ob eine Menge bzw. eine Funktion konvex ist, und begründen, warum das für die Lösbarkeit entscheidend ist.
5. … einen zweidimensionalen Lösungsraum mit matplotlib zeichnen und daraus die optimale Ecke ablesen.
6. … die Konditionszahl \kappa(\mathbf{A}) berechnen, ihren Wert deuten und ein schlecht skaliertes Modell mit Ruiz-Equilibrierung wieder rechenbar machen.
7. … erklären, warum ein Solver-Ergebnis von 0.99999998 niemals mit int() in eine ganze Zahl verwandelt werden darf.

2.3 Warum überhaupt Vektoren und Matrizen?

In Kapitel 1 hatten wir zwei Variablen und drei Nebenbedingungen. Das ließ sich bequem ausschreiben. Reale Modelle haben hunderte bis hunderttausende Variablen — dort wäre Ausschreiben nicht nur unpraktisch, sondern unmöglich.

Die lineare Algebra löst dieses Problem, indem sie viele gleichartige Zahlen zu einem Objekt bündelt. Das ist derselbe Gedanke wie eine Liste in Python: Statt preis_1, preis_2, …, preis_1000 schreibt man preise und arbeitet mit dem Ganzen.

Die Bausteine

Sei n die Anzahl der Entscheidungsvariablen und m die Anzahl der Nebenbedingungen.

Entscheidungsvektor — was wir festlegen:

\mathbf{x} = \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix} \in \mathbb{R}^n

Kosten- bzw. Ertragsvektor — was jede Einheit wert ist:

\mathbf{c} = \begin{pmatrix} c_1 \\ c_2 \\ \vdots \\ c_n \end{pmatrix} \in \mathbb{R}^n

Zielfunktion als Skalarprodukt:

f(\mathbf{x}) = \mathbf{c}^\top \mathbf{x} = \sum_{j=1}^n c_j x_j

📐 Formel-Lesehilfe * \mathbf{c}^\top — der Ertragsvektor, „umgekippt“ zu einer Zeile. * \mathbf{c}^\top \mathbf{x} — Skalarprodukt: Multipliziere jedes c_j mit dem zugehörigen x_j und addiere alles. Das Ergebnis ist eine einzige Zahl. * \sum_{j=1}^n — „addiere für j = 1 bis j = n“.

Ohne Formel gesagt: „Nimm von jedem Produkt die hergestellte Menge mal den Gewinn pro Stück und zähle alles zusammen.“ Im Bot-Beispiel: \mathbf{c} = (150, 250)^\top, \mathbf{x} = (7, 5)^\top, also \mathbf{c}^\top\mathbf{x} = 150 \cdot 7 + 250 \cdot 5 = 2300.

Technologiematrix — wer verbraucht wie viel wovon:

\mathbf{A} = \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} \end{pmatrix} \in \mathbb{R}^{m \times n}

📐 Formel-Lesehilfe Jede Zeile ist eine Ressource, jede Spalte eine Variable. Der Eintrag a_{ij} beantwortet: „Wie viel von Ressource i verbraucht eine Einheit von Variable j?“

Im Bot-Beispiel: \mathbf{A} = \begin{pmatrix} 2 & 5 \\ 4 & 6 \\ 1 & 0 \end{pmatrix} \begin{matrix} \leftarrow \text{vCPU} \\ \leftarrow \text{RAM} \\ \leftarrow \text{Liquidität} \end{matrix} Die erste Zeile (2, 5) heißt: Ein A-Bot braucht 2 vCPUs, ein B-Bot braucht 5. Die dritte Zeile (1, 0) heißt: Die Liquiditätsgrenze zählt nur A-Bots, B-Bots gar nicht.

Kapazitätsvektor (rechte Seite):

\mathbf{b} = \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_m \end{pmatrix} \in \mathbb{R}^m \qquad \text{im Beispiel } \begin{pmatrix} 40 \\ 60 \\ 8 \end{pmatrix}

Die kanonische Standardform

Damit lässt sich jedes lineare Programm in drei Zeilen schreiben:

\min_{\mathbf{x}} \; \mathbf{c}^\top \mathbf{x} \qquad \text{u. d. N.} \qquad \mathbf{A}\mathbf{x} \le \mathbf{b}, \qquad \mathbf{x} \ge \mathbf{0}

📐 Formel-Lesehilfe * \mathbf{A}\mathbf{x} — Matrix-Vektor-Produkt: berechnet auf einen Schlag den Verbrauch aller Ressourcen. Zeile i des Ergebnisses ist a_{i1}x_1 + a_{i2}x_2 + \dots + a_{in}x_n. * \mathbf{A}\mathbf{x} \le \mathbf{b} — die Ungleichung gilt komponentenweise: jede einzelne Zeile muss ihre Kapazität einhalten. * \mathbf{x} \ge \mathbf{0} — alle Variablen sind nichtnegativ.

Ohne Formel gesagt: „Minimiere die Gesamtkosten, ohne bei irgendeiner Ressource über die Kapazität zu gehen, und ohne negative Mengen zu produzieren.“

⚠️ Typische Fehler

• Maximieren statt minimieren. Die Standardform minimiert. Eine Maximierung wird durch Vorzeichenwechsel überführt: \max\ \mathbf{c}^\top\mathbf{x} ist dasselbe wie -\min\ (-\mathbf{c})^\top\mathbf{x}. Wer das vergisst, erhält systematisch die schlechteste statt der besten Lösung — siehe Abschnitt 5.7.
• „\ge“-Bedingungen direkt einsetzen. Die Standardform kennt nur „\le“. Aus 3x_1 + 2x_2 \ge 12 wird durch Multiplikation mit -1: -3x_1 - 2x_2 \le -12. Beim Multiplizieren mit einer negativen Zahl dreht sich das Ungleichheitszeichen um — der häufigste Vorzeichenfehler überhaupt.
• Gleichungen vergessen. h(\mathbf{x}) = b lässt sich als zwei Ungleichungen schreiben: h(\mathbf{x}) \le b und -h(\mathbf{x}) \le -b. Die meisten Solver nehmen Gleichungen aber direkt entgegen — man muss es nicht von Hand machen.

✏️ Handrechnung 2.1: Modell in Matrixform übersetzen

Gegeben: \max\ 3x_1 + 5x_2 \quad \text{u. d. N.}\quad x_1 \le 4,\quad 2x_2 \le 12,\quad 3x_1 + 2x_2 \le 18,\quad x_1, x_2 \ge 0

Schritt 1 — Vektoren und Matrix ablesen: \mathbf{c} = \begin{pmatrix}3\\5\end{pmatrix},\quad \mathbf{A} = \begin{pmatrix}1 & 0\\ 0 & 2\\ 3 & 2\end{pmatrix},\quad \mathbf{b} = \begin{pmatrix}4\\12\\18\end{pmatrix}

Achten Sie auf die Nullen: Die erste Bedingung enthält x_2 gar nicht, also steht dort eine 0. Das Weglassen der Nullen ist der häufigste Anfängerfehler beim Aufstellen von \mathbf{A}.

Schritt 2 — Für die Standardform (Minimierung) negieren: \tilde{\mathbf{c}} = (-3, -5)^\top. Der optimale Zielwert der Minimierung ist dann das Negative des gesuchten Maximums.

Schritt 3 — Probe mit \mathbf{x} = (2,6)^\top: \mathbf{A}\mathbf{x} = \begin{pmatrix}1\cdot2 + 0\cdot6\\ 0\cdot2+2\cdot6\\ 3\cdot2+2\cdot6\end{pmatrix} = \begin{pmatrix}2\\12\\18\end{pmatrix} \le \begin{pmatrix}4\\12\\18\end{pmatrix} \ \checkmark Alle drei Zeilen halten. Zielwert: 3\cdot2 + 5\cdot6 = 36.

Dasselbe in Python

#!/usr/bin/env python3

# Matrixform.py
"""
Kapitel Fundament: Von der ausgeschriebenen Form zur Matrixform - und zurück.
Zeigt, dass beide Schreibweisen dasselbe Modell beschreiben.
"""

import numpy as np
from scipy.optimize import linprog

# --- Modell in Matrixform -------------------------------------------------
# max 3*x1 + 5*x2   u.d.N.   x1 <= 4,  2*x2 <= 12,  3*x1 + 2*x2 <= 18,  x >= 0
c = np.array([3.0, 5.0])                 # Ertragsvektor (Maximierung)
A = np.array([[1.0, 0.0],                # Zeile 1: nur x1 kommt vor
              [0.0, 2.0],                # Zeile 2: nur x2 kommt vor
              [3.0, 2.0]])               # Zeile 3: beide
b = np.array([4.0, 12.0, 18.0])
namen = ["Rohstoff A", "Rohstoff B", "Maschinenzeit"]

# --- Ausgeschriebene Form maschinell erzeugen ------------------------------
def zeige_ausgeschrieben(c, A, b, namen):
    """Druckt die Matrixform als lesbares Ungleichungssystem."""
    terme = " + ".join(f"{c[j]:g}*x{j+1}" for j in range(len(c)))
    print(f"max  {terme}")
    print("u.d.N.")
    for i in range(A.shape[0]):
        summanden = " + ".join(f"{A[i, j]:g}*x{j+1}"
                               for j in range(A.shape[1]) if A[i, j] != 0)
        print(f"   {summanden:<24} <= {b[i]:>5g}   ({namen[i]})")
    print(f"   x1, ..., x{len(c)} >= 0")

zeige_ausgeschrieben(c, A, b, namen)

# --- Zulässigkeit eines Punktes prüfen ------------------------------------
def ist_zulaessig(x, A, b, toleranz=1e-9):
    """Prüft A x <= b und x >= 0 komponentenweise."""
    verbrauch = A @ x                     # Matrix-Vektor-Produkt: alle Zeilen auf einmal
    return bool(np.all(verbrauch <= b + toleranz) and np.all(x >= -toleranz))

for kandidat in [np.array([2.0, 6.0]), np.array([4.0, 3.0]), np.array([4.0, 6.0])]:
    zulaessig = ist_zulaessig(kandidat, A, b)
    zielwert = c @ kandidat
    verbrauch = A @ kandidat
    print(f"\nx = {kandidat}  ->  A x = {verbrauch}  "
          f"{'zulaessig' if zulaessig else 'UNZULAESSIG'},  Z = {zielwert:g}")

# --- Lösen: linprog minimiert, also c negieren -----------------------------
ergebnis = linprog(c=-c, A_ub=A, b_ub=b, bounds=[(0, None)] * len(c), method="highs")
print("\n" + "-" * 60)
print(f"Optimale Loesung:  x* = {np.round(ergebnis.x, 4)}")
print(f"Optimaler Wert:    Z* = {-ergebnis.fun:g}")
print("Hinweis: linprog minimiert, deshalb wurde c negiert und das")
print("         Ergebnis am Ende wieder mit -1 multipliziert.")

Erwartete Ausgabe:

max  3*x1 + 5*x2
u.d.N.
   1*x1                     <=     4   (Rohstoff A)
   2*x2                     <=    12   (Rohstoff B)
   3*x1 + 2*x2              <=    18   (Maschinenzeit)
   x1, ..., x2 >= 0

x = [2. 6.]  ->  A x = [ 2. 12. 18.]  zulaessig,  Z = 36

x = [4. 3.]  ->  A x = [ 4.  6. 18.]  zulaessig,  Z = 27

x = [4. 6.]  ->  A x = [ 4. 12. 24.]  UNZULAESSIG,  Z = 42

------------------------------------------------------------
Optimale Loesung:  x* = [2. 6.]
Optimaler Wert:    Z* = 36
Hinweis: linprog minimiert, deshalb wurde c negiert und das
         Ergebnis am Ende wieder mit -1 multipliziert.

Der dritte Kandidat (4, 6) zeigt die Falle: Er hätte mit Z = 42 den höchsten Zielwert — ist aber unzulässig, weil die Maschinenzeit mit 24 über der Kapazität von 18 liegt. Ein hoher Zielwert allein bedeutet nichts.

2.4 Der zulässige Lösungsraum und das Polyeder

Jede lineare Ungleichung a_{i1}x_1 + a_{i2}x_2 \le b_i definiert geometrisch eine Halbebene (im \mathbb{R}^2) bzw. einen Halbraum (im \mathbb{R}^n): Die zugehörige Gleichung ist eine Gerade (bzw. Hyperebene), und die Ungleichung wählt eine der beiden Seiten aus.

Der Schnitt aller Halbräume bildet den zulässigen Bereich \mathcal{F}:

\mathcal{F} = \{ \mathbf{x} \in \mathbb{R}^n \mid \mathbf{A}\mathbf{x} \le \mathbf{b},\; \mathbf{x} \ge \mathbf{0} \}

📐 Formel-Lesehilfe Die geschweiften Klammern beschreiben eine Menge. Der senkrechte Strich \mid heißt „für die gilt“.

Ohne Formel gesagt: „\mathcal{F} ist die Menge aller Punkte, die gleichzeitig alle Kapazitätsgrenzen einhalten und nicht negativ sind.“ Anschaulich: der Bereich, in dem man überhaupt landen darf.

Geometrisch ist \mathcal{F} ein konvexes Polyeder — ein von ebenen Flächen begrenzter Körper ohne Diese-Eindellungen. Im Zweidimensionalen ein Vieleck, im Dreidimensionalen ein Körper wie ein geschliffener Diamant.

Der Fundamentalsatz der linearen Optimierung

Satz. Besitzt ein lineares Optimierungsproblem eine optimale Lösung, dann liegt mindestens ein optimaler Punkt auf einer Ecke (Extrempunkt) des Polyeders \mathcal{F}.

Warum das gilt — die Anschauung: Die Zielfunktion \mathbf{c}^\top\mathbf{x} hat überall dieselbe Steigungsrichtung; ihre Höhenlinien sind parallele Geraden. Stellen Sie sich vor, Sie schieben eine solche Gerade in Richtung wachsender Zielwerte über das Polyeder. Der letzte Punkt, den sie berührt, bevor sie das Polyeder verlässt, ist eine Ecke — oder, wenn die Gerade zufällig parallel zu einer Kante liegt, eine ganze Kante, deren Endpunkte wiederum Ecken sind.

Warum das praktisch so wichtig ist: Das Polyeder enthält unendlich viele Punkte, aber nur endlich viele Ecken. Der Simplex-Algorithmus (Kapitel 5) muss deshalb nicht unendlich viel absuchen, sondern wandert gezielt von Ecke zu Ecke.

✏️ Handrechnung 2.2: Optimum über Ecken finden

Wir nehmen wieder \max 3x_1 + 5x_2 mit x_1 \le 4, 2x_2 \le 12, 3x_1 + 2x_2 \le 18, x_1, x_2 \ge 0.

Schritt 1 — Ecken bestimmen. Eine Ecke entsteht dort, wo sich zwei Begrenzungsgeraden schneiden (und der Punkt zulässig ist). Die Geraden sind: x_1 = 0, x_2 = 0, x_1 = 4, x_2 = 6, 3x_1 + 2x_2 = 18.

Schnittpunkt von	Ecke	zulässig?	Z = 3x_1 + 5x_2
x_1=0, x_2=0	(0,0)	✓	0
x_1=4, x_2=0	(4,0)	✓	12
x_1=4, 3x_1+2x_2=18	(4,3)	✓	27
x_2=6, 3x_1+2x_2=18	(2,6)	✓	36 ← Optimum
x_1=0, x_2=6	(0,6)	✓	30
x_1=4, x_2=6	(4,6)	✗ (3\cdot4+2\cdot6=24>18)	—

Schritt 2 — Ergebnis ablesen: \mathbf{x}^* = (2,6)^\top, Z^* = 36.

Beachten Sie: Wir haben fünf Punkte geprüft statt unendlich viele. Genau das ist der Gewinn des Fundamentalsatzes. Für n = 2 geht das von Hand; bei n = 50 gibt es zu viele Ecken — dann übernimmt der Simplex-Algorithmus, der nur die verbessernden Ecken besucht.

2.5 Konvexität: die Grenze zwischen leicht und schwer

Hier kommt die vielleicht wichtigste Einsicht des ganzen Buches. In der Optimierungstheorie verläuft die Trennlinie zwischen „zuverlässig lösbar“ und „im Allgemeinen hoffnungslos“ nicht zwischen linear und nichtlinear — sondern zwischen konvex und nicht-konvex.

Konvexe Menge

Eine Menge \mathcal{C} \subseteq \mathbb{R}^n heißt konvex, wenn für alle Punkte \mathbf{x}, \mathbf{y} \in \mathcal{C} und jedes \theta \in [0,1] gilt:

\theta \mathbf{x} + (1-\theta)\mathbf{y} \in \mathcal{C}

📐 Formel-Lesehilfe * \theta („theta“) ist eine Zahl zwischen 0 und 1. * \theta\mathbf{x} + (1-\theta)\mathbf{y} durchläuft für \theta von 1 bis 0 genau die Verbindungsstrecke von \mathbf{x} nach \mathbf{y}. Bei \theta=1 sind wir in \mathbf{x}, bei \theta=0 in \mathbf{y}, bei \theta=0{,}5 genau in der Mitte.

Ohne Formel gesagt: Eine Menge ist konvex, wenn man zwischen zwei beliebigen ihrer Punkte eine gerade Linie ziehen kann, ohne die Menge zu verlassen. Ein Kreis, ein Quadrat und jedes Polyeder sind konvex. Ein Halbmond, ein Ring und ein Stern sind es nicht.

Konvexe Funktion

Eine Funktion f: \mathcal{C} \to \mathbb{R} heißt konvex, wenn für alle \mathbf{x}, \mathbf{y} \in \mathcal{C} und \theta \in [0,1] gilt:

f(\theta\mathbf{x} + (1-\theta)\mathbf{y}) \le \theta f(\mathbf{x}) + (1-\theta)f(\mathbf{y})

📐 Formel-Lesehilfe * Links: der Funktionswert auf der Verbindungsstrecke. * Rechts: der entsprechende Wert auf der Sehne zwischen den beiden Funktionswerten.

Ohne Formel gesagt: Der Funktionsgraph liegt zwischen je zwei Punkten immer unterhalb der geraden Verbindungslinie. Eine konvexe Funktion ist „nach oben offen geschüsselt“ — wie eine Parabel x^2 oder eine Suppenschüssel. Kippen Sie eine Kugel hinein, rollt sie zum tiefsten Punkt, und der ist eindeutig.

Der zentrale Satz

Satz. Bei einem konvexen Optimierungsproblem (konvexe Zielfunktion über konvexem zulässigen Bereich) ist jedes lokale Minimum automatisch auch das globale Minimum.

Warum das den Unterschied macht: Ein Algorithmus, der nur lokal sucht — „gehe bergab, bis es nicht mehr bergab geht“ — findet bei konvexen Problemen garantiert das absolut beste Ergebnis. Bei nicht-konvexen Problemen bleibt derselbe Algorithmus im nächstbesten Tal stecken, und niemand kann ihm ansehen, ob es das tiefste war. Man müsste alle Täler prüfen — und deren Zahl wächst wieder kombinatorisch.

🎯 Merksatz Konvex heißt: Wer bergab geht, kommt am tiefsten Punkt an. Nicht-konvex heißt: Wer bergab geht, kommt irgendwo an — und weiß nicht, ob es der tiefste Punkt war.

Was Sie in der Abbildung sehen. Links sind es fünf Läufe desselben Verfahrens, die in vier verschiedenen Tälern enden — dieselben vier, die Teil 2 von Konvexitaet_Demo.py weiter unten auflistet. In zwei Dimensionen werden aus den Tälern Mulden in einem Gitter, und der Unterschied zwischen Start A und Start B beträgt hier bereits mehr als 24 Einheiten. Keiner der beiden Läufe hat etwas falsch gemacht: Beide enden dort, wo es nicht mehr bergab geht.

Was ist konvex, was nicht?

Konvexität sichtbar machen

#!/usr/bin/env python3

# Konvexitaet_Demo.py
"""
Kapitel Fundament: Konvexität praktisch erfahrbar machen.
Teil 1: Die Sehnen-Bedingung numerisch nachprüfen.
Teil 2: Zeigen, warum lokale Suche bei nicht-konvexen Funktionen scheitert.
"""

import numpy as np
from scipy.optimize import minimize

# --- Teil 1: Sehnen-Test ---------------------------------------------------
def ist_konvex_numerisch(f, unten, oben, proben=2000, seed=0):
    """
    Prüft die Konvexitätsdefinition an zufälligen Punktepaaren:
       f(theta*x + (1-theta)*y)  <=  theta*f(x) + (1-theta)*f(y)  ?
    Findet Gegenbeispiele - beweist aber KEINE Konvexität.
    """
    rng = np.random.default_rng(seed)
    schlimmste_verletzung = 0.0
    for _ in range(proben):
        x, y = rng.uniform(unten, oben, size=2)
        theta = rng.uniform(0.0, 1.0)
        links = f(theta * x + (1 - theta) * y)
        rechts = theta * f(x) + (1 - theta) * f(y)
        schlimmste_verletzung = max(schlimmste_verletzung, links - rechts)
    return schlimmste_verletzung


funktionen = {
    "f(x) = x^2            (konvex)":            lambda x: x ** 2,
    "f(x) = |x|            (konvex)":            lambda x: abs(x),
    "f(x) = e^x            (konvex)":            lambda x: np.exp(x),
    "f(x) = x^3            (NICHT konvex)":      lambda x: x ** 3,
    "f(x) = x^2+3sin(3x)   (NICHT konvex)":      lambda x: x ** 2 + 3 * np.sin(3 * x),
}

print("=" * 68)
print("  TEIL 1: SEHNEN-TEST  (positive Zahl = Konvexität verletzt)")
print("=" * 68)
for name, f in funktionen.items():
    verletzung = ist_konvex_numerisch(f, -3.0, 3.0)
    urteil = "konvex (keine Verletzung gefunden)" if verletzung < 1e-9 \
             else f"NICHT konvex (Verletzung bis {verletzung:.3f})"
    print(f"{name:<34} -> {urteil}")

# --- Teil 2: Lokale Suche von verschiedenen Startpunkten -------------------
def wellige_funktion(x):
    """Nicht-konvex: eine Parabel mit aufmodulierter Welle -> viele lokale Minima."""
    return x[0] ** 2 + 3.0 * np.sin(3.0 * x[0])

print("\n" + "=" * 68)
print("  TEIL 2: LOKALE SUCHE BEI NICHT-KONVEXER FUNKTION")
print("=" * 68)
print(f"{'Startpunkt':>12} | {'gefundenes Minimum':>20} | {'Funktionswert':>14}")
print("-" * 68)

ergebnisse = []
for start in [-3.0, -1.5, 0.0, 1.5, 3.0]:
    res = minimize(wellige_funktion, x0=[start], method="BFGS")
    ergebnisse.append((start, res.x[0], res.fun))
    print(f"{start:>12.1f} | {res.x[0]:>20.4f} | {res.fun:>14.4f}")

bester = min(ergebnisse, key=lambda t: t[2])
print("-" * 68)
print(f"Je nach Startpunkt landet derselbe Algorithmus in "
      f"{len({round(e[1], 3) for e in ergebnisse})} verschiedenen Minima.")
print(f"Das beste gefundene: x = {bester[1]:.4f} mit f = {bester[2]:.4f} "
      f"(Start bei {bester[0]:.1f})")
print("Bei einer KONVEXEN Funktion waeren alle Zeilen identisch --")
print("der Startpunkt waere voellig gleichgueltig.")
print("=" * 68)

Erwartete Ausgabe:

====================================================================
  TEIL 1: SEHNEN-TEST  (positive Zahl = Konvexität verletzt)
====================================================================
f(x) = x^2            (konvex)      -> konvex (keine Verletzung gefunden)
f(x) = |x|            (konvex)      -> konvex (keine Verletzung gefunden)
f(x) = e^x            (konvex)      -> konvex (keine Verletzung gefunden)
f(x) = x^3            (NICHT konvex) -> NICHT konvex (Verletzung bis 12.933)
f(x) = x^2+3sin(3x)   (NICHT konvex) -> NICHT konvex (Verletzung bis 4.800)

====================================================================
  TEIL 2: LOKALE SUCHE BEI NICHT-KONVEXER FUNKTION
====================================================================
  Startpunkt |   gefundenes Minimum |  Funktionswert
--------------------------------------------------------------------
        -3.0 |              -2.4280 |         3.3695
        -1.5 |              -0.4874 |        -2.7448
         0.0 |              -0.4874 |        -2.7448
         1.5 |               1.4606 |        -0.7042
         3.0 |               3.3817 |         9.4570
--------------------------------------------------------------------
Je nach Startpunkt landet derselbe Algorithmus in 4 verschiedenen Minima.
Das beste gefundene: x = -0.4874 mit f = -2.7448 (Start bei -1.5)
Bei einer KONVEXEN Funktion waeren alle Zeilen identisch --
der Startpunkt waere voellig gleichgueltig.
====================================================================

💻 Code-Durchgang

Stelle	Was passiert	Warum wichtig
ist_konvex_numerisch	testet die Definition an Zufallspaaren	Ein einziges Gegenbeispiel widerlegt Konvexität; kein noch so großer Test beweist sie
links - rechts	Abstand zur Sehne	Positiv = Graph liegt über der Sehne = Verletzung
Teil 2, Schleife über Startpunkte	dieselbe Funktion, fünf Starts	Zeigt die praktische Konsequenz: das Ergebnis hängt vom Zufall des Startpunkts ab
method="BFGS"	klassisches lokales Gradientenverfahren	Genau der Verfahrenstyp, den scipy.optimize in Kapitel 11 verwendet

Die Lektion in einer Zeile: Wenn Ihr Optimierungsergebnis vom Startwert abhängt, ist Ihr Problem nicht konvex — und Sie haben keine Garantie, das Beste gefunden zu haben.

Achten Sie besonders auf die erste Zeile der Tabelle: Vom Startpunkt -3{,}0 aus landet das Verfahren bei f = 3{,}37 — einem Wert, der um mehr als 6 Einheiten schlechter ist als das beste gefundene Minimum. Der Algorithmus meldet dabei keinen Fehler und keine Warnung. Er hat korrekt gearbeitet, ein lokales Minimum gefunden und ist stehen geblieben. Genau das ist die Gefahr: Nicht-Konvexität erzeugt keine Fehlermeldungen, sondern stille Fehlentscheidungen.

⚠️ Typische Fehler

• Konvexität mit Linearität verwechseln. f(x) = x^2 ist nichtlinear, aber konvex und damit bequem lösbar. f(x) = x^3 ist ebenfalls nichtlinear, aber nicht konvex und deshalb problematisch. Nichtlinear ist kein Urteil, konvex schon.
• Konvexität nur der Zielfunktion prüfen. Auch der zulässige Bereich muss konvex sein. Ganzzahligkeitsbedingungen zerstören genau diese Eigenschaft — deshalb ist MILP schwer, obwohl alles daran linear ist.
• Eine Matrix ungeprüft als Kovarianzmatrix verwenden. Nur wenn alle Eigenwerte \ge 0 sind, ist \mathbf{x}^\top\mathbf{P}\mathbf{x} konvex. Ein einziger negativer Eigenwert kippt das ganze Problem — mehr dazu in Abschnitt 11.5.

2.6 Geometrische Visualisierung des Lösungsraums

Wir zeichnen nun den Lösungsraum aus Kapitel 1 und lesen das Optimum ab. Zweidimensionale Bilder sind der schnellste Weg, Intuition für höhere Dimensionen aufzubauen — auch wenn man sie dort nicht mehr zeichnen kann.

#!/usr/bin/env python3

# Visualisierung_Loesungsraum.py
"""
Kapitel Fundament: Geometrische Visualisierung eines 2D-Optimierungsraums.

Ecken werden berechnet, auf Zulässigkeit geprüft, bewertet und eingezeichnet.
Zusätzlich wird das ganzzahlige Optimum systematisch bestimmt statt behauptet.
"""

import itertools
import os

import numpy as np
import matplotlib
matplotlib.use("Agg")            # kein Bildschirm nötig
import matplotlib.pyplot as plt

OUTPUT_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), "output")
os.makedirs(OUTPUT_DIR, exist_ok=True)

# --- Modell (identisch zum Kapitel Einfuehrung) ----------------------------
#   max 150*xA + 250*xB
#   u.d.N.  2*xA + 5*xB <= 40   (vCPU)
#           4*xA + 6*xB <= 60   (RAM)
#           1*xA + 0*xB <=  8   (Marktliquidität)
c = np.array([150.0, 250.0])
A = np.array([[2.0, 5.0], [4.0, 6.0], [1.0, 0.0]])
b = np.array([40.0, 60.0, 8.0])
restriktionsnamen = ["vCPU", "RAM", "Marktlimit"]


def ist_zulaessig(punkt, tol=1e-7):
    return np.all(A @ punkt <= b + tol) and np.all(punkt >= -tol)


def berechne_ecken():
    """
    Ecken = Schnittpunkte je zweier Begrenzungsgeraden, die zulässig sind.
    Begrenzungen sind die 3 Restriktionen plus die beiden Achsen xA=0, xB=0.
    """
    geraden = [(A[i], b[i]) for i in range(len(b))]
    geraden.append((np.array([1.0, 0.0]), 0.0))     # xA = 0
    geraden.append((np.array([0.0, 1.0]), 0.0))     # xB = 0

    ecken = []
    for (n1, d1), (n2, d2) in itertools.combinations(geraden, 2):
        M = np.array([n1, n2])
        if abs(np.linalg.det(M)) < 1e-9:            # parallel -> kein Schnittpunkt
            continue
        p = np.linalg.solve(M, np.array([d1, d2]))
        if ist_zulaessig(p) and not any(np.allclose(p, e) for e in ecken):
            ecken.append(p)
    return np.array(ecken)


def bestes_ganzzahliges():
    """Vollständige Suche über das kleine Gitter - hier zulässig, weil winzig."""
    bester_wert, bester_punkt = -np.inf, None
    for xa in range(0, 21):
        for xb in range(0, 21):
            p = np.array([float(xa), float(xb)])
            if ist_zulaessig(p) and c @ p > bester_wert:
                bester_wert, bester_punkt = c @ p, p
    return bester_punkt, bester_wert


# --- Analyse ---------------------------------------------------------------
ecken = berechne_ecken()
werte = ecken @ c
reihenfolge = np.argsort(-werte)

print("=" * 62)
print("  ECKEN DES ZULÄSSIGEN POLYEDERS (nach Zielwert sortiert)")
print("=" * 62)
print(f"{'x_A':>8} {'x_B':>8} {'Z = 150 xA + 250 xB':>24}")
print("-" * 62)
for i in reihenfolge:
    print(f"{ecken[i, 0]:>8.2f} {ecken[i, 1]:>8.2f} {werte[i]:>24,.2f} EUR")

lp_punkt, lp_wert = ecken[reihenfolge[0]], werte[reihenfolge[0]]
ip_punkt, ip_wert = bestes_ganzzahliges()

print("-" * 62)
print(f"Kontinuierliches Optimum (LP): x = ({lp_punkt[0]:.2f}, {lp_punkt[1]:.2f}), "
      f"Z = {lp_wert:,.2f} EUR")
print(f"Ganzzahliges Optimum (IP):     x = ({ip_punkt[0]:.0f}, {ip_punkt[1]:.0f}), "
      f"Z = {ip_wert:,.2f} EUR")
print(f"Preis der Ganzzahligkeit:      {lp_wert - ip_wert:,.2f} EUR "
      f"({(1 - ip_wert / lp_wert) * 100:.2f} %)")
print("=" * 62)

# --- Zeichnung -------------------------------------------------------------
gitter = np.linspace(0, 15, 400)
xa_gitter, xb_gitter = np.meshgrid(gitter, gitter)

plt.figure(figsize=(10, 8))

# Restriktionsgeraden
plt.plot(gitter, (40 - 2 * gitter) / 5, color="tab:blue", lw=2,
         label=r"$2x_A + 5x_B \leq 40$ (vCPU)")
plt.plot(gitter, (60 - 4 * gitter) / 6, color="tab:green", lw=2,
         label=r"$4x_A + 6x_B \leq 60$ (RAM)")
plt.axvline(x=8, color="tab:orange", lw=2, label=r"$x_A \leq 8$ (Marktlimit)")

# Zulässiger Bereich
maske = ((2 * xa_gitter + 5 * xb_gitter <= 40) & (4 * xa_gitter + 6 * xb_gitter <= 60)
         & (xa_gitter <= 8) & (xa_gitter >= 0) & (xb_gitter >= 0))
plt.imshow(maske.astype(int), extent=(0, 15, 0, 15), origin="lower",
           cmap="Greys", alpha=0.25, aspect="auto")

# Höhenlinien der Zielfunktion
Z = 150 * xa_gitter + 250 * xb_gitter
hoehen = plt.contour(xa_gitter, xb_gitter, Z, levels=[500, 1000, 1500, 2000, 2375],
                     colors="purple", linestyles="--", alpha=0.7)
plt.clabel(hoehen, inline=True, fontsize=9, fmt="Z = %1.0f EUR")

# Ecken einzeichnen - jetzt werden sie tatsächlich benutzt
plt.scatter(ecken[:, 0], ecken[:, 1], s=70, facecolors="white",
            edgecolors="black", zorder=4, label="Ecken des Polyeders")
for e, w in zip(ecken, werte):
    plt.annotate(f"({e[0]:.1f}, {e[1]:.1f})\nZ={w:,.0f}", (e[0], e[1]),
                 textcoords="offset points", xytext=(6, 6), fontsize=8)

plt.scatter([lp_punkt[0]], [lp_punkt[1]], color="purple", marker="D", s=110, zorder=5,
            label=f"LP-Optimum ({lp_punkt[0]:.1f}, {lp_punkt[1]:.1f})")
plt.scatter([ip_punkt[0]], [ip_punkt[1]], color="red", s=170, zorder=6,
            label=f"Ganzzahliges Optimum ({ip_punkt[0]:.0f}, {ip_punkt[1]:.0f})")

plt.xlim(0, 12)
plt.ylim(0, 10)
plt.xlabel("Anzahl Arbitrage-Bots ($x_A$)", fontsize=11)
plt.ylabel("Anzahl Trendfolge-Bots ($x_B$)", fontsize=11)
plt.title("Polyeder des zulässigen Bereichs mit Niveaulinien der Zielfunktion", fontsize=13)
plt.grid(True, linestyle=":", alpha=0.6)
plt.legend(loc="upper right", framealpha=0.9)
plt.tight_layout()
ziel = os.path.join(OUTPUT_DIR, "feasible_region_2d.png")
plt.savefig(ziel, dpi=150)
print(f"Visualisierung gespeichert unter '{ziel}'")

Erwartete Ausgabe:

==============================================================
  ECKEN DES ZULÄSSIGEN POLYEDERS (nach Zielwert sortiert)
==============================================================
     x_A      x_B      Z = 150 xA + 250 xB
--------------------------------------------------------------
    7.50     5.00                 2,375.00 EUR
    8.00     4.67                 2,366.67 EUR
    0.00     8.00                 2,000.00 EUR
    8.00     0.00                 1,200.00 EUR
    0.00     0.00                     0.00 EUR
--------------------------------------------------------------
Kontinuierliches Optimum (LP): x = (7.50, 5.00), Z = 2,375.00 EUR
Ganzzahliges Optimum (IP):     x = (7, 5), Z = 2,300.00 EUR
Preis der Ganzzahligkeit:      75.00 EUR (3.16 %)
==============================================================

Zwei Beobachtungen, die Kapitel 6 vorbereiten:

1. Das LP-Optimum liegt bei (7{,}5;\ 5{,}0) — kein zulässiger Betriebszustand, denn ein halber Bot existiert nicht.
2. Das ganzzahlige Optimum (7, 5) ist nicht durch Runden entstanden: Aufrunden auf (8, 5) wäre unzulässig (2\cdot8+5\cdot5=41>40). Hier ging Abrunden gut — Kapitel 6 zeigt Fälle, in denen Runden dramatisch scheitert.

2.7 Wenn die Zahlen nicht zusammenpassen: Kondition und Skalierung

Bis hierher war die lineare Algebra exakt. Ein Rechner rechnet aber nicht exakt, sondern mit rund 16 signifikanten Stellen — und Ihre Betriebsdaten sind ohnehin nur auf wenige Stellen genau bekannt. Dieser Abschnitt beantwortet die Frage, die genau daraus entsteht:

Wenn meine Eingabedaten leicht ungenau sind — wie ungenau ist dann meine Lösung?

Die Antwort trägt einen Namen: die Konditionszahl \kappa(\mathbf{A}). Sie ist der wichtigste Begriff dieses Kapitels für den Produktivbetrieb, weil er erklärt, warum ein mathematisch korrektes Modell trotzdem unbrauchbare Ergebnisse liefern kann — oder sich mit einem sachlich falschen INFEASIBLE verabschiedet.

Die geometrische Vorstellung: zwei fast parallele Linien

Zwei Nebenbedingungen legen im Zweidimensionalen einen Schnittpunkt fest. Stehen die zugehörigen Geraden kreuzweise aufeinander, ist der Schnittpunkt robust: Verschiebt man eine Gerade um einen Millimeter, wandert der Schnittpunkt um einen Millimeter.

Verlaufen die Geraden dagegen fast parallel, entsteht ein spitzer, langgezogener Keil. Derselbe Millimeter Verschiebung schickt den Schnittpunkt jetzt meterweit den Keil entlang. Genau das misst \kappa(\mathbf{A}): den Verstärkungsfaktor zwischen Datenfehler und Lösungsfehler.

\frac{\lVert \Delta \mathbf{x}\rVert}{\lVert \mathbf{x}\rVert} \ \le\ \kappa(\mathbf{A}) \cdot \frac{\lVert \Delta \mathbf{b}\rVert}{\lVert \mathbf{b}\rVert}, \qquad \kappa(\mathbf{A}) = \lVert \mathbf{A}\rVert \cdot \lVert \mathbf{A}^{-1}\rVert = \frac{\sigma_{\max}}{\sigma_{\min}}

🎯 Merksatz \kappa < 10^3 ist unbedenklich, \kappa > 10^6 verdient eine Prüfung, \kappa > 10^{10} heißt: Das Ergebnis ist Rauschen mit Nachkommastellen. Rechnen Sie \kappa mit np.linalg.cond(A) aus — es kostet eine Zeile und beantwortet die Frage, die sonst monatelang niemand stellt.

Woher schlechte Kondition im Alltag kommt

Fast nie aus exotischer Mathematik — fast immer aus Einheiten:

Die dritte Zeile ist der häufigste selbstgemachte Fall: Ein unnötig großes Big-M bläht \kappa auf und lässt Solver bei völlig harmlosen Modellen stundenlang suchen oder falsche Ergebnisse melden.

Ruiz-Equilibrierung: das Modell gesundrechnen

Die Gegenmaßnahme ist erstaunlich schlicht. Man multipliziert jede Zeile und jede Spalte der Matrix mit einem Faktor, sodass die Beträge überall in derselben Größenordnung landen:

\tilde{\mathbf{A}} = \mathbf{D}_r\,\mathbf{A}\,\mathbf{D}_c

mit Diagonalmatrizen \mathbf{D}_r (Zeilen) und \mathbf{D}_c (Spalten). Das Ruiz-Verfahren bestimmt diese Faktoren iterativ: In jedem Durchlauf wird jede Zeile durch die Wurzel ihres größten Betrags geteilt, danach jede Spalte. Nach wenigen Durchläufen liegen alle Zeilen- und Spaltenmaxima bei 1.

Wirtschaftlich passiert dabei nichts: Eine Zeilenskalierung heißt „diese Nebenbedingung in einer anderen Einheit messen“, eine Spaltenskalierung „diese Variable in einer anderen Einheit messen“. Die Lösung rechnet man mit \mathbf{x} = \mathbf{D}_c\,\tilde{\mathbf{x}} zurück.

🎯 Merksatz Jeder ernsthafte Solver skaliert intern selbst — HiGHS, Gurobi, CP-SAT alle. Verlassen Sie sich trotzdem nicht darauf: Die interne Skalierung repariert die Rechnung, nicht die Modellierung. Wenn Ihre Koeffizienten zwölf Größenordnungen überspannen, sagt Ihnen das etwas über Ihr Modell, nicht über den Solver.

Das Experiment

Das folgende Programm führt beides vor — die Fehlerverstärkung und ihre Behebung — und schließt mit der praktischen Konsequenz, die Sie in Kapitel 6 wieder brauchen werden: warum eine Binärvariable mit dem Wert 0.99999998 niemals mit int() gerundet werden darf.

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

Erwartete Ausgabe:

==========================================================================
  1. KONDITIONSZAHL ALS FEHLERVERSTAERKER
==========================================================================
Zwei Gleichungssysteme, beide exakt loesbar mit x = (1, 1).
Die rechte Seite wird um relativ 1e-10 gestoert - so viel Unsicherheit
steckt in JEDER gemessenen Betriebszahl allemal.

Matrix                       kappa(A)  Verst. median   Verst. max
--------------------------------------------------------------------------
gut konditioniert            2.62e+00       1.86e+00     2.56e+00
fast parallele Zeilen        4.00e+08       2.57e+08     4.00e+08

Lesart: Bei der zweiten Matrix wird aus einem Datenfehler in der
10. Nachkommastelle ein Loesungsfehler in der 2. Nachkommastelle.
Das Modell ist mathematisch korrekt - und praktisch wertlos.

==========================================================================
  2. RUIZ-EQUILIBRIERUNG: EINHEITEN GERADERUECKEN
==========================================================================
Ein Produktionsmodell, in dem vier Ressourcen in voellig
verschiedenen Einheiten gemessen werden:
  Zeile 1: Kapitalbindung in Euro          (Groessenordnung 1e7)
  Zeile 2: Katalysatorverbrauch in Tonnen  (Groessenordnung 1e-3)
  Zeile 3: Energie in Wattsekunden         (Groessenordnung 1e5)
  Zeile 4: Ausschussquote als Anteil       (Groessenordnung 1e-2)

                                 kappa(A)  groesster Eintrag  kleinster Eintrag
--------------------------------------------------------------------------
vor der Skalierung               8.67e+10           1.52e+07           1.35e-03
nach Ruiz-Equilibrierung         4.69e+01           1.00e+00           4.24e-01

Ruecktransformation x = c * x_s: groesste Abweichung zur wahren Loesung 2.35e-14
Die Skalierung ist also verlustfrei - sie aendert nur die Zahlen,
nicht das Problem.

==========================================================================
  3. TOLERANZEN: WARUM int() DIE FALSCHE RUNDUNG IST
==========================================================================
Solver liefert x = 1.0

    Solverwert    int()   round() Kommentar
--------------------------------------------------------------------------
    0.99999998        0         1 int() macht aus einer JA- eine NEIN-Entscheidung
    1.00000002        1         1 hier ginge int() zufaellig gut - Verlass ist keiner
    0.49999999        0         0 echt unentschieden: Modell oder Toleranz pruefen!
    2.99999990        2         3 3 Maschinen werden zu 2 - der Plan geht nicht auf

Richtige Vorgehensweise: gegen die Solver-Toleranz pruefen,
dann erst runden - und den Zweifelsfall melden statt still zu raten.
  sichere_ganzzahl(0.99999998) = 1
  sichere_ganzzahl(1.00000002) = 1
  sichere_ganzzahl(0.49999999) -> ValueError: 0.49999999 ist 5.00e-01 von der naechsten ganzen Zahl entfernt - das ist mehr als die Toleranz 1e-06. Ganzzahligkeit im Modell pruefen.
  sichere_ganzzahl(2.9999999) = 3

==========================================================================
Merksatz: Skalieren Sie Ihre Daten, BEVOR der Solver sie sieht -
und runden Sie Solver-Ergebnisse NIE ohne Toleranzpruefung.
==========================================================================

💻 Code-Durchgang

Stelle	Was passiert	Warum das die Kernaussage trägt
Störung in zufälliger Richtung	200 Versuche statt einer festen Störung	Eine proportionale Störung (b \cdot 1{,}0000000001) trifft die empfindliche Richtung meist gar nicht — man sähe fälschlich gar keinen Effekt. Die Schwäche zeigt sich nur, wenn man in alle Richtungen stößt.
Verst. max = 4.00e+08	erreicht exakt \kappa(\mathbf{A})	Die Schranke aus der Formel oben ist scharf: Es gibt eine Richtung, in der die volle Verstärkung eintritt. Sie ist keine pessimistische Abschätzung.
Ruiz: quadratische Matrix	4 \times 4 statt 4 \times 6	Nur bei eindeutiger Lösung beweist die Rückrechnung etwas. Bei unterbestimmten Systemen liefert lstsq die Minimum-Norm-Lösung, und die ist unter Skalierung nicht invariant — eine Falle, in die man leicht tappt.
Abweichung 2.35e-14	Rückrechnung trifft die wahre Lösung	Der Beleg, dass Skalieren nichts kaputt macht: Es ändert die Zahlen, nicht das Problem.
sichere_ganzzahl()	prüft Abstand zur nächsten ganzen Zahl	Der Zweifelsfall (0.49999999) wird gemeldet, nicht stillschweigend geraten. Genau diese Funktion gehört in jedes MILP-Auswertungsskript.

⚠️ Typische Fehler

• int() auf Solver-Ergebnisse. int(0.99999998) == 0 — aus einem „ja, Standort eröffnen“ wird ein „nein“. int() schneidet ab, es rundet nicht. Und selbst round() ist nur mit vorheriger Toleranzprüfung vertretbar.
• Große Zahlen für „Sicherheit“. Ein Big-M von 10^9 „damit es garantiert reicht“ verschlechtert \kappa um neun Größenordnungen und erzeugt genau die Instabilität, die man vermeiden wollte.
• Auf die interne Solver-Skalierung vertrauen. Sie hilft der Rechnung, aber sie kann Ihnen nicht sagen, dass Sie Stück und Tonnen vermischt haben.
• \kappa nie nachschauen. Eine Zeile np.linalg.cond(A) beim Modellaufbau, geloggt neben der Laufzeit, hätte in vielen Projekten Wochen Fehlersuche gespart.

2.8 Übungsaufgaben

Lösungen: Abschnitt A.2.

Aufgabe 2.1 ⭐ — Matrixform lesen. Gegeben \mathbf{c} = (4, 1, 6)^\top, \mathbf{A} = \begin{pmatrix} 1 & 2 & 0 \\ 0 & 1 & 3 \end{pmatrix}, \mathbf{b} = (10, 12)^\top. Schreiben Sie das Modell ausgeschrieben hin. Wie viele Variablen und wie viele Nebenbedingungen hat es?

Aufgabe 2.2 ⭐ — Standardform herstellen. Bringen Sie in die Form \min \mathbf{c}^\top\mathbf{x} u. d. N. \mathbf{A}\mathbf{x} \le \mathbf{b}: \max\ 7x_1 - 2x_2 \quad \text{u. d. N.}\quad 4x_1 + x_2 \ge 20,\quad x_1 - x_2 = 3,\quad x_1, x_2 \ge 0

Aufgabe 2.3 ⭐⭐ — Ecken von Hand. Für \max 2x_1 + 3x_2 u. d. N. x_1 + x_2 \le 8, 2x_1 + x_2 \le 12, x_1, x_2 \ge 0: (a) Zeichnen Sie den zulässigen Bereich auf kariertes Papier. (b) Bestimmen Sie alle Ecken. (c) Werten Sie Z in jeder Ecke aus und geben Sie das Optimum an. (d) Was passiert, wenn die Zielfunktion zu \max 2x_1 + 2x_2 wird? Wie viele optimale Lösungen gibt es dann?

Aufgabe 2.4 ⭐⭐ — Konvexität beurteilen. Konvex oder nicht? Begründen Sie jeweils kurz: (a) f(x) = 5x + 3 (b) f(x) = x^4 (c) f(x) = \sqrt{x} für x \ge 0 (d) f(x) = \ln(x) für x > 0 (e) f(x_1,x_2) = x_1^2 + x_2^2 (f) Die Menge \{(x_1,x_2) : x_1^2 + x_2^2 \le 4\} (g) Die Menge \{(x_1,x_2) : x_1 \cdot x_2 \ge 1,\ x_1,x_2 > 0\}

Aufgabe 2.5 ⭐⭐ — Positive Semidefinitheit prüfen. Prüfen Sie mit NumPy, welche dieser Matrizen als Kovarianzmatrix taugen (alle Eigenwerte \ge 0), und geben Sie bei Nichteignung die implizierte Korrelation an: \mathbf{P}_1=\begin{pmatrix}4&1\\1&9\end{pmatrix},\quad \mathbf{P}_2=\begin{pmatrix}4&7\\7&9\end{pmatrix},\quad \mathbf{P}_3=\begin{pmatrix}1&0.5&0.5\\0.5&1&0.5\\0.5&0.5&1\end{pmatrix}

Aufgabe 2.6 ⭐⭐⭐ — Eigene Visualisierung. Passen Sie Visualisierung_Loesungsraum.py an das Bäckerei-Problem aus Kapitel 1 (die Bäckerei-Handrechnung aus Kapitel 1) an: \max 2{,}5x_1 + 3{,}0x_2 u. d. N. 0{,}5x_1+0{,}6x_2 \le 90, 4x_1+3x_2 \le 600, x_1 \ge 40. Geben Sie alle Ecken mit Zielwerten aus und zeichnen Sie das Polyeder.

Aufgabe 2.7 ⭐⭐⭐ — Konvexität einer Kovarianzmatrix reparieren. Die Matrix \mathbf{P}_2 aus der Aufgabe Positive Semidefinitheit prüfen ist nicht positiv semidefinit. Schreiben Sie eine Funktion repariere(P), die die negativen Eigenwerte auf 0 setzt und die Matrix rekonstruiert (Eigenwert-Clipping): \mathbf{P}_{\text{rep}} = \mathbf{V}\max(\boldsymbol{\Lambda},0)\mathbf{V}^\top. Prüfen Sie das Ergebnis und vergleichen Sie es mit dem Original.

2.9 Finde den Denkfehler

2.10 Micro-Quiz

2.11 Selbsttest

Antworten: Anhang A.

1. Was bedeutet der Ausdruck \mathbf{c}^\top\mathbf{x} inhaltlich, und was für ein Objekt ist das Ergebnis?
2. Warum steht in der Technologiematrix eine 0, wenn eine Variable in einer Bedingung nicht vorkommt — warum lässt man den Eintrag nicht weg?
3. Formulieren Sie den Fundamentalsatz der linearen Optimierung und nennen Sie seinen praktischen Nutzen.
4. Warum ist ein ganzzahliges Problem nicht konvex, obwohl alle Funktionen darin linear sind?
5. Ein Kollege sagt: „Mein Optimierer liefert je nach Startwert unterschiedliche Ergebnisse — das ist wohl ein Bug.“ Was antworten Sie?

2.12 Zusammenfassung

• Matrixform bündelt beliebig viele Variablen und Bedingungen in drei Zeilen: \min \mathbf{c}^\top\mathbf{x} u. d. N. \mathbf{A}\mathbf{x} \le \mathbf{b}, \mathbf{x} \ge \mathbf{0}. Jede Zeile von \mathbf{A} ist eine Ressource, jede Spalte eine Variable.
• Der zulässige Bereich ist der Schnitt aller Halbräume: ein konvexes Polyeder.
• Der Fundamentalsatz garantiert, dass ein Optimum in einer Ecke liegt — endlich viele Kandidaten statt unendlich vieler Punkte.
• Kondition und Skalierung entscheiden darüber, ob aus mathematisch korrekten Zahlen auch brauchbare werden: \kappa(\mathbf{A}) ist der Verstärkungsfaktor zwischen Daten- und Lösungsfehler, gemischte Einheiten sind seine häufigste Ursache, Ruiz-Equilibrierung die Gegenmaßnahme — und Solver-Ergebnisse werden nur nach Toleranzprüfung gerundet.
• Konvexität ist die eigentliche Trennlinie zwischen leicht und schwer. Bei konvexen Problemen ist jedes lokale Optimum global; bei nicht-konvexen hängt das Ergebnis vom Startpunkt ab.
• Prüfen Sie Matrizen, bevor Sie sie als Kovarianzmatrix verwenden: alle Eigenwerte müssen \ge 0 sein.

Ausblick. Kapitel 3 ordnet das Python-Ökosystem: Welche Bibliothek löst welche Problemklasse, und warum gibt es überhaupt mehrere?