Version 04 als eigenes Repository
Erster Commit des Strangs "Optimierte Entscheidungsfindung mit Python" (Version 04). Die Historie der 71 Commits bis zur Trennung bleibt im uebergeordneten Repository OR_mit_Python liegen, das ab jetzt nur noch Version_03 (eingefroren) verwaltet und Version_04/ ignoriert. Bewusst kein "git subtree split": Der Pfad Version_04/ existiert erst seit der Verzeichnistrennung, ein Split braechte daher nur 7 der 41 einschlaegigen Commits - eine Teilhistorie, die vollstaendig aussieht und es nicht ist. Stand: 5 Teile, 23 Kapitel, 5 Anhaenge, 292 Abschnitte, 703 Querverweise, 325 Indexmarken, 73 Beispielprogramme, 32 SVGs, 4 Plotly-Figuren, 25 Notebooks, PDF mit 715 Seiten. Zusaetzlich in diesem Commit: * pyproject.toml mit Abhaengigkeitsgruppen finance, large-scale, api, figures, dev, empfehlungen. Die abgedruckte requirements.txt bleibt unveraendert daneben bestehen. ortools steht in der Grundausstattung, highspy erst in [large-scale] - so kann der HiGHS-Symbolkonflikt bei der schlanken Installation gar nicht erst auftreten. * Dabei zwei Funde: graphviz wird von erzeuge_architektur_diagramme.py importiert, fehlt aber in requirements.txt (jetzt in [figures]); pymoo steht in requirements.txt, wird aber von keinem Programm importiert, sondern nur im Kapitel Metaheuristiken empfohlen (jetzt in [empfehlungen]). * NEUER_TITEL.md nach Kritik_und_Verbesserungsvorschlaege/ verschoben - es ist die Vorlage des Titelblatts, kein Bestandteil des Werks. Die beiden Fundstellen in PROGRESS.md und erzeuge_titelseite.py nachgezogen. * PROGRESS.md nannte noch den Untertitel der ersten Fassung; auf den tatsaechlichen aus erzeuge_titelseite.py korrigiert. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This commit is contained in:
dschlueter
commit
b7af2f1d9a
73 changed files with 262141 additions and 0 deletions
96
OR_HTML_04/programme/Backtest_Fallen.py
Normal file
96
OR_HTML_04/programme/Backtest_Fallen.py
Normal file
|
|
@ -0,0 +1,96 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# Backtest_Fallen.py
|
||||
"""
|
||||
Kapitel Handelsmaschine: Wie leicht man sich selbst betruegt.
|
||||
|
||||
Demonstriert an SIMULIERTEN Daten (also ohne jede echte Prognosekraft),
|
||||
wie gross die scheinbare Ueberrendite durch typische Backtest-Fehler wird.
|
||||
Weil die Daten reines Rauschen sind, ist JEDE positive Ueberrendite ein Artefakt.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
|
||||
def erzeuge_zufallsmarkt(tage=1260, titel=15, seed=7):
|
||||
"""Reines Rauschen: kein Titel hat echte Prognosekraft."""
|
||||
rng = np.random.default_rng(seed)
|
||||
renditen = rng.normal(0.0003, 0.015, size=(tage, titel))
|
||||
index = pd.bdate_range("2020-01-01", periods=tage)
|
||||
return pd.DataFrame(renditen, index=index,
|
||||
columns=[f"T{i:02d}" for i in range(titel)])
|
||||
|
||||
|
||||
def strategie(renditen, lookahead=False, survivorship=False,
|
||||
ohne_kosten=False, gebuehr=0.0015):
|
||||
"""Einfache Momentum-Strategie mit optional eingebauten Fehlern."""
|
||||
fenster = 60
|
||||
kapital = [100.0]
|
||||
gewichte = np.ones(renditen.shape[1]) / renditen.shape[1]
|
||||
|
||||
daten = renditen
|
||||
if survivorship:
|
||||
# FEHLER: nur die im Nachhinein besten Titel ins Universum lassen
|
||||
beste = renditen.sum().nlargest(renditen.shape[1] // 2).index
|
||||
daten = renditen[beste]
|
||||
gewichte = np.ones(len(beste)) / len(beste)
|
||||
|
||||
for t in range(fenster, len(daten) - 1):
|
||||
if t % 21 == 0: # monatlich
|
||||
if lookahead:
|
||||
# FEHLER: Signal nutzt die Rendite von MORGEN
|
||||
signal = daten.iloc[t + 1].values
|
||||
else:
|
||||
signal = daten.iloc[t - fenster:t].mean().values
|
||||
|
||||
neue = (signal > np.median(signal)).astype(float)
|
||||
neue = neue / neue.sum() if neue.sum() > 0 else gewichte
|
||||
if not ohne_kosten:
|
||||
kapital[-1] *= (1 - np.abs(neue - gewichte).sum() * gebuehr)
|
||||
gewichte = neue
|
||||
|
||||
kapital.append(kapital[-1] * (1 + float(gewichte @ daten.iloc[t + 1].values)))
|
||||
return np.array(kapital)
|
||||
|
||||
|
||||
def cagr(verlauf, jahre):
|
||||
return (verlauf[-1] / verlauf[0]) ** (1 / jahre) - 1
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
renditen = erzeuge_zufallsmarkt()
|
||||
jahre = len(renditen) / 252
|
||||
markt = 100 * (1 + renditen.mean(axis=1)).cumprod().values
|
||||
|
||||
print("=" * 84)
|
||||
print(" DIE FUENF SELBSTTAEUSCHUNGEN - GEMESSEN AN REINEM RAUSCHEN")
|
||||
print("=" * 84)
|
||||
print("Die Daten sind zufaellig erzeugt. Es gibt KEINE echte Prognosekraft.")
|
||||
print("Jede Ueberrendite unten ist daher ein reines Artefakt.\n")
|
||||
print(f"{'Variante':<44} {'CAGR':>9} {'ggue. Markt':>13}")
|
||||
print("-" * 84)
|
||||
|
||||
basis = cagr(markt, jahre)
|
||||
print(f"{'Markt (Gleichgewichtung, Referenz)':<44} {basis*100:>8.2f} % "
|
||||
f"{0.0:>12.2f} %")
|
||||
|
||||
varianten = [
|
||||
("Ehrlich: Signal aus Vergangenheit, mit Kosten", {}),
|
||||
("FEHLER 1: Lookahead (Signal kennt morgen)", {"lookahead": True}),
|
||||
("FEHLER 2: Survivorship (nur Gewinner im Universum)", {"survivorship": True}),
|
||||
("FEHLER 3: Transaktionskosten ignoriert", {"ohne_kosten": True}),
|
||||
("FEHLER 1+2+3 kombiniert", {"lookahead": True, "survivorship": True,
|
||||
"ohne_kosten": True}),
|
||||
]
|
||||
for name, argumente in varianten:
|
||||
verlauf = strategie(renditen, **argumente)
|
||||
wert = cagr(verlauf, jahre)
|
||||
print(f"{name:<44} {wert*100:>8.2f} % {(wert-basis)*100:>+12.2f} %")
|
||||
|
||||
print("-" * 84)
|
||||
print("FEHLER 4 (Overfitting) und FEHLER 5 (Data Snooping) lassen sich so nicht")
|
||||
print("zeigen - sie entstehen erst durch WIEDERHOLTES Probieren. Faustregel:")
|
||||
print("Wer 20 Strategien testet, findet auch in reinem Rauschen eine mit")
|
||||
print("Signifikanz auf dem 5-%-Niveau. Genau das ist der Punkt.")
|
||||
print("=" * 84)
|
||||
83
OR_HTML_04/programme/Bausteine_Vorlage.py
Normal file
83
OR_HTML_04/programme/Bausteine_Vorlage.py
Normal file
|
|
@ -0,0 +1,83 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# Bausteine_Vorlage.py
|
||||
"""
|
||||
Kapitel Einfuehrung: Die Bausteine-Vorlage als ausfuellbares Python-Geruest.
|
||||
|
||||
Uebersetzt die Vier-Bausteine-Vorlage aus dem Abschnitt 'Die vier
|
||||
universellen Bausteine jedes OR-Problems' direkt in Code: Wer sie
|
||||
auf Papier ausgefuellt hat, kann sie fast unveraendert in ein loesbares
|
||||
Modell umsetzen. Angewendet auf das Baeckerei-Beispiel der Handrechnung.
|
||||
"""
|
||||
|
||||
from dataclasses import dataclass
|
||||
from ortools.sat.python import cp_model
|
||||
|
||||
|
||||
@dataclass
|
||||
class Bausteine:
|
||||
"""Die vier Bausteine jedes OR-Problems - wortwoertlich aus der Vorlage
|
||||
im Abschnitt 'Die vier universellen Bausteine' uebernommen."""
|
||||
variablen: str
|
||||
parameter: str
|
||||
zielfunktion: str
|
||||
harte_bedingungen: str
|
||||
weiche_bedingungen: str
|
||||
|
||||
def zeige(self) -> None:
|
||||
print("=" * 78)
|
||||
print(" AUSGEFUELLTE BAUSTEINE-VORLAGE")
|
||||
print("=" * 78)
|
||||
for feld, wert in [
|
||||
("Entscheidungsvariablen", self.variablen),
|
||||
("Parameter", self.parameter),
|
||||
("Zielfunktion", self.zielfunktion),
|
||||
("Harte Bedingungen", self.harte_bedingungen),
|
||||
("Weiche Bedingungen", self.weiche_bedingungen),
|
||||
]:
|
||||
print(f"{feld:<24}: {wert}")
|
||||
print("=" * 78)
|
||||
|
||||
|
||||
# --- Baeckerei-Beispiel aus der Handrechnung ---------------------------------
|
||||
VORLAGE = Bausteine(
|
||||
variablen="x1 = Anzahl Brote, x2 = Anzahl Broetchen-Zehnerpackungen "
|
||||
"(beide ganzzahlig >= 0)",
|
||||
parameter="Deckungsbeitraege (2,50; 3,00) EUR; Mehlbedarf (0,5; 0,6) kg; "
|
||||
"Ofenzeit (4; 3) min; Kapazitaeten (90 kg, 600 min); Mindestmenge 40 Brote",
|
||||
zielfunktion="max 2,5*x1 + 3,0*x2 (Einheit: Euro Deckungsbeitrag)",
|
||||
harte_bedingungen="0,5*x1 + 0,6*x2 <= 90 (Mehl); 4*x1 + 3*x2 <= 600 (Ofen); "
|
||||
"x1 >= 40 (Vertrag)",
|
||||
weiche_bedingungen="keine",
|
||||
)
|
||||
|
||||
|
||||
def loese_baeckerei():
|
||||
"""Baut aus der ausgefuellten Vorlage mechanisch ein CP-SAT-Modell."""
|
||||
modell = cp_model.CpModel()
|
||||
x1 = modell.NewIntVar(40, 1000, "Brote") # Vertrag: mindestens 40
|
||||
x2 = modell.NewIntVar(0, 1000, "Broetchen_Zehner")
|
||||
|
||||
modell.Add(5 * x1 + 6 * x2 <= 900) # Mehl, x10 fuer Ganzzahligkeit
|
||||
modell.Add(4 * x1 + 3 * x2 <= 600) # Ofenzeit
|
||||
|
||||
modell.Maximize(25 * x1 + 30 * x2) # Deckungsbeitrag x10
|
||||
|
||||
loeser = cp_model.CpSolver()
|
||||
status = loeser.Solve(modell)
|
||||
if status not in (cp_model.OPTIMAL, cp_model.FEASIBLE):
|
||||
raise SystemExit(f"Keine Loesung: {loeser.StatusName(status)}")
|
||||
return loeser.Value(x1), loeser.Value(x2), loeser.ObjectiveValue() / 10
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
VORLAGE.zeige()
|
||||
|
||||
print("\n--- Aus der Vorlage mechanisch abgeleitetes Modell ---")
|
||||
brote, broetchen, gewinn = loese_baeckerei()
|
||||
print(f"Optimale Brote: {brote}")
|
||||
print(f"Optimale Broetchen-Zehner: {broetchen}")
|
||||
print(f"Maximaler Deckungsbeitrag: {gewinn:.2f} EUR")
|
||||
print()
|
||||
print("Genau diese drei Codezeilen (Variablen, Add, Maximize) entstehen")
|
||||
print("direkt aus den fuenf Zeilen der Vorlage - das ist der ganze Trick.")
|
||||
85
OR_HTML_04/programme/Bellman_Minimalbeispiel.py
Normal file
85
OR_HTML_04/programme/Bellman_Minimalbeispiel.py
Normal file
|
|
@ -0,0 +1,85 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# Bellman_Minimalbeispiel.py
|
||||
"""
|
||||
Kapitel Dynamische Programmierung: Die Handrechnung zur Rueckwaertsinduktion als Code.
|
||||
Zeigt die Wertfunktionstabelle und die optimale Politik Schritt fuer Schritt.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
|
||||
GESAMT = 3 # zu verkaufende Einheiten
|
||||
PERIODEN = 2 # Anzahl Verkaufsperioden
|
||||
|
||||
|
||||
def kosten(menge: int) -> float:
|
||||
"""Ueberproportionale Marktauswirkung: doppelte Menge kostet vierfach."""
|
||||
return float(menge ** 2)
|
||||
|
||||
|
||||
def loese_rueckwaerts():
|
||||
# V[t, x] = minimale Restkosten, wenn zu Beginn von Periode t noch x Stueck offen sind
|
||||
V = np.full((PERIODEN + 1, GESAMT + 1), np.inf)
|
||||
politik = np.zeros((PERIODEN, GESAMT + 1), dtype=int)
|
||||
|
||||
# Endbedingung: nach der letzten Periode darf nichts mehr offen sein
|
||||
V[PERIODEN, 0] = 0.0
|
||||
|
||||
print("=" * 70)
|
||||
print(" RUECKWAERTSINDUKTION SCHRITT FUER SCHRITT")
|
||||
print("=" * 70)
|
||||
|
||||
for t in range(PERIODEN - 1, -1, -1):
|
||||
letzte_periode = (t == PERIODEN - 1)
|
||||
print(f"\nStufe t = {t}" + (" (letzte Periode: alles muss weg)" if letzte_periode
|
||||
else " (freie Wahl der Menge)"))
|
||||
print(f" {'Zustand x':>10} | {'beste Aktion':>12} | {'Sofortkosten':>13} | "
|
||||
f"{'V[t+1]':>9} | {'V[t]':>8}")
|
||||
print(" " + "-" * 62)
|
||||
|
||||
for x in range(GESAMT + 1):
|
||||
aktionen = [x] if letzte_periode else range(x + 1)
|
||||
bester_wert, beste_aktion, beste_teile = np.inf, 0, (0.0, 0.0)
|
||||
|
||||
for n in aktionen:
|
||||
rest = x - n
|
||||
sofort = kosten(n)
|
||||
zukunft = V[t + 1, rest]
|
||||
gesamt = sofort + zukunft
|
||||
if gesamt < bester_wert:
|
||||
bester_wert, beste_aktion = gesamt, n
|
||||
beste_teile = (sofort, zukunft)
|
||||
|
||||
V[t, x] = bester_wert
|
||||
politik[t, x] = beste_aktion
|
||||
print(f" {x:>10} | {beste_aktion:>12} | {beste_teile[0]:>13.1f} | "
|
||||
f"{beste_teile[1]:>9.1f} | {bester_wert:>8.1f}")
|
||||
|
||||
return V, politik
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
V, politik = loese_rueckwaerts()
|
||||
|
||||
# --- Vorwaertspfad: der optimalen Politik folgen ---------------------
|
||||
print("\n" + "=" * 70)
|
||||
print(" OPTIMALER PFAD (Vorwaertssimulation)")
|
||||
print("=" * 70)
|
||||
bestand = GESAMT
|
||||
gesamtkosten = 0.0
|
||||
for t in range(PERIODEN):
|
||||
aktion = politik[t, bestand]
|
||||
gesamtkosten += kosten(aktion)
|
||||
print(f" Periode {t}: Bestand {bestand} -> verkaufe {aktion} "
|
||||
f"(Kosten {kosten(aktion):.1f}) -> Rest {bestand - aktion}")
|
||||
bestand -= aktion
|
||||
|
||||
print(f"\n Gesamtkosten: {gesamtkosten:.1f} (V[0, {GESAMT}] = {V[0, GESAMT]:.1f})")
|
||||
assert abs(gesamtkosten - V[0, GESAMT]) < 1e-9, "Pfadkosten != Wertfunktion!"
|
||||
|
||||
# --- Vergleich mit naiven Strategien ---------------------------------
|
||||
alles_sofort = kosten(GESAMT)
|
||||
print(f"\n Zum Vergleich - alles in Periode 0 verkaufen: {alles_sofort:.1f}")
|
||||
print(f" Ersparnis durch Stueckelung: {alles_sofort - gesamtkosten:.1f} "
|
||||
f"({(1 - gesamtkosten/alles_sofort)*100:.0f} %)")
|
||||
print("=" * 70)
|
||||
202
OR_HTML_04/programme/Benchmark_Skalierung.py
Normal file
202
OR_HTML_04/programme/Benchmark_Skalierung.py
Normal file
|
|
@ -0,0 +1,202 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# Benchmark_Skalierung.py
|
||||
"""
|
||||
Kapitel Testen: Ein Vergleich, dem man glauben kann.
|
||||
|
||||
Solververgleiche stehen in jedem Blog, und die meisten sind wertlos - nicht
|
||||
weil falsch gemessen wurde, sondern weil zu wenig dazugesagt wird. Dieses
|
||||
Programm misst dasselbe Transportproblem in drei Groessen mit vier
|
||||
Bibliotheken und haelt sich dabei an fuenf Regeln, die den Unterschied machen:
|
||||
|
||||
1. EIGENER PROZESS je Bibliothek. Nicht nur wegen des Importkonflikts
|
||||
zwischen ortools und highspy (Kapitel Oekosystem) - auch, damit der
|
||||
Speicherverbrauch der einen nicht in der Messung der anderen auftaucht.
|
||||
2. AUFBAU UND LOESEN GETRENNT messen. Bei grossen Instanzen ist der Aufbau
|
||||
des Modells in Python regelmaessig teurer als das Loesen. Wer nur die
|
||||
Gesamtzeit misst, optimiert am Ende die falsche Haelfte.
|
||||
3. ZIELWERTE GEGENEINANDER PRUEFEN. Eine Bibliothek, die schneller ist und
|
||||
etwas anderes ausrechnet, hat den Vergleich nicht gewonnen. Diese Pruefung
|
||||
ist der wichtigste Teil des Programms.
|
||||
4. SPEICHER MITMESSEN. Bei 10.000 Variablen entscheidet oft er und nicht die
|
||||
Zeit darueber, was auf einer Maschine noch laeuft.
|
||||
5. DIESELBE INSTANZ fuer alle. Feste Saat, kein Neuwuerfeln zwischendurch.
|
||||
|
||||
Und die Einschraenkung, die dazugehoert: Gemessen wird EIN Problemtyp in EINER
|
||||
Formulierung auf EINER Maschine. Das Ergebnis ist keine Rangliste der Solver,
|
||||
sondern eine Entscheidungshilfe fuer genau diesen Fall. Wer es verallgemeinert,
|
||||
macht denselben Fehler wie jemand, der aus einem Backtest auf die Zukunft
|
||||
schliesst (Kapitel Handelsmaschine).
|
||||
|
||||
Benoetigt: numpy; in den Kindprozessen scipy, highspy, ortools, cvxpy
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import json
|
||||
import subprocess
|
||||
import sys
|
||||
import textwrap
|
||||
|
||||
GROESSEN = [(10, 10), (32, 32), (100, 100)] # (Lager, Kunden) -> 100 / 1.024 / 10.000 Variablen
|
||||
|
||||
|
||||
# Jeder Eintrag ist ein eigenstaendiges Programm: Instanz aufbauen, loesen,
|
||||
# Ergebnis als JSON ausgeben. Die Instanz wird in jedem Kindprozess aus
|
||||
# derselben Saat neu erzeugt - so reist nichts ueber die Prozessgrenze,
|
||||
# was das Ergebnis verfaelschen koennte.
|
||||
VORSPANN = """
|
||||
import json, time, resource
|
||||
import numpy as np
|
||||
|
||||
def instanz(m, n):
|
||||
rng = np.random.default_rng(20)
|
||||
kosten = rng.integers(5, 95, (m, n)).astype(float)
|
||||
angebot = rng.integers(50, 150, m).astype(float)
|
||||
bedarf = angebot.sum() * rng.dirichlet(np.ones(n))
|
||||
return kosten, angebot, bedarf
|
||||
|
||||
def speicher_mb():
|
||||
# ru_maxrss ist unter Linux in Kilobyte
|
||||
return resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024
|
||||
|
||||
M, N = {m}, {n}
|
||||
kosten, angebot, bedarf = instanz(M, N)
|
||||
"""
|
||||
|
||||
ANSAETZE = {
|
||||
"scipy.linprog": """
|
||||
from scipy.optimize import linprog
|
||||
t0 = time.perf_counter()
|
||||
c = kosten.reshape(-1)
|
||||
A_ub = np.zeros((M, M * N)); A_eq = np.zeros((N, M * N))
|
||||
for i in range(M):
|
||||
A_ub[i, i * N:(i + 1) * N] = 1.0
|
||||
for j in range(N):
|
||||
A_eq[j, j::N] = 1.0
|
||||
aufbau = time.perf_counter() - t0
|
||||
t0 = time.perf_counter()
|
||||
r = linprog(c=c, A_ub=A_ub, b_ub=angebot, A_eq=A_eq, b_eq=bedarf,
|
||||
bounds=(0, None), method="highs")
|
||||
loesen = time.perf_counter() - t0
|
||||
ausgabe = (float(r.fun), aufbau, loesen, speicher_mb())
|
||||
""",
|
||||
|
||||
"highspy": """
|
||||
import highspy
|
||||
t0 = time.perf_counter()
|
||||
h = highspy.Highs(); h.setOptionValue("output_flag", False)
|
||||
h.addVars(M * N, np.zeros(M * N), np.full(M * N, highspy.kHighsInf))
|
||||
for k in range(M * N):
|
||||
h.changeColCost(k, float(kosten.reshape(-1)[k]))
|
||||
for i in range(M):
|
||||
idx = np.arange(i * N, (i + 1) * N, dtype=np.int32)
|
||||
h.addRow(-highspy.kHighsInf, float(angebot[i]), N, idx, np.ones(N))
|
||||
for j in range(N):
|
||||
idx = np.arange(j, M * N, N, dtype=np.int32)
|
||||
h.addRow(float(bedarf[j]), float(bedarf[j]), M, idx, np.ones(M))
|
||||
aufbau = time.perf_counter() - t0
|
||||
t0 = time.perf_counter(); h.run(); loesen = time.perf_counter() - t0
|
||||
ausgabe = (h.getInfo().objective_function_value, aufbau, loesen, speicher_mb())
|
||||
""",
|
||||
|
||||
"ortools/GLOP": """
|
||||
from ortools.linear_solver import pywraplp
|
||||
t0 = time.perf_counter()
|
||||
s = pywraplp.Solver.CreateSolver("GLOP")
|
||||
x = [[s.NumVar(0, s.infinity(), f"x{i}_{j}") for j in range(N)]
|
||||
for i in range(M)]
|
||||
for i in range(M):
|
||||
s.Add(sum(x[i]) <= float(angebot[i]))
|
||||
for j in range(N):
|
||||
s.Add(sum(x[i][j] for i in range(M)) == float(bedarf[j]))
|
||||
s.Minimize(sum(float(kosten[i, j]) * x[i][j]
|
||||
for i in range(M) for j in range(N)))
|
||||
aufbau = time.perf_counter() - t0
|
||||
t0 = time.perf_counter(); s.Solve(); loesen = time.perf_counter() - t0
|
||||
ausgabe = (s.Objective().Value(), aufbau, loesen, speicher_mb())
|
||||
""",
|
||||
|
||||
"cvxpy": """
|
||||
import cvxpy as cp
|
||||
t0 = time.perf_counter()
|
||||
x = cp.Variable((M, N), nonneg=True)
|
||||
problem = cp.Problem(cp.Minimize(cp.sum(cp.multiply(kosten, x))),
|
||||
[cp.sum(x, axis=1) <= angebot,
|
||||
cp.sum(x, axis=0) == bedarf])
|
||||
aufbau = time.perf_counter() - t0
|
||||
t0 = time.perf_counter(); problem.solve(); loesen = time.perf_counter() - t0
|
||||
ausgabe = (float(problem.value), aufbau, loesen, speicher_mb())
|
||||
""",
|
||||
}
|
||||
|
||||
|
||||
def messe(name: str, quelltext: str, m: int, n: int):
|
||||
"""Fuehrt einen Ansatz in einem eigenen Prozess aus."""
|
||||
programm = (VORSPANN.format(m=m, n=n) + textwrap.dedent(quelltext)
|
||||
+ "\nprint(json.dumps(ausgabe))\n")
|
||||
ergebnis = subprocess.run([sys.executable, "-c", programm],
|
||||
capture_output=True, text=True, timeout=600)
|
||||
if ergebnis.returncode != 0:
|
||||
return None, ergebnis.stderr.strip().splitlines()[-1][:60]
|
||||
return json.loads(ergebnis.stdout.strip().splitlines()[-1]), None
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
print("=" * 92)
|
||||
print(" SKALIERUNGSVERGLEICH: TRANSPORTPROBLEM, VIER BIBLIOTHEKEN")
|
||||
print("=" * 92)
|
||||
print("Jede Zeile ein eigener Prozess. Zeiten und Speicher sind "
|
||||
"hardwareabhaengig,")
|
||||
print("die Zielwerte und ihr Verhaeltnis zueinander nicht.\n")
|
||||
|
||||
for m, n in GROESSEN:
|
||||
kopf = f"--- {m} Lager x {n} Kunden = {m * n:,} Variablen "
|
||||
print(kopf + "-" * max(3, 92 - len(kopf)))
|
||||
print(f" {'Bibliothek':<16} {'Zielwert':>14} {'Aufbau':>9} "
|
||||
f"{'Loesen':>9} {'Anteil':>8} {'Speicher':>10}")
|
||||
print(" " + "-" * 72)
|
||||
zielwerte = {}
|
||||
for name, quelltext in ANSAETZE.items():
|
||||
werte, fehler = messe(name, quelltext, m, n)
|
||||
if werte is None:
|
||||
print(f" {name:<16} nicht verfuegbar: {fehler}")
|
||||
continue
|
||||
ziel, aufbau, loesen, speicher = werte
|
||||
zielwerte[name] = ziel
|
||||
anteil = aufbau / (aufbau + loesen) * 100
|
||||
print(f" {name:<16} {ziel:>14,.2f} {aufbau:>8.3f}s "
|
||||
f"{loesen:>8.3f}s {anteil:>7.0f}% {speicher:>9.0f} MB")
|
||||
|
||||
# Die wichtigste Zeile: Rechnen alle dasselbe aus?
|
||||
spanne = max(zielwerte.values()) - min(zielwerte.values())
|
||||
bezug = max(abs(v) for v in zielwerte.values())
|
||||
print(f" {'':16} Spannweite der Zielwerte: {spanne:.2e} "
|
||||
f"(relativ {spanne / bezug:.1e})")
|
||||
if spanne / bezug > 1e-6:
|
||||
print(" ACHTUNG: Die Bibliotheken widersprechen sich - "
|
||||
"der Zeitvergleich ist wertlos.")
|
||||
print()
|
||||
|
||||
print("=" * 92)
|
||||
print(" WAS MAN AUS SO EINER TABELLE ABLESEN DARF - UND WAS NICHT")
|
||||
print("=" * 92)
|
||||
print("DARF man ablesen:")
|
||||
print(" * Die Spalte 'Anteil' - wie viel der Zeit in den AUFBAU geht statt")
|
||||
print(" ins Loesen. Wenn dort 80 % stehen, ist ein schnellerer Solver die")
|
||||
print(" falsche Antwort; dann gehoert das Modell vektorisiert aufgebaut")
|
||||
print(" (Kapitel Oekosystem).")
|
||||
print(" * Die Groessenordnung des Speicherbedarfs. Sie entscheidet, was auf")
|
||||
print(" einer bestimmten Maschine ueberhaupt laeuft.")
|
||||
print(" * Wie sich beides mit der Groesse ENTWICKELT. Der Trend ist")
|
||||
print(" uebertragbarer als der Absolutwert.")
|
||||
print()
|
||||
print("NICHT ablesen darf man:")
|
||||
print(" * 'Bibliothek X ist schneller als Y.' Gemessen wurde EIN")
|
||||
print(" Problemtyp in EINER Formulierung. Ein MILP, ein QP oder eine")
|
||||
print(" andere Modellierung desselben Problems koennen die Reihenfolge")
|
||||
print(" umdrehen.")
|
||||
print(" * Etwas ueber Ihre Maschine. Diese Zahlen stammen von einer")
|
||||
print(" anderen. Der Sinn des Programms ist, dass Sie es auf Ihrer")
|
||||
print(" laufen lassen.")
|
||||
print("=" * 92)
|
||||
180
OR_HTML_04/programme/Betriebsueberwachung.py
Normal file
180
OR_HTML_04/programme/Betriebsueberwachung.py
Normal file
|
|
@ -0,0 +1,180 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# Betriebsueberwachung.py
|
||||
"""
|
||||
Kapitel Praxisfallen: Das Modell, das still aufgehoert hat zu optimieren.
|
||||
|
||||
Ein Optimierungsmodell geht in Betrieb und laeuft jede Nacht. Es liefert
|
||||
zuverlaessig eine Zahl, die Job-Dauer bleibt konstant, es gibt keine
|
||||
Fehlermeldung. Alles sieht gut aus.
|
||||
|
||||
Waehrenddessen waechst das Unternehmen. Aus 15 Lagern werden 40, dann 80. Das
|
||||
Zeitlimit von drei Sekunden - vor zwei Jahren grosszuegig bemessen - reicht
|
||||
irgendwann nicht mehr. Der Solver bricht ab und liefert die beste Loesung, die
|
||||
er bis dahin gefunden hat.
|
||||
|
||||
Das ist nicht falsch. Es ist sogar richtig so. Aber es ist etwas anderes als
|
||||
das, was das Modell einmal geleistet hat - und niemand merkt es, weil im
|
||||
Protokoll nur die Kosten stehen.
|
||||
|
||||
Dieses Programm zeigt die Verschlechterung und die Ueberwachung, die sie
|
||||
sichtbar macht.
|
||||
|
||||
WICHTIG: ortools wird hier nicht importiert (Konflikt mit highspy, siehe
|
||||
Kapitel Oekosystem).
|
||||
|
||||
Benoetigt: numpy, highspy
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import time
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import highspy
|
||||
|
||||
ZEITLIMIT = 3.0 # das Budget, das der Nachtjob hat
|
||||
GAP_WARNSCHWELLE = 0.02 # ab 2 % Gap wollen wir es wissen
|
||||
ZEITAUSLASTUNG_WARNSCHWELLE = 0.9
|
||||
|
||||
|
||||
@dataclass
|
||||
class Laufprotokoll:
|
||||
"""Was nach JEDEM Produktivlauf protokolliert gehoert.
|
||||
|
||||
Die meisten Systeme schreiben nur 'kosten'. Genau deshalb faellt eine
|
||||
schleichende Verschlechterung ueber Monate nicht auf.
|
||||
"""
|
||||
zeitpunkt: str
|
||||
lager: int
|
||||
kunden: int
|
||||
status: str
|
||||
kosten: float
|
||||
gap: float
|
||||
laufzeit: float
|
||||
knoten: int
|
||||
|
||||
@property
|
||||
def warnungen(self) -> list[str]:
|
||||
meldungen = []
|
||||
if self.status != "Optimal":
|
||||
meldungen.append(f"nicht beweisbar optimal ({self.status})")
|
||||
if self.gap > GAP_WARNSCHWELLE:
|
||||
meldungen.append(f"Gap {self.gap:.1%} ueber Schwelle "
|
||||
f"{GAP_WARNSCHWELLE:.0%}")
|
||||
if self.laufzeit > ZEITAUSLASTUNG_WARNSCHWELLE * ZEITLIMIT:
|
||||
meldungen.append(f"Zeitbudget zu {self.laufzeit / ZEITLIMIT:.0%} "
|
||||
f"ausgeschoepft")
|
||||
return meldungen
|
||||
|
||||
|
||||
def plane_netzwerk(n_lager: int, n_kunden: int, zeitlimit: float,
|
||||
saat: int = 7) -> tuple[str, float, float, float, int]:
|
||||
"""Standortplanung wie in Kapitel MILP, nur mit wachsender Groesse."""
|
||||
rng = np.random.default_rng(saat)
|
||||
fixkosten = rng.uniform(3000, 9000, n_lager)
|
||||
transport = rng.uniform(5, 60, (n_lager, n_kunden))
|
||||
bedarf = rng.uniform(10, 60, n_kunden)
|
||||
kapazitaet = np.full(n_lager, bedarf.sum() * 0.22)
|
||||
|
||||
modell = highspy.Highs()
|
||||
modell.setOptionValue("output_flag", False)
|
||||
modell.setOptionValue("time_limit", zeitlimit)
|
||||
|
||||
anzahl_x = n_lager * n_kunden
|
||||
unendlich = highspy.kHighsInf
|
||||
modell.addVars(anzahl_x, np.zeros(anzahl_x), np.full(anzahl_x, unendlich))
|
||||
modell.addVars(n_lager, np.zeros(n_lager), np.ones(n_lager))
|
||||
for i in range(n_lager):
|
||||
modell.changeColIntegrality(anzahl_x + i, highspy.HighsVarType.kInteger)
|
||||
modell.changeColCost(anzahl_x + i, fixkosten[i])
|
||||
for j in range(n_kunden):
|
||||
modell.changeColCost(i * n_kunden + j, transport[i, j])
|
||||
|
||||
for j in range(n_kunden):
|
||||
index = np.array([i * n_kunden + j for i in range(n_lager)], dtype=np.int32)
|
||||
modell.addRow(bedarf[j], bedarf[j], len(index), index, np.ones(len(index)))
|
||||
for i in range(n_lager):
|
||||
index = np.array([i * n_kunden + j for j in range(n_kunden)] + [anzahl_x + i],
|
||||
dtype=np.int32)
|
||||
werte = np.concatenate([np.ones(n_kunden), [-kapazitaet[i]]])
|
||||
modell.addRow(-unendlich, 0.0, len(index), index, werte)
|
||||
|
||||
t0 = time.perf_counter()
|
||||
modell.run()
|
||||
laufzeit = time.perf_counter() - t0
|
||||
info = modell.getInfo()
|
||||
return (modell.modelStatusToString(modell.getModelStatus()),
|
||||
info.objective_function_value, info.mip_gap, laufzeit,
|
||||
info.mip_node_count)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
# So ist das Unternehmen ueber zwei Jahre gewachsen.
|
||||
entwicklung = [("2024-Q1", 15, 40), ("2024-Q3", 25, 70),
|
||||
("2025-Q1", 40, 110), ("2025-Q3", 60, 160),
|
||||
("2026-Q1", 80, 220)]
|
||||
|
||||
print("=" * 84)
|
||||
print(" WAS DAS PROTOKOLL ZEIGT - UND WAS ES ZEIGEN SOLLTE")
|
||||
print("=" * 84)
|
||||
print(f"Derselbe Nachtjob, unveraendert, mit {ZEITLIMIT:.0f} Sekunden Zeitlimit.\n")
|
||||
|
||||
protokolle = []
|
||||
for zeitpunkt, lager, kunden in entwicklung:
|
||||
status, kosten, gap, laufzeit, knoten = plane_netzwerk(
|
||||
lager, kunden, ZEITLIMIT)
|
||||
protokolle.append(Laufprotokoll(zeitpunkt, lager, kunden, status,
|
||||
kosten, gap, laufzeit, knoten))
|
||||
|
||||
print("So sieht das ueblich gefuehrte Protokoll aus:\n")
|
||||
print(f" {'Zeitpunkt':<10} {'Kosten':>12}")
|
||||
print(" " + "-" * 24)
|
||||
for p in protokolle:
|
||||
print(f" {p.zeitpunkt:<10} {p.kosten:>12,.0f}")
|
||||
print("\n -> Die Kosten steigen. Das Unternehmen waechst ja auch.")
|
||||
print(" Nichts an dieser Tabelle deutet auf ein Problem hin.")
|
||||
|
||||
print("\n" + "-" * 84)
|
||||
print("So sieht ein vollstaendiges Protokoll aus:\n")
|
||||
print(f" {'Zeitpunkt':<10} {'Groesse':>10} {'Status':>20} {'Kosten':>11} "
|
||||
f"{'Gap':>8} {'Zeit':>7}")
|
||||
print(" " + "-" * 80)
|
||||
for p in protokolle:
|
||||
groesse = f"{p.lager}x{p.kunden}"
|
||||
print(f" {p.zeitpunkt:<10} {groesse:>10} {p.status:>20} "
|
||||
f"{p.kosten:>11,.0f} {p.gap * 100:7.2f} % {p.laufzeit:6.2f}s")
|
||||
|
||||
print("\n" + "-" * 84)
|
||||
print("Und so sehen die Warnungen aus, die daraus folgen:\n")
|
||||
for p in protokolle:
|
||||
if p.warnungen:
|
||||
print(f" {p.zeitpunkt}: " + "; ".join(p.warnungen))
|
||||
else:
|
||||
print(f" {p.zeitpunkt}: in Ordnung")
|
||||
|
||||
letzte = protokolle[-1]
|
||||
erste_warnung = next(p for p in protokolle if p.warnungen)
|
||||
|
||||
print("\n" + "=" * 84)
|
||||
print(" WAS DA PASSIERT IST")
|
||||
print("=" * 84)
|
||||
print(f"Seit {erste_warnung.zeitpunkt} erreicht der Job das Zeitlimit und liefert")
|
||||
print("nicht mehr das Optimum, sondern die beste bis dahin gefundene Loesung.")
|
||||
print(f"Im letzten Lauf betraegt der Abstand zum Bestmoeglichen "
|
||||
f"{letzte.gap:.1%}.")
|
||||
print()
|
||||
print("Der Job ist nicht abgestuerzt. Er hat keine Fehlermeldung erzeugt.")
|
||||
print("Die Laufzeit ist sogar bemerkenswert STABIL geblieben - genau deshalb,")
|
||||
print("weil das Zeitlimit greift. Ein Ueberwachungssystem, das auf Abstuerze")
|
||||
print("und Laufzeitspitzen achtet, sieht hier nichts.")
|
||||
print()
|
||||
print("Drei Zahlen gehoeren deshalb in jedes Protokoll eines Optimierungsjobs:")
|
||||
print(" * der STATUS - 'Optimal' oder etwas anderes?")
|
||||
print(" * der GAP - wie weit ist die Loesung vom Bestmoeglichen entfernt?")
|
||||
print(" * die LAUFZEIT im Verhaeltnis zum Limit - wie nah am Anschlag?")
|
||||
print()
|
||||
print("Die dritte ist die frueheste Warnung: Sie steigt, lange bevor der Gap")
|
||||
print("sichtbar wird, und gibt Zeit zum Handeln.")
|
||||
print("=" * 84)
|
||||
194
OR_HTML_04/programme/Big_M_Falle.py
Normal file
194
OR_HTML_04/programme/Big_M_Falle.py
Normal file
|
|
@ -0,0 +1,194 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# Big_M_Falle.py
|
||||
"""
|
||||
Kapitel MILP: Was ein zu grosses Big-M wirklich anrichtet.
|
||||
|
||||
Lehrbuecher warnen vor grossem M mit dem Hinweis "die Relaxation wird
|
||||
schwach, der Solver braucht mehr Knoten". Das stimmt - ist aber die
|
||||
harmlosere Haelfte der Wahrheit. Die gefaehrlichere: Bei sehr grossem M
|
||||
kann der Solver eine Loesung als ganzzahlig ANNEHMEN, in der die
|
||||
Binaervariablen bei 1e-8 stehen. Dann liefern zugeschaltete Anlagen Ware
|
||||
aus, waehrend das Modell ihre Fixkosten mit 0 verbucht.
|
||||
|
||||
Beispiel: Standortplanung, 12 moegliche Lager, 40 Kunden.
|
||||
min sum_i fix_i * y_i + sum_ij kosten_ij * x_ij
|
||||
u.d.N. sum_i x_ij = bedarf_j (jeder Kunde beliefert)
|
||||
sum_j x_ij <= M_i * y_i (Lager offen, wenn es liefert)
|
||||
y_i binaer, x_ij >= 0
|
||||
|
||||
Zwei Laeufe:
|
||||
1. M knapp gewaehlt (= tatsaechliche Lagerkapazitaet) -> richtig
|
||||
2. M = 1e7-fach zu gross, Presolve abgeschaltet -> stilles Desaster
|
||||
Beide werden mit derselben Pruefung kontrolliert, die den Fall auffliegen laesst.
|
||||
|
||||
Benoetigt: numpy, highspy
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import time
|
||||
|
||||
import numpy as np
|
||||
import highspy
|
||||
|
||||
RNG = np.random.default_rng(11)
|
||||
|
||||
N_LAGER, N_KUNDE = 12, 40
|
||||
FIXKOSTEN = RNG.uniform(2000, 5000, N_LAGER)
|
||||
TRANSPORT = RNG.uniform(5, 40, (N_LAGER, N_KUNDE))
|
||||
BEDARF = RNG.uniform(10, 60, N_KUNDE)
|
||||
KAPAZITAET = BEDARF.sum() * 0.45 # jedes Lager schafft 45 % des Gesamtbedarfs
|
||||
|
||||
# Toleranz, ab der ein Solver eine Variable als ganzzahlig durchgehen laesst.
|
||||
# HiGHS und die meisten anderen verwenden 1e-6 als Standard.
|
||||
GANZZAHL_TOLERANZ = 1e-6
|
||||
|
||||
|
||||
def loese(big_m: float, presolve: str = "on") -> dict:
|
||||
"""Baut und loest das Standortmodell. Liefert Loesung und Solverkennzahlen."""
|
||||
hochschule = highspy.Highs()
|
||||
hochschule.setOptionValue("output_flag", False)
|
||||
hochschule.setOptionValue("presolve", presolve)
|
||||
hochschule.setOptionValue("time_limit", 300.0)
|
||||
|
||||
anzahl_x = N_LAGER * N_KUNDE
|
||||
unendlich = highspy.kHighsInf
|
||||
|
||||
# Spalten: erst alle x_ij, dann die y_i
|
||||
hochschule.addVars(anzahl_x, np.zeros(anzahl_x), np.full(anzahl_x, unendlich))
|
||||
hochschule.addVars(N_LAGER, np.zeros(N_LAGER), np.ones(N_LAGER))
|
||||
for i in range(N_LAGER):
|
||||
hochschule.changeColIntegrality(anzahl_x + i,
|
||||
highspy.HighsVarType.kInteger)
|
||||
hochschule.changeColCost(anzahl_x + i, FIXKOSTEN[i])
|
||||
for j in range(N_KUNDE):
|
||||
hochschule.changeColCost(i * N_KUNDE + j, TRANSPORT[i, j])
|
||||
|
||||
# Jeder Kunde wird genau beliefert
|
||||
for j in range(N_KUNDE):
|
||||
index = np.array([i * N_KUNDE + j for i in range(N_LAGER)], dtype=np.int32)
|
||||
hochschule.addRow(BEDARF[j], BEDARF[j], len(index), index,
|
||||
np.ones(len(index)))
|
||||
|
||||
# Die Kopplung: sum_j x_ij - M * y_i <= 0
|
||||
for i in range(N_LAGER):
|
||||
index = np.array([i * N_KUNDE + j for j in range(N_KUNDE)] + [anzahl_x + i],
|
||||
dtype=np.int32)
|
||||
werte = np.concatenate([np.ones(N_KUNDE), [-big_m]])
|
||||
hochschule.addRow(-unendlich, 0.0, len(index), index, werte)
|
||||
|
||||
t0 = time.perf_counter()
|
||||
hochschule.run()
|
||||
dauer = time.perf_counter() - t0
|
||||
|
||||
loesung = np.array(hochschule.getSolution().col_value)
|
||||
info = hochschule.getInfo()
|
||||
return {
|
||||
"x": loesung[:anzahl_x].reshape(N_LAGER, N_KUNDE),
|
||||
"y": loesung[anzahl_x:],
|
||||
"zielwert": info.objective_function_value,
|
||||
"knoten": info.mip_node_count,
|
||||
"dauer": dauer,
|
||||
}
|
||||
|
||||
|
||||
def pruefe(ergebnis: dict) -> tuple[bool, list[str]]:
|
||||
"""Die Pruefung, die in jedes MILP-Auswertungsskript gehoert.
|
||||
|
||||
Sie rechnet die Kosten AUS DER LOESUNG neu aus, statt dem Zielwert des
|
||||
Solvers zu glauben - und vergleicht beide. Genau diese Gegenrechnung
|
||||
entlarvt eine Loesung, in der Binaervariablen bei 1e-8 haengengeblieben
|
||||
sind.
|
||||
"""
|
||||
beanstandungen = []
|
||||
x, y = ergebnis["x"], ergebnis["y"]
|
||||
|
||||
# 1. Sind die Binaervariablen wirklich binaer?
|
||||
abstand = np.abs(y - np.round(y))
|
||||
if abstand.max() > GANZZAHL_TOLERANZ:
|
||||
beanstandungen.append(
|
||||
f"y ist nicht ganzzahlig: groesster Abstand {abstand.max():.2e}")
|
||||
|
||||
# 2. Liefert ein Lager, dessen Schalter aus ist?
|
||||
liefert = x.sum(axis=1) > 1e-6
|
||||
geschlossen_aber_aktiv = np.where(liefert & (y < 0.5))[0]
|
||||
if geschlossen_aber_aktiv.size:
|
||||
beanstandungen.append(
|
||||
f"Lager {geschlossen_aber_aktiv.tolist()} liefern Ware, "
|
||||
f"gelten im Modell aber als geschlossen")
|
||||
|
||||
# 3. Stimmt der Zielwert mit den echten Kosten ueberein?
|
||||
echte_fixkosten = FIXKOSTEN[liefert].sum()
|
||||
echte_transportkosten = float((TRANSPORT * x).sum())
|
||||
echte_kosten = echte_fixkosten + echte_transportkosten
|
||||
if abs(echte_kosten - ergebnis["zielwert"]) > 1e-4 * max(1.0, echte_kosten):
|
||||
beanstandungen.append(
|
||||
f"Zielwert {ergebnis['zielwert']:,.2f} weicht von den echten Kosten "
|
||||
f"{echte_kosten:,.2f} ab (Differenz {echte_kosten - ergebnis['zielwert']:,.2f})")
|
||||
|
||||
return not beanstandungen, beanstandungen
|
||||
|
||||
|
||||
def zeige(titel: str, ergebnis: dict) -> None:
|
||||
x, y = ergebnis["x"], ergebnis["y"]
|
||||
liefert = x.sum(axis=1) > 1e-6
|
||||
print(f"\n{titel}")
|
||||
print(f" Zielwert laut Solver {ergebnis['zielwert']:>14,.2f} EUR")
|
||||
print(f" Knoten / Zeit {ergebnis['knoten']:>14,} "
|
||||
f"/ {ergebnis['dauer']:.3f} s")
|
||||
print(f" Lager mit y = 1 {int((y > 0.5).sum()):>14}")
|
||||
print(f" Lager, die tatsaechlich liefern {int(liefert.sum()):>10}")
|
||||
unter = y[y < 0.5]
|
||||
print(f" groesster y-Wert unter 0.5 "
|
||||
f"{(unter.max() if unter.size else 0.0):>14.3e}")
|
||||
print(f" Fixkosten real / verbucht {FIXKOSTEN[liefert].sum():>14,.2f} "
|
||||
f"/ {float((FIXKOSTEN * y).sum()):,.2f} EUR")
|
||||
|
||||
in_ordnung, beanstandungen = pruefe(ergebnis)
|
||||
if in_ordnung:
|
||||
print(" PRUEFUNG: bestanden")
|
||||
else:
|
||||
print(" PRUEFUNG: DURCHGEFALLEN")
|
||||
for text in beanstandungen:
|
||||
print(f" - {text}")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
print("=" * 78)
|
||||
print(" DIE BIG-M-FALLE: STANDORTPLANUNG, 12 LAGER, 40 KUNDEN")
|
||||
print("=" * 78)
|
||||
print(f"Tatsaechliche Lagerkapazitaet: {KAPAZITAET:,.1f} Einheiten.")
|
||||
print("Genau das ist das kleinstmoegliche gueltige M - mehr kann ein Lager")
|
||||
print("ohnehin nicht ausliefern.")
|
||||
|
||||
knapp = loese(KAPAZITAET)
|
||||
zeige("[1] M = Kapazitaet (richtig gewaehlt)", knapp)
|
||||
|
||||
gross = loese(1e7 * KAPAZITAET, presolve="off")
|
||||
zeige("[2] M = 10 Millionen mal Kapazitaet, Presolve abgeschaltet", gross)
|
||||
|
||||
print("\n" + "=" * 78)
|
||||
print(" WAS DA PASSIERT IST")
|
||||
print("=" * 78)
|
||||
fehlbetrag = knapp["zielwert"] - gross["zielwert"]
|
||||
print(f"Lauf [2] meldet {gross['zielwert']:,.2f} EUR und sieht damit um")
|
||||
print(f"{fehlbetrag:,.2f} EUR BESSER aus als die richtige Loesung - ein Ergebnis,")
|
||||
print("ueber das sich jeder Auftraggeber freuen wuerde.")
|
||||
print()
|
||||
print("Der Grund steht in der Zeile 'groesster y-Wert unter 0.5': Die")
|
||||
print(f"Schaltervariablen stehen bei rund "
|
||||
f"{gross['y'][gross['y'] < 0.5].max():.0e}.")
|
||||
print(f"Das ist kleiner als die Ganzzahltoleranz {GANZZAHL_TOLERANZ:.0e}, also gilt")
|
||||
print("y = 0 - 'Lager geschlossen'. Zugleich ist M so gross, dass")
|
||||
print(" sum_j x_ij <= M * 1e-8")
|
||||
print("immer noch reichlich Liefermenge erlaubt. Die Lager liefern also,")
|
||||
print("ohne dass ihre Fixkosten je bezahlt werden. Der Fachbegriff dafuer")
|
||||
print("ist 'trickle flow'.")
|
||||
print()
|
||||
print("WICHTIG: Mit eingeschaltetem Presolve (Standard) faellt HiGHS hier")
|
||||
print("nicht darauf herein - es zieht M selbst zurecht. Verlassen Sie sich")
|
||||
print("nicht darauf: Presolve kann das nur, wenn eine implizite Schranke")
|
||||
print("herleitbar ist. Die Pruefung aus pruefe() kostet Millisekunden und")
|
||||
print("funktioniert immer.")
|
||||
print("=" * 78)
|
||||
96
OR_HTML_04/programme/Bot_Allokation.py
Normal file
96
OR_HTML_04/programme/Bot_Allokation.py
Normal file
|
|
@ -0,0 +1,96 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# Bot_Allokation.py
|
||||
"""
|
||||
Kapitel Einfuehrung: Erstes Optimierungsmodell mit Google OR-Tools (CP-SAT).
|
||||
Problem: Server-Allokation für Trading-Bots.
|
||||
|
||||
Modell mit expliziter Nebenbedingung statt versteckter Variablengrenze,
|
||||
vollständiger Statusauswertung und Ergebnisprüfung.
|
||||
"""
|
||||
|
||||
from ortools.sat.python import cp_model
|
||||
|
||||
# --- Parameter (die "Fakten" des Problems) ---------------------------------
|
||||
GEWINN_A, GEWINN_B = 150, 250 # Euro pro Bot und Tag
|
||||
CPU_A, CPU_B = 2, 5 # vCPU-Bedarf je Bot
|
||||
RAM_A, RAM_B = 4, 6 # GB Arbeitsspeicher je Bot
|
||||
CPU_GESAMT, RAM_GESAMT = 40, 60 # verfügbare Kapazitäten
|
||||
MAX_ARBITRAGE = 8 # Marktliquiditätsgrenze
|
||||
|
||||
|
||||
def loese_bot_allokation(zeitlimit_s: float = 10.0):
|
||||
"""Baut das Modell, löst es und gibt eine Auswertung aus."""
|
||||
|
||||
# 1. Modell instanziieren
|
||||
modell = cp_model.CpModel()
|
||||
|
||||
# 2. Entscheidungsvariablen anlegen
|
||||
# NewIntVar(untere_grenze, obere_grenze, name)
|
||||
# Die obere Grenze ist bewusst großzügig; die echte Beschränkung
|
||||
# formulieren wir unten als Nebenbedingung, damit das Modell die
|
||||
# mathematische Formulierung 1:1 abbildet.
|
||||
x_a = modell.NewIntVar(0, 100, "Arbitrage_Bots")
|
||||
x_b = modell.NewIntVar(0, 100, "Trend_Bots")
|
||||
|
||||
# 3. Nebenbedingungen definieren
|
||||
c_cpu = modell.Add(CPU_A * x_a + CPU_B * x_b <= CPU_GESAMT) # vCPU-Limit
|
||||
c_ram = modell.Add(RAM_A * x_a + RAM_B * x_b <= RAM_GESAMT) # RAM-Limit
|
||||
c_liq = modell.Add(x_a <= MAX_ARBITRAGE) # Marktliquidität
|
||||
|
||||
# 4. Zielfunktion: Maximiere den Tagesgewinn
|
||||
modell.Maximize(GEWINN_A * x_a + GEWINN_B * x_b)
|
||||
|
||||
# 5. Solver konfigurieren und ausführen
|
||||
loeser = cp_model.CpSolver()
|
||||
loeser.parameters.max_time_in_seconds = zeitlimit_s
|
||||
status = loeser.Solve(modell)
|
||||
|
||||
# 6. Status auswerten -- IMMER alle Fälle behandeln
|
||||
status_text = loeser.StatusName(status)
|
||||
if status not in (cp_model.OPTIMAL, cp_model.FEASIBLE):
|
||||
print(f"Keine verwertbare Lösung. Solver-Status: {status_text}")
|
||||
if status == cp_model.INFEASIBLE:
|
||||
print("Das Modell ist unlösbar - die Bedingungen widersprechen sich.")
|
||||
return None
|
||||
|
||||
a, b = loeser.Value(x_a), loeser.Value(x_b)
|
||||
gewinn = loeser.ObjectiveValue()
|
||||
|
||||
# 7. Ergebnis ausgeben
|
||||
print("=" * 58)
|
||||
print(" OPTIMALE BOT-ALLOKATION")
|
||||
print("=" * 58)
|
||||
print(f"Solver-Status: {status_text}"
|
||||
f"{' (beweisbar optimal)' if status == cp_model.OPTIMAL else ' (zulässig, nicht bewiesen)'}")
|
||||
print(f"Arbitrage-Bots (x_A): {a} Instanzen")
|
||||
print(f"Trendfolge-Bots (x_B): {b} Instanzen")
|
||||
print(f"Täglicher Max-Gewinn: {gewinn:,.2f} EUR")
|
||||
|
||||
genutzt_cpu = CPU_A * a + CPU_B * b
|
||||
genutzt_ram = RAM_A * a + RAM_B * b
|
||||
print("-" * 58)
|
||||
print(f"CPU-Auslastung: {genutzt_cpu} / {CPU_GESAMT} vCPUs "
|
||||
f"({genutzt_cpu / CPU_GESAMT * 100:.1f} %)")
|
||||
print(f"RAM-Auslastung: {genutzt_ram} / {RAM_GESAMT} GB "
|
||||
f"({genutzt_ram / RAM_GESAMT * 100:.1f} %)")
|
||||
print(f"Arbitrage-Limit: {a} / {MAX_ARBITRAGE}")
|
||||
|
||||
# 8. Selbstkontrolle: Sind wirklich alle Bedingungen eingehalten?
|
||||
assert genutzt_cpu <= CPU_GESAMT, "vCPU-Limit verletzt!"
|
||||
assert genutzt_ram <= RAM_GESAMT, "RAM-Limit verletzt!"
|
||||
assert a <= MAX_ARBITRAGE, "Arbitrage-Limit verletzt!"
|
||||
|
||||
# 9. Zum Vergleich: die naive Strategie "nur der lukrativere Bot"
|
||||
naiv_b = min(CPU_GESAMT // CPU_B, RAM_GESAMT // RAM_B)
|
||||
naiv_gewinn = GEWINN_B * naiv_b
|
||||
print("-" * 58)
|
||||
print(f"Naive Strategie (nur Trendfolge): {naiv_b} Bots -> {naiv_gewinn:,.2f} EUR")
|
||||
print(f"Vorteil der Optimierung: {gewinn - naiv_gewinn:,.2f} EUR pro Tag "
|
||||
f"({(gewinn / naiv_gewinn - 1) * 100:.1f} %)")
|
||||
print("=" * 58)
|
||||
return a, b, gewinn
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
loese_bot_allokation()
|
||||
103
OR_HTML_04/programme/Brute_Force_Vergleich.py
Normal file
103
OR_HTML_04/programme/Brute_Force_Vergleich.py
Normal file
|
|
@ -0,0 +1,103 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# Brute_Force_Vergleich.py
|
||||
"""
|
||||
Kapitel Einfuehrung: Warum Ausprobieren scheitert.
|
||||
Vergleicht vollständige Enumeration mit einem Constraint-Solver an einem
|
||||
Zuordnungsproblem wachsender Größe (n Mitarbeiter auf n Aufgaben).
|
||||
"""
|
||||
|
||||
import itertools
|
||||
import math
|
||||
import time
|
||||
|
||||
import numpy as np
|
||||
from ortools.sat.python import cp_model
|
||||
|
||||
|
||||
def erzeuge_kostenmatrix(n: int, seed: int = 7) -> np.ndarray:
|
||||
"""Zufällige, aber reproduzierbare Kosten: Wer bearbeitet welche Aufgabe wie teuer?"""
|
||||
rng = np.random.default_rng(seed)
|
||||
return rng.integers(low=10, high=100, size=(n, n))
|
||||
|
||||
|
||||
def loese_brute_force(kosten: np.ndarray) -> tuple[float, tuple, int]:
|
||||
"""
|
||||
Probiert ALLE n! Zuordnungen durch und behält die beste.
|
||||
Rückgabe: (bester Kostenwert, beste Permutation, Anzahl geprüfter Kombinationen)
|
||||
"""
|
||||
n = len(kosten)
|
||||
bester_wert = math.inf
|
||||
beste_zuordnung = None
|
||||
geprueft = 0
|
||||
|
||||
# itertools.permutations(range(n)) liefert nacheinander jede Reihenfolge
|
||||
for zuordnung in itertools.permutations(range(n)):
|
||||
# zuordnung[i] = Aufgabe, die Mitarbeiter i übernimmt
|
||||
wert = sum(kosten[i][zuordnung[i]] for i in range(n))
|
||||
geprueft += 1
|
||||
if wert < bester_wert:
|
||||
bester_wert = wert
|
||||
beste_zuordnung = zuordnung
|
||||
|
||||
return bester_wert, beste_zuordnung, geprueft
|
||||
|
||||
|
||||
def loese_mit_solver(kosten: np.ndarray) -> tuple[float, tuple]:
|
||||
"""Dasselbe Problem als Constraint-Programm — der Solver probiert NICHT alles durch."""
|
||||
n = len(kosten)
|
||||
modell = cp_model.CpModel()
|
||||
|
||||
# x[i][j] = 1 <=> Mitarbeiter i übernimmt Aufgabe j
|
||||
x = [[modell.NewBoolVar(f"x_{i}_{j}") for j in range(n)] for i in range(n)]
|
||||
|
||||
for i in range(n):
|
||||
modell.AddExactlyOne(x[i][j] for j in range(n)) # jeder genau eine Aufgabe
|
||||
for j in range(n):
|
||||
modell.AddExactlyOne(x[i][j] for i in range(n)) # jede Aufgabe genau einmal
|
||||
|
||||
modell.Minimize(sum(int(kosten[i][j]) * x[i][j] for i in range(n) for j in range(n)))
|
||||
|
||||
loeser = cp_model.CpSolver()
|
||||
loeser.parameters.max_time_in_seconds = 30.0
|
||||
status = loeser.Solve(modell)
|
||||
if status not in (cp_model.OPTIMAL, cp_model.FEASIBLE):
|
||||
raise RuntimeError("Solver fand keine Lösung.")
|
||||
|
||||
zuordnung = tuple(
|
||||
next(j for j in range(n) if loeser.Value(x[i][j]) == 1) for i in range(n)
|
||||
)
|
||||
return loeser.ObjectiveValue(), zuordnung
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
print("=" * 72)
|
||||
print(" VOLLSTÄNDIGES DURCHPROBIEREN vs. OPERATIONS RESEARCH")
|
||||
print("=" * 72)
|
||||
print(f"{'n':>3} | {'Kombinationen':>18} | {'Brute Force':>14} | {'Solver':>10} | gleich?")
|
||||
print("-" * 72)
|
||||
|
||||
for n in [4, 6, 8, 9, 10]: # n = 11 dauert bereits ~2 Minuten
|
||||
kosten = erzeuge_kostenmatrix(n)
|
||||
|
||||
t0 = time.perf_counter()
|
||||
wert_bf, zuord_bf, geprueft = loese_brute_force(kosten)
|
||||
zeit_bf = time.perf_counter() - t0
|
||||
|
||||
t0 = time.perf_counter()
|
||||
wert_or, zuord_or = loese_mit_solver(kosten)
|
||||
zeit_or = time.perf_counter() - t0
|
||||
|
||||
gleich = "ja" if abs(wert_bf - wert_or) < 1e-9 else "NEIN!"
|
||||
print(f"{n:>3} | {geprueft:>18,} | {zeit_bf:>11.4f} s | {zeit_or:>7.4f} s | {gleich}")
|
||||
|
||||
print("-" * 72)
|
||||
# Hochrechnung: Wie lange bräuchte Brute Force bei n = 20?
|
||||
kombis_20 = math.factorial(20)
|
||||
# Messbasis: Prüfungen pro Sekunde aus dem letzten Lauf schätzen
|
||||
pro_sekunde = geprueft / max(zeit_bf, 1e-9)
|
||||
jahre = kombis_20 / pro_sekunde / (60 * 60 * 24 * 365.25)
|
||||
print(f"Hochrechnung für n = 20: {kombis_20:,} Kombinationen")
|
||||
print(f"Bei gemessenen {pro_sekunde:,.0f} Prüfungen/s wären das {jahre:,.0f} Jahre.")
|
||||
print("Der Solver löst dieselbe Instanz in Sekundenbruchteilen.")
|
||||
print("=" * 72)
|
||||
223
OR_HTML_04/programme/CP_SAT_Statusfaelle.py
Normal file
223
OR_HTML_04/programme/CP_SAT_Statusfaelle.py
Normal file
|
|
@ -0,0 +1,223 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# CP_SAT_Statusfaelle.py
|
||||
"""
|
||||
Kapitel CP-SAT: Alle fuenf Antworten von CP-SAT - und was jede bedeutet.
|
||||
|
||||
Der haeufigste Fehler beim Einsatz von CP-SAT im Betrieb ist eine Zeile wie
|
||||
|
||||
if status == cp_model.OPTIMAL:
|
||||
...verwerte die Loesung...
|
||||
|
||||
Sie wirft in drei von fuenf Faellen etwas weg, das man haette gebrauchen
|
||||
koennen, und verschweigt in einem weiteren Fall einen Modellierungsfehler.
|
||||
|
||||
Dieses Programm erzeugt jeden Statusfall absichtlich und zeigt die passende
|
||||
Reaktion:
|
||||
|
||||
OPTIMAL beweisbar bestmoeglich -> ausfuehren
|
||||
FEASIBLE zulaessig, Zeit war um -> Gap pruefen, dann entscheiden
|
||||
INFEASIBLE kein Plan existiert -> harte Regeln lockern
|
||||
MODEL_INVALID das Modell ist fehlerhaft -> Programmierfehler beheben
|
||||
UNKNOWN nichts gefunden, Zeit war um -> mehr Zeit oder Heuristik
|
||||
|
||||
Benoetigt: numpy, ortools
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
from ortools.sat.python import cp_model
|
||||
|
||||
|
||||
@dataclass
|
||||
class Auswertung:
|
||||
"""Was nach einem CP-SAT-Lauf tatsaechlich feststeht."""
|
||||
status: str
|
||||
hat_loesung: bool
|
||||
beweisbar_optimal: bool
|
||||
zielwert: float | None
|
||||
schranke: float | None
|
||||
gap: float | None
|
||||
empfehlung: str
|
||||
|
||||
|
||||
def werte_aus(loeser: cp_model.CpSolver, status: int,
|
||||
mit_zielfunktion: bool) -> Auswertung:
|
||||
"""Uebersetzt einen CP-SAT-Status in eine Handlungsempfehlung.
|
||||
|
||||
Genau diese Funktion fehlt in den meisten Projekten. Sie ist kurz, aber
|
||||
sie ist der Unterschied zwischen einem Prototyp und einem System, das
|
||||
nachts ohne Aufsicht laeuft.
|
||||
"""
|
||||
name = loeser.StatusName(status)
|
||||
hat_loesung = status in (cp_model.OPTIMAL, cp_model.FEASIBLE)
|
||||
|
||||
zielwert = schranke = gap = None
|
||||
if hat_loesung and mit_zielfunktion:
|
||||
zielwert = loeser.ObjectiveValue()
|
||||
schranke = loeser.BestObjectiveBound()
|
||||
nenner = max(abs(zielwert), 1e-9)
|
||||
gap = abs(zielwert - schranke) / nenner
|
||||
|
||||
if status == cp_model.OPTIMAL:
|
||||
empfehlung = "Ausfuehren. Besser geht es nachweislich nicht."
|
||||
elif status == cp_model.FEASIBLE:
|
||||
empfehlung = (f"Brauchbar. Der Plan ist hoechstens {gap * 100:.2f} % "
|
||||
f"schlechter als das theoretisch Bestmoegliche."
|
||||
if gap is not None else
|
||||
"Brauchbar, aber nicht als optimal bewiesen.")
|
||||
elif status == cp_model.INFEASIBLE:
|
||||
empfehlung = ("KEIN Plan existiert. Das ist eine Aussage ueber Ihre "
|
||||
"harten Regeln, nicht ueber die Rechenzeit: Auch "
|
||||
"unendlich viel Zeit wuerde nichts aendern. Regeln "
|
||||
"lockern oder weich machen (Anhang C).")
|
||||
elif status == cp_model.MODEL_INVALID:
|
||||
empfehlung = ("Das MODELL ist fehlerhaft, nicht das Problem. "
|
||||
"Typisch: eine Variable mit leerem Wertebereich. "
|
||||
"model.Validate() nennt die Ursache.")
|
||||
else: # UNKNOWN
|
||||
empfehlung = ("Nichts gefunden - aber auch nicht widerlegt. Mehr Zeit "
|
||||
"geben, das Modell vereinfachen oder mit einer "
|
||||
"Heuristik vorbelegen (AddHint).")
|
||||
|
||||
return Auswertung(
|
||||
status=name,
|
||||
hat_loesung=hat_loesung,
|
||||
beweisbar_optimal=status == cp_model.OPTIMAL,
|
||||
zielwert=zielwert, schranke=schranke, gap=gap,
|
||||
empfehlung=empfehlung,
|
||||
)
|
||||
|
||||
|
||||
def loese(baue, zeitlimit: float, mit_zielfunktion: bool) -> Auswertung:
|
||||
modell = cp_model.CpModel()
|
||||
baue(modell)
|
||||
loeser = cp_model.CpSolver()
|
||||
loeser.parameters.max_time_in_seconds = zeitlimit
|
||||
# Fuer reproduzierbare Buchausgaben; im Betrieb weglassen.
|
||||
loeser.parameters.num_workers = 1
|
||||
loeser.parameters.random_seed = 1
|
||||
status = loeser.Solve(modell)
|
||||
return werte_aus(loeser, status, mit_zielfunktion)
|
||||
|
||||
|
||||
# --- Die fuenf Faelle --------------------------------------------------------
|
||||
|
||||
def dienstplan_loesbar(modell: cp_model.CpModel) -> None:
|
||||
"""Vier Personen, fuenf Dienste, hoechstens zwei je Person - loesbar."""
|
||||
personen, dienste = 4, 5
|
||||
x = [[modell.NewBoolVar(f"x_{i}_{j}") for j in range(dienste)]
|
||||
for i in range(personen)]
|
||||
for j in range(dienste):
|
||||
modell.AddExactlyOne(x[i][j] for i in range(personen))
|
||||
for i in range(personen):
|
||||
modell.Add(sum(x[i]) <= 2)
|
||||
|
||||
|
||||
def dienstplan_unloesbar(modell: cp_model.CpModel) -> None:
|
||||
"""Dieselben fuenf Dienste, aber hoechstens EIN Dienst je Person.
|
||||
|
||||
Vier Personen koennen zusammen hoechstens vier Dienste uebernehmen - fuer
|
||||
fuenf Dienste reicht das nicht. Kein Solver der Welt findet hier etwas.
|
||||
"""
|
||||
personen, dienste = 4, 5
|
||||
x = [[modell.NewBoolVar(f"x_{i}_{j}") for j in range(dienste)]
|
||||
for i in range(personen)]
|
||||
for j in range(dienste):
|
||||
modell.AddExactlyOne(x[i][j] for i in range(personen))
|
||||
for i in range(personen):
|
||||
modell.Add(sum(x[i]) <= 1)
|
||||
|
||||
|
||||
def modell_fehlerhaft(modell: cp_model.CpModel) -> None:
|
||||
"""Eine Variable mit leerem Wertebereich: untere Schranke > obere.
|
||||
|
||||
Im Alltag entsteht das durch einen Rechenfehler in den Grenzen, etwa
|
||||
NewIntVar(mindestbesetzung, kapazitaet, ...) bei falsch sortierten Daten.
|
||||
"""
|
||||
kaputt = modell.NewIntVar(5, 2, "leerer_bereich")
|
||||
modell.Add(kaputt >= 0)
|
||||
|
||||
|
||||
def lastverteilung(modell: cp_model.CpModel) -> None:
|
||||
"""120 Auftraege auf 11 Maschinen - gross genug, dass der Optimalitaets-
|
||||
beweis laenger dauert als das Finden einer sehr guten Loesung."""
|
||||
rng = np.random.default_rng(3)
|
||||
dauer = rng.integers(100, 9000, 120)
|
||||
n, maschinen = len(dauer), 11
|
||||
obergrenze = int(dauer.sum())
|
||||
|
||||
x = [[modell.NewBoolVar(f"x_{i}_{k}") for k in range(maschinen)]
|
||||
for i in range(n)]
|
||||
for i in range(n):
|
||||
modell.AddExactlyOne(x[i])
|
||||
belegung = [modell.NewIntVar(0, obergrenze, f"l_{k}") for k in range(maschinen)]
|
||||
for k in range(maschinen):
|
||||
modell.Add(belegung[k] == sum(int(dauer[i]) * x[i][k] for i in range(n)))
|
||||
makespan = modell.NewIntVar(0, obergrenze, "makespan")
|
||||
modell.AddMaxEquality(makespan, belegung)
|
||||
modell.Minimize(makespan)
|
||||
|
||||
|
||||
def zahlpartition(modell: cp_model.CpModel) -> None:
|
||||
"""46 grosse Zahlen exakt in zwei gleich schwere Haelften teilen.
|
||||
|
||||
Ein beruehmt schweres Problem: Es gibt keine Zielfunktion, entweder es
|
||||
geht auf oder nicht - und beides zu entscheiden dauert lange.
|
||||
"""
|
||||
rng = np.random.default_rng(1)
|
||||
gewichte = rng.integers(10**6, 2 * 10**6, 46)
|
||||
x = [modell.NewBoolVar(f"x_{i}") for i in range(len(gewichte))]
|
||||
modell.Add(sum(int(gewichte[i]) * x[i] for i in range(len(gewichte)))
|
||||
== int(gewichte.sum() // 2))
|
||||
|
||||
|
||||
def zeige(titel: str, a: Auswertung) -> None:
|
||||
print(f"\n{titel}")
|
||||
print(f" Status {a.status}")
|
||||
print(f" Loesung vorhanden {'ja' if a.hat_loesung else 'nein'}")
|
||||
if a.zielwert is not None:
|
||||
print(f" Zielwert / Schranke {a.zielwert:,.0f} / {a.schranke:,.0f}"
|
||||
f" (Gap {a.gap * 100:.3f} %)")
|
||||
print(f" -> {a.empfehlung}")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
print("=" * 78)
|
||||
print(" DIE FUENF ANTWORTEN VON CP-SAT")
|
||||
print("=" * 78)
|
||||
|
||||
zeige("[1] Dienstplan, hoechstens 2 Dienste je Person",
|
||||
loese(dienstplan_loesbar, 10.0, mit_zielfunktion=False))
|
||||
|
||||
zeige("[2] Lastverteilung 120 Auftraege / 11 Maschinen, 0,5 s Limit",
|
||||
loese(lastverteilung, 0.5, mit_zielfunktion=True))
|
||||
|
||||
zeige("[3] Derselbe Fall mit 10 s Limit",
|
||||
loese(lastverteilung, 10.0, mit_zielfunktion=True))
|
||||
|
||||
zeige("[4] Dienstplan, hoechstens 1 Dienst je Person",
|
||||
loese(dienstplan_unloesbar, 10.0, mit_zielfunktion=False))
|
||||
|
||||
zeige("[5] Variable mit leerem Wertebereich",
|
||||
loese(modell_fehlerhaft, 10.0, mit_zielfunktion=False))
|
||||
|
||||
zeige("[6] Zahlpartition mit 46 grossen Zahlen, 2 s Limit",
|
||||
loese(zahlpartition, 2.0, mit_zielfunktion=False))
|
||||
|
||||
print("\n" + "=" * 78)
|
||||
print(" DIE DREI, DIE MAN NICHT VERWECHSELN DARF")
|
||||
print("=" * 78)
|
||||
print("INFEASIBLE Es gibt keine Loesung. Eine Aussage ueber Ihr MODELL.")
|
||||
print(" Mehr Rechenzeit aendert daran nichts.")
|
||||
print("UNKNOWN Es wurde keine gefunden. Eine Aussage ueber die ZEIT.")
|
||||
print(" Ob es eine gibt, ist offen.")
|
||||
print("MODEL_INVALID Das Modell ist gar kein gueltiges Modell. Eine")
|
||||
print(" Aussage ueber Ihren CODE - immer ein Programmierfehler.")
|
||||
print()
|
||||
print("Wer diese drei in einem 'else: return None' zusammenfasst, verliert")
|
||||
print("genau die Information, die zur Fehlersuche noetig waere.")
|
||||
print("=" * 78)
|
||||
242
OR_HTML_04/programme/CP_SAT_Vertretungssystem.py
Normal file
242
OR_HTML_04/programme/CP_SAT_Vertretungssystem.py
Normal file
|
|
@ -0,0 +1,242 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# CP_SAT_Vertretungssystem.py
|
||||
"""
|
||||
Kapitel CP-SAT: Dynamisches Vertretungs- und Einsatzplanungssystem mit CP-SAT.
|
||||
|
||||
Eigenschaften:
|
||||
* Solver-Status wird ueber SolverStatus aus or_kern.py ausgewertet - dieselbe
|
||||
Fallunterscheidung wie bei jedem anderen Solver im Buch
|
||||
* Strafkosten-Zerlegung wird ausgewiesen (Erklaerbarkeit, Kapitel Praxisfallen)
|
||||
* Fairness-Kriterium ergaenzt (Spannweite der Belastung minimieren)
|
||||
* Vorabdiagnose auf offensichtliche Unloesbarkeit
|
||||
* Abnahmepruefung des fertigen Plans OHNE den Solver zu fragen
|
||||
|
||||
Benoetigt: ortools, pandas, pydantic (ueber or_kern)
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from ortools.sat.python import cp_model
|
||||
import pandas as pd
|
||||
|
||||
from or_kern import Loesung, SolverStatus, status_von_cpsat
|
||||
|
||||
# --- 1. Datenbasis --------------------------------------------------------
|
||||
SLOTS = [1, 2, 3, 4]
|
||||
FAECHER = ["Mathematik", "Physik", "Mathematik", "Informatik"]
|
||||
KLASSEN = ["Klasse 8a", "Klasse 10b", "Klasse 7a", "Klasse 9c"]
|
||||
|
||||
PERSONAL = ["Frau_Mueller", "Herr_Schmidt", "Frau_Albrecht", "Herr_Bauer", "Frau_Koch"]
|
||||
|
||||
QUALIFIKATION = {
|
||||
"Frau_Mueller": {"Mathematik", "Physik"},
|
||||
"Herr_Schmidt": {"Informatik", "Mathematik"},
|
||||
"Frau_Albrecht": {"Physik", "Informatik"},
|
||||
"Herr_Bauer": {"Mathematik"},
|
||||
"Frau_Koch": {"Mathematik", "Physik", "Informatik"},
|
||||
}
|
||||
|
||||
VORBELASTUNG = { # bereits geleistete Stunden heute
|
||||
"Frau_Mueller": 1, "Herr_Schmidt": 0, "Frau_Albrecht": 2,
|
||||
"Herr_Bauer": 0, "Frau_Koch": 1,
|
||||
}
|
||||
|
||||
MAX_VERTRETUNGEN = 2
|
||||
STRAFE_VORBELASTUNG = 50 # je Stunde Vorbelastung und Zuweisung
|
||||
STRAFE_LOCH = 80 # je Freistundenloch
|
||||
STRAFE_UNFAIRNESS = 30 # je Einheit Spannweite der Gesamtbelastung
|
||||
|
||||
|
||||
def pruefe_grundsaetzliche_loesbarkeit() -> bool:
|
||||
"""Vorabdiagnose: Gibt es fuer jeden Slot ueberhaupt qualifiziertes Personal?"""
|
||||
ok = True
|
||||
for s, fach in enumerate(FAECHER):
|
||||
kandidaten = [p for p in PERSONAL if fach in QUALIFIKATION[p]]
|
||||
if not kandidaten:
|
||||
print(f" UNLOESBAR: Fuer Slot {SLOTS[s]} ({fach}) gibt es niemanden.")
|
||||
ok = False
|
||||
kapazitaet = len(PERSONAL) * MAX_VERTRETUNGEN
|
||||
if kapazitaet < len(SLOTS):
|
||||
print(f" UNLOESBAR: Kapazitaet {kapazitaet} < {len(SLOTS)} offene Stunden.")
|
||||
ok = False
|
||||
return ok
|
||||
|
||||
|
||||
def pruefe_plan(plan: dict[int, str]) -> list[str]:
|
||||
"""Prueft den fertigen Plan gegen die harten Regeln - ohne den Solver.
|
||||
|
||||
Der Solver kann nur pruefen, was ihm gesagt wurde. Diese Funktion prueft
|
||||
gegen die ANFORDERUNG und benutzt dafuer bewusst keine Modellvariable
|
||||
(siehe Kapitel Praxisfallen). Leere Liste heisst bestanden.
|
||||
"""
|
||||
beanstandungen: list[str] = []
|
||||
if sorted(plan) != list(range(len(SLOTS))):
|
||||
beanstandungen.append(f"nicht jede Stunde genau einmal besetzt: {sorted(plan)}")
|
||||
return beanstandungen
|
||||
for s, person in plan.items():
|
||||
if FAECHER[s] not in QUALIFIKATION[person]:
|
||||
beanstandungen.append(f"{person} ist nicht fuer {FAECHER[s]} qualifiziert")
|
||||
for person in PERSONAL:
|
||||
anzahl = sum(1 for p in plan.values() if p == person)
|
||||
if anzahl > MAX_VERTRETUNGEN:
|
||||
beanstandungen.append(f"{person} hat {anzahl} Vertretungen "
|
||||
f"(hoechstens {MAX_VERTRETUNGEN})")
|
||||
return beanstandungen
|
||||
|
||||
|
||||
def plane_vertretung() -> list[dict[str, str]] | None:
|
||||
if not pruefe_grundsaetzliche_loesbarkeit():
|
||||
return None
|
||||
|
||||
modell = cp_model.CpModel()
|
||||
|
||||
# Entscheidungsvariablen: x[person, slot] = 1 <=> Person uebernimmt Slot
|
||||
x = {(p, s): modell.NewBoolVar(f"zuweisung_{p}_slot{s+1}")
|
||||
for p in PERSONAL for s in range(len(SLOTS))}
|
||||
|
||||
# --- HARTE NEBENBEDINGUNGEN ------------------------------------------
|
||||
# H1: Jede Stunde genau einmal besetzen
|
||||
for s in range(len(SLOTS)):
|
||||
modell.AddExactlyOne(x[p, s] for p in PERSONAL)
|
||||
|
||||
# H2: Qualifikation - unqualifizierte Zuweisung ausschliessen
|
||||
for s, fach in enumerate(FAECHER):
|
||||
for p in PERSONAL:
|
||||
if fach not in QUALIFIKATION[p]:
|
||||
modell.Add(x[p, s] == 0)
|
||||
|
||||
# H3: Hoechstens MAX_VERTRETUNGEN Stunden je Person
|
||||
for p in PERSONAL:
|
||||
modell.Add(sum(x[p, s] for s in range(len(SLOTS))) <= MAX_VERTRETUNGEN)
|
||||
|
||||
# (H4 "nicht an zwei Orten gleichzeitig" ist durch H1 bereits erfuellt:
|
||||
# jeder Slot hat genau eine Person, und eine Person kann pro Slot nur
|
||||
# eine Variable auf 1 setzen.)
|
||||
|
||||
# --- WEICHE ZIELE ----------------------------------------------------
|
||||
strafterme = []
|
||||
|
||||
# W1: Vorbelastete Personen schonen
|
||||
strafe_last = []
|
||||
for p in PERSONAL:
|
||||
zuweisungen = sum(x[p, s] for s in range(len(SLOTS)))
|
||||
strafe_last.append(VORBELASTUNG[p] * STRAFE_VORBELASTUNG * zuweisungen)
|
||||
strafterme.extend(strafe_last)
|
||||
|
||||
# W2: Freistundenloecher vermeiden (Muster: arbeitet, frei, arbeitet)
|
||||
loch_variablen = []
|
||||
for p in PERSONAL:
|
||||
for s in range(len(SLOTS) - 2):
|
||||
loch = modell.NewBoolVar(f"loch_{p}_{s}")
|
||||
# Reifizierung: loch == 1 <=> (x_s AND NOT x_{s+1} AND x_{s+2})
|
||||
modell.AddBoolAnd([x[p, s], x[p, s + 1].Not(), x[p, s + 2]]).OnlyEnforceIf(loch)
|
||||
modell.AddBoolOr([x[p, s].Not(), x[p, s + 1], x[p, s + 2].Not()]) \
|
||||
.OnlyEnforceIf(loch.Not())
|
||||
loch_variablen.append(loch)
|
||||
strafterme.append(loch * STRAFE_LOCH)
|
||||
|
||||
# W3: Fairness - Spannweite der Gesamtbelastung minimieren
|
||||
gesamtlast = {}
|
||||
for p in PERSONAL:
|
||||
last = modell.NewIntVar(0, len(SLOTS) + max(VORBELASTUNG.values()), f"last_{p}")
|
||||
modell.Add(last == VORBELASTUNG[p] + sum(x[p, s] for s in range(len(SLOTS))))
|
||||
gesamtlast[p] = last
|
||||
max_last = modell.NewIntVar(0, 10, "max_last")
|
||||
min_last = modell.NewIntVar(0, 10, "min_last")
|
||||
modell.AddMaxEquality(max_last, list(gesamtlast.values()))
|
||||
modell.AddMinEquality(min_last, list(gesamtlast.values()))
|
||||
spannweite = modell.NewIntVar(0, 10, "spannweite")
|
||||
modell.Add(spannweite == max_last - min_last)
|
||||
strafterme.append(spannweite * STRAFE_UNFAIRNESS)
|
||||
|
||||
modell.Minimize(sum(strafterme))
|
||||
|
||||
# --- LOESEN -----------------------------------------------------------
|
||||
loeser = cp_model.CpSolver()
|
||||
loeser.parameters.max_time_in_seconds = 5.0
|
||||
# Reproduzierbarkeit: ein Arbeiter, fester Startwert. Ohne das kann CP-SAT
|
||||
# bei mehreren gleich guten Plaenen von Lauf zu Lauf einen anderen liefern
|
||||
# (siehe JobShop_Intervalle.py im naechsten Abschnitt).
|
||||
loeser.parameters.num_workers = 1
|
||||
loeser.parameters.random_seed = 1
|
||||
rohstatus = loeser.Solve(modell)
|
||||
|
||||
# Der Rueckgabewert wird in die gemeinsame Sprache uebersetzt. Ab hier
|
||||
# sieht die Auswertung genauso aus wie bei HiGHS, GLOP oder CVXPY.
|
||||
status = status_von_cpsat(rohstatus)
|
||||
loesung = Loesung(
|
||||
status=status,
|
||||
zielwert=loeser.ObjectiveValue() if status.brauchbar else None,
|
||||
schranke=loeser.BestObjectiveBound() if status.brauchbar else None,
|
||||
laufzeit=loeser.WallTime())
|
||||
|
||||
print("=" * 78)
|
||||
print(" DYNAMISCHER VERTRETUNGSPLAN (CP-SAT OPTIMIERT)")
|
||||
print("=" * 78)
|
||||
|
||||
if status.modellfehler:
|
||||
print(f"Kein zulaessiger Plan moeglich (Status: {status.value}).")
|
||||
print("Die harten Regeln widersprechen sich. Naechster Schritt: eine Regel")
|
||||
print("weich machen (Kapitel Praxisfallen) oder den Konflikt eingrenzen")
|
||||
print("(Anhang Fehlerdiagnose).")
|
||||
return None
|
||||
if not status.brauchbar:
|
||||
print(f"Solver lieferte kein Ergebnis (Status: {status.value}).")
|
||||
print("Zeitlimit erhoehen oder das Modell vereinfachen.")
|
||||
return None
|
||||
|
||||
guete = ("beweisbar optimal" if status is SolverStatus.OPTIMAL
|
||||
else "zulaessig, aber nicht bewiesen - Gap beachten")
|
||||
print(f"Solver-Status: {loeser.StatusName(rohstatus)} -> {status.value} ({guete})")
|
||||
print(f"{loesung.als_bericht()}\n")
|
||||
|
||||
# --- Plan ausgeben ----------------------------------------------------
|
||||
zuweisung = {s: next(p for p in PERSONAL if loeser.Value(x[p, s]) == 1)
|
||||
for s in range(len(SLOTS))}
|
||||
plan = [{
|
||||
"Stunde": f"Std. {SLOTS[s]}",
|
||||
"Klasse": KLASSEN[s],
|
||||
"Fach": FAECHER[s],
|
||||
"Vertretung": zuweisung[s],
|
||||
"Vorbelastung": f"{VORBELASTUNG[zuweisung[s]]} Std.",
|
||||
} for s in range(len(SLOTS))]
|
||||
print(pd.DataFrame(plan).to_string(index=False))
|
||||
|
||||
# --- Abnahmepruefung gegen die Anforderung ----------------------------
|
||||
beanstandungen = pruefe_plan(zuweisung)
|
||||
if beanstandungen:
|
||||
print("\nABNAHMEPRUEFUNG FEHLGESCHLAGEN:")
|
||||
for beanstandung in beanstandungen:
|
||||
print(f" - {beanstandung}")
|
||||
return None
|
||||
print("\nAbnahmepruefung: bestanden (Besetzung, Qualifikation, Hoechstzahl).")
|
||||
|
||||
# --- Erklaerbarkeit: Strafkosten aufschluesseln -----------------------
|
||||
kosten_last = sum(VORBELASTUNG[p] * STRAFE_VORBELASTUNG
|
||||
* sum(loeser.Value(x[p, s]) for s in range(len(SLOTS)))
|
||||
for p in PERSONAL)
|
||||
kosten_loecher = sum(loeser.Value(v) for v in loch_variablen) * STRAFE_LOCH
|
||||
kosten_fairness = loeser.Value(spannweite) * STRAFE_UNFAIRNESS
|
||||
|
||||
print("\n--- Woraus bestehen die Strafkosten? ---")
|
||||
print(f" Vorbelastung geschont: {kosten_last:5.0f} Punkte")
|
||||
print(f" Freistundenloecher: {kosten_loecher:5.0f} Punkte "
|
||||
f"({sum(loeser.Value(v) for v in loch_variablen)} Loch/Loecher)")
|
||||
print(f" Fairness (Spannweite {loeser.Value(spannweite)}): {kosten_fairness:5.0f} Punkte")
|
||||
print(f" {'Summe':<23} {loeser.ObjectiveValue():5.0f} Punkte")
|
||||
|
||||
print("\n--- Auslastung nach Optimierung ---")
|
||||
for p in PERSONAL:
|
||||
heute = sum(loeser.Value(x[p, s]) for s in range(len(SLOTS)))
|
||||
print(f" * {p:<15}: {VORBELASTUNG[p]} vorher + {heute} Vertretung "
|
||||
f"= {VORBELASTUNG[p] + heute} Stunden")
|
||||
|
||||
print(f"\nSuchstatistik: {loeser.NumBranches()} Verzweigungen, "
|
||||
f"{loeser.NumConflicts()} Konflikte, {loeser.WallTime():.3f} s")
|
||||
print("=" * 78)
|
||||
return plan
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
plane_vertretung()
|
||||
188
OR_HTML_04/programme/CVaR_Portfolio.py
Normal file
188
OR_HTML_04/programme/CVaR_Portfolio.py
Normal file
|
|
@ -0,0 +1,188 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# CVaR_Portfolio.py
|
||||
"""
|
||||
Kapitel CVaR: CVaR-Optimierung mit L1-Transaktionskosten via CVXPY.
|
||||
|
||||
Eigenschaften:
|
||||
* Spaltenreihenfolge erzwungen
|
||||
* Einheiten konsistent (alles taeglich, Annualisierung nur in der Ausgabe)
|
||||
* Szenario-Nebenbedingungen VEKTORISIERT statt in einer Python-Schleife
|
||||
(eine Matrixbedingung statt S einzelner Constraints - deutlich schneller)
|
||||
* Vergleich CVaR- gegen Varianz-Optimierung
|
||||
* Nachrechnung von VaR/CVaR aus den realisierten Szenarien
|
||||
"""
|
||||
|
||||
import time
|
||||
|
||||
import cvxpy as cp
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import yfinance as yf
|
||||
|
||||
TICKER = ["AAPL", "MSFT", "NVDA", "AMZN", "JNJ", "PFE", "JPM", "GS", "XOM", "CVX"]
|
||||
ALPHA = 0.95 # Konfidenzniveau: schlechteste 5 % der Tage
|
||||
MAX_GEWICHT = 0.25
|
||||
GEBUEHRENSATZ = 0.002 # 0,2 % je Einheit Turnover (Spread + Brokerage)
|
||||
RISIKOAVERSION = 1.5 # bezogen auf TAEGLICHE Groessen
|
||||
HANDELSTAGE = 252
|
||||
|
||||
|
||||
def lade_renditen():
|
||||
ende = pd.Timestamp.today().normalize()
|
||||
start = ende - pd.DateOffset(years=2)
|
||||
roh = yf.download(TICKER, start=start, end=ende, auto_adjust=True, progress=False)
|
||||
if roh.empty:
|
||||
raise SystemExit("Download fehlgeschlagen (Netz, Ticker oder Rate-Limit pruefen).")
|
||||
|
||||
if isinstance(roh.columns, pd.MultiIndex):
|
||||
kurse = roh["Close"][TICKER].dropna() # erzwingt eigene Spaltenreihenfolge
|
||||
else:
|
||||
kurse = roh[["Close"]].dropna()
|
||||
kurse.columns = TICKER
|
||||
assert list(kurse.columns) == TICKER, "Spaltenreihenfolge weicht ab!"
|
||||
return kurse.pct_change().dropna()
|
||||
|
||||
|
||||
def optimiere_cvar(R, w_alt, vektorisiert=True):
|
||||
"""
|
||||
Maximiere: taegliche Rendite - lambda * CVaR - Transaktionskosten
|
||||
Alle Groessen TAEGLICH.
|
||||
"""
|
||||
S, N = R.shape
|
||||
mu_taeglich = R.mean(axis=0)
|
||||
|
||||
w = cp.Variable(N, nonneg=True)
|
||||
gamma = cp.Variable() # wird im Optimum zum VaR
|
||||
u = cp.Variable(S, nonneg=True) # Ueberschuss ueber die Schwelle
|
||||
|
||||
cvar = gamma + (1.0 / (S * (1.0 - ALPHA))) * cp.sum(u)
|
||||
turnover = cp.norm1(w - w_alt)
|
||||
kosten = GEBUEHRENSATZ * turnover
|
||||
|
||||
ziel = cp.Maximize(mu_taeglich @ w - RISIKOAVERSION * cvar - kosten)
|
||||
|
||||
bedingungen = [cp.sum(w) == 1, w <= MAX_GEWICHT]
|
||||
if vektorisiert:
|
||||
# EINE Matrixbedingung statt S einzelner - deutlich schneller
|
||||
bedingungen.append(u >= -(R @ w) - gamma)
|
||||
else:
|
||||
for s in range(S): # Alternative: S einzelne Constraints (langsamer)
|
||||
bedingungen.append(u[s] >= -R[s] @ w - gamma)
|
||||
|
||||
problem = cp.Problem(ziel, bedingungen)
|
||||
problem.solve()
|
||||
if problem.status not in ("optimal", "optimal_inaccurate"):
|
||||
raise SystemExit(f"CVaR-Problem nicht loesbar: {problem.status}")
|
||||
return w.value, float(gamma.value), float(cvar.value), float(turnover.value)
|
||||
|
||||
|
||||
def optimiere_varianz(R, w_alt):
|
||||
"""Klassisches Mean-Variance zum Vergleich - ebenfalls taeglich gerechnet.
|
||||
|
||||
ACHTUNG, DCP-Falle: Die Standardabweichung ist hier NICHT als
|
||||
cp.sqrt(cp.quad_form(w, sigma)) formulierbar. cp.sqrt ist konkav und
|
||||
verlangt ein konkaves Argument; quad_form ist konvex - CVXPY lehnt den
|
||||
Ausdruck mit einem DCPError ab, und zwar voellig zu Recht (Anhang
|
||||
Fehlerdiagnose). cp.psd_wrap() hilft dagegen nicht: Es behebt eine
|
||||
NUMERISCHE Beanstandung an sigma, keine Regelverletzung im Aufbau.
|
||||
|
||||
Der Standardweg ist die Cholesky-Zerlegung sigma = L L^T. Damit gilt
|
||||
w' sigma w = ||L^T w||^2, also ist die Standardabweichung die 2-Norm
|
||||
eines AFFINEN Ausdrucks - konvex und damit regelkonform.
|
||||
"""
|
||||
N = R.shape[1]
|
||||
mu_taeglich = R.mean(axis=0)
|
||||
sigma = np.cov(R, rowvar=False, ddof=1)
|
||||
# Der winzige Diagonalzuschlag faengt den Fall ab, dass sigma numerisch
|
||||
# nur halbdefinit ist (mehr Titel als Handelstage, doppelte Spalten).
|
||||
L = np.linalg.cholesky(sigma + 1e-12 * np.eye(N))
|
||||
|
||||
w = cp.Variable(N, nonneg=True)
|
||||
ziel = cp.Maximize(mu_taeglich @ w
|
||||
- RISIKOAVERSION * cp.norm2(L.T @ w)
|
||||
- GEBUEHRENSATZ * cp.norm1(w - w_alt))
|
||||
problem = cp.Problem(ziel, [cp.sum(w) == 1, w <= MAX_GEWICHT])
|
||||
problem.solve()
|
||||
if problem.status not in ("optimal", "optimal_inaccurate"):
|
||||
raise SystemExit(f"Varianz-Problem nicht loesbar: {problem.status}")
|
||||
return w.value
|
||||
|
||||
|
||||
def realisierte_kennzahlen(w, R):
|
||||
"""VaR und CVaR direkt aus den Szenarien - unabhaengige Gegenprobe."""
|
||||
portfoliorenditen = R @ w
|
||||
verluste = -portfoliorenditen
|
||||
var = float(np.quantile(verluste, ALPHA))
|
||||
cvar = float(verluste[verluste >= var].mean())
|
||||
return var, cvar, float(portfoliorenditen.mean()), float(portfoliorenditen.std(ddof=1))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
renditen = lade_renditen()
|
||||
R = renditen.values
|
||||
S, N = R.shape
|
||||
w_alt = np.ones(N) / N # Ausgangslage: Gleichgewichtung
|
||||
|
||||
t0 = time.perf_counter()
|
||||
w_cvar, var_modell, cvar_modell, turnover = optimiere_cvar(R, w_alt, vektorisiert=True)
|
||||
dauer_vektor = time.perf_counter() - t0
|
||||
|
||||
w_var = optimiere_varianz(R, w_alt)
|
||||
|
||||
print("=" * 90)
|
||||
print(" CVaR-PORTFOLIO-OPTIMIERUNG MIT TRANSAKTIONSKOSTEN")
|
||||
print("=" * 90)
|
||||
print(f"Datenbasis: {S} Handelstage, {N} Titel | Konfidenzniveau "
|
||||
f"{ALPHA*100:.0f} % | Loesungszeit {dauer_vektor:.2f} s\n")
|
||||
|
||||
# --- Gegenprobe: Modellwerte gegen realisierte Szenariowerte ---------
|
||||
var_real, cvar_real, mu_real, sd_real = realisierte_kennzahlen(w_cvar, R)
|
||||
print("--- Gegenprobe: stimmen Modell und Szenarien ueberein? ---")
|
||||
print(f" VaR aus dem Modell (gamma): {var_modell*100:7.4f} % | "
|
||||
f"aus den Szenarien: {var_real*100:7.4f} %")
|
||||
print(f" CVaR aus dem Modell: {cvar_modell*100:7.4f} % | "
|
||||
f"aus den Szenarien: {cvar_real*100:7.4f} %")
|
||||
assert abs(cvar_modell - cvar_real) < 1e-4, "CVaR stimmt nicht mit den Szenarien!"
|
||||
print(" -> Der Rockafellar-Uryasev-Trick liefert exakt den empirischen CVaR.")
|
||||
|
||||
# --- Kennzahlen beider Portfolios ------------------------------------
|
||||
print(f"\n{'Portfolio':<24} {'Rendite p.a.':>13} {'Vola p.a.':>11} "
|
||||
f"{'VaR 95% (Tag)':>15} {'CVaR 95% (Tag)':>16} {'Turnover':>10}")
|
||||
print("-" * 90)
|
||||
for name, w in [("CVaR-optimiert", w_cvar), ("Varianz-optimiert", w_var),
|
||||
("Gleichgewichtung", w_alt)]:
|
||||
v, c, m, s = realisierte_kennzahlen(w, R)
|
||||
to = float(np.abs(w - w_alt).sum())
|
||||
print(f"{name:<24} {m*HANDELSTAGE*100:>12.2f} % "
|
||||
f"{s*np.sqrt(HANDELSTAGE)*100:>10.2f} % {v*100:>14.3f} % "
|
||||
f"{c*100:>15.3f} % {to*100:>9.1f} %")
|
||||
|
||||
print("\nHinweis zur Annualisierung: Renditen werden mit 252 skaliert,")
|
||||
print("Volatilitaeten mit sqrt(252). Fuer VaR/CVaR ist eine solche Skalierung")
|
||||
print("nur unter starken Annahmen (Unabhaengigkeit, kein Drift) zulaessig -")
|
||||
print("sie werden hier deshalb bewusst als TAGESwerte ausgewiesen.")
|
||||
|
||||
# --- Allokationstabelle ------------------------------------------------
|
||||
print("\n--- Allokation ---")
|
||||
tabelle = pd.DataFrame({
|
||||
"Ticker": TICKER,
|
||||
"vorher": [f"{v*100:5.1f} %" for v in w_alt],
|
||||
"CVaR-opt.": [f"{v*100:5.1f} %" for v in w_cvar],
|
||||
"Handel": [f"{(w_cvar[i]-w_alt[i])*100:+6.1f} %" for i in range(N)],
|
||||
"Varianz-opt.": [f"{v*100:5.1f} %" for v in w_var],
|
||||
})
|
||||
print(tabelle.to_string(index=False))
|
||||
print(f"\nTurnover {turnover*100:.1f} % -> Transaktionskosten "
|
||||
f"{GEBUEHRENSATZ*turnover*100:.3f} % des Portfoliowerts")
|
||||
|
||||
# --- Laufzeitvergleich vektorisiert vs. Schleife ----------------------
|
||||
if S <= 600: # bei sehr vielen Szenarien zu langsam
|
||||
t0 = time.perf_counter()
|
||||
optimiere_cvar(R, w_alt, vektorisiert=False)
|
||||
dauer_schleife = time.perf_counter() - t0
|
||||
print(f"\n--- Laufzeit: {S} Nebenbedingungen aufbauen ---")
|
||||
print(f" vektorisiert (u >= -(R @ w) - gamma): {dauer_vektor:6.2f} s")
|
||||
print(f" Schleife ueber Szenarien (V01): {dauer_schleife:6.2f} s "
|
||||
f"({dauer_schleife/dauer_vektor:.1f}x langsamer)")
|
||||
print("=" * 90)
|
||||
133
OR_HTML_04/programme/Data_Snooping.py
Normal file
133
OR_HTML_04/programme/Data_Snooping.py
Normal file
|
|
@ -0,0 +1,133 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# Data_Snooping.py
|
||||
"""
|
||||
Kapitel Handelsmaschine: Die Selbsttaeuschung, die man einem Backtest nicht ansieht.
|
||||
|
||||
Backtest_Fallen.py zeigt drei Fehler, die IM Backtest stecken: Lookahead,
|
||||
Survivorship, vergessene Kosten. Alle drei kann man einem Programm ansehen,
|
||||
wenn man es liest.
|
||||
|
||||
Die beiden anderen - Overfitting und Data Snooping - entstehen nicht im
|
||||
Backtest, sondern in der ARBEITSWEISE davor. Sie hinterlassen im Code keine
|
||||
Spur. Ein Backtest, der zwanzigmal angepasst wurde, sieht am Ende genauso
|
||||
sauber aus wie einer, der beim ersten Versuch funktioniert hat.
|
||||
|
||||
Dieses Programm macht den Effekt messbar. Es testet Strategien, die
|
||||
NACHWEISLICH keinerlei Prognosekraft haben (ihre taeglichen Ergebnisse sind
|
||||
Zufallszahlen mit Erwartungswert exakt null), und zeigt, wie gut die jeweils
|
||||
BESTE davon aussieht - allein als Funktion der Anzahl der Versuche.
|
||||
|
||||
Benoetigt: numpy, scipy
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import numpy as np
|
||||
from scipy import stats
|
||||
|
||||
HANDELSTAGE = 1260 # fuenf Jahre
|
||||
TAGESSCHWANKUNG = 0.010 # 1 % taeglich
|
||||
WIEDERHOLUNGEN = 300 # damit die Tabelle nicht vom Zufall abhaengt
|
||||
NIVEAU = 0.05
|
||||
|
||||
RNG = np.random.default_rng(42)
|
||||
|
||||
|
||||
def bewerte(ergebnisse: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
|
||||
"""Sharpe-Kennzahl und einseitiger p-Wert je Strategie.
|
||||
|
||||
ergebnisse hat die Form (Strategien, Handelstage). Der p-Wert prueft die
|
||||
Nullhypothese 'mittleres Tagesergebnis ist null' - also genau die Frage,
|
||||
die ein Backtest beantworten soll.
|
||||
"""
|
||||
mittel = ergebnisse.mean(axis=1)
|
||||
streuung = ergebnisse.std(axis=1, ddof=1)
|
||||
sharpe = mittel / streuung * np.sqrt(252)
|
||||
t_wert = mittel / (streuung / np.sqrt(ergebnisse.shape[1]))
|
||||
return sharpe, stats.t.sf(t_wert, ergebnisse.shape[1] - 1)
|
||||
|
||||
|
||||
def probiere(anzahl_strategien: int) -> tuple[float, float, float, float]:
|
||||
"""Liefert (beste Sharpe, ihr p-Wert, Anteil scheinbar signifikanter,
|
||||
Anteil, der auch die Bonferroni-Huerde nimmt) - gemittelt ueber
|
||||
WIEDERHOLUNGEN unabhaengige Durchgaenge."""
|
||||
beste_sharpe, beste_p, anteil_signifikant, ueberlebt_bonferroni = [], [], [], []
|
||||
|
||||
for _ in range(WIEDERHOLUNGEN):
|
||||
# Der springende Punkt: Erwartungswert exakt 0. Keine dieser
|
||||
# Strategien hat einen echten Vorteil - per Konstruktion.
|
||||
ergebnisse = RNG.normal(0.0, TAGESSCHWANKUNG,
|
||||
(anzahl_strategien, HANDELSTAGE))
|
||||
sharpe, p_werte = bewerte(ergebnisse)
|
||||
|
||||
beste = int(sharpe.argmax())
|
||||
beste_sharpe.append(float(sharpe[beste]))
|
||||
beste_p.append(float(p_werte[beste]))
|
||||
anteil_signifikant.append(float((p_werte < NIVEAU).mean()))
|
||||
# Bonferroni: Wer n Tests macht, muss die Huerde durch n teilen.
|
||||
ueberlebt_bonferroni.append(
|
||||
float((p_werte < NIVEAU / anzahl_strategien).mean()))
|
||||
|
||||
return (float(np.mean(beste_sharpe)), float(np.mean(beste_p)),
|
||||
float(np.mean(anteil_signifikant)),
|
||||
float(np.mean(ueberlebt_bonferroni)))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
print("=" * 78)
|
||||
print(" WIE GUT SIEHT DIE BESTE VON N STRATEGIEN AUS,")
|
||||
print(" WENN KEINE EINZIGE ETWAS TAUGT?")
|
||||
print("=" * 78)
|
||||
print(f"{HANDELSTAGE} Handelstage, taegliche Ergebnisse mit Erwartungswert "
|
||||
f"EXAKT null.")
|
||||
print(f"Jede Zeile ist der Mittelwert aus {WIEDERHOLUNGEN} unabhaengigen "
|
||||
f"Durchgaengen.\n")
|
||||
|
||||
print(f"{'getestet':>9} {'beste Sharpe':>14} {'ihr p-Wert':>12} "
|
||||
f"{'scheinbar sign.':>16}")
|
||||
print("-" * 78)
|
||||
ergebnisse = {}
|
||||
for anzahl in (1, 10, 100, 1000):
|
||||
sharpe, p_wert, anteil, _ = probiere(anzahl)
|
||||
ergebnisse[anzahl] = (sharpe, p_wert, anteil)
|
||||
print(f"{anzahl:9d} {sharpe:14.2f} {p_wert:12.4f} {anteil * 100:15.1f} %")
|
||||
print("-" * 78)
|
||||
|
||||
sharpe_100, p_100, _ = ergebnisse[100]
|
||||
sharpe_1000, p_1000, _ = ergebnisse[1000]
|
||||
|
||||
print("\nLesen Sie die dritte Spalte von unten nach oben:")
|
||||
print(f" Wer {1000} Varianten durchprobiert, findet eine mit p = {p_1000:.4f} -")
|
||||
print(" ein Wert, den jede Zeitschrift als 'hochsignifikant' druckt.")
|
||||
print(f" Wer {100} probiert, findet p = {p_100:.4f} und eine Sharpe-Kennzahl")
|
||||
print(f" von {sharpe_100:.2f} - damit geht man zum Vorgesetzten.")
|
||||
print(" Und keine dieser Strategien hat auch nur den Hauch eines Vorteils.")
|
||||
print()
|
||||
print("Die LETZTE Spalte bleibt dagegen konstant bei rund 5 %. Das ist kein")
|
||||
print("Widerspruch, sondern die Definition: Ein Test zum 5-%-Niveau irrt in")
|
||||
print("5 % der Faelle. Was sich aendert, ist nicht der Anteil, sondern die")
|
||||
print("ANZAHL - und weil man nur die beste Strategie praesentiert, sieht man")
|
||||
print("von allen Fehlversuchen genau einen: den erfolgreichsten.")
|
||||
|
||||
print("\n" + "=" * 78)
|
||||
print(" WAS MAN DAGEGEN TUN KANN")
|
||||
print("=" * 78)
|
||||
print("1. ZAEHLEN. Notieren Sie, wie viele Varianten Sie ausprobiert haben -")
|
||||
print(" auch die verworfenen, auch die 'nur mal schnell'. Ohne diese Zahl")
|
||||
print(" ist kein p-Wert und keine Sharpe-Kennzahl interpretierbar.")
|
||||
print()
|
||||
print("2. KORRIGIEREN. Die Bonferroni-Korrektur teilt die Huerde durch die")
|
||||
print(" Zahl der Versuche:")
|
||||
for anzahl in (10, 100, 1000):
|
||||
print(f" bei {anzahl:5d} Versuchen: p < {NIVEAU / anzahl:.5f} statt "
|
||||
f"p < {NIVEAU:.2f}")
|
||||
print(" Sie ist konservativ, aber sie ist ehrlich. Fuer Handelsstrategien")
|
||||
print(" gibt es verfeinerte Varianten (Deflated Sharpe Ratio nach Bailey")
|
||||
print(" und Lopez de Prado), die dasselbe Prinzip verfolgen.")
|
||||
print()
|
||||
print("3. ZURUECKHALTEN. Legen Sie einen Zeitraum beiseite, BEVOR Sie")
|
||||
print(" anfangen, und ruehren Sie ihn bis zum Schluss nicht an. Er ist")
|
||||
print(" der einzige Test, den Sie nicht durch Probieren verderben koennen -")
|
||||
print(" und er ist genau einmal verwendbar.")
|
||||
print("=" * 78)
|
||||
47
OR_HTML_04/programme/Diversifikation_Demo.py
Normal file
47
OR_HTML_04/programme/Diversifikation_Demo.py
Normal file
|
|
@ -0,0 +1,47 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# Diversifikation_Demo.py
|
||||
"""
|
||||
Kapitel Markowitz: Portfoliorisiko in Abhaengigkeit von der Korrelation - die
|
||||
Handrechnung zum Gratis-Mittagessen, geprueft ueber die analytische Formel und
|
||||
ueber simulierte, korrelierte Renditen.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
|
||||
MU = 0.08
|
||||
SIGMA = 0.20
|
||||
|
||||
|
||||
def sigma_portfolio(rho: float, w1: float = 0.5) -> float:
|
||||
"""Analytische Formel: sigma_p^2 = w1^2 s1^2 + w2^2 s2^2 + 2 w1 w2 rho s1 s2."""
|
||||
w2 = 1 - w1
|
||||
varianz = w1**2 * SIGMA**2 + w2**2 * SIGMA**2 + 2 * w1 * w2 * rho * SIGMA**2
|
||||
return float(np.sqrt(varianz))
|
||||
|
||||
|
||||
def simuliere(rho: float, n: int = 500_000, seed: int = 7) -> float:
|
||||
"""Erzeugt n korrelierte Renditepaare und misst die Portfolio-Volatilitaet empirisch."""
|
||||
rng = np.random.default_rng(seed)
|
||||
kovarianz = np.array([[SIGMA**2, rho * SIGMA**2],
|
||||
[rho * SIGMA**2, SIGMA**2]])
|
||||
renditen = rng.multivariate_normal([MU, MU], kovarianz, size=n)
|
||||
portfolio = 0.5 * renditen[:, 0] + 0.5 * renditen[:, 1]
|
||||
return float(portfolio.std(ddof=1))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
print("=" * 78)
|
||||
print(" PORTFOLIORISIKO IN ABHAENGIGKEIT VON DER KORRELATION")
|
||||
print("=" * 78)
|
||||
print(f"{'rho':>6} | {'sigma_p (Formel)':>18} | {'sigma_p (Simulation)':>20} | {'Risikoreduktion':>16}")
|
||||
print("-" * 78)
|
||||
for rho in [1.0, 0.5, 0.0, -0.5, -1.0]:
|
||||
formel = sigma_portfolio(rho)
|
||||
sim = simuliere(rho)
|
||||
reduktion = (1 - formel / SIGMA) * 100
|
||||
print(f"{rho:6.1f} | {formel*100:16.2f} % | {sim*100:18.2f} % | {reduktion:14.0f} %")
|
||||
|
||||
print("\nDie simulierten Werte (500.000 gezogene Renditepaare je rho) bestaetigen")
|
||||
print("die Formel aus der Handrechnung bis auf statistisches Rauschen - und das,")
|
||||
print("obwohl Formel und Simulation nichts voneinander wissen.")
|
||||
78
OR_HTML_04/programme/Dualitaet_Nachweis.py
Normal file
78
OR_HTML_04/programme/Dualitaet_Nachweis.py
Normal file
|
|
@ -0,0 +1,78 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# Dualitaet_Nachweis.py
|
||||
"""
|
||||
Kapitel LP: Primales und duales Problem unabhaengig loesen und den starken
|
||||
Dualitaetssatz sowie den komplementaeren Schlupf numerisch nachweisen.
|
||||
|
||||
Modell aus der Simplex-Handrechnung (Bot-Allokation, LP-Relaxation):
|
||||
max 150*x1 + 250*x2 u.d.N. 2*x1+5*x2<=40, 4*x1+6*x2<=60, x1<=8, x>=0
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
from scipy.optimize import linprog
|
||||
|
||||
# --- Primales Problem --------------------------------------------------------
|
||||
C_PRIMAL = np.array([150.0, 250.0])
|
||||
A_PRIMAL = np.array([[2.0, 5.0], [4.0, 6.0], [1.0, 0.0]])
|
||||
B_PRIMAL = np.array([40.0, 60.0, 8.0])
|
||||
|
||||
# --- Duales Problem: min b^T y u.d.N. A^T y >= c, y >= 0 ------------------
|
||||
# linprog kennt nur <=, also A^T y >= c <=> -A^T y <= -c
|
||||
C_DUAL = B_PRIMAL
|
||||
A_DUAL = -A_PRIMAL.T
|
||||
B_DUAL = -C_PRIMAL
|
||||
|
||||
|
||||
def loese():
|
||||
primal = linprog(c=-C_PRIMAL, A_ub=A_PRIMAL, b_ub=B_PRIMAL,
|
||||
bounds=[(0, None)] * 2, method="highs")
|
||||
dual = linprog(c=C_DUAL, A_ub=A_DUAL, b_ub=B_DUAL,
|
||||
bounds=[(0, None)] * 3, method="highs")
|
||||
if not primal.success or not dual.success:
|
||||
raise SystemExit("Primal oder Dual nicht loesbar.")
|
||||
return primal, dual
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
primal, dual = loese()
|
||||
x = primal.x
|
||||
y = dual.x
|
||||
z_primal = -primal.fun
|
||||
z_dual = dual.fun
|
||||
|
||||
print("=" * 78)
|
||||
print(" PRIMALES PROBLEM")
|
||||
print("=" * 78)
|
||||
print(f"x* = ({x[0]:.4f}, {x[1]:.4f})")
|
||||
print(f"Z* = {z_primal:.4f}")
|
||||
|
||||
print("\n" + "=" * 78)
|
||||
print(" DUALES PROBLEM")
|
||||
print("=" * 78)
|
||||
print(f"y* = ({y[0]:.4f}, {y[1]:.4f}, {y[2]:.4f})")
|
||||
print(f"W* = {z_dual:.4f}")
|
||||
|
||||
print("\n" + "=" * 78)
|
||||
print(" STARKER DUALITAETSSATZ: c^T x* == b^T y* ?")
|
||||
print("=" * 78)
|
||||
print(f" Primal Z* = {z_primal:.6f}")
|
||||
print(f" Dual W* = {z_dual:.6f}")
|
||||
differenz = abs(z_primal - z_dual)
|
||||
print(f" Differenz = {differenz:.2e} -> "
|
||||
f"{'BESTAETIGT' if differenz < 1e-6 else 'VERLETZT!'}")
|
||||
|
||||
print("\n" + "=" * 78)
|
||||
print(" KOMPLEMENTAERER SCHLUPF: s_i * y_i == 0 fuer alle i ?")
|
||||
print("=" * 78)
|
||||
schlupf = B_PRIMAL - A_PRIMAL @ x
|
||||
ressourcen = ["vCPU (s1)", "RAM (s2)", "Marktlimit (s3)"]
|
||||
for name, s, yi in zip(ressourcen, schlupf, y):
|
||||
produkt = s * yi
|
||||
print(f" {name:<16} Schlupf s={s:6.4f} Schattenpreis y={yi:6.4f} "
|
||||
f"s*y={produkt:.2e} {'OK' if abs(produkt) < 1e-6 else 'VERLETZT!'}")
|
||||
|
||||
print("\nFazit: Das dual geloeste y* stimmt exakt mit den Schattenpreisen")
|
||||
print("überein, die die Simplex-Rechnung von Hand in der Z-Zeile")
|
||||
print("ablas - unabhaengig voneinander berechnet, identisches Ergebnis.")
|
||||
print("=" * 78)
|
||||
120
OR_HTML_04/programme/Ein_System_Vier_Ansaetze.py
Normal file
120
OR_HTML_04/programme/Ein_System_Vier_Ansaetze.py
Normal file
|
|
@ -0,0 +1,120 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# Ein_System_Vier_Ansaetze.py
|
||||
"""
|
||||
Kapitel Oekosystem: Dasselbe LP in vier Bibliotheken.
|
||||
max 10*x1 + 15*x2 + 25*x3
|
||||
u.d.N. x1 + x2 + 2*x3 <= 40
|
||||
2*x1 + 3*x2 + x3 <= 50
|
||||
x >= 0
|
||||
|
||||
Deckt scipy.optimize, highspy, CVXPY und OR-Tools/GLOP ab, mit Kreuzvergleich am Ende.
|
||||
|
||||
WICHTIG: Jeder Solver läuft in einem EIGENEN Prozess, weil sich ortools und
|
||||
highspy auf vielen Systemen nicht gemeinsam importieren lassen (beide bringen
|
||||
eine eigene HiGHS-Kopie mit -> Symbolkonflikt).
|
||||
"""
|
||||
|
||||
import json
|
||||
import subprocess
|
||||
import sys
|
||||
import textwrap
|
||||
import time
|
||||
|
||||
ERWARTET = 530.0 # Ergebnis der Handrechnung zum Produktionsprogramm
|
||||
|
||||
# Jeder Eintrag ist ein eigenständiges Miniprogramm, das sein Ergebnis als
|
||||
# JSON auf stdout ausgibt. So bleibt jeder Import in seinem eigenen Prozess.
|
||||
ANSAETZE: dict[str, str] = {
|
||||
|
||||
"scipy.optimize.linprog": """
|
||||
from scipy.optimize import linprog
|
||||
res = linprog(c=[-10.0, -15.0, -25.0], # linprog MINIMIERT -> negieren
|
||||
A_ub=[[1, 1, 2], [2, 3, 1]], b_ub=[40, 50],
|
||||
bounds=[(0, None)] * 3, method="highs")
|
||||
ausgabe = (-res.fun, list(res.x))
|
||||
""",
|
||||
|
||||
"highspy (natives HiGHS)": """
|
||||
import numpy as np, highspy
|
||||
h = highspy.Highs()
|
||||
h.setOptionValue("output_flag", False)
|
||||
h.addVars(3, np.zeros(3), np.full(3, highspy.kHighsInf))
|
||||
h.changeObjectiveSense(highspy.ObjSense.kMaximize)
|
||||
for j, wert in enumerate([10.0, 15.0, 25.0]):
|
||||
h.changeColCost(j, wert)
|
||||
# CSR-Format: starts[i] = Beginn von Zeile i in indices/values
|
||||
h.addRows(2, np.full(2, -highspy.kHighsInf), np.array([40.0, 50.0]), 6,
|
||||
np.array([0, 3], dtype=np.int32),
|
||||
np.array([0, 1, 2, 0, 1, 2], dtype=np.int32),
|
||||
np.array([1.0, 1.0, 2.0, 2.0, 3.0, 1.0]))
|
||||
h.run()
|
||||
ausgabe = (h.getInfo().objective_function_value,
|
||||
list(h.getSolution().col_value[:3]))
|
||||
""",
|
||||
|
||||
"cvxpy": """
|
||||
import numpy as np, cvxpy as cp
|
||||
x = cp.Variable(3, nonneg=True)
|
||||
problem = cp.Problem(cp.Maximize(np.array([10.0, 15.0, 25.0]) @ x),
|
||||
[np.array([[1, 1, 2], [2, 3, 1]]) @ x <= np.array([40, 50])])
|
||||
problem.solve()
|
||||
ausgabe = (float(problem.value), [float(v) for v in x.value])
|
||||
""",
|
||||
|
||||
"ortools / GLOP": """
|
||||
from ortools.linear_solver import pywraplp
|
||||
s = pywraplp.Solver.CreateSolver("GLOP")
|
||||
x = [s.NumVar(0, s.infinity(), f"x{j+1}") for j in range(3)]
|
||||
A = [[1, 1, 2], [2, 3, 1]]
|
||||
for i, kap in enumerate([40, 50]):
|
||||
s.Add(sum(A[i][j] * x[j] for j in range(3)) <= kap)
|
||||
s.Maximize(10 * x[0] + 15 * x[1] + 25 * x[2])
|
||||
s.Solve()
|
||||
ausgabe = (s.Objective().Value(), [v.solution_value() for v in x])
|
||||
""",
|
||||
}
|
||||
|
||||
|
||||
def fuehre_in_eigenem_prozess_aus(quelltext: str) -> tuple[float, list[float]]:
|
||||
"""Startet den Codeschnipsel als separaten Python-Prozess und liest das Ergebnis."""
|
||||
programm = textwrap.dedent(quelltext) + "\nimport json; print(json.dumps(ausgabe))\n"
|
||||
ergebnis = subprocess.run([sys.executable, "-c", programm],
|
||||
capture_output=True, text=True, timeout=120)
|
||||
if ergebnis.returncode != 0:
|
||||
raise RuntimeError(ergebnis.stderr.strip().splitlines()[-1])
|
||||
wert, loesung = json.loads(ergebnis.stdout.strip().splitlines()[-1])
|
||||
return wert, loesung
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
print("=" * 78)
|
||||
print(" EIN SYSTEM - VIER ANSAETZE (je eigener Prozess)")
|
||||
print("=" * 78)
|
||||
print(f"{'Bibliothek':<26} {'Z*':>10} {'x1':>7} {'x2':>7} {'x3':>7} {'Zeit':>10}")
|
||||
print("-" * 78)
|
||||
|
||||
werte = []
|
||||
for name, quelltext in ANSAETZE.items():
|
||||
t0 = time.perf_counter()
|
||||
try:
|
||||
wert, x = fuehre_in_eigenem_prozess_aus(quelltext)
|
||||
except RuntimeError as fehler:
|
||||
print(f"{name:<26} nicht verfuegbar: {fehler[:40]}")
|
||||
continue
|
||||
dauer = time.perf_counter() - t0
|
||||
werte.append(wert)
|
||||
print(f"{name:<26} {wert:>10.2f} {x[0]:>7.2f} {x[1]:>7.2f} {x[2]:>7.2f} "
|
||||
f"{dauer:>8.2f} s")
|
||||
|
||||
print("-" * 78)
|
||||
spanne = max(werte) - min(werte)
|
||||
print(f"Spannweite zwischen den Bibliotheken: {spanne:.2e}")
|
||||
print(f"Abweichung zur Handrechnung ({ERWARTET:.0f}): "
|
||||
f"{abs(werte[0] - ERWARTET):.2e}")
|
||||
assert spanne < 1e-6, "Die Bibliotheken widersprechen sich!"
|
||||
assert abs(werte[0] - ERWARTET) < 1e-6, "Ergebnis weicht von der Handrechnung ab!"
|
||||
print("Alle Wege fuehren zum selben, von Hand bestaetigten Optimum.")
|
||||
print("(Die Zeiten enthalten den Prozessstart und den Import - sie messen")
|
||||
print(" NICHT die reine Solverleistung, siehe Uebung 3.5.)")
|
||||
print("=" * 78)
|
||||
138
OR_HTML_04/programme/Entropie_Maximierte_Allokation.py
Normal file
138
OR_HTML_04/programme/Entropie_Maximierte_Allokation.py
Normal file
|
|
@ -0,0 +1,138 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# Entropie_Maximierte_Allokation.py
|
||||
"""
|
||||
Kapitel QP/NLP: Nichtlineare Optimierung (NLP) mit scipy SLSQP.
|
||||
Modell: Mean-Variance-Portfolio mit Shannon-Entropie-Diversifikation.
|
||||
|
||||
Achtung, haeufiger Fehler: Sigma per A@A.T zu erzeugen und anschliessend die
|
||||
Diagonale zu ueberschreiben zerstoert die positive Semidefinitheit - die
|
||||
Matrix kann dadurch einen negativen Eigenwert bekommen und unmoegliche
|
||||
Korrelationen (> 1) implizieren. Dieses Programm baut Sigma stattdessen als
|
||||
D * C * D aus Volatilitaeten und einer echten Korrelationsmatrix -
|
||||
konstruktionsbedingt immer PSD.
|
||||
|
||||
Zusaetzlich: Multistart, weil SLSQP nur lokale Optima findet.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
from scipy.optimize import minimize
|
||||
|
||||
ASSETS = ["Tech-Aktien", "Rohstoffe", "US-Treasuries", "Krypto"]
|
||||
MU = np.array([0.14, 0.07, 0.03, 0.22]) # erwartete Jahresrenditen
|
||||
VOLA = np.array([0.20, 0.14, 0.07, 0.40]) # Volatilitaeten p.a.
|
||||
KORRELATION = np.array([
|
||||
[1.00, 0.25, -0.10, 0.55],
|
||||
[0.25, 1.00, 0.05, 0.20],
|
||||
[-0.10, 0.05, 1.00, -0.15],
|
||||
[0.55, 0.20, -0.15, 1.00],
|
||||
])
|
||||
|
||||
ALPHA = 1.0 # Gewicht des Ertrags
|
||||
BETA = 0.015 # Gewicht der Entropie-Diversifikation
|
||||
UNTERGRENZE = 0.001
|
||||
OBERGRENZE = 0.80
|
||||
|
||||
|
||||
def baue_kovarianz() -> np.ndarray:
|
||||
"""Sigma = D * C * D. Immer PSD, wenn C eine gueltige Korrelationsmatrix ist."""
|
||||
eig_c = np.linalg.eigvalsh(KORRELATION)
|
||||
if eig_c.min() < -1e-10:
|
||||
raise ValueError(f"Korrelationsmatrix unmoeglich (Eigenwert {eig_c.min():.4f}).")
|
||||
D = np.diag(VOLA)
|
||||
Sigma = D @ KORRELATION @ D
|
||||
eig_s = np.linalg.eigvalsh(Sigma)
|
||||
assert eig_s.min() > 0, "Sigma nicht positiv definit!"
|
||||
return Sigma
|
||||
|
||||
|
||||
SIGMA = baue_kovarianz()
|
||||
|
||||
|
||||
def zielfunktion(w, alpha=ALPHA, beta=BETA):
|
||||
"""Risiko - Ertrag - Entropiepraemie (wird minimiert)."""
|
||||
varianz = 0.5 * w @ SIGMA @ w
|
||||
ertrag = MU @ w
|
||||
entropie = -np.sum(w * np.log(w + 1e-12))
|
||||
return varianz - alpha * ertrag - beta * entropie
|
||||
|
||||
|
||||
def gradient(w, alpha=ALPHA, beta=BETA):
|
||||
"""Exakter analytischer Gradient - beschleunigt und stabilisiert SLSQP."""
|
||||
grad_varianz = SIGMA @ w # d/dw von 0.5 w'Sigma w
|
||||
grad_ertrag = -alpha * MU
|
||||
grad_entropie = beta * (np.log(w + 1e-12) + 1.0)
|
||||
return grad_varianz + grad_ertrag + grad_entropie
|
||||
|
||||
|
||||
def optimiere(startpunkt, alpha=ALPHA, beta=BETA):
|
||||
"""alpha und beta werden durchgereicht - so bleibt die Funktion seiteneffektfrei."""
|
||||
return minimize(
|
||||
zielfunktion, startpunkt, args=(alpha, beta), jac=gradient,
|
||||
method="SLSQP", bounds=[(UNTERGRENZE, OBERGRENZE)] * len(MU),
|
||||
constraints=({"type": "eq",
|
||||
"fun": lambda w: np.sum(w) - 1.0,
|
||||
"jac": lambda w: np.ones(len(MU))}),
|
||||
options={"ftol": 1e-12, "maxiter": 300})
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
n = len(MU)
|
||||
eig = np.linalg.eigvalsh(SIGMA)
|
||||
|
||||
print("=" * 74)
|
||||
print(" NICHTLINEARE ENTROPIE-OPTIMIERTE ASSET-ALLOKATION")
|
||||
print("=" * 74)
|
||||
print("Pruefung der Kovarianzmatrix:")
|
||||
print(f" Eigenwerte: {np.round(eig, 6)}")
|
||||
print(f" Positiv definit: {bool(eig.min() > 0)}")
|
||||
print(f" Groesste Korrelation ausserhalb der Diagonale: "
|
||||
f"{np.abs(KORRELATION - np.eye(n)).max():.2f} (muss <= 1 sein)")
|
||||
|
||||
# --- Multistart: SLSQP findet nur lokale Optima -----------------------
|
||||
rng = np.random.default_rng(7)
|
||||
startpunkte = [np.ones(n) / n] # Gleichgewichtung
|
||||
for _ in range(9):
|
||||
z = rng.random(n) + 0.05
|
||||
startpunkte.append(z / z.sum())
|
||||
|
||||
ergebnisse = [optimiere(s) for s in startpunkte]
|
||||
erfolgreich = [r for r in ergebnisse if r.success]
|
||||
bestes = min(erfolgreich, key=lambda r: r.fun)
|
||||
zielwerte = np.array([r.fun for r in erfolgreich])
|
||||
|
||||
print(f"\nMultistart mit {len(startpunkte)} Startpunkten:")
|
||||
print(f" Erfolgreich konvergiert: {len(erfolgreich)}")
|
||||
print(f" Spannweite der Zielwerte: {zielwerte.max() - zielwerte.min():.2e}")
|
||||
print(" -> " + ("alle Startpunkte fuehren zum selben Optimum (Indiz fuer Konvexitaet)"
|
||||
if zielwerte.max() - zielwerte.min() < 1e-6
|
||||
else "ACHTUNG: verschiedene lokale Optima gefunden!"))
|
||||
|
||||
w_opt = bestes.x
|
||||
rendite = MU @ w_opt
|
||||
volatilitaet = np.sqrt(w_opt @ SIGMA @ w_opt)
|
||||
entropie = -np.sum(w_opt * np.log(w_opt))
|
||||
|
||||
print(f"\nKonvergenz: {bestes.message}")
|
||||
print(f"Iterationen: {bestes.nit}")
|
||||
print(f"Erwartete Jahresrendite: {rendite * 100:6.2f} %")
|
||||
print(f"Erwartete Volatilitaet: {volatilitaet * 100:6.2f} %")
|
||||
print(f"Diversifikations-Entropie: {entropie:.4f} "
|
||||
f"(Maximum bei Gleichgewichtung: {np.log(n):.4f})")
|
||||
|
||||
print("\n" + pd.DataFrame({
|
||||
"Asset": ASSETS,
|
||||
"Erw. Rendite": [f"{r*100:.1f} %" for r in MU],
|
||||
"Volatilitaet": [f"{v*100:.1f} %" for v in VOLA],
|
||||
"Gewicht": [f"{w*100:6.2f} %" for w in w_opt],
|
||||
}).to_string(index=False))
|
||||
|
||||
# --- Vergleich: was passiert ohne Entropieterm? ----------------------
|
||||
ohne = optimiere(np.ones(n) / n, beta=0.0)
|
||||
print("\n--- Wirkung des Entropieterms ---")
|
||||
print(f" Mit Entropie (beta={BETA}): Gewichte {np.round(w_opt * 100, 1)}")
|
||||
print(f" Ohne Entropie (beta=0): Gewichte {np.round(ohne.x * 100, 1)}")
|
||||
print(" Der Entropieterm zieht Kapital aus der Spitzenposition heraus,")
|
||||
print(" ohne dass eine harte Obergrenze noetig waere.")
|
||||
print("=" * 74)
|
||||
85
OR_HTML_04/programme/Erklaerbarkeit.py
Normal file
85
OR_HTML_04/programme/Erklaerbarkeit.py
Normal file
|
|
@ -0,0 +1,85 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# Erklaerbarkeit.py
|
||||
"""
|
||||
Kapitel Praxisfallen: Ein Constraint-Trace, der die drei Fragen der Faustregel fuer
|
||||
Erklaerbarkeit (Abschnitt 'Die fuenf typischen Praxisfallen') fuer eine
|
||||
einzelne Zuweisung beantwortet:
|
||||
Warum diese Loesung? Warum nicht die Alternative? Was wuerde sie verbessern?
|
||||
"""
|
||||
|
||||
from ortools.sat.python import cp_model
|
||||
|
||||
PERSONAL = ["Anna", "Ben", "Clara"]
|
||||
BEREITS_STUNDEN = {"Anna": 1, "Ben": 2, "Clara": 0} # diese Woche schon geleistet
|
||||
MAX_STUNDEN = 2 # gesetzliche Obergrenze
|
||||
STRAFE_VORBELASTUNG = 50 # pro bereits geleisteter Stunde
|
||||
STRAFKOSTEN_WUNSCH = {"Anna": 0, "Ben": 30, "Clara": 100}
|
||||
|
||||
|
||||
def loese(max_stunden: int):
|
||||
modell = cp_model.CpModel()
|
||||
x = {p: modell.NewBoolVar(f"x_{p}") for p in PERSONAL}
|
||||
modell.AddExactlyOne(x.values())
|
||||
|
||||
kosten = []
|
||||
for p in PERSONAL:
|
||||
if BEREITS_STUNDEN[p] + 1 > max_stunden:
|
||||
modell.Add(x[p] == 0) # hartes Ausschlusskriterium
|
||||
strafe = BEREITS_STUNDEN[p] * STRAFE_VORBELASTUNG + STRAFKOSTEN_WUNSCH[p]
|
||||
kosten.append(strafe * x[p])
|
||||
|
||||
modell.Minimize(sum(kosten))
|
||||
loeser = cp_model.CpSolver()
|
||||
status = loeser.Solve(modell)
|
||||
if status not in (cp_model.OPTIMAL, cp_model.FEASIBLE):
|
||||
return None
|
||||
gewinner = next(p for p in PERSONAL if loeser.Value(x[p]))
|
||||
return gewinner, loeser.ObjectiveValue()
|
||||
|
||||
|
||||
def gesamtkosten(p: str) -> int:
|
||||
return BEREITS_STUNDEN[p] * STRAFE_VORBELASTUNG + STRAFKOSTEN_WUNSCH[p]
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
gewinner, kosten_opt = loese(MAX_STUNDEN)
|
||||
|
||||
print("=" * 78)
|
||||
print(" CONSTRAINT-TRACE: Warum bekommt", gewinner, "die Vertretungsstunde?")
|
||||
print("=" * 78)
|
||||
|
||||
print("\n1. WARUM DIESE LOESUNG?")
|
||||
v, w = BEREITS_STUNDEN[gewinner] * STRAFE_VORBELASTUNG, STRAFKOSTEN_WUNSCH[gewinner]
|
||||
print(f" {gewinner}: Vorbelastung {v} Strafpunkte + Wunsch {w} Strafpunkte "
|
||||
f"= {v + w} Strafpunkte insgesamt (Minimum).")
|
||||
|
||||
print("\n2. WARUM NICHT DIE NAHELIEGENDEN ALTERNATIVEN?")
|
||||
for p in PERSONAL:
|
||||
if p == gewinner:
|
||||
continue
|
||||
if BEREITS_STUNDEN[p] + 1 > MAX_STUNDEN:
|
||||
print(f" {p}: AUSGESCHLOSSEN - hat bereits {BEREITS_STUNDEN[p]} Stunden, "
|
||||
f"eine dritte würde die Obergrenze von {MAX_STUNDEN} Stunden verletzen "
|
||||
f"(harte Bedingung, nicht verhandelbar).")
|
||||
else:
|
||||
differenz = gesamtkosten(p) - kosten_opt
|
||||
print(f" {p}: zulässig, aber {differenz:.0f} Strafpunkte teurer als "
|
||||
f"{gewinner} ({gesamtkosten(p):.0f} statt {kosten_opt:.0f}).")
|
||||
|
||||
print("\n3. WAS WUERDE DIE LOESUNG VERBESSERN?")
|
||||
gewinner_gelockert, kosten_gelockert = loese(MAX_STUNDEN + 1)
|
||||
if kosten_gelockert < kosten_opt:
|
||||
print(f" Obergrenze auf {MAX_STUNDEN + 1} Stunden gelockert: {gewinner_gelockert} "
|
||||
f"uebernimmt jetzt zu {kosten_gelockert:.0f} Strafpunkten "
|
||||
f"({kosten_opt - kosten_gelockert:.0f} weniger als heute).")
|
||||
print(" -> Diese Bedingung hat tatsaechlich einen Preis.")
|
||||
else:
|
||||
print(f" Obergrenze auf {MAX_STUNDEN + 1} Stunden gelockert: Ergebnis bleibt bei "
|
||||
f"{gewinner} mit {kosten_gelockert:.0f} Strafpunkten (unveraendert).")
|
||||
blockiert = ", ".join(p for p in PERSONAL
|
||||
if BEREITS_STUNDEN[p] + 1 > MAX_STUNDEN and p != gewinner)
|
||||
print(f" -> Die Obergrenze blockiert zwar {blockiert}, kostet aber gerade NICHTS:")
|
||||
print(" Selbst ohne sie waere die Loesung dieselbe.")
|
||||
print(" Blockierend und kostenrelevant sind zwei verschiedene Dinge.")
|
||||
print("=" * 78)
|
||||
96
OR_HTML_04/programme/Excel_Bruecke.py
Normal file
96
OR_HTML_04/programme/Excel_Bruecke.py
Normal file
|
|
@ -0,0 +1,96 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# Excel_Bruecke.py
|
||||
"""
|
||||
Kapitel Einfuehrung: Vom Excel-Solver zu Python - dieselbe Rechnung, dieselbe Datei.
|
||||
|
||||
Das Programm bildet den vollstaendigen Arbeitsweg ab, den ein Excel-Modell im
|
||||
Betrieb sonst von Hand geht:
|
||||
|
||||
Tabellenblatt lesen -> Modell bauen -> loesen -> Ergebnis zurueckschreiben
|
||||
|
||||
Damit es ohne Vorbereitung laeuft, legt es die Eingabedatei beim ersten Start
|
||||
selbst an. Ersetzen Sie sie durch Ihre eigene Datei - solange die Spalten
|
||||
gleich heissen, aendert sich am Code nichts.
|
||||
|
||||
Die vier Arbeitsschritte kommen aus or_kern.py (Kapitel Praxisfallen). Was
|
||||
hier steht, ist nur noch das, was diesen Fall ausmacht - der Rest ist
|
||||
gemeinsamer Unterbau. Wer das Modul nicht kennt: Es ist im Buch vollstaendig
|
||||
abgedruckt und liegt im selben Verzeichnis.
|
||||
|
||||
Erzeugte Dateien (im Arbeitsverzeichnis):
|
||||
produktionsmix.xlsx Eingabe: Blaetter 'Produkte' und 'Kapazitaeten'
|
||||
produktionsmix_ergebnis.xlsx Ausgabe: Blaetter 'Plan' und 'Kennzahlen'
|
||||
|
||||
Benoetigt: pandas, openpyxl, ortools, pydantic (ueber or_kern)
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import os
|
||||
|
||||
import pandas as pd
|
||||
|
||||
from or_kern import (lade_produktionsproblem, loese_mit_glop, pruefe_loesung,
|
||||
schreibe_ergebnis)
|
||||
|
||||
EINGABE = "produktionsmix.xlsx"
|
||||
AUSGABE = "produktionsmix_ergebnis.xlsx"
|
||||
|
||||
|
||||
def lege_beispieldatei_an(pfad: str) -> None:
|
||||
"""Schreibt die Eingabemappe, wie sie auch aus dem Controlling kaeme."""
|
||||
produkte = pd.DataFrame({
|
||||
"Produkt": ["Tisch", "Stuhl", "Regal"],
|
||||
"Deckungsbeitrag": [240.0, 60.0, 130.0],
|
||||
"Montagestunden": [3.0, 1.0, 2.0],
|
||||
"Plattenmaterial": [6.0, 1.0, 4.0],
|
||||
})
|
||||
kapazitaeten = pd.DataFrame({
|
||||
"Ressource": ["Montagestunden", "Plattenmaterial"],
|
||||
"Verfuegbar": [150.0, 240.0],
|
||||
})
|
||||
with pd.ExcelWriter(pfad, engine="openpyxl") as mappe:
|
||||
produkte.to_excel(mappe, sheet_name="Produkte", index=False)
|
||||
kapazitaeten.to_excel(mappe, sheet_name="Kapazitaeten", index=False)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
if not os.path.exists(EINGABE):
|
||||
lege_beispieldatei_an(EINGABE)
|
||||
print(f"Beispiel-Eingabedatei angelegt: {EINGABE}")
|
||||
|
||||
# 1. Lesen und pruefen. Faellt hier etwas auf, ist der Fehler noch
|
||||
# zuzuordnen - spaeter waere er nur noch eine merkwuerdige Zahl.
|
||||
problem = lade_produktionsproblem(EINGABE)
|
||||
|
||||
# 2. Loesen. Einzige solverabhaengige Zeile des ganzen Programms -
|
||||
# 'loese_mit_scipy' waere ein Einzeiler-Wechsel.
|
||||
loesung = loese_mit_glop(problem)
|
||||
|
||||
# 3. Abnahmepruefung gegen die Anforderungen, ohne den Solver zu fragen.
|
||||
beanstandungen = pruefe_loesung(problem, loesung)
|
||||
if beanstandungen:
|
||||
raise SystemExit("Abnahmepruefung fehlgeschlagen:\n - "
|
||||
+ "\n - ".join(beanstandungen))
|
||||
|
||||
# 4. Zurueckschreiben - im Format, das die Fachabteilung ohnehin benutzt.
|
||||
schreibe_ergebnis(AUSGABE, problem, loesung)
|
||||
|
||||
print("=" * 66)
|
||||
print(" PRODUKTIONSPLAN AUS DER EXCEL-MAPPE")
|
||||
print("=" * 66)
|
||||
for produkt in problem.produkte:
|
||||
menge = loesung.werte[produkt.name]
|
||||
print(f"{produkt.name:<18} {menge:8.1f} Stueck "
|
||||
f"{menge * produkt.deckungsbeitrag:12,.2f} EUR")
|
||||
print("-" * 66)
|
||||
print(loesung.als_bericht())
|
||||
print("Abnahmepruefung: bestanden\n")
|
||||
|
||||
# Die Schattenpreise beantworten die Frage, die der Plan nicht stellt:
|
||||
# Was waere eine zusaetzliche Einheit dieser Ressource wert?
|
||||
for ressource, preis in loesung.schattenpreise.items():
|
||||
print(f" eine Einheit {ressource:<18} mehr waere wert: {preis:7.2f} EUR")
|
||||
print("=" * 66)
|
||||
print(f"Ergebnis geschrieben nach: {AUSGABE}")
|
||||
147
OR_HTML_04/programme/Finanzdaten_Ledoit_Wolf.py
Normal file
147
OR_HTML_04/programme/Finanzdaten_Ledoit_Wolf.py
Normal file
|
|
@ -0,0 +1,147 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# Finanzdaten_Ledoit_Wolf.py
|
||||
"""
|
||||
Kapitel Finanzdaten: Automatisierte Finanzdaten-Pipeline und Kovarianz-Shrinkage.
|
||||
|
||||
Eigenschaften:
|
||||
* Spaltenreihenfolge wird erzwungen und per assert geprueft
|
||||
* beide Shrinkage-Ziele werden berechnet und verglichen
|
||||
* Datenqualitaetspruefungen (Luecken, Ausreisser, Mindestlaenge)
|
||||
* Eigenwertspektrum wird ausgewiesen, nicht nur die Konditionszahl
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import yfinance as yf
|
||||
from sklearn.covariance import LedoitWolf
|
||||
|
||||
HANDELSTAGE = 252
|
||||
|
||||
|
||||
def lade_kurse(tickers: list[str], start, ende) -> pd.DataFrame:
|
||||
"""Laedt bereinigte Schlusskurse und garantiert die Spaltenreihenfolge."""
|
||||
print(f"Lade {len(tickers)} Ticker von {start} bis {ende} ...")
|
||||
roh = yf.download(tickers, start=start, end=ende, auto_adjust=True, progress=False)
|
||||
if roh.empty:
|
||||
raise SystemExit("Download fehlgeschlagen (Netz, Ticker oder Rate-Limit pruefen).")
|
||||
|
||||
if isinstance(roh.columns, pd.MultiIndex):
|
||||
# [tickers] erzwingt die gewuenschte Spaltenreihenfolge
|
||||
kurse = roh["Close"][tickers].dropna()
|
||||
else:
|
||||
kurse = roh[["Close"]].dropna()
|
||||
kurse.columns = tickers
|
||||
|
||||
assert list(kurse.columns) == tickers, (
|
||||
f"Spaltenreihenfolge weicht ab!\n erwartet: {tickers}\n"
|
||||
f" erhalten: {list(kurse.columns)}")
|
||||
return kurse
|
||||
|
||||
|
||||
def pruefe_datenqualitaet(kurse: pd.DataFrame, mindest_tage: int = 200) -> None:
|
||||
"""Faengt die haeufigsten Datenprobleme ab, bevor sie ins Modell gelangen."""
|
||||
T, N = kurse.shape
|
||||
print(f"\nDatenqualitaet: {T} Handelstage, {N} Titel (T/N = {T/N:.1f})")
|
||||
if T < mindest_tage:
|
||||
raise SystemExit(f"Zu wenige Beobachtungen ({T} < {mindest_tage}).")
|
||||
if T < N:
|
||||
print(" WARNUNG: T < N - die Stichprobenkovarianz ist singulaer!")
|
||||
elif T < 3 * N:
|
||||
print(" WARNUNG: T < 3N - erhebliches Schaetzrauschen zu erwarten.")
|
||||
|
||||
renditen = kurse.pct_change().dropna()
|
||||
extreme = (renditen.abs() > 0.25).sum().sum()
|
||||
if extreme:
|
||||
print(f" Hinweis: {extreme} Tagesrenditen ueber 25 % "
|
||||
f"(Splits, Sondersituationen oder Datenfehler pruefen).")
|
||||
luecken = kurse.isna().sum().sum()
|
||||
print(f" Fehlende Werte nach dropna: {luecken}")
|
||||
|
||||
|
||||
def ziel_konstante_korrelation(renditen: np.ndarray) -> np.ndarray:
|
||||
"""
|
||||
Shrinkage-Ziel nach Ledoit/Wolf 2003: Einzelvarianzen behalten,
|
||||
alle Korrelationen durch ihren Mittelwert ersetzen.
|
||||
(sklearn nutzt stattdessen die skalierte Einheitsmatrix - siehe oben.)
|
||||
"""
|
||||
S = np.cov(renditen, rowvar=False, ddof=1)
|
||||
d = np.sqrt(np.diag(S))
|
||||
korrelation = S / np.outer(d, d)
|
||||
n = len(S)
|
||||
r_quer = (korrelation.sum() - n) / (n * (n - 1)) # Mittel ohne Diagonale
|
||||
F = r_quer * np.outer(d, d)
|
||||
np.fill_diagonal(F, np.diag(S))
|
||||
return F
|
||||
|
||||
|
||||
def analysiere(kurse: pd.DataFrame):
|
||||
renditen = kurse.pct_change().dropna()
|
||||
T, N = renditen.shape
|
||||
|
||||
mu = renditen.mean().values * HANDELSTAGE
|
||||
sigma_stichprobe = renditen.cov().values * HANDELSTAGE
|
||||
|
||||
lw = LedoitWolf(assume_centered=False).fit(renditen.values)
|
||||
sigma_lw = lw.covariance_ * HANDELSTAGE
|
||||
delta = lw.shrinkage_
|
||||
|
||||
print("\n" + "=" * 84)
|
||||
print(" FINANZDATEN-PIPELINE UND MATRIX-KONDITIONIERUNG")
|
||||
print("=" * 84)
|
||||
print(f"Beobachtungen T: {T}")
|
||||
print(f"Titel N: {N}")
|
||||
print(f"Optimales Shrinkage delta: {delta:.4f} "
|
||||
f"({delta*100:.1f} % Gewicht auf dem strukturierten Ziel)")
|
||||
|
||||
eig_stich = np.linalg.eigvalsh(sigma_stichprobe)
|
||||
eig_lw = np.linalg.eigvalsh(sigma_lw)
|
||||
|
||||
print("\n--- Eigenwertspektrum (annualisiert) ---")
|
||||
print(f"{'':<14} {'kleinster':>12} {'Median':>12} {'groesster':>12} "
|
||||
f"{'Kondition':>12}")
|
||||
for name, eig in [("Stichprobe", eig_stich), ("Ledoit-Wolf", eig_lw)]:
|
||||
kondition = eig.max() / max(eig.min(), 1e-12)
|
||||
print(f"{name:<14} {eig.min():>12.6f} {np.median(eig):>12.6f} "
|
||||
f"{eig.max():>12.6f} {kondition:>12.1f}")
|
||||
|
||||
print("\nDie Konditionszahl misst, wie stark sich kleine Datenaenderungen auf")
|
||||
print("die Inverse auswirken. Je kleiner, desto stabiler die Portfoliogewichte.")
|
||||
|
||||
# --- Vergleich der beiden Shrinkage-Ziele ----------------------------
|
||||
S_taeglich = np.cov(renditen.values, rowvar=False, ddof=1)
|
||||
F_identitaet = np.eye(N) * np.trace(S_taeglich) / N
|
||||
F_korrelation = ziel_konstante_korrelation(renditen.values)
|
||||
|
||||
print("\n--- Die beiden Shrinkage-Ziele im Vergleich ---")
|
||||
for name, F in [("skalierte Einheitsmatrix (sklearn)", F_identitaet),
|
||||
("konstante Korrelation (LW 2003)", F_korrelation)]:
|
||||
gemischt = ((1 - delta) * S_taeglich + delta * F) * HANDELSTAGE
|
||||
eig = np.linalg.eigvalsh(gemischt)
|
||||
print(f" {name:<36}: Kondition {eig.max()/max(eig.min(),1e-12):8.1f}, "
|
||||
f"kleinster EW {eig.min():.6f}")
|
||||
|
||||
print("\n--- Erwartete Renditen (annualisiert) ---")
|
||||
uebersicht = pd.DataFrame({
|
||||
"Ticker": kurse.columns,
|
||||
"Rendite p.a.": [f"{r*100:+7.2f} %" for r in mu],
|
||||
"Volatilitaet p.a.": [f"{np.sqrt(sigma_lw[i, i])*100:6.2f} %" for i in range(N)],
|
||||
})
|
||||
print(uebersicht.to_string(index=False))
|
||||
print("=" * 84)
|
||||
return renditen, mu, sigma_stichprobe, sigma_lw
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
UNIVERSUM = ["AAPL", "MSFT", "NVDA", "AMZN", # Technologie
|
||||
"JNJ", "PFE", # Gesundheit
|
||||
"JPM", "GS", # Banken
|
||||
"XOM", "CVX"] # Energie
|
||||
|
||||
# Relatives 2-Jahres-Fenster, damit das Beispiel nicht "altert"
|
||||
ende = pd.Timestamp.today().normalize()
|
||||
start = ende - pd.DateOffset(years=2)
|
||||
|
||||
kurse = lade_kurse(UNIVERSUM, start.strftime("%Y-%m-%d"), ende.strftime("%Y-%m-%d"))
|
||||
pruefe_datenqualitaet(kurse)
|
||||
analysiere(kurse)
|
||||
53
OR_HTML_04/programme/Fluch_des_Durchschnitts.py
Normal file
53
OR_HTML_04/programme/Fluch_des_Durchschnitts.py
Normal file
|
|
@ -0,0 +1,53 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
# Fluch_des_Durchschnitts.py
|
||||
"""
|
||||
Kapitel Unsicherheit: Der Fluch des Durchschnitts (die Handrechnung dazu) als Scan ueber alle
|
||||
moeglichen Kapazitaeten - zeigt, dass das Optimum nicht beim Mittelwert liegt.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
|
||||
SZENARIEN = np.array([100, 250, 500])
|
||||
WAHRSCHEINLICHKEITEN = np.array([0.5, 0.3, 0.2])
|
||||
PREIS_VORAB = 40.0
|
||||
PREIS_ZUKAUF = 120.0
|
||||
PREIS_VERWALTUNG = 5.0
|
||||
|
||||
ERWARTUNGSWERT = float(SZENARIEN @ WAHRSCHEINLICHKEITEN)
|
||||
|
||||
|
||||
def erwartete_kosten(x):
|
||||
kosten = np.where(
|
||||
SZENARIEN >= x,
|
||||
PREIS_VORAB * x + PREIS_ZUKAUF * (SZENARIEN - x),
|
||||
PREIS_VORAB * x + PREIS_VERWALTUNG * (x - SZENARIEN),
|
||||
)
|
||||
return float(kosten @ WAHRSCHEINLICHKEITEN)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
print("=" * 70)
|
||||
print(" DER FLUCH DES DURCHSCHNITTS")
|
||||
print("=" * 70)
|
||||
print(f"Erwarteter Bedarf (Mittelwert): {ERWARTUNGSWERT:.1f}")
|
||||
print(f"Kosten bei naiver Planung x={ERWARTUNGSWERT:.0f}: "
|
||||
f"{erwartete_kosten(ERWARTUNGSWERT):,.2f} EUR\n")
|
||||
|
||||
x_werte = np.arange(0, 501, 1)
|
||||
kosten_werte = np.array([erwartete_kosten(x) for x in x_werte])
|
||||
x_optimal = x_werte[np.argmin(kosten_werte)]
|
||||
kosten_optimal = kosten_werte.min()
|
||||
|
||||
print(f"Optimales x (durch Scan gefunden): {x_optimal}")
|
||||
print(f"Kosten beim Optimum: {kosten_optimal:,.2f} EUR")
|
||||
print(f"Ersparnis gegenueber naiver Planung: "
|
||||
f"{erwartete_kosten(ERWARTUNGSWERT) - kosten_optimal:,.2f} EUR\n")
|
||||
|
||||
print("Kosten fuer ausgewaehlte x zum Vergleich:")
|
||||
for x in [100, 225, 250, 300, 500]:
|
||||
markierung = " <- Mittelwert" if x == 225 else (" <- Optimum" if x == 250 else "")
|
||||
print(f" x={x:4d}: {erwartete_kosten(x):>12,.2f} EUR{markierung}")
|
||||
|
||||
print("\nDas Optimum liegt exakt auf einem Szenariowert (250 = 'Volatil'),")
|
||||
print("nicht beim Mittelwert 225 - typisch fuer asymmetrische Kostenfunktionen.")
|
||||
Some files were not shown because too many files have changed in this diff Show more
Loading…
Reference in a new issue