Anhang D: Spickzettel der Solver

Wofür dieser Anhang gedacht ist: Sie wissen, was Sie modellieren wollen, und suchen nur noch, wie die gewählte Bibliothek es schreibt. Jede Seite hat denselben Aufbau — Modell, Variablen, Nebenbedingungen, Lösen, alle Statusfälle, Lösung auslesen, Stolpersteine. So lassen sich die Seiten nebeneinanderlegen.

Alle Schnipsel lösen dasselbe Problem — das Produktionsprogramm aus Abschnitt 3.1 mit dem bekannten Optimum 530 und den Schattenpreisen 12 und 1:

\max\; 10x_1 + 15x_2 + 25x_3 \quad\text{u.d.N.}\quad x_1 + x_2 + 2x_3 \le 40,\;\; 2x_1 + 3x_2 + x_3 \le 50,\;\; x \ge 0

Damit ist jeder Schnipsel selbstprüfend: Kommt bei Ihnen etwas anderes als 530 heraus, liegt es an der Übertragung, nicht am Modell. (Die Ausnahme ist CP-SAT — ein rein stetiges LP ist dort das falsche Werkzeug; die Seite zeigt stattdessen die CP-SAT-eigenen Bausteine.)

Was dieser Anhang nicht ist. Kein Vergleich und keine Empfehlung. Welche Bibliothek für welche Aufgabe taugt, steht in Abschnitt 3.6; denselben Fall in vier Bibliotheken nebeneinander zeigt Ein_System_Vier_Ansaetze.py, die beiden Modellierungssprachen Modellierungsschichten.py. Hier geht es allein ums Nachschlagen.

D1 — SciPy: `linprog` und `milp`

Wofür. Die Einstiegsschicht: keine zusätzliche Installation, HiGHS als Unterbau, ideal für lineare und gemischt-ganzzahlige Probleme in Matrixform. Wofür nicht: alles, was sich nicht als Matrix schreiben lässt, und jede nichtlineare Zielfunktion.

Lineares Programm

import numpy as np
from scipy.optimize import linprog

# linprog MINIMIERT immer -> zum Maximieren die Zielfunktion negieren
ergebnis = linprog(
    c=[-10.0, -15.0, -25.0],                       # Zielkoeffizienten
    A_ub=[[1, 1, 2], [2, 3, 1]], b_ub=[40, 50],    # A_ub @ x <= b_ub
    A_eq=None, b_eq=None,                          # Gleichungen, falls vorhanden
    bounds=[(0, None)] * 3,                        # je Variable (unten, oben)
    method="highs")

if ergebnis.status == 0:
    print(f"optimal: {-ergebnis.fun:.2f}")         # Vorzeichen zuruecknehmen!
    print(f"x = {np.round(ergebnis.x, 4)}")
    print(f"Schattenpreise: {-ergebnis.ineqlin.marginals}")
elif ergebnis.status == 2:
    print("INFEASIBLE - kein zulaessiger Punkt")
elif ergebnis.status == 3:
    print("UNBOUNDED - Zielfunktion unbeschraenkt")
elif ergebnis.status == 1:
    print("Iterations- oder Zeitlimit erreicht")
else:
    print(f"numerisches Problem (status {ergebnis.status}): {ergebnis.message}")

Ganzzahlig: `milp`

import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds

# Ganzzahlig: dasselbe Problem, aber x muss ganzzahlig sein
ergebnis = milp(
    c=[-10.0, -15.0, -25.0],                       # auch milp MINIMIERT
    constraints=LinearConstraint([[1, 1, 2], [2, 3, 1]], -np.inf, [40, 50]),
    integrality=[1, 1, 1],                         # 0 = stetig, 1 = ganzzahlig
    bounds=Bounds(0, np.inf))

if ergebnis.status == 0:
    print(f"optimal: {-ergebnis.fun:.2f}  x = {np.round(ergebnis.x).astype(int)}")
    print(f"MIP-Gap: {ergebnis.mip_gap:.4f}")
elif ergebnis.status == 1:
    print("Zeitlimit - beste gefundene Loesung nutzen, Gap pruefen")
elif ergebnis.status == 2:
    print("INFEASIBLE")
elif ergebnis.status == 3:
    print("UNBOUNDED")
else:
    print(f"kein Ergebnis: {ergebnis.message}")

