From f70a786d7c7648151c7f81af7bea8abc1c36bcaa Mon Sep 17 00:00:00 2001 From: Andreas Date: Wed, 16 Sep 2026 12:17:26 +0200 Subject: [PATCH] feat(optimization): port tested terminal and tail value primitives Source d2e2d58237339454dd8bb226f92677c5987f8b27. 22 primitive tests pass; integration with the optimizer, forecast horizon and API is still pending. Co-authored-by: Andreas Co-authored-by: Christin --- .../optimization/genetic/forecast.py | 44 +++ .../optimization/genetic/tailvalue.py | 303 +++++++++++++++ .../optimization/genetic/terminalvalue.py | 359 ++++++++++++++++++ tests/test_tailvalue_physics.py | 200 ++++++++++ tests/test_terminalvalue.py | 127 +++++++ 5 files changed, 1033 insertions(+) create mode 100644 src/akkudoktoreos/optimization/genetic/forecast.py create mode 100644 src/akkudoktoreos/optimization/genetic/tailvalue.py create mode 100644 src/akkudoktoreos/optimization/genetic/terminalvalue.py create mode 100644 tests/test_tailvalue_physics.py create mode 100644 tests/test_terminalvalue.py diff --git a/src/akkudoktoreos/optimization/genetic/forecast.py b/src/akkudoktoreos/optimization/genetic/forecast.py new file mode 100644 index 00000000..8d7ad7d9 --- /dev/null +++ b/src/akkudoktoreos/optimization/genetic/forecast.py @@ -0,0 +1,44 @@ +"""Resample actual forecast intervals without extrapolating missing provider data.""" + +from typing import Any + +import numpy as np +import pandas as pd + + +def bounded_forecast_array( + prediction: Any, + *, + key: str, + start_datetime: Any, + end_datetime: Any, + interval: Any, + **kwargs: Any, +) -> np.ndarray: + """Hold interval averages only within their source interval. + + EOS forecasts carry interval starts. Infer source cadence from timestamps, + conservatively bounded to one hour; never extend the last value indefinitely. + Explicit NaNs and holes remain missing. Downsampling requires full coverage. + """ + target_seconds = int(interval.total_seconds()) + try: + series = prediction.key_to_series(key, dropna=False) + except KeyError: + series = pd.Series(dtype=float, index=pd.DatetimeIndex([], tz="UTC")) + series = pd.to_numeric(series, errors="coerce").sort_index() + series = series[~series.index.duplicated(keep="last")] + cadence = 3600 + if len(series) > 1: + gaps = np.diff(series.index.as_unit("ns").asi8) / 1e9 + cadence = int(min(3600, np.min(gaps[gaps > 0]))) + step = min(cadence, target_seconds) + index = pd.date_range(start=start_datetime, end=end_datetime, freq=f"{step}s", inclusive="left") + if series.empty: + sampled = pd.Series(np.nan, index=index) + else: + sampled = series.reindex(index, method="ffill", tolerance=pd.Timedelta(seconds=cadence - 1)) + groups = sampled.resample(f"{target_seconds}s", origin=start_datetime) + result = groups.mean() + result[groups.count().to_numpy() < target_seconds / step] = np.nan + return result.to_numpy(dtype=float) diff --git a/src/akkudoktoreos/optimization/genetic/tailvalue.py b/src/akkudoktoreos/optimization/genetic/tailvalue.py new file mode 100644 index 00000000..91a6f371 --- /dev/null +++ b/src/akkudoktoreos/optimization/genetic/tailvalue.py @@ -0,0 +1,303 @@ +"""Chronological, run-local battery lookahead using the production device physics. + +The Bellman recursion interpolates continuation values on a stored-energy grid. +Tail actions never enter the executable control arrays. The selected path may +be returned separately as diagnostics so users can inspect the lookahead. +""" + +import numpy as np +from pydantic import PrivateAttr + +from akkudoktoreos.devices.genetic.battery import Battery +from akkudoktoreos.devices.genetic.inverter import Inverter +from akkudoktoreos.optimization.genetic.terminalvalue import TailPlanSlot, TerminalValueCurve + + +TailAction = tuple[int, int, float, float] + + +def _action_name(action: TailAction) -> str: + dc, discharge, ac_rate, export_rate = action + if export_rate > 0: + return "BATTERY_EXPORT" + if ac_rate > 0: + return "GRID_CHARGE" + if dc and discharge: + return "SELF_CONSUMPTION" + if dc: + return "PV_CHARGE_ONLY" + if discharge: + return "DISCHARGE_ONLY" + return "HOLD" + + +def _simulate_action( + *, + bat: Battery, + inv: Inverter, + energy_wh: float, + action: TailAction, + price: float, + load: float, + pv: float, + tariff: float, + direct_marketing: bool, +) -> dict[str, float]: + """Apply one tail action from one stored-energy state.""" + dc, discharge, ac_rate, export = action + bat.soc_wh = float(energy_wh) + bat._charged_raw_wh_per_slot.fill(0) + bat._discharged_raw_wh_per_slot.fill(0) + ac_enabled = inv.ac_to_dc_efficiency > 0 and ( + inv.max_ac_charge_power_w is None or inv.max_ac_charge_power_w > 0 + ) + bat.charge_array[0] = ac_rate if ac_rate > 0 and ac_enabled else dc + bat.discharge_array[0] = discharge if export == 0 or tariff > 0 else 0 + sold, bought, losses, _ = inv.process_energy( + pv, + load, + 0, + allow_battery_grid_export=direct_marketing and export > 0 and tariff > 0, + battery_grid_export_factor=export, + ) + ac_grid_charge_wh = 0.0 + if ac_rate > 0 and inv.ac_to_dc_efficiency > 0: + rate = ac_rate + if inv.max_ac_charge_power_w is not None and bat.max_charge_power_w > 0: + rate = min( + rate, + inv.max_ac_charge_power_w * inv.ac_to_dc_efficiency / bat.max_charge_power_w, + ) + bat.charge_array[0] = rate + if rate > 0: + stored, loss = bat.charge_energy(None, 0, charge_factor=rate) + ac_grid_charge_wh = (stored + loss) / inv.ac_to_dc_efficiency + bought += ac_grid_charge_wh + losses += loss + max(ac_grid_charge_wh - stored - loss, 0.0) + if direct_marketing and tariff < 0: + sold = 0.0 + discharged_wh = bat.discharged_energy_wh(0) + reward = ( + sold * tariff - bought * price - (discharged_wh * bat.levelized_cost_of_storage_kwh / 1000) + ) + return { + "next_state_wh": bat.soc_wh, + "reward_euro": reward, + "grid_export_wh": sold, + "grid_import_wh": bought, + "battery_charge_wh": (bat._charged_raw_wh_per_slot[0] * bat.charging_efficiency), + "battery_discharge_wh": discharged_wh, + "losses_wh": losses, + "ac_grid_charge_wh": ac_grid_charge_wh, + } + + +class TailValueCurve(TerminalValueCurve): + """Value of usable AC battery energy, including the value of empty capacity. + + Neither values nor marginal values are constrained to be monotone. The empty + state can earn money by charging at negative prices and selling later. + """ + + _trace_context: dict = PrivateAttr(default_factory=dict) + + def value(self, energy_wh: float) -> float: + return ( + float(np.interp(energy_wh, self.energy_wh, self.value_euro)) if self.energy_wh else 0.0 + ) + + def component_values(self, energy_wh: float) -> tuple[float, float]: + """Return tail operating cash flow and continuation credit separately.""" + if not self.energy_wh: + return 0.0, 0.0 + operating = float(np.interp(energy_wh, self.energy_wh, self.operating_value_euro)) + continuation = float(np.interp(energy_wh, self.energy_wh, self.continuation_value_euro)) + return operating, continuation + + def diagnostic_plan(self, energy_wh: float, control_horizon_hours: float) -> list[TailPlanSlot]: + """Replay the optimal tail path for one control-end battery state.""" + context = self._trace_context + if not context: + return [] + bat = Battery( + context["battery_parameters"], + prediction_hours=1, + slot_duration_h=context["slot_duration_h"], + ) + inv = Inverter( + context["inverter_parameters"], + battery=bat, + slot_duration_h=context["slot_duration_h"], + ) + conversion = bat.discharging_efficiency * inv.dc_to_ac_efficiency + state_wh = bat.min_soc_wh + (energy_wh / conversion if conversion > 0 else 0.0) + state_wh = float(np.clip(state_wh, bat.min_soc_wh, bat.max_soc_wh)) + plan: list[TailPlanSlot] = [] + arrays = zip( + context["prices"], + context["load"], + context["pv"], + context["tariffs"], + ) + for slot, (price, load, pv, tariff) in enumerate(arrays): + candidates: list[tuple[float, TailAction, dict[str, float], float]] = [] + for action in context["actions"]: + result = _simulate_action( + bat=bat, + inv=inv, + energy_wh=state_wh, + action=action, + price=price, + load=load, + pv=pv, + tariff=tariff, + direct_marketing=context["direct_marketing"], + ) + remaining = float( + np.interp( + result["next_state_wh"], + context["states"], + context["future_values"][slot + 1], + ) + ) + candidates.append((result["reward_euro"] + remaining, action, result, remaining)) + candidates.sort(key=lambda candidate: candidate[0], reverse=True) + chosen_value, chosen_action, chosen_result, remaining = candidates[0] + alternative = next( + ( + candidate + for candidate in candidates[1:] + if _action_name(candidate[1]) != _action_name(chosen_action) + ), + candidates[1] if len(candidates) > 1 else candidates[0], + ) + dc, discharge, ac_rate, export_rate = chosen_action + start_soc = state_wh / bat.capacity_wh * 100 + state_wh = chosen_result["next_state_wh"] + plan.append( + TailPlanSlot( + slot=slot, + hour_from_start=control_horizon_hours + slot * context["slot_duration_h"], + action=_action_name(chosen_action), + alternative_action=_action_name(alternative[1]), + decision_margin_euro=max(chosen_value - alternative[0], 0.0), + soc_start_percentage=start_soc, + soc_end_percentage=state_wh / bat.capacity_wh * 100, + pv_wh=pv, + load_wh=load, + grid_import_wh=chosen_result["grid_import_wh"], + grid_export_wh=chosen_result["grid_export_wh"], + battery_charge_wh=chosen_result["battery_charge_wh"], + battery_discharge_wh=chosen_result["battery_discharge_wh"], + import_price_euro_per_kwh=price * 1000, + feed_in_tariff_euro_per_kwh=tariff * 1000, + slot_value_euro=chosen_result["reward_euro"], + remaining_value_euro=remaining, + ac_charge_factor=ac_rate, + dc_charge_allowed=dc, + discharge_allowed=discharge, + battery_grid_export_factor=export_rate, + ) + ) + return plan + + +def build_tail_value_curve( + *, + battery: Battery, + inverter: Inverter, + prices_euro_per_wh: np.ndarray, + load_wh: np.ndarray, + pv_wh: np.ndarray, + feed_in_euro_per_wh: np.ndarray, + continuation: TerminalValueCurve, + charge_rates: list[float], + export_rates: list[float], + direct_marketing: bool, + grid_points: int = 101, +) -> TailValueCurve: + """Solve the finite tail once, backwards in time, without mutating run devices. + + LCOS uses delivered DC energy, exactly as in GeneticSimulation. State grid + endpoints include battery minimum and maximum SOC. Continuous next states are + interpolated rather than rounded (which would invent or destroy energy). + """ + arrays = [ + np.asarray(a, dtype=float) + for a in (prices_euro_per_wh, load_wh, pv_wh, feed_in_euro_per_wh) + ] + if len({len(a) for a in arrays}) != 1 or any(not np.isfinite(a).all() for a in arrays): + raise ValueError("Tail forecasts must have equal lengths and contain only finite values") + bat = Battery(battery.parameters, prediction_hours=1, slot_duration_h=battery.slot_duration_h) + bat.charge_array = np.zeros(1, dtype=float) + inv = Inverter(inverter.parameters, battery=bat, slot_duration_h=battery.slot_duration_h) + states = np.linspace(bat.min_soc_wh, bat.max_soc_wh, grid_points) + usable = (states - bat.min_soc_wh) * bat.discharging_efficiency * inv.dc_to_ac_efficiency + continuation_values = np.array([continuation.value(e) for e in usable]) + operating_values = np.zeros(len(states)) + values = continuation_values.copy() + # (DC charge, local discharge, AC rate, export rate). Preserve production + # modes; direct marketing permits disabling DC charge to make headroom. + actions: list[TailAction] = [(1, 0, 0.0, 0.0), (1, 1, 0.0, 0.0)] + actions += [(1, 0, rate, 0.0) for rate in charge_rates if rate > 0] + if direct_marketing: + actions += [(0, 0, 0.0, 0.0), (0, 1, 0.0, 0.0)] + actions += [(0, 0, rate, 0.0) for rate in charge_rates if rate > 0] + actions += [(dc, 1, 0.0, rate) for dc in (0, 1) for rate in export_rates if rate > 0] + future_values: list[np.ndarray] = [np.empty(0)] * (len(arrays[0]) + 1) + future_values[-1] = continuation_values.copy() + for slot in reversed(range(len(arrays[0]))): + price, load, pv, tariff = (array[slot] for array in arrays) + best = np.full(len(states), -np.inf) + best_operating = np.zeros(len(states)) + best_continuation = np.zeros(len(states)) + for action in actions: + next_states = np.empty(len(states)) + rewards = np.empty(len(states)) + for i, energy in enumerate(states): + result = _simulate_action( + bat=bat, + inv=inv, + energy_wh=energy, + action=action, + price=price, + load=load, + pv=pv, + tariff=tariff, + direct_marketing=direct_marketing, + ) + rewards[i] = result["reward_euro"] + next_states[i] = result["next_state_wh"] + candidate_operating = rewards + np.interp(next_states, states, operating_values) + candidate_continuation = np.interp(next_states, states, continuation_values) + candidate = candidate_operating + candidate_continuation + better = candidate > best + best[better] = candidate[better] + best_operating[better] = candidate_operating[better] + best_continuation[better] = candidate_continuation[better] + values = best + operating_values = best_operating + continuation_values = best_continuation + future_values[slot] = best.copy() + result = TailValueCurve( + energy_wh=usable.tolist(), + value_euro=values.tolist(), + operating_value_euro=operating_values.tolist(), + continuation_value_euro=continuation_values.tolist(), + marginal_euro_per_kwh=(np.diff(values) / np.maximum(np.diff(usable), 1e-9) * 1000).tolist(), + window_slots=len(arrays[0]), + ) + result._trace_context = { + "battery_parameters": battery.parameters, + "inverter_parameters": inverter.parameters, + "slot_duration_h": battery.slot_duration_h, + "states": states, + "future_values": future_values, + "actions": actions, + "prices": arrays[0], + "load": arrays[1], + "pv": arrays[2], + "tariffs": arrays[3], + "direct_marketing": direct_marketing, + } + return result diff --git a/src/akkudoktoreos/optimization/genetic/terminalvalue.py b/src/akkudoktoreos/optimization/genetic/terminalvalue.py new file mode 100644 index 00000000..f8c45990 --- /dev/null +++ b/src/akkudoktoreos/optimization/genetic/terminalvalue.py @@ -0,0 +1,359 @@ +"""Terminal value of the energy left in the battery at the end of the horizon. + +The optimizer stops at the horizon, but the energy still stored in the battery +keeps its worth: it replaces grid imports that would otherwise be paid for +afterwards. Crediting that worth with a single price per kWh - the historical +``preis_euro_pro_wh_akku`` - cannot describe it, because the value of stored +energy is **not linear in the amount stored**: + +- The first kWh replaces the most expensive hour after the horizon. +- The next one replaces the second most expensive hour, and so on. +- Once every hour that PV cannot cover is served, further energy replaces + nothing; it is worth an export at best, and nothing at worst. + +The resulting value function is monotone and concave. A scalar has to pick one +slope: high enough for the first kWh means hoarding a full battery, low enough +for the last kWh means running it empty by midnight. This module builds the +curve instead. + +There is no forecast beyond the horizon, so the trailing window of the horizon +itself stands in for the day that follows: same season, same household rhythm, +same tariff structure. That approximation is the reason the curve is a planning +aid, not a prediction - which is also why the marginal values are deliberately +conservative wherever a choice exists. +""" + +from typing import Optional + +import numpy as np +from loguru import logger +from pydantic import Field + +from akkudoktoreos.core.pydantic import PydanticBaseModel + + +class TerminalValueCurve(PydanticBaseModel): + """Piecewise linear, concave value of battery energy left at the horizon. + + ``energy_wh`` and ``value_euro`` are the breakpoints of the cumulative + value, ``marginal_euro_per_kwh`` the slope of each segment. Both arrays + start at the origin; the curve is flat beyond its last breakpoint. + """ + + energy_wh: list[float] = Field( + default_factory=list, + json_schema_extra={ + "description": "Breakpoints of usable AC energy left in the battery [Wh]." + }, + ) + value_euro: list[float] = Field( + default_factory=list, + json_schema_extra={"description": "Cumulative credit at each breakpoint [EUR]."}, + ) + operating_value_euro: list[float] = Field( + default_factory=list, + json_schema_extra={ + "description": "Tail operating component at each breakpoint [EUR]; empty for a proxy curve." + }, + ) + continuation_value_euro: list[float] = Field( + default_factory=list, + json_schema_extra={ + "description": "Continuation component at each breakpoint [EUR]; empty for a proxy curve." + }, + ) + marginal_euro_per_kwh: list[float] = Field( + default_factory=list, + json_schema_extra={ + "description": ( + "Marginal value of the segment that starts at each breakpoint " + "[EUR/kWh]. May be negative or non-monotone in TAIL mode." + ) + }, + ) + residual_energy_wh: float = Field( + default=0.0, + json_schema_extra={ + "description": ( + "Energy up to which the curve is backed by residual load - the " + "knee. Everything beyond it is only worth an export." + ) + }, + ) + window_slots: int = Field( + default=0, + json_schema_extra={ + "description": ( + "Number of trailing horizon slots the curve was derived from. " + "Fewer slots than a full day mean a shorter proxy period." + ) + }, + ) + + def value(self, energy_wh: float) -> float: + """Return the credit for ``energy_wh`` of usable AC energy [EUR]. + + Args: + energy_wh: Usable AC energy left in the battery. + + Returns: + Interpolated value of the curve; 0.0 for an empty curve. + """ + if not self.energy_wh or energy_wh <= 0.0: + return 0.0 + return float(np.interp(energy_wh, self.energy_wh, self.value_euro)) + + +class TailDiagnostics(PydanticBaseModel): + """Forecast summary used by the deterministic tail optimization.""" + + slots: int = 0 + slot_hours: float = 0.0 + soc_grid_points: int = 0 + min_import_price_euro_per_kwh: float = 0.0 + max_import_price_euro_per_kwh: float = 0.0 + min_feed_in_tariff_euro_per_kwh: float = 0.0 + max_feed_in_tariff_euro_per_kwh: float = 0.0 + negative_import_price_slots: int = 0 + positive_battery_export_slots: int = 0 + + +class TailPlanSlot(PydanticBaseModel): + """One diagnostic slot of the optimal tail path. + + These values explain the lookahead used for fitness. They are diagnostics + only and are never copied into the executable control arrays. + """ + + slot: int + hour_from_start: float + action: str + alternative_action: str = "" + decision_margin_euro: float = 0.0 + soc_start_percentage: float + soc_end_percentage: float + pv_wh: float + load_wh: float + grid_import_wh: float + grid_export_wh: float + battery_charge_wh: float + battery_discharge_wh: float + import_price_euro_per_kwh: float + feed_in_tariff_euro_per_kwh: float + slot_value_euro: float + remaining_value_euro: float + ac_charge_factor: float + dc_charge_allowed: int + discharge_allowed: int + battery_grid_export_factor: float + + +class TerminalValueResult(PydanticBaseModel): + """What the optimizer credited for the energy left in the battery.""" + + control_horizon_hours: float = 0 + requested_tail_hours: float = 0 + effective_tail_hours: float = 0 + tail_end_hour: float = 0 + continuation_mode: str = "FIXED" + + mode: str = Field( + json_schema_extra={ + "description": "Terminal value mode the run used: TAIL, AUTO or FIXED.", + "examples": ["TAIL", "AUTO", "FIXED"], + } + ) + battery_energy_wh: float = Field( + default=0.0, + json_schema_extra={ + "description": "Usable AC energy left in the battery at the end of the horizon [Wh]." + }, + ) + credited_euro: float = Field( + default=0.0, + json_schema_extra={"description": "Credit applied to the total balance [EUR]."}, + ) + tail_operating_euro: float = Field( + default=0.0, + json_schema_extra={ + "description": ( + "Optimal net cash flow within the effective tail for the selected " + "control-end battery state [EUR]." + ) + }, + ) + continuation_value_euro: float = Field( + default=0.0, + json_schema_extra={ + "description": ( + "Continuation credit remaining at the end of the optimal tail path [EUR]." + ) + }, + ) + curve: Optional[TerminalValueCurve] = Field( + default=None, + json_schema_extra={ + "description": ( + "Combined tail value curve (tail operation plus continuation) read by fitness; " + "None in FIXED mode." + ) + }, + ) + continuation_curve: Optional[TerminalValueCurve] = Field( + default=None, + json_schema_extra={ + "description": "Conservative AUTO proxy constructed at the effective tail end." + }, + ) + tail_diagnostics: Optional[TailDiagnostics] = None + tail_plan: list[TailPlanSlot] = Field( + default_factory=list, + json_schema_extra={ + "description": ( + "Diagnostic optimal battery path inside the tail. It explains the " + "lookahead but is never an executable control plan." + ) + }, + ) + reason: str = Field( + default="", + json_schema_extra={ + "description": ( + "Why this mode applied. Empty in AUTO mode; in FIXED mode it " + "says whether FIXED was configured or whether AUTO fell back " + "because no curve could be derived." + ), + "examples": ["", "terminal_value_mode is FIXED"], + }, + ) + + +def build_terminal_value_curve( + *, + prices_euro_per_wh: np.ndarray, + load_wh: np.ndarray, + pv_wh: np.ndarray, + feed_in_euro_per_wh: np.ndarray, + max_energy_wh: float, + lcos_euro_per_kwh: float = 0.0, + dc_to_ac_efficiency: float = 1.0, + grid_export_allowed: bool = False, +) -> TerminalValueCurve: + """Build the terminal value curve from the trailing horizon window. + + Every slot of the window contributes its residual load - the part of the + load that PV does not cover - at its import price. Sorting those slots by + price and accumulating them yields the marginal value of the first, second, + ... kWh in the battery. Energy beyond the residual load can only be + exported, and only when direct marketing allows it. + + Args: + prices_euro_per_wh: Import prices of the window [EUR/Wh]. + load_wh: Load per slot of the window [Wh]. + pv_wh: PV generation per slot of the window [Wh]. + feed_in_euro_per_wh: Feed-in tariff of the window [EUR/Wh]. + max_energy_wh: Usable AC energy of a full battery [Wh]; the curve ends here. + lcos_euro_per_kwh: Levelized cost of storage, already charged per + delivered DC energy in the simulation and therefore subtracted here + so stored energy is not credited twice. + dc_to_ac_efficiency: Inverter efficiency, used to convert the LCOS from + delivered DC energy to the AC energy of the curve. + grid_export_allowed: Whether the battery may feed the grid (direct + marketing). Without it, energy beyond the residual load gets no + credit: it can neither be exported nor is its use covered by the + proxy window. + + Returns: + The curve; empty when the window carries no usable information. + """ + window = min(len(prices_euro_per_wh), len(load_wh), len(pv_wh)) + if window <= 0 or max_energy_wh <= 0.0: + return TerminalValueCurve() + + residual = np.maximum(load_wh[:window] - pv_wh[:window], 0.0) + prices = np.asarray(prices_euro_per_wh[:window], dtype=float) + + # LCOS is charged on delivered DC energy; the curve is in AC energy. + lcos_per_wh_ac = (lcos_euro_per_kwh / 1000.0) / max(dc_to_ac_efficiency, 1e-9) + + order = np.argsort(-prices) + energy_points: list[float] = [0.0] + value_points: list[float] = [0.0] + marginals: list[float] = [] + + cumulative_energy = 0.0 + cumulative_value = 0.0 + for index in order: + slot_energy = float(residual[index]) + if slot_energy <= 0.0: + continue + # Negative or very cheap hours are not worth storing energy for. + marginal = max(float(prices[index]) - lcos_per_wh_ac, 0.0) + if marginal <= 0.0: + continue + slot_energy = min(slot_energy, max_energy_wh - cumulative_energy) + if slot_energy <= 0.0: + break + cumulative_energy += slot_energy + cumulative_value += slot_energy * marginal + energy_points.append(cumulative_energy) + value_points.append(cumulative_value) + marginals.append(marginal * 1000.0) + + # Everything beyond the residual load can only be sold. A median feed-in + # tariff rather than the best one: exporting all of it in the single best + # slot is not something the horizon can promise. + residual_energy_wh = cumulative_energy + if grid_export_allowed and cumulative_energy < max_energy_wh: + positive_feed_in = [ + float(value) for value in feed_in_euro_per_wh[:window] if float(value) > 0.0 + ] + export_marginal = max( + (float(np.median(positive_feed_in)) if positive_feed_in else 0.0) - lcos_per_wh_ac, + 0.0, + ) + if export_marginal > 0.0: + remaining = max_energy_wh - cumulative_energy + cumulative_energy += remaining + cumulative_value += remaining * export_marginal + energy_points.append(cumulative_energy) + value_points.append(cumulative_value) + marginals.append(export_marginal * 1000.0) + + if len(energy_points) <= 1: + logger.debug("Terminal value curve is empty - no priced residual load in the window.") + return TerminalValueCurve(window_slots=window) + + # The segment slopes are decreasing by construction (prices were sorted), + # so the curve is concave; the export tail is the flattest segment. + return TerminalValueCurve( + energy_wh=energy_points, + value_euro=value_points, + marginal_euro_per_kwh=marginals, + residual_energy_wh=residual_energy_wh, + window_slots=window, + ) + + +def trailing_window( + values: Optional[np.ndarray], + end_slot: int, + window_slots: int, +) -> np.ndarray: + """Return the ``window_slots`` values in front of ``end_slot``. + + Args: + values: Full slot array, or None. + end_slot: Exclusive end of the window (end of the optimization horizon). + window_slots: Desired window length; a shorter horizon yields less. + + Returns: + The window as a float array, empty when no data is available. + """ + if values is None: + return np.zeros(0, dtype=float) + end = min(int(end_slot), len(values)) + start = max(end - int(window_slots), 0) + if end <= start: + return np.zeros(0, dtype=float) + return np.asarray(values[start:end], dtype=float) diff --git a/tests/test_tailvalue_physics.py b/tests/test_tailvalue_physics.py new file mode 100644 index 00000000..0d7a222a --- /dev/null +++ b/tests/test_tailvalue_physics.py @@ -0,0 +1,200 @@ +"""Economic tail scenarios and hard control/forecast boundaries.""" + +from unittest.mock import patch + +import numpy as np +import pandas as pd +import pytest + +from akkudoktoreos.config.config import SettingsEOSDefaults +from akkudoktoreos.core.coreabc import get_ems +from akkudoktoreos.devices.genetic.battery import Battery +from akkudoktoreos.devices.genetic.inverter import Inverter +from akkudoktoreos.optimization.genetic.forecast import bounded_forecast_array +from akkudoktoreos.optimization.genetic.genetic import GeneticOptimization +from akkudoktoreos.devices.genetic.inverter import InverterParameters +from akkudoktoreos.devices.genetic.battery import SolarPanelBatteryParameters +from akkudoktoreos.optimization.genetic.geneticparams import GeneticOptimizationParameters +from akkudoktoreos.optimization.genetic.tailvalue import build_tail_value_curve +from akkudoktoreos.optimization.genetic.terminalvalue import TerminalValueCurve +from akkudoktoreos.utils.datetimeutil import to_datetime, to_duration + + +def devices(power=1000, efficiency=1.0, lcos=0, ac_limit=None, export_power=5000): + bat = Battery( + SolarPanelBatteryParameters( + device_id="battery1", + capacity_wh=1000, + max_charge_power_w=power, + charging_efficiency=efficiency, + discharging_efficiency=efficiency, + initial_soc_percentage=50, + levelized_cost_of_storage_kwh=lcos, + charge_rates=[0, 0.5, 1], + ), + prediction_hours=1, + ) + inv = Inverter( + InverterParameters( + device_id="inverter1", + battery_id="battery1", + max_power_wh=export_power, + dc_to_ac_efficiency=1, + ac_to_dc_efficiency=1, + max_ac_charge_power_w=ac_limit, + ), + battery=bat, + ) + return bat, inv + + +def curve( + prices=(-0.1, 0.3), + tariffs=(0, 0.3), + direct=True, + continuation=None, + load=None, + pv=None, + **kwargs, +): + bat, inv = devices(**kwargs) + return build_tail_value_curve( + battery=bat, + inverter=inv, + prices_euro_per_wh=np.array(prices) / 1000, + feed_in_euro_per_wh=np.array(tariffs) / 1000, + load_wh=np.zeros(len(prices)) if load is None else np.array(load), + pv_wh=np.zeros(len(prices)) if pv is None else np.array(pv), + continuation=continuation or TerminalValueCurve(), + charge_rates=[0.5, 1], + export_rates=[1], + direct_marketing=direct, + ) + + +def test_headroom_has_value_and_empty_state_can_earn(): + c = curve() + assert c.value(0) == pytest.approx(0.4) + tail, continuation = c.component_values(0) + assert tail == pytest.approx(0.4) + assert continuation == pytest.approx(0.0) + assert c.value(0) == pytest.approx(tail + continuation) + assert c.value(500) > c.value(1000) + assert any(v < 0 for v in c.marginal_euro_per_kwh) + + +def test_chronology_changes_arbitrage(): + forward = curve() + reverse = curve(prices=(0.3, -0.1), tariffs=(0.3, 0)) + assert forward.value(0) > reverse.value(0) + + +def test_discharge_and_ac_power_limits(): + limited = curve(prices=(1,), tariffs=(1,), power=100) + assert limited.value(1000) == pytest.approx(0.1) + limited_ac = curve(ac_limit=100) + assert limited_ac.value(0) == pytest.approx(0.04) + limited_inverter = curve(prices=(1,), tariffs=(1,), export_power=50) + assert limited_inverter.value(1000) == pytest.approx(0.05) + + +def test_losses_and_lcos_reduce_arbitrage(): + ideal = curve(prices=(0.1, 0.3)) + lossy = curve(prices=(0.1, 0.3), efficiency=0.8) + assert 0 < lossy.value(0) < ideal.value(0) + assert curve(prices=(0.1, 0.3), lcos=0.25).value(0) == pytest.approx(0) + + +def test_no_battery_export_without_permission(): + assert curve(prices=(0.1, 0.3), direct=False).value(0) == pytest.approx(0) + + +def test_pv_surplus_can_be_stored_for_local_load(): + c = curve(prices=(0.2, 0.3), tariffs=(0, 0), pv=[1000, 0], load=[0, 1000], direct=False) + assert c.value(0) == pytest.approx(0) + # Without PV the same empty battery must buy energy to serve the load. + assert curve(prices=(0.2, 0.3), tariffs=(0, 0), load=[0, 1000], direct=False).value( + 0 + ) < c.value(0) + + +def test_continuation_survives_tail_end(): + continuation = TerminalValueCurve(energy_wh=[0, 1000], value_euro=[0, 0.2]) + c = curve(prices=(0.5,), tariffs=(0,), continuation=continuation) + assert c.value(1000) == pytest.approx(0.2) + tail, continuation_credit = c.component_values(1000) + assert tail == pytest.approx(0.0) + assert continuation_credit == pytest.approx(0.2) + + +def test_tail_diagnostic_plan_explains_the_selected_path(): + c = curve() + plan = c.diagnostic_plan(0, control_horizon_hours=24) + assert len(plan) == 2 + assert plan[0].hour_from_start == 24 + assert plan[0].action == "GRID_CHARGE" + assert plan[0].soc_end_percentage > plan[0].soc_start_percentage + assert plan[0].grid_import_wh > 0 + assert plan[1].action == "BATTERY_EXPORT" + assert plan[1].soc_end_percentage < plan[1].soc_start_percentage + assert plan[1].grid_export_wh > 0 + assert sum(slot.slot_value_euro for slot in plan) == pytest.approx(c.value(0)) + + + + +def test_provider_values_are_not_extrapolated(): + from types import SimpleNamespace + + start = to_datetime("2026-09-05T00:00:00Z") + series = pd.Series([1.0, 2.0], index=pd.date_range(start=start, periods=2, freq="h")) + provider = SimpleNamespace(key_to_series=lambda *a, **kw: series) + result = bounded_forecast_array( + provider, + key="price", + start_datetime=start, + end_datetime=start.add(hours=3), + interval=to_duration("15 minutes"), + ) + assert result[:8].tolist() == [1.0] * 4 + [2.0] * 4 + assert np.isnan(result[8:]).all() + + + +def test_disabled_ac_conversion_cannot_earn_negative_price_revenue(): + bat, inv = devices() + inv.parameters.ac_to_dc_efficiency = 0 + c = build_tail_value_curve( + battery=bat, + inverter=inv, + prices_euro_per_wh=np.array([-0.001, 0.001]), + feed_in_euro_per_wh=np.array([0.0, 0.001]), + load_wh=np.zeros(2), + pv_wh=np.zeros(2), + continuation=TerminalValueCurve(), + charge_rates=[1], + export_rates=[1], + direct_marketing=True, + ) + assert c.value(0) == pytest.approx(0) + assert bat.soc_wh == 500 # Building the tail never mutates the real battery. + + + +def test_missing_provider_key_stays_missing(): + from types import SimpleNamespace + + def unavailable(*a, **kw): + raise KeyError("price unavailable") + + start = to_datetime("2026-09-05T00:00:00Z") + result = bounded_forecast_array( + SimpleNamespace(key_to_series=unavailable), + key="price", + start_datetime=start, + end_datetime=start.add(hours=2), + interval=to_duration("1 hour"), + ) + assert np.isnan(result).all() + assert len(result) == 2 + diff --git a/tests/test_terminalvalue.py b/tests/test_terminalvalue.py new file mode 100644 index 00000000..cdcce97f --- /dev/null +++ b/tests/test_terminalvalue.py @@ -0,0 +1,127 @@ +"""Tests for the concave terminal value of the energy left in the battery.""" + +import numpy as np +import pytest + +from akkudoktoreos.optimization.genetic.terminalvalue import ( + build_terminal_value_curve, + trailing_window, +) + + +def _curve(**overrides): + """Two expensive slots, one cheap one, no PV, 10 kWh of usable battery.""" + params = dict( + prices_euro_per_wh=np.array([0.0004, 0.0003, 0.0001]), + load_wh=np.array([1000.0, 1000.0, 1000.0]), + pv_wh=np.array([0.0, 0.0, 0.0]), + feed_in_euro_per_wh=np.array([0.00008, 0.00008, 0.00008]), + max_energy_wh=10000.0, + lcos_euro_per_kwh=0.0, + dc_to_ac_efficiency=1.0, + grid_export_allowed=False, + ) + params.update(overrides) + return build_terminal_value_curve(**params) + + +def test_marginal_value_follows_the_most_expensive_hours_first(): + """The first stored kWh replaces the most expensive slot, then the next.""" + curve = _curve() + + # 0.40, 0.30 and 0.10 EUR/kWh, in that order. + assert curve.marginal_euro_per_kwh == pytest.approx([0.4, 0.3, 0.1]) + assert curve.energy_wh == pytest.approx([0.0, 1000.0, 2000.0, 3000.0]) + assert curve.value_euro == pytest.approx([0.0, 0.4, 0.7, 0.8]) + + +def test_curve_is_concave_and_saturates(): + """Marginal values only decrease, and beyond the last breakpoint nothing is added.""" + curve = _curve() + marginals = curve.marginal_euro_per_kwh + + assert all(a >= b for a, b in zip(marginals, marginals[1:])) + # The residual load of the window is 3 kWh - more energy replaces nothing. + assert curve.value(3000.0) == pytest.approx(0.8) + assert curve.value(9000.0) == pytest.approx(0.8) + + +def test_value_interpolates_within_a_segment(): + """Half of the first slot is worth half of the first segment.""" + curve = _curve() + assert curve.value(500.0) == pytest.approx(0.2) + + +def test_pv_reduces_the_residual_load(): + """Only load that PV cannot cover can be replaced by stored energy.""" + curve = _curve(pv_wh=np.array([600.0, 1000.0, 0.0])) + + # Slot 0 keeps 400 Wh, slot 1 is fully covered by PV, slot 2 keeps 1000 Wh. + assert curve.energy_wh == pytest.approx([0.0, 400.0, 1400.0]) + assert curve.marginal_euro_per_kwh == pytest.approx([0.4, 0.1]) + + +def test_lcos_is_subtracted_from_the_marginal_value(): + """Storage cost is already charged on discharge and must not be credited twice.""" + curve = _curve(lcos_euro_per_kwh=0.05, dc_to_ac_efficiency=1.0) + assert curve.marginal_euro_per_kwh == pytest.approx([0.35, 0.25, 0.05]) + + +def test_negative_prices_do_not_create_value(): + """Storing energy for an hour that pays nothing is not worth anything.""" + curve = _curve(prices_euro_per_wh=np.array([0.0004, -0.0001, 0.0])) + assert curve.marginal_euro_per_kwh == pytest.approx([0.4]) + assert curve.value(5000.0) == pytest.approx(0.4) + + +def test_export_tail_only_with_direct_marketing(): + """Surplus beyond the residual load is worth an export - if export is allowed.""" + without = _curve(grid_export_allowed=False) + with_export = _curve(grid_export_allowed=True) + + assert without.value(10000.0) == pytest.approx(0.8) + # 7 kWh beyond the residual load at the median feed-in tariff of 0.08 EUR/kWh. + assert with_export.value(10000.0) == pytest.approx(0.8 + 7.0 * 0.08) + assert with_export.marginal_euro_per_kwh[-1] == pytest.approx(0.08) + + +def test_residual_energy_marks_the_knee(): + """The knee separates load-backed value from the export tail.""" + without = _curve(grid_export_allowed=False) + with_export = _curve(grid_export_allowed=True) + + # 3 kWh of residual load in the window, whether or not export is allowed. + assert without.residual_energy_wh == pytest.approx(3000.0) + assert with_export.residual_energy_wh == pytest.approx(3000.0) + # Only the export tail reaches beyond it. + assert without.energy_wh[-1] == pytest.approx(3000.0) + assert with_export.energy_wh[-1] == pytest.approx(10000.0) + + +def test_curve_is_capped_by_the_usable_battery_energy(): + """A battery smaller than the residual load ends the curve early.""" + curve = _curve(max_energy_wh=1500.0) + assert curve.energy_wh[-1] == pytest.approx(1500.0) + assert curve.value(5000.0) == pytest.approx(0.4 + 0.5 * 0.3) + + +def test_empty_window_yields_an_empty_curve(): + """Without data there is no curve, and no credit.""" + curve = build_terminal_value_curve( + prices_euro_per_wh=np.zeros(0), + load_wh=np.zeros(0), + pv_wh=np.zeros(0), + feed_in_euro_per_wh=np.zeros(0), + max_energy_wh=10000.0, + ) + assert curve.energy_wh == [] + assert curve.value(5000.0) == 0.0 + + +def test_trailing_window_takes_the_end_of_the_horizon(): + values = np.arange(10, dtype=float) + + assert list(trailing_window(values, end_slot=8, window_slots=3)) == [5.0, 6.0, 7.0] + # A window longer than the horizon yields what there is. + assert list(trailing_window(values, end_slot=2, window_slots=5)) == [0.0, 1.0] + assert list(trailing_window(None, end_slot=8, window_slots=3)) == []