bunkerfix

src/optimization/model_builder.py

@@ -120,6 +120,33 @@ def build_model(
 ) -> pyo.ConcreteModel:
     if step_size_tonnes <= 0:
         raise ValueError("step_size_tonnes must be positive")
+
+    # =========================================================================
+    # INHALTSVERZEICHNIS build_model()
+    # =========================================================================
+    #   1  Datenvorbereitung          Bedarf, Bounds, Wochenstruktur
+    #   2  Sets & Parameter           I, J, W, D, S; Bedarfs-/Toleranzparameter
+    #   3  Entscheidungsvariablen     k (ganzzahlig), x = step_size * k
+    #   4  Bunkerlogik                Bunker-Vars und Bilanzen (optional)
+    #   5  Aggregationen              y_delivery_shift, y_delivery, y_week, y_month
+    #   6  Liefertoleranzen KW        täglich, wöchentlich, monatlich, gesamt
+    #   7  Veredlungstoleranzen       täglich, wöchentlich, monatlich (NO/WZ/ALL)
+    #   8  Mischungsgrenzen           alpha hard (mix_lower/upper), soft (mix_target)
+    #   9  Abweichungsverfolgung      dev, z_max, no_three_in_a_row, devV
+    #  10  Schichtglättung            tagesübergreifend: Spätschicht Mi–Di, Früh-/Nacht Mo–Do / Fr–So
+    #  11  Zielfunktion               Strafterme und Lambda-Gewichte
+    #  12  Routenverbote              Reichwalde->V, Welzow->B4
+    #  13  Förder-/Verladekapazitäten monatlich, Schicht, Tag
+    #  14  Verfügbarkeiten            dynamische Caps aus Verfuegbarkeiten.parquet
+    #  15  Zugdurchlass               Streckencaps aus zugdurchlass.parquet
+    # =========================================================================
+
+    # =========================================================================
+    # 1  DATENVORBEREITUNG
+    #    Bedarf-DataFrame aus Excel aufbereiten, Bergbauwoche berechnen,
+    #    Bounds für Kraftwerke und Veredlung aus den Parquet-Tabellen laden.
+    # =========================================================================
+
     bedarf = tables["rohkohlebedarf"][
         [
             "datum",
@@ -153,6 +180,12 @@ def build_model(
     bounds_power_plants = tables["bounds_power_plants"]
     bounds_day = bounds_power_plants[bounds_power_plants["zeitraum"] == "pro Tag"].copy()
 
+    # =========================================================================
+    # 2  SETS & PARAMETER
+    #    Pyomo-Modell anlegen, Mengen I/J/W/D/S definieren, Bedarfs- und
+    #    Tagestoleranzparameter initialisieren und aus den Tabellen befüllen.
+    # =========================================================================
+
     model = pyo.ConcreteModel()
 
     J_POWER = ["J", "SP", "B3", "B4"]
@@ -244,6 +277,12 @@ def build_model(
             model.a_min_day[j, w, d] = a_min
             model.a_max_day[j, w, d] = a_max
 
+    # =========================================================================
+    # 3  ENTSCHEIDUNGSVARIABLEN
+    #    k[i,j,w,d,s] als ganzzahlige Lieferschritte, x = step_size * k
+    #    als abgeleiteter kontinuierlicher Fluss in Tonnen.
+    # =========================================================================
+
     model.step_size_tonnes = pyo.Param(initialize=float(step_size_tonnes), mutable=False)
     model.k = pyo.Var(model.I, model.J, model.W, model.D, model.S, domain=pyo.NonNegativeIntegers)
 
@@ -252,6 +291,12 @@ def build_model(
 
     model.x = pyo.Expression(model.I, model.J, model.W, model.D, model.S, rule=x_rule)
 
+    # =========================================================================
+    # 4  BUNKERLOGIK
+    #    Falls eine Bunkertabelle vorliegt: Bunker-Vars, Ausfluss-Expressions
+    #    und tägliche Bilanzbedingungen anlegen; sonst out = x direkt.
+    # =========================================================================
+
     bunker_df = tables.get("bunker")
     if bunker_df is not None:
         j_bunker_map = {
@@ -357,6 +402,21 @@ def build_model(
                         + sum(model.x[i, j, w, d, s] for s in model.S)
                         - sum(model.bunker_out[i, j, w, d, s] for s in model.S)
                     )
+
+        # Nebenbedingung: Täglicher Bunkerabfluss entspricht der Kraftwerksnachfrage.
+        # Damit wird bunker_out physikalisch korrekt an den Verbrauch gebunden.
+        model.bunker_out_demand = pyo.ConstraintList()
+        for j in model.J_BUNKER:
+            for w, d in time_points_day:
+                demand_expr = (
+                    model.dV_N[w, d] + model.dV_W[w, d]
+                    if j == "V"
+                    else model.d[j, w, d]
+                )
+                model.bunker_out_demand.add(
+                    sum(model.bunker_out[i, j, w, d, s] for i in model.I for s in model.S)
+                    == demand_expr
+                )
     else:
         def out_rule(m, i, j, w, d, s):
             return m.x[i, j, w, d, s]
@@ -373,6 +433,12 @@ def build_model(
 
         model.y = pyo.Expression(model.J, model.W, model.D, rule=y_rule)
 
+    # =========================================================================
+    # 5  AGGREGATIONEN
+    #    Expressions für Schicht-, Tages-, Wochen- und Monatsflüsse aufbauen,
+    #    die in Constraints und Zielfunktion wiederverwendet werden.
+    # =========================================================================
+
     def y_delivery_shift_rule(model, j, w, d, s):
         return sum(model.x[i, j, w, d, s] for i in model.I)
 
@@ -383,6 +449,12 @@ def build_model(
 
     model.y_delivery = pyo.Expression(model.J, model.W, model.D, rule=y_delivery_rule)
 
+    # =========================================================================
+    # 6  LIEFERTOLERANZEN KRAFTWERKE
+    #    Tägliche, wöchentliche und monatliche Toleranzbänder je Kraftwerk
+    #    sowie Gesamtsystemgrenzen über alle Werke.
+    # =========================================================================
+
     J_delivery = [j for j in model.J if j != "V"]
 
     def delivery_tolerance_rule(model, j, w, d):
@@ -392,24 +464,9 @@ def build_model(
             model.d[j, w, d] + model.a_max_day[j, w, d],
         )
 
-    # Constraint: Daily delivery must stay within demand tolerance for power plants.
+    # Nebenbedingung: Tägliche Lieferung muss innerhalb des Toleranzbandes der Kraftwerke liegen.
     model.delivery_tolerance = pyo.Constraint(J_delivery, model.W, model.D, rule=delivery_tolerance_rule)
 
-    model.shift_balance_dev = pyo.Var(model.J, model.W, model.D, model.S, domain=pyo.NonNegativeReals)
-    model.shift_balance_abs = pyo.ConstraintList()
-    for j in model.J:
-        for w in model.W:
-            for d in model.D:
-                for s in model.S:
-                    model.shift_balance_abs.add(
-                        model.shift_balance_dev[j, w, d, s]
-                        >= model.y_delivery_shift[j, w, d, s] - model.y_delivery[j, w, d] / 3
-                    )
-                    model.shift_balance_abs.add(
-                        model.shift_balance_dev[j, w, d, s]
-                        >= model.y_delivery[j, w, d] / 3 - model.y_delivery_shift[j, w, d, s]
-                    )
-
     week_to_month = bedarf.groupby("week")["datum"].min().dt.month.to_dict()
     model.M = pyo.Set(initialize=sorted(set(week_to_month.values())))
     J_power = [j for j in model.J if j != "V"]
@@ -490,7 +547,7 @@ def build_model(
             model.d_week[j, w] + model.a_max_week[j, w],
         )
 
-    # Constraint: Weekly deliveries per plant stay within weekly tolerance band.
+    # Nebenbedingung: Wöchentliche Lieferung je Kraftwerk innerhalb des Wochentoleranzbandes.
     model.week_tolerance = pyo.Constraint(J_power, model.W, rule=week_tolerance_rule)
 
     def month_tolerance_rule(model, j, m):
@@ -500,7 +557,7 @@ def build_model(
             model.d_month[j, m] + model.a_max_month[j, m],
         )
 
-    # Constraint: Monthly deliveries per plant stay within monthly tolerance band.
+    # Nebenbedingung: Monatliche Lieferung je Kraftwerk innerhalb des Monatstoleranzbandes.
     model.month_tolerance = pyo.Constraint(J_power, model.M, rule=month_tolerance_rule)
 
     def week_total_rule(model, w):
@@ -510,7 +567,7 @@ def build_model(
             model.d_week_total[w] + model.a_max_week_total[w],
         )
 
-    # Constraint: Total weekly deliveries across all plants stay within tolerance band.
+    # Nebenbedingung: Wöchentliche Gesamtlieferung über alle Kraftwerke innerhalb des Toleranzbandes.
     model.week_total_tolerance = pyo.Constraint(model.W, rule=week_total_rule)
 
     def month_total_rule(model, m):
@@ -520,9 +577,15 @@ def build_model(
             model.d_month_total[m] + model.a_max_month_total[m],
         )
 
-    # Constraint: Total monthly deliveries across all plants stay within tolerance band.
+    # Nebenbedingung: Monatliche Gesamtlieferung über alle Kraftwerke innerhalb des Toleranzbandes.
     model.month_total_tolerance = pyo.Constraint(model.M, rule=month_total_rule)
 
+    # =========================================================================
+    # 7  VEREDLUNGSTOLERANZEN
+    #    Tages-, Wochen- und Monatsgrenzen für Nochten-, Welzow- und
+    #    Gesamtkohle an die Veredlung (ISP).
+    # =========================================================================
+
     model.a_min_V_day = pyo.Param(["N", "W", "ALL"], model.W, model.D, initialize=0.0, mutable=True)
     model.a_max_V_day = pyo.Param(["N", "W", "ALL"], model.W, model.D, initialize=0.0, mutable=True)
     model.a_min_V_week = pyo.Param(["N", "W", "ALL"], model.W, initialize=0.0, mutable=True)
@@ -581,7 +644,7 @@ def build_model(
             m.dV_W[w, d] + m.a_max_V_day["W", w, d],
         )
 
-    # Constraint: Veredlung day tolerance for Welzow coal.
+    # Nebenbedingung: Tagestoleranz Veredlung für Welzower Kohle.
     model.v_w_day_tol = pyo.Constraint(model.W, model.D, rule=v_w_day_rule)
 
     def v_n_day_rule(m, w, d):
@@ -591,7 +654,7 @@ def build_model(
             m.dV_N[w, d] + m.a_max_V_day["N", w, d],
         )
 
-    # Constraint: Veredlung day tolerance for Nochten coal.
+    # Nebenbedingung: Tagestoleranz Veredlung für Nochtener Kohle.
     model.v_n_day_tol = pyo.Constraint(model.W, model.D, rule=v_n_day_rule)
 
     def v_all_day_rule(m, w, d):
@@ -601,7 +664,7 @@ def build_model(
             m.dV_W[w, d] + m.dV_N[w, d] + m.a_max_V_day["ALL", w, d],
         )
 
-    # Constraint: Veredlung day tolerance for combined coal types.
+    # Nebenbedingung: Tagestoleranz Veredlung für kombinierten Bedarf (NO+WZ).
     model.v_all_day_tol = pyo.Constraint(model.W, model.D, rule=v_all_day_rule)
 
     def v_w_week_rule(m, w):
@@ -611,7 +674,7 @@ def build_model(
             sum(m.dV_W[w, d] for d in m.D) + m.a_max_V_week["W", w],
         )
 
-    # Constraint: Veredlung weekly tolerance for Welzow coal.
+    # Nebenbedingung: Wochentoleranz Veredlung für Welzower Kohle.
     model.v_w_week_tol = pyo.Constraint(model.W, rule=v_w_week_rule)
 
     def v_n_week_rule(m, w):
@@ -621,7 +684,7 @@ def build_model(
             sum(m.dV_N[w, d] for d in m.D) + m.a_max_V_week["N", w],
         )
 
-    # Constraint: Veredlung weekly tolerance for Nochten coal.
+    # Nebenbedingung: Wochentoleranz Veredlung für Nochtener Kohle.
     model.v_n_week_tol = pyo.Constraint(model.W, rule=v_n_week_rule)
 
     def v_all_week_rule(m, w):
@@ -631,7 +694,7 @@ def build_model(
             sum(m.dV_W[w, d] + m.dV_N[w, d] for d in m.D) + m.a_max_V_week["ALL", w],
         )
 
-    # Constraint: Veredlung weekly tolerance for combined coal types.
+    # Nebenbedingung: Wochentoleranz Veredlung für kombinierten Bedarf (NO+WZ).
     model.v_all_week_tol = pyo.Constraint(model.W, rule=v_all_week_rule)
 
     def v_w_month_rule(m, mm):
@@ -642,7 +705,7 @@ def build_model(
             sum(m.dV_W[w, d] for w in weeks_in_m for d in m.D) + m.a_max_V_month["W", mm],
         )
 
-    # Constraint: Veredlung monthly tolerance for Welzow coal.
+    # Nebenbedingung: Monatstoleranz Veredlung für Welzower Kohle.
     model.v_w_month_tol = pyo.Constraint(model.M, rule=v_w_month_rule)
 
     def v_n_month_rule(m, mm):
@@ -653,7 +716,7 @@ def build_model(
             sum(m.dV_N[w, d] for w in weeks_in_m for d in m.D) + m.a_max_V_month["N", mm],
         )
 
-    # Constraint: Veredlung monthly tolerance for Nochten coal.
+    # Nebenbedingung: Monatstoleranz Veredlung für Nochtener Kohle.
     model.v_n_month_tol = pyo.Constraint(model.M, rule=v_n_month_rule)
 
     def v_all_month_rule(m, mm):
@@ -664,7 +727,7 @@ def build_model(
             sum(m.dV_W[w, d] + m.dV_N[w, d] for w in weeks_in_m for d in m.D) + m.a_max_V_month["ALL", mm],
         )
 
-    # Constraint: Veredlung monthly tolerance for combined coal types.
+    # Nebenbedingung: Monatstoleranz Veredlung für kombinierten Bedarf (NO+WZ).
     model.v_all_month_tol = pyo.Constraint(model.M, rule=v_all_month_rule)
 
     j_map = {
@@ -682,6 +745,12 @@ def build_model(
 
     mix_df = tables["kohle_mix"]
 
+    # =========================================================================
+    # 8  MISCHUNGSGRENZEN
+    #    Harte Mindest- und Maximalanteile je Quelle/Ziel (alpha_min/max),
+    #    plus weiche Zielbandabweichungen (mix_target_low/high/mid).
+    # =========================================================================
+
     model.alpha_min = pyo.Param(model.I, model.J, initialize=0.0, mutable=True)
     model.alpha_max = pyo.Param(model.I, model.J, initialize=1.0, mutable=True)
 
@@ -714,7 +783,7 @@ def build_model(
             return pyo.Constraint.Skip
         return m.x_day[i, j, w, d] <= m.alpha_max[i, j] * m.y_day[j, w, d]
 
-    # Constraint: Mix lower/upper bound per source and plant (hard), per day.
+    # Nebenbedingung: Harte Mischungsunter- und -obergrenzen je Quelle und Kraftwerk, täglich.
     model.mix_lower = pyo.Constraint(model.I, model.J, model.W, model.D, rule=mix_lower_rule)
     model.mix_upper = pyo.Constraint(model.I, model.J, model.W, model.D, rule=mix_upper_rule)
 
@@ -746,7 +815,7 @@ def build_model(
             return pyo.Constraint.Skip
         return m.x_day[i, j, w, d] <= m.alpha_target_high[i, j] * m.y_day[j, w, d] + m.mix_target_high_dev[i, j, w, d]
 
-    # Constraint: Mix target lower/upper band (soft via deviation variable).
+    # Nebenbedingung: Weiches Mischungszielband je Quelle und Kraftwerk (Abweichungsvariable).
     model.mix_target_low = pyo.Constraint(model.I, model.J, model.W, model.D, rule=mix_target_low_rule)
     model.mix_target_high = pyo.Constraint(model.I, model.J, model.W, model.D, rule=mix_target_high_rule)
 
@@ -835,6 +904,12 @@ def build_model(
                         >= model.bunker_target[j] - model.bunker_total[j, w_prev, d_prev]
                     )
 
+    # =========================================================================
+    # 9  ABWEICHUNGSVERFOLGUNG
+    #    Tagesabweichung in pos/neg aufteilen, Binärvariable für Maximalabweichung
+    #    setzen und die Drei-Tage-Regel in Folge begrenzen.
+    # =========================================================================
+
     def demand_rule(m, j, w, d):
         return m.dV_N[w, d] + m.dV_W[w, d] if j == "V" else m.d[j, w, d]
 
@@ -850,7 +925,7 @@ def build_model(
     def dev_balance_rule(model, j, w, d):
         return model.dev_pos[j, w, d] - model.dev_neg[j, w, d] == model.dev[j, w, d]
 
-    # Constraint: Split demand deviation into positive/negative components.
+    # Nebenbedingung: Nachfrageabweichung in positiven und negativen Anteil aufteilen.
     model.dev_balance = pyo.Constraint(J_dev, model.W, model.D, rule=dev_balance_rule)
 
     model.z_max = pyo.Var(J_power, model.W, model.D, domain=pyo.Binary)
@@ -858,7 +933,7 @@ def build_model(
     def max_reached_rule(m, j, w, d):
         return m.dev_pos[j, w, d] <= m.a_max_day[j, w, d] * m.z_max[j, w, d]
 
-    # Constraint: Flag when daily deviation hits the maximum tolerance.
+    # Nebenbedingung: Binärvariable setzen, wenn Tagesabweichung die Maximumtoleranz erreicht.
     model.max_reached = pyo.Constraint(J_power, model.W, model.D, rule=max_reached_rule)
 
     D_list = list(model.D.ordered_data())
@@ -867,7 +942,7 @@ def build_model(
         d1, d2, d3 = D_list[idx : idx + 3]
         return m.z_max[j, w, d1] + m.z_max[j, w, d2] + m.z_max[j, w, d3] <= 3
 
-    # Constraint: Prevent three consecutive days with max deviation.
+    # Nebenbedingung: Mehr als drei aufeinanderfolgende Tage mit Maximalabweichung verhindern.
     model.no_three_in_a_row = pyo.Constraint(J_power, model.W, range(len(D_list) - 2), rule=no_three_in_a_row_rule)
 
     def y_V_nochten(m, w, d):
@@ -884,52 +959,95 @@ def build_model(
     model.devV_W_pos = pyo.Var(model.W, model.D, domain=pyo.NonNegativeReals)
     model.devV_W_neg = pyo.Var(model.W, model.D, domain=pyo.NonNegativeReals)
 
-    # Constraint: Split Nochten Veredlung deviation into positive/negative parts.
+    # Nebenbedingung: Veredlungsabweichung Nochten in positiven und negativen Anteil aufteilen.
     model.devV_N_balance = pyo.Constraint(
         model.W, model.D, rule=lambda m, w, d: m.devV_N_pos[w, d] - m.devV_N_neg[w, d] == m.devV_N[w, d]
     )
-    # Constraint: Split Welzow Veredlung deviation into positive/negative parts.
+    # Nebenbedingung: Veredlungsabweichung Welzow in positiven und negativen Anteil aufteilen.
     model.devV_W_balance = pyo.Constraint(
         model.W, model.D, rule=lambda m, w, d: m.devV_W_pos[w, d] - m.devV_W_neg[w, d] == m.devV_W[w, d]
     )
 
+    # =========================================================================
+    # 10  SCHICHTGLÄTTUNG
+    #     Tagesübergreifende Gleichmäßigkeit je Schicht innerhalb fixer Zeitfenster:
+    #       Spätschicht:       Mi–Di  (7 Tage, wochenübergreifend)
+    #       Früh-/Nachtschicht: Mo–Do (innerhalb Bergbauwoche) und Fr–So (Wochenübergang)
+    #     Soft-Constraint: Paarweise Tagesdifferenz wird ab SHIFT_SMOOTH_TOL bestraft.
+    # =========================================================================
 
-    lambda_dev = 100_000
+    SHIFT_SMOOTH_TOL = step_size_tonnes  # tolerierte Tagesdifferenz je Schicht in Tonnen
+
+    lambda_dev = 100_000_000
     lambda_v = 100_000
-    lambda_shift = 100_000_000
+    lambda_shift_smooth = 100_000_000
     lambda_mix = 4_000_000_000
     lambda_mix_mid = 4_000_000
     lambda_b3_bunker_mix = 100_000_000
-    lambda_shift_balance = 5_000_000
     lambda_bunker_target = 50_000_000
-    SHIFT_TOL = step_size_tonnes
-    SHIFT_PAIRS = [("F", "S"), ("S", "N"), ("F", "N")]
-    model.SHIFT_PAIRS = pyo.Set(initialize=SHIFT_PAIRS, dimen=2)
 
-    model.shift_dev_soft = pyo.Var(model.I, model.J, model.W, model.D, model.SHIFT_PAIRS, domain=pyo.NonNegativeReals)
+    # Nachfolgewoche je Bergbauwoche (für wochenübergreifende Fenster)
-    model.shift_excess = pyo.Var(model.I, model.J, model.W, model.D, model.SHIFT_PAIRS, domain=pyo.NonNegativeReals)
+    nw_map = {weeks[i]: weeks[i + 1] for i in range(len(weeks) - 1)}
+
+    # Alle paarweisen Tages-Kombinationen je Schicht-Fenster sammeln.
+    # Einträge: (schicht, woche_a, tag_a, woche_b, tag_b)
+    smooth_pairs: list[tuple] = []
+
+    for w in weeks:
+        nw = nw_map.get(w)
+
+        # Spätschicht: Mittwoch bis Dienstag der Folgewoche (7 Tage)
+        sp = [(w, "Mi"), (w, "Do"), (w, "Fr")]
+        if nw is not None:
+            sp += [(nw, "Sa"), (nw, "So"), (nw, "Mo"), (nw, "Di")]
+        for a_idx, (wa, da) in enumerate(sp):
+            for wb, db in sp[a_idx + 1 :]:
+                smooth_pairs.append(("S", wa, da, wb, db))
+
+        for s in ("F", "N"):
+            # Gruppe 1: Montag bis Donnerstag (innerhalb einer Bergbauwoche)
+            g1 = [(w, "Mo"), (w, "Di"), (w, "Mi"), (w, "Do")]
+            for a_idx, (wa, da) in enumerate(g1):
+                for wb, db in g1[a_idx + 1 :]:
+                    smooth_pairs.append((s, wa, da, wb, db))
+
+            # Gruppe 2: Freitag bis Sonntag (Wochenübergang)
+            g2 = [(w, "Fr")]
+            if nw is not None:
+                g2 += [(nw, "Sa"), (nw, "So")]
+            for a_idx, (wa, da) in enumerate(g2):
+                for wb, db in g2[a_idx + 1 :]:
+                    smooth_pairs.append((s, wa, da, wb, db))
+
+    model.SMOOTH_PAIR_IDX = pyo.Set(initialize=range(len(smooth_pairs)), dimen=1)
+    model.shift_smooth_dev = pyo.Var(model.I, model.J, model.SMOOTH_PAIR_IDX, domain=pyo.NonNegativeReals)
+    model.shift_smooth_excess = pyo.Var(model.I, model.J, model.SMOOTH_PAIR_IDX, domain=pyo.NonNegativeReals)
+
+    # Nebenbedingung: Absolute Tagesdifferenz je Schicht-Fenster-Paar und Toleranzüberschreitung.
+    model.shift_smooth_abs = pyo.ConstraintList()
+    model.shift_smooth_excess_def = pyo.ConstraintList()
 
-    # ConstraintList: Enforce absolute shift differences and excess over tolerance.
-    model.shift_dev_abs = pyo.ConstraintList()
-    model.shift_excess_def = pyo.ConstraintList()
     for i in model.I:
         for j in model.J:
-            for w in model.W:
+            for p_idx, (s, w_a, d_a, w_b, d_b) in enumerate(smooth_pairs):
-                for d in model.D:
+                model.shift_smooth_abs.add(
-                    for s1, s2 in SHIFT_PAIRS:
+                    model.shift_smooth_dev[i, j, p_idx]
-                        if (i, j, w, d, s1) in model.x and (i, j, w, d, s2) in model.x:
+                    >= model.x[i, j, w_a, d_a, s] - model.x[i, j, w_b, d_b, s]
-                            model.shift_dev_abs.add(
+                )
-                                model.shift_dev_soft[i, j, w, d, s1, s2]
+                model.shift_smooth_abs.add(
-                                >= model.x[i, j, w, d, s1] - model.x[i, j, w, d, s2]
+                    model.shift_smooth_dev[i, j, p_idx]
-                            )
+                    >= model.x[i, j, w_b, d_b, s] - model.x[i, j, w_a, d_a, s]
-                            model.shift_dev_abs.add(
+                )
-                                model.shift_dev_soft[i, j, w, d, s1, s2]
+                model.shift_smooth_excess_def.add(
-                                >= model.x[i, j, w, d, s2] - model.x[i, j, w, d, s1]
+                    model.shift_smooth_excess[i, j, p_idx]
-                            )
+                    >= model.shift_smooth_dev[i, j, p_idx] - SHIFT_SMOOTH_TOL
-                            model.shift_excess_def.add(
+                )
-                                model.shift_excess[i, j, w, d, s1, s2]
+
-                                >= model.shift_dev_soft[i, j, w, d, s1, s2] - SHIFT_TOL
+    # =========================================================================
-                            )
+    # 11  ZIELFUNKTION
+    #    Alle Strafterme mit Lambda-Gewichten zusammenführen und als
+    #    Minimierungsziel im Modell registrieren.
+    # =========================================================================
 
     def objective_rule(model):
         deviation_penalty = (
@@ -943,18 +1061,11 @@ def build_model(
             )
         )
 
-        smoothness_penalty = lambda_shift * sum(
+        shift_smooth_penalty = lambda_shift_smooth * sum(
-            model.shift_excess[i, j, w, d, s1, s2]
+            model.shift_smooth_excess[i, j, p_idx]
             for i in model.I
             for j in model.J
-            for w in model.W
+            for p_idx in model.SMOOTH_PAIR_IDX
-            for d in model.D
-            for s1, s2 in model.SHIFT_PAIRS
-            if (i, j, w, d, s1) in model.x and (i, j, w, d, s2) in model.x
-        )
-
-        shift_balance_penalty = lambda_shift_balance * sum(
-            model.shift_balance_dev[j, w, d, s] for j in model.J for w in model.W for d in model.D for s in model.S
         )
 
         mix_target_penalty = lambda_mix * sum(
@@ -985,27 +1096,32 @@ def build_model(
 
         return (
             deviation_penalty
-            + smoothness_penalty
+            + shift_smooth_penalty
-            + shift_balance_penalty
             + mix_target_penalty
             + mix_target_mid_penalty
             + bunker_penalty
             + b3_bunker_mix_penalty
         )
 
-    # Objective: penalize deviations/smoothness, apply route bonuses and penalties.
+    # Zielfunktion: Abweichungen und Ungleichmäßigkeit bestrafen, Routenanreize anwenden.
     model.obj = pyo.Objective(rule=objective_rule, sense=pyo.minimize)
 
+    # =========================================================================
+    # 12  ROUTENVERBOTE
+    #    Fachlich unzulässige Routen hart auf null setzen:
+    #    Reichwalde→Veredlung und Welzow→Boxberg Werk 4.
+    # =========================================================================
+
     def forbid_reichwalde_V_rule(m, w, d, s):
         return m.x["Reichwalde", "V", w, d, s] == 0
 
-    # Constraint: Disallow Reichwalde coal to Veredlung.
+    # Nebenbedingung: Reichwalde-Kohle an Veredlung verboten.
     model.forbid_reichwalde_V = pyo.Constraint(model.W, model.D, model.S, rule=forbid_reichwalde_V_rule)
 
     def forbid_welzow_B4_rule(m, w, d, s):
         return m.x["Welzow", "B4", w, d, s] == 0
 
-    # Constraint: Disallow Welzow coal to Boxberg Werk 4.
+    # Nebenbedingung: Welzow-Kohle an Boxberg Werk 4 verboten.
     model.forbid_welzow_B4 = pyo.Constraint(model.W, model.D, model.S, rule=forbid_welzow_B4_rule)
 
     cap_df = tables["foerderkapaz"]
@@ -1017,6 +1133,12 @@ def build_model(
         "Welzow-Süd": "Welzow",
     }
 
+    # =========================================================================
+    # 13  FÖRDER- UND VERLADEKAPAZITÄTEN
+    #    Monatliche Förderobergrenzen je Tagebau, Verladecaps für Boxberg
+    #    und Welzow-Süd auf Schicht- und Tagesebene.
+    # =========================================================================
+
     model.F_max_month = pyo.Param(model.I, initialize=1e12, mutable=True)
     for _, row in cap_month.iterrows():
         i = i_map.get(row["tagebau"])
@@ -1029,10 +1151,10 @@ def build_model(
 
     model.F_month = pyo.Expression(model.I, model.M, rule=F_month_rule)
 
-    # Constraint: Monthly production cap per mine.
+    # Nebenbedingung: Monatliche Förderobergrenze je Tagebau.
     model.cap_month = pyo.Constraint(model.I, model.M, rule=lambda m, i, mm: m.F_month[i, mm] <= m.F_max_month[i])
 
-    # Constraint: Joint monthly cap for Reichwalde + Nochten.
+    # Nebenbedingung: Gemeinsame monatliche Förderobergrenze für Reichwalde und Nochten.
     model.cap_month_RWNO = pyo.Constraint(
         model.M, rule=lambda m, mm: m.F_month["Reichwalde", mm] + m.F_month["Nochten", mm] <= 3_000_000.0
     )
@@ -1060,7 +1182,7 @@ def build_model(
     def boxberg_shift_cap_rule(m, w, d, s):
         return sum(m.x[i, j, w, d, s] for i in BOXBERG_SOURCES if i in m.I for j in BOXBERG_DEST if j in m.J) <= BOXBERG_SHIFT_CAP
 
-    # Constraint: Boxberg loading capacity per shift.
+    # Nebenbedingung: Verladungskapazität Boxberg je Schicht.
     model.cap_boxberg_shift = pyo.Constraint(model.W, model.D, model.S, rule=boxberg_shift_cap_rule)
 
     def boxberg_day_cap_rule(m, w, d):
@@ -1073,13 +1195,13 @@ def build_model(
             for s in m.S
         ) <= BOXBERG_DAY_CAP
 
-    # Constraint: Boxberg loading capacity per day.
+    # Nebenbedingung: Verladungskapazität Boxberg je Tag.
     model.cap_boxberg_day = pyo.Constraint(model.W, model.D, rule=boxberg_day_cap_rule)
 
     def welzow_shift_cap_rule(m, w, d, s):
         return sum(m.x[i, j, w, d, s] for i in WELZOW_SOURCES if i in m.I for j in WELZOW_DEST if j in m.J) <= WELZOW_SHIFT_CAP
 
-    # Constraint: Welzow loading capacity per shift.
+    # Nebenbedingung: Verladungskapazität Welzow-Süd je Schicht.
     model.cap_welzow_shift = pyo.Constraint(model.W, model.D, model.S, rule=welzow_shift_cap_rule)
 
     def welzow_day_cap_rule(m, w, d):
@@ -1092,9 +1214,15 @@ def build_model(
             for s in m.S
         ) <= WELZOW_DAY_CAP
 
-    # Constraint: Welzow loading capacity per day.
+    # Nebenbedingung: Verladungskapazität Welzow-Süd je Tag.
     model.cap_welzow_day = pyo.Constraint(model.W, model.D, rule=welzow_day_cap_rule)
 
+    # =========================================================================
+    # 14  VERFÜGBARKEITEN
+    #    Schichtweise dynamische Obergrenzen aus Verfuegbarkeiten.parquet
+    #    für Welzow-Süd und den Boxberg-Verbund (RW + NO).
+    # =========================================================================
+
     avail = tables.get("Verfuegbarkeiten")
     if avail is not None:
         day_map = {"Sat": "Sa", "Sun": "So", "Mon": "Mo", "Tue": "Di", "Wed": "Mi", "Thu": "Do", "Fri": "Fr"}
@@ -1142,7 +1270,7 @@ def build_model(
                 <= v
             )
 
-        # Constraint: Dynamic availability cap for Welzow (per shift).
+        # Nebenbedingung: Dynamische Verfügbarkeitsgrenze Welzow-Süd je Schicht.
         model.cap_welzow_con = pyo.Constraint(model.W, model.D, model.S, rule=cap_welzow_rule)
 
         def cap_rw_n_rule(m, w, d, s):
@@ -1159,9 +1287,15 @@ def build_model(
                 <= v
             )
 
-        # Constraint: Dynamic availability cap for Reichwalde + Nochten (per shift).
+        # Nebenbedingung: Dynamische Verfügbarkeitsgrenze Reichwalde + Nochten je Schicht.
         model.cap_rw_n_con = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_rule)
 
+    # =========================================================================
+    # 15  ZUGDURCHLASS
+    #    Streckenbezogene Kapazitätsgrenzen aus zugdurchlass.parquet für
+    #    KLP- und KUP-Routen, inkl. kombinierter Mengenbeschränkungen.
+    # =========================================================================
+
     zug = tables.get("zugdurchlass")
     if zug is not None:
         df = zug.copy()
@@ -1200,31 +1334,31 @@ def build_model(
         def cap_rw_n_to_j_rule(m, w, d, s):
             return m.x["Reichwalde", "J", w, d, s] + m.x["Nochten", "J", w, d, s] <= CAP_RW_N_TO_J
 
-        # Constraint: KLP capacity from RW+NO to J.
+        # Nebenbedingung: KLP-Kapazität Reichwalde/Nochten nach Jänschwalde.
         model.cap_rw_n_to_j = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_to_j_rule)
 
         def cap_rw_n_to_sp_rule(m, w, d, s):
             return m.x["Reichwalde", "SP", w, d, s] + m.x["Nochten", "SP", w, d, s] <= CAP_RW_N_TO_SP
 
-        # Constraint: KLP capacity from RW+NO to SP.
+        # Nebenbedingung: KLP-Kapazität Reichwalde/Nochten nach Schwarze Pumpe.
         model.cap_rw_n_to_sp = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_to_sp_rule)
 
         def cap_rw_n_to_v_rule(m, w, d, s):
             return m.x["Reichwalde", "V", w, d, s] + m.x["Nochten", "V", w, d, s] <= CAP_RW_N_TO_V
 
-        # Constraint: KLP capacity from RW+NO to Veredlung.
+        # Nebenbedingung: KLP-Kapazität Reichwalde/Nochten zur Veredlung.
         model.cap_rw_n_to_v = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_to_v_rule)
 
         def cap_rw_n_to_bw3_rule(m, w, d, s):
             return m.x["Reichwalde", "B3", w, d, s] + m.x["Nochten", "B3", w, d, s] <= CAP_RW_N_TO_B3
 
-        # Constraint: KLP capacity from RW+NO to B3.
+        # Nebenbedingung: KLP-Kapazität Reichwalde/Nochten nach Boxberg Werk 3.
         model.cap_rw_n_to_bw3 = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_to_bw3_rule)
 
         def cap_w_to_j_rule(m, w, d, s):
             return m.x["Welzow", "J", w, d, s] <= CAP_W_TO_J
 
-        # Constraint: KUP capacity from Welzow to J.
+        # Nebenbedingung: KUP-Kapazität Welzow nach Jänschwalde.
         model.cap_w_to_j = pyo.Constraint(model.W, model.D, model.S, rule=cap_w_to_j_rule)
 
         model.k_w_to_j_steps = pyo.Var(model.W, model.D, model.S, domain=pyo.NonNegativeIntegers)
@@ -1232,25 +1366,25 @@ def build_model(
         def welzow_j_multiple_rule(m, w, d, s):
             return m.x["Welzow", "J", w, d, s] == (2 * m.step_size_tonnes) * m.k_w_to_j_steps[w, d, s]
 
-        # Constraint: Enforce 2kt multiples for Welzow->J flows.
+        # Nebenbedingung: Welzow→Jänschwalde-Fluss auf 2.000-t-Vielfache beschränken.
         model.welzow_j_multiple_2kt = pyo.Constraint(model.W, model.D, model.S, rule=welzow_j_multiple_rule)
 
         def cap_w_to_sp_rule(m, w, d, s):
             return m.x["Welzow", "SP", w, d, s] <= CAP_W_TO_SP
 
-        # Constraint: KUP capacity from Welzow to SP.
+        # Nebenbedingung: KUP-Kapazität Welzow nach Schwarze Pumpe.
         model.cap_w_to_sp = pyo.Constraint(model.W, model.D, model.S, rule=cap_w_to_sp_rule)
 
         def cap_w_to_v_rule(m, w, d, s):
             return m.x["Welzow", "V", w, d, s] <= CAP_W_TO_V
 
-        # Constraint: KUP capacity from Welzow to Veredlung.
+        # Nebenbedingung: KUP-Kapazität Welzow zur Veredlung.
         model.cap_w_to_v = pyo.Constraint(model.W, model.D, model.S, rule=cap_w_to_v_rule)
 
         def cap_w_to_bw3_rule(m, w, d, s):
             return m.x["Welzow", "B3", w, d, s] <= CAP_W_TO_B3
 
-        # Constraint: KUP capacity from Welzow to B3.
+        # Nebenbedingung: KUP-Kapazität Welzow nach Boxberg Werk 3.
         model.cap_w_to_bw3 = pyo.Constraint(model.W, model.D, model.S, rule=cap_w_to_bw3_rule)
 
         def cap_rw_n_to_sp_v_rule(m, w, d, s):
@@ -1262,7 +1396,7 @@ def build_model(
                 <= CAP_RW_N_TO_SP_V
             )
 
-        # Constraint: KLP combined capacity to SP + Veredlung.
+        # Nebenbedingung: KLP-Gesamtkapazität nach Schwarze Pumpe und Veredlung.
         model.cap_rw_n_to_sp_v = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_to_sp_v_rule)
 
         def cap_rw_n_to_j_sp_v_b3_rule(m, w, d, s):
@@ -1278,13 +1412,13 @@ def build_model(
                 <= CAP_RW_N_TO_ALL
             )
 
-        # Constraint: KLP combined capacity to J + SP + V + B3.
+        # Nebenbedingung: KLP-Gesamtkapazität nach Jänschwalde, SP, Veredlung und B3.
         model.cap_rw_n_to_j_sp_v_b3 = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_to_j_sp_v_b3_rule)
 
         def cap_w_to_sp_v_rule(m, w, d, s):
             return m.x["Welzow", "SP", w, d, s] + m.x["Welzow", "V", w, d, s] <= CAP_W_TO_SP_V
 
-        # Constraint: KUP combined capacity to SP + V.
+        # Nebenbedingung: KUP-Gesamtkapazität nach Schwarze Pumpe und Veredlung.
         model.cap_w_to_sp_v = pyo.Constraint(model.W, model.D, model.S, rule=cap_w_to_sp_v_rule)
 
         def cap_w_to_sp_v_b3_rule(m, w, d, s):
@@ -1295,7 +1429,7 @@ def build_model(
                 <= CAP_W_TO_SP_V_B3
             )
 
-        # Constraint: KUP combined capacity to SP + V + B3.
+        # Nebenbedingung: KUP-Gesamtkapazität nach Schwarze Pumpe, Veredlung und B3.
         model.cap_w_to_sp_v_b3 = pyo.Constraint(model.W, model.D, model.S, rule=cap_w_to_sp_v_b3_rule)
 
         def cap_rw_n_w_to_j_rule(m, w, d, s):
@@ -1306,11 +1440,10 @@ def build_model(
                 <= CAP_RW_N_W_TO_J
             )
 
-        # Constraint: Combined KLP/KUP capacity to J.
+        # Nebenbedingung: Kombinierte KLP/KUP-Kapazität nach Jänschwalde.
         model.cap_rw_n_w_to_j = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_w_to_j_rule)
 
         kvb_nord = tables.get("zugdurchlass_kvb_nord")
-        print(kvb_nord)
         if kvb_nord is not None:
             day_map = {"Sat": "Sa", "Sun": "So", "Mon": "Mo", "Tue": "Di", "Wed": "Mi", "Thu": "Do", "Fri": "Fr"}
             schicht_order = ["F", "S", "N"]
@@ -1346,7 +1479,7 @@ def build_model(
                     <= v
                 )
 
-            # Constraint: KVB Nord short-term capacity to J (per day/shift if provided).
+            # Nebenbedingung: KVB-Nord-Kurzfristkapazität nach Jänschwalde je Schicht (sofern vorhanden).
             model.cap_kvb_nord = pyo.Constraint(model.W, model.D, model.S, rule=cap_kvb_nord_rule)
 
         def cap_rw_n_w_to_bw3_rule(m, w, d, s):
@@ -1357,7 +1490,7 @@ def build_model(
                 <= CAP_RW_N_W_TO_B3
             )
 
-        # Constraint: Combined KLP/KUP capacity to B3.
+        # Nebenbedingung: Kombinierte KLP/KUP-Kapazität nach Boxberg Werk 3.
         model.cap_rw_n_w_to_bw3 = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_w_to_bw3_rule)
 
         def cap_rw_n_j_and_w_b3_rule(m, w, d, s):
@@ -1368,7 +1501,7 @@ def build_model(
                 <= CAP_RW_N_J_AND_W_B3
             )
 
-        # Constraint: Combined capacity (KLP->J, KUP->B3).
+        # Nebenbedingung: Kombinierte Kapazität KLP→Jänschwalde und KUP→Boxberg Werk 3.
         model.cap_rw_n_j_and_w_b3 = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_j_and_w_b3_rule)
 
         def cap_rw_n_to_j_sp_v_and_w_to_b3_rule(m, w, d, s):
@@ -1383,7 +1516,7 @@ def build_model(
                 <= CAP_RW_N_JSPV_AND_W_B3
             )
 
-        # Constraint: Combined capacity (KLP->J/SP/V, KUP->B3).
+        # Nebenbedingung: Kombinierte Kapazität KLP→J/SP/Veredlung und KUP→Boxberg Werk 3.
         model.cap_rw_n_to_j_sp_v_and_w_to_b3 = pyo.Constraint(
             model.W, model.D, model.S, rule=cap_rw_n_to_j_sp_v_and_w_to_b3_rule
         )
@@ -1433,10 +1566,6 @@ def solve_model(
         opt.options["TimeLimit"] = time_limit
         opt.options["MIPGap"] = mip_gap
         opt.options["LogFile"] = "gurobi_first_opt_model.log"
-    elif solver_name == "scip":
-        # SCIP uses hierarchical parameter names.
-        opt.options["limits/time"] = time_limit
-        opt.options["limits/gap"] = mip_gap
     solve_kwargs = {"tee": True, "symbolic_solver_labels": True}
     if use_warmstart and hasattr(opt, "warm_start_capable"):
         try:

webapp/backend/main.py

@@ -39,7 +39,7 @@ def _get_solver_availability() -> dict[str, bool]:
         return {"highs": False, "gurobi": False}
 
     availability: dict[str, bool] = {}
-    for solver_name in ("highs", "gurobi", "scip"):
+    for solver_name in ("highs", "gurobi"):
         try:
             solver = pyo.SolverFactory(solver_name)
             availability[solver_name] = bool(solver) and bool(

webapp/frontend/src/App.jsx

@@ -81,44 +81,7 @@ function parseSolverProgress(logText, solverName) {
     }
   };
 
-  if (solver === "scip") {
+  if (solver === "gurobi") {
-    for (const line of reversed) {
-      if (!metrics.gapPct) {
-        const m = line.match(/gap[^0-9-]*([0-9]+(?:\.[0-9]+)?)\s*%/i);
-        if (m) assignIfNumber("gapPct", m[1]);
-      }
-      if (!metrics.nodes) {
-        const m = line.match(/Solving Nodes\s*:\s*([0-9]+)/i);
-        if (m) assignIfNumber("nodes", m[1]);
-      }
-      if (!metrics.bestPrimal) {
-        const m = line.match(/Primal Bound\s*:\s*([+-]?[0-9.eE]+)/i);
-        if (m) assignIfNumber("bestPrimal", m[1]);
-      }
-      if (!metrics.bestDual) {
-        const m = line.match(/Dual Bound\s*:\s*([+-]?[0-9.eE]+)/i);
-        if (m) assignIfNumber("bestDual", m[1]);
-      }
-
-      // Best-effort parsing of SCIP progress table rows: time|node|...|gap
-      if ((!metrics.gapPct || !metrics.nodes) && line.includes("|")) {
-        if (!metrics.nodes) {
-          const nodeCol = line.match(/^\s*[^|]+\|\s*([0-9]+)\s*\|/);
-          if (nodeCol) assignIfNumber("nodes", nodeCol[1]);
-        }
-        if (!metrics.gapPct) {
-          const percents = [...line.matchAll(/([0-9]+(?:\.[0-9]+)?)\s*%/g)];
-          if (percents.length > 0) {
-            assignIfNumber("gapPct", percents[percents.length - 1][1]);
-          }
-        }
-      }
-
-      if (metrics.gapPct != null && metrics.nodes != null && metrics.bestPrimal != null) {
-        break;
-      }
-    }
-  } else if (solver === "gurobi") {
     for (const line of reversed) {
       if (!metrics.gapPct) {
         const m = line.match(/\bgap\s+([0-9]+(?:\.[0-9]+)?)%/i);
@@ -141,7 +104,6 @@ function parseSolverProgress(logText, solverName) {
   }
 
   const statusPatterns = [
-    /SCIP Status\s*:\s*(.+)/i,
     /Model status\s*:\s*(.+)/i,
     /Status\s*:\s*(.+)/i,
   ];