Die drei häufigsten Stolpersteine

linprog minimiert immer. Zum Maximieren c negieren — und beim Ausgeben des Zielwerts das Vorzeichen wieder zurücknehmen. Dieselbe Negation dreht auch die Schattenpreise (Abschnitt 5.7).
status == 0 prüfen, nicht res.success allein. success ist bei status == 1 (Limit erreicht) False, obwohl eine brauchbare Lösung vorliegen kann.
bounds gilt je Variable. bounds=(0, None) setzt alle Variablen gleich; bounds=[(0, None), (0, 10), ...] einzeln. Wer die Liste vergisst, bekommt stillschweigend überall dieselbe Schranke.

D2 — HiGHS über `highspy`

Wofür. Derselbe Solver wie unter SciPy, aber direkt gesteuert: Optionen, Warm-Starts, inkrementelles Ändern eines bestehenden Modells. Wofür nicht: schnelles Hinschreiben — die CSR-Matrixübergabe ist fehleranfällig.

import numpy as np
import highspy

h = highspy.Highs()
h.setOptionValue("output_flag", False)             # Solverprotokoll abschalten
h.setOptionValue("time_limit", 60.0)
h.setOptionValue("mip_rel_gap", 0.01)              # 1 % Gap genuegt

# Variablen: Anzahl, Untergrenzen, Obergrenzen
h.addVars(3, np.zeros(3), np.full(3, highspy.kHighsInf))
h.changeObjectiveSense(highspy.ObjSense.kMaximize)
for spalte, wert in enumerate([10.0, 15.0, 25.0]):
    h.changeColCost(spalte, wert)

# Zeilen im CSR-Format: starts[i] = Beginn von Zeile i in indices/values
h.addRows(2, np.full(2, -highspy.kHighsInf), np.array([40.0, 50.0]), 6,
          np.array([0, 3], dtype=np.int32),                  # starts
          np.array([0, 1, 2, 0, 1, 2], dtype=np.int32),      # Spaltenindizes
          np.array([1.0, 1.0, 2.0, 2.0, 3.0, 1.0]))          # Koeffizienten
# Ganzzahligkeit: h.changeColsIntegrality(...) mit highspy.HighsVarType.kInteger

h.run()
status = h.getModelStatus()
if status == highspy.HighsModelStatus.kOptimal:
    print(f"optimal: {h.getInfo().objective_function_value:.2f}")
    print(f"x = {np.round(h.getSolution().col_value[:3], 4)}")
    print(f"Schattenpreise: {np.round(h.getSolution().row_dual[:2], 4)}")
elif status == highspy.HighsModelStatus.kInfeasible:
    print("INFEASIBLE")
elif status == highspy.HighsModelStatus.kUnbounded:
    print("UNBOUNDED")
elif status == highspy.HighsModelStatus.kTimeLimit:
    print(f"Zeitlimit, Gap {h.getInfo().mip_gap:.3f}")
else:
    print("kein Optimum:", h.modelStatusToString(status))

Die drei häufigsten Stolpersteine

Nicht zusammen mit ortools importieren. Beide bringen eine eigene HiGHS-Kopie mit; im selben Prozess endet das in undefined symbol (Anhang C, C9). Auch cvxpy zieht highspy bei der Solver-Erkennung mit hinein.
Das CSR-Format stimmt oder es stimmt still nicht. starts hat so viele Einträge wie Zeilen, indices und values so viele wie Nichtnullen. Ein falscher starts-Eintrag erzeugt ein anderes, aber lösbares Modell — es fällt nur durch ein falsches Ergebnis auf.
output_flag abschalten, sonst überschwemmt das Solverprotokoll jede Ausgabe.

D3 — OR-Tools: `pywraplp`

Wofür. Bequeme algebraische Schreibweise für LP und MILP mit umschaltbarem Backend (GLOP, SCIP, CBC, SAT). Wofür nicht: Scheduling und kombinatorische Bedingungen — dafür ist CP-SAT (D4) da.

from ortools.linear_solver import pywraplp

# "GLOP" = LP, "SCIP" oder "CBC" = MILP, "SAT" = CP-SAT als MILP-Backend
loeser = pywraplp.Solver.CreateSolver("GLOP")
if loeser is None:
    raise SystemExit("Solver nicht verfuegbar")
