241 lines
9.7 KiB
Python
241 lines
9.7 KiB
Python
import numpy as np
|
|
import pandas as pd
|
|
from scipy.optimize import Bounds, LinearConstraint, milp
|
|
from scipy.sparse import lil_matrix
|
|
|
|
DT_H = 5.0 / 60.0
|
|
|
|
|
|
def train_artifact(kind):
|
|
return {"trained": True, "type": "battery_48h_cost_milp_v3", "source": kind}
|
|
|
|
|
|
def _cfg_float(config, key, default):
|
|
try:
|
|
value = config.get(key, default)
|
|
return float(default if value is None or value == "" else value)
|
|
except Exception:
|
|
return float(default)
|
|
|
|
|
|
def _cfg_bool(config, key, default=False):
|
|
value = config.get(key, default)
|
|
if value is None or value == "":
|
|
return bool(default)
|
|
if isinstance(value, bool):
|
|
return value
|
|
return str(value).strip().lower() in {"1", "true", "yes", "ja", "on"}
|
|
|
|
|
|
def _use_dynamic(config, key):
|
|
value = str(config.get(key, "") or "").strip().lower()
|
|
return any(token in value for token in ("dynam", "marktpreis", "referenzmarktpreis", "market"))
|
|
|
|
|
|
def _price(data_obj, timestamp, column, fallback, dynamic_enabled):
|
|
if not dynamic_enabled:
|
|
return fallback
|
|
frame = data_obj["df_fut"]
|
|
if column in frame.columns and timestamp in frame.index:
|
|
try:
|
|
value = float(frame.at[timestamp, column])
|
|
if np.isfinite(value):
|
|
return value
|
|
except Exception:
|
|
pass
|
|
return fallback
|
|
|
|
|
|
def _battery_meta(config, data_obj):
|
|
cap_kwh = _cfg_float(config, "batt_capacity_kwh", 0.0)
|
|
max_power_w = _cfg_float(config, "batt_power_kw", 0.0) * 1000.0
|
|
min_soc = _cfg_float(config, "batt_min_soc", _cfg_float(config, "batt_min_soc_percent", 0.0))
|
|
max_soc = _cfg_float(config, "batt_max_soc", _cfg_float(config, "batt_max_soc_percent", 100.0))
|
|
start_soc_value = data_obj.get("current_soc_perc")
|
|
if start_soc_value is None:
|
|
start_soc_value = min_soc
|
|
try:
|
|
start_soc = float(start_soc_value)
|
|
except (TypeError, ValueError):
|
|
start_soc = min_soc
|
|
min_soc = max(0.0, min(100.0, min_soc))
|
|
max_soc = max(min_soc, min(100.0, max_soc))
|
|
start_soc = max(min_soc, min(max_soc, start_soc))
|
|
charge_eff = max(0.01, min(1.0, _cfg_float(config, "batt_charge_efficiency", 0.95)))
|
|
discharge_eff = max(0.01, min(1.0, _cfg_float(config, "batt_discharge_efficiency", 0.95)))
|
|
return cap_kwh, max_power_w, min_soc, max_soc, start_soc, charge_eff, discharge_eff
|
|
|
|
|
|
def _quarter_groups(index):
|
|
groups = {}
|
|
for position, timestamp in enumerate(index):
|
|
quarter = timestamp.floor("15min") if hasattr(timestamp, "floor") else position // 3
|
|
groups.setdefault(quarter, []).append(position)
|
|
return list(groups.values())
|
|
|
|
|
|
def _fallback_plan(index, residual_w):
|
|
grid = {timestamp: float(value) for timestamp, value in zip(index, residual_w)}
|
|
battery = {timestamp: 0.0 for timestamp in index}
|
|
return {"grid": grid, "battery": battery, "solver": "fallback"}
|
|
|
|
|
|
def optimize_battery_plan(data_obj, pv_dict, load_dict):
|
|
config = data_obj["config"]
|
|
cap_kwh, max_power_w, min_soc, max_soc, start_soc, charge_eff, discharge_eff = _battery_meta(config, data_obj)
|
|
index = list(data_obj["df_fut"].index)
|
|
if not index:
|
|
return {"grid": {}, "battery": {}, "solver": "empty"}
|
|
|
|
load_w = np.array([max(0.0, float(load_dict.get(t, 0.0))) for t in index])
|
|
pv_w = np.array([max(0.0, float(pv_dict.get(t, 0.0))) for t in index])
|
|
residual_w = load_w - pv_w
|
|
if cap_kwh <= 0.0 or max_power_w <= 0.0:
|
|
return _fallback_plan(index, residual_w)
|
|
|
|
import_fixed = _cfg_float(config, "tarif_bezug_fest", 0.30)
|
|
export_fixed = _cfg_float(config, "tarif_einspeisung_fest", 0.10)
|
|
import_dynamic = _use_dynamic(config, "tarif_bezug")
|
|
export_dynamic = _use_dynamic(config, "tarif_einspeisung")
|
|
import_price = np.array([
|
|
_price(data_obj, t, "import_price", import_fixed, import_dynamic) for t in index
|
|
])
|
|
export_price = np.array([
|
|
_price(data_obj, t, "export_price", export_fixed, export_dynamic) for t in index
|
|
])
|
|
|
|
n = len(index)
|
|
imp, exp, charge, discharge, curtail, soc = 0, n, 2 * n, 3 * n, 4 * n, 5 * n
|
|
peak = 6 * n + 1
|
|
battery_mode = peak + 1
|
|
grid_mode = battery_mode + n
|
|
variable_count = grid_mode + n
|
|
|
|
max_import_w = max(float(load_w.max(initial=0.0)) + max_power_w, max_power_w, 1.0)
|
|
configured_import_limit = _cfg_float(config, "grid_import_limit_w", 0.0)
|
|
if configured_import_limit > 0.0:
|
|
max_import_w = min(max_import_w, configured_import_limit)
|
|
max_export_w = max(float(pv_w.max(initial=0.0)) + max_power_w, max_power_w, 1.0)
|
|
configured_export_limit = _cfg_float(config, "grid_export_limit_w", 0.0)
|
|
if configured_export_limit > 0.0:
|
|
max_export_w = min(max_export_w, configured_export_limit)
|
|
|
|
lower = np.zeros(variable_count)
|
|
upper = np.full(variable_count, np.inf)
|
|
upper[imp:imp + n] = max_import_w
|
|
upper[exp:exp + n] = max_export_w
|
|
upper[charge:charge + n] = max_power_w
|
|
if not _cfg_bool(config, "batt_grid_charging_enabled", False):
|
|
upper[charge:charge + n] = np.minimum(max_power_w, np.maximum(0.0, pv_w - load_w))
|
|
upper[discharge:discharge + n] = max_power_w
|
|
upper[curtail:curtail + n] = pv_w
|
|
reserve_soc = max(
|
|
min_soc,
|
|
min(100.0, _cfg_float(config, "batt_economic_reserve_soc_percent", 10.0)),
|
|
)
|
|
economic_min_soc = max(min_soc, min(start_soc, reserve_soc))
|
|
min_kwh = cap_kwh * economic_min_soc / 100.0
|
|
max_kwh = cap_kwh * max_soc / 100.0
|
|
lower[soc:soc + n + 1] = min_kwh
|
|
upper[soc:soc + n + 1] = max_kwh
|
|
current_peak_kw = max(0.0, float(data_obj.get("current_month_peak_kw", 0.0) or 0.0))
|
|
lower[peak] = current_peak_kw
|
|
upper[peak] = max(current_peak_kw, max_import_w / 1000.0)
|
|
upper[battery_mode:battery_mode + n] = 1.0
|
|
upper[grid_mode:grid_mode + n] = 1.0
|
|
|
|
objective = np.zeros(variable_count)
|
|
objective[imp:imp + n] = import_price * DT_H / 1000.0
|
|
objective[exp:exp + n] = -export_price * DT_H / 1000.0
|
|
degradation = max(0.0, _cfg_float(config, "batt_degradation_chf_kwh", 0.03))
|
|
objective[charge:charge + n] = (degradation / 2.0 + 1e-7) * DT_H / 1000.0
|
|
objective[discharge:discharge + n] = (degradation / 2.0 + 1e-7) * DT_H / 1000.0
|
|
objective[curtail:curtail + n] = 1e-9 * DT_H / 1000.0
|
|
objective[peak] = max(0.0, _cfg_float(config, "tarif_peak_fest", 0.0))
|
|
terminal_value = _cfg_float(config, "batt_terminal_value_chf_kwh", np.median(import_price))
|
|
objective[soc + n] = -max(0.0, terminal_value) * discharge_eff
|
|
|
|
equality_rows = 2 * n + 1
|
|
equality = lil_matrix((equality_rows, variable_count), dtype=float)
|
|
equality_rhs = np.zeros(equality_rows)
|
|
for i in range(n):
|
|
equality[i, imp + i] = 1.0
|
|
equality[i, exp + i] = -1.0
|
|
equality[i, charge + i] = -1.0
|
|
equality[i, discharge + i] = 1.0
|
|
equality[i, curtail + i] = -1.0
|
|
equality_rhs[i] = residual_w[i]
|
|
|
|
row = n + i
|
|
equality[row, soc + i] = -1.0
|
|
equality[row, soc + i + 1] = 1.0
|
|
equality[row, charge + i] = -charge_eff * DT_H / 1000.0
|
|
equality[row, discharge + i] = DT_H / (1000.0 * discharge_eff)
|
|
equality[2 * n, soc] = 1.0
|
|
equality_rhs[2 * n] = cap_kwh * start_soc / 100.0
|
|
|
|
quarter_groups = _quarter_groups(pd.Index(index))
|
|
inequality_rows = 4 * n + len(quarter_groups)
|
|
inequality = lil_matrix((inequality_rows, variable_count), dtype=float)
|
|
inequality_upper = np.zeros(inequality_rows)
|
|
row = 0
|
|
for i in range(n):
|
|
inequality[row, charge + i] = 1.0
|
|
inequality[row, battery_mode + i] = -max_power_w
|
|
row += 1
|
|
inequality[row, discharge + i] = 1.0
|
|
inequality[row, battery_mode + i] = max_power_w
|
|
inequality_upper[row] = max_power_w
|
|
row += 1
|
|
inequality[row, imp + i] = 1.0
|
|
inequality[row, grid_mode + i] = -max_import_w
|
|
row += 1
|
|
inequality[row, exp + i] = 1.0
|
|
inequality[row, grid_mode + i] = max_export_w
|
|
inequality_upper[row] = max_export_w
|
|
row += 1
|
|
for group in quarter_groups:
|
|
for i in group:
|
|
inequality[row, imp + i] = 1.0 / (len(group) * 1000.0)
|
|
inequality[row, peak] = -1.0
|
|
row += 1
|
|
|
|
integrality = np.zeros(variable_count, dtype=int)
|
|
integrality[battery_mode:battery_mode + n] = 1
|
|
integrality[grid_mode:grid_mode + n] = 1
|
|
constraints = [
|
|
LinearConstraint(equality.tocsr(), equality_rhs, equality_rhs),
|
|
LinearConstraint(inequality.tocsr(), -np.inf, inequality_upper),
|
|
]
|
|
result = milp(
|
|
objective,
|
|
integrality=integrality,
|
|
bounds=Bounds(lower, upper),
|
|
constraints=constraints,
|
|
options={"time_limit": max(5.0, _cfg_float(config, "batt_optimizer_timeout_seconds", 30.0))},
|
|
)
|
|
if not result.success or result.x is None:
|
|
return _fallback_plan(index, residual_w)
|
|
|
|
grid_values = result.x[imp:imp + n] - result.x[exp:exp + n]
|
|
battery_values = result.x[charge:charge + n] - result.x[discharge:discharge + n]
|
|
threshold_w = max(25.0, max_power_w * 0.005)
|
|
grid_values[np.abs(grid_values) < threshold_w] = 0.0
|
|
battery_values[np.abs(battery_values) < threshold_w] = 0.0
|
|
return {
|
|
"grid": {t: float(v) for t, v in zip(index, grid_values)},
|
|
"battery": {t: float(v) for t, v in zip(index, battery_values)},
|
|
"solver": "scipy-milp",
|
|
"objective_chf": float(result.fun),
|
|
"planned_peak_kw": float(result.x[peak]),
|
|
"start_soc_percent": float(start_soc),
|
|
"economic_min_soc_percent": float(economic_min_soc),
|
|
}
|
|
|
|
|
|
def optimize_grid_setpoint(data_obj, pv_dict, load_dict, forecast_id=None):
|
|
plan = optimize_battery_plan(data_obj, pv_dict, load_dict)
|
|
if forecast_id is not None:
|
|
data_obj.setdefault("battery_plans", {})[int(forecast_id)] = plan
|
|
return plan["grid"]
|