operations_research/bilder_04/erzeuge_lastabwurf_verteilung.py

#!/usr/bin/env python3

# erzeuge_lastabwurf_verteilung.py
"""
Erzeugt das Verteilungsdiagramm zum Supply-Chain-Kapitel:

    bilder_04/kap_supplychain_lastabwurf.svg

Das vorhandene Bild kap_supplychain_wirkung.svg zerlegt die Kostendifferenz der
beiden Plaene nach Kostenarten. Dieses hier zeigt etwas anderes: die 40
Szenarien einzeln.

Warum das noetig ist: Ein Mittelwert von 874 870 EUR gegen 351 356 EUR sagt
nicht, WIE die Differenz zustande kommt. Die Verteilung sagt es. Der
Erwartungswert-Plan ist in den ruhigen Szenarien tatsaechlich billiger - und
faellt in den windarmen Szenarien so weit durch, dass der Mittelwert kippt. Das
ist die Gestalt, die eine einzelne Kennzahl verbirgt.

Die Szenarien mit Lastabwurf sind eigens markiert; ihre Zahl muss die im
Kapitel abgedruckte sein (28 von 40 gegen 0 von 40), sonst bricht das Skript
ab.

RECHENZEIT: Das Skript loest 2 x 40 Fahrplaene zuzueglich zweier
Commitment-Probleme. Das dauert einige Minuten - anders ist die Verteilung
nicht ehrlich zu bekommen.

Aufruf (aus dem Repository-Wurzelverzeichnis):
    python3 bilder_04/erzeuge_lastabwurf_verteilung.py

Benoetigt: numpy, matplotlib, ortools
"""

from __future__ import annotations

import os
import sys

import numpy as np

HIER = os.path.dirname(os.path.abspath(__file__))
BASIS = os.path.dirname(HIER)
sys.path.insert(0, HIER)
sys.path.insert(0, os.path.join(
    BASIS, "Operations_Research_mit_Python_Version_04_Programme"))

from stil_04 import FARBEN, neue_figur, speichere          # noqa: E402
from Kraftwerkseinsatz import (SZENARIEN, WIND_ERWARTUNG,  # noqa: E402
                               bewerte, erzeuge_windszenarien, plane)

# Aus der Befundtabelle des Kapitels - nur zum Gegenpruefen.
ABWURF_IM_BUCH = {"Erwartungswert": 28, "alle Szenarien": 0}


def pruefe_abwurf(zaehlung: dict[str, int]) -> None:
    for name, erwartet in ABWURF_IM_BUCH.items():
        if zaehlung[name] != erwartet:
            raise SystemExit(
                f"Lastabwurf-Bild: Der Plan '{name}' wirft in "
                f"{zaehlung[name]} von {SZENARIEN} Szenarien Last ab, die "
                f"Befundtabelle des Kapitels nennt {erwartet}.")


def zeichne(kosten: dict, fehlt: dict) -> None:
    figur, achse = neue_figur(breite=8.2, hoehe=4.6)

    for name, farbe in (("Erwartungswert", FARBEN["fehler"]),
                        ("alle Szenarien", FARBEN["gut"])):
        werte = np.sort(kosten[name]) / 1000.0
        rang = np.arange(1, len(werte) + 1)
        achse.step(werte, rang, where="post", color=farbe, linewidth=2.2,
                   zorder=4,
                   label=f"Plan „{name}“ — Mittel "
                         f"{kosten[name].mean() / 1000:,.0f} Tsd. €"
                         .replace(",", "."))

        betroffen = int((fehlt[name] > 1e-6).sum())
        if betroffen:
            mit_abwurf = np.sort(kosten[name][fehlt[name] > 1e-6]) / 1000.0
            achse.scatter(mit_abwurf,
                          np.searchsorted(werte, mit_abwurf, side="right"),
                          s=26, color=farbe, zorder=5, alpha=0.85)

    achse.set_xlabel("Kosten des Tages in Tausend €")
    achse.set_ylabel("Szenarien mit höchstens diesen Kosten")
    achse.set_title(f"Dieselben {SZENARIEN} Windszenarien, zwei Pläne — "
                    f"Punkte: Szenarien mit Lastabwurf",
                    fontsize=10.5, color=FARBEN["text"])
    achse.set_ylim(0, SZENARIEN + 1)
    achse.legend(fontsize=9, loc="lower right", frameon=False)

    achse.annotate(
        f"Der Erwartungswert-Plan ist links günstiger —\n"
        f"und wirft in {int((fehlt['Erwartungswert'] > 1e-6).sum())} von "
        f"{SZENARIEN} Szenarien Last ab.\n"
        f"Der zweistufige Plan in keinem.",
        xy=(0.03, 0.95), xycoords="axes fraction", fontsize=9,
        color=FARBEN["text"], va="top")

    speichere(figur, "kap_supplychain_lastabwurf")


if __name__ == "__main__":
    # erzeuge_windszenarien() liefert (wind, ist_flaute), plane() liefert
    # (plan, dauer, optimal) - beide Rueckgaben werden hier nur teilweise
    # gebraucht.
    wind, _ist_flaute = erzeuge_windszenarien()
    gleich = np.full(SZENARIEN, 1.0 / SZENARIEN)

    kosten, fehlt, zaehlung = {}, {}, {}
    for name, plan in (
            ("Erwartungswert", plane(WIND_ERWARTUNG.reshape(1, -1), [1.0])[0]),
            ("alle Szenarien", plane(wind, gleich)[0])):
        kosten[name], fehlt[name] = bewerte(plan, wind)
        zaehlung[name] = int((fehlt[name] > 1e-6).sum())
        print(f"  Plan „{name}“: Mittel {kosten[name].mean():>12,.0f} EUR, "
              f"schlimmster Tag {kosten[name].max():>12,.0f} EUR, "
              f"Lastabwurf in {zaehlung[name]:>2} von {SZENARIEN}")

    pruefe_abwurf(zaehlung)
    zeichne(kosten, fehlt)