loeser.SetTimeLimit(60_000)                        # Millisekunden!

unendlich = loeser.infinity()
x = [loeser.NumVar(0, unendlich, f"x{j}") for j in range(3)]
# ganzzahlig: loeser.IntVar(0, unendlich, "n") | binaer: loeser.BoolVar("b")

loeser.Add(x[0] + x[1] + 2 * x[2] <= 40)
loeser.Add(2 * x[0] + 3 * x[1] + x[2] <= 50)
loeser.Maximize(10 * x[0] + 15 * x[1] + 25 * x[2])

status = loeser.Solve()
if status == pywraplp.Solver.OPTIMAL:
    print(f"optimal: {loeser.Objective().Value():.2f}")
    print(f"x = {[round(v.solution_value(), 4) for v in x]}")
elif status == pywraplp.Solver.FEASIBLE:
    print("zulaessig, nicht bewiesen optimal (Zeitlimit)")
elif status == pywraplp.Solver.INFEASIBLE:
    print("INFEASIBLE")
elif status == pywraplp.Solver.UNBOUNDED:
    print("UNBOUNDED")
else:
    print("ABNORMAL / NOT_SOLVED - Modell oder Solver pruefen")

Die drei häufigsten Stolpersteine

CreateSolver liefert None, wenn der Backend-Name unbekannt oder nicht gebaut ist — immer prüfen, statt am None später zu scheitern.
SetTimeLimit erwartet Millisekunden, nicht Sekunden. Ein SetTimeLimit(60) bricht nach einer sechzigstel Sekunde ab.
FEASIBLE ist kein OPTIMAL. Bei Zeitlimit liefert der Solver eine gültige, aber möglicherweise schlechte Lösung — den Gap mitberichten, nicht die Zahl allein.

D4 — CP-SAT: `cp_model`

Wofür. Scheduling, Zuordnung, Reihenfolgen, alles Kombinatorische mit globalen Bedingungen. Wofür nicht: stetige Größen — CP-SAT rechnet ausschließlich ganzzahlig. Wer Nachkommastellen braucht, skaliert (Cent statt Euro, Promille statt Anteil).

from ortools.sat.python import cp_model

modell = cp_model.CpModel()

# Variablen - CP-SAT rechnet ausschliesslich mit GANZEN Zahlen
x = modell.NewIntVar(0, 100, "x")                  # untere, obere Schranke, Name
y = modell.NewIntVar(0, 100, "y")
b = modell.NewBoolVar("b")                         # 0/1

# Nebenbedingungen
modell.Add(2 * x + 3 * y <= 50)
modell.Add(x >= 5).OnlyEnforceIf(b)                # gilt nur, wenn b wahr ist
modell.AddAllDifferent([x, y])                     # globale Bedingung
modell.AddMaxEquality(z := modell.NewIntVar(0, 100, "z"), [x, y])

modell.Maximize(10 * x + 15 * y - 3 * z)

loeser = cp_model.CpSolver()
loeser.parameters.max_time_in_seconds = 10.0
loeser.parameters.num_workers = 1                  # 1 = reproduzierbar
loeser.parameters.random_seed = 1
status = loeser.Solve(modell)

if status == cp_model.OPTIMAL:
    print(f"optimal: {loeser.ObjectiveValue():.0f}  x={loeser.Value(x)} y={loeser.Value(y)}")
elif status == cp_model.FEASIBLE:
    print(f"zulaessig, nicht bewiesen optimal - Gap-Schranke: {loeser.BestObjectiveBound()}")
elif status == cp_model.INFEASIBLE:
    print("INFEASIBLE - Bedingungen widersprechen sich")
elif status == cp_model.MODEL_INVALID:
    print("Modellfehler:", modell.Validate())
else:
    print("UNKNOWN - Zeit abgelaufen, ohne eine Loesung zu finden")
print(f"Laufzeit {loeser.WallTime():.3f}s, {loeser.NumBranches()} Verzweigungen")

Die drei häufigsten Stolpersteine

