operations_research/OR_HTML_04/programme/Spaltengenerierung.py

#!/usr/bin/env python3

# Spaltengenerierung.py
"""
Kapitel Dekomposition: Das Modell umbauen, statt die Loesung zu raten.

Eine Papierfabrik liefert Rollen von 5.600 mm Breite. Die Kunden bestellen
schmalere Breiten, und aus jeder Mutterrolle werden mehrere davon geschnitten.
Gesucht ist der Schnittplan, der mit den wenigsten Mutterrollen auskommt.

Der naheliegende Modellansatz - "welches Stueck kommt auf welche Rolle?" -
fuehrt in eine Sackgasse: Er braucht eine Binaervariable je Stueck-und-Rolle,
und weil alle Rollen gleich sind, ist er hochgradig symmetrisch. Der Solver
probiert dieselbe Loesung in tausend Umbenennungen durch.

Die Spaltengenerierung dreht die Frage um:

    Nicht "welches Stueck auf welche Rolle",
    sondern "wie oft schneide ich welches MUSTER".

Damit verschwindet die Symmetrie - aber es entsteht ein neues Problem: Die Zahl
der moeglichen Muster waechst kombinatorisch. Man kann sie nicht aufschreiben.
Der Trick besteht darin, sie auch nicht aufzuschreiben, sondern nur die wenigen
zu erzeugen, die tatsaechlich gebraucht werden - und zwar mit Hilfe der
Schattenpreise aus dem Kapitel LP.

Das Programm zeigt vier Dinge:

  1. Die beiden Formulierungen desselben Problems und ihre Groesse.
  2. Wie schnell die Musterzahl waechst - und wie wenige davon genuegen.
  3. Die Schleife selbst: Master-LP, Dualwerte, Rucksack, neue Spalte.
  4. Wann sich der Aufwand lohnt und wann nicht. Die Antwort haengt an einer
     einzigen Kennzahl der Instanz, und sie erklaert nebenbei, warum eine
     Faustregel bei manchen Zuschnittproblemen schon optimal ist.

Benoetigt: numpy, scipy
"""

from __future__ import annotations

import numpy as np
from scipy.optimize import linprog

ROLLENBREITE = 5600

# Zwei Instanzen mit demselben Bedarf, aber verschiedenen Breiten. Der
# Unterschied ist die eine Kennzahl, um die es in Teil 4 geht.
BEDARF = [22, 25, 12, 14, 18, 18, 20, 10, 12, 14, 16, 18, 20]
BREIT = [1380, 1520, 1560, 1710, 1820, 1880, 1930, 2000, 2050, 2100, 2140,
         2150, 2200]
SCHMAL = [380, 520, 560, 710, 820, 880, 930, 1000, 1050, 1100, 1140,
          1150, 1200]


# --- Die Musterzahl ---------------------------------------------------------

