feat: add reduce_to_radial_network for looped systems

docs/new_features.ipynb

@@ -11,7 +11,56 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Version 2.3.7\n",
+    "## Version 2.3.8\n",
+    "\n",
+    "### New `reduce_to_radial_network` utility\n",
+    "\n",
+    "New function `reduce_to_radial_network(dist_system, name)` in `gdm.distribution.model_reduction.reducer` (also exported from `gdm.distribution`) that reduces a looped distribution system to a radial network:\n",
+    "\n",
+    "- Every electrical loop is broken deterministically: if the loop contains a closed switch, that switch is opened; otherwise the shortest line on the loop (ties broken by component name) is replaced with an open `MatrixImpedanceSwitch` between the same buses.\n",
+    "- The input system is not modified — a reduced deep copy is returned under the given name.\n",
+    "- Geometry-based lines are converted to their matrix representation at 60 Hz when selected for conversion.\n",
+    "- Only cycle edges are touched, so the reduced network stays connected and `get_directed_graph` prunes no edges or emits no warnings.\n",
+    "\n",
+    "### Multiple switches in a loop\n",
+    "\n",
+    "Fixed cycle breaking in `DistributionSystem.get_directed_graph` when producing a radial network and a loop contains more than one switch edge:\n",
+    "\n",
+    "- `find_switch_buses_in_cycle` now returns `(bus_1, bus_2)` tuples for each switch edge found in the cycle instead of a flat list of bus names.\n",
+    "- The first switch edge is removed with explicit direction checking (`has_edge` in both directions); cycles whose switch edge is not present in the DFS tree are skipped with a warning instead of raising.\n",
+    "\n",
+    "### Faster power DataFrame assembly\n",
+    "\n",
+    "Refactored the internals of `get_combined_load_time_series_df` and `get_combined_solar_time_series_df`:\n",
+    "\n",
+    "- Timestamps for `SingleTimeSeries` are generated with `pd.date_range` instead of a Python list comprehension.\n",
+    "- Each time series slice is built as a dict of ndarrays (`value`, `timestamp`) and scalars; the final DataFrame is assembled in one pass with `np.concatenate`/`np.repeat`, avoiding per-slice DataFrames and repeated `pd.concat` calls.\n",
+    "- The check that `per_phase_function` is required when `aggregate_phases=False` now runs up front.\n",
+    "\n",
+    "### MCP SDK 2.0 support\n",
+    "\n",
+    "The MCP server was updated to the new `mcp>=2.x` SDK API:\n",
+    "\n",
+    "- The `Server` instance is created with `on_list_tools`/`on_call_tool` callbacks and an explicit version instead of using decorators.\n",
+    "- Tool handlers now receive a `ServerRequestContext` and return `ListToolsResult` / `CallToolResult`, flagging errors via the `is_error` field.\n",
+    "\n",
+    "### Upgrade handler\n",
+    "\n",
+    "- Added the `from__2_3_7__to__2_3_8` upgrade handler, which runs component metadata migration on existing data.\n",
+    "- Fixed the 2.3.2 to 2.3.3 handler so metadata migration only runs when a component actually needs it (fixes an upgrade failure).\n",
+    "\n",
+    "### Dependency updates\n",
+    "\n",
+    "- Updated `pandas` requirement from `~=3.0.3` to `~=3.0.5`.\n",
+    "- Updated the optional `mcp` requirement from `>=1.28.1` to `>=2.1.0`.\n",
+    "- Pinned `infrasys` to `==1.2.1`."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+     "## Version 2.3.7\n",
     "\n",
     "### Parallel single-phase transformer aggregation: 2-unit groups\n",
     "\n",

pyproject.toml