Alles ist ganzzahlig. 0.5 * x gibt es nicht. Skalieren Sie das ganze Modell auf eine feinere Einheit, statt zu runden.
UNKNOWN heißt nicht INFEASIBLE. Es heißt: Die Zeit war zu knapp, um überhaupt etwas zu finden. Die beiden zu verwechseln ist einer der teuersten Fehler in Produktion — das Modell wird für widersprüchlich erklärt, obwohl es lösbar ist.
Ohne num_workers = 1 ist der Lauf nicht reproduzierbar. Mehrere Suchstränge finden je nach Zeitverlauf verschiedene, gleich gute Lösungen. Für Tests und für abgedruckte Ausgaben Worker und Seed festnageln.

D5 — CVXPY

Wofür. Konvexe Probleme: quadratische Ziele, Normen, CVaR, alles mit Regularisierungstermen. Wofür nicht: große kombinatorische Modelle — der Aufbau der Ausdrücke wird dann selbst zum Engpass.

import numpy as np
import cvxpy as cp

x = cp.Variable(3, nonneg=True)                    # nonneg=True statt x >= 0
gewichte = cp.Variable(3)
ganzzahlig = cp.Variable(3, integer=True)          # macht daraus ein MIP

ziel = cp.Maximize(np.array([10.0, 15.0, 25.0]) @ x)
bedingungen = [np.array([[1, 1, 2], [2, 3, 1]]) @ x <= np.array([40.0, 50.0])]

problem = cp.Problem(ziel, bedingungen)
if not problem.is_dcp():                           # VOR dem Loesen pruefen
    print("nicht DCP:", [c for c in bedingungen if not c.is_dcp()])
problem.solve()

if problem.status == cp.OPTIMAL:
    print(f"optimal: {problem.value:.2f}  x = {np.round(x.value, 4)}")
    print(f"Schattenpreis: {np.round(bedingungen[0].dual_value, 4)}")
elif problem.status == cp.OPTIMAL_INACCURATE:
    print("Loesung numerisch unsicher - Skalierung pruefen, anderen Solver testen")
elif problem.status == cp.INFEASIBLE:
    print("INFEASIBLE")
elif problem.status == cp.UNBOUNDED:
    print("UNBOUNDED")
else:
    print("Solverfehler:", problem.status)

Die drei häufigsten Stolpersteine

DCP-Regeln vor dem Lösen prüfen. problem.is_dcp() und dann die einzelnen Bedingungen — das nennt die Schuldige, statt einen DCPError ohne Ort zu werfen (Anhang C, C8).
OPTIMAL_INACCURATE ist kein Erfolg. Der Solver hat aufgegeben und meldet das leise. Diesen Fall immer eigens behandeln.
Der Aufbau kann teurer sein als das Lösen. Schleifen über Szenarien durch Vektorausdrücke ersetzen; Parameter statt Neuaufbau, wenn sich nur Zahlen ändern.

D6 — Modellierungssprachen: Pyomo und Linopy

Wofür. Beide trennen Modell von Solver: dasselbe Modell läuft ohne Änderung unter HiGHS, CBC, Gurobi. Pyomo denkt in Mengen und Indizes wie eine mathematische Formulierung; Linopy denkt in beschrifteten Arrays und baut Nebenbedingungen als Matrixoperation statt in Python-Schleifen. Wofür nicht: ein Modell mit zehn Nebenbedingungen — dort ist der Aufwand höher als der Nutzen.

Pyomo

import pyomo.environ as pyo

modell = pyo.ConcreteModel()
modell.J = pyo.RangeSet(0, 2)                      # Indexmenge
modell.x = pyo.Var(modell.J, domain=pyo.NonNegativeReals)
# ganzzahlig: domain=pyo.NonNegativeIntegers | binaer: domain=pyo.Binary

ertrag = {0: 10.0, 1: 15.0, 2: 25.0}
modell.ziel = pyo.Objective(expr=sum(ertrag[j] * modell.x[j] for j in modell.J),
                            sense=pyo.maximize)
modell.montage = pyo.Constraint(
    expr=modell.x[0] + modell.x[1] + 2 * modell.x[2] <= 40)
modell.pruefung = pyo.Constraint(
    expr=2 * modell.x[0] + 3 * modell.x[1] + modell.x[2] <= 50)
modell.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT)   # fuer Schattenpreise

ergebnis = pyo.SolverFactory("appsi_highs").solve(modell)
zustand = ergebnis.solver.termination_condition
if zustand == pyo.TerminationCondition.optimal:
    print(f"optimal: {pyo.value(modell.ziel):.2f}")
    print(f"x = {[round(pyo.value(modell.x[j]), 4) for j in modell.J]}")
    print(f"Schattenpreis Montage: {modell.dual[modell.montage]:.4f}")