def zaehle_muster(breiten, rollenbreite: int = ROLLENBREITE) -> int:
    """Wie viele zulaessige Schnittmuster gibt es ueberhaupt?

    Nur zum Zeigen - in einem produktiven Programm wuerde man das nie
    ausrechnen, weil die Zahl bei realistischen Instanzen jede Vorstellung
    sprengt. Genau das ist der Punkt.
    """
    speicher: dict[tuple[int, int], int] = {}

    def ab(i: int, rest: int) -> int:
        if i == len(breiten):
            return 1
        if (i, rest) in speicher:
            return speicher[(i, rest)]
        summe = sum(ab(i + 1, rest - anzahl * breiten[i])
                    for anzahl in range(rest // breiten[i] + 1))
        speicher[(i, rest)] = summe
        return summe

    return ab(0, rollenbreite)


# --- Die beiden Bausteine der Schleife --------------------------------------

def loese_master(muster, bedarf, ganzzahlig: bool = False):
    """Restringiertes Master: Wie oft schneide ich jedes bekannte Muster?

        min  sum_p x_p           unter   sum_p a_ip * x_p >= bedarf_i

    Als LP gerechnet liefert es zusaetzlich die DUALWERTE - und die sind der
    eigentliche Ertrag: Der Dualwert zu Breite i sagt, was ein zusaetzliches
    Stueck dieser Breite an Rollen kostet. Es ist derselbe Schattenpreis wie
    im Kapitel LP, nur dass er hier nicht berichtet, sondern weiterverarbeitet
    wird.
    """
    matrix = -np.array(muster, dtype=float).T        # >= wird zu <= mit Minus
    ergebnis = linprog(np.ones(len(muster)),
                       A_ub=matrix, b_ub=-np.array(bedarf, dtype=float),
                       bounds=(0, None),
                       integrality=(1 if ganzzahlig else 0), method="highs")
    if not ergebnis.success:
        raise SystemExit(f"Master nicht loesbar: {ergebnis.message}")
    dual = None if ganzzahlig else -ergebnis.ineqlin.marginals
    return ergebnis.fun, ergebnis.x, dual


def bestes_neues_muster(dual, breiten, rollenbreite: int = ROLLENBREITE):
    """Pricing: Gibt es ein Muster, das sich noch lohnt?

    Gesucht ist das Muster mit dem groessten Gesamtwert zu den aktuellen
    Schattenpreisen - unter der Bedingung, dass es in die Rolle passt. Das ist
    ein RUCKSACKPROBLEM (Kapitel MILP, Rucksack.py): Preise sind der Nutzen,
    Breiten das Gewicht, die Rollenbreite die Kapazitaet. Hier mit dynamischer
    Programmierung geloest, weil die Breiten ganzzahlig sind.

    Bei Gleichstand wird bewusst der KLEINSTE Index bevorzugt (>-Vergleich mit
    Toleranz statt >=). Sonst haengt das erzeugte Muster von der Reihenfolge
    der Gleitkommaoperationen ab, und das Programm liefert von Lauf zu Lauf
    verschiedene Ausgaben.
    """
    wert = np.zeros(rollenbreite + 1)
    herkunft = [-1] * (rollenbreite + 1)
    for platz in range(rollenbreite + 1):
        for i, breite in enumerate(breiten):
            if breite <= platz and wert[platz - breite] + dual[i] > wert[platz] + 1e-9:
                wert[platz] = wert[platz - breite] + dual[i]
                herkunft[platz] = i

    muster = [0] * len(breiten)
    platz = rollenbreite
    while herkunft[platz] >= 0:
        i = herkunft[platz]
        muster[i] += 1
        platz -= breiten[i]
    return float(wert[rollenbreite]), muster


def spaltengenerierung(breiten, bedarf, rollenbreite: int = ROLLENBREITE):
    """Die Schleife. Startbasis: je ein Muster mit nur einer Breite.

    Abbruch, wenn kein Muster mehr einen Wert ueber 1 hat: Eine zusaetzliche
    Rolle kostet 1, also lohnt sich ein neues Muster nur, wenn es zu den
    aktuellen Preisen mehr als 1 wert ist. Das ist das Kriterium der
    reduzierten Kosten.
    """
    anzahl = len(breiten)
    muster = [[0] * anzahl for _ in range(anzahl)]
    for i in range(anzahl):
        muster[i][i] = rollenbreite // breiten[i]

    runden = 0
    while True:
        zielwert, _, dual = loese_master(muster, bedarf)
        wert, neu = bestes_neues_muster(dual, breiten, rollenbreite)
        runden += 1
        if wert <= 1 + 1e-6:
            return muster, zielwert, runden
        muster.append(neu)


def first_fit(breiten, bedarf, rollenbreite: int = ROLLENBREITE) -> int:
    """Die Faustregel: groesstes Stueck zuerst, auf die erste passende Rolle."""
    stuecke = sorted([b for b, menge in zip(breiten, bedarf)
                      for _ in range(menge)], reverse=True)
    rollen: list[list[int]] = []
    for stueck in stuecke:
        for rolle in rollen:
            if sum(rolle) + stueck <= rollenbreite:
                rolle.append(stueck)
                break
        else:
            rollen.append([stueck])
    return len(rollen)


def auswerten(name: str, breiten, bedarf) -> dict:
    muster, schranke, runden = spaltengenerierung(breiten, bedarf)
    ganz, _, _ = loese_master(muster, bedarf, ganzzahlig=True)
    return {"name": name, "muster_gesamt": zaehle_muster(breiten),
            "muster_erzeugt": len(muster), "runden": runden,
            "schranke": schranke, "ganzzahlig": int(round(ganz)),
            "first_fit": first_fit(breiten, bedarf),
            "stuecke_je_rolle": ROLLENBREITE / np.mean(breiten)}


if __name__ == "__main__":
    print("=" * 84)
    print("  SPALTENGENERIERUNG: NUR DIE MUSTER ERZEUGEN, DIE MAN BRAUCHT")
    print("=" * 84)
    print(f"Mutterrolle {ROLLENBREITE:,} mm, {len(BREIT)} bestellte Breiten, "
          f"{sum(BEDARF)} Zuschnitte.\n")

    # --- 1. Zwei Formulierungen ------------------------------------------
    stuecke = sum(BEDARF)
    # Wie viele Rollen braucht man hoechstens? Die Faustregel liefert eine
    # brauchbare Obergrenze - mehr Rollen als das wird niemand benoetigen.
    rollen_obergrenze = first_fit(BREIT, BEDARF)
    print("1. Zwei Modelle fuer dieselbe Aufgabe\n")
    print(f"  {'Formulierung':<34} {'Binaervariablen':>16} {'symmetrisch?':>16}")
    print("  " + "-" * 70)
    print(f"  {'Stueck -> Rolle (naheliegend)':<34} "
          f"{stuecke * rollen_obergrenze:>16,} {'ja':>16}")
    print(f"  {'wie oft welches Muster':<34} "
          f"{'eine je Muster':>16} {'nein':>16}")
    print(f"\n  Der naheliegende Ansatz braucht eine Variable je Stueck und Rolle:")
    print(f"  {stuecke} Zuschnitte x {rollen_obergrenze} Rollen (Obergrenze aus der Faustregel)")
    print(f"  = {stuecke * rollen_obergrenze:,} Binaervariablen.")
    print("\n  Schlimmer als die Zahl ist die Symmetrie: Alle Mutterrollen sind")
    print("  gleich, also beschreibt jede Loesung dieselbe Schnittvorschrift in")
    print("  unzaehligen Umbenennungen. Branch-and-Bound probiert sie einzeln")
    print("  durch und kommt nicht voran - der Suchbaum waechst, ohne dass sich")
    print("  die Schranke bewegt.")
    print("\n  Das Mustermodell hat je Muster genau eine Variable und kennt keine")
    print("  einzelnen Rollen mehr. Damit ist die Symmetrie weg. Sein Problem ist")
    print("  ein anderes - und zwar das folgende.")

    # --- 2. Die Musterzahl ------------------------------------------------
    print("\n" + "-" * 84)
    print("2. Warum man die Muster nicht aufschreiben kann\n")
    gross = auswerten("breite Zuschnitte (1.380-2.200 mm)", BREIT, BEDARF)
    schmal = auswerten("schmale Zuschnitte (380-1.200 mm)", SCHMAL, BEDARF)

    print(f"  {'Instanz':<36} {'Muster':>12} {'davon erzeugt':>14} {'Anteil':>9}")
    print("  " + "-" * 76)
    for fall in (gross, schmal):
        anteil = fall["muster_erzeugt"] / fall["muster_gesamt"]
        print(f"  {fall['name']:<36} {fall['muster_gesamt']:>12,} "
              f"{fall['muster_erzeugt']:>14} {anteil:>8.1%}")

    print(f"\n  Schon das Halbieren der Breiten laesst die Musterzahl von "
          f"{gross['muster_gesamt']:,} auf")
    print(f"  {schmal['muster_gesamt']:,} springen - Faktor "
          f"{schmal['muster_gesamt'] / gross['muster_gesamt']:.0f}. Bei einer echten "
          f"Papierfabrik mit")
    print("  vierzig Breiten und Millimeterschritten sind es mehr, als sich")
    print("  speichern liesse.")
    print(f"\n  Gebraucht werden davon {schmal['muster_erzeugt']} - "
          f"{schmal['muster_erzeugt'] / schmal['muster_gesamt']:.2%} der Gesamtzahl.")

    # --- 3. Die Schleife --------------------------------------------------
    print("\n" + "-" * 84)
    print("3. Die Schleife an der breiten Instanz\n")
    print(f"  Startbasis: {len(BREIT)} triviale Muster (je Rolle nur eine Breite)")
    print(f"  Runden bis kein Muster mehr lohnt: {gross['runden']}")
    print(f"  Muster am Ende: {gross['muster_erzeugt']}")
    print(f"  LP-Schranke: {gross['schranke']:.4f} Rollen")
    print(f"\n  Die LP-Schranke ist eine ZUSAGE: Weniger als "
          f"{np.ceil(gross['schranke'] - 1e-9):.0f} Rollen sind")
    print("  nicht moeglich - unabhaengig davon, wie clever man weiterschneidet.")
    print("  Genau diese Aussage fehlt einer Heuristik (Kapitel Metaheuristiken).")

    # --- 4. Wann es sich lohnt --------------------------------------------
    print("\n" + "-" * 84)
    print("4. Wann sich der Aufwand lohnt - und wann nicht\n")
    print(f"  {'Instanz':<36} {'Faustregel':>11} {'exakt':>8} {'Schranke':>10} "
          f"{'Ersparnis':>11}")
    print("  " + "-" * 80)
    for fall in (gross, schmal):
        ersparnis = (fall["first_fit"] - fall["ganzzahlig"]) / fall["first_fit"]
        print(f"  {fall['name']:<36} {fall['first_fit']:>11} "
              f"{fall['ganzzahlig']:>8} {fall['schranke']:>10.2f} "
              f"{ersparnis:>10.0%}")

    print(f"\n  Bei den breiten Zuschnitten spart die Spaltengenerierung "
          f"{gross['first_fit'] - gross['ganzzahlig']} von")
    print(f"  {gross['first_fit']} Rollen. Bei den schmalen spart sie NICHTS - dort ist die")
    print("  Faustregel bereits optimal.")
    print(f"\n  Der Unterschied haengt an einer einzigen Kennzahl:\n")
    print(f"  {'Instanz':<36} {'Stuecke je Rolle (etwa)':>24}")
    print("  " + "-" * 62)
    for fall in (gross, schmal):
        print(f"  {fall['name']:<36} {fall['stuecke_je_rolle']:>24.1f}")

    print("\n  Passen nur zwei bis drei Stuecke auf eine Rolle, entscheidet jede")
    print("  einzelne Zuordnung viel, und eine kurzsichtige Regel verschenkt")
    print("  ganze Rollen. Passen sechs oder mehr darauf, gleichen sich die")
    print("  Fehler aus - die Reste sind klein gegen die Rollenbreite, und die")
    print("  Faustregel trifft es fast immer.")

    print("\n" + "=" * 84)
    print("  WAS MAN DARAUS MITNIMMT")
    print("=" * 84)
    print("1. Der Perspektivwechsel ist die eigentliche Arbeit: nicht 'welches")
    print("   Stueck wohin', sondern 'wie oft welches Muster'. Damit verschwindet")
    print("   die Symmetrie, die das naheliegende Modell unloesbar macht.")
    print("2. Der Preis dafuer ist eine unaufschreibbare Zahl von Variablen. Die")
    print("   Spaltengenerierung zahlt ihn nicht, sondern erzeugt nur die wenigen")
    print("   Spalten, die die Dualwerte als lohnend ausweisen.")
    print("3. Das Teilproblem ist ein Rucksack - ein Verfahren, das die Leser")
    print("   dieses Buchs schon kennen. Neu ist allein die Schleife.")
    print("4. Und die unbequeme Erkenntnis: Ob sich das alles lohnt, entscheidet")
    print("   die Instanz, nicht die Methode. Bei sechs Stuecken je Rolle ist die")
    print("   Faustregel so gut wie das Optimum - und der Ertrag der")
    print("   Spaltengenerierung liegt dann allein in der SCHRANKE, die beweist,")
    print("   dass man aufhoeren kann zu suchen.")
    print("=" * 84)