#!/usr/bin/env python3

# Lokale_Optima_Multistart.py
"""
Kapitel QP/NLP: Was passiert, wenn die Konvexitaet fehlt.

CVXPY verweigert nicht-konvexe Probleme - das ist sein Schutzmechanismus.
scipy.optimize.minimize verweigert nichts. Es rechnet, meldet 'success: True'
und liefert ein Ergebnis. Nur ist das Ergebnis dann kein Optimum, sondern
irgendein lokales Minimum, das vom Startpunkt abhaengt.

Beispiel aus dem Einkauf: 700 Tonnen Rohstoff werden auf vier Lieferanten
verteilt. Jeder gewaehrt einen MENGENRABATT - der Stueckpreis faellt, je mehr
man bei ihm bestellt:

    preis_i(q) = basis_i * (1 - rabatt_i * (1 - exp(-q / skala_i)))

Genau das macht die Zielfunktion nicht-konvex: Grosse Bestellungen lohnen
sich ueberproportional, es gibt also mehrere sinnvolle "Cluster"-Loesungen -
und dazwischen schlechtere Taeler.

Das Programm zeigt drei Dinge:
  1. Ein einzelner Lauf liefert ein plausibles Ergebnis - ohne jede Warnung.
  2. 200 Startpunkte foerdern mehrere verschiedene lokale Optima zutage.
  3. Der Unterschied zwischen bestem und schlechtestem betraegt hier 8 %.

Benoetigt: numpy, scipy
"""

from __future__ import annotations

import numpy as np
from scipy.optimize import minimize

# Vier Lieferanten
NAMEN = ["Nord AG", "Ost GmbH", "Sued KG", "West SE"]
BASISPREIS = np.array([50.0, 47.0, 53.0, 45.0])     # EUR je Tonne ohne Rabatt
MAX_RABATT = np.array([0.30, 0.18, 0.35, 0.12])     # hoechstmoeglicher Rabatt
RABATT_SKALA = np.array([120.0, 260.0, 90.0, 300.0])  # wie schnell er greift
KAPAZITAET = np.array([400.0, 400.0, 400.0, 400.0])
BEDARF = 700.0

RNG = np.random.default_rng(0)


def gesamtkosten(menge: np.ndarray) -> float:
    """Einkaufskosten bei mengenabhaengigem Stueckpreis.

    Der Rabatt waechst mit der Bestellmenge und laeuft gegen MAX_RABATT.
    Dadurch ist der Stueckpreis fallend - und die Gesamtkostenfunktion
    nicht mehr konvex.
    """
    stueckpreis = BASISPREIS * (1 - MAX_RABATT * (1 - np.exp(-menge / RABATT_SKALA)))
    return float(stueckpreis @ menge)


NEBENBEDINGUNGEN = [{"type": "eq", "fun": lambda q: q.sum() - BEDARF}]
GRENZEN = [(0.0, k) for k in KAPAZITAET]


def optimiere_von(startpunkt: np.ndarray):
    """Ein Lauf von einem gegebenen Startpunkt aus."""
    return minimize(gesamtkosten, startpunkt, method="SLSQP",
                    bounds=GRENZEN, constraints=NEBENBEDINGUNGEN)


def zufaelliger_start() -> np.ndarray:
    """Zufaellige Aufteilung, die den Bedarf bereits erfuellt."""
    anteil = RNG.random(len(NAMEN))
    return anteil / anteil.sum() * BEDARF


def zeige_plan(titel: str, menge: np.ndarray, kosten: float) -> None:
    print(f"\n{titel}")
    for name, m in zip(NAMEN, menge):
        anteil = m / BEDARF * 100
        print(f"    {name:<10} {m:7.1f} t  ({anteil:4.1f} %)")
    print(f"    {'Gesamtkosten':<10} {kosten:9,.2f} EUR")


