231 lines
8.8 KiB
Python
231 lines
8.8 KiB
Python
|
#!/usr/bin/env python3
|
|||
|
|
|
|||
|
|
# erzeuge_runden.py
|
|||
|
|
"""
|
|||
|
|
Erzeugt das Diagramm "Warum Runden scheitert" zum MILP-Kapitel:
|
|||
|
|
|
|||
|
|
bilder_04/kap_milp_runden.svg statisch, fuer das PDF
|
|||
|
|
bilder_04/kap_milp_runden.png Rasterfassung
|
|||
|
|
|
|||
|
|
Zwei Bilder nebeneinander, beide aus dem Kapitel selbst:
|
|||
|
|
|
|||
|
|
Links - die Handrechnung des Kapitels, gezeichnet:
|
|||
|
|
max x1 + x2 u.d.N. 2 x1 + 2 x2 <= 3, x ganzzahlig
|
|||
|
|
Der Witz ist die Entartung: Die Zielfunktion laeuft PARALLEL zur
|
|||
|
|
Restriktion, also ist die ganze Kante optimal, nicht eine Ecke.
|
|||
|
|
Aufrunden auf (1,1) fuehrt aus dem zulaessigen Bereich heraus,
|
|||
|
|
Abrunden auf (0,0) auf den Wert null. Beides sieht man sofort.
|
|||
|
|
Rechts - die Messung aus Runden_Gegenbeispiel.py, mit derselben Saat
|
|||
|
|
nachgerechnet: Wie oft ist die aufgerundete Loesung unzulaessig,
|
|||
|
|
und was kostet Abrunden - mit und ohne gieriges Auffuellen?
|
|||
|
|
|
|||
|
|
Warum ueberhaupt neu: Das Vorgaengerbild zeigte ein LP, das im Buch nicht
|
|||
|
|
vorkommt (LP-Optimum bei (1,8; 3,4), Z = 42,6), ohne die Zielfunktion zu nennen
|
|||
|
|
- unmittelbar ueber der Handrechnung, die etwas voellig anderes rechnet.
|
|||
|
|
|
|||
|
|
Aufruf (aus dem Repository-Wurzelverzeichnis, dauert rund eine Minute):
|
|||
|
|
python3 bilder_04/erzeuge_runden.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 linprog
|
|||
|
|
|
|||
|
|
# 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")
|
|||
|
|
|
|||
|
|
GROESSEN = [(5, 2), (10, 3), (20, 5), (40, 8)]
|
|||
|
|
INSTANZEN = 200
|
|||
|
|
SAAT = 2026
|
|||
|
|
|
|||
|
|
|
|||
|
|
# --- Messung: wortgleich zu Runden_Gegenbeispiel.py -------------------------
|
|||
|
|
|
|||
|
|
def erzeuge_instanz(n, m, rng):
|
|||
|
|
c = rng.integers(3, 20, size=n).astype(float)
|
|||
|
|
A = rng.integers(1, 9, size=(m, n)).astype(float)
|
|||
|
|
b = (A.sum(axis=1) * rng.uniform(0.25, 0.45)).round()
|
|||
|
|
return c, A, b
|
|||
|
|
|
|||
|
|
|
|||
|
|
def loese_lp(c, A, b, ganzzahlig=False):
|
|||
|
|
n = len(c)
|
|||
|
|
ergebnis = linprog(
|
|||
|
|
c=-c, A_ub=A, b_ub=b, bounds=[(0, None)] * n,
|
|||
|
|
integrality=np.ones(n) if ganzzahlig else None,
|
|||
|
|
method="highs")
|
|||
|
|
return (-ergebnis.fun, ergebnis.x) if ergebnis.success else (None, None)
|
|||
|
|
|
|||
|
|
|
|||
|
|
def abrunden_und_reparieren(x_lp, c, A, b):
|
|||
|
|
x = np.floor(x_lp + 1e-9)
|
|||
|
|
verbessert = True
|
|||
|
|
while verbessert:
|
|||
|
|
verbessert = False
|
|||
|
|
for j in np.argsort(-c):
|
|||
|
|
kandidat = x.copy()
|
|||
|
|
kandidat[j] += 1
|
|||
|
|
if np.all(A @ kandidat <= b + 1e-9):
|
|||
|
|
x = kandidat
|
|||
|
|
verbessert = True
|
|||
|
|
break
|
|||
|
|
return c @ x, x
|
|||
|
|
|
|||
|
|
|
|||
|
|
def messreihe():
|
|||
|
|
"""Dieselbe Schleife wie im Buchprogramm, mit derselben Saat.
|
|||
|
|
|
|||
|
|
Die Saat sitzt VOR der Schleife ueber die Groessen, nicht darin - so wie
|
|||
|
|
dort. Nur deshalb kommen hier dieselben Zahlen heraus wie in der
|
|||
|
|
abgedruckten Tabelle.
|
|||
|
|
"""
|
|||
|
|
rng = np.random.default_rng(SAAT)
|
|||
|
|
zeilen = []
|
|||
|
|
for n, m in GROESSEN:
|
|||
|
|
unzulaessig = 0
|
|||
|
|
verluste_ab, verluste_gierig = [], []
|
|||
|
|
gezaehlt = 0
|
|||
|
|
for _ in range(INSTANZEN):
|
|||
|
|
c, A, b = erzeuge_instanz(n, m, rng)
|
|||
|
|
z_lp, x_lp = loese_lp(c, A, b, ganzzahlig=False)
|
|||
|
|
z_ip, _ = loese_lp(c, A, b, ganzzahlig=True)
|
|||
|
|
if z_lp is None or z_ip is None or z_ip <= 0:
|
|||
|
|
continue
|
|||
|
|
gezaehlt += 1
|
|||
|
|
x_auf = np.ceil(x_lp - 1e-9)
|
|||
|
|
if np.any(A @ x_auf > b + 1e-9):
|
|||
|
|
unzulaessig += 1
|
|||
|
|
x_ab = np.floor(x_lp + 1e-9)
|
|||
|
|
verluste_ab.append(1.0 - (c @ x_ab) / z_ip)
|
|||
|
|
z_gierig, _ = abrunden_und_reparieren(x_lp, c, A, b)
|
|||
|
|
verluste_gierig.append(1.0 - z_gierig / z_ip)
|
|||
|
|
zeilen.append({
|
|||
|
|
"name": f"{n}×{m}",
|
|||
|
|
"unzulaessig": unzulaessig / gezaehlt,
|
|||
|
|
"ab": float(np.mean(verluste_ab)),
|
|||
|
|
"gierig": float(np.mean(verluste_gierig)),
|
|||
|
|
})
|
|||
|
|
return zeilen
|
|||
|
|
|
|||
|
|
|
|||
|
|
# --- Zeichnen ---------------------------------------------------------------
|
|||
|
|
|
|||
|
|
def zeichne_handrechnung(achse) -> None:
|
|||
|
|
"""Die Handrechnung des Kapitels: max x1+x2 u.d.N. 2x1+2x2 <= 3."""
|
|||
|
|
ecken = np.array([[0.0, 0.0], [1.5, 0.0], [0.0, 1.5]])
|
|||
|
|
achse.fill(ecken[:, 0], ecken[:, 1], color=CYAN, alpha=0.13, zorder=1)
|
|||
|
|
achse.plot([1.5, 0.0], [0.0, 1.5], color=CYAN, linewidth=2.0, zorder=3,
|
|||
|
|
label="$2x_1+2x_2 = 3$ — zugleich das LP-Optimum")
|
|||
|
|
|
|||
|
|
# Ganzzahlige Punkte im Fenster, zulaessig und unzulaessig getrennt
|
|||
|
|
for x1 in range(0, 3):
|
|||
|
|
for x2 in range(0, 3):
|
|||
|
|
zulaessig = 2 * x1 + 2 * x2 <= 3
|
|||
|
|
achse.scatter(x1, x2, s=48, zorder=4,
|
|||
|
|
color=INDIGO if zulaessig else "#cbd5e1")
|
|||
|
|
|
|||
|
|
achse.scatter([0.75], [0.75], marker="*", s=260, color=ROT, zorder=6)
|
|||
|
|
achse.annotate("LP-Lösung $(0{,}75;\\ 0{,}75)$\n$Z_{LP} = 1{,}5$",
|
|||
|
|
(0.75, 0.75), textcoords="offset points", xytext=(-12, -6),
|
|||
|
|
ha="right", va="top", fontsize=9, color=ROT,
|
|||
|
|
fontweight="bold")
|
|||
|
|
|
|||
|
|
achse.scatter([1], [1], s=150, facecolors="none", edgecolors=AMBER,
|
|||
|
|
linewidths=2.4, zorder=6)
|
|||
|
|
achse.annotate("aufgerundet $(1;1)$\n$2+2=4>3$ — unzulässig",
|
|||
|
|
(1, 1), textcoords="offset points", xytext=(12, -6),
|
|||
|
|
fontsize=8.5, color=AMBER, va="top")
|
|||
|
|
|
|||
|
|
achse.scatter([0], [0], s=150, facecolors="none", edgecolors="#64748b",
|
|||
|
|
linewidths=2.2, zorder=6)
|
|||
|
|
achse.annotate("abgerundet $(0;0)$\n$Z = 0$ — alles verschenkt",
|
|||
|
|
(0, 0), textcoords="offset points", xytext=(14, 8),
|
|||
|
|
fontsize=8.5, color="#475569")
|
|||
|
|
|
|||
|
|
for punkt in ((1, 0), (0, 1)):
|
|||
|
|
achse.scatter(*punkt, s=140, color=GRUEN, zorder=7)
|
|||
|
|
achse.annotate("wahres Optimum $Z_{IP}=1$", (1, 0),
|
|||
|
|
textcoords="offset points", xytext=(12, -4),
|
|||
|
|
va="top", fontsize=8.5, color=GRUEN, fontweight="bold")
|
|||
|
|
|
|||
|
|
achse.set_xlim(-0.35, 2.35)
|
|||
|
|
achse.set_ylim(-0.35, 2.35)
|
|||
|
|
achse.set_xlabel("$x_1$")
|
|||
|
|
achse.set_ylabel("$x_2$")
|
|||
|
|
achse.set_title("Die Handrechnung als Bild:\nrunden geht in beide Richtungen schief",
|
|||
|
|
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)
|
|||
|
|
|
|||
|
|
|
|||
|
|
def zeichne_messung(achse, zeilen) -> None:
|
|||
|
|
namen = [z["name"] for z in zeilen]
|
|||
|
|
stellen = np.arange(len(zeilen))
|
|||
|
|
breite = 0.26
|
|||
|
|
|
|||
|
|
achse.bar(stellen - breite, [z["unzulaessig"] * 100 for z in zeilen],
|
|||
|
|
breite, color=AMBER, label="Aufrunden unzulässig")
|
|||
|
|
achse.bar(stellen, [z["ab"] * 100 for z in zeilen],
|
|||
|
|
breite, color=ROT, label="Abrunden: mittlerer Verlust")
|
|||
|
|
achse.bar(stellen + breite, [z["gierig"] * 100 for z in zeilen],
|
|||
|
|
breite, color=GRUEN, label="abrunden + gierig auffüllen")
|
|||
|
|
|
|||
|
|
for stelle, zeile in zip(stellen, zeilen):
|
|||
|
|
for versatz, wert in ((-breite, zeile["unzulaessig"]),
|
|||
|
|
(0.0, zeile["ab"]), (breite, zeile["gierig"])):
|
|||
|
|
achse.annotate(f"{wert * 100:.0f}", (stelle + versatz, wert * 100),
|
|||
|
|
textcoords="offset points", xytext=(0, 3),
|
|||
|
|
ha="center", fontsize=7.5, color="#334155")
|
|||
|
|
|
|||
|
|
achse.set_xticks(stellen)
|
|||
|
|
achse.set_xticklabels(namen)
|
|||
|
|
achse.set_xlabel("Instanzgröße (Variablen × Nebenbedingungen)")
|
|||
|
|
achse.set_ylabel("Anteil bzw. Verlust (%)")
|
|||
|
|
achse.set_ylim(0, 118)
|
|||
|
|
achse.set_title(f"Gemessen an je {INSTANZEN} Zufallsinstanzen:\n"
|
|||
|
|
"Aufrunden wird mit der Größe hoffnungslos", fontsize=10.5)
|
|||
|
|
achse.grid(axis="y", linestyle=":", alpha=0.45)
|
|||
|
|
achse.legend(fontsize=8, loc="upper center", frameon=False, ncol=1)
|
|||
|
|
for seite in ("top", "right"):
|
|||
|
|
achse.spines[seite].set_visible(False)
|
|||
|
|
|
|||
|
|
|
|||
|
|
def zeichne(zeilen) -> None:
|
|||
|
|
figur, (links, rechts) = plt.subplots(1, 2, figsize=(11.0, 4.8))
|
|||
|
|
zeichne_handrechnung(links)
|
|||
|
|
zeichne_messung(rechts, zeilen)
|
|||
|
|
figur.tight_layout()
|
|||
|
|
|
|||
|
|
os.makedirs(BILDER, exist_ok=True)
|
|||
|
|
for endung in ("svg", "png"):
|
|||
|
|
pfad = os.path.join(BILDER, f"kap_milp_runden.{endung}")
|
|||
|
|
figur.savefig(pfad, format=endung, dpi=160,
|
|||
|
|
metadata={"Date": None} if endung == "svg" else None)
|
|||
|
|
print(f"geschrieben: {pfad}")
|
|||
|
|
plt.close(figur)
|
|||
|
|
|
|||
|
|
|
|||
|
|
if __name__ == "__main__":
|
|||
|
|
zeilen = messreihe()
|
|||
|
|
print(f"{'n x m':>8} | {'Aufrunden unzul.':>17} | {'Abrunden':>10} | {'gierig':>8}")
|
|||
|
|
print("-" * 54)
|
|||
|
|
for zeile in zeilen:
|
|||
|
|
print(f"{zeile['name']:>8} | {zeile['unzulaessig'] * 100:>16.1f} % | "
|
|||
|
|
f"{zeile['ab'] * 100:>9.1f} % | {zeile['gierig'] * 100:>7.1f} %")
|
|||
|
|
zeichne(zeilen)
|