{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Kapitel 16: Die Strukturbrücke — dieselbe Mathematik, zwei Welten\n", "\n", "Begleitnotebook zu *Optimierte Entscheidungsfindung mit Python*. Die Codezellen sind identisch mit den im Buch abgedruckten Programmen.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Einmalig ausfuehren: installiert alle im Buch verwendeten Pakete.\n", "# Lokal in einer virtuellen Umgebung genauso gueltig wie in Google Colab.\n", "%pip install --quiet ortools highspy cvxpy scipy numpy pandas polars \\\n", " scikit-learn matplotlib plotly pyomo linopy pymoo pydantic openpyxl" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Das Programm\n", "\n", "`Strukturbruecke.py`\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#!/usr/bin/env python3\n", "\n", "# Strukturbruecke.py\n", "\"\"\"\n", "Kapitel Bruecke: Derselbe Code, zwei Welten - der Beweis statt der Behauptung.\n", "\n", "Der Teil-Auftakt stellt eine Tabelle auf: \"Ressourcen auf Produkte verteilen\"\n", "entspreche \"Kapital auf Anlagen verteilen\", \"gegen den Worst Case absichern\"\n", "entspreche \"Absicherung gegen Kursabstuerze\". Solche Tabellen stehen in vielen\n", "Buechern. Sie sind billig - und man kann sie pruefen.\n", "\n", "Dieses Programm prueft sie. Es fuettert zweimal DENSELBEN Code mit Daten aus\n", "zwei Welten und zeigt die Ergebnisse nebeneinander:\n", "\n", " 1. Die Allokation als LP. Das Domaenenmodell aus or_kern.py, einmal mit einer\n", " Schreinerei (Montagestunden, Plattenmaterial) und einmal mit einem Depot\n", " (Kapital, Risikobudget). Gleiche Klasse, gleicher Modellbauer, gleiche\n", " Abnahmepruefung - und der Schattenpreis heisst in der einen Welt\n", " \"Wert einer zusaetzlichen Montagestunde\" und in der anderen \"Preis des\n", " Risikos\".\n", " 2. Die Absicherung gegen den schlechtesten Fall als CVaR. EINE Funktion,\n", " einmal mit Kursrenditen und einmal mit Lieferverzuegen.\n", " 3. Was NICHT hinueberreicht. Der ehrliche Teil: Drei Unterschiede, die die\n", " Analogie begrenzen - und die man kennen muss, bevor man Methoden aus\n", " der einen Welt in die andere traegt.\n", "\n", "ZUR SOLVERWAHL: Teil 1 benutzt 'loese_mit_scipy' und nicht 'loese_mit_glop'.\n", "Das ist kein Zufall - CVXPY laedt fuer Teil 2 highspy, und ortools vertraegt\n", "sich damit nicht im selben Prozess (Kapitel Oekosystem). Wer hier GLOP nimmt,\n", "bekommt beim cvxpy-Import eine Fehlermeldung ueber ein 'undefined symbol'.\n", "Genau deshalb laedt or_kern.py seine Solver erst beim Aufruf.\n", "\n", "Benoetigt: numpy, cvxpy, scipy und pydantic (ueber or_kern)\n", "\"\"\"\n", "\n", "from __future__ import annotations\n", "\n", "import numpy as np\n", "\n", "from or_kern import (Produkt, Produktionsproblem, loese_mit_scipy,\n", " pruefe_loesung)\n", "\n", "SAAT = 7\n", "SZENARIEN = 500\n", "ALPHA = 0.95 # die schlechtesten 5 % der Faelle\n", "MAX_ANTEIL = 0.40 # Streuungsgebot in beiden Welten\n", "\n", "\n", "# --- Teil 1: Dieselbe Klasse, zwei Welten ----------------------------------\n", "\n", "def schreinerei() -> Produktionsproblem:\n", " \"\"\"Montagestunden und Plattenmaterial auf Tische und Stuehle verteilen.\"\"\"\n", " return Produktionsproblem(\n", " produkte=[\n", " Produkt(name=\"Tisch\", deckungsbeitrag=240.0,\n", " verbrauch={\"Montagestunden\": 3.0, \"Plattenmaterial\": 6.0}),\n", " Produkt(name=\"Stuhl\", deckungsbeitrag=60.0,\n", " verbrauch={\"Montagestunden\": 1.0, \"Plattenmaterial\": 1.0}),\n", " Produkt(name=\"Regal\", deckungsbeitrag=130.0,\n", " verbrauch={\"Montagestunden\": 2.0, \"Plattenmaterial\": 4.0}),\n", " ],\n", " kapazitaeten={\"Montagestunden\": 150.0, \"Plattenmaterial\": 240.0})\n", "\n", "\n", "def depot() -> Produktionsproblem:\n", " \"\"\"Kapital und Risikobudget auf Anlageklassen verteilen.\n", "\n", " Eine \"Einheit\" ist hier 1.000 EUR Anlagesumme. Der Deckungsbeitrag ist der\n", " erwartete Jahresertrag dieser Einheit, der Verbrauch die beanspruchte\n", " Kapital- und Risikomenge. Das ist keine Analogie, sondern buchstaeblich\n", " dasselbe Modell - deshalb passt es in dieselbe Klasse.\n", " \"\"\"\n", " return Produktionsproblem(\n", " produkte=[\n", " Produkt(name=\"Aktien Welt\", deckungsbeitrag=75.0,\n", " verbrauch={\"Kapital (Tsd. EUR)\": 1.0, \"Risikobudget\": 1.00}),\n", " Produkt(name=\"Anleihen\", deckungsbeitrag=28.0,\n", " verbrauch={\"Kapital (Tsd. EUR)\": 1.0, \"Risikobudget\": 0.22}),\n", " Produkt(name=\"Immobilienfonds\", deckungsbeitrag=46.0,\n", " verbrauch={\"Kapital (Tsd. EUR)\": 1.0, \"Risikobudget\": 0.55}),\n", " ],\n", " kapazitaeten={\"Kapital (Tsd. EUR)\": 150.0, \"Risikobudget\": 90.0})\n", "\n", "\n", "def berichte_allokation(titel: str, problem: Produktionsproblem,\n", " einheit: str, ertragsname: str) -> None:\n", " \"\"\"Ein Bericht fuer beide Welten - nur die Beschriftung wechselt.\"\"\"\n", " loesung = loese_mit_scipy(problem)\n", " beanstandungen = pruefe_loesung(problem, loesung)\n", "\n", " print(f\" {titel}\")\n", " print(f\" {loesung.als_bericht()}\")\n", " for produkt in problem.produkte:\n", " print(f\" {produkt.name:<18} {loesung.werte[produkt.name]:8.2f} {einheit}\")\n", " print(f\" {ertragsname:<18} {loesung.zielwert:8.2f} EUR\")\n", " for ressource, preis in loesung.schattenpreise.items():\n", " print(f\" Schattenpreis {ressource:<24} {preis:7.2f} EUR\")\n", " print(f\" Abnahmepruefung: \"\n", " f\"{'bestanden' if not beanstandungen else beanstandungen}\")\n", "\n", "\n", "# --- Teil 2: Eine CVaR-Funktion, zwei Welten -------------------------------\n", "\n", "def optimiere_cvar(verluste: np.ndarray, ertrag: np.ndarray,\n", " mindestertrag: float, alpha: float = ALPHA,\n", " max_anteil: float = MAX_ANTEIL):\n", " \"\"\"Minimiert den CVaR der Verluste unter einer Mindestertragsbedingung.\n", "\n", " 'verluste' hat die Form (Szenarien x Optionen) und enthaelt, was in\n", " Szenario s passiert, wenn eine Einheit in Option j steckt. Was ein\n", " \"Verlust\" ist, entscheidet allein die Einheit der Matrix: Prozentpunkte\n", " Kursverlust oder Tage Lieferverzug - die Formel sieht keinen Unterschied.\n", "\n", " Das ist die Rockafellar-Uryasev-Formulierung aus dem Kapitel CVaR, hier\n", " ohne jede Aenderung wiederverwendet.\n", " \"\"\"\n", " import cvxpy as cp\n", "\n", " anzahl_szenarien, anzahl_optionen = verluste.shape\n", " anteil = cp.Variable(anzahl_optionen, nonneg=True)\n", " schwelle = cp.Variable() # wird im Optimum zum VaR\n", " ueberschuss = cp.Variable(anzahl_szenarien, nonneg=True)\n", "\n", " cvar = schwelle + (1.0 / (anzahl_szenarien * (1 - alpha))) * cp.sum(ueberschuss)\n", " problem = cp.Problem(\n", " cp.Minimize(cvar),\n", " [ueberschuss >= verluste @ anteil - schwelle,\n", " cp.sum(anteil) == 1,\n", " anteil <= max_anteil,\n", " ertrag @ anteil >= mindestertrag])\n", " problem.solve()\n", " if problem.status not in (\"optimal\", \"optimal_inaccurate\"):\n", " raise SystemExit(f\"CVaR-Problem nicht loesbar: {problem.status}\")\n", " return anteil.value, float(cvar.value), float(schwelle.value)\n", "\n", "\n", "def kursszenarien(rng) -> tuple[np.ndarray, np.ndarray, list[str]]:\n", " \"\"\"Taegliche Verluste (negative Renditen) von sechs Anlageklassen.\"\"\"\n", " namen = [\"Aktien Welt\", \"Aktien EU\", \"Schwellenlaender\",\n", " \"Staatsanleihen\", \"Unternehmensanl.\", \"Rohstoffe\"]\n", " rendite_pa = np.array([0.080, 0.065, 0.110, 0.025, 0.045, 0.070])\n", " schwankung = np.array([0.180, 0.160, 0.260, 0.040, 0.075, 0.210])\n", " taeglich = rng.normal(rendite_pa / 252, schwankung / np.sqrt(252),\n", " (SZENARIEN, len(namen)))\n", " return -taeglich * 100.0, rendite_pa * 100.0, namen\n", "\n", "\n", "def lieferszenarien(rng) -> tuple[np.ndarray, np.ndarray, list[str]]:\n", " \"\"\"Lieferverzug in Tagen bei sechs Lieferanten.\n", "\n", " Der Aufbau ist bewusst anders als bei den Kursen: Hier gibt es einen\n", " normalen Verzug UND seltene Totalausfaelle, die 14 Tage kosten. Das ist\n", " genau die Art fetter Raender, wegen der man in beiden Welten CVaR statt\n", " Standardabweichung benutzt.\n", " \"\"\"\n", " namen = [\"Nordwerk\", \"Sued-Metall\", \"Fernost A\", \"Lokalzulieferer\",\n", " \"Fernost B\", \"Osteuropa\"]\n", " zuverlaessigkeit = np.array([0.94, 0.97, 0.90, 0.99, 0.92, 0.95])\n", " verzugsstreuung = np.array([3.5, 1.8, 5.0, 0.9, 4.2, 2.8])\n", " normaler_verzug = np.maximum(0.0, rng.normal(0.5, 1.0, (SZENARIEN, len(namen)))\n", " * verzugsstreuung)\n", " totalausfall = (rng.random((SZENARIEN, len(namen))) > zuverlaessigkeit) * 14.0\n", " return normaler_verzug + totalausfall, zuverlaessigkeit, namen\n", "\n", "\n", "def berichte_cvar(titel: str, verluste, ertrag, namen, mindestertrag,\n", " einheit: str, ertragsname: str):\n", " anteile, cvar, var = optimiere_cvar(verluste, ertrag, mindestertrag)\n", " mittel = float((verluste @ anteile).mean())\n", " print(f\" {titel}\")\n", " print(f\" VaR {ALPHA:.0%}: {var:8.3f} {einheit}\")\n", " print(f\" CVaR {ALPHA:.0%}: {cvar:8.3f} {einheit} \"\n", " f\"(Mittelwert ueber alle Szenarien: {mittel:.3f})\")\n", " print(f\" {ertragsname}: {float(ertrag @ anteile):.3f}\")\n", " print(\" Aufteilung:\")\n", " for name, anteil in zip(namen, anteile):\n", " balken = \"#\" * int(round(anteil * 40))\n", " grenze = \" <- an der Streuungsgrenze\" if anteil > MAX_ANTEIL - 1e-4 else \"\"\n", " print(f\" {name:<18} {anteil:6.1%} {balken}{grenze}\")\n", " return anteile\n", "\n", "\n", "if __name__ == \"__main__\":\n", " rng = np.random.default_rng(SAAT)\n", "\n", " print(\"=\" * 80)\n", " print(\" DIESELBE STRUKTUR, ZWEI WELTEN\")\n", " print(\"=\" * 80)\n", "\n", " # --- Teil 1 ----------------------------------------------------------\n", " print(\"\\n1. Allokation als LP - EINE Klasse, EIN Modellbauer\\n\")\n", " berichte_allokation(\"Werkstatt: Produktionsprogramm\", schreinerei(),\n", " \"Stueck\", \"Deckungsbeitrag\")\n", " print()\n", " berichte_allokation(\"Depot: Anlageaufteilung\", depot(),\n", " \"Tsd. \", \"Erwarteter Ertrag\")\n", "\n", " print(\"\\n Beide Ausgaben stammen aus derselben Funktion \"\n", " \"'berichte_allokation'.\")\n", " print(\" Ausgetauscht wurden nur die Daten und die Beschriftungen -\")\n", " print(\" keine Zeile Modellcode.\")\n", " print()\n", " print(\" Lesen Sie die Schattenpreise nebeneinander: In der Werkstatt sagt\")\n", " print(\" er, was eine zusaetzliche Montagestunde wert waere. Im Depot sagt\")\n", " print(\" dieselbe Zahl, was eine zusaetzliche Einheit Risikobudget wert\")\n", " print(\" waere - der PREIS DES RISIKOS. Das ist kein Sprachbild, sondern\")\n", " print(\" derselbe Dualwert derselben Nebenbedingung.\")\n", "\n", " # --- Teil 2 ----------------------------------------------------------\n", " print(\"\\n\" + \"-\" * 80)\n", " print(\"2. Absicherung gegen den schlechtesten Fall - EINE CVaR-Funktion\\n\")\n", "\n", " kurse, rendite, anlagen = kursszenarien(rng)\n", " anteile_depot = berichte_cvar(\n", " \"Depot: die schlechtesten 5 % der Handelstage\",\n", " kurse, rendite, anlagen, 5.5,\n", " \"% je Tag \", \"Erwartete Jahresrendite (%)\")\n", " print()\n", " lieferung, zuverlaessig, lieferanten = lieferszenarien(rng)\n", " anteile_einkauf = berichte_cvar(\"Einkauf: die schlechtesten 5 % der Bestellungen\",\n", " lieferung, zuverlaessig, lieferanten, 0.945,\n", " \"Tage Verzug \", \"Mittlere Zuverlaessigkeit\")\n", "\n", " print(\"\\n Auch hier: eine Funktion, zwei Aufrufe. Die Zielfunktion\")\n", " print(\" interessiert sich nicht dafuer, ob in der Matrix Prozentpunkte\")\n", " print(\" oder Tage stehen.\")\n", " print()\n", " print(f\" Interessant ist, WO die Streuungsgrenze von {MAX_ANTEIL:.0%} bindet:\")\n", " print(f\" Depot : {(anteile_depot > MAX_ANTEIL - 1e-4).sum()} von \"\n", " f\"{len(anteile_depot)} Posten am Anschlag\")\n", " print(f\" Einkauf : {(anteile_einkauf > MAX_ANTEIL - 1e-4).sum()} von \"\n", " f\"{len(anteile_einkauf)} Posten am Anschlag\")\n", " print()\n", " print(\" Im Einkauf zieht es die Loesung an den sicheren Lokalzulieferer,\")\n", " print(\" bis die Grenze sie stoppt. Im Depot nicht - dort verhindert die\")\n", " print(\" Mindestrendite, dass alles in Anleihen wandert. Zwei verschiedene\")\n", " print(\" Bremsen also, und sie stehen an verschiedenen Stellen des Modells:\")\n", " print(\" einmal in einer Nebenbedingung ueber die Anteile, einmal in einer\")\n", " print(\" ueber den Ertrag. Wer eine Struktur uebertraegt, uebertraegt eben\")\n", " print(\" nicht automatisch mit, WELCHE Bedingung am Ende bindet.\")\n", "\n", " # --- Teil 3: Der ehrliche Teil ---------------------------------------\n", " print(\"\\n\" + \"=\" * 80)\n", " print(\" WAS NICHT HINUEBERREICHT\")\n", " print(\"=\" * 80)\n", " print(\"Die Struktur traegt. Drei Unterschiede tragen NICHT mit, und wer sie\")\n", " print(\"uebersieht, macht aus einer nuetzlichen Analogie einen Fehler:\")\n", " print()\n", " print(\"1. WOHER DIE ZAHLEN KOMMEN.\")\n", " print(\" In der Werkstatt ist der Verbrauch je Tisch gemessen - drei\")\n", " print(\" Montagestunden sind drei Montagestunden. Im Depot ist die\")\n", " print(\" erwartete Rendite GESCHAETZT, und zwar mit einem Fehler, der\")\n", " print(\" groesser sein kann als die Unterschiede zwischen den Anlagen\")\n", " print(\" (Kapitel Markowitz, Renditeschaetzung_Falle.py). Dieselbe\")\n", " print(\" Optimierung ist im einen Fall Planung und im anderen\")\n", " print(\" Fehlerverstaerkung.\")\n", " print()\n", " print(\"2. OB DIE VERGANGENHEIT ETWAS UEBER DIE ZUKUNFT SAGT.\")\n", " print(\" Lieferzeiten haben physikalische Ursachen: Entfernung, Zoll,\")\n", " print(\" Kapazitaet. Sie aendern sich langsam und nachvollziehbar.\")\n", " print(\" Kursrenditen entstehen aus dem Verhalten von Marktteilnehmern,\")\n", " print(\" die selbst auf Modelle reagieren - dort verschwindet ein\")\n", " print(\" erkanntes Muster oft genau deshalb, weil es erkannt wurde\")\n", " print(\" (Kapitel Handelsmaschine, Data_Snooping.py).\")\n", " print()\n", " print(\"3. OB TEILBARKEIT ERLAUBT IST.\")\n", " print(\" 37,4 % eines Aktienfonds sind ein normaler Auftrag. 37,4 % eines\")\n", " print(\" Lieferanten sind es nicht - Vertraege, Mindestabnahmen und\")\n", " print(\" Ruestzeiten machen Einkaufsentscheidungen ganzzahlig. Genau\")\n", " print(\" deshalb ist Teil II voller MILP und Teil V fast frei davon.\")\n", " print()\n", " print(\"Die Bruecke traegt also die MODELLE, nicht die Annahmen. Wer sie\")\n", " print(\"benutzt, spart sich das Lernen der Methoden - nicht das Nachdenken\")\n", " print(\"ueber die Daten.\")\n", " print(\"=\" * 80)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.11" } }, "nbformat": 4, "nbformat_minor": 5 }