operations_research/bilder_04/erzeuge_kkt.py

#!/usr/bin/env python3

# erzeuge_kkt.py
"""
Erzeugt das KKT-Diagramm zum QP/NLP-Kapitel:

    bilder_04/kap_qp_kkt.svg

Zwei Bilder nebeneinander:

  Links  - die Handrechnung des Kapitels, eindimensional:
               min (x-5)^2   u.d.N.  x - 3 <= 0   ->  x* = 3, lambda* = 4
           Sichtbar wird, was lambda BEDEUTET: Es ist die Steigung von f am
           Rand, also der Preis dafuer, dass die Grenze nicht weiter darf.
  Rechts - dieselbe Aussage in zwei Dimensionen: Am Optimum zeigen der
           Zielgradient und der mit lambda gewichtete Restriktionsgradient
           GENAU GEGENEINANDER und heben sich auf.

Warum neu: Im Vorgaengerbild war der Zielpfeil mit "-grad f" beschriftet und
zeigte entgegen lambda*grad g. Das ist die Konstellation einer MAXIMIERUNG. Das
Kapitel formuliert KKT aber fuer die Minimierung - dort gilt
grad f + lambda*grad g = 0, also zeigen grad f (ohne Minus) und lambda*grad g
gegeneinander. Eine Beschriftung, ein Vorzeichen, und die Aussage kippt.

Alle Werte werden gerechnet, nicht eingetragen: Das Optimum kommt aus
scipy.optimize, die Gradienten aus der Ableitung, und das Skript prueft alle
vier KKT-Bedingungen numerisch, bevor es zeichnet.

Aufruf (aus dem Repository-Wurzelverzeichnis):
    python3 bilder_04/erzeuge_kkt.py

Benoetigt: numpy, scipy, matplotlib
"""

from __future__ import annotations

import os

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import minimize

# Reproduzierbare SVG-Ausgabe (siehe erzeuge_titelseite.py)
plt.rcParams["svg.hashsalt"] = "or-mit-python-v04"

BASIS = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
BILDER = os.path.join(BASIS, "bilder_04")

INDIGO, CYAN, AMBER, GRUEN, ROT = (
    "#4338ca", "#0891b2", "#b45309", "#059669", "#be123c")

# --- Eindimensional: die Handrechnung des Kapitels --------------------------
ZIEL_1D = 5.0          # unbeschraenktes Minimum von (x-5)^2
GRENZE_1D = 3.0        # x <= 3


# --- Zweidimensional: dieselbe Struktur, ein Ziel mehr ----------------------
ZENTRUM = np.array([3.5, 3.5])     # unbeschraenktes Minimum
BUDGET = 4.0                       # x1 + x2 <= 4


def f(x):
    return float(np.sum((np.asarray(x) - ZENTRUM) ** 2))


def grad_f(x):
    return 2.0 * (np.asarray(x) - ZENTRUM)


def g(x):
    return float(np.sum(x) - BUDGET)


GRAD_G = np.array([1.0, 1.0])


def loese_2d():
    ergebnis = minimize(f, x0=np.zeros(2), method="SLSQP",
                        constraints=[{"type": "ineq", "fun": lambda x: -g(x)}])
    if not ergebnis.success:
        raise SystemExit(f"SLSQP gescheitert: {ergebnis.message}")
    x = np.asarray(ergebnis.x)
    # lambda aus der Stationaritaet: grad f + lambda * grad g = 0
    lam = float(-(grad_f(x) @ GRAD_G) / (GRAD_G @ GRAD_G))
    return x, lam


def pruefe_kkt(x, lam, toleranz=1e-5) -> None:
    """Alle vier Bedingungen, bevor gezeichnet wird - nicht danach."""
    stationaritaet = grad_f(x) + lam * GRAD_G
    pruefungen = {
        "Stationarität grad f + lambda*grad g = 0":
            float(np.linalg.norm(stationaritaet)) < toleranz,
        "primale Zulässigkeit g(x) <= 0": g(x) < toleranz,
        "duale Zulässigkeit lambda >= 0": lam > -toleranz,
        "komplementärer Schlupf lambda*g(x) = 0": abs(lam * g(x)) < toleranz,
    }
    for name, erfuellt in pruefungen.items():
        print(f"  {'OK ' if erfuellt else 'NEIN'}  {name}")
    if not all(pruefungen.values()):
        raise SystemExit("FEHLER: KKT nicht erfuellt - es wird nichts gezeichnet.")


def zahl(wert: float) -> str:
    text = f"{wert:.2f}".rstrip("0").rstrip(".")
    return text.replace(".", ",")


