diff --git a/src/optimization/__pycache__/model_builder.cpython-313.pyc b/src/optimization/__pycache__/model_builder.cpython-313.pyc index e3cdeb7..e0dc93a 100644 Binary files a/src/optimization/__pycache__/model_builder.cpython-313.pyc and b/src/optimization/__pycache__/model_builder.cpython-313.pyc differ diff --git a/src/optimization/model_builder.py b/src/optimization/model_builder.py index b6b560c..7091272 100644 --- a/src/optimization/model_builder.py +++ b/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) - model.shift_excess = pyo.Var(model.I, model.J, model.W, model.D, model.SHIFT_PAIRS, domain=pyo.NonNegativeReals) + # Nachfolgewoche je Bergbauwoche (für wochenübergreifende Fenster) + 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 d in model.D: - for s1, s2 in SHIFT_PAIRS: - if (i, j, w, d, s1) in model.x and (i, j, w, d, s2) in model.x: - model.shift_dev_abs.add( - model.shift_dev_soft[i, j, w, d, s1, s2] - >= model.x[i, j, w, d, s1] - model.x[i, j, w, d, s2] - ) - model.shift_dev_abs.add( - model.shift_dev_soft[i, j, w, d, s1, s2] - >= model.x[i, j, w, d, s2] - model.x[i, j, w, d, s1] - ) - 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 - ) + for p_idx, (s, w_a, d_a, w_b, d_b) in enumerate(smooth_pairs): + model.shift_smooth_abs.add( + model.shift_smooth_dev[i, j, p_idx] + >= model.x[i, j, w_a, d_a, s] - model.x[i, j, w_b, d_b, s] + ) + model.shift_smooth_abs.add( + 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_smooth_excess_def.add( + model.shift_smooth_excess[i, j, p_idx] + >= model.shift_smooth_dev[i, j, p_idx] - SHIFT_SMOOTH_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( - model.shift_excess[i, j, w, d, s1, s2] + shift_smooth_penalty = lambda_shift_smooth * sum( + model.shift_smooth_excess[i, j, p_idx] for i in model.I for j in model.J - for w in model.W - 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 + for p_idx in model.SMOOTH_PAIR_IDX ) mix_target_penalty = lambda_mix * sum( @@ -985,27 +1096,32 @@ def build_model( return ( deviation_penalty - + smoothness_penalty - + shift_balance_penalty + + shift_smooth_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: diff --git a/webapp/backend/main.py b/webapp/backend/main.py index f05ba32..9c442aa 100644 --- a/webapp/backend/main.py +++ b/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( diff --git a/webapp/frontend/src/App.jsx b/webapp/frontend/src/App.jsx index 405c87e..a388110 100644 --- a/webapp/frontend/src/App.jsx +++ b/webapp/frontend/src/App.jsx @@ -81,44 +81,7 @@ function parseSolverProgress(logText, solverName) { } }; - if (solver === "scip") { - 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") { + 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, ];