operations_research/bilder_04/erzeuge_big_m_vergleich.py

#!/usr/bin/env python3

# erzeuge_big_m_vergleich.py
"""
Erzeugt das Big-M-Diagramm zum MILP-Kapitel:

    bilder_04/kap_milp_big_m.svg

Zwei Bilder derselben Bedingung x <= M*y mit y in [0,1], einmal mit straffem
und einmal mit ueberdimensioniertem M. Gezeigt wird die LP-RELAXATION - also
das, was der Solver sieht, bevor er y auf 0 oder 1 zwingt.

Der Punkt: Beide Modelle sind fuer y in {0,1} GLEICHWERTIG. Sie beschreiben
denselben zulaessigen Bereich, sobald y ganzzahlig ist. Verschieden ist nur die
Relaxation - und die entscheidet, wie gut die Schranke ist und wie lange der
Solver sucht. Bei zu grossem M kommt der Trickle Flow dazu: Ein y knapp ueber
null gilt dem Solver als geschlossen und laesst trotzdem grossen Fluss durch.

Die Trickle-Flow-Zahlen in der Anmerkung stammen aus der Loesung im Anhang des
Buchs (M = 7,2e9, y = 4e-8, rund 289 Einheiten Fluss). Die beiden M-Werte der
Geometrie sind Anschauungswerte; das Skript rechnet aus, ab welchem y die
Relaxation die volle Kapazitaet freigibt, statt eine Zahl zu behaupten. Die
Flaeche des Relaxationsgebiets waere der falsche Massstab - sie ist wegen der
Kapazitaetsgrenze in beiden Faellen aehnlich gross.

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

Benoetigt: numpy, matplotlib
"""

from __future__ import annotations

import os
import sys

import numpy as np

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
from stil_04 import FARBEN, speichere      # noqa: E402
import matplotlib.pyplot as plt            # noqa: E402

KAPAZITAET = 40.0        # der tatsaechlich moegliche Hoechstwert von x
M_STRAFF = 40.0          # so klein wie moeglich: M = Kapazitaet
M_GROSS = 500.0          # "sicherheitshalber gross" - der haeufige Fehler

# Aus der Loesung im Anhang: ein verschwindend kleines y laesst grossen Fluss
# durch, weil es mit einem riesigen M multipliziert wird.
TRICKLE_M = 7.2e9
TRICKLE_Y = 4e-8


def schwelle(m: float) -> float:
    """Kleinstes y, bei dem die Relaxation die volle Kapazitaet freigibt.

    Das ist der aussagekraeftige Vergleich - nicht die Flaeche des
    Relaxationsgebiets. Die ist wegen der Kapazitaetsgrenze in beiden Faellen
    aehnlich gross; der Unterschied liegt darin, wie WENIG y genuegt, um alles
    durchzulassen."""
    return min(1.0, KAPAZITAET / m)


def zeichne(faktor: float) -> None:
    figur, achsen = plt.subplots(1, 2, figsize=(9.2, 4.0), sharey=True)
    y = np.linspace(0, 1, 300)

    for achse, m, titel, farbe in (
            (achsen[0], M_STRAFF, f"straff: $M = {M_STRAFF:.0f}$",
             FARBEN["gut"]),
            (achsen[1], M_GROSS, f"zu groß: $M = {M_GROSS:.0f}$",
             FARBEN["fehler"])):
        obergrenze = np.minimum(m * y, KAPAZITAET)
        achse.fill_between(y, 0, obergrenze, color=farbe, alpha=0.16, zorder=2)
        achse.plot(y, obergrenze, "-", color=farbe, linewidth=2.2, zorder=4,
                   label=r"$x \le M\,y$ und $x \le$ Kapazität")
        achse.axhline(KAPAZITAET, color=FARBEN["linie"], linewidth=1.0,
                      linestyle="--", zorder=3)
        achse.text(0.02, KAPAZITAET + 1.5, "Kapazität", fontsize=8.5,
                   color=FARBEN["gedaempft"])

        # Die beiden Punkte, die das ganzzahlige Modell zulaesst.
        achse.plot([0, 1], [0, KAPAZITAET], "o", color=FARBEN["haupt"],
                   markersize=8, zorder=5)
        achse.annotate(r"$y=0$: zu", xy=(0, 0), xytext=(8, 10),
                       textcoords="offset points", fontsize=8.5,
                       color=FARBEN["haupt"])
        achse.annotate(r"$y=1$: offen", xy=(1, KAPAZITAET), xytext=(-64, -18),
                       textcoords="offset points", fontsize=8.5,
                       color=FARBEN["haupt"])

        achse.set_xlim(0, 1.02)
        achse.set_ylim(0, KAPAZITAET * 1.25)
        achse.set_xlabel(r"Schaltvariable $y \in [0,1]$ (relaxiert)")
        achse.set_title(titel, fontsize=10.5, color=farbe)

    achsen[0].set_ylabel(r"Menge $x$")
    achsen[1].annotate(
        f"Schon $y = {schwelle(M_GROSS):.2f}$ gibt die volle\n"
        f"Kapazität frei — links braucht es dafür $y = 1$.\n"
        f"Der Solver darf das Lager also fast\n"
        f"geschlossen lassen und trotzdem alles liefern.",
        xy=(0.16, KAPAZITAET * 0.42), fontsize=9, color=FARBEN["fehler"],
        ha="left")
    achsen[1].plot([schwelle(M_GROSS)], [KAPAZITAET], "o",
                   color=FARBEN["fehler"], markersize=8, zorder=6)

    figur.suptitle("Dieselbe Bedingung, zwei Zahlen für $M$ — gleich für "
                   "$y \\in \\{0,1\\}$, verschieden für den Solver",
                   fontsize=10.5, color=FARBEN["text"], y=1.02)
    speichere(figur, "kap_milp_big_m")


if __name__ == "__main__":
    faktor = schwelle(M_STRAFF) / schwelle(M_GROSS)
    print(f"  volle Kapazität ab y = {schwelle(M_STRAFF):.3f} (M = "
          f"{M_STRAFF:.0f})")
    print(f"  volle Kapazität ab y = {schwelle(M_GROSS):.3f} (M = "
          f"{M_GROSS:.0f})  → {faktor:.0f}-mal weniger y genügt")
    print(f"  Trickle Flow im Anhang: y = {TRICKLE_Y:g} gilt als 0, "
          f"laesst aber M*y = {TRICKLE_M * TRICKLE_Y:.0f} Einheiten durch.")
    zeichne(faktor)