def zeichne_1d(achse, lam_1d: float) -> None:
    x = np.linspace(0.5, 6.5, 400)
    achse.plot(x, (x - ZIEL_1D) ** 2, color=INDIGO, linewidth=2.0,
               label="$f(x) = (x-5)^2$")
    achse.axvspan(GRENZE_1D, 6.5, color=ROT, alpha=0.07)
    achse.axvline(GRENZE_1D, color=ROT, linewidth=1.8, linestyle="--")
    achse.annotate("unzulässig:\n$x > 3$", (GRENZE_1D, 20),
                   textcoords="offset points", xytext=(10, 0),
                   fontsize=8.5, color=ROT, va="center")

    achse.scatter([ZIEL_1D], [0.0], s=90, facecolors="none",
                  edgecolors="#94a3b8", linewidths=2.0, zorder=4)
    achse.annotate("freies Minimum $x=5$\n— nicht erlaubt", (ZIEL_1D, 0.0),
                   textcoords="offset points", xytext=(10, 16),
                   ha="left", fontsize=8.5, color="#64748b")

    wert = (GRENZE_1D - ZIEL_1D) ** 2
    achse.scatter([GRENZE_1D], [wert], s=130, color=GRUEN, zorder=5)
    achse.annotate(f"$x^*=3$, $f=4$", (GRENZE_1D, wert),
                   textcoords="offset points", xytext=(-10, 12),
                   ha="right", fontsize=9, color=GRUEN, fontweight="bold")

    # Die Tangente am Rand: ihre Steigung ist -lambda*
    steigung = 2 * (GRENZE_1D - ZIEL_1D)
    nahe = np.linspace(GRENZE_1D - 1.2, GRENZE_1D + 0.8, 20)
    achse.plot(nahe, wert + steigung * (nahe - GRENZE_1D), color=AMBER,
               linewidth=1.6, linestyle="-",
               label=f"Steigung $-\\lambda^* = {zahl(steigung)}$")

    achse.set_xlabel("$x$")
    achse.set_ylabel("$f(x)$")
    achse.set_title("Die Handrechnung: $\\lambda^*$ ist die Steigung am Rand",
                    fontsize=10.5)
    achse.set_ylim(-1.5, 22)
    achse.grid(linestyle=":", alpha=0.45)
    achse.legend(fontsize=8, loc="upper left", frameon=False)
    for seite in ("top", "right"):
        achse.spines[seite].set_visible(False)


def zeichne_2d(achse, x_stern, lam) -> None:
    gitter = np.linspace(-0.2, 5.4, 300)
    X1, X2 = np.meshgrid(gitter, gitter)
    werte = (X1 - ZENTRUM[0]) ** 2 + (X2 - ZENTRUM[1]) ** 2
    achse.contour(X1, X2, werte, levels=[0.5, 2, 4.5, 8, 12.5, 18],
                  colors=[INDIGO], linewidths=0.8, alpha=0.55)

    rand = np.linspace(-0.2, 5.4, 10)
    achse.plot(rand, BUDGET - rand, color=ROT, linewidth=2.0,
               label="$g(x)=x_1+x_2-4=0$")
    achse.fill_between(rand, -0.2, np.minimum(BUDGET - rand, 5.4),
                       color=CYAN, alpha=0.08)

    achse.scatter(*ZENTRUM, s=90, facecolors="none", edgecolors="#94a3b8",
                  linewidths=2.0, zorder=4)
    achse.annotate("freies Minimum\n— unzulässig", ZENTRUM,
                   textcoords="offset points", xytext=(12, -8),
                   va="top", fontsize=8.5, color="#64748b")

    achse.scatter(*x_stern, s=140, color=GRUEN, zorder=6)
    achse.annotate(f"$x^* = ({zahl(x_stern[0])}; {zahl(x_stern[1])})$",
                   x_stern, textcoords="offset points", xytext=(-12, -20),
                   ha="right", fontsize=9, color=GRUEN, fontweight="bold")

    zielgradient = grad_f(x_stern)
    achse.annotate("", xy=x_stern + zielgradient * 0.45, xytext=x_stern,
                   arrowprops=dict(arrowstyle="-|>", color=INDIGO, lw=2.4))
    achse.annotate("$\\nabla f$ — zeigt bergauf,\nweg vom freien Minimum",
                   x_stern + zielgradient * 0.45,
                   textcoords="offset points", xytext=(-8, -8),
                   ha="right", va="top", fontsize=8.5, color=INDIGO)

    restriktionsgradient = lam * GRAD_G
    achse.annotate("", xy=x_stern + restriktionsgradient * 0.45, xytext=x_stern,
                   arrowprops=dict(arrowstyle="-|>", color=AMBER, lw=2.4))
    achse.annotate(f"$\\lambda^*\\nabla g$ mit $\\lambda^* = {zahl(lam)}$",
                   x_stern + restriktionsgradient * 0.45,
                   textcoords="offset points", xytext=(-10, 8),
                   ha="right", fontsize=8.5, color=AMBER)

    achse.set_xlim(-0.2, 5.4)
    achse.set_ylim(-0.2, 5.4)
    achse.set_xlabel("$x_1$")
    achse.set_ylabel("$x_2$")
    achse.set_title("Am Optimum heben sich beide Gradienten auf:\n"
                    "$\\nabla f + \\lambda^*\\nabla g = 0$", fontsize=10.5)
    achse.grid(linestyle=":", alpha=0.45)
    achse.legend(fontsize=8, loc="upper right", frameon=False)
    for seite in ("top", "right"):
        achse.spines[seite].set_visible(False)


if __name__ == "__main__":
    x_stern, lam = loese_2d()
    print(f"2D: x* = ({x_stern[0]:.4f}; {x_stern[1]:.4f}), lambda* = {lam:.4f}")
    pruefe_kkt(x_stern, lam)

    lam_1d = -2 * (GRENZE_1D - ZIEL_1D)
    print(f"1D: x* = {GRENZE_1D:g}, lambda* = {lam_1d:g} "
          f"(Handrechnung: 3 und 4)")
    if abs(lam_1d - 4.0) > 1e-9:
        raise SystemExit("FEHLER: 1D-Lambda weicht von der Handrechnung ab.")

    figur, (links, rechts) = plt.subplots(1, 2, figsize=(11.0, 4.8))
    zeichne_1d(links, lam_1d)
    zeichne_2d(rechts, x_stern, lam)
    figur.tight_layout()
    os.makedirs(BILDER, exist_ok=True)
    for endung in ("svg",):
        pfad = os.path.join(BILDER, f"kap_qp_kkt.{endung}")
        figur.savefig(pfad, format=endung, dpi=160,
                      metadata={"Date": None} if endung == "svg" else None)
        print(f"geschrieben: {pfad}")
    plt.close(figur)