@@ -24,7 +24,7 @@ dependencies = [
   "semver",
   "networkx", 
   "pydantic",
-  "infrasys~=1.2", 
+  "infrasys==1.2.1", 
   "importlib_metadata",
   "typer",
   "pandas~=3.0.5",

src/gdm/distribution/__init__.py

@@ -5,6 +5,7 @@
 # from gdm.tracked_changes import apply_tracked_changes, apply_update_scenario, get_distribution_system_on_date
 from gdm.distribution.model_reduction.reducer import (
     reduce_to_primary_system,
+    reduce_to_radial_network,
     reduce_to_three_phase_system,
 )
 from gdm.distribution.distribution_graph import build_graph_from_system

src/gdm/distribution/model_reduction/__init__.py

@@ -1,2 +1,3 @@
 from gdm.distribution.model_reduction.reducer import reduce_to_primary_system
+from gdm.distribution.model_reduction.reducer import reduce_to_radial_network
 from gdm.distribution.model_reduction.reducer import reduce_to_three_phase_system

src/gdm/distribution/model_reduction/reducer.py

@@ -4,14 +4,22 @@
 from infrasys.time_series_models import SingleTimeSeries, TimeSeriesData
 import networkx as nx
 
+from gdm.distribution.components.base.distribution_branch_base import DistributionBranchBase
+from gdm.distribution.components.base.distribution_switch_base import DistributionSwitchBase
 from gdm.distribution.components.distribution_bus import DistributionBus
 from gdm.distribution.components.distribution_load import DistributionLoad
 from gdm.distribution.components.distribution_solar import DistributionSolar
 from gdm.distribution.components.distribution_battery import DistributionBattery
+from gdm.distribution.components.geometry_branch import GeometryBranch
+from gdm.distribution.components.matrix_impedance_branch import MatrixImpedanceBranch
+from gdm.distribution.components.matrix_impedance_switch import MatrixImpedanceSwitch
 from gdm.distribution.distribution_system import (
     DistributionSystem,
     UserAttributes,
 )
+from gdm.distribution.equipment.matrix_impedance_switch_equipment import (
+    MatrixImpedanceSwitchEquipment,
+)
 from gdm.distribution.enums import Phase
 from gdm.distribution.sys_functools import (
     get_aggregated_load_time_series,
@@ -174,3 +182,174 @@ def reduce_to_primary_system(
         agg_timeseries,
         time_series_type,
     )
+
+
+def _normalize_cycle(cycle: list[str]) -> tuple[str, ...]:
+    """Rotate a cycle so that it starts at its lexicographically smallest bus."""
+    start = min(range(len(cycle)), key=cycle.__getitem__)
+    return tuple(cycle[start:] + cycle[:start])
+
+
+def _closed_edge_graph(system: DistributionSystem) -> nx.MultiGraph:
+    """Build an undirected graph containing only electrically closed edges.
+
+    Mirrors the edge filtering performed by ``get_directed_graph``: open switch
+    edges (any phase open) are excluded, while all other branches and
+    transformers count as closed.
+    """
+    graph = system.get_undirected_graph()
+    closed_graph = nx.MultiGraph()
+    for node in graph.nodes():
+        closed_graph.add_node(node)
+    for u, v, key, data in graph.edges(data=True, keys=True):
+        if data.get("is_closed"):
+            closed_graph.add_edge(u, v, key=key, **data)
+    return closed_graph
+
+
+def _cycle_components(
+    system: DistributionSystem,
+    cycle: tuple[str, ...],
+    graph: nx.MultiGraph,
+    is_candidate: Callable[[type], bool],
+) -> list[DistributionBranchBase]:
+    """Resolve the branch components on a cycle that match ``is_candidate``.
+
+    Components are deduplicated and returned sorted by name so that any
+    downstream selection is deterministic.
+    """
+    seen: dict[uuid.UUID, DistributionBranchBase] = {}
+    for i in range(len(cycle)):
+        bus_1, bus_2 = cycle[i], cycle[(i + 1) % len(cycle)]
+        if not graph.has_edge(bus_1, bus_2):
+            continue
+        for data in graph[bus_1][bus_2].values():
+            component_type = data["type"]
+            if is_candidate(component_type):
+                component = system.get_component(component_type, data["name"])
+                seen.setdefault(component.uuid, component)
+    return sorted(seen.values(), key=lambda component: component.name)
+
+
+def _switch_equipment_from_branch(
+    line: DistributionBranchBase,
+) -> MatrixImpedanceSwitchEquipment:
+    """Build switch equipment from a branch's impedance data."""
+    if isinstance(line, GeometryBranch):
+        line = line.to_matrix_representation(frequency_hz=60)
+    if not isinstance(line, MatrixImpedanceBranch):
+        msg = (
+            f"Cannot convert {line.__class__.__name__} to a switch. Only "
+            "MatrixImpedanceBranch and GeometryBranch lines are supported."
+        )
+        raise ValueError(msg)
+    return MatrixImpedanceSwitchEquipment(**line.equipment.model_dump(exclude_none=True))
+
+
+def _convert_line_to_open_switch(
+    system: DistributionSystem, line: DistributionBranchBase
+) -> MatrixImpedanceSwitch:
+    """Replace ``line`` with an open switch between the same buses."""
+    name = f"{line.name}_switch"
+    existing_names = {switch.name for switch in system.get_components(MatrixImpedanceSwitch)}
+    suffix = 1
+    while name in existing_names:
+        suffix += 1
+        name = f"{line.name}_switch_{suffix}"
+
+    switch = MatrixImpedanceSwitch(
+        buses=line.buses,
+        length=line.length,
+        phases=list(line.phases),
+        substation=line.substation,
+        feeder=line.feeder,
+        name=name,
+        is_closed=[False] * len(line.phases),
+        equipment=_switch_equipment_from_branch(line),
+        uuid=uuid.uuid5(uuid.NAMESPACE_URL, f"gdm.radial_switch.{name}"),
+    )
+    system.remove_component(line)
+    system.add_component(switch)
+    return switch
+
+
+def reduce_to_radial_network(dist_system: DistributionSystem, name: str) -> DistributionSystem:
+    """Reduce a looped distribution system to a radial network.
+
+    Every electrical loop in the input is broken deterministically: if a loop
+    contains a closed switch, that switch is opened; otherwise the shortest
+    line on the loop (ties broken by component name) is replaced with an open
+    ``MatrixImpedanceSwitch`` between the same buses. The input system is not
+    modified; a reduced copy is returned.
+
+    Parameters
+    ----------
+    dist_system : DistributionSystem
+        Source distribution system, possibly containing loops. Must be a single
+        connected component fed from one voltage source (the same requirement
+        as ``get_source_bus``).
+    name : str
+        Name of the reduced system.
+
+    Returns
+    -------
+    DistributionSystem
+        A radial copy of the input system in which no closed loop remains, so
+        that ``get_directed_graph`` prunes no edges and emits no "pruned from
+        DFS tree" warnings.
+
+    Notes
+    -----
+    - Geometry-based lines are converted to their matrix representation at
+      60 Hz when they are selected for conversion.
+    - Only cycle edges are touched, so the reduced network stays connected.
+    """
+    system = dist_system.deepcopy()
+    system.name = name
+    closed_graph = _closed_edge_graph(system)
+
+    while True:
+        cycles = sorted(
+            _normalize_cycle(cycle) for cycle in DistributionSystem.get_cycles(closed_graph)
+        )
+        if not cycles:
+            break
+        cycle = cycles[0]
+
+        switches = _cycle_components(
+            system, cycle, closed_graph, lambda t: issubclass(t, DistributionSwitchBase)
+        )
+        if switches:
+            switch = switches[0]
+            switch.is_closed = [False] * len(switch.phases)
+            bus_1, bus_2 = switch.buses[0].name, switch.buses[1].name
+            for key in [
+                k for k, data in closed_graph[bus_1][bus_2].items() if data["name"] == switch.name
+            ]:
+                closed_graph.remove_edge(bus_1, bus_2, key=key)
+            continue
+
+        lines = _cycle_components(
+            system,
+            cycle,
+            closed_graph,
+            lambda t: (
+                issubclass(t, DistributionBranchBase) and not issubclass(t, DistributionSwitchBase)
+            ),
+        )
+        if not lines:
+            msg = (
+                f"Cycle {list(cycle)} contains no branches or switches that can be "
+                "opened to break the loop; cannot reduce the system to a radial network."
+            )
+            raise ValueError(msg)
+
+        line = min(lines, key=lambda b: (b.length.to("meter").magnitude, b.name))
+        _convert_line_to_open_switch(system, line)
+        bus_1, bus_2 = line.buses[0].name, line.buses[1].name
+        for key in [
+            k for k, data in closed_graph[bus_1][bus_2].items() if data["name"] == line.name
+        ]:
+            closed_graph.remove_edge(bus_1, bus_2, key=key)
+
+    return system