#!/usr/bin/env python3 # MILP_Portfolio_Fixgebuehren.py """ Kapitel MILP: MILP-Portfolio-Selektion mit Fixkosten und Kardinalitaet. CP-SAT-freie, gut lesbare Formulierung ueber scipy/HiGHS (kein manueller CSR-Matrixaufbau, deutlich leichter nachvollziehbar), mit Vergleich gegen die Loesung OHNE Restriktionen. Die Auswertung laeuft ueber SolverStatus und das Loesung-Objekt aus or_kern.py: Alle Statusfaelle werden behandelt, und der MIP-Gap steht im Bericht - statt eines blossen 'success', das nicht verraet, ob der Solver fertig geworden ist oder nur aufgegeben hat. Benoetigt: numpy, pandas, scipy, pydantic (ueber or_kern) """ from __future__ import annotations import time import numpy as np import pandas as pd from scipy.optimize import linprog from or_kern import Loesung, SolverStatus, status_von_scipy ANLAGEN = ["US-Aktien", "EU-Aktien", "Emerging-Markets", "Staatsanleihen", "Unternehmensanleihen", "Rohstoffe"] RENDITE = np.array([0.11, 0.08, 0.13, 0.03, 0.05, 0.07]) # erwartet, p.a. BUDGET = 100_000.0 MIN_POSITION = 10_000.0 # L MAX_POSITION = 40_000.0 # U (dient zugleich als Big-M!) GEBUEHR = 50.0 # F, je aktivierter Position MAX_POSITIONEN = 3 # K N = len(ANLAGEN) def baue_und_loese(zeitlimit: float | None = None) -> Loesung: """ Variablenreihenfolge: [x_0..x_{N-1}, y_0..y_{N-1}] Zielfunktion (Maximierung -> fuer linprog negiert): max sum(rendite_i * x_i) - GEBUEHR * sum(y_i) """ c = np.concatenate([-RENDITE, np.full(N, GEBUEHR)]) # negiert = Minimierung # --- Gleichungsnebenbedingung: Budget vollstaendig investiert ----------- A_eq = np.zeros((1, 2 * N)) A_eq[0, :N] = 1.0 b_eq = np.array([BUDGET]) zeilen, grenzen = [], [] # --- Kardinalitaet: sum(y_i) <= K -------------------------------------- zeile = np.zeros(2 * N) zeile[N:] = 1.0 zeilen.append(zeile) grenzen.append(MAX_POSITIONEN) # --- Obergrenze (Big-M): x_i - U*y_i <= 0 ----------------------------- for i in range(N): zeile = np.zeros(2 * N) zeile[i] = 1.0 zeile[N + i] = -MAX_POSITION zeilen.append(zeile) grenzen.append(0.0) # --- Untergrenze: L*y_i - x_i <= 0 (entspricht x_i >= L*y_i) -------- for i in range(N): zeile = np.zeros(2 * N) zeile[i] = -1.0 zeile[N + i] = MIN_POSITION zeilen.append(zeile) grenzen.append(0.0) A_ub = np.array(zeilen) b_ub = np.array(grenzen) schranken = [(0.0, MAX_POSITION)] * N + [(0.0, 1.0)] * N ganzzahligkeit = np.concatenate([np.zeros(N), np.ones(N)]) # y binaer t0 = time.perf_counter() ergebnis = linprog(c=c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=schranken, integrality=ganzzahligkeit, method="highs", options={"time_limit": zeitlimit} if zeitlimit else None) laufzeit = time.perf_counter() - t0 status = status_von_scipy(ergebnis) if not status.brauchbar: return Loesung(status=status, laufzeit=laufzeit) # Zurueck in die Maximierungswelt: Zielwert UND Schranke negieren. # Aus beiden zusammen rechnet das Loesung-Objekt den MIP-Gap aus. x, y = ergebnis.x[:N], np.round(ergebnis.x[N:]) return Loesung( status=status, werte={**{name: float(w) for name, w in zip(ANLAGEN, x)}, **{f"aktiv:{name}": float(w) for name, w in zip(ANLAGEN, y)}}, zielwert=float(-ergebnis.fun), schranke=float(-ergebnis.mip_dual_bound), laufzeit=laufzeit) def ohne_restriktionen() -> tuple[float, np.ndarray]: """Vergleichsfall: nur Budget, keine Gebuehren/Kardinalitaet/Mindestgroesse.""" ergebnis = linprog(c=-RENDITE, A_eq=[np.ones(N)], b_eq=[BUDGET], bounds=[(0, None)] * N, method="highs") status = status_von_scipy(ergebnis) if not status.brauchbar: raise SystemExit(f"Vergleichsfall nicht loesbar: {status.value}") return -ergebnis.fun, ergebnis.x def pruefe_portfolio(loesung: Loesung, toleranz: float = 1e-6) -> list[str]: """Prueft die Loesung gegen die Anforderungen - ohne den Solver zu fragen. Bewusst kein assert: Eine Beanstandungsliste laesst sich protokollieren, weiterreichen und testen. Ein assert verschwindet ausserdem, sobald jemand Python mit -O startet. """ x = np.array([loesung.werte[name] for name in ANLAGEN]) y = np.array([loesung.werte[f"aktiv:{name}"] for name in ANLAGEN]) beanstandungen: list[str] = [] if abs(x.sum() - BUDGET) > 1e-4: beanstandungen.append(f"Budget nicht exakt investiert: {x.sum():,.2f}") if y.sum() > MAX_POSITIONEN + toleranz: beanstandungen.append(f"{y.sum():.0f} Positionen statt hoechstens " f"{MAX_POSITIONEN}") for i, name in enumerate(ANLAGEN): if abs(y[i] - round(y[i])) > toleranz: beanstandungen.append(f"{name}: y = {y[i]!r} ist nicht ganzzahlig") elif y[i] and not (MIN_POSITION - toleranz <= x[i] <= MAX_POSITION + toleranz): beanstandungen.append(f"{name}: {x[i]:,.2f} EUR verletzt die " f"Groessengrenzen") elif not y[i] and x[i] > toleranz: beanstandungen.append(f"{name}: inaktiv, aber {x[i]:,.2f} EUR " f"investiert") return beanstandungen if __name__ == "__main__": loesung = baue_und_loese() print("=" * 88) print(" OPTIMALE MILP-PORTFOLIO-ALLOKATION MIT FIXGEBUEHREN") print("=" * 88) print(f"Budget: {BUDGET:,.0f} EUR | max. {MAX_POSITIONEN} Positionen | " f"je {MIN_POSITION:,.0f}-{MAX_POSITION:,.0f} EUR | Gebuehr {GEBUEHR:.0f} EUR\n") # Zuerst der Status - erst danach interessieren die Zahlen. if loesung.status.modellfehler: raise SystemExit(f"Das Modell ist nicht loesbar ({loesung.status.value}). " f"Naechster Schritt: Anhang Fehlerdiagnose.") if not loesung.status.brauchbar: raise SystemExit(f"Keine Loesung erhalten ({loesung.status.value}). " f"Zeitlimit erhoehen oder Modell vereinfachen.") if loesung.status is SolverStatus.ZULAESSIG: print(f"ACHTUNG: nicht beweisbar optimal - Gap {loesung.gap:.2%}\n") x = np.array([loesung.werte[name] for name in ANLAGEN]) y = np.array([loesung.werte[f"aktiv:{name}"] for name in ANLAGEN], dtype=int) tabelle = pd.DataFrame({ "Anlage": ANLAGEN, "Aktiv": ["JA" if y[i] else "-" for i in range(N)], "Investition (EUR)": [f"{x[i]:,.0f}" for i in range(N)], "Anteil": [f"{x[i]/BUDGET*100:5.1f} %" for i in range(N)], "Erw. Rendite": [f"{RENDITE[i]*100:4.1f} %" for i in range(N)], "Erw. Ertrag (EUR)": [f"{x[i]*RENDITE[i]:,.0f}" for i in range(N)], }) print(tabelle.to_string(index=False)) brutto = float(RENDITE @ x) gebuehren = float(GEBUEHR * y.sum()) print("-" * 88) print(f"Erwarteter Bruttoertrag: {brutto:>12,.2f} EUR") print(f"Ordergebuehren: {-gebuehren:>12,.2f} EUR ({y.sum()} Positionen)") print(f"Netto-Erwartungswert: {loesung.zielwert:>12,.2f} EUR") print(f"\n{loesung.als_bericht()}") # --- Vergleich mit dem unbeschraenkten Fall --------------------------- z_frei, x_frei = ohne_restriktionen() print("-" * 88) print(f"Zum Vergleich ohne jede Restriktion (alles in den Bestwert): " f"{z_frei:,.2f} EUR") print(f"Kosten der Realitaet (Gebuehren, Streuung, Mindestgroessen): " f"{z_frei - loesung.zielwert:,.2f} EUR " f"({(1 - loesung.zielwert/z_frei)*100:.2f} %)") # --- Alle Nebenbedingungen nachpruefen -------------------------------- beanstandungen = pruefe_portfolio(loesung) print("-" * 88) if beanstandungen: raise SystemExit("Abnahmepruefung fehlgeschlagen:\n - " + "\n - ".join(beanstandungen)) print("Abnahmepruefung: alle Nebenbedingungen geprueft und eingehalten.") print("=" * 88)