#!/usr/bin/env python3 # Modellierungsschichten.py """ Kapitel Oekosystem: Pyomo und Linopy - zwei Modellierungsschichten fuer grosse Modelle. Geloest wird dasselbe Produktionsproblem wie im Vierfach-Vergleich: max 10*x1 + 15*x2 + 25*x3 u.d.N. x1 + x2 + 2*x3 <= 40 2*x1 + 3*x2 + x3 <= 50 x >= 0 Handrechnung: Z* = 530 bei x = (0, 12, 14). Der Vergleich zeigt die beiden Denkweisen: * Pyomo - algebraisch, indexbasiert, Industriestandard fuer Grossmodelle, trennt Modellstruktur sauber von den Daten (AbstractModel). * Linopy - beschriftete Arrays (xarray): eine Zeile Code erzeugt Tausende Nebenbedingungen auf einmal, ohne Python-Schleife. Beide bringen KEINEN eigenen Solver mit; hier rechnet in beiden Faellen HiGHS. WICHTIG: ortools wird in diesem Programm bewusst NICHT importiert - es vertraegt sich nicht mit der HiGHS-Kopie, die Pyomo und Linopy laden (siehe die Stolperfalle im Abschnitt 'Ein System - vier Programmieransaetze'). Benoetigt: pyomo, linopy, xarray, pandas, highspy, numpy """ from __future__ import annotations import time import numpy as np import pandas as pd import pyomo.environ as pyo import xarray as xr import linopy ERWARTET = 530.0 # Ergebnis der Handrechnung PRODUKTE = ["Standard", "Komfort", "Premium"] RESSOURCEN = ["Material", "Montage"] DECKUNGSBEITRAG = np.array([10.0, 15.0, 25.0]) VERBRAUCH = np.array([[1.0, 1.0, 2.0], # Material je Produkt [2.0, 3.0, 1.0]]) # Montage je Produkt VORRAT = np.array([40.0, 50.0]) # --- Pyomo: algebraisch und indexbasiert ------------------------------------ def loese_mit_pyomo() -> tuple[float, list[float], float]: """Pyomo denkt in Mengen und Indizes, wie ein Mathematiker es aufschreibt. `Constraint(RESSOURCEN, rule=...)` erzeugt fuer JEDES Element der Menge eine Nebenbedingung - das ist das 'fuer alle i' der Formelsprache, unmittelbar in Code uebersetzt. """ t0 = time.perf_counter() modell = pyo.ConcreteModel(name="Produktionsprogramm") modell.P = pyo.Set(initialize=PRODUKTE) modell.R = pyo.Set(initialize=RESSOURCEN) modell.db = pyo.Param(modell.P, initialize=dict(zip(PRODUKTE, DECKUNGSBEITRAG))) modell.a = pyo.Param(modell.R, modell.P, initialize={ (r, p): VERBRAUCH[i, j] for i, r in enumerate(RESSOURCEN) for j, p in enumerate(PRODUKTE)}) modell.vorrat = pyo.Param(modell.R, initialize=dict(zip(RESSOURCEN, VORRAT))) modell.x = pyo.Var(modell.P, domain=pyo.NonNegativeReals) modell.ziel = pyo.Objective( expr=sum(modell.db[p] * modell.x[p] for p in modell.P), sense=pyo.maximize) def kapazitaet(m, r): return sum(m.a[r, p] * m.x[p] for p in m.P) <= m.vorrat[r] modell.kapazitaet = pyo.Constraint(modell.R, rule=kapazitaet) ergebnis = pyo.SolverFactory("appsi_highs").solve(modell) dauer = time.perf_counter() - t0 status = ergebnis.solver.termination_condition if status != pyo.TerminationCondition.optimal: raise RuntimeError(f"Pyomo meldet Status: {status}") return (float(pyo.value(modell.ziel)), [float(pyo.value(modell.x[p])) for p in PRODUKTE], dauer) # --- Linopy: beschriftete Arrays -------------------------------------------- def loese_mit_linopy() -> tuple[float, list[float], float]: """Linopy denkt in beschrifteten Arrays (xarray). Der entscheidende Unterschied: `(verbrauch * x).sum("produkt") <= vorrat` ist EINE Zeile und erzeugt so viele Nebenbedingungen, wie die Dimension 'ressource' Eintraege hat. Bei 2 Ressourcen faellt das nicht auf, bei 200 000 schon - dort entstehen sie als Matrixoperation statt in einer Python-Schleife. """ t0 = time.perf_counter() # Benannte Indizes statt blosser Listen: Dadurch heissen die Achsen # 'produkt' und 'ressource', und xarray fuehrt sie beim Rechnen von allein # richtig zusammen. Ohne Namen vergibt linopy 'dim_0', und man muss # spaeter umbenennen - eine haeufige Stolperstelle. produkt = pd.Index(PRODUKTE, name="produkt") ressource = pd.Index(RESSOURCEN, name="ressource") modell = linopy.Model() modell.add_variables(lower=0, coords=[produkt], name="menge") x = modell.variables["menge"] db = xr.DataArray(DECKUNGSBEITRAG, coords=[produkt]) verbrauch = xr.DataArray(VERBRAUCH, coords=[ressource, produkt]) vorrat = xr.DataArray(VORRAT, coords=[ressource]) # EINE Zeile - sie erzeugt so viele Nebenbedingungen, wie die Achse # 'ressource' Eintraege hat. Genau das ist der Punkt. modell.add_constraints((verbrauch * x).sum("produkt") <= vorrat, name="kapazitaet") modell.add_objective((db * x).sum(), sense="max") modell.solve(solver_name="highs", output_flag=False) dauer = time.perf_counter() - t0 if modell.termination_condition != "optimal": raise RuntimeError(f"Linopy meldet Status: {modell.termination_condition}") loesung = modell.variables["menge"].solution.to_series() return (float(modell.objective.value), [float(loesung[p]) for p in PRODUKTE], dauer) if __name__ == "__main__": print("=" * 80) print(" MODELLIERUNGSSCHICHTEN FUER GROSSE MODELLE") print("=" * 80) print(f"{'Schicht':<14} {'Z*':>10} {'Standard':>10} {'Komfort':>10} " f"{'Premium':>10} {'Zeit':>10}") print("-" * 80) ergebnisse = [] for name, loeser in [("Pyomo", loese_mit_pyomo), ("Linopy", loese_mit_linopy)]: ziel, mengen, dauer = loeser() ergebnisse.append(ziel) print(f"{name:<14} {ziel:>10.2f} {mengen[0]:>10.2f} {mengen[1]:>10.2f} " f"{mengen[2]:>10.2f} {dauer:>8.2f} s") print("-" * 80) for ziel in ergebnisse: assert abs(ziel - ERWARTET) < 1e-6, \ f"Abweichung von der Handrechnung: {ziel} statt {ERWARTET}" print(f"Beide stimmen mit der Handrechnung ueberein (Z* = {ERWARTET:.0f}).") print() print("Wann welche Schicht?") print(" Pyomo -> wenn Modellstruktur und Daten getrennt bleiben sollen,") print(" wenn nichtlineare Terme oder MINLP dazukommen koennen,") print(" wenn spaeter ein kommerzieller Solver angebunden wird.") print(" Linopy -> wenn die Daten ohnehin als beschriftete Arrays vorliegen") print(" (Energiesystem-, Netz- und Zeitreihenmodelle) und das") print(" Modell zehntausende gleichartige Nebenbedingungen hat.") print("=" * 80)