elif zustand == pyo.TerminationCondition.infeasible:
    print("INFEASIBLE")
elif zustand == pyo.TerminationCondition.unbounded:
    print("UNBOUNDED")
elif zustand == pyo.TerminationCondition.maxTimeLimit:
    print("Zeitlimit")
else:
    print("kein Optimum:", zustand)

Linopy

import linopy
import numpy as np
import pandas as pd
import xarray as xr

# Benannte Indizes statt blosser Listen - sonst heissen die Achsen "dim_0"
produkt = pd.Index(["Rahmen", "Gehaeuse", "Deckel"], name="produkt")
ressource = pd.Index(["Montage", "Pruefung"], name="ressource")

modell = linopy.Model()
modell.add_variables(lower=0, coords=[produkt], name="menge")
x = modell.variables["menge"]
# ganzzahlig: integer=True | binaer: binary=True

ertrag = xr.DataArray([10.0, 15.0, 25.0], coords=[produkt])
verbrauch = xr.DataArray([[1.0, 1.0, 2.0], [2.0, 3.0, 1.0]],
                         coords=[ressource, produkt])
vorrat = xr.DataArray([40.0, 50.0], coords=[ressource])

modell.add_objective((ertrag * x).sum(), sense="max")
# EINE Zeile, so viele Nebenbedingungen wie 'ressource' Eintraege hat:
modell.add_constraints((verbrauch * x).sum("produkt") <= vorrat, name="kapazitaet")

modell.solve(solver_name="highs", output_flag=False)

if modell.termination_condition == "optimal":
    print(f"optimal: {modell.objective.value:.2f}")
    print(x.solution.to_series().round(4).to_dict())
    print("Schattenpreise:",
          modell.constraints["kapazitaet"].dual.to_series().round(4).to_dict())
elif modell.termination_condition == "infeasible":
    print("INFEASIBLE")
elif modell.termination_condition == "unbounded":
    print("UNBOUNDED")
else:
    print("kein Optimum:", modell.status, modell.termination_condition)

Die drei häufigsten Stolpersteine

Der Solver ist ein eigenes Programm. SolverFactory("appsi_highs") scheitert, wenn HiGHS nicht auffindbar ist — die Fehlermeldung nennt dann das Modell, nicht die fehlende Installation.
Schattenpreise kommen nur auf Anforderung. Bei Pyomo braucht es Suffix(direction=IMPORT) vor dem Lösen; wer ihn vergisst, bekommt einen KeyError statt einer Warnung.
Bei Linopy die Indizes benennen (pd.Index(..., name="produkt")). Ohne Namen heißen die Achsen dim_0, und jede spätere Zuordnung wird zum Ratespiel.

Die gemeinsame Regel

Alle sechs Seiten haben denselben längsten Abschnitt: die Statusauswertung. Das ist kein Zufall. Ein Solveraufruf hat nie zwei Ausgänge, sondern mindestens fünf — optimal, zulässig ohne Beweis, unlösbar, unbeschränkt, abgebrochen. Code, der nur if erfolgreich: prüft, verwechselt früher oder später „keine Lösung gefunden“ mit „es gibt keine Lösung“, und diese Verwechslung merkt niemand, bis sie teuer wird.

Anhang D: Spickzettel der Solver

D1 — SciPy: linprog und milp

Lineares Programm

Ganzzahlig: milp

Die drei häufigsten Stolpersteine

D2 — HiGHS über highspy

Die drei häufigsten Stolpersteine

D3 — OR-Tools: pywraplp

Die drei häufigsten Stolpersteine

D4 — CP-SAT: cp_model

Die drei häufigsten Stolpersteine

D5 — CVXPY

Die drei häufigsten Stolpersteine

D6 — Modellierungssprachen: Pyomo und Linopy

Pyomo

Linopy

Die drei häufigsten Stolpersteine

Die gemeinsame Regel

D1 — SciPy: `linprog` und `milp`

Ganzzahlig: `milp`

D2 — HiGHS über `highspy`

D3 — OR-Tools: `pywraplp`

D4 — CP-SAT: `cp_model`