#!/usr/bin/env python3 # or_kern.py """ Kapitel Praxisfallen: Der gemeinsame Unterbau fuer produktionsreife OR-Modelle. Alle Beispiele dieses Buchs loesen dieselben vier Aufgaben immer wieder: Daten einlesen und pruefen, ein Modell bauen, den Solverstatus auswerten, die Loesung gegen die Wirklichkeit kontrollieren. Dieses Modul zieht diese vier Aufgaben aus den Einzelprogrammen heraus. Rohdaten (Excel/CSV) -> Domaenenmodell (Pydantic, geprueft) -> Modellbauer (solverabhaengig) -> Loesung (DTO, solverunabhaengig) -> Abnahmepruefung Der Sinn dieser Trennung: Die Domaenenschicht weiss nichts von Solvern, die Loesungsschicht weiss nichts vom Modellaufbau. Wer den Solver wechselt, tauscht genau EINEN Baustein aus - der Rest bleibt. WICHTIG ZUM IMPORT: Dieses Modul importiert von sich aus KEINE Solver- bibliothek. Der Grund steht im Kapitel Oekosystem: ortools und highspy vertragen sich nicht im selben Prozess. Die Statusuebersetzer laden ihre Bibliothek erst, wenn sie aufgerufen werden - so bleibt or_kern.py mit jedem Solver benutzbar. Benoetigt: pydantic (v2), numpy; optional pandas und openpyxl fuer Excel. """ from __future__ import annotations from enum import Enum from typing import Annotated, Any, Sequence import numpy as np from pydantic import BaseModel, Field, model_validator # Wiederverwendbare Feldtypen. Sie tragen die Pruefung im Typ, nicht im Code - # damit gilt sie ueberall, wo der Typ verwendet wird. PositiveZahl = Annotated[float, Field(gt=0)] NichtNegativ = Annotated[float, Field(ge=0)] # --- 1. Solverstatus: eine Sprache fuer fuenf Bibliotheken ------------------- class SolverStatus(str, Enum): """Was ein Solverlauf ergeben hat - unabhaengig davon, wer gerechnet hat. Die fuenf im Buch verwendeten Bibliotheken benennen dasselbe Ergebnis unterschiedlich: 'Optimal', 'OPTIMAL', 2, 'optimal'. Wer darauf direkt prueft, bindet seinen Auswertungscode an eine Bibliothek. """ OPTIMAL = "optimal" # beweisbar bestmoeglich ZULAESSIG = "zulaessig" # brauchbare Loesung, Beweis fehlt UNZULAESSIG = "unzulaessig" # es gibt keine Loesung (Aussage ueber das Modell) UNBESCHRAENKT = "unbeschraenkt" # Ziel waechst ins Unendliche ZEITLIMIT = "zeitlimit" # abgebrochen, nichts Brauchbares gefunden FEHLERHAFT = "fehlerhaft" # das Modell ist kein gueltiges Modell UNBEKANNT = "unbekannt" # alles Uebrige @property def brauchbar(self) -> bool: """Habe ich etwas in der Hand, das ich ausfuehren kann?""" return self in (SolverStatus.OPTIMAL, SolverStatus.ZULAESSIG) @property def modellfehler(self) -> bool: """Liegt die Ursache im Modell (und nicht in der Rechenzeit)?""" return self in (SolverStatus.UNZULAESSIG, SolverStatus.UNBESCHRAENKT, SolverStatus.FEHLERHAFT) def status_von_pywraplp(rohstatus: int) -> SolverStatus: """OR-Tools linear_solver (GLOP, SCIP, CBC).""" from ortools.linear_solver import pywraplp zuordnung = { pywraplp.Solver.OPTIMAL: SolverStatus.OPTIMAL, pywraplp.Solver.FEASIBLE: SolverStatus.ZULAESSIG, pywraplp.Solver.INFEASIBLE: SolverStatus.UNZULAESSIG, pywraplp.Solver.UNBOUNDED: SolverStatus.UNBESCHRAENKT, pywraplp.Solver.ABNORMAL: SolverStatus.FEHLERHAFT, pywraplp.Solver.NOT_SOLVED: SolverStatus.ZEITLIMIT, } return zuordnung.get(rohstatus, SolverStatus.UNBEKANNT) def status_von_cpsat(rohstatus: int) -> SolverStatus: """OR-Tools CP-SAT (siehe Kapitel CP-SAT, alle fuenf Faelle).""" from ortools.sat.python import cp_model zuordnung = { cp_model.OPTIMAL: SolverStatus.OPTIMAL, cp_model.FEASIBLE: SolverStatus.ZULAESSIG, cp_model.INFEASIBLE: SolverStatus.UNZULAESSIG, cp_model.MODEL_INVALID: SolverStatus.FEHLERHAFT, cp_model.UNKNOWN: SolverStatus.ZEITLIMIT, } return zuordnung.get(rohstatus, SolverStatus.UNBEKANNT) def status_von_highs(statustext: str) -> SolverStatus: """HiGHS ueber highspy - hier kommt der Status als Text.""" zuordnung = { "Optimal": SolverStatus.OPTIMAL, "Time limit reached": SolverStatus.ZULAESSIG, # meist mit Loesung "Solution limit reached": SolverStatus.ZULAESSIG, "Infeasible": SolverStatus.UNZULAESSIG, "Unbounded": SolverStatus.UNBESCHRAENKT, "Primal infeasible or unbounded": SolverStatus.UNZULAESSIG, "Model error": SolverStatus.FEHLERHAFT, } return zuordnung.get(statustext, SolverStatus.UNBEKANNT) def status_von_scipy(ergebnis: Any) -> SolverStatus: """scipy.optimize.linprog und milp liefern ein Ergebnisobjekt mit .status.""" zuordnung = { 0: SolverStatus.OPTIMAL, 1: SolverStatus.ZEITLIMIT, # Iterations-/Zeitgrenze 2: SolverStatus.UNZULAESSIG, 3: SolverStatus.UNBESCHRAENKT, 4: SolverStatus.FEHLERHAFT, # numerische Schwierigkeiten } return zuordnung.get(int(ergebnis.status), SolverStatus.UNBEKANNT) def status_von_cvxpy(statustext: str) -> SolverStatus: """CVXPY - Statuszeichenketten wie 'optimal' oder 'infeasible'.""" zuordnung = { "optimal": SolverStatus.OPTIMAL, "optimal_inaccurate": SolverStatus.ZULAESSIG, "infeasible": SolverStatus.UNZULAESSIG, "infeasible_inaccurate": SolverStatus.UNZULAESSIG, "unbounded": SolverStatus.UNBESCHRAENKT, "unbounded_inaccurate": SolverStatus.UNBESCHRAENKT, } return zuordnung.get(statustext, SolverStatus.UNBEKANNT) # --- 2. Die Loesung als solverunabhaengiges Datenobjekt ---------------------- class Loesung(BaseModel): """Das Ergebnis eines Solverlaufs, so wie es weiterverarbeitet wird. Bewusst OHNE Referenz auf Solver, Modell oder Variablen: Ein Bericht, ein Test oder eine Weiterverarbeitung soll nicht wissen muessen, womit gerechnet wurde. """ status: SolverStatus werte: dict[str, float] = Field(default_factory=dict) zielwert: float | None = None schranke: float | None = None laufzeit: float | None = None # Schattenpreise je Nebenbedingung, sofern der Solver sie liefert # (nur bei kontinuierlichen Problemen, siehe Kapitel LP). Auch das ist # ein solverunabhaengiger Begriff und gehoert deshalb hierher. schattenpreise: dict[str, float] = Field(default_factory=dict) @property def gap(self) -> float | None: """Relativer Abstand zwischen gefundener Loesung und bewiesener Schranke. Das ist die Zahl fuer den Bericht: keine Schaetzung, sondern eine Zusage - schlechter als das kann die Loesung nicht sein. """ if self.zielwert is None or self.schranke is None: return None nenner = max(abs(self.zielwert), 1e-9) return abs(self.zielwert - self.schranke) / nenner def als_bericht(self) -> str: """Eine Zeile fuers Betriebsprotokoll - siehe Betriebsueberwachung.py.""" teile = [f"Status: {self.status.value}"] if self.zielwert is not None: teile.append(f"Zielwert: {self.zielwert:,.2f}") if self.gap is not None: teile.append(f"Gap: {self.gap:.2%}") if self.laufzeit is not None: teile.append(f"Zeit: {self.laufzeit:.2f}s") return " | ".join(teile) # --- 3. Domaenenmodell: die Fakten, bevor ein Solver sie sieht --------------- class Produkt(BaseModel): """Ein Produkt mit Deckungsbeitrag und Ressourcenverbrauch.""" name: str = Field(min_length=1) deckungsbeitrag: float verbrauch: dict[str, NichtNegativ] model_config = {"frozen": True} # Stammdaten aendern sich nicht im Lauf class Produktionsproblem(BaseModel): """Produktionsprogrammplanung - das Modell aus Kapitel Einfuehrung. Die Pruefungen stehen hier und nicht im Solvercode. Damit gelten sie unabhaengig davon, mit welcher Bibliothek spaeter gerechnet wird - und sie schlagen beim EINLESEN zu, wo der Fehler noch zuzuordnen ist. """ produkte: list[Produkt] = Field(min_length=1) kapazitaeten: dict[str, PositiveZahl] @model_validator(mode="after") def pruefe_ressourcen(self) -> "Produktionsproblem": namen = [p.name for p in self.produkte] if len(set(namen)) != len(namen): doppelt = sorted({n for n in namen if namen.count(n) > 1}) raise ValueError(f"Produktnamen kommen mehrfach vor: {doppelt}") for produkt in self.produkte: unbekannt = set(produkt.verbrauch) - set(self.kapazitaeten) if unbekannt: raise ValueError( f"Produkt '{produkt.name}' verbraucht Ressourcen ohne " f"Kapazitaetsangabe: {sorted(unbekannt)}") return self @property def ressourcen(self) -> list[str]: """Feste Reihenfolge - siehe die Spaltenfalle in Kapitel Finanzdaten.""" return sorted(self.kapazitaeten) def verbrauchsmatrix(self) -> np.ndarray: """Zeilen = Ressourcen, Spalten = Produkte (Matrixform, Kapitel Fundament).""" return np.array([[p.verbrauch.get(r, 0.0) for p in self.produkte] for r in self.ressourcen]) def kapazitaetsvektor(self) -> np.ndarray: return np.array([self.kapazitaeten[r] for r in self.ressourcen]) def deckungsbeitragsvektor(self) -> np.ndarray: return np.array([p.deckungsbeitrag for p in self.produkte]) # --- 4. Abnahmepruefung: Loesung gegen Anforderung, ohne Solver -------------- def pruefe_loesung(problem: Produktionsproblem, loesung: Loesung, ganzzahlig: Sequence[str] = (), toleranz: float = 1e-6) -> list[str]: """Prueft eine Loesung gegen die Anforderungen - OHNE den Solver zu fragen. Diese Funktion darf ausdruecklich KEIN Solverobjekt und keine Modellvariable benutzen. Eine Pruefung aus denselben Bausteinen wie das Modell prueft das Modell gegen sich selbst und findet nichts (siehe die Denkfehler in den Kapiteln Graphen und MILP). Liefert eine Liste von Beanstandungen; leer heisst bestanden. """ beanstandungen: list[str] = [] if not loesung.status.brauchbar: return [f"Kein verwertbares Ergebnis (Status: {loesung.status.value})"] fehlend = [p.name for p in problem.produkte if p.name not in loesung.werte] if fehlend: return [f"Loesung enthaelt keine Werte fuer: {fehlend}"] mengen = np.array([loesung.werte[p.name] for p in problem.produkte]) if (mengen < -toleranz).any(): negativ = [p.name for p, m in zip(problem.produkte, mengen) if m < -toleranz] beanstandungen.append(f"negative Mengen bei: {negativ}") verbrauch = problem.verbrauchsmatrix() @ mengen for ressource, ist, grenze in zip(problem.ressourcen, verbrauch, problem.kapazitaetsvektor()): if ist > grenze + toleranz: beanstandungen.append( f"{ressource}: Verbrauch {ist:,.3f} ueber Kapazitaet {grenze:,.3f}") for name in ganzzahlig: wert = loesung.werte[name] if abs(wert - round(wert)) > toleranz: beanstandungen.append(f"'{name}' = {wert!r} ist nicht ganzzahlig") if loesung.zielwert is not None: nachgerechnet = float(problem.deckungsbeitragsvektor() @ mengen) if abs(nachgerechnet - loesung.zielwert) > toleranz * max(1.0, abs(nachgerechnet)): beanstandungen.append( f"Zielwert {loesung.zielwert:,.4f} passt nicht zu den Mengen " f"(nachgerechnet {nachgerechnet:,.4f})") return beanstandungen # --- 5. Ein Modellbauer je Solver - austauschbar ---------------------------- def loese_mit_glop(problem: Produktionsproblem) -> Loesung: """Modellbauer fuer OR-Tools GLOP (kontinuierlich).""" import time from ortools.linear_solver import pywraplp solver = pywraplp.Solver.CreateSolver("GLOP") menge = {p.name: solver.NumVar(0, solver.infinity(), p.name) for p in problem.produkte} # Die Nebenbedingung wird gemerkt - ohne die Referenz gibt es spaeter # keinen Schattenpreis (siehe Kapitel LP). bedingung = {} for ressource in problem.ressourcen: bedingung[ressource] = solver.Add( sum(menge[p.name] * p.verbrauch.get(ressource, 0.0) for p in problem.produkte) <= problem.kapazitaeten[ressource], name=ressource) solver.Maximize(sum(menge[p.name] * p.deckungsbeitrag for p in problem.produkte)) t0 = time.perf_counter() rohstatus = solver.Solve() laufzeit = time.perf_counter() - t0 status = status_von_pywraplp(rohstatus) if not status.brauchbar: return Loesung(status=status, laufzeit=laufzeit) return Loesung(status=status, werte={name: v.solution_value() for name, v in menge.items()}, zielwert=solver.Objective().Value(), schranke=solver.Objective().Value(), schattenpreise={r: b.dual_value() for r, b in bedingung.items()}, laufzeit=laufzeit) def loese_mit_scipy(problem: Produktionsproblem) -> Loesung: """Derselbe Fall mit scipy.optimize.linprog - anderer Solver, gleiche Loesung. Beachten Sie, wie wenig hier steht: Die Daten kommen fertig geprueft aus dem Domaenenmodell, die Auswertung geht ans gemeinsame Loesung-Objekt. Genau das ist der Ertrag der Trennung. """ import time from scipy.optimize import linprog t0 = time.perf_counter() ergebnis = linprog(c=-problem.deckungsbeitragsvektor(), # linprog minimiert A_ub=problem.verbrauchsmatrix(), b_ub=problem.kapazitaetsvektor(), bounds=(0, None), method="highs") laufzeit = time.perf_counter() - t0 status = status_von_scipy(ergebnis) if not status.brauchbar: return Loesung(status=status, laufzeit=laufzeit) # linprog rechnet minimierend mit negierten Kosten - die Dualwerte # muessen entsprechend zurueckgedreht werden. return Loesung(status=status, werte={p.name: float(x) for p, x in zip(problem.produkte, ergebnis.x)}, zielwert=float(-ergebnis.fun), schranke=float(-ergebnis.fun), schattenpreise={r: float(-m) for r, m in zip(problem.ressourcen, ergebnis.ineqlin.marginals)}, laufzeit=laufzeit) # --- 6. Ein- und Ausgabe: Tabellen sind die Wirklichkeit --------------------- # # In den meisten Betrieben liegen die Zahlen in einer Tabellenkalkulation # (Kapitel Einfuehrung). Diese beiden Funktionen sind die einzige Stelle im # Modul, die das weiss - alles andere arbeitet mit dem Domaenenmodell. # # pandas und openpyxl werden bewusst erst hier importiert: Wer or_kern nur # fuer Statusauswertung und Pruefung benutzt, soll sie nicht installieren # muessen. STAMMSPALTEN = ("Produkt", "Deckungsbeitrag") def lade_produktionsproblem(pfad: str, blatt_produkte: str = "Produkte", blatt_kapazitaeten: str = "Kapazitaeten" ) -> Produktionsproblem: """Liest eine Excel-Mappe und liefert ein geprueftes Domaenenmodell. Alle Spalten ausser 'Produkt' und 'Deckungsbeitrag' gelten als Verbrauchsspalten - eine neue Ressource ist damit eine neue SPALTE in der Mappe und erfordert keine Codeaenderung. Die eigentliche Pruefung passiert nicht hier, sondern im Konstruktor von Produktionsproblem. Diese Funktion uebersetzt nur Tabelle in Objekt; was ein gueltiges Problem ist, steht an genau einer Stelle. """ import pandas as pd produkte_tabelle = pd.read_excel(pfad, sheet_name=blatt_produkte) kapazitaeten_tabelle = pd.read_excel(pfad, sheet_name=blatt_kapazitaeten) fehlende = [s for s in STAMMSPALTEN if s not in produkte_tabelle.columns] if fehlende: raise ValueError(f"Blatt '{blatt_produkte}': Spalten fehlen: {fehlende}") if produkte_tabelle.isna().any().any(): zeilen = produkte_tabelle[produkte_tabelle.isna().any(axis=1)].index.tolist() raise ValueError(f"Blatt '{blatt_produkte}': leere Zellen in Zeile(n) {zeilen}") ressourcen = [s for s in produkte_tabelle.columns if s not in STAMMSPALTEN] produkte = [ Produkt(name=str(zeile.Produkt), deckungsbeitrag=float(zeile.Deckungsbeitrag), verbrauch={r: float(getattr(zeile, r)) for r in ressourcen}) for zeile in produkte_tabelle.itertuples() ] kapazitaeten = {str(z.Ressource): float(z.Verfuegbar) for z in kapazitaeten_tabelle.itertuples()} return Produktionsproblem(produkte=produkte, kapazitaeten=kapazitaeten) def schreibe_ergebnis(pfad: str, problem: Produktionsproblem, loesung: Loesung) -> None: """Schreibt Plan und Kennzahlen als zwei Blaetter zurueck nach Excel. Die Fachabteilung bekommt das Ergebnis im Format, das sie ohnehin benutzt. Akzeptanz entscheidet ueber Projekterfolg (Kapitel Praxisfallen). """ import pandas as pd if not loesung.status.brauchbar: raise ValueError(f"Nichts zu schreiben - Status: {loesung.status.value}") mengen = np.array([loesung.werte[p.name] for p in problem.produkte]) plan = pd.DataFrame({ "Produkt": [p.name for p in problem.produkte], "Menge": mengen, "Deckungsbeitrag": mengen * problem.deckungsbeitragsvektor(), }) verbrauch = problem.verbrauchsmatrix() @ mengen kennzahlen: dict[str, float | str] = { "Status": loesung.status.value, "Gesamtdeckungsbeitrag": loesung.zielwert or 0.0, } for ressource, ist, grenze in zip(problem.ressourcen, verbrauch, problem.kapazitaetsvektor()): kennzahlen[f"{ressource}: Verbrauch"] = float(ist) kennzahlen[f"{ressource}: Auslastung %"] = float(100.0 * ist / grenze) if ressource in loesung.schattenpreise: kennzahlen[f"{ressource}: Schattenpreis"] = \ loesung.schattenpreise[ressource] kennzahlen_tabelle = pd.DataFrame({ "Kennzahl": list(kennzahlen), "Wert": [round(w, 2) if isinstance(w, float) else w for w in kennzahlen.values()], }) with pd.ExcelWriter(pfad, engine="openpyxl") as mappe: plan.round(2).to_excel(mappe, sheet_name="Plan", index=False) kennzahlen_tabelle.to_excel(mappe, sheet_name="Kennzahlen", index=False) if __name__ == "__main__": # Die Schreinerei aus Kapitel Einfuehrung - jetzt als Domaenenmodell. schreinerei = Produktionsproblem( produkte=[ Produkt(name="Tisch", deckungsbeitrag=240.0, verbrauch={"Montagestunden": 3.0, "Plattenmaterial": 6.0}), Produkt(name="Stuhl", deckungsbeitrag=60.0, verbrauch={"Montagestunden": 1.0, "Plattenmaterial": 1.0}), ], kapazitaeten={"Montagestunden": 150.0, "Plattenmaterial": 240.0}, ) print("=" * 78) print(" DERSELBE FALL, ZWEI SOLVER, EIN AUSWERTUNGSCODE") print("=" * 78) for name, bauer in [("OR-Tools GLOP", loese_mit_glop), ("scipy / HiGHS", loese_mit_scipy)]: loesung = bauer(schreinerei) beanstandungen = pruefe_loesung(schreinerei, loesung) print(f"\n{name}") print(f" {loesung.als_bericht()}") for produkt, menge in loesung.werte.items(): print(f" {produkt:<10} {menge:8.2f}") print(f" Abnahmepruefung: " f"{'bestanden' if not beanstandungen else beanstandungen}") print("\n" + "=" * 78) print(" WAS DIE PRUEFUNGEN ABFANGEN") print("=" * 78) faelle = [ ("Kapazitaet 0", dict( produkte=schreinerei.produkte, kapazitaeten={"Montagestunden": 0.0, "Plattenmaterial": 240.0})), ("Ressource ohne Kapazitaet", dict( produkte=[Produkt(name="Regal", deckungsbeitrag=130.0, verbrauch={"Lackieren": 2.0})], kapazitaeten={"Montagestunden": 150.0})), ("doppelter Produktname", dict( produkte=[schreinerei.produkte[0], schreinerei.produkte[0]], kapazitaeten=schreinerei.kapazitaeten)), ] for beschreibung, daten in faelle: try: Produktionsproblem(**daten) print(f" {beschreibung:<28} NICHT erkannt (!)") except Exception as fehler: meldung = str(fehler).splitlines() kern = next((z.strip() for z in meldung if "Value error" in z or "greater than" in z), meldung[-1]) print(f" {beschreibung:<28} abgefangen: {kern[:44]}") print("\nAlle drei scheitern beim EINLESEN - nicht erst beim Loesen und") print("schon gar nicht erst im Bericht. Das ist der ganze Zweck der") print("Domaenenschicht.") print("=" * 78)