if __name__ == "__main__":
    print("=" * 78)
    print("  NICHT-KONVEX: DERSELBE CODE, VERSCHIEDENE ERGEBNISSE")
    print("=" * 78)
    print(f"{BEDARF:.0f} t Rohstoff auf {len(NAMEN)} Lieferanten mit Mengenrabatt.")

    # --- 1. Ein einziger Lauf, so wie man es zuerst schreibt --------------
    erster = optimiere_von(np.full(len(NAMEN), BEDARF / len(NAMEN)))
    print(f"\n[1] EIN Lauf, Startpunkt 'gleichmaessig verteilt'")
    print(f"    scipy meldet: success={erster.success}, "
          f"'{erster.message}'")
    zeige_plan("    Ergebnis:", erster.x, erster.fun)
    print("\n    Nichts an dieser Ausgabe deutet darauf hin, dass etwas fehlt.")

    # --- 1b. Der kaufmaennisch naheliegende Startpunkt --------------------
    # "Kaufe bei den beiden Lieferanten mit dem guenstigsten Basispreis" -
    # West SE (45) und Ost GmbH (47), jeweils bis zur Kapazitaetsgrenze.
    guenstigste = np.argsort(BASISPREIS)[:2]
    kaufmaennisch = np.zeros(len(NAMEN))
    rest = BEDARF
    for i in guenstigste:
        kaufmaennisch[i] = min(KAPAZITAET[i], rest)
        rest -= kaufmaennisch[i]
    zweiter = optimiere_von(kaufmaennisch)
    print(f"\n[1b] EIN Lauf, Startpunkt 'die zwei mit dem guenstigsten Basispreis'")
    print(f"    scipy meldet: success={zweiter.success}")
    zeige_plan("    Ergebnis:", zweiter.x, zweiter.fun)
    print(f"\n    Dasselbe Programm, derselbe Aufruf, ein anderer Startpunkt -")
    print(f"    und {zweiter.fun - erster.fun:,.2f} EUR Unterschied "
          f"({(zweiter.fun / erster.fun - 1) * 100:.1f} %).")

    # --- 2. Multistart: dasselbe Problem, viele Startpunkte ---------------
    laeufe = []
    for _ in range(200):
        ergebnis = optimiere_von(zufaelliger_start())
        if ergebnis.success:
            laeufe.append((float(ergebnis.fun), ergebnis.x))

    # Ergebnisse, die sich um weniger als 1 Cent unterscheiden, sind dasselbe
    # lokale Optimum - zusammenfassen, sonst zaehlt man Rundungsrauschen.
    optima: list[tuple[float, np.ndarray]] = []
    for wert, plan in sorted(laeufe, key=lambda t: t[0]):
        if not optima or abs(wert - optima[-1][0]) > 0.01:
            optima.append((wert, plan))

    print("\n" + "=" * 78)
    print(f"[2] 200 zufaellige Startpunkte -> {len(laeufe)} erfolgreiche Laeufe")
    print(f"    darunter {len(optima)} VERSCHIEDENE lokale Optima:")
    print()
    print(f"    {'Rang':>5} {'Kosten':>13} {'Abstand zum besten':>20}   Aufteilung (t)")
    print("    " + "-" * 70)
    bester = optima[0][0]
    for rang, (wert, plan) in enumerate(optima, start=1):
        abstand = (wert / bester - 1) * 100
        aufteilung = " ".join(f"{m:5.0f}" for m in plan)
        print(f"    {rang:>5} {wert:>13,.2f} {abstand:>19.2f} %   {aufteilung}")

    # --- 3. Was das kostet ------------------------------------------------
    schlechtester = optima[-1]
    print("\n" + "=" * 78)
    print("  WAS AUF DEM SPIEL STEHT")
    print("=" * 78)
    zeige_plan("Bester gefundener Plan:", optima[0][1], optima[0][0])
    zeige_plan("Schlechtestes lokales Optimum:", schlechtester[1], schlechtester[0])
    unterschied = schlechtester[0] - optima[0][0]
    print(f"\n  Unterschied: {unterschied:,.2f} EUR "
          f"({unterschied / optima[0][0] * 100:.1f} %)")
    print(f"\n  Lauf [1]  (gleichmaessiger Start):     {erster.fun:>10,.2f} EUR")
    print(f"  Lauf [1b] (kaufmaennischer Start):    {zweiter.fun:>10,.2f} EUR"
          f"   <- {(zweiter.fun / optima[0][0] - 1) * 100:.1f} % ueber dem besten")
    print()
    print("  Bemerkenswert: Der kaufmaennisch NAHELIEGENDE Startpunkt fuehrt in")
    print("  das schlechteste Ergebnis von allen - schlechter als jedes der 200")
    print("  zufaellig gefundenen lokalen Optima. Wer beim guenstigsten")
    print("  Basispreis anfaengt, uebersieht, dass hier der Mengenrabatt")
    print("  entscheidet und nicht der Listenpreis.")

    print()
    print("Drei Konsequenzen fuer die Praxis:")
    print("  1. 'success: True' heisst bei nicht-konvexen Problemen NICHT 'optimal'.")
    print("     Es heisst nur: 'Ich bin an einer Stelle angekommen, an der es in")
    print("     keine Richtung mehr bergab geht.'")
    print("  2. Ein einzelner Lauf ist wertlos. Nehmen Sie viele Startpunkte und")
    print("     berichten Sie die STREUUNG mit - sie ist Ihre einzige Auskunft")
    print("     darueber, wie zerklueftet die Landschaft ist.")
    print("  3. Auch Multistart liefert KEINE Garantie. Dass hier nichts unter")
    print(f"     {bester:,.2f} EUR gefunden wurde, beweist nicht, dass es nichts gibt.")
    print("=" * 78)