commit f984bc4cf84ee7120c85f0672001d8ee7c6a3cba Author: Nicolai Date: Tue Mar 17 11:52:27 2026 +0100 Initial clean commit diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..f3b872a --- /dev/null +++ b/.gitignore @@ -0,0 +1,14 @@ +data/ +*.parquet +*.xlsx +*.pdf +.venv/ +node_modules/ +.DS_Store +latex/*.pdf +*.log +data/ +*.parquet +*.xlsx +latex/ +latex/*.pdf diff --git a/README.md b/README.md new file mode 100644 index 0000000..089a491 --- /dev/null +++ b/README.md @@ -0,0 +1,2 @@ +# LEAG-COALLOG + diff --git a/src/optimization/__init__.py b/src/optimization/__init__.py new file mode 100644 index 0000000..2c68ca6 --- /dev/null +++ b/src/optimization/__init__.py @@ -0,0 +1,3 @@ +from .model_builder import build_model, load_tables, solve_model + +__all__ = ["build_model", "load_tables", "solve_model"] diff --git a/src/optimization/__pycache__/__init__.cpython-313.pyc b/src/optimization/__pycache__/__init__.cpython-313.pyc new file mode 100644 index 0000000..d65bd4b Binary files /dev/null and b/src/optimization/__pycache__/__init__.cpython-313.pyc differ diff --git a/src/optimization/__pycache__/model_builder.cpython-313.pyc b/src/optimization/__pycache__/model_builder.cpython-313.pyc new file mode 100644 index 0000000..4a9a064 Binary files /dev/null and b/src/optimization/__pycache__/model_builder.cpython-313.pyc differ diff --git a/src/optimization/__pycache__/run_optimization.cpython-313.pyc b/src/optimization/__pycache__/run_optimization.cpython-313.pyc new file mode 100644 index 0000000..0b42911 Binary files /dev/null and b/src/optimization/__pycache__/run_optimization.cpython-313.pyc differ diff --git a/src/optimization/model_builder.py b/src/optimization/model_builder.py new file mode 100644 index 0000000..4d3434b --- /dev/null +++ b/src/optimization/model_builder.py @@ -0,0 +1,1479 @@ +from __future__ import annotations + +import math +from pathlib import Path + +import numpy as np +import pandas as pd +import pyomo.environ as pyo +from pyomo.environ import value +from pyomo.opt import TerminationCondition + +DEFAULT_STEP_SIZE_TONNES = 1000 + +# ----------------------------------------------------------------------------- +# OVERVIEW +# ----------------------------------------------------------------------------- +# This module contains the full Pyomo optimization model that was previously +# implemented in the notebook. The workflow is: +# 1) load_tables(...) loads all parquet inputs from a directory. +# 2) build_model(...) constructs a ConcreteModel with: +# - Sets: sources (I), consumers (J), weeks (W), days (D), shifts (S), months (M) +# - Parameters: demands, tolerances, capacities, mix targets +# - Decision variable: k (integer steps of 1,000 t), with x as derived flow +# - Constraints (in order of creation): +# 1) Daily delivery tolerances (power plants) +# 2) Weekly/monthly plant tolerances + total system tolerances +# 3) Veredlung day/week/month tolerances (Nochten/Welzow/Total) +# 4) Mix bounds (hard min/max) and target band (soft) +# 5) Deviation split & max-deviation pattern controls +# 6) Veredlung deviation splits (Nochten/Welzow) +# 7) Smoothness constraints across shifts +# 8) Route exclusions and capacity limits (monthly, loading, availability) +# 9) Transport caps from zugdurchlass (if table provided) +# - Objective: penalizes deviations & smoothness, applies route incentives +# 3) solve_model(...) runs the solver with configured limits. +# ----------------------------------------------------------------------------- + + +def load_tables(data_dir: Path) -> dict[str, pd.DataFrame]: + tables: dict[str, pd.DataFrame] = {} + for parquet_file in sorted(data_dir.glob("*.parquet")): + tables[parquet_file.stem] = pd.read_parquet(parquet_file) + return tables + + +def _match_kraftwerk_row(df: pd.DataFrame, j_key: str) -> pd.Series: + if j_key == "J": + mask = df["kraftwerk"].str.contains("JW") + elif j_key == "SP": + mask = df["kraftwerk"].str.contains("SP") + elif j_key == "B3": + mask = df["kraftwerk"].str.contains("BW3") + elif j_key == "B4": + mask = df["kraftwerk"].str.contains("BW4") + else: + mask = np.zeros(len(df), dtype=bool) + + row = df.loc[mask] + if row.empty: + raise ValueError(f"Keine Bounds Zeile für Kraftwerk {j_key} gefunden") + return row.iloc[0] + + +def _match_ver_row(df: pd.DataFrame, label_exact: str) -> pd.Series: + row = df[df["kohlesorte"] == label_exact] + if row.empty: + raise ValueError(f"Keine Veredlung-Bounds für {label_exact}") + return row.iloc[0] + + +def _match_bounds_row(df: pd.DataFrame, j_key: str, zeitraum: str) -> pd.Series: + label_map = { + "J": "Jänschwalde", + "SP": "Schwarze Pumpe", + "B3": "Boxberg Werk 3", + "B4": "Boxberg Werk 4", + "ALL": "Gesamt", + } + label = label_map[j_key] + row = df[(df["zeitraum"] == zeitraum) & df["kraftwerk"].str.contains(label)] + if row.empty: + raise ValueError(f"Keine Bounds-Zeile für {j_key} ({zeitraum})") + return row.iloc[0] + + +def _tol_from_row( + row: pd.Series, + d_nom: float, + step: int = DEFAULT_STEP_SIZE_TONNES, + min_steps: int = 1, +) -> tuple[float, float]: + lowers: list[float] = [] + uppers: list[float] = [] + if not pd.isna(row.get("minus")): + lowers.append(float(row["minus"])) + if not pd.isna(row.get("plus")): + uppers.append(float(row["plus"])) + if not pd.isna(row.get("minus_pct")): + lowers.append(float(row["minus_pct"]) * d_nom) + if not pd.isna(row.get("plus_pct")): + uppers.append(float(row["plus_pct"]) * d_nom) + if not lowers or not uppers: + raise ValueError("Toleranz nicht definiert") + + a_min = max(lowers) + a_max = min(uppers) + a_min = min(a_min, -min_steps * step) + a_max = max(a_max, min_steps * step) + return a_min, a_max + + +def _mining_week_from_date(datum: pd.Timestamp) -> int: + shifted = datum + pd.Timedelta(days=2) + iso = shifted.isocalendar() + return int(iso.year) * 100 + int(iso.week) + + +def build_model( + tables: dict[str, pd.DataFrame], step_size_tonnes: int = DEFAULT_STEP_SIZE_TONNES +) -> pyo.ConcreteModel: + if step_size_tonnes <= 0: + raise ValueError("step_size_tonnes must be positive") + bedarf = tables["rohkohlebedarf"][ + [ + "datum", + "JW", + "SP", + "BW3", + "BW4", + "Veredel_Nochtener", + "Veredel_Welzower", + ] + ].copy() + + bedarf["datum"] = pd.to_datetime(bedarf["datum"]) + bedarf["weekday"] = bedarf["datum"].dt.weekday + shifted = bedarf["datum"] + pd.Timedelta(days=2) + iso = shifted.dt.isocalendar() + bedarf["week"] = iso.year.astype(int) * 100 + iso.week.astype(int) + + weekday_map = {5: "Sa", 6: "So", 0: "Mo", 1: "Di", 2: "Mi", 3: "Do", 4: "Fr"} + bedarf["day"] = bedarf["weekday"].map(weekday_map) + + wd_to_date = ( + bedarf.sort_values("datum") + .drop_duplicates(subset=["week", "day"]) + .set_index(["week", "day"])["datum"] + .to_dict() + ) + ih_dates_welzow = set() + ih_dates_boxberg = set() + + bounds_power_plants = tables["bounds_power_plants"] + bounds_day = bounds_power_plants[bounds_power_plants["zeitraum"] == "pro Tag"].copy() + + model = pyo.ConcreteModel() + + J_POWER = ["J", "SP", "B3", "B4"] + model.J = pyo.Set(initialize=J_POWER + ["V"]) + + weeks = sorted(bedarf["week"].unique().tolist()) + model.W = pyo.Set(initialize=weeks) + + model.D = pyo.Set(initialize=["Sa", "So", "Mo", "Di", "Mi", "Do", "Fr"]) + model.S = pyo.Set(initialize=["F", "S", "N"]) + model.I = pyo.Set(initialize=["Reichwalde", "Nochten", "Welzow"]) + + I_W = ["Welzow"] + I_N = ["Nochten"] + + model.d = pyo.Param(model.J, model.W, model.D, initialize=0.0, mutable=True, within=pyo.NonNegativeReals) + model.dV_N = pyo.Param(model.W, model.D, initialize=0.0, mutable=True) + model.dV_W = pyo.Param(model.W, model.D, initialize=0.0, mutable=True) + model.a_min_day = pyo.Param(model.J, model.W, model.D, initialize=0.0, mutable=True) + model.a_max_day = pyo.Param(model.J, model.W, model.D, initialize=0.0, mutable=True) + + column_map = {"JW": "J", "SP": "SP", "BW3": "B3", "BW4": "B4"} + bounds_by_j = {j: _match_kraftwerk_row(bounds_day, j) for j in model.J if j != "V"} + + ver_bounds = tables["veredelung_bounds"] + v_day = ver_bounds[ver_bounds["zeitraum"] == "pro Tag"] + v_week = ver_bounds[ver_bounds["zeitraum"] == "pro Woche"] + v_month = ver_bounds[ver_bounds["zeitraum"] == "pro Monat"] + + ver_row_day_N = _match_ver_row(v_day, "Nochtener-Kohle") + ver_row_day_W = _match_ver_row(v_day, "Welzower-Kohle") + ver_row_day_ALL = _match_ver_row(v_day, "Gesamt (NO+WZ)") + + ver_row_week_N = _match_ver_row(v_week, "Nochtener-Kohle") + ver_row_week_W = _match_ver_row(v_week, "Welzower-Kohle") + ver_row_week_ALL = _match_ver_row(v_week, "Gesamt (NO+WZ)") + + ver_row_month_N = _match_ver_row(v_month, "Nochtener-Kohle") + ver_row_month_W = _match_ver_row(v_month, "Welzower-Kohle") + ver_row_month_ALL = _match_ver_row(v_month, "Gesamt (NO+WZ)") + + for _, row in bedarf.iterrows(): + w = int(row["week"]) + d = row["day"] + + if (w in model.W) and (d in model.D): + vN = row.get("Veredel_Nochtener", np.nan) + if not pd.isna(vN): + model.dV_N[w, d] = float(vN) + vW = row.get("Veredel_Welzower", np.nan) + if not pd.isna(vW): + model.dV_W[w, d] = float(vW) + + for col, j in column_map.items(): + if not ((j in model.J) and (w in model.W) and (d in model.D)): + continue + + val = row.get(col, np.nan) + if pd.isna(val): + continue + + d_nom = float(val) + model.d[j, w, d] = d_nom + + b = bounds_by_j[j] + lower_candidates: list[float] = [] + upper_candidates: list[float] = [] + + minus_abs = b.get("minus") + plus_abs = b.get("plus") + minus_pct = b.get("minus_pct") + plus_pct = b.get("plus_pct") + + if not pd.isna(minus_abs): + lower_candidates.append(float(minus_abs)) + if not pd.isna(plus_abs): + upper_candidates.append(float(plus_abs)) + if not pd.isna(minus_pct): + lower_candidates.append(float(minus_pct) * d_nom) + if not pd.isna(plus_pct): + upper_candidates.append(float(plus_pct) * d_nom) + + if not lower_candidates or not upper_candidates: + raise ValueError(f"Keine gültigen Toleranzen für {j}, Woche {w}, Tag {d}") + + a_min = max(lower_candidates) + a_max = min(upper_candidates) + + model.a_min_day[j, w, d] = a_min + model.a_max_day[j, w, d] = a_max + + 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) + + def x_rule(model, i, j, w, d, s): + return model.step_size_tonnes * model.k[i, j, w, d, s] + + model.x = pyo.Expression(model.I, model.J, model.W, model.D, model.S, rule=x_rule) + + bunker_df = tables.get("bunker") + if bunker_df is not None: + j_bunker_map = { + "Jänschwalde": "J", + "Schwarze Pumpe": "SP", + "Boxberg Werk 3": "B3", + "ISP": "V", + } + bunker_rows = [] + for _, row in bunker_df.iterrows(): + j = j_bunker_map.get(row["anlage"]) + if j is None: + continue + bunker_rows.append((j, row)) + + J_BUNKER = sorted({j for j, _ in bunker_rows}) + model.J_BUNKER = pyo.Set(initialize=J_BUNKER) + + bunker_init = {j: 0.0 for j in J_BUNKER} + bunker_target = {j: 0.0 for j in J_BUNKER} + bunker_max = {j: 0.0 for j in J_BUNKER} + for j, row in bunker_rows: + bunker_init[j] = float(row["anfang_mo_di_t"]) # placeholder, refined below + bunker_target[j] = float(row["zielbestand_t"]) + bunker_max[j] = float(row["maximal_t"]) + + # Determine initial date for the month to pick correct anfang_* column. + first_date = min(wd_to_date.values()) + first_weekday = first_date.strftime("%a") + use_mo_di = first_weekday in {"Mon", "Tue"} + for j, row in bunker_rows: + bunker_init[j] = float(row["anfang_mo_di_t" if use_mo_di else "anfang_rest_t"]) + + model.bunker_init = pyo.Param(model.J_BUNKER, initialize=bunker_init, mutable=True, within=pyo.NonNegativeReals) + model.bunker_target = pyo.Param( + model.J_BUNKER, initialize=bunker_target, mutable=True, within=pyo.NonNegativeReals + ) + model.bunker_max = pyo.Param(model.J_BUNKER, initialize=bunker_max, mutable=True, within=pyo.NonNegativeReals) + + model.bunker = pyo.Var(model.I, model.J_BUNKER, model.W, model.D, domain=pyo.NonNegativeReals) + model.bunker_out = pyo.Var(model.I, model.J_BUNKER, model.W, model.D, model.S, domain=pyo.NonNegativeReals) + + def bunker_total_rule(m, j, w, d): + return sum(m.bunker[i, j, w, d] for i in m.I) + + model.bunker_total = pyo.Expression(model.J_BUNKER, model.W, model.D, rule=bunker_total_rule) + + def out_rule(m, i, j, w, d, s): + if j in m.J_BUNKER: + return m.bunker_out[i, j, w, d, s] + return m.x[i, j, w, d, s] + + model.out = pyo.Expression(model.I, model.J, model.W, model.D, model.S, rule=out_rule) + + def y_shift_rule(model, j, w, d, s): + return sum(model.out[i, j, w, d, s] for i in model.I) + + model.y_shift = pyo.Expression(model.J, model.W, model.D, model.S, rule=y_shift_rule) + + def y_rule(model, j, w, d): + return sum(model.y_shift[j, w, d, s] for s in model.S) + + model.y = pyo.Expression(model.J, model.W, model.D, rule=y_rule) + + time_points = [] + for w in model.W: + for d in model.D: + date = wd_to_date.get((w, d)) + if date is not None: + time_points.append((date, w, d)) + time_points.sort() + time_points_day = [(w, d) for _, w, d in time_points] + prev_map = {} + for idx, (w, d) in enumerate(time_points_day): + if idx == 0: + continue + prev_map[(w, d)] = time_points_day[idx - 1] + + model.bunker_balance = pyo.ConstraintList() + model.bunker_init_con = pyo.ConstraintList() + + first_point = time_points_day[0] if time_points_day else None + if first_point is not None: + w0, d0 = first_point + for j in model.J_BUNKER: + model.bunker_init_con.add( + sum(model.bunker[i, j, w0, d0] for i in model.I) + == model.bunker_init[j] + + sum(model.x[i, j, w0, d0, s] for i in model.I for s in model.S) + - sum(model.bunker_out[i, j, w0, d0, s] for i in model.I for s in model.S) + ) + + for j in model.J_BUNKER: + for w, d in time_points_day: + prev = prev_map.get((w, d)) + if prev is None: + continue + wp, dp = prev + for i in model.I: + model.bunker_balance.add( + model.bunker[i, j, w, d] + == model.bunker[i, j, wp, dp] + + 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) + ) + else: + def out_rule(m, i, j, w, d, s): + return m.x[i, j, w, d, s] + + model.out = pyo.Expression(model.I, model.J, model.W, model.D, model.S, rule=out_rule) + + def y_shift_rule(model, j, w, d, s): + return sum(model.x[i, j, w, d, s] for i in model.I) + + model.y_shift = pyo.Expression(model.J, model.W, model.D, model.S, rule=y_shift_rule) + + def y_rule(model, j, w, d): + return sum(model.x[i, j, w, d, s] for i in model.I for s in model.S) + + model.y = pyo.Expression(model.J, model.W, model.D, rule=y_rule) + + def y_delivery_shift_rule(model, j, w, d, s): + return sum(model.x[i, j, w, d, s] for i in model.I) + + model.y_delivery_shift = pyo.Expression(model.J, model.W, model.D, model.S, rule=y_delivery_shift_rule) + + def y_delivery_rule(model, j, w, d): + return sum(model.y_delivery_shift[j, w, d, s] for s in model.S) + + model.y_delivery = pyo.Expression(model.J, model.W, model.D, rule=y_delivery_rule) + + J_delivery = [j for j in model.J if j != "V"] + + def delivery_tolerance_rule(model, j, w, d): + return ( + model.d[j, w, d] + model.a_min_day[j, w, d], + model.y_delivery[j, w, d], + model.d[j, w, d] + model.a_max_day[j, w, d], + ) + + # Constraint: Daily delivery must stay within demand tolerance for power plants. + 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"] + week_tolerance_debug = 0 # extra weekly slack for debugging infeasibility + + model.d_week = pyo.Param(model.J, model.W, initialize=0.0, mutable=True) + model.d_month = pyo.Param(model.J, model.M, initialize=0.0, mutable=True) + + model.a_min_week = pyo.Param(model.J, model.W, initialize=0.0, mutable=True) + model.a_max_week = pyo.Param(model.J, model.W, initialize=0.0, mutable=True) + model.a_min_month = pyo.Param(model.J, model.M, initialize=0.0, mutable=True) + model.a_max_month = pyo.Param(model.J, model.M, initialize=0.0, mutable=True) + + model.d_week_total = pyo.Param(model.W, initialize=0.0, mutable=True) + model.d_month_total = pyo.Param(model.M, initialize=0.0, mutable=True) + model.a_min_week_total = pyo.Param(model.W, initialize=0.0, mutable=True) + model.a_max_week_total = pyo.Param(model.W, initialize=0.0, mutable=True) + model.a_min_month_total = pyo.Param(model.M, initialize=0.0, mutable=True) + model.a_max_month_total = pyo.Param(model.M, initialize=0.0, mutable=True) + + def y_week_rule(model, j, w): + return sum(model.y_delivery[j, w, d] for d in model.D) + + model.y_week = pyo.Expression(model.J, model.W, rule=y_week_rule) + + def y_month_rule(model, j, m): + weeks_in_m = [w for w, mm in week_to_month.items() if mm == m and w in model.W] + return sum(model.y_week[j, w] for w in weeks_in_m) + + model.y_month = pyo.Expression(model.J, model.M, rule=y_month_rule) + + model.y_week_total = pyo.Expression(model.W, rule=lambda m, w: sum(m.y_week[j, w] for j in J_power)) + model.y_month_total = pyo.Expression(model.M, rule=lambda m, mm: sum(m.y_month[j, mm] for j in J_power)) + + bounds_week = bounds_power_plants[bounds_power_plants["zeitraum"] == "pro Woche"] + bounds_month = bounds_power_plants[bounds_power_plants["zeitraum"] == "pro Monat"] + + for j in J_power: + for w in model.W: + d_w = sum(model.d[j, w, d] for d in model.D) + model.d_week[j, w] = d_w + row = _match_bounds_row(bounds_week, j, "pro Woche") + a_min, a_max = _tol_from_row(row, d_w, step=step_size_tonnes) + model.a_min_week[j, w] = a_min - week_tolerance_debug + model.a_max_week[j, w] = a_max + week_tolerance_debug + + for j in J_power: + for m in model.M: + weeks_in_m = [w for w, mm in week_to_month.items() if mm == m and w in model.W] + d_m = sum(model.d_week[j, w] for w in weeks_in_m) + model.d_month[j, m] = d_m + row = _match_bounds_row(bounds_month, j, "pro Monat") + a_min, a_max = _tol_from_row(row, d_m, step=step_size_tonnes) + model.a_min_month[j, m] = a_min + model.a_max_month[j, m] = a_max + + row_week_total = _match_bounds_row(bounds_week, "ALL", "pro Woche") + row_month_total = _match_bounds_row(bounds_month, "ALL", "pro Monat") + for w in model.W: + d_tot = sum(model.d_week[j, w] for j in J_power) + model.d_week_total[w] = d_tot + a_min, a_max = _tol_from_row(row_week_total, d_tot, step=step_size_tonnes) + model.a_min_week_total[w] = a_min - week_tolerance_debug + model.a_max_week_total[w] = a_max + week_tolerance_debug + + for m in model.M: + weeks_in_m = [w for w, mm in week_to_month.items() if mm == m and w in model.W] + d_tot = sum(model.d_month[j, m] for j in J_power) + model.d_month_total[m] = d_tot + a_min, a_max = _tol_from_row(row_month_total, d_tot, step=step_size_tonnes) + model.a_min_month_total[m] = a_min + model.a_max_month_total[m] = a_max + + def week_tolerance_rule(model, j, w): + return ( + model.d_week[j, w] + model.a_min_week[j, w], + model.y_week[j, w], + model.d_week[j, w] + model.a_max_week[j, w], + ) + + # Constraint: Weekly deliveries per plant stay within weekly tolerance band. + model.week_tolerance = pyo.Constraint(J_power, model.W, rule=week_tolerance_rule) + + def month_tolerance_rule(model, j, m): + return ( + model.d_month[j, m] + model.a_min_month[j, m], + model.y_month[j, m], + model.d_month[j, m] + model.a_max_month[j, m], + ) + + # Constraint: Monthly deliveries per plant stay within monthly tolerance band. + model.month_tolerance = pyo.Constraint(J_power, model.M, rule=month_tolerance_rule) + + def week_total_rule(model, w): + return ( + model.d_week_total[w] + model.a_min_week_total[w], + model.y_week_total[w], + model.d_week_total[w] + model.a_max_week_total[w], + ) + + # Constraint: Total weekly deliveries across all plants stay within tolerance band. + model.week_total_tolerance = pyo.Constraint(model.W, rule=week_total_rule) + + def month_total_rule(model, m): + return ( + model.d_month_total[m] + model.a_min_month_total[m], + model.y_month_total[m], + model.d_month_total[m] + model.a_max_month_total[m], + ) + + # Constraint: Total monthly deliveries across all plants stay within tolerance band. + model.month_total_tolerance = pyo.Constraint(model.M, rule=month_total_rule) + + 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) + model.a_max_V_week = pyo.Param(["N", "W", "ALL"], model.W, initialize=0.0, mutable=True) + model.a_min_V_month = pyo.Param(["N", "W", "ALL"], model.M, initialize=0.0, mutable=True) + model.a_max_V_month = pyo.Param(["N", "W", "ALL"], model.M, initialize=0.0, mutable=True) + + for w in model.W: + for d in model.D: + dN = model.dV_N[w, d] + dW = model.dV_W[w, d] + dA = dN + dW + model.a_min_V_day["N", w, d], model.a_max_V_day["N", w, d] = _tol_from_row( + ver_row_day_N, dN, step=step_size_tonnes + ) + model.a_min_V_day["W", w, d], model.a_max_V_day["W", w, d] = _tol_from_row( + ver_row_day_W, dW, step=step_size_tonnes + ) + model.a_min_V_day["ALL", w, d], model.a_max_V_day["ALL", w, d] = _tol_from_row( + ver_row_day_ALL, dA, step=step_size_tonnes + ) + + for w in model.W: + dN = sum(model.dV_N[w, dd] for dd in model.D) + dW = sum(model.dV_W[w, dd] for dd in model.D) + dA = dN + dW + a_min_n, a_max_n = _tol_from_row(ver_row_week_N, dN, step=step_size_tonnes) + a_min_w, a_max_w = _tol_from_row(ver_row_week_W, dW, step=step_size_tonnes) + a_min_a, a_max_a = _tol_from_row(ver_row_week_ALL, dA, step=step_size_tonnes) + model.a_min_V_week["N", w] = a_min_n - week_tolerance_debug + model.a_max_V_week["N", w] = a_max_n + week_tolerance_debug + model.a_min_V_week["W", w] = a_min_w - week_tolerance_debug + model.a_max_V_week["W", w] = a_max_w + week_tolerance_debug + model.a_min_V_week["ALL", w] = a_min_a - week_tolerance_debug + model.a_max_V_week["ALL", w] = a_max_a + week_tolerance_debug + + for m in model.M: + weeks_in_m = [w for w, mm in week_to_month.items() if mm == m and w in model.W] + dN = sum(model.dV_N[w, dd] for w in weeks_in_m for dd in model.D) + dW = sum(model.dV_W[w, dd] for w in weeks_in_m for dd in model.D) + dA = dN + dW + model.a_min_V_month["N", m], model.a_max_V_month["N", m] = _tol_from_row( + ver_row_month_N, dN, step=step_size_tonnes + ) + model.a_min_V_month["W", m], model.a_max_V_month["W", m] = _tol_from_row( + ver_row_month_W, dW, step=step_size_tonnes + ) + model.a_min_V_month["ALL", m], model.a_max_V_month["ALL", m] = _tol_from_row( + ver_row_month_ALL, dA, step=step_size_tonnes + ) + + def v_w_day_rule(m, w, d): + return ( + m.dV_W[w, d] + m.a_min_V_day["W", w, d], + sum(m.x[i, "V", w, d, s] for i in I_W for s in m.S), + m.dV_W[w, d] + m.a_max_V_day["W", w, d], + ) + + # Constraint: Veredlung day tolerance for Welzow coal. + model.v_w_day_tol = pyo.Constraint(model.W, model.D, rule=v_w_day_rule) + + def v_n_day_rule(m, w, d): + return ( + m.dV_N[w, d] + m.a_min_V_day["N", w, d], + sum(m.x[i, "V", w, d, s] for i in I_N for s in m.S), + m.dV_N[w, d] + m.a_max_V_day["N", w, d], + ) + + # Constraint: Veredlung day tolerance for Nochten coal. + model.v_n_day_tol = pyo.Constraint(model.W, model.D, rule=v_n_day_rule) + + def v_all_day_rule(m, w, d): + return ( + m.dV_W[w, d] + m.dV_N[w, d] + m.a_min_V_day["ALL", w, d], + sum(m.x[i, "V", w, d, s] for i in I_W + I_N for s in m.S), + 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. + model.v_all_day_tol = pyo.Constraint(model.W, model.D, rule=v_all_day_rule) + + def v_w_week_rule(m, w): + return ( + sum(m.dV_W[w, d] for d in m.D) + m.a_min_V_week["W", w], + sum(m.x[i, "V", w, d, s] for i in I_W for d in m.D for s in m.S), + 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. + model.v_w_week_tol = pyo.Constraint(model.W, rule=v_w_week_rule) + + def v_n_week_rule(m, w): + return ( + sum(m.dV_N[w, d] for d in m.D) + m.a_min_V_week["N", w], + sum(m.x[i, "V", w, d, s] for i in I_N for d in m.D for s in m.S), + 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. + model.v_n_week_tol = pyo.Constraint(model.W, rule=v_n_week_rule) + + def v_all_week_rule(m, w): + return ( + sum(m.dV_W[w, d] + m.dV_N[w, d] for d in m.D) + m.a_min_V_week["ALL", w], + sum(m.x[i, "V", w, d, s] for i in I_W + I_N for d in m.D for s in m.S), + 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. + model.v_all_week_tol = pyo.Constraint(model.W, rule=v_all_week_rule) + + def v_w_month_rule(m, mm): + weeks_in_m = [w for w, month in week_to_month.items() if month == mm and w in m.W] + return ( + sum(m.dV_W[w, d] for w in weeks_in_m for d in m.D) + m.a_min_V_month["W", mm], + sum(m.x[i, "V", w, d, s] for i in I_W for w in weeks_in_m for d in m.D for s in m.S), + 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. + model.v_w_month_tol = pyo.Constraint(model.M, rule=v_w_month_rule) + + def v_n_month_rule(m, mm): + weeks_in_m = [w for w, month in week_to_month.items() if month == mm and w in m.W] + return ( + sum(m.dV_N[w, d] for w in weeks_in_m for d in m.D) + m.a_min_V_month["N", mm], + sum(m.x[i, "V", w, d, s] for i in I_N for w in weeks_in_m for d in m.D for s in m.S), + 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. + model.v_n_month_tol = pyo.Constraint(model.M, rule=v_n_month_rule) + + def v_all_month_rule(m, mm): + weeks_in_m = [w for w, month in week_to_month.items() if month == mm and w in m.W] + return ( + sum(m.dV_W[w, d] + m.dV_N[w, d] for w in weeks_in_m for d in m.D) + m.a_min_V_month["ALL", mm], + sum(m.x[i, "V", w, d, s] for i in I_W + I_N for w in weeks_in_m for d in m.D for s in m.S), + 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. + model.v_all_month_tol = pyo.Constraint(model.M, rule=v_all_month_rule) + + j_map = { + "Jänschwalde": "J", + "Schwarze Pumpe": "SP", + "Boxberg Werk 3": "B3", + "Boxberg Werk 4": "B4", + } + + i_map = { + "Reichwalder-Kohle": "Reichwalde", + "Nochtener-Kohle": "Nochten", + "Welzower-Kohle": "Welzow", + } + + mix_df = tables["kohle_mix"] + + 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) + + for _, row in mix_df.iterrows(): + j = j_map.get(row["kraftwerk"]) + i = i_map.get(row["kohlesorte"]) + if (i in model.I) and (j in model.J): + model.alpha_min[i, j] = float(row["minimal"]) + model.alpha_max[i, j] = float(row["maximal"]) + + J_MIX = ["J", "SP", "B3", "B4", "V"] + + def x_day_rule(m, i, j, w, d): + return sum(m.x[i, j, w, d, s] for s in m.S) + + model.x_day = pyo.Expression(model.I, model.J, model.W, model.D, rule=x_day_rule) + + def y_day_rule(m, j, w, d): + return sum(m.x_day[i, j, w, d] for i in m.I) + + model.y_day = pyo.Expression(model.J, model.W, model.D, rule=y_day_rule) + + def mix_lower_rule(m, i, j, w, d): + if bunker_df is not None and (j not in J_MIX or j not in m.J_BUNKER): + return pyo.Constraint.Skip + return m.alpha_min[i, j] * m.y_day[j, w, d] <= m.x_day[i, j, w, d] + + def mix_upper_rule(m, i, j, w, d): + if bunker_df is not None and (j not in J_MIX or j not in m.J_BUNKER): + 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. + 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) + + model.alpha_target_low = pyo.Param(model.I, model.J, initialize=0.0, mutable=True) + model.alpha_target_high = pyo.Param(model.I, model.J, initialize=1.0, mutable=True) + model.alpha_target_mid = pyo.Param(model.I, model.J, initialize=0.5, mutable=True) + + for _, row in mix_df.iterrows(): + j = j_map.get(row["kraftwerk"]) + i = i_map.get(row["kohlesorte"]) + if (i in model.I) and (j in model.J): + low = float(row["ziel_low"]) + high = float(row["ziel_high"]) + model.alpha_target_low[i, j] = low + model.alpha_target_high[i, j] = high + model.alpha_target_mid[i, j] = 0.5 * (low + high) + + model.mix_target_low_dev = pyo.Var(model.I, model.J, model.W, model.D, domain=pyo.NonNegativeReals) + model.mix_target_high_dev = pyo.Var(model.I, model.J, model.W, model.D, domain=pyo.NonNegativeReals) + model.mix_target_mid_dev = pyo.Var(model.I, model.J, model.W, model.D, domain=pyo.NonNegativeReals) + + def mix_target_low_rule(m, i, j, w, d): + if bunker_df is not None and (j not in J_MIX or j not in m.J_BUNKER): + return pyo.Constraint.Skip + return m.x_day[i, j, w, d] + m.mix_target_low_dev[i, j, w, d] >= m.alpha_target_low[i, j] * m.y_day[j, w, d] + + def mix_target_high_rule(m, i, j, w, d): + if bunker_df is not None and (j not in J_MIX or j not in m.J_BUNKER): + 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). + 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) + + def mix_target_mid_hi(m, i, j, w, d): + if bunker_df is not None and (j not in J_MIX or j not in m.J_BUNKER): + return pyo.Constraint.Skip + return m.x_day[i, j, w, d] - m.alpha_target_mid[i, j] * m.y_day[j, w, d] <= m.mix_target_mid_dev[i, j, w, d] + + def mix_target_mid_lo(m, i, j, w, d): + if bunker_df is not None and (j not in J_MIX or j not in m.J_BUNKER): + return pyo.Constraint.Skip + return m.alpha_target_mid[i, j] * m.y_day[j, w, d] - m.x_day[i, j, w, d] <= m.mix_target_mid_dev[i, j, w, d] + + model.mix_target_mid_hi = pyo.Constraint(model.I, model.J, model.W, model.D, rule=mix_target_mid_hi) + model.mix_target_mid_lo = pyo.Constraint(model.I, model.J, model.W, model.D, rule=mix_target_mid_lo) + + if bunker_df is not None: + vorfahr_df = tables.get("bunker_vorfahrfenster") + if vorfahr_df is not None and not vorfahr_df.empty: + vorfahrfenster_tage = int(vorfahr_df.iloc[0]["vorfahrfenster_tage"]) + else: + vorfahrfenster_tage = 0 + + date_to_wd = {pd.to_datetime(v).normalize(): k for k, v in wd_to_date.items()} + + model.bunker_target_short = pyo.Var(model.J_BUNKER, model.W, model.D, domain=pyo.NonNegativeReals) + + model.bunker_target_con = pyo.ConstraintList() + model.bunker_max_con = pyo.ConstraintList() + + for j in model.J_BUNKER: + for w in model.W: + for d in model.D: + model.bunker_max_con.add(model.bunker_total[j, w, d] <= model.bunker_max[j]) + + # B3 bunker mix targets (soft) to discourage Welzow and favor RW/NO in bunker. + model.b3_bunker_mix_low_dev = pyo.Var(model.I, model.W, model.D, domain=pyo.NonNegativeReals) + model.b3_bunker_mix_high_dev = pyo.Var(model.I, model.W, model.D, domain=pyo.NonNegativeReals) + + def b3_bunker_mix_low_rule(m, i, w, d): + if "B3" not in m.J_BUNKER: + return pyo.Constraint.Skip + return ( + m.bunker[i, "B3", w, d] + m.b3_bunker_mix_low_dev[i, w, d] + >= m.alpha_target_low[i, "B3"] * m.bunker_total["B3", w, d] + ) + + def b3_bunker_mix_high_rule(m, i, w, d): + if "B3" not in m.J_BUNKER: + return pyo.Constraint.Skip + return ( + m.bunker[i, "B3", w, d] + <= m.alpha_target_high[i, "B3"] * m.bunker_total["B3", w, d] + m.b3_bunker_mix_high_dev[i, w, d] + ) + + model.b3_bunker_mix_low = pyo.Constraint(model.I, model.W, model.D, rule=b3_bunker_mix_low_rule) + model.b3_bunker_mix_high = pyo.Constraint(model.I, model.W, model.D, rule=b3_bunker_mix_high_rule) + + # Zielbestand vor Montag (Sonntag) als soft constraint. + for date, (w, d) in date_to_wd.items(): + if date.strftime("%a") != "Mon": + continue + prev_date = date - pd.Timedelta(days=1) + if prev_date not in date_to_wd: + continue + w_prev, d_prev = date_to_wd[prev_date] + for j in model.J_BUNKER: + model.bunker_target_con.add( + model.bunker_target_short[j, w_prev, d_prev] + >= model.bunker_target[j] - model.bunker_total[j, w_prev, d_prev] + ) + + ih_map = {"welzow": ih_dates_welzow, "boxberg": ih_dates_boxberg} + + for src, ih_dates in ih_map.items(): + for ih_date in ih_dates: + target_date = ih_date - pd.Timedelta(days=vorfahrfenster_tage) + if target_date not in date_to_wd: + continue + w_prev, d_prev = date_to_wd[target_date] + for j in model.J_BUNKER: + if j not in J_MIX: + continue + model.bunker_target_con.add( + model.bunker_target_short[j, w_prev, d_prev] + >= model.bunker_target[j] - model.bunker_total[j, w_prev, d_prev] + ) + + model.lambda_J = pyo.Param(initialize=100_000, mutable=True, within=pyo.NonNegativeReals) + + def welzow_to_J(m): + return sum( + m.x["Welzow", "J", w, d, s] + for w in m.W + for d in m.D + for s in m.S + if ("Welzow", "J", w, d, s) in m.x + ) + + model.bonus_welzow_J = pyo.Expression(rule=welzow_to_J) + + 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] + + model.demand = pyo.Expression(model.J, model.W, model.D, rule=demand_rule) + model.dev = pyo.Expression( + model.J, model.W, model.D, rule=lambda m, j, w, d: m.y_delivery[j, w, d] - m.demand[j, w, d] + ) + + J_dev = list(J_power) + model.dev_pos = pyo.Var(J_dev, model.W, model.D, domain=pyo.NonNegativeReals) + model.dev_neg = pyo.Var(J_dev, model.W, model.D, domain=pyo.NonNegativeReals) + + 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. + 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) + + 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. + model.max_reached = pyo.Constraint(J_power, model.W, model.D, rule=max_reached_rule) + + D_list = list(model.D.ordered_data()) + + def no_three_in_a_row_rule(m, j, w, idx): + 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. + 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): + return sum(m.x["Nochten", "V", w, d, s] for s in m.S) + + def y_V_welzow(m, w, d): + return sum(m.x["Welzow", "V", w, d, s] for s in m.S) + + model.devV_N = pyo.Expression(model.W, model.D, rule=lambda m, w, d: y_V_nochten(m, w, d) - m.dV_N[w, d]) + model.devV_W = pyo.Expression(model.W, model.D, rule=lambda m, w, d: y_V_welzow(m, w, d) - m.dV_W[w, d]) + + model.devV_N_pos = pyo.Var(model.W, model.D, domain=pyo.NonNegativeReals) + model.devV_N_neg = pyo.Var(model.W, model.D, domain=pyo.NonNegativeReals) + 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. + 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. + 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] + ) + + model.lambda_glatt = pyo.Param(initialize=1e5, mutable=True, within=pyo.NonNegativeReals) + + model.m_ijwd = pyo.Expression( + model.I, + model.J, + model.W, + model.D, + rule=lambda m, i, j, w, d: (1 / len(m.S)) + * sum(m.x[i, j, w, d, s] for s in m.S if (i, j, w, d, s) in m.x), + ) + + model.t_glatt = pyo.Var(model.I, model.J, model.W, model.D, model.S, domain=pyo.NonNegativeReals) + + def glatt_hi(m, i, j, w, d, s): + return m.t_glatt[i, j, w, d, s] >= m.x[i, j, w, d, s] - m.m_ijwd[i, j, w, d] + + def glatt_lo(m, i, j, w, d, s): + return m.t_glatt[i, j, w, d, s] >= -(m.x[i, j, w, d, s] - m.m_ijwd[i, j, w, d]) + + # Constraint: Upper absolute deviation from average shift flow (smoothness). + model.glatt_hi = pyo.Constraint(model.I, model.J, model.W, model.D, model.S, rule=glatt_hi) + # Constraint: Lower absolute deviation from average shift flow (smoothness). + model.glatt_lo = pyo.Constraint(model.I, model.J, model.W, model.D, model.S, rule=glatt_lo) + + lambda_dev = 100_000 + lambda_v = 100_000 + lambda_shift = 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) + + # 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 + ) + + def objective_rule(model): + deviation_penalty = ( + lambda_dev + * sum(model.dev_pos[j, w, d] + model.dev_neg[j, w, d] for j in J_dev for w in model.W for d in model.D) + + lambda_v + * sum( + model.devV_N_pos[w, d] + model.devV_N_neg[w, d] + model.devV_W_pos[w, d] + model.devV_W_neg[w, d] + for w in model.W + for d in model.D + ) + ) + + smoothness_penalty = lambda_shift * sum( + model.shift_excess[i, j, w, d, s1, s2] + 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 + ) + + mix_target_penalty = lambda_mix * sum( + model.mix_target_low_dev[i, j, w, d] + model.mix_target_high_dev[i, j, w, d] + for i in model.I + for j in model.J + for w in model.W + for d in model.D + ) + + mix_target_mid_penalty = lambda_mix_mid * sum( + model.mix_target_mid_dev[i, j, w, d] for i in model.I for j in model.J for w in model.W for d in model.D + ) + + bunker_penalty = 0 + if bunker_df is not None: + bunker_penalty = lambda_bunker_target * sum( + model.bunker_target_short[j, w, d] for j in model.J_BUNKER for w in model.W for d in model.D + ) + b3_bunker_mix_penalty = lambda_b3_bunker_mix * sum( + model.b3_bunker_mix_low_dev[i, w, d] + model.b3_bunker_mix_high_dev[i, w, d] + for i in model.I + for w in model.W + for d in model.D + ) + else: + b3_bunker_mix_penalty = 0 + + return ( + deviation_penalty + + smoothness_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. + model.obj = pyo.Objective(rule=objective_rule, sense=pyo.minimize) + + def forbid_reichwalde_V_rule(m, w, d, s): + return m.x["Reichwalde", "V", w, d, s] == 0 + + # Constraint: Disallow Reichwalde coal to Veredlung. + 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. + model.forbid_welzow_B4 = pyo.Constraint(model.W, model.D, model.S, rule=forbid_welzow_B4_rule) + + cap_df = tables["foerderkapaz"] + cap_month = cap_df[cap_df["zeitraum"] == "pro Monat"] + + i_map = { + "Reichwalde (RW)": "Reichwalde", + "Nochten (NO)": "Nochten", + "Welzow-Süd": "Welzow", + } + + model.F_max_month = pyo.Param(model.I, initialize=1e12, mutable=True) + for _, row in cap_month.iterrows(): + i = i_map.get(row["tagebau"]) + if i in model.I: + model.F_max_month[i] = float(row["maximal"]) + + def F_month_rule(m, i, mm): + weeks_in_m = [w for w, mth in week_to_month.items() if mth == mm and w in m.W] + return sum(m.x[i, j, w, d, s] for j in m.J for w in weeks_in_m for d in m.D for s in m.S) + + model.F_month = pyo.Expression(model.I, model.M, rule=F_month_rule) + + # Constraint: Monthly production cap per mine. + 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. + 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 + ) + + verladung_cap = tables["verladungskap"] + + def cap_lookup(verladung_label: str, zeitraum_label: str) -> float: + row = verladung_cap[ + (verladung_cap["verladung"] == verladung_label) & (verladung_cap["zeitraum"] == zeitraum_label) + ] + if row.empty: + raise ValueError(f"Keine Verladungskap für {verladung_label} / {zeitraum_label}") + return float(row.iloc[0]["maximal"]) + + BOXBERG_SOURCES = ["Reichwalde", "Nochten"] + BOXBERG_DEST = ["J", "SP", "B3", "V"] + BOXBERG_SHIFT_CAP = cap_lookup("Boxberg (RW+NO)", "pro Schicht") + BOXBERG_DAY_CAP = cap_lookup("Boxberg (RW+NO)", "pro Tag") + + WELZOW_SOURCES = ["Welzow"] + WELZOW_DEST = ["J", "SP", "B3", "V"] + WELZOW_SHIFT_CAP = cap_lookup("Welzow-Süd", "pro Schicht") + WELZOW_DAY_CAP = cap_lookup("Welzow-Süd", "pro Tag") + + 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. + 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): + 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 + for s in m.S + ) <= BOXBERG_DAY_CAP + + # Constraint: Boxberg loading capacity per day. + 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. + 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): + 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 + for s in m.S + ) <= WELZOW_DAY_CAP + + # Constraint: Welzow loading capacity per day. + model.cap_welzow_day = pyo.Constraint(model.W, model.D, rule=welzow_day_cap_rule) + + 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"} + schicht_order = ["F", "S", "N"] + + model.cap_Welzow = pyo.Param( + model.W, model.D, model.S, initialize=math.nan, mutable=True, within=pyo.Any + ) + model.cap_RW_N = pyo.Param( + model.W, model.D, model.S, initialize=math.nan, mutable=True, within=pyo.Any + ) + + for _, row in avail.iterrows(): + datum = row["datum"] + w = _mining_week_from_date(datum) + d = day_map[datum.strftime("%a")] + + caps_w = [row["Welzow_Sued_S1_t"], row["Welzow_Sued_S2_t"], row["Welzow_Sued_S3_t"]] + caps_rn = [row["Boxberg_NO_RW_S1_t"], row["Boxberg_NO_RW_S2_t"], row["Boxberg_NO_RW_S3_t"]] + + if all((not math.isnan(v)) and v == 0 for v in caps_w): + ih_dates_welzow.add(pd.to_datetime(datum).normalize()) + if all((not math.isnan(v)) and v == 0 for v in caps_rn): + ih_dates_boxberg.add(pd.to_datetime(datum).normalize()) + + for idx, s in enumerate(schicht_order): + val_w = caps_w[idx] + if not math.isnan(val_w): + model.cap_Welzow[w, d, s] = val_w + val_rn = caps_rn[idx] + if not math.isnan(val_rn): + model.cap_RW_N[w, d, s] = val_rn + + def cap_welzow_rule(m, w, d, s): + v = value(m.cap_Welzow[w, d, s]) + if math.isnan(v): + return pyo.Constraint.Skip + return ( + sum( + m.x[i, j, w, d, s] + for i in ["Welzow"] + for j in m.J + if (i, j, w, d, s) in m.x.index_set() + ) + <= v + ) + + # Constraint: Dynamic availability cap for Welzow (per shift). + 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): + v = value(m.cap_RW_N[w, d, s]) + if math.isnan(v): + return pyo.Constraint.Skip + return ( + sum( + m.x[i, j, w, d, s] + for i in ["Reichwalde", "Nochten"] + for j in m.J + if (i, j, w, d, s) in m.x.index_set() + ) + <= v + ) + + # Constraint: Dynamic availability cap for Reichwalde + Nochten (per shift). + model.cap_rw_n_con = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_rule) + + zug = tables.get("zugdurchlass") + if zug is not None: + df = zug.copy() + + def cap_from_table(start_exact: str, ziel_exact: str) -> float: + s = df["start"].astype(str).str.strip() + z = df["ziel"].astype(str).str.strip() + mask = (s == start_exact) & (z == ziel_exact) + row = df[mask] + if row.empty: + raise ValueError(f"cap not found for start={start_exact!r}, ziel={ziel_exact!r}") + return float(row.iloc[0]["maximal"]) + + CAP_RW_N_TO_J = cap_from_table("Verladung Boxberg (KLP)", "KW Jänschwalde") + CAP_RW_N_TO_SP = cap_from_table("Verladung Boxberg (KLP)", "KW Schwarze Pumpe") + CAP_RW_N_TO_V = cap_from_table("Verladung Boxberg (KLP)", "Veredlung ISP") + CAP_RW_N_TO_B3 = cap_from_table("Verladung Boxberg (KLP)", "KW Boxberg Werk 3") + + CAP_W_TO_J = cap_from_table("Verladung Welzow-Süd (KUP)", "KW Jänschwalde") + CAP_W_TO_SP = cap_from_table("Verladung Welzow-Süd (KUP)", "KW Schwarze Pumpe") + CAP_W_TO_V = cap_from_table("Verladung Welzow-Süd (KUP)", "Veredlung ISP") + CAP_W_TO_B3 = cap_from_table("Verladung Welzow-Süd (KUP)", "KW Boxberg Werk 3") + + CAP_RW_N_TO_SP_V = cap_from_table("Verladung Boxberg (KLP)", "KW SP + V ISP") + CAP_RW_N_TO_ALL = cap_from_table("Verladung Boxberg (KLP)", "KW JW + KW SP + V ISP + KW B3") + + CAP_W_TO_SP_V = cap_from_table("Verladung Welzow-Süd (KUP)", "KW SP + V ISP") + CAP_W_TO_SP_V_B3 = cap_from_table("Verladung Welzow-Süd (KUP)", "KW SP + V ISP + KW B3") + + CAP_RW_N_W_TO_J = cap_from_table("KUP + KLP", "KW Jänschwalde") + CAP_RW_N_W_TO_B3 = cap_from_table("KUP + KLP", "KW Boxberg Werk 3") + + CAP_RW_N_J_AND_W_B3 = cap_from_table("KLP zum KW JW", "KUP zum KW B3") + CAP_RW_N_JSPV_AND_W_B3 = cap_from_table("KLP zum KW JW + KW SP + V ISP", "KUP zum KW B3") + + 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. + 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. + 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. + 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. + 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. + 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) + + 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. + 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. + 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. + 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. + 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): + return ( + m.x["Reichwalde", "SP", w, d, s] + + m.x["Nochten", "SP", w, d, s] + + m.x["Reichwalde", "V", w, d, s] + + m.x["Nochten", "V", w, d, s] + <= CAP_RW_N_TO_SP_V + ) + + # Constraint: KLP combined capacity to SP + 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): + return ( + m.x["Reichwalde", "J", w, d, s] + + m.x["Reichwalde", "SP", w, d, s] + + m.x["Reichwalde", "V", w, d, s] + + m.x["Reichwalde", "B3", w, d, s] + + m.x["Nochten", "J", w, d, s] + + m.x["Nochten", "SP", w, d, s] + + m.x["Nochten", "V", w, d, s] + + m.x["Nochten", "B3", w, d, s] + <= CAP_RW_N_TO_ALL + ) + + # Constraint: KLP combined capacity to J + SP + V + 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. + 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): + return ( + m.x["Welzow", "SP", w, d, s] + + m.x["Welzow", "V", w, d, s] + + m.x["Welzow", "B3", w, d, s] + <= CAP_W_TO_SP_V_B3 + ) + + # Constraint: KUP combined capacity to SP + V + 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): + return ( + m.x["Reichwalde", "J", w, d, s] + + m.x["Nochten", "J", w, d, s] + + m.x["Welzow", "J", w, d, s] + <= CAP_RW_N_W_TO_J + ) + + # Constraint: Combined KLP/KUP capacity to J. + 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"] + + model.cap_KVB_Nord = pyo.Param( + model.W, + model.D, + model.S, + initialize=math.nan, + mutable=True, + within=pyo.Any, + ) + + for _, row in kvb_nord.iterrows(): + datum = row["datum"] + w = _mining_week_from_date(datum) + d = day_map[datum.strftime("%a")] + + caps = [row["KVB_Nord_S1_t"], row["KVB_Nord_S2_t"], row["KVB_Nord_S3_t"]] + for idx, s in enumerate(schicht_order): + val = caps[idx] + if not math.isnan(val): + model.cap_KVB_Nord[w, d, s] = val + + def cap_kvb_nord_rule(m, w, d, s): + v = value(m.cap_KVB_Nord[w, d, s]) + if math.isnan(v): + return pyo.Constraint.Skip + return ( + m.x["Reichwalde", "J", w, d, s] + + m.x["Nochten", "J", w, d, s] + + m.x["Welzow", "J", w, d, s] + <= v + ) + + # Constraint: KVB Nord short-term capacity to J (per day/shift if provided). + 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): + return ( + m.x["Reichwalde", "B3", w, d, s] + + m.x["Nochten", "B3", w, d, s] + + m.x["Welzow", "B3", w, d, s] + <= CAP_RW_N_W_TO_B3 + ) + + # Constraint: Combined KLP/KUP capacity to B3. + 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): + return ( + m.x["Reichwalde", "J", w, d, s] + + m.x["Nochten", "J", w, d, s] + + m.x["Welzow", "B3", w, d, s] + <= CAP_RW_N_J_AND_W_B3 + ) + + # Constraint: Combined capacity (KLP->J, KUP->B3). + 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): + return ( + m.x["Reichwalde", "J", w, d, s] + + m.x["Reichwalde", "SP", w, d, s] + + m.x["Reichwalde", "V", w, d, s] + + m.x["Nochten", "J", w, d, s] + + m.x["Nochten", "SP", w, d, s] + + m.x["Nochten", "V", w, d, s] + + m.x["Welzow", "B3", w, d, s] + <= CAP_RW_N_JSPV_AND_W_B3 + ) + + # Constraint: Combined capacity (KLP->J/SP/V, KUP->B3). + 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 + ) + + model._wd_to_date = wd_to_date + model._week_to_month = week_to_month + + return model + + +def solve_model( + model: pyo.ConcreteModel, + solver_name: str = "gurobi", + time_limit: int = 600, + mip_gap: float = 0.05, + iis_path: Path = Path("results/iis.ilp"), +) -> pyo.SolverResults: + if solver_name == "highs": + solver = pyo.SolverFactory("highs") + solver.options.update( + { + "time_limit": time_limit, + "mip_rel_gap": mip_gap, + "threads": 8, + "output_flag": True, + } + ) + return solver.solve(model, tee=True) + + opt = pyo.SolverFactory(solver_name) + if solver_name == "gurobi": + 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 + results = opt.solve(model, tee=True, symbolic_solver_labels=True) + + if solver_name == "gurobi" and results.solver.termination_condition in { + TerminationCondition.infeasible, + TerminationCondition.infeasibleOrUnbounded, + }: + try: + iis_path.parent.mkdir(parents=True, exist_ok=True) + persistent = pyo.SolverFactory("gurobi_persistent") + persistent.set_instance(model, symbolic_solver_labels=True) + persistent.options["TimeLimit"] = time_limit + persistent.options["MIPGap"] = mip_gap + persistent.solve(tee=False) + persistent._solver_model.computeIIS() + persistent._solver_model.write(str(iis_path)) + except Exception as exc: + print(f"WARNING: Failed to write IIS file: {exc}") + + return results diff --git a/src/optimization/run_optimization.py b/src/optimization/run_optimization.py new file mode 100644 index 0000000..e15588d --- /dev/null +++ b/src/optimization/run_optimization.py @@ -0,0 +1,731 @@ +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +import pandas as pd +import pyomo.environ as pyo +from pyomo.environ import value + +SRC_ROOT = Path(__file__).resolve().parents[1] +if str(SRC_ROOT) not in sys.path: + sys.path.insert(0, str(SRC_ROOT)) + +from optimization.model_builder import build_model, load_tables, solve_model + + + + +def report_results(model: pyo.ConcreteModel, max_rows: int) -> None: + print("Objective value:", value(model.obj)) + print("Non-zero production decisions (k):") + printed = 0 + for i in model.I: + for j in model.J: + for w in model.W: + for d in model.D: + for s in model.S: + qty = value(model.k[i, j, w, d, s]) + if qty > 1e-6: + print(f" {i} -> {j} (W{w} D{d} S{s}): {qty:.0f}") + printed += 1 + if printed >= max_rows: + print(" ... output truncated ...") + return + + +def export_results(model: pyo.ConcreteModel, output_path: Path) -> None: + output_path.parent.mkdir(parents=True, exist_ok=True) + wd_to_date = getattr(model, "_wd_to_date", {}) + + def safe_value(var) -> float: + val = pyo.value(var, exception=False) + return float(val) if val is not None else 0.0 + + def autosize_worksheet(ws, df, index_cols=None, max_width=25): + if index_cols is None: + index_cols = list(df.index.names) + idx_names = list(index_cols) + col_widths = [max(10, max(len(str(n)) for n in idx_names if n is not None) + 2) if idx_names else 10] + for col_idx, col in enumerate(df.columns): + header = " / ".join([str(c) for c in col]) if isinstance(col, tuple) else str(col) + max_len = max(len(header), 8) + sample = df.iloc[:200, col_idx] + max_len = max(max_len, sample.astype(str).str.len().max()) + col_widths.append(min(max_width, int(max_len) + 2)) + return col_widths + + def adjust_widths_for_labels(df, widths, label_scale, index_scale=None): + adjusted = widths[:] + if index_scale is not None and adjusted: + adjusted[0] = max(10, int(adjusted[0] * index_scale)) + if hasattr(df.columns, "get_level_values"): + top_level = df.columns.get_level_values(0) + for idx, label in enumerate(top_level, start=1): + if label in label_scale: + adjusted[idx] = max(6, int(adjusted[idx] * label_scale[label])) + return adjusted + lieferungen_schicht = [] + for j in model.J: + for w in model.W: + for d in model.D: + for s in model.S: + if j == "V": + nachfrage = pyo.value(model.dV_N[w, d] + model.dV_W[w, d]) + else: + nachfrage = pyo.value(model.d[j, w, d]) + use_bunker_out = hasattr(model, "bunker_out") and j in getattr(model, "J_BUNKER", []) + delivery_sum = sum(safe_value(model.x[i, j, w, d, s]) for i in model.I) + out_sum = ( + sum(safe_value(model.bunker_out[i, j, w, d, s]) for i in model.I) + if use_bunker_out + else delivery_sum + ) + bunker_inflow = 0.0 + if use_bunker_out: + x_sum = sum(safe_value(model.x[i, j, w, d, s]) for i in model.I) + out_sum = sum(safe_value(model.bunker_out[i, j, w, d, s]) for i in model.I) + bunker_inflow = round(x_sum - out_sum, 2) + def flow_val(i_name: str) -> float: + return safe_value(model.x[i_name, j, w, d, s]) + lieferungen_schicht.append( + { + "kraftwerk": j, + "woche": w, + "tag": d, + "datum": wd_to_date.get((w, d)), + "schicht": s, + "nachfrage_tonnen": nachfrage, + "lieferung_tonnen": delivery_sum, + "lieferungsabweichung_tonnen": round(delivery_sum - nachfrage, 2), + "bunkerzufluss_tonnen": bunker_inflow, + "Nochten": flow_val("Nochten"), + "Reichwalde": flow_val("Reichwalde"), + "Welzow": flow_val("Welzow"), + } + ) + + order_k_pw = ["J", "SP", "B3", "B4"] + order_k_v = ["V"] + order_sources = ["Reichwalde", "Nochten", "Welzow"] + order_s = ["F", "S", "N"] + + df_raw = pd.DataFrame(lieferungen_schicht).copy() + df_raw["datum"] = pd.to_datetime(df_raw["datum"]) + + v_demand_map = { + (int(w), d): { + "welzow": float(pyo.value(model.dV_W[w, d])), + "nochten": float(pyo.value(model.dV_N[w, d])), + } + for w in model.W + for d in model.D + } + + df_raw["nachfrage_welzow"] = df_raw.get("nachfrage_welzow", pd.Series(index=df_raw.index)) + df_raw["nachfrage_nochten"] = df_raw.get("nachfrage_nochten", pd.Series(index=df_raw.index)) + df_raw["nachfrage_welzow"] = df_raw.apply( + lambda r: ( + r["nachfrage_welzow"] + if pd.notna(r["nachfrage_welzow"]) + else v_demand_map.get((int(r["woche"]), r["tag"]), {}).get("welzow", 0) + ), + axis=1, + ) + df_raw["nachfrage_nochten"] = df_raw.apply( + lambda r: ( + r["nachfrage_nochten"] + if pd.notna(r["nachfrage_nochten"]) + else v_demand_map.get((int(r["woche"]), r["tag"]), {}).get("nochten", 0) + ), + axis=1, + ) + + df = df_raw.rename(columns={"lieferungen_tonnen": "lieferung_tonnen"}).copy() + if "lieferung_tonnen" not in df.columns: + df["lieferung_tonnen"] = df[order_sources].sum(axis=1) + + present_sources = [c for c in order_sources if c in df.columns] + + df_src = ( + df.pivot_table( + index=["datum", "woche", "tag"], + columns=["kraftwerk", "schicht"], + values=present_sources, + aggfunc="sum", + ) + .fillna(0) + ) + df_src.columns = df_src.columns.reorder_levels([1, 0, 2]) + df_src = df_src.reindex( + columns=pd.MultiIndex.from_product([order_k_pw + order_k_v, present_sources, order_s]), + fill_value=0, + ) + + df_demand = ( + df.groupby(["datum", "woche", "tag", "kraftwerk"])["nachfrage_tonnen"] + .first() + .unstack("kraftwerk") + .reindex(columns=order_k_pw + order_k_v, fill_value=0) + .reindex(df_src.index, fill_value=0) + ) + + df_v_demand_split = ( + df[df["kraftwerk"] == "V"] + .groupby(["datum", "woche", "tag"])[["nachfrage_welzow", "nachfrage_nochten"]] + .first() + .reindex(df_src.index, fill_value=0) + ) + + totals_plain = df_src.T.groupby(level=0).sum().T + totals_plain = totals_plain.reindex(columns=order_k_pw + order_k_v, fill_value=0) + + day_diff_rows = [] + for w in model.W: + for d in model.D: + date = wd_to_date.get((w, d)) + if date is None: + continue + row = {"datum": date, "woche": w, "tag": d} + for j in order_k_pw + order_k_v: + if j == "V": + demand = pyo.value(model.dV_N[w, d] + model.dV_W[w, d]) + else: + demand = pyo.value(model.d[j, w, d]) + delivered = pyo.value(model.y_delivery[j, w, d]) + row[j] = delivered - demand + day_diff_rows.append(row) + + day_diff_plain = ( + pd.DataFrame(day_diff_rows) + .set_index(["datum", "woche", "tag"]) + .reindex(totals_plain.index, fill_value=0) + .reindex(columns=order_k_pw + order_k_v, fill_value=0) + ) + + totals = totals_plain.copy() + totals.columns = pd.MultiIndex.from_tuples([(k, "Gesamt", "") for k in totals.columns]) + + demand_cols = df_demand.copy() + demand_cols.columns = pd.MultiIndex.from_tuples([(k, "Nachfrage", "") for k in demand_cols.columns]) + + day_diff_cols = day_diff_plain.copy() + day_diff_cols.columns = pd.MultiIndex.from_tuples( + [(k, "Lieferungstagesabweichung", "") for k in day_diff_cols.columns] + ) + + v_demand_cols = df_v_demand_split.copy() + v_demand_cols.columns = pd.MultiIndex.from_tuples( + [ + ("V", "Nachfrage_Welzow", ""), + ("V", "Nachfrage_Nochtener", ""), + ] + ) + + has_bunker = hasattr(model, "bunker") + col_order = [] + for k in order_k_pw: + for src in present_sources: + for sch in order_s: + col_order.append((k, src, sch)) + col_order.append((k, "Gesamt", "")) + col_order.append((k, "Nachfrage", "")) + col_order.append((k, "Lieferungstagesabweichung", "")) + + col_order += [( "V", src, sch) for src in present_sources for sch in order_s] + col_order += [ + ("V", "Nachfrage_Welzow", ""), + ("V", "Nachfrage_Nochtener", ""), + ("V", "Gesamt", ""), + ("V", "Nachfrage", ""), + ("V", "Lieferungstagesabweichung", ""), + ] + + df_out = pd.concat([df_src, v_demand_cols, totals, demand_cols, day_diff_cols], axis=1) + + df_out = df_out.reindex( + columns=col_order, + fill_value=0, + ) + + weekday_order = ["Mo", "Di", "Mi", "Do", "Fr", "Sa", "So"] + idx = df_out.index + df_out = df_out.copy() + df_out.index = pd.MultiIndex.from_arrays( + [ + pd.to_datetime(idx.get_level_values("datum")), + idx.get_level_values("woche"), + pd.Categorical(idx.get_level_values("tag"), categories=weekday_order, ordered=True), + ], + names=["datum", "woche", "tag"], + ) + df_out = df_out.sort_index(level=["datum", "woche", "tag"]) + + df = df_out.copy() + df_out = df_out / 1000 + + bunker_sheet = None + if has_bunker: + bunker_rows = [] + for j in getattr(model, "J_BUNKER", []): + for w in model.W: + for d in model.D: + date = wd_to_date.get((w, d)) + if date is None: + continue + bunker_total = sum(safe_value(model.bunker[i, j, w, d]) for i in model.I) + bunker_rows.append( + { + "kraftwerk": j, + "woche": w, + "tag": d, + "datum": date, + "bunkerbestand_tonnen": bunker_total, + } + ) + + if bunker_rows: + bunker_df = pd.DataFrame(bunker_rows) + bunker_df["datum"] = pd.to_datetime(bunker_df["datum"]) + bunker_df = bunker_df.sort_values(["kraftwerk", "datum", "woche", "tag"]) + bunker_df["vortags_bunkerbestand_tonnen"] = bunker_df.groupby("kraftwerk")[ + "bunkerbestand_tonnen" + ].shift(1) + bunker_df["vortags_bunkerbestand_tonnen"] = bunker_df["vortags_bunkerbestand_tonnen"].fillna(0) + bunker_pivot = ( + bunker_df.pivot_table( + index=["datum", "woche", "tag"], + columns=["kraftwerk"], + values=["bunkerbestand_tonnen"], + aggfunc="first", + ) + .fillna(0) + ) + bunker_pivot = bunker_pivot.reindex(columns=order_k_pw + order_k_v, level=1, fill_value=0) + bunker_pivot.columns = pd.MultiIndex.from_tuples( + [(k, "Bunkerbestand", "") for k in bunker_pivot.columns.get_level_values(1)] + ) + bunker_prev_pivot = ( + bunker_df.pivot_table( + index=["datum", "woche", "tag"], + columns=["kraftwerk"], + values=["vortags_bunkerbestand_tonnen"], + aggfunc="first", + ) + .fillna(0) + ) + bunker_prev_pivot = bunker_prev_pivot.reindex(columns=order_k_pw + order_k_v, level=1, fill_value=0) + bunker_prev_pivot.columns = pd.MultiIndex.from_tuples( + [(k, "Vortags_Bunkerbestand", "") for k in bunker_prev_pivot.columns.get_level_values(1)] + ) + inflow_pivot = None + if "bunkerzufluss_tonnen" in df_raw.columns: + inflow_pivot = ( + df_raw.pivot_table( + index=["datum", "woche", "tag"], + columns=["kraftwerk"], + values=["bunkerzufluss_tonnen"], + aggfunc="sum", + ) + .fillna(0) + ) + inflow_pivot.columns = pd.MultiIndex.from_tuples( + [(k, "Bunkerzufluss", "") for k in inflow_pivot.columns.get_level_values(1)] + ) + + frames = [df] + frames.append(bunker_prev_pivot.reindex(df.index, fill_value=0)) + if inflow_pivot is not None: + frames.append(inflow_pivot.reindex(df.index, fill_value=0)) + frames.append(bunker_pivot.reindex(df.index, fill_value=0)) + bunker_sheet = pd.concat(frames, axis=1) + + col_order_bunker = [] + for k in order_k_pw: + for src in present_sources: + for sch in order_s: + col_order_bunker.append((k, src, sch)) + col_order_bunker.append((k, "Gesamt", "")) + col_order_bunker.append((k, "Nachfrage", "")) + col_order_bunker.append((k, "Lieferungstagesabweichung", "")) + if (k, "Vortags_Bunkerbestand", "") in bunker_sheet.columns: + col_order_bunker.append((k, "Vortags_Bunkerbestand", "")) + if (k, "Bunkerzufluss", "") in bunker_sheet.columns: + col_order_bunker.append((k, "Bunkerzufluss", "")) + if (k, "Bunkerbestand", "") in bunker_sheet.columns: + col_order_bunker.append((k, "Bunkerbestand", "")) + + col_order_bunker += [("V", src, sch) for src in present_sources for sch in order_s] + col_order_bunker += [ + ("V", "Nachfrage_Welzow", ""), + ("V", "Nachfrage_Nochtener", ""), + ("V", "Gesamt", ""), + ("V", "Nachfrage", ""), + ("V", "Lieferungstagesabweichung", ""), + ] + if ("V", "Vortags_Bunkerbestand", "") in bunker_sheet.columns: + col_order_bunker.append(("V", "Vortags_Bunkerbestand", "")) + if ("V", "Bunkerzufluss", "") in bunker_sheet.columns: + col_order_bunker.append(("V", "Bunkerzufluss", "")) + if ("V", "Bunkerbestand", "") in bunker_sheet.columns: + col_order_bunker.append(("V", "Bunkerbestand", "")) + + bunker_sheet = bunker_sheet.reindex(columns=col_order_bunker, fill_value=0) + + try: + import xlsxwriter # type: ignore + + excel_engine = "xlsxwriter" + except Exception: + excel_engine = "openpyxl" + + with pd.ExcelWriter(output_path, engine=excel_engine) as writer: + df.to_excel(writer, sheet_name="Sheet1") + if excel_engine == "xlsxwriter": + ws1 = writer.sheets["Sheet1"] + widths = autosize_worksheet(ws1, df) + widths = adjust_widths_for_labels( + df, + widths, + { + "Reichwalde": 0.5, + "Nochten": 0.5, + "Welzow": 0.5, + "Gesamt": 0.5, + "Nachfrage": 0.5, + "Bunkerbestand": 0.5, + "Bunkerzufluss": 0.5, + "Vortags_Bunkerbestand": 0.5, + }, + index_scale=1.2, + ) + for i, w in enumerate(widths): + ws1.set_column(i, i, w) + else: + ws1 = writer.sheets["Sheet1"] + widths = autosize_worksheet(ws1, df) + widths = adjust_widths_for_labels( + df, + widths, + { + "Reichwalde": 0.5, + "Nochten": 0.5, + "Welzow": 0.5, + "Gesamt": 0.5, + "Nachfrage": 0.5, + "Bunkerbestand": 0.5, + "Bunkerzufluss": 0.5, + "Vortags_Bunkerbestand": 0.5, + }, + index_scale=1.2, + ) + for i, w in enumerate(widths): + ws1.column_dimensions[chr(65 + i)].width = w + + if bunker_sheet is not None: + bunker_sheet.to_excel(writer, sheet_name="mit_Bunkerbestand") + if excel_engine == "xlsxwriter": + workbook = writer.book + worksheet = writer.sheets["mit_Bunkerbestand"] + widths = autosize_worksheet(worksheet, bunker_sheet) + widths = adjust_widths_for_labels( + bunker_sheet, + widths, + { + "Reichwalde": 0.5, + "Nochten": 0.5, + "Welzow": 0.5, + "Gesamt": 0.5, + "Nachfrage": 0.5, + "Bunkerbestand": 0.5, + "Bunkerzufluss": 0.5, + "Vortags_Bunkerbestand": 0.5, + }, + index_scale=1.2, + ) + for i, w in enumerate(widths): + worksheet.set_column(i, i, w) + + header_fmt = workbook.add_format({"bold": True, "bg_color": "#E6E6E6", "border": 1}) + block_colors = ["#DCEFFE", "#FDEBD0", "#E8F8F5", "#FADBD8", "#E8DAEF", "#FEF9E7"] + block_formats = [ + workbook.add_format({"bold": True, "bg_color": color, "border": 1}) for color in block_colors + ] + + index_cols = len(bunker_sheet.index.names) + n_header_rows = bunker_sheet.columns.nlevels + + # Base header formatting for index columns. + for r in range(n_header_rows): + for c in range(index_cols): + worksheet.write(r, c, "", header_fmt) + + # Apply block colors per Kraftwerk on header rows. + top_level = bunker_sheet.columns.get_level_values(0) + blocks = {} + for idx, label in enumerate(top_level): + blocks.setdefault(label, []).append(idx) + + for b_idx, (label, cols) in enumerate(blocks.items()): + fmt = block_formats[b_idx % len(block_formats)] + for r in range(n_header_rows): + for c in cols: + value = bunker_sheet.columns.get_level_values(r)[c] + worksheet.write(r, index_cols + c, value, fmt) + else: + from openpyxl.styles import Border, Font, PatternFill, Side + + worksheet = writer.sheets["mit_Bunkerbestand"] + widths = autosize_worksheet(worksheet, bunker_sheet) + widths = adjust_widths_for_labels( + bunker_sheet, + widths, + { + "Reichwalde": 0.5, + "Nochten": 0.5, + "Welzow": 0.5, + "Gesamt": 0.5, + "Nachfrage": 0.5, + "Bunkerbestand": 0.5, + "Bunkerzufluss": 0.5, + "Vortags_Bunkerbestand": 0.5, + }, + index_scale=1.2, + ) + for i, w in enumerate(widths): + worksheet.column_dimensions[chr(65 + i)].width = w + header_fill = PatternFill("solid", fgColor="E6E6E6") + block_colors = ["DCEFFE", "FDEBD0", "E8F8F5", "FADBD8", "E8DAEF", "FEF9E7"] + block_fills = [PatternFill("solid", fgColor=c) for c in block_colors] + bold_font = Font(bold=True) + border = Border( + left=Side(style="thin"), + right=Side(style="thin"), + top=Side(style="thin"), + bottom=Side(style="thin"), + ) + + index_cols = len(bunker_sheet.index.names) + n_header_rows = bunker_sheet.columns.nlevels + top_level = bunker_sheet.columns.get_level_values(0) + blocks = {} + for idx, label in enumerate(top_level): + blocks.setdefault(label, []).append(idx) + + for r in range(n_header_rows): + for c in range(index_cols): + cell = worksheet.cell(row=r + 1, column=c + 1) + cell.fill = header_fill + cell.font = bold_font + cell.border = border + + for b_idx, (label, cols) in enumerate(blocks.items()): + fill = block_fills[b_idx % len(block_fills)] + for r in range(n_header_rows): + for c in cols: + cell = worksheet.cell(row=r + 1, column=index_cols + c + 1) + cell.fill = fill + cell.font = bold_font + cell.border = border + + # Kohlesorten-Mischverhaeltnis (gesamter Zeitraum) + j_name_map = { + "J": "Jänschwalde", + "SP": "Schwarze Pumpe", + "B3": "Boxberg Werk 3", + "B4": "Boxberg Werk 4", + } + i_name_map = { + "Reichwalde": "Reichwalder-Kohle", + "Nochten": "Nochtener-Kohle", + "Welzow": "Welzower-Kohle", + } + + # Empirical mix over full horizon based on delivered quantities (x). + total_delivered_by_j = {} + for j_code in j_name_map: + if j_code not in model.J: + continue + total_delivered_by_j[j_code] = sum( + safe_value(model.x[i, j_code, w, d, s]) + for i in model.I + for w in model.W + for d in model.D + for s in model.S + if (i, j_code, w, d, s) in model.x + ) + + # Empirical bunker mix over full horizon based on bunker stock. + total_bunker_by_j = {} + if hasattr(model, "bunker"): + for j_code in j_name_map: + if j_code not in getattr(model, "J_BUNKER", []): + continue + total_bunker_by_j[j_code] = sum( + safe_value(model.bunker[i, j_code, w, d]) + for i in model.I + for w in model.W + for d in model.D + if (i, j_code, w, d) in model.bunker + ) + + mix_rows = [] + for j_code, j_name in j_name_map.items(): + if j_code not in model.J: + continue + for i_code, i_name in i_name_map.items(): + if i_code not in model.I: + continue + denom = total_delivered_by_j.get(j_code, 0.0) + num = sum( + safe_value(model.x[i_code, j_code, w, d, s]) + for w in model.W + for d in model.D + for s in model.S + if (i_code, j_code, w, d, s) in model.x + ) + empirisch = round(100 * num / denom, 2) if denom > 0 else 0.0 + bunker_empirisch = 0.0 + if hasattr(model, "bunker") and j_code in total_bunker_by_j: + denom_b = total_bunker_by_j.get(j_code, 0.0) + num_b = sum( + safe_value(model.bunker[i_code, j_code, w, d]) + for w in model.W + for d in model.D + if (i_code, j_code, w, d) in model.bunker + ) + bunker_empirisch = round(100 * num_b / denom_b, 2) if denom_b > 0 else 0.0 + mix_rows.append( + { + "kraftwerk": j_name, + "kohlesorte": i_name, + "ziel_low": round(100 * pyo.value(model.alpha_target_low[i_code, j_code]), 2), + "ziel_high": round(100 * pyo.value(model.alpha_target_high[i_code, j_code]), 2), + "maximal": round(100 * pyo.value(model.alpha_max[i_code, j_code]), 2), + "minimal": round(100 * pyo.value(model.alpha_min[i_code, j_code]), 2), + "empirisch": empirisch, + "bunker_empirisch": bunker_empirisch, + } + ) + + mix_df = pd.DataFrame(mix_rows) + mix_df.to_excel(writer, sheet_name="Kohlemischverhältnis", index=False) + if excel_engine == "xlsxwriter": + ws_mix = writer.sheets["Kohlemischverhältnis"] + workbook = writer.book + mix_block_colors = ["#E8F8F5", "#FDEBD0", "#DCEFFE", "#FADBD8"] + mix_formats = [workbook.add_format({"bg_color": c}) for c in mix_block_colors] + red_fill = workbook.add_format({"bg_color": "#F5B7B1"}) + green_fill = workbook.add_format({"bg_color": "#D4EFDF"}) + if not mix_df.empty: + emp_col = mix_df.columns.get_loc("empirisch") + bunker_emp_col = mix_df.columns.get_loc("bunker_empirisch") + current = None + block_idx = -1 + for r, row in mix_df.iterrows(): + if row["kraftwerk"] != current: + block_idx += 1 + current = row["kraftwerk"] + fmt = mix_formats[block_idx % len(mix_formats)] + for c in range(0, min(6, mix_df.shape[1])): + ws_mix.write(r + 1, c, row.iloc[c], fmt) + emp_fmt = ( + red_fill + if row["empirisch"] < row["ziel_low"] or row["empirisch"] > row["ziel_high"] + else green_fill + ) + ws_mix.write(r + 1, emp_col, row["empirisch"], emp_fmt) + ws_mix.write(r + 1, bunker_emp_col, row["bunker_empirisch"]) + widths = autosize_worksheet(ws_mix, mix_df, index_cols=[]) + for i, w in enumerate(widths[1:]): + ws_mix.set_column(i, i, w) + else: + ws_mix = writer.sheets["Kohlemischverhältnis"] + if not mix_df.empty: + from openpyxl.styles import PatternFill + + mix_block_colors = ["E8F8F5", "FDEBD0", "DCEFFE", "FADBD8"] + mix_fills = [PatternFill("solid", fgColor=c) for c in mix_block_colors] + red_fill = PatternFill("solid", fgColor="F5B7B1") + green_fill = PatternFill("solid", fgColor="D4EFDF") + current = None + block_idx = -1 + for row_idx, kraftwerk in enumerate(mix_df["kraftwerk"], start=2): + if kraftwerk != current: + block_idx += 1 + current = kraftwerk + fill = mix_fills[block_idx % len(mix_fills)] + for col_idx in range(1, min(7, mix_df.shape[1]) + 1): + ws_mix.cell(row=row_idx, column=col_idx).fill = fill + emp_col = mix_df.columns.get_loc("empirisch") + 1 + bunker_emp_col = mix_df.columns.get_loc("bunker_empirisch") + 1 + for r, row in mix_df.iterrows(): + fill = red_fill if row["empirisch"] < row["ziel_low"] or row["empirisch"] > row["ziel_high"] else green_fill + ws_mix.cell(row=r + 2, column=emp_col).fill = fill + # No group fill for bunker_empirisch (column H). + widths = autosize_worksheet(ws_mix, mix_df, index_cols=[]) + for i, w in enumerate(widths[1:]): + ws_mix.column_dimensions[chr(65 + i)].width = w + + +def main() -> None: + parser = argparse.ArgumentParser(description="Run the Pyomo optimization model.") + parser.add_argument( + "--data-dir", + type=Path, + default=Path("data/processed"), + help="Directory containing input parquet files.", + ) + parser.add_argument( + "--solver", + default="gurobi", + help="Solver name passed to Pyomo (default: gurobi).", + ) + parser.add_argument( + "--max-rows", + type=int, + default=50, + help="Maximum number of non-zero decision variable rows to print.", + ) + parser.add_argument( + "--time-limit", + type=int, + default=600, + help="Time limit (seconds) for the solver (default: 600).", + ) + parser.add_argument( + "--output-xlsx", + type=Path, + default=Path("data/out/output.xlsx"), + help="Excel output file for deliveries by plant/week/day/shift.", + ) + parser.add_argument( + "--mip-gap", + type=float, + default=0.03, + help="MIP gap tolerance (default: 0.03).", + ) + parser.add_argument( + "--step-size-tonnes", + type=int, + default=1000, + choices=[960, 1000], + help="Discrete train step size in tonnes (default: 1000).", + ) + args = parser.parse_args() + + tables = load_tables(args.data_dir) + model = build_model(tables, step_size_tonnes=args.step_size_tonnes) + + solve_model(model, args.solver, args.time_limit, args.mip_gap) + # report_results(model, args.max_rows) + export_results(model, args.output_xlsx) + + +if __name__ == "__main__": + main() + + +# uv run python src/optimization/run_optimization.py --solver gurobi --mip-gap 0.05 +# uv run python src/optimization/run_optimization.py --solver highs diff --git a/src/preprocessing/exploration_preprocess.py b/src/preprocessing/exploration_preprocess.py new file mode 100644 index 0000000..3604515 --- /dev/null +++ b/src/preprocessing/exploration_preprocess.py @@ -0,0 +1,505 @@ +# Generated from exploration.ipynb +# %% +import os +from pathlib import Path + +import numpy as np +import pandas as pd + +# %% +PROJECT_ROOT = Path(__file__).resolve().parents[2] +DEFAULT_INPUT = PROJECT_ROOT / "data/input/PoC1_Rohkohleverteilung_Input_Parameter.xlsx" +DEFAULT_OUTPUT = PROJECT_ROOT / "data/processed" +INPUT_XLSX = Path(os.environ.get("POC1_INPUT_XLSX", str(DEFAULT_INPUT))) +OUTPUT_DIR = Path(os.environ.get("POC1_OUTPUT_DIR", str(DEFAULT_OUTPUT))) +OUTPUT_DIR.mkdir(parents=True, exist_ok=True) +path = INPUT_XLSX +# %% [markdown] +# # Mappe Parameter +# %% [markdown] +# ## Erlaubte Abweichung +# ### Kraftwerke +# %% +raw = pd.read_excel(path, sheet_name="Parameter", header=None) + +t1 = raw.iloc[0:18, 0:8].copy() +# 1. komplett leere Zeilen raus +df = t1.dropna(how="all").reset_index(drop=True) + +# 2. Einheitenzeile und Labelzeile identifizieren +unit_idx = df[df.apply(lambda r: r.astype(str).str.contains(r"\[kt\]").any(), axis=1)].index[0] +label_idx = unit_idx + 1 + +unit_row = df.loc[unit_idx] +label_row = df.loc[label_idx] + +# 3. Spaltennamen bauen +headers = [] +for u, l in zip(unit_row, label_row): + if pd.isna(l): + headers.append(str(u)) + elif pd.isna(u): + headers.append(str(l)) + else: + headers.append(f"{l} ({u})") + +# 4. Daten unterhalb der Headerzeilen +data = df.loc[label_idx+1:].reset_index(drop=True) +data.columns = headers + +# Neue Kopie mit eindeutigen Spaltennamen +data2 = data.copy() +data2.columns = [ + "col0", + "titel", + "col2", + "zeitraum", + "minus_kt", + "plus_kt", + "minus_pct", + "plus_pct", +] + +# Titelzeilen entfernen +mask_title = data2["titel"].isin([ + "Erlaubte Abweichungen der Bedarfserfüllung", + "Kraftwerk", + "Kraftwerke", + "Veredlung ISP" +]) +data3 = data2[~mask_title].reset_index(drop=True) + +# kraftwerk steht in col2, nach unten füllen +data3["kraftwerk"] = data3["col2"].ffill() + +# Zeilen ohne Zeitraum entfernen +data3 = data3.dropna(subset=["zeitraum"]).reset_index(drop=True) + +# numerische Spalten konvertieren +for col in ["minus_kt", "plus_kt", "minus_pct", "plus_pct"]: + data3[col] = pd.to_numeric(data3[col], errors="coerce") + +# Endauswahl +result = data3[["kraftwerk", "zeitraum", "minus_kt", "plus_kt", "minus_pct", "plus_pct"]] +data3 = data2.copy() + +data3["kraftwerk"] = data3["col2"].ffill() +data3 = data3.dropna(subset=["zeitraum"]).reset_index(drop=True) + +for col in ["minus_kt", "plus_kt", "minus_pct", "plus_pct"]: + data3[col] = pd.to_numeric(data3[col], errors="coerce") + +bounds_power_plants = data3[["kraftwerk", "zeitraum", "minus_kt", "plus_kt", "minus_pct", "plus_pct"]] +bounds_power_plants[["minus"]] = bounds_power_plants[["minus_kt"]]*1000 +bounds_power_plants[["plus"]] = bounds_power_plants[["plus_kt"]]*1000 +bounds_power_plants.drop(columns=["minus_kt", "plus_kt"],inplace=True) +bounds_power_plants.to_parquet(OUTPUT_DIR / "bounds_power_plants.parquet") +print("Saved bounds_power_plants.parquet") +bounds_power_plants +# %% [markdown] +# ### Veredlung +# %% +raw = pd.read_excel(path, sheet_name="Parameter", header=None) + +t1_upper = raw.iloc[0:4, 0:8].copy() +t1_upper +t2 = raw.iloc[18:28, 0:8].copy() + +t2 = pd.concat([t1_upper,t2], axis=0) +# 1. komplett leere Zeilen raus +df = t2.dropna(how="all").reset_index(drop=True) + +# 2. Einheitenzeile und Labelzeile identifizieren +unit_idx = df[df.apply(lambda r: r.astype(str).str.contains(r"\[kt\]").any(), axis=1)].index[0] +label_idx = unit_idx + 1 + +unit_row = df.loc[unit_idx] +label_row = df.loc[label_idx] + +# 3. Spaltennamen bauen +headers = [] +for u, l in zip(unit_row, label_row): + if pd.isna(l): + headers.append(str(u)) + elif pd.isna(u): + headers.append(str(l)) + else: + headers.append(f"{l} ({u})") + +# 4. Daten unterhalb der Headerzeilen +data = df.loc[label_idx+1:].reset_index(drop=True) +data.columns = headers + +# Neue Kopie mit eindeutigen Spaltennamen +data2 = data.copy() +data2.columns = [ + "col0", + "titel", + "col2", + "zeitraum", + "minus_kt", + "plus_kt", + "minus_pct", + "plus_pct", +] + +# Titelzeilen entfernen +mask_title = data2["titel"].isin([ + "Erlaubte Abweichungen der Bedarfserfüllung", + "Kraftwerk", + "Kraftwerke", + "Veredlung ISP" +]) +# kraftwerk steht in col2, nach unten füllen +data3["kraftwerk"] = data3["col2"].ffill() + +# Zeilen ohne Zeitraum entfernen +data3 = data3.dropna(subset=["zeitraum"]).reset_index(drop=True) + +# numerische Spalten konvertieren +for col in ["minus_kt", "plus_kt", "minus_pct", "plus_pct"]: + data3[col] = pd.to_numeric(data3[col], errors="coerce") + +# Endauswahl +result = data3[["kraftwerk", "zeitraum", "minus_kt", "plus_kt", "minus_pct", "plus_pct"]] +data3 = data2.copy() + +data3["kraftwerk"] = data3["col2"].ffill() +data3 = data3.dropna(subset=["zeitraum"]).reset_index(drop=True) + +for col in ["minus_kt", "plus_kt", "minus_pct", "plus_pct"]: + data3[col] = pd.to_numeric(data3[col], errors="coerce") + +veredelung_bounds = data3[["kraftwerk", "zeitraum", "minus_kt", "plus_kt", "minus_pct", "plus_pct"]] +# veredelung_bounds.rename({"kraftwerk":"kohleart"}, inplace=True) +veredelung_bounds = veredelung_bounds.rename(columns={"kraftwerk":"kohlesorte"}) +veredelung_bounds[["minus"]] = veredelung_bounds[["minus_kt"]]*1000 +veredelung_bounds[["plus"]] = veredelung_bounds[["plus_kt"]]*1000 +veredelung_bounds.drop(columns=["minus_kt", "plus_kt"],inplace=True) +veredelung_bounds.to_parquet(OUTPUT_DIR / "veredelung_bounds.parquet") +print("Saved veredelung_bounds.parquet") +veredelung_bounds +# %% [markdown] +# ## Kohlesorten-Mischverhältnis +# %% + + +raw = pd.read_excel(path, sheet_name="Parameter", header=None) + +# J3:P16 -> Zeilen 2:16, Spalten 9:16 (0-basiert, rechte Grenze exklusiv) +block = raw.iloc[2:16, 9:16].copy() + +# Leere Zeilen entfernen +df = block.dropna(how="all").reset_index(drop=True) + +# Header finden: Einheitenzeile enthält "[%]" +unit_idx = df[df[12].astype(str).str.contains(r"\[%\]", regex=True)].index[0] +label_idx = unit_idx + 1 + +unit_row = df.loc[unit_idx] +label_row = df.loc[label_idx] + +# Spaltenüberschriften konstruieren +headers = [] +for u, l in zip(unit_row, label_row): + if pd.isna(l): + headers.append(str(u)) + elif pd.isna(u): + headers.append(str(l)) + else: + headers.append(f"{l} ({u})") + +# Daten unterhalb der Headerzeilen +data = df.loc[label_idx+1:].reset_index(drop=True) +data.columns = headers + +# Spalten sinnvoll benennen +data2 = data.copy() +data2.columns = [ + "titel", + "kraftwerk", + "kohlesorte", + "ziel_low", + "ziel_high", + "maximal", + "minimal", +] + +# Kraftwerk nach unten füllen +data2["kraftwerk"] = data2["kraftwerk"].ffill() + +# *** WICHTIG: nach kohlesorte filtern, nicht nach titel *** +data3 = data2[data2["kohlesorte"].notna()].reset_index(drop=True) + +# Finaler DataFrame +kohle_mix = data3[[ + "kraftwerk", + "kohlesorte", + "ziel_low", + "ziel_high", + "maximal", + "minimal", +]] + +# numerische Spalten in konsistente floats umwandeln +num_cols = ["ziel_low", "ziel_high", "maximal", "minimal"] +kohle_mix[num_cols] = kohle_mix[num_cols].apply(pd.to_numeric, errors="coerce") + +kohle_mix.to_parquet(OUTPUT_DIR / "kohle_mix.parquet") +print("Saved kohle_mix.parquet") +kohle_mix + +# %% [markdown] +# ## Förderkapazitäten +# %% + +# J19:M23 -> Zeilen 18:23, Spalten 9:13 (0-basiert) +block = pd.read_excel(path, sheet_name="Parameter", header=None).iloc[18:23, 9:13].copy() + +# leere Zeilen entfernen +df = block.dropna(how="all").reset_index(drop=True) + +# Titelzeile "Förderkapazität" rausfiltern +df = df[df[9] != "Förderkapazität"].reset_index(drop=True) + +# Spalten benennen +df.columns = ["kategorie", "tagebau", "zeitraum", "maximal"] + +# "Tagebau" nach unten auffüllen (für Nochten, Gesamt, Welzow-Süd) +df["kategorie"] = df["kategorie"].ffill() + +# Maximalwert in Zahl konvertieren +df["maximal"] = pd.to_numeric(df["maximal"], errors="coerce") + +# falls du nur die wesentlichen Infos brauchst: +foerderkap = df[["tagebau", "zeitraum", "maximal"]] + +foerderkap["maximal"] = foerderkap["maximal"]*1000 +foerderkap.to_parquet(OUTPUT_DIR / "foerderkapaz.parquet") +print("Saved foerderkap.parquet") +foerderkap +# %% [markdown] +# ## Verladungskapazitäten +# %% +raw = pd.read_excel(path, sheet_name="Parameter", header=None) + +# J26:M30 -> rows 25:30, cols 9:13 +block = raw.iloc[25:30, 9:13].copy() + +# vollständig leere Zeilen entfernen +df = block.dropna(how="all").reset_index(drop=True) + +# Titelzeile entfernen +df = df[df[9] != "Verladungskapazität"].reset_index(drop=True) + +# Spalten sinnvoll benennen +df.columns = ["kategorie", "verladung", "zeitraum", "maximal"] + +# Verladung nach unten auffüllen +df["verladung"] = df["verladung"].ffill() + +# Maximalwert numerisch +df["maximal"] = pd.to_numeric(df["maximal"], errors="coerce") + +# finale Auswahl +verladung = df[["verladung", "zeitraum", "maximal"]] +verladung["maximal"] = verladung.maximal*1000 +verladung.to_parquet(OUTPUT_DIR / "verladungskap.parquet") +print("Saved verladungskap.parquet") +verladung +# %% [markdown] +# ## Zugdurchlass +# %% +raw = pd.read_excel(path, sheet_name="Parameter", header=None) +block = raw.iloc[3:21, 17:24].copy() +block.replace("unlimitiert", np.inf, inplace=True) +# komplett leere Zeilen entfernen +df = block.dropna(how="all").reset_index(drop=True) + +# Kopfzeilen / Überschrift entfernen +mask_header = df.apply( + lambda r: r.astype(str).str.contains( + "Zugdurchlasskapazität|Maximal|Vielfaches von", regex=True + ).any(), + axis=1, +) +df = df[~mask_header].reset_index(drop=True) + +# Spalten benennen +df.columns = ["von", "start", "zum", "ziel", "zeitraum", "maximal", "vielfaches_von"] + +# numerische Spalten nach float konvertieren +for col in ["maximal", "vielfaches_von"]: + df[col] = pd.to_numeric(df[col], errors="coerce") + +# alle numerischen Werte * 1000 +df["maximal"] = df["maximal"] * 1000 + +# fertiger DataFrame +zugdurchlass = df.copy() +zugdurchlass.to_parquet(OUTPUT_DIR / "zugdurchlass.parquet") +zugdurchlass +# %% [markdown] +# # Mappe Rohkohlebedarf +# %% +raw = pd.read_excel(path, sheet_name="Rohkohlebedarf", header=None) + +df = raw.iloc[2:36, 1:15].copy().reset_index(drop=True) + +jahr = int(df.loc[0, 2]) +monat = str(df.loc[1, 2]) + +kw_header = df.loc[1, 4:10].tolist() +ver_header = df.loc[1, 12:14].tolist() + +kw_names = kw_header[:-1] + ["Gesamt_KW"] +ver_names = ver_header[:-1] + ["Gesamt_Veredlung"] + +data = df.loc[3:].reset_index(drop=True) + +out = pd.DataFrame() +out["jahr"] = jahr +out["monat"] = monat +out["datum"] = pd.to_datetime(data[1]) + +for idx, name in zip(range(4, 4+len(kw_names)), kw_names): + out[name] = pd.to_numeric(data[idx], errors="coerce") + +for idx, name in zip(range(12, 12+len(ver_names)), ver_names): + out[name] = pd.to_numeric(data[idx], errors="coerce") + +# units: convert kt to t +kw_cols = kw_names +out[kw_cols] = out[kw_cols] * 1000 + +# convert Gesamt_Veredlung from kt to t +welz_col, nocht_col, ges_ver_col = ver_names +out[ges_ver_col] = out[ges_ver_col] * 1000 + +out = pd.DataFrame() +out["datum"] = pd.to_datetime(data[1]) +out["jahr"] = jahr +out["monat"] = monat + +for idx, name in zip(range(4, 4+len(kw_names)), kw_names): + out[name] = pd.to_numeric(data[idx], errors="coerce") + +for idx, name in zip(range(12, 12+len(ver_names)), ver_names): + out[name] = pd.to_numeric(data[idx], errors="coerce") + +kw_cols = kw_names +out[kw_cols] = out[kw_cols] * 1000 +welz_col, nocht_col, ges_ver_col = ver_names +out[ges_ver_col] = out[ges_ver_col] * 1000 + +out.rename(columns={"Welzower Kohle":"Veredel_Welzower", "Nochtener Kohle": "Veredel_Nochtener"}, inplace=True) +out.to_parquet(OUTPUT_DIR / "rohkohlebedarf.parquet") +print("Saved rohkohlebedarf.parquet") +out.round(5) + +# %% [markdown] +# # Mappe Verfügbarkeit +# %% +import pandas as pd + + +raw = pd.read_excel(path, sheet_name="Verfügbarkeit", header=None) + +# Jahr/Monat bleiben statisch in C3/C4 +jahr = int(raw.iloc[2, 2]) +monat = str(raw.iloc[3, 2]) + +# B8:J38 (Datum in Spalte B, Wochentag in Spalte C) +df = raw.iloc[7:38, 1:10].copy().reset_index(drop=True) + +# "Datum"-Zeile entfernen +data = df.copy() +data = data[data[1] != "Datum"].reset_index(drop=True) +data = data[data[1].notna()].reset_index(drop=True) + +ver = pd.DataFrame() +ver["datum"] = pd.to_datetime(data[1]) +ver["jahr"] = jahr +ver["monat"] = monat + +# Welzow-Süd Schicht 1 2 3 in cols 4,5,6 +ver["Welzow_Sued_S1_t"] = pd.to_numeric(data[4], errors="coerce") * 1000 +ver["Welzow_Sued_S2_t"] = pd.to_numeric(data[5], errors="coerce") * 1000 +ver["Welzow_Sued_S3_t"] = pd.to_numeric(data[6], errors="coerce") * 1000 + +# Boxberg (NO+RW) Schicht 1 2 3 in cols 7,8,9 +ver["Boxberg_NO_RW_S1_t"] = pd.to_numeric(data[7], errors="coerce") * 1000 +ver["Boxberg_NO_RW_S2_t"] = pd.to_numeric(data[8], errors="coerce") * 1000 +ver["Boxberg_NO_RW_S3_t"] = pd.to_numeric(data[9], errors="coerce") * 1000 + +# ver.to_parquet +ver.to_parquet(OUTPUT_DIR / "Verfuegbarkeiten.parquet") +print("Saved Verfuegbarkeiten.parquet") +ver.round(5) + +# KVB Nord Zugdurchlasskapazitäten (L8:N38) +kvb_block = raw.iloc[7:38, 11:14].copy().reset_index(drop=True) +kvb_block = kvb_block.iloc[: len(data)].reset_index(drop=True) + +kvb = pd.DataFrame() +kvb["datum"] = pd.to_datetime(data[1]) +kvb["jahr"] = jahr +kvb["monat"] = monat +kvb["KVB_Nord_S1_t"] = pd.to_numeric(kvb_block[11], errors="coerce") * 1000 +kvb["KVB_Nord_S2_t"] = pd.to_numeric(kvb_block[12], errors="coerce") * 1000 +kvb["KVB_Nord_S3_t"] = pd.to_numeric(kvb_block[13], errors="coerce") * 1000 + +kvb.to_parquet(OUTPUT_DIR / "zugdurchlass_kvb_nord.parquet") +print("Saved zugdurchlass_kvb_nord.parquet") + + +# %% [markdown] +# # Bunker + +raw = pd.read_excel(path, sheet_name="Parameter", header=None) + +# Ausschnitt R23:W30 (0-basiert: Zeilen 22–29, Spalten 17–22) +bunker = raw.iloc[22:30, 17:23].copy() + +# Kraftwerks- und Veredlungsbunker (Jänschwalde, SP, BW3, ISP) +plants = bunker.iloc[25-22:29-22].reset_index(drop=True) +plants.columns = [ + "typ", + "anlage", + "anfang_mo_di_kt", + "anfang_rest_kt", + "zielbestand_kt", + "maximal_kt", +] + +# Typ auffüllen +plants["typ"] = plants["typ"].ffill() + +# numerische Werte konvertieren +for col in ["anfang_mo_di_kt", "anfang_rest_kt", "zielbestand_kt", "maximal_kt"]: + plants[col] = pd.to_numeric(plants[col], errors="coerce") + +# *** Umrechnung in t *** +plants = plants.rename(columns={ + "anfang_mo_di_kt": "anfang_mo_di_t", + "anfang_rest_kt": "anfang_rest_t", + "zielbestand_kt": "zielbestand_t", + "maximal_kt": "maximal_t", +}) + +plants[["anfang_mo_di_t", "anfang_rest_t", "zielbestand_t", "maximal_t"]] *= 1000 + +# Vorfahrfenster +vorfahrfenster_tage = pd.to_numeric(bunker.loc[29, 19], errors="coerce") + +plants.to_parquet(OUTPUT_DIR / "bunker.parquet") +print("Saved bunker.parquet") +plants + +pd.DataFrame([{"vorfahrfenster_tage": vorfahrfenster_tage}]).to_parquet( + OUTPUT_DIR / "bunker_vorfahrfenster.parquet" +) +print("Saved bunker_vorfahrfenster.parquet") + +# %% +print("\n ####################### Done preprocessing exploration_preprocess.py ####################### \n")