feat: configurable timestep selection for the run pipeline (#663)

edisgo/io/powermodels_io.py

@@ -370,7 +370,7 @@ def from_powermodels(
         # calculate relative error
         df2 = deepcopy(df)
         for flex in df2.columns:
-            if type(hv_flex_dict[flex]) == pd.Series:
+            if isinstance(hv_flex_dict[flex], pd.Series):
                 abs_error = abs(df2[flex].values - hv_flex_dict[flex].values)
                 rel_error = [
                     abs_error[i] / hv_flex_dict[flex].iloc[i]
@@ -379,10 +379,15 @@ def from_powermodels(
                     for i in range(len(abs_error))
                 ]
             else:
-                abs_error = abs(df2[flex].values - hv_flex_dict[flex].sum(axis=1).values)
+                abs_error = abs(
+                    df2[flex].values - hv_flex_dict[flex].sum(axis=1).values
+                )
                 rel_error = [
                     abs_error[i] / hv_flex_dict[flex].sum(axis=1).iloc[i]
-                    if ((abs_error > 0.01)[i] & (hv_flex_dict[flex].sum(axis=1).iloc[i] != 0))
+                    if (
+                        (abs_error > 0.01)[i]
+                        & (hv_flex_dict[flex].sum(axis=1).iloc[i] != 0)
+                    )
                     else 0
                     for i in range(len(abs_error))
                 ]
@@ -1061,6 +1066,18 @@ def _build_battery_storage(
         * edisgo_obj.topology.storage_units_df.max_hours
     )
 
+    # The end-of-period SoC step (timeindex[-1] + freq) is only used as the OPF
+    # boundary (soc_end) and is not an optimized time step. When the time index
+    # is a reduced, non-contiguous selection, that step can fall in a gap and be
+    # missing from the source SoC series (which only carried a trailing step for
+    # the very last interval), leaving it NaN. A NaN boundary makes the Julia OPF
+    # fail with "Inf - Inf". Forward-fill (then back-fill) so the boundary takes
+    # the interval's last valid SoC — a harmless approximation for a throwaway
+    # scaffolding step.
+    edisgo_obj.overlying_grid.storage_units_soc = (
+        edisgo_obj.overlying_grid.storage_units_soc.ffill().bfill()
+    )
+
     for stor_i in np.arange(len(flexible_storage_units)):
         idx_bus = _mapping(
             psa_net,
@@ -1357,6 +1374,12 @@ def _build_heat_storage(psa_net, pm, edisgo_obj, s_base, flexible_hps, opf_versi
     edisgo_obj.overlying_grid.heat_storage_units_soc = pd.concat(
         [df_decentral, df_central], axis=1
     )
+    # Fill the end-of-period boundary SoC step (see storage note above) so a
+    # reduced, non-contiguous time index does not leave a NaN boundary that
+    # breaks the Julia OPF.
+    edisgo_obj.overlying_grid.heat_storage_units_soc = (
+        edisgo_obj.overlying_grid.heat_storage_units_soc.ffill().bfill()
+    )
 
     heat_storage_df = heat_storage_df.loc[flexible_hps]
     for stor_i in np.arange(len(flexible_hps)):
@@ -1616,7 +1639,7 @@ def _build_hv_requirements(
     )
 
     for i in np.arange(len(opf_flex)):
-        if type(hv_flex_dict[opf_flex[i]]) == pd.DataFrame:
+        if isinstance(hv_flex_dict[opf_flex[i]], pd.DataFrame):
             pm["HV_requirements"][str(i + 1)] = {
                 "P": hv_flex_dict[opf_flex[i]].sum(axis=1).iloc[0],
                 "name": opf_flex[i],
@@ -1957,7 +1980,7 @@ def _build_component_timeseries(
 
     if (kind == "HV_requirements") & (pm["opf_version"] in [3, 4]):
         for i in np.arange(len(opf_flex)):
-            if type(hv_flex_dict[opf_flex[i]])==pd.DataFrame:
+            if isinstance(hv_flex_dict[opf_flex[i]], pd.DataFrame):
                 pm_comp[(str(i + 1))] = {
                     "P": hv_flex_dict[opf_flex[i]].sum(axis=1).round(20).tolist(),
                 }

edisgo/opf/powermodels_opf.py

@@ -9,20 +9,136 @@
 #
 # SPDX-License-Identifier: AGPL-3.0-or-later
 
+import copy
 import json
 import logging
 import os
 import subprocess
 import sys
 
 import numpy as np
+import pandas as pd
 
 from edisgo.flex_opt import exceptions
 from edisgo.io.powermodels_io import from_powermodels
 from edisgo.network.topology import Topology
 
 logger = logging.getLogger(__name__)
 
+# Time-indexed opf_results attributes that from_powermodels overwrites on each
+# call. When the OPF is run separately per interval, these must be concatenated
+# across intervals so opf_results covers the full (reduced) time index. Nested
+# containers (LineVariables etc.) are listed via their sub-frame attribute names.
+_OPF_FLAT_TIME_FRAMES = (
+    "slack_generator_t",
+    "hv_requirement_slacks_t",
+)
+_OPF_NESTED_TIME_FRAMES = (
+    "lines_t",
+    "heat_storage_t",
+    "grid_slacks_t",
+    "battery_storage_t",
+)
+
+
+def _with_freq(index):
+    """Return the DatetimeIndex with its frequency inferred/attached if regular.
+
+    Reducing/uniting time indices drops the ``freq`` attribute; several
+    downstream consumers (notably the powermodels OPF) do
+    ``timeindex[-1] + timeindex.freq`` and break on ``freq is None``. This
+    re-attaches the freq when the index is regularly spaced (a no-op otherwise).
+    """
+    if index.freq is not None or len(index) < 2:
+        return index
+    inferred = pd.infer_freq(index)
+    if inferred is not None:
+        try:
+            return pd.DatetimeIndex(index, freq=inferred)
+        except (ValueError, TypeError):
+            return index
+    return index
+
+
+def _contiguous_intervals(timeindex):
+    """
+    Split a time index into contiguous intervals.
+
+    Automatic timestep selection can reduce the time index to disconnected
+    intervals (e.g. one load-case and one feed-in-case week). This helper detects
+    the gap(s) so the OPF can be run separately per interval — storage/heat state
+    does not carry across a gap, so a single OPF over the concatenated steps would
+    be wrong.
+
+    A boundary is placed wherever the spacing between two consecutive time steps
+    exceeds the regular step (the smallest spacing in the index). A contiguous
+    index therefore yields a single interval. Each returned interval has its
+    ``freq`` restored (set operations that produced the reduced index drop it).
+
+    Parameters
+    ----------
+    timeindex : pandas.DatetimeIndex
+
+    Returns
+    -------
+    list of pandas.DatetimeIndex
+        One entry per contiguous interval, in chronological order. Empty index
+        in -> empty list out; a single time step -> one interval.
+    """
+    timeindex = timeindex.sort_values()
+    if len(timeindex) <= 1:
+        return [timeindex] if len(timeindex) else []
+    diffs = timeindex[1:] - timeindex[:-1]
+    step = diffs.min()
+    breaks = [i + 1 for i, d in enumerate(diffs) if d > step]
+    starts = [0] + breaks
+    ends = breaks + [len(timeindex)]
+    return [_with_freq(timeindex[s:e]) for s, e in zip(starts, ends)]
+
+
+def _snapshot_opf_time_frames(opf_results):
+    """Copy the time-indexed opf_results frames produced by one interval's OPF."""
+    snap = {}
+    for attr in _OPF_FLAT_TIME_FRAMES:
+        snap[attr] = getattr(opf_results, attr).copy()
+    for attr in _OPF_NESTED_TIME_FRAMES:
+        container = getattr(opf_results, attr)
+        snap[attr] = {
+            sub: getattr(container, sub).copy() for sub in container._attributes()
+        }
+    return snap
+
+
+def _merge_opf_time_frames(opf_results, snapshots):
+    """
+    Concatenate per-interval opf_results snapshots by time index and write them
+    back onto ``opf_results``, so its detailed frames cover the full reduced
+    index rather than only the last interval's.
+    """
+
+    def _concat(frames):
+        frames = [f for f in frames if f is not None and not f.empty]
+        if not frames:
+            return pd.DataFrame()
+        return pd.concat(frames).sort_index()
+
+    for attr in _OPF_FLAT_TIME_FRAMES:
+        setattr(opf_results, attr, _concat([s[attr] for s in snapshots]))
+    for attr in _OPF_NESTED_TIME_FRAMES:
+        container = getattr(opf_results, attr)
+        for sub in container._attributes():
+            setattr(container, sub, _concat([s[attr][sub] for s in snapshots]))
+
+    # Recompute the overlying_grid summary (opf_version 3/4) from the merged HV
+    # requirement slacks, since it is a reduction over the whole time index.
+    hv = opf_results.hv_requirement_slacks_t
+    if not hv.empty:
+        opf_results.overlying_grid = pd.DataFrame(
+            columns=["Highest error", "Mean error", "Sum error"],
+            index=hv.columns,
+            data=pd.concat([hv.max(), hv.mean(), hv.sum()], axis=1).values,
+        )
+
 
 def pm_optimize(
     edisgo_obj,
@@ -37,8 +153,162 @@ def pm_optimize(
     silence_moi: bool = False,
 ) -> None:
     """
-    Run OPF for edisgo object in julia subprocess and write results of OPF to edisgo
-    object. Results of OPF are time series of operation schedules of flexibilities.
+    Run OPF for the edisgo object and write results back to it.
+
+    If the time index is a single contiguous interval, this runs one OPF
+    (:func:`_pm_optimize_single`). If the time index is NON-contiguous
+    (disconnected intervals, e.g. from automatic timestep selection), each
+    contiguous interval is optimized separately and independently — storage/heat
+    state does not carry across the gap — and the results are combined:
+
+    * per-interval operation schedules accumulate in ``edisgo.timeseries``;
+    * the detailed ``edisgo.opf_results`` frames are merged by time index;
+    * a per-interval solve report is stored in
+      ``edisgo.opf_results.interval_results``;
+    * if any interval was infeasible, an
+      :class:`~.flex_opt.exceptions.InfeasibleModelError` is raised after the
+      feasible intervals' results have been stored.
+
+    The overlying-grid SOC attributes and reactive-power time series (which
+    ``to_powermodels`` / ``from_powermodels`` mutate or replace on the current
+    interval) are snapshotted and restored pristine before each interval so a
+    later interval sees intact input. Parameters are as for
+    :func:`_pm_optimize_single`.
+    """
+    opf_kwargs = dict(
+        s_base=s_base,
+        flexible_cps=flexible_cps,
+        flexible_hps=flexible_hps,
+        flexible_loads=flexible_loads,
+        flexible_storage_units=flexible_storage_units,
+        opf_version=opf_version,
+        method=method,
+        warm_start=warm_start,
+        silence_moi=silence_moi,
+    )
+
+    intervals = _contiguous_intervals(edisgo_obj.timeseries.timeindex)
+    if len(intervals) <= 1:
+        # single contiguous optimization. Re-set the (freq-restored) interval so
+        # the OPF sees a time index with a frequency — set operations upstream
+        # (e.g. timestep selection) drop it, and the OPF needs timeindex.freq.
+        if intervals:
+            edisgo_obj.set_timeindex(intervals[0])
+        _pm_optimize_single(edisgo_obj, **opf_kwargs)
+        return
+
+    logger.info(
+        f"pm_optimize: time index has {len(intervals)} disconnected intervals; "
+        f"running a separate OPF per interval."
+    )
+    full_timeindex = edisgo_obj.timeseries.timeindex
+
+    # Snapshot the shared input state that per-interval OPF runs mutate:
+    #  * overlying-grid SOC attributes are rewritten in place by to_powermodels;
+    #  * the reactive-power time series are fully REPLACED (not .loc-updated) by
+    #    the set_time_series_reactive_power_control() call inside from_powermodels.
+    # Reactive power was set on the full reduced index before this call; restore
+    # this input pristine before each interval. Active-power frames are NOT
+    # restored — they accumulate each interval's OPF results via .loc.
+    og = edisgo_obj.overlying_grid
+    og_snapshot = {attr: copy.deepcopy(getattr(og, attr)) for attr in og._attributes}
+    reactive_attrs = [
+        "_generators_reactive_power",
+        "_loads_reactive_power",
+        "_storage_units_reactive_power",
+    ]
+    reactive_snapshot = {
+        attr: copy.deepcopy(getattr(edisgo_obj.timeseries, attr, None))
+        for attr in reactive_attrs
+    }
+
+    def _restore_pristine_inputs():
+        for attr, value in og_snapshot.items():
+            setattr(og, attr, copy.deepcopy(value))
+        for attr, value in reactive_snapshot.items():
+            if value is not None:
+                setattr(edisgo_obj.timeseries, attr, copy.deepcopy(value))
+
+    # Pre-allocate the storage active-power schedule over the FULL reduced index
+    # so from_powermodels .loc-accumulates each interval's storage result instead
+    # of replacing the frame with an interval-only one (which would drop earlier
+    # intervals' storage schedules).
+    su_names = edisgo_obj.topology.storage_units_df.index
+    if len(su_names) > 0 and edisgo_obj.timeseries.storage_units_active_power.empty:
+        edisgo_obj.timeseries.storage_units_active_power = pd.DataFrame(
+            0.0, index=full_timeindex, columns=su_names
+        )
+
+    snapshots = []
+    report = []
+    try:
+        for interval in intervals:
+            _restore_pristine_inputs()
+            edisgo_obj.set_timeindex(interval)
+            entry = {
+                "start": interval[0],
+                "end": interval[-1],
+                "status": None,
+                "solver": None,
+                "solution_time": None,
+            }
+            try:
+                _pm_optimize_single(edisgo_obj, **opf_kwargs)
+                entry["status"] = edisgo_obj.opf_results.status
+                entry["solver"] = edisgo_obj.opf_results.solver
+                entry["solution_time"] = edisgo_obj.opf_results.solution_time
+                snapshots.append(_snapshot_opf_time_frames(edisgo_obj.opf_results))
+            except exceptions.InfeasibleModelError:
+                entry["status"] = "infeasible"
+                logger.warning(
+                    f"pm_optimize: OPF infeasible for interval "
+                    f"{interval[0]}..{interval[-1]}."
+                )
+            report.append(entry)
+    finally:
+        # restore the full (reduced) index so all intervals' schedules are exposed
+        # and undo the per-interval mutations of the overlying-grid/reactive input.
+        _restore_pristine_inputs()
+        edisgo_obj.set_timeindex(full_timeindex)
+
+    _merge_opf_time_frames(edisgo_obj.opf_results, snapshots)
+    edisgo_obj.opf_results.interval_results = report
+    solution_times = [
+        e["solution_time"] for e in report if e["solution_time"] is not None
+    ]
+    edisgo_obj.opf_results.solution_time = (
+        sum(solution_times) if solution_times else None
+    )
+    statuses = [e["status"] for e in report]
+    infeasible = [e for e in report if e["status"] == "infeasible"]
+    edisgo_obj.opf_results.status = (
+        "infeasible" if infeasible else (statuses[0] if statuses else None)
+    )
+    if infeasible:
+        raise exceptions.InfeasibleModelError(
+            f"OPF infeasible for {len(infeasible)} of {len(intervals)} time "
+            f"intervals; see edisgo.opf_results.interval_results. Results for "
+            f"feasible intervals have been stored."
+        )
+
+
+def _pm_optimize_single(
+    edisgo_obj,
+    s_base: int = 1,
+    flexible_cps: np.ndarray | None = None,
+    flexible_hps: np.ndarray | None = None,
+    flexible_loads: np.ndarray | None = None,
+    flexible_storage_units: np.ndarray | None = None,
+    opf_version: int = 1,
+    method: str = "soc",
+    warm_start: bool = False,
+    silence_moi: bool = False,
+) -> None:
+    """
+    Run a single-interval OPF for the edisgo object in a julia subprocess and
+    write results back to the edisgo object. Assumes the time index is a single
+    contiguous interval; :func:`pm_optimize` is the public entry point that
+    handles non-contiguous indices by calling this per interval.
 
     Parameters
     ----------