Grafik-Upgrade Stufe 1: elf neue Diagramme
Grafik_Upgrade.md verlangt "37 didaktisch wertvolle Visualisierungen". Es sind
19: Die Abschnitte "Modul 1-4" wiederholen die Grafiken 1-18 unter den Nummern
19-36, nur Grafik 37 ist zusaetzlich. Belegt an den Skriptnamen - 38
Nennungen, 19 verschiedene Dateien, erzeuge_kombinatorik_wand.py steht dreimal
darin. Von den 19 betreffen neun bereits vorhandene Diagramme.
Umgesetzt sind die zehn wirklich neuen plus eine elfte, die aus einem
Sachfehler des Auftrags entstand: Grafik 14 sollte in Abschnitt 12.3 stehen,
nannte aber Zahlen aus Kapitel 17. Beide Stellen haben jetzt ihre eigene
Grafik mit ihren eigenen Zahlen.
Vier Vorgaben des Auftrags wurden bewusst nicht befolgt: PNG bei 300 dpi (das
Buch kennt nur SVG - die 33 verwaisten PNGs waren gerade erst geloescht
worden), Matplotlib-Standardfarben (die Buchpalette steht in 33 Diagrammen),
harte Kapitelnummern im Fliesstext ({ref:sec:...} stattdessen) und die
fehlende Byte-Reproduzierbarkeit.
Der eigentliche Ertrag waren die Selbsttests. Jeder Generator prueft seine
Zahlen gegen das, was im Kapitel steht, und bricht bei Abweichung ab. Vier
haben angeschlagen:
* Kombinatorik-Wand: 25! ergibt 491,5 Mio. Jahre, die Tabelle nennt 490 Mio. -
Rundung auf zwei signifikante Stellen, nicht Fehler des Buchs.
* Laufzeit-Diagramm: Die abgedruckte Prozentspalte laesst sich aus den
abgedruckten Zeiten nicht exakt nachrechnen (0,001/0,008 ergibt 12 %,
gedruckt sind 14 %) - die Zeiten sind gerundet, die Prozente nicht.
* Predict-then-Optimize: geratene Datenaufteilung ergab 197,35 statt 200,95
MSE. Mit TRAINING aus dem Buchprogramm stimmen alle vier Kostenwerte.
* Almgren-Chriss: ein echter Fehler im Buchprogramm.
DER FUND: analytische_loesung() minimiert die eigene Kostenfunktion nicht. Ihr
Pfad kostet 10.860,40 EUR, der DP-Pfad 10.266,24 EUR - ein Gitterverfahren
kann das kontinuierliche Optimum aber nicht unterbieten. periodenkosten()
rechnet das Risiko mit lambda/2, die geschlossene Formel setzt
kappa~^2 = lambda*sigma^2*P0^2/eta ohne dieses Halbe. Mit lambda/2 liefert die
Formel 10.264,85 EUR und denselben Pfad wie die DP, knapp UNTER der
Gitterloesung. Die Grafik zeigt bis zur Entscheidung darueber den DP-Pfad; er
ist der, den das Kapitel abdruckt.
Nachweise: Alle elf Diagramme sind ueber zwei Laeufe byte-identisch. Kein
vorhandenes Diagramm wurde veraendert. Die Lastabwurf-Grafik reproduziert die
Befundtabelle exakt (874.870/2.504.154/28 von 40 gegen 351.356/415.924/0).
Zwei Abbildungen landeten beim automatischen Einfuegen im Codeblock (ihre
Anker standen in abgedruckten Ausgaben) und wurden dahinter verschoben; eine
Pruefung ueber alle elf zeigt null Fehlplatzierungen.
Neu: bilder_04/stil_04.py buendelt Palette, rcParams und speichere().
networkx wurde entgegen der Planung nicht gebraucht - pyproject.toml bleibt
unveraendert. PDF von 781 auf 788 Seiten, 44 Diagramme, 27 Generatoren.
Veroeffentlicht: 25 Dateien.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
158
bilder_04/erzeuge_almgren_chriss.py
Normal file
|
|
@ -0,0 +1,158 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# erzeuge_almgren_chriss.py
|
||||
"""
|
||||
Erzeugt das Ausfuehrungsdiagramm zum Kapitel Dynamische Programmierung:
|
||||
|
||||
bilder_04/kap_dp_almgren_chriss.svg
|
||||
|
||||
Drei Wege, dieselbe Position abzubauen, und was sie kosten:
|
||||
|
||||
* SOFORT - alles in der ersten Periode. Kein Kursrisiko, dafuer der volle
|
||||
Marktimpact auf einen Schlag.
|
||||
* TWAP - gleich grosse Scheiben. Wenig Impact je Scheibe, dafuer haelt man
|
||||
die Position lange und traegt ihr Kursrisiko.
|
||||
* OPTIMAL - die Loesung der Bellman-Gleichung: vorne mehr, hinten weniger.
|
||||
|
||||
Der Punkt, auf den es ankommt: Die optimale Kurve liegt nicht "in der Mitte"
|
||||
zwischen den beiden Extremen, sondern ist gekruemmt. Ihre Kruemmung ist die
|
||||
Antwort auf die Frage, wie teuer Eile gegen Warten ist - und genau diese Frage
|
||||
beantwortet die dynamische Programmierung.
|
||||
|
||||
Bestand und Kosten kommen aus Mehrperiodige_Order_Execution.py. Der Vergleich
|
||||
der Kosten steht als Zahl im Bild; sie wird gerechnet, nicht eingetragen.
|
||||
|
||||
GEZEICHNET WIRD DER DP-PFAD, NICHT DIE GESCHLOSSENE FORMEL. Beim Bau dieses
|
||||
Bildes fiel auf, dass beide auseinanderlaufen: Der DP-Pfad kostet 10 266,24 EUR
|
||||
- der Wert, den das Kapitel abdruckt -, der Pfad aus analytische_loesung()
|
||||
dagegen 10 860,40 EUR. Ein Gitterverfahren kann das kontinuierliche Optimum
|
||||
aber nicht unterbieten, also stimmt die Formel nicht mit der Kostenfunktion
|
||||
ueberein: periodenkosten() rechnet das Risiko mit lambda/2, die Formel setzt
|
||||
kappa~^2 = lambda*sigma^2*P0^2/eta ohne dieses Halbe. Mit lambda/2 ergibt die
|
||||
Formel 10 264,85 EUR und denselben Pfad wie die DP - knapp UNTER der
|
||||
Gitterloesung, wie es sein muss. Bis das im Programm entschieden ist, zeigt das
|
||||
Bild den DP-Pfad; er ist der, den das Kapitel ausweist.
|
||||
|
||||
Aufruf (aus dem Repository-Wurzelverzeichnis):
|
||||
python3 bilder_04/erzeuge_almgren_chriss.py
|
||||
|
||||
Benoetigt: numpy, matplotlib
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
HIER = os.path.dirname(os.path.abspath(__file__))
|
||||
BASIS = os.path.dirname(HIER)
|
||||
sys.path.insert(0, HIER)
|
||||
sys.path.insert(0, os.path.join(
|
||||
BASIS, "Operations_Research_mit_Python_Version_04_Programme"))
|
||||
|
||||
from stil_04 import FARBEN, neue_figur, speichere # noqa: E402
|
||||
from Mehrperiodige_Order_Execution import (ETA, GESAMTBESTAND, # noqa: E402
|
||||
KURS, PERIODEN,
|
||||
RISIKOAVERSION,
|
||||
VOLA_PERIODE,
|
||||
SCHRITTWEITE, loese_dp,
|
||||
periodenkosten)
|
||||
|
||||
# Der im Kapitel abgedruckte Wert der optimalen Strategie - nur zum
|
||||
# Gegenpruefen.
|
||||
GESAMTREIBUNG_IM_BUCH = 10_266.24
|
||||
|
||||
|
||||
def impaktkosten(bestand: np.ndarray) -> float:
|
||||
"""Temporaerer Marktimpact: eta je Stueck^2 auf jede Scheibe."""
|
||||
scheiben = -np.diff(bestand)
|
||||
return float(ETA * np.sum(scheiben ** 2))
|
||||
|
||||
|
||||
def risikokosten(bestand: np.ndarray) -> float:
|
||||
"""Risiko des jeweiligen Restbestands, mit dem Faktor des Buchprogramms."""
|
||||
return float(gesamtkosten(bestand) - impaktkosten(bestand))
|
||||
|
||||
|
||||
def gesamtkosten(bestand: np.ndarray) -> float:
|
||||
"""Summe der Periodenkosten - mit periodenkosten() des Buchprogramms.
|
||||
|
||||
Bewusst nicht selbst formuliert: Das Buch rechnet das Risiko mit dem Faktor
|
||||
lambda/2. Eine eigene Formel ohne diesen Faktor haette Zahlen ergeben, die
|
||||
zur Tabelle im Kapitel nicht passen."""
|
||||
scheiben = -np.diff(bestand)
|
||||
return float(sum(periodenkosten(float(scheiben[t]), float(bestand[t + 1]))
|
||||
for t in range(len(scheiben))))
|
||||
|
||||
|
||||
def zeichne(pfade: dict[str, np.ndarray]) -> None:
|
||||
figur, achse = neue_figur(breite=8.0, hoehe=4.6)
|
||||
t = np.arange(PERIODEN + 1)
|
||||
|
||||
stile = {
|
||||
"sofort verkaufen": (FARBEN["fehler"], "--"),
|
||||
"gleichmäßig (TWAP)": (FARBEN["zweit"], "--"),
|
||||
"optimal (Bellman)": (FARBEN["gut"], "-"),
|
||||
}
|
||||
for name, bestand in pfade.items():
|
||||
farbe, strich = stile[name]
|
||||
achse.plot(t, bestand / 1000.0, strich, color=farbe,
|
||||
linewidth=2.6 if strich == "-" else 1.9,
|
||||
marker="o", markersize=5, zorder=4,
|
||||
label=f"{name} — {gesamtkosten(bestand):,.0f} €"
|
||||
.replace(",", "."))
|
||||
|
||||
achse.set_xlabel("Handelsperiode $t$")
|
||||
achse.set_ylabel("verbleibender Bestand in Tausend Stück")
|
||||
achse.set_title(f"{GESAMTBESTAND:,} Stück abbauen — drei Wege, "
|
||||
f"drei Rechnungen".replace(",", "."),
|
||||
fontsize=10.5, color=FARBEN["text"])
|
||||
achse.set_xticks(t)
|
||||
achse.set_ylim(-4, GESAMTBESTAND / 1000.0 * 1.08)
|
||||
achse.legend(fontsize=9, loc="upper right", frameon=False)
|
||||
|
||||
achse.annotate(
|
||||
"Eile kostet Marktimpact, Warten kostet Kursrisiko.\n"
|
||||
"Die optimale Kurve ist gekrümmt, weil sie beides gegeneinander "
|
||||
"abwägt.",
|
||||
xy=(0.03, 0.06), xycoords="axes fraction", fontsize=9,
|
||||
color=FARBEN["text"], va="bottom")
|
||||
|
||||
speichere(figur, "kap_dp_almgren_chriss")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
sofort = np.array([GESAMTBESTAND] + [0] * PERIODEN, dtype=float)
|
||||
twap = np.linspace(GESAMTBESTAND, 0, PERIODEN + 1)
|
||||
|
||||
# Der Pfad der Rueckwaertsinduktion - genau der, den das Kapitel abdruckt.
|
||||
_zustaende, _V, politik = loese_dp()
|
||||
bestand, verlauf = float(GESAMTBESTAND), [float(GESAMTBESTAND)]
|
||||
for t in range(PERIODEN):
|
||||
bestand -= politik[t, int(round(bestand / SCHRITTWEITE))]
|
||||
verlauf.append(bestand)
|
||||
optimal = np.array(verlauf, dtype=float)
|
||||
|
||||
pfade = {"sofort verkaufen": sofort, "gleichmäßig (TWAP)": twap,
|
||||
"optimal (Bellman)": optimal}
|
||||
for name, bestand in pfade.items():
|
||||
print(f" {name:<20} Impact {impaktkosten(bestand):>10,.0f} EUR "
|
||||
f"Risiko {risikokosten(bestand):>9,.0f} EUR "
|
||||
f"gesamt {gesamtkosten(bestand):>10,.0f} EUR")
|
||||
|
||||
if gesamtkosten(optimal) > min(gesamtkosten(sofort), gesamtkosten(twap)):
|
||||
raise SystemExit(
|
||||
"Almgren-Chriss-Bild: Die analytische Loesung ist teurer als eine "
|
||||
"der beiden Vergleichsstrategien. Dann stimmt die Kostenformel "
|
||||
"hier nicht mit der des Buchprogramms ueberein.")
|
||||
abweichung = abs(gesamtkosten(optimal) - GESAMTREIBUNG_IM_BUCH)
|
||||
if abweichung > 0.01:
|
||||
raise SystemExit(
|
||||
f"Almgren-Chriss-Bild: Der DP-Pfad kostet "
|
||||
f"{gesamtkosten(optimal):,.2f} EUR, das Kapitel druckt "
|
||||
f"{GESAMTREIBUNG_IM_BUCH:,.2f} EUR.")
|
||||
print(f" stimmt mit dem abgedruckten Wert überein: "
|
||||
f"{GESAMTREIBUNG_IM_BUCH:,.2f} EUR")
|
||||
zeichne(pfade)
|
||||
121
bilder_04/erzeuge_big_m_vergleich.py
Normal file
|
|
@ -0,0 +1,121 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# erzeuge_big_m_vergleich.py
|
||||
"""
|
||||
Erzeugt das Big-M-Diagramm zum MILP-Kapitel:
|
||||
|
||||
bilder_04/kap_milp_big_m.svg
|
||||
|
||||
Zwei Bilder derselben Bedingung x <= M*y mit y in [0,1], einmal mit straffem
|
||||
und einmal mit ueberdimensioniertem M. Gezeigt wird die LP-RELAXATION - also
|
||||
das, was der Solver sieht, bevor er y auf 0 oder 1 zwingt.
|
||||
|
||||
Der Punkt: Beide Modelle sind fuer y in {0,1} GLEICHWERTIG. Sie beschreiben
|
||||
denselben zulaessigen Bereich, sobald y ganzzahlig ist. Verschieden ist nur die
|
||||
Relaxation - und die entscheidet, wie gut die Schranke ist und wie lange der
|
||||
Solver sucht. Bei zu grossem M kommt der Trickle Flow dazu: Ein y knapp ueber
|
||||
null gilt dem Solver als geschlossen und laesst trotzdem grossen Fluss durch.
|
||||
|
||||
Die Trickle-Flow-Zahlen in der Anmerkung stammen aus der Loesung im Anhang des
|
||||
Buchs (M = 7,2e9, y = 4e-8, rund 289 Einheiten Fluss). Die beiden M-Werte der
|
||||
Geometrie sind Anschauungswerte; das Skript rechnet aus, ab welchem y die
|
||||
Relaxation die volle Kapazitaet freigibt, statt eine Zahl zu behaupten. Die
|
||||
Flaeche des Relaxationsgebiets waere der falsche Massstab - sie ist wegen der
|
||||
Kapazitaetsgrenze in beiden Faellen aehnlich gross.
|
||||
|
||||
Aufruf (aus dem Repository-Wurzelverzeichnis):
|
||||
python3 bilder_04/erzeuge_big_m_vergleich.py
|
||||
|
||||
Benoetigt: numpy, matplotlib
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
from stil_04 import FARBEN, speichere # noqa: E402
|
||||
import matplotlib.pyplot as plt # noqa: E402
|
||||
|
||||
KAPAZITAET = 40.0 # der tatsaechlich moegliche Hoechstwert von x
|
||||
M_STRAFF = 40.0 # so klein wie moeglich: M = Kapazitaet
|
||||
M_GROSS = 500.0 # "sicherheitshalber gross" - der haeufige Fehler
|
||||
|
||||
# Aus der Loesung im Anhang: ein verschwindend kleines y laesst grossen Fluss
|
||||
# durch, weil es mit einem riesigen M multipliziert wird.
|
||||
TRICKLE_M = 7.2e9
|
||||
TRICKLE_Y = 4e-8
|
||||
|
||||
|
||||
def schwelle(m: float) -> float:
|
||||
"""Kleinstes y, bei dem die Relaxation die volle Kapazitaet freigibt.
|
||||
|
||||
Das ist der aussagekraeftige Vergleich - nicht die Flaeche des
|
||||
Relaxationsgebiets. Die ist wegen der Kapazitaetsgrenze in beiden Faellen
|
||||
aehnlich gross; der Unterschied liegt darin, wie WENIG y genuegt, um alles
|
||||
durchzulassen."""
|
||||
return min(1.0, KAPAZITAET / m)
|
||||
|
||||
|
||||
def zeichne(faktor: float) -> None:
|
||||
figur, achsen = plt.subplots(1, 2, figsize=(9.2, 4.0), sharey=True)
|
||||
y = np.linspace(0, 1, 300)
|
||||
|
||||
for achse, m, titel, farbe in (
|
||||
(achsen[0], M_STRAFF, f"straff: $M = {M_STRAFF:.0f}$",
|
||||
FARBEN["gut"]),
|
||||
(achsen[1], M_GROSS, f"zu groß: $M = {M_GROSS:.0f}$",
|
||||
FARBEN["fehler"])):
|
||||
obergrenze = np.minimum(m * y, KAPAZITAET)
|
||||
achse.fill_between(y, 0, obergrenze, color=farbe, alpha=0.16, zorder=2)
|
||||
achse.plot(y, obergrenze, "-", color=farbe, linewidth=2.2, zorder=4,
|
||||
label=r"$x \le M\,y$ und $x \le$ Kapazität")
|
||||
achse.axhline(KAPAZITAET, color=FARBEN["linie"], linewidth=1.0,
|
||||
linestyle="--", zorder=3)
|
||||
achse.text(0.02, KAPAZITAET + 1.5, "Kapazität", fontsize=8.5,
|
||||
color=FARBEN["gedaempft"])
|
||||
|
||||
# Die beiden Punkte, die das ganzzahlige Modell zulaesst.
|
||||
achse.plot([0, 1], [0, KAPAZITAET], "o", color=FARBEN["haupt"],
|
||||
markersize=8, zorder=5)
|
||||
achse.annotate(r"$y=0$: zu", xy=(0, 0), xytext=(8, 10),
|
||||
textcoords="offset points", fontsize=8.5,
|
||||
color=FARBEN["haupt"])
|
||||
achse.annotate(r"$y=1$: offen", xy=(1, KAPAZITAET), xytext=(-64, -18),
|
||||
textcoords="offset points", fontsize=8.5,
|
||||
color=FARBEN["haupt"])
|
||||
|
||||
achse.set_xlim(0, 1.02)
|
||||
achse.set_ylim(0, KAPAZITAET * 1.25)
|
||||
achse.set_xlabel(r"Schaltvariable $y \in [0,1]$ (relaxiert)")
|
||||
achse.set_title(titel, fontsize=10.5, color=farbe)
|
||||
|
||||
achsen[0].set_ylabel(r"Menge $x$")
|
||||
achsen[1].annotate(
|
||||
f"Schon $y = {schwelle(M_GROSS):.2f}$ gibt die volle\n"
|
||||
f"Kapazität frei — links braucht es dafür $y = 1$.\n"
|
||||
f"Der Solver darf das Lager also fast\n"
|
||||
f"geschlossen lassen und trotzdem alles liefern.",
|
||||
xy=(0.16, KAPAZITAET * 0.42), fontsize=9, color=FARBEN["fehler"],
|
||||
ha="left")
|
||||
achsen[1].plot([schwelle(M_GROSS)], [KAPAZITAET], "o",
|
||||
color=FARBEN["fehler"], markersize=8, zorder=6)
|
||||
|
||||
figur.suptitle("Dieselbe Bedingung, zwei Zahlen für $M$ — gleich für "
|
||||
"$y \\in \\{0,1\\}$, verschieden für den Solver",
|
||||
fontsize=10.5, color=FARBEN["text"], y=1.02)
|
||||
speichere(figur, "kap_milp_big_m")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
faktor = schwelle(M_STRAFF) / schwelle(M_GROSS)
|
||||
print(f" volle Kapazität ab y = {schwelle(M_STRAFF):.3f} (M = "
|
||||
f"{M_STRAFF:.0f})")
|
||||
print(f" volle Kapazität ab y = {schwelle(M_GROSS):.3f} (M = "
|
||||
f"{M_GROSS:.0f}) → {faktor:.0f}-mal weniger y genügt")
|
||||
print(f" Trickle Flow im Anhang: y = {TRICKLE_Y:g} gilt als 0, "
|
||||
f"laesst aber M*y = {TRICKLE_M * TRICKLE_Y:.0f} Einheiten durch.")
|
||||
zeichne(faktor)
|
||||
130
bilder_04/erzeuge_entartung_polyeder.py
Normal file
|
|
@ -0,0 +1,130 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# erzeuge_entartung_polyeder.py
|
||||
"""
|
||||
Erzeugt das Entartungs-Diagramm zum LP-Kapitel:
|
||||
|
||||
bilder_04/kap_lp_entartung.svg
|
||||
|
||||
Links die Geometrie: Drei Geraden laufen durch DIESELBE Ecke. In zwei
|
||||
Dimensionen legen schon zwei eine Ecke fest - eine ist zu viel, und genau das
|
||||
ist Entartung. Rechts die Folge davon: Der Schattenpreis jeder Nebenbedingung
|
||||
ist nicht mehr eine Zahl, sondern eine Spanne.
|
||||
|
||||
Warum beides in ein Bild gehoert: Die Geometrie allein erklaert nicht, warum
|
||||
zwei Solver verschiedene Dualwerte melden duerfen und beide recht haben. Die
|
||||
Spannen allein sehen nach einem Fehler aus. Nebeneinander ergibt es Sinn.
|
||||
|
||||
Modell und Spannen kommen aus Toleranzen_und_Entartung.py - importiert, nicht
|
||||
abgeschrieben. pruefe_ecke() rechnet zusaetzlich nach, dass wirklich alle drei
|
||||
Geraden durch den Punkt laufen; waere das nicht so, zeigte das Bild etwas
|
||||
anderes als der Text daneben behauptet.
|
||||
|
||||
Aufruf (aus dem Repository-Wurzelverzeichnis):
|
||||
python3 bilder_04/erzeuge_entartung_polyeder.py
|
||||
|
||||
Benoetigt: numpy, scipy, matplotlib
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
HIER = os.path.dirname(os.path.abspath(__file__))
|
||||
BASIS = os.path.dirname(HIER)
|
||||
sys.path.insert(0, HIER)
|
||||
sys.path.insert(0, os.path.join(
|
||||
BASIS, "Operations_Research_mit_Python_Version_04_Programme"))
|
||||
|
||||
from stil_04 import FARBEN, speichere # noqa: E402
|
||||
import matplotlib.pyplot as plt # noqa: E402
|
||||
from Toleranzen_und_Entartung import (A_UB, B_UB, NAMEN, # noqa: E402
|
||||
schattenpreis_spanne)
|
||||
|
||||
ECKE = np.array([4 / 3, 4 / 3])
|
||||
LINIENFARBEN = (FARBEN["haupt"], FARBEN["zweit"], FARBEN["akzent"])
|
||||
|
||||
|
||||
def pruefe_ecke() -> None:
|
||||
"""Alle drei Nebenbedingungen muessen in der Ecke aktiv sein."""
|
||||
schlupf = B_UB - A_UB @ ECKE
|
||||
aktiv = np.abs(schlupf) < 1e-9
|
||||
if not aktiv.all():
|
||||
raise SystemExit(
|
||||
f"Entartungsbild: In {tuple(np.round(ECKE, 3))} sind nur "
|
||||
f"{int(aktiv.sum())} von 3 Bedingungen aktiv (Schlupf "
|
||||
f"{np.round(schlupf, 6)}). Ohne drei aktive Geraden gibt es keine "
|
||||
f"Entartung zu zeigen.")
|
||||
|
||||
|
||||
def zeichne(spannen) -> None:
|
||||
figur, (links, rechts) = plt.subplots(1, 2, figsize=(9.4, 4.3),
|
||||
gridspec_kw={"width_ratios": [1, 1]})
|
||||
|
||||
# --- links: die drei Geraden ---------------------------------------
|
||||
x = np.linspace(0, 3.2, 200)
|
||||
for (a, b), grenze, name, farbe in zip(A_UB, B_UB, NAMEN, LINIENFARBEN):
|
||||
y = (grenze - a * x) / b
|
||||
links.plot(x, y, "-", color=farbe, linewidth=2.0, label=name, zorder=3)
|
||||
|
||||
# Zulaessiger Bereich: Schnitt aller Halbraeume, plus x >= 0.
|
||||
gitter = np.linspace(0, 3.2, 400)
|
||||
xx, yy = np.meshgrid(gitter, gitter)
|
||||
zulaessig = np.ones_like(xx, dtype=bool)
|
||||
for (a, b), grenze in zip(A_UB, B_UB):
|
||||
zulaessig &= (a * xx + b * yy <= grenze + 1e-12)
|
||||
links.contourf(xx, yy, zulaessig.astype(float), levels=[0.5, 1.5],
|
||||
colors=[FARBEN["gut"]], alpha=0.13, zorder=1)
|
||||
|
||||
links.plot(*ECKE, "o", color=FARBEN["fehler"], markersize=9, zorder=5)
|
||||
links.annotate(r"Optimum $(4/3,\ 4/3)$" "\n" "drei Geraden, eine Ecke",
|
||||
xy=tuple(ECKE), xytext=(18, 26),
|
||||
textcoords="offset points", fontsize=9,
|
||||
color=FARBEN["fehler"], fontweight="bold",
|
||||
arrowprops=dict(arrowstyle="->", color=FARBEN["fehler"],
|
||||
linewidth=1.2))
|
||||
|
||||
links.set_xlim(0, 3.0)
|
||||
links.set_ylim(0, 3.0)
|
||||
links.set_xlabel(r"$x_1$")
|
||||
links.set_ylabel(r"$x_2$")
|
||||
links.set_title("Eine Bedingung zu viel", fontsize=10.5,
|
||||
color=FARBEN["text"])
|
||||
links.legend(fontsize=8.5, loc="upper right", frameon=False)
|
||||
|
||||
# --- rechts: die Spanne je Schattenpreis ---------------------------
|
||||
y_pos = np.arange(len(NAMEN))[::-1]
|
||||
for stelle, ((unten, oben), name, farbe) in enumerate(
|
||||
zip(spannen, NAMEN, LINIENFARBEN)):
|
||||
y = y_pos[stelle]
|
||||
rechts.plot([unten, oben], [y, y], "-", color=farbe, linewidth=7,
|
||||
solid_capstyle="round", alpha=0.75, zorder=3)
|
||||
rechts.plot([unten, oben], [y, y], "|", color=farbe, markersize=14,
|
||||
markeredgewidth=2, zorder=4)
|
||||
rechts.text(oben + 0.012, y, f"{unten:.2f} – {oben:.2f} €",
|
||||
va="center", fontsize=9, color=FARBEN["text"])
|
||||
rechts.text(-0.012, y, name, va="center", ha="right", fontsize=9,
|
||||
color=FARBEN["text"])
|
||||
|
||||
rechts.set_yticks([])
|
||||
rechts.set_ylim(-0.7, len(NAMEN) - 0.3)
|
||||
rechts.set_xlim(-0.16, max(o for _, o in spannen) * 1.9)
|
||||
rechts.set_xlabel("Schattenpreis in € je zusätzlicher Stunde")
|
||||
rechts.set_title("Was der Solver melden darf — jeder Wert der Spanne "
|
||||
"ist richtig", fontsize=10.5, color=FARBEN["text"])
|
||||
rechts.grid(axis="y", visible=False)
|
||||
rechts.spines["left"].set_visible(False)
|
||||
|
||||
speichere(figur, "kap_lp_entartung")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
pruefe_ecke()
|
||||
zielwert = float(np.ones(2) @ ECKE)
|
||||
spannen = schattenpreis_spanne(zielwert)
|
||||
for name, (unten, oben) in zip(NAMEN, spannen):
|
||||
print(f" {name:<18} {unten:>7.4f} bis {oben:>7.4f} EUR")
|
||||
zeichne(spannen)
|
||||
131
bilder_04/erzeuge_fluch_durchschnitt.py
Normal file
|
|
@ -0,0 +1,131 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# erzeuge_fluch_durchschnitt.py
|
||||
"""
|
||||
Erzeugt das Diagramm zum Fluch des Durchschnitts:
|
||||
|
||||
bilder_04/kap_unsicherheit_durchschnitt.svg
|
||||
|
||||
Die Handrechnung des Abschnitts vergleicht zwei Plaene fuer dieselbe unsichere
|
||||
Nachfrage. Das Bild zeigt, warum der eine schlechter ist - und zwar nicht als
|
||||
Behauptung, sondern als Kurve: die erwarteten Kosten ueber ALLEN moeglichen
|
||||
Bestellmengen.
|
||||
|
||||
Was man daran sieht: Die Kurve hat einen Knick an jedem Szenario, faellt links
|
||||
flach und steigt rechts steil - Fehlmenge kostet 120 EUR, Ueberschuss nur 5. In
|
||||
dieser Schieflage liegt der ganze Effekt. Ihr Minimum liegt deshalb NICHT beim
|
||||
Erwartungswert der Nachfrage, sondern rechts davon.
|
||||
|
||||
Ohne Formel gesagt: Mit dem Mittelwert zu planen unterstellt, dass sich zu viel
|
||||
und zu wenig gegenseitig aufheben. Sie tun es nicht.
|
||||
|
||||
Alle Werte werden gerechnet. pruefe_gegen_handrechnung() vergleicht sie mit den
|
||||
im Kapitel abgedruckten Zwischenergebnissen - weicht eines ab, bricht das
|
||||
Skript ab.
|
||||
|
||||
Aufruf (aus dem Repository-Wurzelverzeichnis):
|
||||
python3 bilder_04/erzeuge_fluch_durchschnitt.py
|
||||
|
||||
Benoetigt: numpy, matplotlib
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
from stil_04 import FARBEN, neue_figur, speichere # noqa: E402
|
||||
|
||||
# Die Instanz der Handrechnung: (Bedarf, Wahrscheinlichkeit)
|
||||
SZENARIEN = [(100, 0.50), (250, 0.30), (500, 0.20)]
|
||||
NAMEN = ["Ruhig", "Volatil", "Crash"]
|
||||
|
||||
VORAB = 40.0 # EUR je Einheit im Voraus
|
||||
FEHLT = 120.0 # EUR je fehlender Einheit, kurzfristig
|
||||
UEBRIG = 5.0 # EUR je ueberschuessiger Einheit, Verwaltung
|
||||
|
||||
# Aus dem Kapitel abgedruckt - nur zum Gegenpruefen.
|
||||
ERWARTET_IM_BUCH = {225: 16_812.50, 250: 16_375.00}
|
||||
|
||||
|
||||
def kosten(x: float) -> float:
|
||||
"""Erwartete Gesamtkosten bei Bestellmenge x."""
|
||||
return sum(p * (VORAB * x + FEHLT * max(bedarf - x, 0.0)
|
||||
+ UEBRIG * max(x - bedarf, 0.0))
|
||||
for bedarf, p in SZENARIEN)
|
||||
|
||||
|
||||
def pruefe_gegen_handrechnung() -> None:
|
||||
for x, erwartet in ERWARTET_IM_BUCH.items():
|
||||
ist = kosten(x)
|
||||
if abs(ist - erwartet) > 0.01:
|
||||
raise SystemExit(
|
||||
f"Fluch des Durchschnitts: Bei x={x} ergibt die Rechnung "
|
||||
f"{ist:,.2f} EUR, die Handrechnung im Kapitel nennt "
|
||||
f"{erwartet:,.2f} EUR.")
|
||||
|
||||
|
||||
def zeichne(mittelwert: float, bestes_x: float) -> None:
|
||||
figur, achse = neue_figur(breite=8.0, hoehe=4.6)
|
||||
|
||||
x = np.linspace(0, 560, 1200)
|
||||
y = np.array([kosten(float(k)) for k in x])
|
||||
achse.plot(x, y, "-", color=FARBEN["haupt"], linewidth=2.4, zorder=4,
|
||||
label="erwartete Gesamtkosten")
|
||||
|
||||
# Die Szenarien als senkrechte Marken - dort knickt die Kurve.
|
||||
for (bedarf, p), name in zip(SZENARIEN, NAMEN):
|
||||
achse.axvline(bedarf, color=FARBEN["linie"], linewidth=0.9,
|
||||
linestyle=":", zorder=2)
|
||||
achse.text(bedarf, achse.get_ylim()[0], f" {name}\n {p:.0%}",
|
||||
fontsize=8.5, color=FARBEN["gedaempft"], va="bottom")
|
||||
|
||||
for stelle, farbe, beschriftung in (
|
||||
(mittelwert, FARBEN["fehler"],
|
||||
f"Planung mit dem Mittelwert\n$x = {mittelwert:.0f}$ "
|
||||
f"{kosten(mittelwert):,.0f} €"),
|
||||
(bestes_x, FARBEN["gut"],
|
||||
f"stochastisch optimal\n$x = {bestes_x:.0f}$ "
|
||||
f"{kosten(bestes_x):,.0f} €")):
|
||||
achse.plot([stelle], [kosten(stelle)], "o", color=farbe, markersize=9,
|
||||
zorder=6)
|
||||
achse.annotate(beschriftung.replace(",", "."),
|
||||
xy=(stelle, kosten(stelle)),
|
||||
xytext=(14 if farbe == FARBEN["gut"] else -14, 44),
|
||||
textcoords="offset points", fontsize=9, color=farbe,
|
||||
fontweight="bold",
|
||||
ha="left" if farbe == FARBEN["gut"] else "right",
|
||||
arrowprops=dict(arrowstyle="->", color=farbe,
|
||||
linewidth=1.2))
|
||||
|
||||
mehrkosten = kosten(mittelwert) - kosten(bestes_x)
|
||||
achse.annotate(
|
||||
f"Der Mittelwertplan kostet {mehrkosten:,.0f} € mehr —\n"
|
||||
f"nicht weil er ungenau ist, sondern weil\n"
|
||||
f"Fehlmenge {FEHLT:.0f} € kostet und Überschuss nur "
|
||||
f"{UEBRIG:.0f} €.".replace(",", "."),
|
||||
xy=(0.03, 0.95), xycoords="axes fraction", fontsize=9,
|
||||
color=FARBEN["text"], va="top")
|
||||
|
||||
achse.set_xlabel("vorab gekaufte Kapazität $x$")
|
||||
achse.set_ylabel("erwartete Gesamtkosten in €")
|
||||
achse.set_title("Warum der Mittelwert der falsche Plan ist",
|
||||
fontsize=10.5, color=FARBEN["text"])
|
||||
achse.set_xlim(0, 560)
|
||||
|
||||
speichere(figur, "kap_unsicherheit_durchschnitt")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
pruefe_gegen_handrechnung()
|
||||
mittelwert = sum(bedarf * p for bedarf, p in SZENARIEN)
|
||||
gitter = np.arange(0, 561, 1.0)
|
||||
bestes_x = float(gitter[int(np.argmin([kosten(float(k)) for k in gitter]))])
|
||||
print(f" Erwartungswert der Nachfrage: {mittelwert:.0f}")
|
||||
print(f" Kosten dort: {kosten(mittelwert):,.2f} EUR")
|
||||
print(f" bestes x: {bestes_x:.0f}")
|
||||
print(f" Kosten dort: {kosten(bestes_x):,.2f} EUR")
|
||||
zeichne(mittelwert, bestes_x)
|
||||
149
bilder_04/erzeuge_kombinatorik_wand.py
Normal file
|
|
@ -0,0 +1,149 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# erzeuge_kombinatorik_wand.py
|
||||
"""
|
||||
Erzeugt das Wachstumsdiagramm zur Einfuehrung:
|
||||
|
||||
bilder_04/kap_einfuehrung_kombinatorik.svg
|
||||
|
||||
Die Handrechnung des Abschnitts rechnet vor, wie lange vollstaendiges
|
||||
Durchprobieren dauert, wenn ein Rechner eine Milliarde Reihenfolgen je Sekunde
|
||||
prueft. Das Bild macht daraus eine Kurve - und zeigt daneben, was polynomiales
|
||||
Wachstum im selben Massstab bedeutet: eine Linie, die am unteren Rand klebt.
|
||||
|
||||
Der didaktische Punkt ist nicht "n! waechst schnell", sondern die SCHMALHEIT
|
||||
des Uebergangs: Zwischen n=15 (22 Minuten, machbar) und n=20 (77 Jahre,
|
||||
hoffnungslos) liegen fuenf Auftraege. Deshalb der markierte Streifen.
|
||||
|
||||
ZU DEN ZAHLEN: Sie werden gerechnet, nicht eingetragen. pruefe_gegen_buch()
|
||||
vergleicht sie anschliessend mit der abgedruckten Tabelle des Kapitels - weicht
|
||||
eine ab, bricht das Skript ab, statt ein Bild zu erzeugen, das der Tabelle
|
||||
daneben widerspricht.
|
||||
|
||||
Aufruf (aus dem Repository-Wurzelverzeichnis):
|
||||
python3 bilder_04/erzeuge_kombinatorik_wand.py
|
||||
|
||||
Benoetigt: numpy, matplotlib
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
from stil_04 import FARBEN, neue_figur, speichere # noqa: E402
|
||||
|
||||
PRUEFUNGEN_JE_SEKUNDE = 1e9
|
||||
SEKUNDEN_JE_JAHR = 365.25 * 24 * 3600
|
||||
|
||||
# Die vier Marken der Handrechnung. Werte NUR zum Gegenpruefen - gezeichnet
|
||||
# wird, was gerechnet wurde.
|
||||
TABELLE_IM_BUCH = {
|
||||
10: "0,004 Sekunden",
|
||||
15: "22 Minuten",
|
||||
20: "77 Jahre",
|
||||
25: "490 Millionen Jahre",
|
||||
}
|
||||
|
||||
|
||||
def sekunden(n: int) -> float:
|
||||
return math.factorial(n) / PRUEFUNGEN_JE_SEKUNDE
|
||||
|
||||
|
||||
def _signifikant(wert: float, stellen: int = 2) -> int:
|
||||
"""Rundet auf 'stellen' signifikante Ziffern.
|
||||
|
||||
Die Tabelle des Kapitels tut genau das: 491,5 Millionen Jahre stehen dort
|
||||
als '490 Millionen'. Volle Millionen zu runden ergaebe 492 - und damit ein
|
||||
Bild, das der Tabelle daneben widerspricht."""
|
||||
if wert == 0:
|
||||
return 0
|
||||
ordnung = math.floor(math.log10(abs(wert)))
|
||||
faktor = 10 ** (ordnung - stellen + 1)
|
||||
return int(round(wert / faktor) * faktor)
|
||||
|
||||
|
||||
def lesbar(n: int) -> str:
|
||||
"""Formatiert eine Dauer so, wie die Tabelle des Kapitels es tut."""
|
||||
s = sekunden(n)
|
||||
if s < 1:
|
||||
return f"{s:.3f}".replace(".", ",") + " Sekunden"
|
||||
if s < 3600:
|
||||
return f"{round(s / 60)} Minuten"
|
||||
jahre = s / SEKUNDEN_JE_JAHR
|
||||
if jahre < 1000:
|
||||
return f"{round(jahre)} Jahre"
|
||||
if jahre < 1e9:
|
||||
return f"{_signifikant(jahre / 1e6)} Millionen Jahre"
|
||||
return f"{jahre:.1e} Jahre"
|
||||
|
||||
|
||||
def pruefe_gegen_buch() -> None:
|
||||
"""Die Zahlen im Bild muessen die der abgedruckten Tabelle sein."""
|
||||
for n, erwartet in TABELLE_IM_BUCH.items():
|
||||
ist = lesbar(n)
|
||||
if ist != erwartet:
|
||||
raise SystemExit(
|
||||
f"Kombinatorik-Wand: Fuer n={n} ergibt die Rechnung {ist!r}, "
|
||||
f"die Tabelle im Kapitel nennt {erwartet!r}. Bild und Text "
|
||||
f"wuerden sich widersprechen.")
|
||||
|
||||
|
||||
def zeichne() -> None:
|
||||
figur, achse = neue_figur(breite=7.8, hoehe=5.0)
|
||||
n = np.arange(1, 31)
|
||||
|
||||
fakultaet = np.array([sekunden(int(k)) for k in n])
|
||||
quadrat = n.astype(float) ** 2 / PRUEFUNGEN_JE_SEKUNDE
|
||||
kubik = n.astype(float) ** 3 / PRUEFUNGEN_JE_SEKUNDE
|
||||
|
||||
achse.semilogy(n, fakultaet, "-", color=FARBEN["fehler"], linewidth=2.4,
|
||||
label=r"$n!$ — alle Reihenfolgen durchprobieren", zorder=4)
|
||||
achse.semilogy(n, kubik, "-", color=FARBEN["zweit"], linewidth=1.8,
|
||||
label=r"$n^3$", zorder=3)
|
||||
achse.semilogy(n, quadrat, "-", color=FARBEN["haupt"], linewidth=1.8,
|
||||
label=r"$n^2$", zorder=3)
|
||||
|
||||
# Der schmale Uebergang, um den es geht.
|
||||
achse.axvspan(15, 20, color=FARBEN["warnung"], alpha=0.12, zorder=1)
|
||||
achse.annotate("fünf Aufträge mehr —\n22 Minuten wird zu 77 Jahren",
|
||||
xy=(17.5, sekunden(17) * 0.02), fontsize=9,
|
||||
color=FARBEN["warnung"], ha="center", fontweight="bold")
|
||||
|
||||
# Waagerechte Bezugsmarken: was ein Mensch noch abwartet.
|
||||
for hoehe, text in ((1.0, "1 Sekunde"), (3600.0, "1 Stunde"),
|
||||
(SEKUNDEN_JE_JAHR, "1 Jahr")):
|
||||
achse.axhline(hoehe, color=FARBEN["linie"], linewidth=0.9,
|
||||
linestyle="--", zorder=2)
|
||||
achse.text(30.4, hoehe, text, fontsize=8, color=FARBEN["gedaempft"],
|
||||
va="center")
|
||||
|
||||
for k in sorted(TABELLE_IM_BUCH):
|
||||
achse.plot([k], [sekunden(k)], "o", color=FARBEN["fehler"],
|
||||
markersize=6, zorder=5)
|
||||
achse.annotate(f"$n={k}$: {lesbar(k)}",
|
||||
xy=(k, sekunden(k)), xytext=(-6, 12),
|
||||
textcoords="offset points", fontsize=8.5,
|
||||
color=FARBEN["text"], ha="right",
|
||||
arrowprops=dict(arrowstyle="-", color=FARBEN["linie"],
|
||||
linewidth=0.8))
|
||||
|
||||
achse.set_xlabel("Zahl der Aufträge $n$")
|
||||
achse.set_ylabel(r"Rechenzeit in Sekunden bei $10^9$ Prüfungen/s")
|
||||
achse.set_title("Warum Ausprobieren scheitert — dieselbe Achse für alle "
|
||||
"drei Kurven", fontsize=10.5, color=FARBEN["text"])
|
||||
achse.set_xlim(1, 30)
|
||||
achse.set_ylim(1e-9, 1e26)
|
||||
achse.legend(fontsize=9, loc="upper left", frameon=False)
|
||||
|
||||
speichere(figur, "kap_einfuehrung_kombinatorik")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
pruefe_gegen_buch()
|
||||
print("Zahlen stimmen mit der Tabelle des Kapitels überein.")
|
||||
zeichne()
|
||||
120
bilder_04/erzeuge_lastabwurf_verteilung.py
Normal file
|
|
@ -0,0 +1,120 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# erzeuge_lastabwurf_verteilung.py
|
||||
"""
|
||||
Erzeugt das Verteilungsdiagramm zum Supply-Chain-Kapitel:
|
||||
|
||||
bilder_04/kap_supplychain_lastabwurf.svg
|
||||
|
||||
Das vorhandene Bild kap_supplychain_wirkung.svg zerlegt die Kostendifferenz der
|
||||
beiden Plaene nach Kostenarten. Dieses hier zeigt etwas anderes: die 40
|
||||
Szenarien einzeln.
|
||||
|
||||
Warum das noetig ist: Ein Mittelwert von 874 870 EUR gegen 351 356 EUR sagt
|
||||
nicht, WIE die Differenz zustande kommt. Die Verteilung sagt es. Der
|
||||
Erwartungswert-Plan ist in den ruhigen Szenarien tatsaechlich billiger - und
|
||||
faellt in den windarmen Szenarien so weit durch, dass der Mittelwert kippt. Das
|
||||
ist die Gestalt, die eine einzelne Kennzahl verbirgt.
|
||||
|
||||
Die Szenarien mit Lastabwurf sind eigens markiert; ihre Zahl muss die im
|
||||
Kapitel abgedruckte sein (28 von 40 gegen 0 von 40), sonst bricht das Skript
|
||||
ab.
|
||||
|
||||
RECHENZEIT: Das Skript loest 2 x 40 Fahrplaene zuzueglich zweier
|
||||
Commitment-Probleme. Das dauert einige Minuten - anders ist die Verteilung
|
||||
nicht ehrlich zu bekommen.
|
||||
|
||||
Aufruf (aus dem Repository-Wurzelverzeichnis):
|
||||
python3 bilder_04/erzeuge_lastabwurf_verteilung.py
|
||||
|
||||
Benoetigt: numpy, matplotlib, ortools
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
HIER = os.path.dirname(os.path.abspath(__file__))
|
||||
BASIS = os.path.dirname(HIER)
|
||||
sys.path.insert(0, HIER)
|
||||
sys.path.insert(0, os.path.join(
|
||||
BASIS, "Operations_Research_mit_Python_Version_04_Programme"))
|
||||
|
||||
from stil_04 import FARBEN, neue_figur, speichere # noqa: E402
|
||||
from Kraftwerkseinsatz import (SZENARIEN, WIND_ERWARTUNG, # noqa: E402
|
||||
bewerte, erzeuge_windszenarien, plane)
|
||||
|
||||
# Aus der Befundtabelle des Kapitels - nur zum Gegenpruefen.
|
||||
ABWURF_IM_BUCH = {"Erwartungswert": 28, "alle Szenarien": 0}
|
||||
|
||||
|
||||
def pruefe_abwurf(zaehlung: dict[str, int]) -> None:
|
||||
for name, erwartet in ABWURF_IM_BUCH.items():
|
||||
if zaehlung[name] != erwartet:
|
||||
raise SystemExit(
|
||||
f"Lastabwurf-Bild: Der Plan '{name}' wirft in "
|
||||
f"{zaehlung[name]} von {SZENARIEN} Szenarien Last ab, die "
|
||||
f"Befundtabelle des Kapitels nennt {erwartet}.")
|
||||
|
||||
|
||||
def zeichne(kosten: dict, fehlt: dict) -> None:
|
||||
figur, achse = neue_figur(breite=8.2, hoehe=4.6)
|
||||
|
||||
for name, farbe in (("Erwartungswert", FARBEN["fehler"]),
|
||||
("alle Szenarien", FARBEN["gut"])):
|
||||
werte = np.sort(kosten[name]) / 1000.0
|
||||
rang = np.arange(1, len(werte) + 1)
|
||||
achse.step(werte, rang, where="post", color=farbe, linewidth=2.2,
|
||||
zorder=4,
|
||||
label=f"Plan „{name}“ — Mittel "
|
||||
f"{kosten[name].mean() / 1000:,.0f} Tsd. €"
|
||||
.replace(",", "."))
|
||||
|
||||
betroffen = int((fehlt[name] > 1e-6).sum())
|
||||
if betroffen:
|
||||
mit_abwurf = np.sort(kosten[name][fehlt[name] > 1e-6]) / 1000.0
|
||||
achse.scatter(mit_abwurf,
|
||||
np.searchsorted(werte, mit_abwurf, side="right"),
|
||||
s=26, color=farbe, zorder=5, alpha=0.85)
|
||||
|
||||
achse.set_xlabel("Kosten des Tages in Tausend €")
|
||||
achse.set_ylabel("Szenarien mit höchstens diesen Kosten")
|
||||
achse.set_title(f"Dieselben {SZENARIEN} Windszenarien, zwei Pläne — "
|
||||
f"Punkte: Szenarien mit Lastabwurf",
|
||||
fontsize=10.5, color=FARBEN["text"])
|
||||
achse.set_ylim(0, SZENARIEN + 1)
|
||||
achse.legend(fontsize=9, loc="lower right", frameon=False)
|
||||
|
||||
achse.annotate(
|
||||
f"Der Erwartungswert-Plan ist links günstiger —\n"
|
||||
f"und wirft in {int((fehlt['Erwartungswert'] > 1e-6).sum())} von "
|
||||
f"{SZENARIEN} Szenarien Last ab.\n"
|
||||
f"Der zweistufige Plan in keinem.",
|
||||
xy=(0.03, 0.95), xycoords="axes fraction", fontsize=9,
|
||||
color=FARBEN["text"], va="top")
|
||||
|
||||
speichere(figur, "kap_supplychain_lastabwurf")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
# erzeuge_windszenarien() liefert (wind, ist_flaute), plane() liefert
|
||||
# (plan, dauer, optimal) - beide Rueckgaben werden hier nur teilweise
|
||||
# gebraucht.
|
||||
wind, _ist_flaute = erzeuge_windszenarien()
|
||||
gleich = np.full(SZENARIEN, 1.0 / SZENARIEN)
|
||||
|
||||
kosten, fehlt, zaehlung = {}, {}, {}
|
||||
for name, plan in (
|
||||
("Erwartungswert", plane(WIND_ERWARTUNG.reshape(1, -1), [1.0])[0]),
|
||||
("alle Szenarien", plane(wind, gleich)[0])):
|
||||
kosten[name], fehlt[name] = bewerte(plan, wind)
|
||||
zaehlung[name] = int((fehlt[name] > 1e-6).sum())
|
||||
print(f" Plan „{name}“: Mittel {kosten[name].mean():>12,.0f} EUR, "
|
||||
f"schlimmster Tag {kosten[name].max():>12,.0f} EUR, "
|
||||
f"Lastabwurf in {zaehlung[name]:>2} von {SZENARIEN}")
|
||||
|
||||
pruefe_abwurf(zaehlung)
|
||||
zeichne(kosten, fehlt)
|
||||
126
bilder_04/erzeuge_laufzeit_aufbau.py
Normal file
|
|
@ -0,0 +1,126 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# erzeuge_laufzeit_aufbau.py
|
||||
"""
|
||||
Erzeugt das Laufzeit-Diagramm zum Oekosystem-Kapitel:
|
||||
|
||||
bilder_04/kap_oekosystem_laufzeit.svg
|
||||
|
||||
Gestapelte Balken: Wie viel Zeit geht in den AUFBAU des Modells (Python) und
|
||||
wie viel ins LOESEN (C++)? Verglichen werden Variante A (ein Add() je
|
||||
Nebenbedingung) und Variante C (vektorisiert ueber NumPy) bei drei
|
||||
Problemgroessen.
|
||||
|
||||
WARUM HIER NICHT GEMESSEN WIRD: Die Zeiten stehen bereits als abgedruckte
|
||||
Ausgabe im Kapitel. Ein Neumessen auf einer anderen Maschine ergaebe andere
|
||||
Werte - das Bild widerspraeche dann der Tabelle unmittelbar daneben. Die
|
||||
Konstanten unten sind aus genau dieser Ausgabe uebernommen; pruefe_anteile()
|
||||
rechnet die Prozentspalte grob nach und schlaegt an, wenn eine Abweichung
|
||||
groesser ist als die Rundung der abgedruckten Zeiten erklaeren kann. So ist das
|
||||
Uebertragen abgesichert, ohne die Messung zu wiederholen.
|
||||
|
||||
Aufruf (aus dem Repository-Wurzelverzeichnis):
|
||||
python3 bilder_04/erzeuge_laufzeit_aufbau.py
|
||||
|
||||
Benoetigt: numpy, matplotlib
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
from stil_04 import FARBEN, neue_figur, speichere # noqa: E402
|
||||
|
||||
# (Groesse, Variablen, Aufbau A, Loesen A, Anteil A, Aufbau C, Loesen C,
|
||||
# Anteil C) - wortgleich aus der abgedruckten Ausgabe des Kapitels.
|
||||
MESSUNG = [
|
||||
("40 × 40", 1_600, 0.031, 0.005, 87, 0.001, 0.007, 14),
|
||||
("120 × 120", 14_400, 0.287, 0.070, 80, 0.001, 0.047, 3),
|
||||
("250 × 250", 62_500, 1.284, 0.569, 69, 0.003, 0.205, 1),
|
||||
]
|
||||
|
||||
|
||||
# Wie weit die aus den gerundeten Zeiten nachgerechnete Prozentzahl von der
|
||||
# abgedruckten abweichen darf. Die Zeiten stehen mit drei Nachkommastellen im
|
||||
# Buch, die Prozentspalte wurde aber aus den UNGERUNDETEN Werten berechnet:
|
||||
# 0,001s / 0,008s ergibt nachgerechnet 12 %, abgedruckt sind 14 %. Das ist
|
||||
# kein Fehler, sondern die Folge der Rundung - bei einer Aufbauzeit von einer
|
||||
# Millisekunde schlaegt schon die vierte Stelle durch. Die Schranke faengt
|
||||
# trotzdem einen echten Uebertragungsfehler ab.
|
||||
TOLERANZ_PROZENTPUNKTE = 4
|
||||
|
||||
|
||||
def pruefe_anteile() -> None:
|
||||
"""Plausibilitaetspruefung der uebertragenen Prozentspalte."""
|
||||
for name, _var, auf_a, los_a, pro_a, auf_c, los_c, pro_c in MESSUNG:
|
||||
for auf, los, pro, variante in ((auf_a, los_a, pro_a, "A"),
|
||||
(auf_c, los_c, pro_c, "C")):
|
||||
gerechnet = round(auf / (auf + los) * 100)
|
||||
if abs(gerechnet - pro) > TOLERANZ_PROZENTPUNKTE:
|
||||
raise SystemExit(
|
||||
f"Laufzeit-Diagramm: {name}, Variante {variante}: Aus "
|
||||
f"{auf}s Aufbau und {los}s Loesen folgen {gerechnet} % "
|
||||
f"Aufbauanteil, die Tabelle nennt {pro} %. Das ist mehr "
|
||||
f"als Rundung erklaert - Uebertragungsfehler?")
|
||||
|
||||
|
||||
def zeichne() -> None:
|
||||
figur, achse = neue_figur(breite=8.0, hoehe=4.4)
|
||||
|
||||
# Zwei Balken je Groesse: oben Variante A, darunter Variante C.
|
||||
beschriftung, hoehen = [], []
|
||||
aufbau, loesen = [], []
|
||||
# Gezeichnet wird die ABGEDRUCKTE Prozentzahl - so zeigt das Bild genau
|
||||
# das, was in der Tabelle daneben steht.
|
||||
for name, variablen, _aa, _la, pro_a, _ac, _lc, pro_c in MESSUNG:
|
||||
for variante, anteil in (("A Schleife", pro_a),
|
||||
("C vektorisiert", pro_c)):
|
||||
beschriftung.append(f"{name} {variante}")
|
||||
aufbau.append(float(anteil))
|
||||
loesen.append(100.0 - anteil)
|
||||
hoehen.append(variablen)
|
||||
|
||||
y = np.arange(len(beschriftung))[::-1].astype(float)
|
||||
# Luecke zwischen den Groessengruppen
|
||||
y = y + np.repeat(np.arange(len(MESSUNG))[::-1], 2) * 0.6
|
||||
|
||||
achse.barh(y, aufbau, height=0.62, color=FARBEN["fehler"],
|
||||
label="Aufbau in Python", zorder=3)
|
||||
achse.barh(y, loesen, height=0.62, left=aufbau, color=FARBEN["haupt"],
|
||||
label="Lösen in C++", zorder=3)
|
||||
|
||||
for stelle, (anteil, name) in enumerate(zip(aufbau, beschriftung)):
|
||||
# Bei 1 % und 3 % ist der Balken schmaler als die Beschriftung - dann
|
||||
# steht sie rechts daneben, in der Farbe des Balkens.
|
||||
if anteil >= 10:
|
||||
achse.text(anteil / 2, y[stelle], f"{round(anteil)} %",
|
||||
ha="center", va="center", fontsize=8.5, color="white",
|
||||
fontweight="bold", zorder=4)
|
||||
else:
|
||||
achse.text(anteil + 1.5, y[stelle], f"{round(anteil)} %",
|
||||
ha="left", va="center", fontsize=8.5,
|
||||
color=FARBEN["fehler"], fontweight="bold", zorder=4)
|
||||
achse.text(101, y[stelle], name, va="center", fontsize=8.5,
|
||||
color=FARBEN["text"])
|
||||
|
||||
achse.set_yticks([])
|
||||
achse.set_xlim(0, 100)
|
||||
achse.set_xlabel("Anteil an der Gesamtzeit in Prozent")
|
||||
achse.set_title("Wo die Zeit wirklich hingeht — dieselbe Aufgabe, zwei "
|
||||
"Arten sie aufzuschreiben", fontsize=10.5,
|
||||
color=FARBEN["text"])
|
||||
achse.grid(axis="y", visible=False)
|
||||
achse.legend(fontsize=9, loc="lower center", bbox_to_anchor=(0.5, -0.28),
|
||||
ncol=2, frameon=False)
|
||||
|
||||
speichere(figur, "kap_oekosystem_laufzeit")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
pruefe_anteile()
|
||||
print("Prozentspalte stimmt mit den abgedruckten Zeiten überein.")
|
||||
zeichne()
|
||||
139
bilder_04/erzeuge_predict_vs_optimize.py
Normal file
|
|
@ -0,0 +1,139 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# erzeuge_predict_vs_optimize.py
|
||||
"""
|
||||
Erzeugt das Diagramm zum Predict-then-Optimize-Kapitel:
|
||||
|
||||
bilder_04/kap_prognose_entscheidungsguete.svg
|
||||
|
||||
Vier Verfahren, zwei Massstaebe: links der Prognosefehler (MSE), rechts die
|
||||
Kosten je Tag. Die Balken sind absichtlich nebeneinander gestellt, denn der
|
||||
Befund des Kapitels ist genau ihr Auseinanderlaufen: Das Verfahren mit dem
|
||||
BESTEN MSE hat die HOECHSTEN Kosten.
|
||||
|
||||
Ohne Formel gesagt: Der Prognostiker optimiert seinen Fehler, der Planer traegt
|
||||
die Kosten - und beide messen etwas anderes. Wer die Prognose verbessert, ohne
|
||||
zu fragen, wofuer sie gebraucht wird, verbessert die falsche Zahl.
|
||||
|
||||
Alle Werte kommen aus Predict_then_Optimize.py: dieselben Daten, dieselbe Saat,
|
||||
dieselben vier Verfahren. pruefe_befund() stellt sicher, dass die Aussage des
|
||||
Kapitels im Bild auch wirklich zu sehen ist - sonst waere das Bild ein
|
||||
Widerspruch zum Text daneben.
|
||||
|
||||
Aufruf (aus dem Repository-Wurzelverzeichnis):
|
||||
python3 bilder_04/erzeuge_predict_vs_optimize.py
|
||||
|
||||
Benoetigt: numpy, scipy, scikit-learn, matplotlib
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
HIER = os.path.dirname(os.path.abspath(__file__))
|
||||
BASIS = os.path.dirname(HIER)
|
||||
sys.path.insert(0, HIER)
|
||||
sys.path.insert(0, os.path.join(
|
||||
BASIS, "Operations_Research_mit_Python_Version_04_Programme"))
|
||||
|
||||
from stil_04 import FARBEN, speichere # noqa: E402
|
||||
import matplotlib.pyplot as plt # noqa: E402
|
||||
from scipy.stats import norm # noqa: E402
|
||||
from sklearn.linear_model import LinearRegression, QuantileRegressor # noqa: E402
|
||||
from Predict_then_Optimize import (KRITISCHES_VERHAELTNIS, TAGE, # noqa: E402
|
||||
TRAINING, erzeuge_daten, tageskosten)
|
||||
|
||||
KURZ = ["Punkt-\nprognose", "+ Zuschlag aus der\nNormalverteilung",
|
||||
"+ Zuschlag auf\nKosten trainiert", "Quantil-\nregression"]
|
||||
|
||||
|
||||
def rechne() -> list[tuple[str, float, float]]:
|
||||
"""Dieselben vier Verfahren wie im Buchprogramm."""
|
||||
merkmale, nachfrage = erzeuge_daten()[:2]
|
||||
# Dieselbe Aufteilung wie im Buchprogramm - eine selbst gewaehlte Grenze
|
||||
# ergaebe andere Zahlen als die abgedruckte Tabelle.
|
||||
lernen = slice(0, TRAINING)
|
||||
pruefen = slice(TRAINING, TAGE)
|
||||
|
||||
kq = LinearRegression().fit(merkmale[lernen], nachfrage[lernen])
|
||||
punkt = kq.predict(merkmale[pruefen])
|
||||
rest = nachfrage[lernen] - kq.predict(merkmale[lernen])
|
||||
pauschal = norm.ppf(KRITISCHES_VERHAELTNIS) * rest.std()
|
||||
|
||||
kandidaten = np.linspace(-10.0, 30.0, 401)
|
||||
trainingsprognose = kq.predict(merkmale[lernen])
|
||||
kostenzuschlag = float(kandidaten[np.argmin(
|
||||
[tageskosten(trainingsprognose + z, nachfrage[lernen])
|
||||
for z in kandidaten])])
|
||||
|
||||
quantil = QuantileRegressor(quantile=KRITISCHES_VERHAELTNIS, alpha=0.0,
|
||||
solver="highs").fit(merkmale[lernen],
|
||||
nachfrage[lernen])
|
||||
quantilprognose = quantil.predict(merkmale[pruefen])
|
||||
|
||||
verfahren = [
|
||||
(KURZ[0], punkt, punkt),
|
||||
(KURZ[1], punkt, punkt + pauschal),
|
||||
(KURZ[2], punkt, punkt + kostenzuschlag),
|
||||
(KURZ[3], quantilprognose, quantilprognose),
|
||||
]
|
||||
aus = []
|
||||
for name, prognose, bestellung in verfahren:
|
||||
mse = float(((prognose - nachfrage[pruefen]) ** 2).mean())
|
||||
aus.append((name, mse, tageskosten(bestellung, nachfrage[pruefen])))
|
||||
return aus
|
||||
|
||||
|
||||
def pruefe_befund(werte) -> None:
|
||||
"""Die Aussage des Kapitels muss im Bild sichtbar sein."""
|
||||
bester_mse = min(werte, key=lambda z: z[1])
|
||||
teuerster = max(werte, key=lambda z: z[2])
|
||||
if bester_mse[0] != teuerster[0]:
|
||||
raise SystemExit(
|
||||
f"Predict-then-Optimize-Bild: Das Verfahren mit dem besten MSE "
|
||||
f"({bester_mse[0].replace(chr(10), ' ')}) ist nicht das teuerste "
|
||||
f"({teuerster[0].replace(chr(10), ' ')}). Das Kapitel behauptet "
|
||||
f"genau das - Bild und Text wuerden sich widersprechen.")
|
||||
|
||||
|
||||
def zeichne(werte) -> None:
|
||||
figur, achsen = plt.subplots(1, 2, figsize=(9.4, 4.4))
|
||||
namen = [n for n, _, _ in werte]
|
||||
y = np.arange(len(namen))[::-1]
|
||||
|
||||
for achse, index, titel, einheit, farbe in (
|
||||
(achsen[0], 1, "Prognosefehler (MSE)", "", FARBEN["zweit"]),
|
||||
(achsen[1], 2, "Kosten je Tag", " €", FARBEN["warnung"])):
|
||||
wert = [w[index] for w in werte]
|
||||
bester = min(wert)
|
||||
farben = [FARBEN["gut"] if abs(w - bester) < 1e-9 else farbe
|
||||
for w in wert]
|
||||
achse.barh(y, wert, height=0.62, color=farben, zorder=3)
|
||||
for stelle, w in enumerate(wert):
|
||||
achse.text(w * 1.01, y[stelle], f"{w:,.2f}{einheit}"
|
||||
.replace(",", "."), va="center", fontsize=9,
|
||||
color=FARBEN["text"])
|
||||
achse.set_yticks(y)
|
||||
achse.set_yticklabels(namen if index == 1 else [], fontsize=8.5,
|
||||
color=FARBEN["text"])
|
||||
achse.set_xlim(0, max(wert) * 1.22)
|
||||
achse.set_title(titel, fontsize=10.5, color=FARBEN["text"])
|
||||
achse.grid(axis="y", visible=False)
|
||||
|
||||
figur.suptitle("Grün ist jeweils der beste Wert — und es ist nicht "
|
||||
"dasselbe Verfahren",
|
||||
fontsize=11, color=FARBEN["text"], y=1.02)
|
||||
speichere(figur, "kap_prognose_entscheidungsguete")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
werte = rechne()
|
||||
for name, mse, kosten in werte:
|
||||
print(f" {name.replace(chr(10), ' '):<40} MSE {mse:>8.2f} "
|
||||
f"Kosten {kosten:>7.2f} EUR")
|
||||
pruefe_befund(werte)
|
||||
print(" Befund bestätigt: bester MSE = höchste Kosten.")
|
||||
zeichne(werte)
|
||||
100
bilder_04/erzeuge_spaltengenerierung_flow.py
Normal file
|
|
@ -0,0 +1,100 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# erzeuge_spaltengenerierung_flow.py
|
||||
"""
|
||||
Erzeugt das Regelkreis-Diagramm zum Dekompositions-Kapitel:
|
||||
|
||||
bilder_04/kap_dekomposition_regelkreis.svg
|
||||
|
||||
Spaltengenerierung ist kein Algorithmus, den man am Stueck liest, sondern ein
|
||||
Kreislauf zwischen zwei Problemen, die einander abwechselnd Zahlen zuwerfen:
|
||||
|
||||
Master-LP -- Dualwerte pi_i --> Pricing (Rucksack)
|
||||
Master-LP <-- neue Spalte a -- Pricing
|
||||
|
||||
Das Bild macht daraus einen Kreis mit der Abbruchbedingung als Ausgang. Es ist
|
||||
bewusst schematisch: Zahlen stehen keine darin, denn die Schleife gilt fuer
|
||||
jede Instanz gleichermassen. Was sie zeigt, ist die Richtung der Information -
|
||||
und dass das Verfahren erst endet, wenn das Pricing NICHTS mehr findet.
|
||||
|
||||
Gezeichnet mit matplotlib.patches statt Graphviz: Das Bild hat fuenf Kaesten
|
||||
und vier Pfeile: dafuer eine weitere Abhaengigkeit einzufuehren, waere
|
||||
unverhaeltnismaessig - und die Positionen sollen zu den uebrigen Diagrammen des
|
||||
Buchs passen, nicht zu Graphviz' Layoutgeschmack.
|
||||
|
||||
Aufruf (aus dem Repository-Wurzelverzeichnis):
|
||||
python3 bilder_04/erzeuge_spaltengenerierung_flow.py
|
||||
|
||||
Benoetigt: matplotlib
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import os
|
||||
import sys
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
from stil_04 import FARBEN, speichere # noqa: E402
|
||||
import matplotlib.pyplot as plt # noqa: E402
|
||||
from matplotlib.patches import FancyArrowPatch, FancyBboxPatch # noqa: E402
|
||||
|
||||
|
||||
def kasten(achse, x, y, breite, hoehe, titel, zeilen, farbe, fuellung):
|
||||
achse.add_patch(FancyBboxPatch(
|
||||
(x, y), breite, hoehe,
|
||||
boxstyle="round,pad=0,rounding_size=0.09",
|
||||
facecolor=fuellung, edgecolor=farbe, linewidth=1.8, zorder=2))
|
||||
achse.text(x + breite / 2, y + hoehe - 0.20, titel, ha="center",
|
||||
va="center", fontsize=10, fontweight="bold", color=farbe,
|
||||
zorder=3)
|
||||
achse.text(x + breite / 2, y + hoehe / 2 - 0.16, zeilen, ha="center",
|
||||
va="center", fontsize=8.8, color=FARBEN["text"], zorder=3)
|
||||
|
||||
|
||||
def pfeil(achse, von, nach, text, farbe, bogen=0.0, versatz=(0, 0)):
|
||||
achse.add_patch(FancyArrowPatch(
|
||||
von, nach, arrowstyle="-|>", mutation_scale=16, linewidth=1.8,
|
||||
color=farbe, connectionstyle=f"arc3,rad={bogen}", zorder=4,
|
||||
shrinkA=2, shrinkB=2))
|
||||
mitte = ((von[0] + nach[0]) / 2 + versatz[0],
|
||||
(von[1] + nach[1]) / 2 + versatz[1])
|
||||
achse.text(*mitte, text, ha="center", va="center", fontsize=9,
|
||||
color=farbe, fontweight="bold", zorder=5)
|
||||
|
||||
|
||||
def zeichne() -> None:
|
||||
figur, achse = plt.subplots(figsize=(8.6, 4.3))
|
||||
achse.set_xlim(0, 10)
|
||||
achse.set_ylim(0, 5)
|
||||
achse.axis("off")
|
||||
|
||||
kasten(achse, 0.3, 2.5, 3.5, 1.9, "Master-LP (RMP)",
|
||||
"rechnet mit den bisher\nbekannten Mustern\nund liefert Dualwerte",
|
||||
FARBEN["haupt"], "#eef2ff")
|
||||
kasten(achse, 6.2, 2.5, 3.5, 1.9, "Pricing: Rucksackproblem",
|
||||
"sucht das EINE Muster mit\nden kleinsten reduzierten\nKosten",
|
||||
FARBEN["zweit"], "#ecfeff")
|
||||
kasten(achse, 3.1, 0.15, 3.8, 1.35, "Optimum bewiesen",
|
||||
"kein Muster mit negativen\nreduzierten Kosten mehr",
|
||||
FARBEN["gut"], "#ecfdf5")
|
||||
|
||||
pfeil(achse, (3.9, 3.9), (6.1, 3.9),
|
||||
r"Dualwerte $\pi_i$", FARBEN["haupt"], bogen=-0.28,
|
||||
versatz=(0, 0.62))
|
||||
pfeil(achse, (6.1, 3.0), (3.9, 3.0),
|
||||
r"neue Spalte $a$, wenn $1-\sum_i \pi_i a_i < 0$",
|
||||
FARBEN["zweit"], bogen=-0.28, versatz=(0, -0.72))
|
||||
pfeil(achse, (7.4, 2.4), (6.0, 1.5), "sonst", FARBEN["gut"], bogen=0.25,
|
||||
versatz=(0.55, 0.12))
|
||||
|
||||
achse.text(5.0, 4.85,
|
||||
"Die Schleife endet nicht, wenn sie müde wird, sondern wenn das "
|
||||
"Pricing nichts mehr findet.",
|
||||
ha="center", fontsize=9.5, color=FARBEN["gedaempft"],
|
||||
style="italic")
|
||||
|
||||
speichere(figur, "kap_dekomposition_regelkreis")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
zeichne()
|
||||
115
bilder_04/erzeuge_vrp_subtours.py
Normal file
|
|
@ -0,0 +1,115 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# erzeuge_vrp_subtours.py
|
||||
"""
|
||||
Erzeugt das Kurzzyklus-Diagramm zum Graphen-Kapitel:
|
||||
|
||||
bilder_04/kap_graphen_subtouren.svg
|
||||
|
||||
Links, was ein Tourenmodell OHNE Kurzzyklus-Bedingungen zurueckgeben darf:
|
||||
Jeder Kunde hat genau einen Vorgaenger und einen Nachfolger - alle
|
||||
Gradbedingungen sind erfuellt -, und trotzdem ist die Loesung unbrauchbar, weil
|
||||
die Fahrzeuge nie am Depot vorbeikommen. Rechts dieselbe Instanz mit gueltiger
|
||||
Rundreise ab und bis Depot.
|
||||
|
||||
Der Punkt: Das linke Bild ist kein Rechenfehler. Es erfuellt jede Bedingung,
|
||||
die naiv formuliert wurde. Was fehlt, ist eine Bedingung, an die man beim
|
||||
Aufschreiben nicht denkt - und genau dafuer gibt es MTZ.
|
||||
|
||||
Die Instanz stammt aus erzeuge_vrp_touren.py, also aus demselben Datensatz wie
|
||||
das Buchprogramm. Die Kurzzyklen links sind daraus konstruiert: Es sind die
|
||||
Kundenfolgen der gueltigen Loesung, nur ohne das Depot geschlossen. So zeigt
|
||||
das Bild denselben Fall und nicht ein zweites, erfundenes Beispiel.
|
||||
|
||||
Aufruf (aus dem Repository-Wurzelverzeichnis):
|
||||
python3 bilder_04/erzeuge_vrp_subtours.py
|
||||
|
||||
Benoetigt: numpy, matplotlib, ortools
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
from stil_04 import FARBEN, speichere # noqa: E402
|
||||
import matplotlib.pyplot as plt # noqa: E402
|
||||
from erzeuge_vrp_touren import (TOURFARBEN, erzeuge_daten, # noqa: E402
|
||||
loese)
|
||||
|
||||
|
||||
def kurzzyklen(touren) -> list[list[int]]:
|
||||
"""Aus jeder gueltigen Tour einen Kreis ohne Depot machen.
|
||||
|
||||
Genau das darf ein Modell liefern, das nur 'jeder Kunde einmal angefahren'
|
||||
fordert: In sich geschlossene Kreise, die das Depot nie beruehren."""
|
||||
aus = []
|
||||
for tour in touren:
|
||||
kunden = [k for k in tour["stationen"] if k != 0]
|
||||
if len(kunden) >= 2:
|
||||
aus.append(kunden + [kunden[0]])
|
||||
return aus
|
||||
|
||||
|
||||
def pruefe_grade(zyklen, anzahl_kunden: int) -> None:
|
||||
"""Die unzulaessige Loesung muss die Gradbedingungen trotzdem erfuellen."""
|
||||
besucht = [k for z in zyklen for k in z[:-1]]
|
||||
if sorted(besucht) != list(range(1, anzahl_kunden + 1)):
|
||||
raise SystemExit(
|
||||
"Subtour-Bild: Die konstruierten Kurzzyklen besuchen nicht jeden "
|
||||
"Kunden genau einmal. Dann waere das linke Bild kein gueltiges "
|
||||
"Gegenbeispiel, sondern einfach falsch.")
|
||||
|
||||
|
||||
def zeichne(koordinaten, touren, zyklen) -> None:
|
||||
figur, achsen = plt.subplots(1, 2, figsize=(9.6, 4.6), sharex=True,
|
||||
sharey=True)
|
||||
|
||||
for achse, titel in ((achsen[0], "ohne Kurzzyklus-Bedingungen"),
|
||||
(achsen[1], "mit MTZ — jede Tour über das Depot")):
|
||||
for ort in range(1, len(koordinaten)):
|
||||
achse.scatter(*koordinaten[ort], s=34, color=FARBEN["text"],
|
||||
zorder=3)
|
||||
achse.scatter(*koordinaten[0], marker="s", s=150,
|
||||
color=FARBEN["fehler"], zorder=4)
|
||||
achse.annotate("Depot", koordinaten[0], textcoords="offset points",
|
||||
xytext=(9, 7), fontsize=9, color=FARBEN["fehler"],
|
||||
fontweight="bold")
|
||||
achse.set_xlabel("km")
|
||||
achse.set_title(titel, fontsize=10.5, color=FARBEN["text"])
|
||||
|
||||
# links: die Kurzzyklen
|
||||
for nummer, (zyklus, farbe) in enumerate(zip(zyklen, TOURFARBEN), start=1):
|
||||
weg = koordinaten[zyklus]
|
||||
achsen[0].plot(weg[:, 0], weg[:, 1], "-", color=farbe, linewidth=1.8,
|
||||
alpha=0.9, zorder=2)
|
||||
achsen[0].annotate(
|
||||
"Jeder Kunde hat genau einen Vor- und\n"
|
||||
"einen Nachfolger — und kein Fahrzeug\n"
|
||||
"kommt je am Depot vorbei.",
|
||||
xy=(0.03, 0.03), xycoords="axes fraction", fontsize=9,
|
||||
color=FARBEN["fehler"], va="bottom")
|
||||
|
||||
# rechts: die gueltige Loesung
|
||||
for tour, farbe in zip(touren, TOURFARBEN):
|
||||
weg = koordinaten[tour["stationen"]]
|
||||
achsen[1].plot(weg[:, 0], weg[:, 1], "-", color=farbe, linewidth=1.8,
|
||||
alpha=0.9, zorder=2)
|
||||
|
||||
achsen[0].set_ylabel("km")
|
||||
figur.suptitle("Warum Gradbedingungen allein nicht genügen",
|
||||
fontsize=11, color=FARBEN["text"], y=1.01)
|
||||
speichere(figur, "kap_graphen_subtouren")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
koordinaten, distanz = erzeuge_daten()
|
||||
touren = loese(distanz)
|
||||
zyklen = kurzzyklen(touren)
|
||||
pruefe_grade(zyklen, len(koordinaten) - 1)
|
||||
print(f" {len(zyklen)} Kurzzyklen konstruiert, "
|
||||
f"{sum(len(z) - 1 for z in zyklen)} Kunden — jeder genau einmal.")
|
||||
zeichne(koordinaten, touren, zyklen)
|
||||
115
bilder_04/erzeuge_zielgroessen_radar.py
Normal file
|
|
@ -0,0 +1,115 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# erzeuge_zielgroessen_radar.py
|
||||
"""
|
||||
Erzeugt das Netzdiagramm zum Kapitel "Vom Management-Wunsch zum Modell":
|
||||
|
||||
bilder_04/kap_modellierung_zielgroessen.svg
|
||||
|
||||
Drei Saetze aus derselben Besprechung - "Umsatz maximieren", "Maschinen
|
||||
auslasten", "Deckungsbeitrag maximieren" - ergeben drei verschiedene Plaene aus
|
||||
denselben 14 Anfragen. Das Netzdiagramm legt sie uebereinander. Der Punkt, auf
|
||||
den es ankommt: Der Auslastungsplan sieht auf seiner eigenen Achse perfekt aus
|
||||
und faellt auf der Deckungsbeitragsachse durch.
|
||||
|
||||
KEINE ZAHL IST HIER EINGETRAGEN. Das Skript importiert Vom_Wunsch_zum_Modell.py
|
||||
und loest dieselben Modelle mit demselben Solver. Weicht das Buchprogramm
|
||||
spaeter ab, aendert sich das Bild mit - nicht der Text neben einem alten Bild.
|
||||
|
||||
Aufruf (aus dem Repository-Wurzelverzeichnis):
|
||||
python3 bilder_04/erzeuge_zielgroessen_radar.py
|
||||
|
||||
Benoetigt: numpy, scipy, matplotlib
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
HIER = os.path.dirname(os.path.abspath(__file__))
|
||||
BASIS = os.path.dirname(HIER)
|
||||
sys.path.insert(0, HIER)
|
||||
# Das Buchprogramm ist die Quelle der Instanz - importiert, nicht kopiert.
|
||||
sys.path.insert(0, os.path.join(
|
||||
BASIS, "Operations_Research_mit_Python_Version_04_Programme"))
|
||||
|
||||
from stil_04 import FARBEN, speichere # noqa: E402
|
||||
import matplotlib.pyplot as plt # noqa: E402
|
||||
from Vom_Wunsch_zum_Modell import (AUFTRAEGE, DECKUNGSBEITRAG, # noqa: E402
|
||||
STUNDEN, UMSATZ, kennzahlen,
|
||||
waehle_aus)
|
||||
|
||||
# Die fuenf Achsen des Netzes und wie sie aus kennzahlen() entstehen. Jede
|
||||
# Achse wird auf ihr eigenes Maximum ueber die drei Plaene skaliert - sonst
|
||||
# erschluege der Umsatz in Euro alle uebrigen Groessen.
|
||||
ACHSEN = [
|
||||
("Deckungsbeitrag", lambda k: k["db"]),
|
||||
("Umsatz", lambda k: k["umsatz"]),
|
||||
("Maschinen-\nauslastung", lambda k: k["auslastung"]),
|
||||
("Stammkunden\nbedient", lambda k: len({a[1] for a in AUFTRAEGE if a[2]})
|
||||
- k["kunden_leer"]),
|
||||
("angenommene\nAufträge", lambda k: k["angenommen"]),
|
||||
]
|
||||
|
||||
PLAENE = [
|
||||
("„Umsatz maximieren“", UMSATZ, FARBEN["zweit"]),
|
||||
("„Maschinen auslasten“", STUNDEN, FARBEN["warnung"]),
|
||||
("„Deckungsbeitrag maximieren“", DECKUNGSBEITRAG, FARBEN["gut"]),
|
||||
]
|
||||
|
||||
|
||||
def loese_alle() -> list[tuple[str, dict, str]]:
|
||||
ergebnis = []
|
||||
for name, ziel, farbe in PLAENE:
|
||||
plan = waehle_aus(ziel)
|
||||
if plan is None:
|
||||
raise SystemExit(f"Radar: Fuer {name} liefert das Modell keinen Plan.")
|
||||
ergebnis.append((name, kennzahlen(plan), farbe))
|
||||
return ergebnis
|
||||
|
||||
|
||||
def zeichne(plaene) -> None:
|
||||
werte = {name: [f(k) for _, f in ACHSEN] for name, k, _ in plaene}
|
||||
maxima = [max(werte[name][i] for name in werte) or 1.0
|
||||
for i in range(len(ACHSEN))]
|
||||
|
||||
winkel = np.linspace(0, 2 * np.pi, len(ACHSEN), endpoint=False)
|
||||
geschlossen = np.concatenate([winkel, winkel[:1]])
|
||||
|
||||
figur, achse = plt.subplots(figsize=(6.6, 5.8),
|
||||
subplot_kw={"projection": "polar"})
|
||||
for name, _k, farbe in plaene:
|
||||
anteile = [w / m * 100 for w, m in zip(werte[name], maxima)]
|
||||
anteile += anteile[:1]
|
||||
achse.plot(geschlossen, anteile, "-", color=farbe, linewidth=2.0,
|
||||
label=name, zorder=3)
|
||||
achse.fill(geschlossen, anteile, color=farbe, alpha=0.12, zorder=2)
|
||||
|
||||
achse.set_xticks(winkel)
|
||||
achse.set_xticklabels([bez for bez, _ in ACHSEN], fontsize=9,
|
||||
color=FARBEN["text"])
|
||||
achse.set_yticks([25, 50, 75, 100])
|
||||
achse.set_yticklabels(["25 %", "50 %", "75 %", "100 %"], fontsize=8,
|
||||
color=FARBEN["gedaempft"])
|
||||
achse.set_ylim(0, 108)
|
||||
achse.grid(color=FARBEN["linie"], linewidth=0.7, alpha=0.6)
|
||||
achse.spines["polar"].set_color(FARBEN["linie"])
|
||||
achse.set_title("Drei Sätze aus derselben Besprechung, drei Pläne\n"
|
||||
"je Achse auf den besten der drei Pläne bezogen",
|
||||
fontsize=10.5, color=FARBEN["text"], pad=24)
|
||||
achse.legend(fontsize=9, loc="lower center", bbox_to_anchor=(0.5, -0.20),
|
||||
frameon=False)
|
||||
|
||||
speichere(figur, "kap_modellierung_zielgroessen")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
plaene = loese_alle()
|
||||
for name, k, _ in plaene:
|
||||
print(f" {name:<32} Umsatz {k['umsatz']:>8,.0f} "
|
||||
f"DB {k['db']:>8,.0f} Auslastung {k['auslastung']:>3.0f} % "
|
||||
f"angenommen {k['angenommen']:>2}")
|
||||
zeichne(plaene)
|
||||
2284
bilder_04/kap_dekomposition_regelkreis.svg
Normal file
|
After Width: | Height: | Size: 66 KiB |
1909
bilder_04/kap_dp_almgren_chriss.svg
Normal file
|
After Width: | Height: | Size: 60 KiB |
2440
bilder_04/kap_einfuehrung_kombinatorik.svg
Normal file
|
After Width: | Height: | Size: 74 KiB |
2026
bilder_04/kap_graphen_subtouren.svg
Normal file
|
After Width: | Height: | Size: 62 KiB |
3309
bilder_04/kap_lp_entartung.svg
Normal file
|
After Width: | Height: | Size: 95 KiB |
3483
bilder_04/kap_milp_big_m.svg
Normal file
|
After Width: | Height: | Size: 106 KiB |
1556
bilder_04/kap_modellierung_zielgroessen.svg
Normal file
|
After Width: | Height: | Size: 53 KiB |
1820
bilder_04/kap_oekosystem_laufzeit.svg
Normal file
|
After Width: | Height: | Size: 55 KiB |
1888
bilder_04/kap_prognose_entscheidungsguete.svg
Normal file
|
After Width: | Height: | Size: 55 KiB |
1972
bilder_04/kap_supplychain_lastabwurf.svg
Normal file
|
After Width: | Height: | Size: 66 KiB |
2542
bilder_04/kap_unsicherheit_durchschnitt.svg
Normal file
|
After Width: | Height: | Size: 74 KiB |
92
bilder_04/stil_04.py
Normal file
|
|
@ -0,0 +1,92 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# stil_04.py
|
||||
"""
|
||||
Gemeinsamer Unterbau der Diagramm-Generatoren: Farben, Grundeinstellungen und
|
||||
das Schreiben der SVG-Datei.
|
||||
|
||||
Warum es das gibt: Elf neue Generatoren trugen sonst elfmal denselben Vorspann.
|
||||
Eine Zeile daraus ist besonders heikel -
|
||||
|
||||
plt.rcParams["svg.hashsalt"] = "or-mit-python-v04"
|
||||
|
||||
Ohne sie wuerfelt matplotlib die clip-path-IDs je Figure neu; die Datei aendert
|
||||
sich dann bei JEDEM Lauf, ohne dass sich am Bild etwas geaendert haette. Das
|
||||
faellt nicht sofort auf, sondern erst als unlesbarer git-diff - und dann sucht
|
||||
man die Ursache lange. Dasselbe gilt fuer metadata={"Date": None} beim
|
||||
Speichern: Sonst steht der Schreibzeitstempel auf die Mikrosekunde im SVG.
|
||||
|
||||
Die 33 aelteren Generatoren setzen beides jeweils selbst. Sie bleiben
|
||||
unveraendert - eine Umstellung waere ein eigener Schritt und braeuchte den
|
||||
Nachweis, dass jedes einzelne Bild byte-identisch bleibt.
|
||||
|
||||
Benutzung:
|
||||
|
||||
from stil_04 import FARBEN, neue_figur, speichere
|
||||
|
||||
figur, achse = neue_figur(breite=7.6, hoehe=4.6)
|
||||
...
|
||||
speichere(figur, "kap_einfuehrung_kombinatorik")
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import os
|
||||
|
||||
import matplotlib
|
||||
matplotlib.use("Agg")
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
# Reproduzierbare SVG-Ausgabe - siehe Modulkommentar.
|
||||
plt.rcParams["svg.hashsalt"] = "or-mit-python-v04"
|
||||
|
||||
# Serifenlos und ueberall vorhanden. Bewusst kein Stylesheet von seaborn: Das
|
||||
# Buch hat eine eigene Palette, und zwei Stilsysteme nebeneinander sehen nach
|
||||
# zwei Buechern aus.
|
||||
plt.rcParams["font.family"] = "DejaVu Sans"
|
||||
plt.rcParams["axes.grid"] = True
|
||||
plt.rcParams["grid.linestyle"] = ":"
|
||||
plt.rcParams["grid.alpha"] = 0.4
|
||||
plt.rcParams["axes.spines.top"] = False
|
||||
plt.rcParams["axes.spines.right"] = False
|
||||
|
||||
BASIS = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
|
||||
BILDER = os.path.join(BASIS, "bilder_04")
|
||||
|
||||
# Die Palette des Buchs - dieselben Werte wie in den 33 vorhandenen Diagrammen
|
||||
# und auf der Titelseite. Die Namen sagen, wofuer die Farbe steht, nicht wie
|
||||
# sie aussieht: So bleibt die Bedeutung erhalten, falls ein Ton je gewechselt
|
||||
# wird.
|
||||
FARBEN = {
|
||||
"haupt": "#4338ca", # Indigo - die primaere Groesse, Zielfunktion
|
||||
"zweit": "#0891b2", # Teal - die Vergleichsgroesse
|
||||
"warnung": "#b45309", # Amber - teuer, aber zulaessig
|
||||
"gut": "#059669", # Emerald - Optimum, zulaessige Menge
|
||||
"fehler": "#be123c", # Rose - unzulaessig, Verlust, Bruch
|
||||
"akzent": "#7c3aed", # Violett - Schranken, Niveaulinien
|
||||
"text": "#334155", # Slate - Beschriftung
|
||||
"gedaempft": "#64748b", # Slate hell - Nebeninformation
|
||||
"linie": "#94a3b8", # Hilfslinien
|
||||
"flaeche": "#f1f5f9", # Hinterlegte Flaechen
|
||||
}
|
||||
|
||||
|
||||
def neue_figur(breite: float = 7.6, hoehe: float = 4.6, **kwargs):
|
||||
"""Figure und Achse in der Standardgroesse des Buchs."""
|
||||
return plt.subplots(figsize=(breite, hoehe), **kwargs)
|
||||
|
||||
|
||||
def speichere(figur, name: str) -> str:
|
||||
"""Schreibt bilder_04/<name>.svg und meldet den Pfad.
|
||||
|
||||
Nur SVG: Das PDF bindet die Bilder ueber \\includesvg (inkscape) ein, und
|
||||
PNG-Zweitfassungen gab es in diesem Projekt schon einmal - sie wurden von
|
||||
nirgendwo referenziert und mussten wieder geloescht werden."""
|
||||
figur.tight_layout()
|
||||
os.makedirs(BILDER, exist_ok=True)
|
||||
pfad = os.path.join(BILDER, f"{name}.svg")
|
||||
figur.savefig(pfad, format="svg", dpi=160, bbox_inches="tight",
|
||||
metadata={"Date": None})
|
||||
plt.close(figur)
|
||||
print(f"geschrieben: {pfad}")
|
||||
return pfad
|
||||