operations_research/Notebooks_04/prognose.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kapitel 15: Predict-then-Optimize{idx:Predict-then-Optimize} — die bessere Prognose, die schlechtere Entscheidung\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",
    "`Predict_then_Optimize.py`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "\n",
    "# Predict_then_Optimize.py\n",
    "\"\"\"\n",
    "Kapitel Prognose: Die bessere Prognose trifft die schlechtere Entscheidung.\n",
    "\n",
    "Eine Baeckerei muss jeden Abend entscheiden, wie viel sie fuer den naechsten Tag\n",
    "ansetzt. Zu wenig kostet die Marge des entgangenen Verkaufs, zu viel kostet den\n",
    "Einkaufspreis der Retoure. Das ist das Newsvendor-Problem aus dem Kapitel\n",
    "Unsicherheit - nur dass die Nachfrage diesmal nicht aus einer Verteilung kommt,\n",
    "sondern PROGNOSTIZIERT werden muss: aus Wochentag, Temperatur und Aktionstagen.\n",
    "\n",
    "Damit zerfaellt die Aufgabe in zwei Schritte, und genau an der Naht entsteht der\n",
    "Fehler, um den es hier geht:\n",
    "\n",
    "    PREDICT     ein Modell schaetzt die Nachfrage\n",
    "    OPTIMIZE    daraus wird eine Bestellmenge\n",
    "\n",
    "Der Prognostiker optimiert seinen Modellfehler, meist den MSE. Der Planer traegt\n",
    "die Kosten. Beide messen etwas anderes - und die beiden Masse widersprechen\n",
    "einander. Das Programm zeigt:\n",
    "\n",
    "  1. Vier Verfahren, verglichen nach MSE UND nach Entscheidungskosten. Das\n",
    "     Verfahren mit dem BESTEN MSE hat die HOECHSTEN Kosten.\n",
    "  2. Warum ein pauschaler Sicherheitszuschlag zu kurz greift - die Streuung der\n",
    "     Nachfrage haengt selbst von den Merkmalen ab.\n",
    "  3. Eine Messfalle, in die der Autor dieses Programms zuerst selbst getappt\n",
    "     ist: Bei kurzen Testzeitraeumen ist der MSE-Vergleich nicht stabil.\n",
    "\n",
    "Benoetigt: numpy, scipy, scikit-learn\n",
    "\"\"\"\n",
    "\n",
    "from __future__ import annotations\n",
    "\n",
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "from sklearn.linear_model import LinearRegression, QuantileRegressor\n",
    "\n",
    "VERKAUFSPREIS = 9.0\n",
    "EINKAUFSPREIS = 3.0\n",
    "KOSTEN_FEHLMENGE = VERKAUFSPREIS - EINKAUFSPREIS      # entgangene Marge: 6 EUR\n",
    "KOSTEN_UEBERHANG = EINKAUFSPREIS                      # Retoure:           3 EUR\n",
    "KRITISCHES_VERHAELTNIS = KOSTEN_FEHLMENGE / (KOSTEN_FEHLMENGE + KOSTEN_UEBERHANG)\n",
    "\n",
    "TAGE = 5000               # Simulation, siehe Hinweis unten\n",
    "TRAINING = 1000\n",
    "SAAT = 11\n",
    "\n",
    "\n",
    "def erzeuge_daten(tage: int = TAGE, saat: int = SAAT):\n",
    "    \"\"\"Taegliche Nachfrage mit Wochentag, Temperatur und Aktionstagen.\n",
    "\n",
    "    Die entscheidende Eigenschaft steckt in 'streuung': An Aktionstagen ist die\n",
    "    Nachfrage nicht nur hoeher, sondern auch viel UNSICHERER. Solche\n",
    "    heteroskedastischen Daten sind der Normalfall - und der Grund, warum ein\n",
    "    pauschaler Sicherheitszuschlag nicht genuegt.\n",
    "    \"\"\"\n",
    "    rng = np.random.default_rng(saat)\n",
    "    wochentag = np.arange(tage) % 7\n",
    "    temperatur = (12 + 10 * np.sin(2 * np.pi * np.arange(tage) / 365)\n",
    "                  + rng.normal(0, 3, tage))\n",
    "    aktion = (rng.random(tage) < 0.15).astype(float)\n",
    "\n",
    "    merkmale = np.column_stack([np.eye(7)[wochentag][:, 1:], temperatur, aktion])\n",
    "    erwartung = (120\n",
    "                 + np.eye(7)[wochentag] @ np.array([0, 10, 12, 14, 18, 35, -40])\n",
    "                 + 1.8 * temperatur + 45 * aktion)\n",
    "    streuung = 8 + 22 * aktion\n",
    "    nachfrage = np.maximum(0.0, erwartung + rng.normal(0, 1, tage) * streuung)\n",
    "    return merkmale, nachfrage, aktion\n",
    "\n",
    "\n",
    "def tageskosten(bestellung: np.ndarray, nachfrage: np.ndarray) -> float:\n",
    "    \"\"\"Die Zahl, auf die es ankommt - und die kein Prognosemass kennt.\"\"\"\n",
    "    fehlmenge = np.maximum(0.0, nachfrage - bestellung)\n",
    "    ueberhang = np.maximum(0.0, bestellung - nachfrage)\n",
    "    return float((KOSTEN_FEHLMENGE * fehlmenge\n",
    "                  + KOSTEN_UEBERHANG * ueberhang).mean())\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    merkmale, nachfrage, aktion = erzeuge_daten()\n",
    "    lernen = slice(0, TRAINING)\n",
    "    pruefen = slice(TRAINING, TAGE)\n",
    "\n",
    "    print(\"=\" * 84)\n",
    "    print(\"  DIE BESSERE PROGNOSE TRIFFT DIE SCHLECHTERE ENTSCHEIDUNG\")\n",
    "    print(\"=\" * 84)\n",
    "    print(f\"Verkaufspreis {VERKAUFSPREIS:.0f} EUR, Einkauf {EINKAUFSPREIS:.0f} EUR.\")\n",
    "    print(f\"Fehlmenge kostet {KOSTEN_FEHLMENGE:.0f} EUR, Ueberhang \"\n",
    "          f\"{KOSTEN_UEBERHANG:.0f} EUR je Stueck.\")\n",
    "    print(f\"Kritisches Verhaeltnis: {KRITISCHES_VERHAELTNIS:.3f} - der Planer sollte \"\n",
    "          f\"also das\")\n",
    "    print(f\"{KRITISCHES_VERHAELTNIS:.1%}-Quantil der Nachfrage bestellen, nicht ihren \"\n",
    "          f\"Erwartungswert.\")\n",
    "    print(f\"\\nTraining: Tag 1 bis {TRAINING}. Bewertung: die restlichen \"\n",
    "          f\"{TAGE - TRAINING} Tage.\\n\")\n",
    "\n",
    "    # --- Die vier Verfahren ----------------------------------------------\n",
    "    kleinste_quadrate = LinearRegression().fit(merkmale[lernen], nachfrage[lernen])\n",
    "    punktprognose = kleinste_quadrate.predict(merkmale[pruefen])\n",
    "\n",
    "    restfehler = nachfrage[lernen] - kleinste_quadrate.predict(merkmale[lernen])\n",
    "    pauschalzuschlag = norm.ppf(KRITISCHES_VERHAELTNIS) * restfehler.std()\n",
    "\n",
    "    # Ein Zuschlag, der nicht aus der Normalverteilung kommt, sondern direkt\n",
    "    # auf den Trainingsdaten die Kosten minimiert.\n",
    "    kandidaten = np.linspace(-10.0, 30.0, 401)\n",
    "    trainingsprognose = kleinste_quadrate.predict(merkmale[lernen])\n",
    "    kostenzuschlag = float(kandidaten[np.argmin(\n",
    "        [tageskosten(trainingsprognose + z, nachfrage[lernen]) for z in kandidaten])])\n",
    "\n",
    "    # Und das Verfahren, das von vornherein das richtige Quantil schaetzt.\n",
    "    quantilmodell = QuantileRegressor(quantile=KRITISCHES_VERHAELTNIS,\n",
    "                                      alpha=0.0, solver=\"highs\")\n",
    "    quantilmodell.fit(merkmale[lernen], nachfrage[lernen])\n",
    "    quantilprognose = quantilmodell.predict(merkmale[pruefen])\n",
    "\n",
    "    verfahren = [\n",
    "        (\"bestelle die Punktprognose\", punktprognose, punktprognose),\n",
    "        (\"+ Zuschlag aus der Normalverteilung\",\n",
    "         punktprognose, punktprognose + pauschalzuschlag),\n",
    "        (\"+ Zuschlag auf Kosten trainiert\",\n",
    "         punktprognose, punktprognose + kostenzuschlag),\n",
    "        (\"Quantilregression aufs kritische Quantil\",\n",
    "         quantilprognose, quantilprognose),\n",
    "    ]\n",
    "\n",
    "    print(f\"  {'Verfahren':<42} {'MSE':>9} {'Kosten/Tag':>12} {'gegen Zeile 1':>14}\")\n",
    "    print(\"  \" + \"-\" * 80)\n",
    "    ergebnisse = {}\n",
    "    for name, prognose, bestellung in verfahren:\n",
    "        mse = float(((prognose - nachfrage[pruefen]) ** 2).mean())\n",
    "        kosten = tageskosten(bestellung, nachfrage[pruefen])\n",
    "        ergebnisse[name] = (mse, kosten)\n",
    "        basis = ergebnisse[verfahren[0][0]][1]\n",
    "        vergleich = \"\" if name == verfahren[0][0] else f\"{(kosten - basis) / basis:+13.1%}\"\n",
    "        print(f\"  {name:<42} {mse:>9.1f} {kosten:>10.2f} EUR {vergleich:>14}\")\n",
    "\n",
    "    bester_mse = min(ergebnisse, key=lambda k: ergebnisse[k][0])\n",
    "    beste_kosten = min(ergebnisse, key=lambda k: ergebnisse[k][1])\n",
    "    print(f\"\\n  bester MSE:      {bester_mse}\")\n",
    "    print(f\"  beste Kosten:    {beste_kosten}\")\n",
    "    print(f\"\\n  Das Verfahren mit dem besten MSE hat die HOECHSTEN Kosten, und das\")\n",
    "    print(f\"  Verfahren mit den besten Kosten hat einen um \"\n",
    "          f\"{(ergebnisse[beste_kosten][0] / ergebnisse[bester_mse][0] - 1):.0%} SCHLECHTEREN MSE.\")\n",
    "    print(\"  Wer Prognosemodelle nach MSE auswaehlt, waehlt hier das falsche.\")\n",
    "\n",
    "    # --- Warum der pauschale Zuschlag zu kurz greift ---------------------\n",
    "    print(\"\\n\" + \"-\" * 84)\n",
    "    print(\"Warum ein pauschaler Zuschlag nicht genuegt\\n\")\n",
    "    ist_aktion = aktion[pruefen] > 0.5\n",
    "    print(f\"  {'Verfahren':<42} {'normale Tage':>14} {'Aktionstage':>14}\")\n",
    "    print(\"  \" + \"-\" * 74)\n",
    "    for name, _, bestellung in verfahren:\n",
    "        normal = tageskosten(bestellung[~ist_aktion], nachfrage[pruefen][~ist_aktion])\n",
    "        aktionstag = tageskosten(bestellung[ist_aktion], nachfrage[pruefen][ist_aktion])\n",
    "        print(f\"  {name:<42} {normal:>10.2f} EUR {aktionstag:>10.2f} EUR\")\n",
    "\n",
    "    # Nachgerechnet statt behauptet: Welcher Zuschlag waere je Tagesart richtig?\n",
    "    aktion_training = aktion[lernen] > 0.5\n",
    "    z = norm.ppf(KRITISCHES_VERHAELTNIS)\n",
    "    richtig_normal = z * restfehler[~aktion_training].std()\n",
    "    richtig_aktion = z * restfehler[aktion_training].std()\n",
    "    print(f\"\\n  Der pauschale Zuschlag betraegt {pauschalzuschlag:.1f} Stueck. Aus den\")\n",
    "    print(f\"  Trainingsresten getrennt nach Tagesart waere richtig:\")\n",
    "    print(f\"    normale Tage : {richtig_normal:5.1f} Stueck\")\n",
    "    print(f\"    Aktionstage  : {richtig_aktion:5.1f} Stueck\")\n",
    "    print(f\"  Ein Zuschlag fuer alle Tage kann nur einen Mittelweg treffen - hier\")\n",
    "    print(f\"  ist er an normalen Tagen {pauschalzuschlag / richtig_normal:.1f}-mal zu gross und an\")\n",
    "    print(f\"  Aktionstagen nur {pauschalzuschlag / richtig_aktion:.0%} dessen, was noetig waere.\")\n",
    "    print(f\"\\n  Die Quantilregression schaetzt das {KRITISCHES_VERHAELTNIS:.1%}-Quantil \"\n",
    "          f\"direkt aus den\")\n",
    "    print(f\"  Merkmalen und darf deshalb an verschiedenen Tagen verschieden weit\")\n",
    "    print(f\"  ueber dem Erwartungswert liegen. Genau das ist der Unterschied\")\n",
    "    print(f\"  zwischen 'ein Modell und danach eine Formel' und 'ein Modell, das\")\n",
    "    print(f\"  weiss, wofuer es gebraucht wird'.\")\n",
    "\n",
    "    # --- Die Messfalle ---------------------------------------------------\n",
    "    print(\"\\n\" + \"-\" * 84)\n",
    "    print(\"Eine Messfalle, in die der Autor zuerst selbst getappt ist\\n\")\n",
    "    print(\"  Der erste Entwurf dieses Programms bewertete auf 230 Testtagen - ein\")\n",
    "    print(\"  realistischer Zeitraum. Dort hatte die Quantilregression den BESSEREN\")\n",
    "    print(\"  MSE, und die ganze Aussage des Kapitels stand auf dem Kopf.\")\n",
    "    print(\"\\n  Wie oft das passiert, laesst sich ausmessen:\\n\")\n",
    "    rng = np.random.default_rng(0)\n",
    "    print(f\"  {'Testfenster':>14} {'QR sieht MSE-besser aus':>26}\")\n",
    "    print(\"  \" + \"-\" * 42)\n",
    "    for fenster in (180, 365, 730, 2000):\n",
    "        treffer = 0\n",
    "        versuche = 400\n",
    "        for _ in range(versuche):\n",
    "            start = int(rng.integers(TRAINING, TAGE - fenster))\n",
    "            ausschnitt = slice(start, start + fenster)\n",
    "            mse_punkt = ((kleinste_quadrate.predict(merkmale[ausschnitt])\n",
    "                          - nachfrage[ausschnitt]) ** 2).mean()\n",
    "            mse_quantil = ((quantilmodell.predict(merkmale[ausschnitt])\n",
    "                            - nachfrage[ausschnitt]) ** 2).mean()\n",
    "            treffer += mse_quantil < mse_punkt\n",
    "        print(f\"  {fenster:>10} Tage {treffer / versuche:>24.1%}\")\n",
    "\n",
    "    print(\"\\n  Bei einem halben Jahr Testdaten sieht das schlechtere Modell in gut\")\n",
    "    print(\"  jedem zehnten Fall besser aus. Das ist keine grosse Zahl - aber wer\")\n",
    "    print(\"  EINMAL misst, hat genau eine Ziehung aus dieser Verteilung.\")\n",
    "    print(\"\\n  Die Lehre ist nicht 'nimm 4.000 Testtage' - die hat niemand. Sie\")\n",
    "    print(\"  lautet: Ein Kennzahlenvergleich ohne Angabe seiner Streuung ist keine\")\n",
    "    print(\"  Aussage. Bei kurzen Zeitraeumen gehoert eine Kreuzvalidierung dazu.\")\n",
    "\n",
    "    print(\"\\n\" + \"=\" * 84)\n",
    "    print(\"  WAS MAN DARAUS MITNIMMT\")\n",
    "    print(\"=\" * 84)\n",
    "    print(\"Der Prognostiker optimiert den MSE, der Planer traegt die Kosten - und\")\n",
    "    print(\"die beiden Masse zeigen hier in verschiedene Richtungen. Drei Saetze:\")\n",
    "    print()\n",
    "    print(\"  1. Sagen Sie nicht den Erwartungswert vorher, sondern die Groesse, die\")\n",
    "    print(\"     in die Entscheidung eingeht. Beim Newsvendor ist das das kritische\")\n",
    "    print(\"     Quantil - und das kann man direkt schaetzen.\")\n",
    "    print(\"  2. Bewerten Sie Prognosemodelle an den ENTSCHEIDUNGSKOSTEN. Die sind\")\n",
    "    print(\"     in Euro und damit vergleichbar; ein MSE ist es nicht.\")\n",
    "    print(\"  3. Ein pauschaler Sicherheitszuschlag ist besser als nichts und\")\n",
    "    print(\"     schlechter als ein Modell, das die Unsicherheit selbst aus den\")\n",
    "    print(\"     Merkmalen liest.\")\n",
    "    print()\n",
    "    print(\"Der naechste Schritt - Prognosemodelle so zu trainieren, dass sie die\")\n",
    "    print(\"Entscheidungskosten direkt minimieren (Smart Predict-then-Optimize,\")\n",
    "    print(\"differenzierbare Optimierungsschichten) - ist Forschungsstand und\")\n",
    "    print(\"erfordert Bibliotheken wie cvxpylayers. Die dritte Zeile der Tabelle\")\n",
    "    print(\"oben ist seine einfachste denkbare Form: ein einziger Parameter, auf\")\n",
    "    print(\"Kosten statt auf Fehler trainiert.\")\n",
    "    print(\"=\" * 84)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}