From c9fc75ae389877d696f36a427c397271e1627a6d Mon Sep 17 00:00:00 2001 From: Bpetrick <170433522+BenjPetr@users.noreply.github.com> Date: Fri, 17 Apr 2026 09:29:08 +0000 Subject: [PATCH 1/7] refactor: prepare separation slack and reference bus --- src/dc_plus/importing/import_schema.py | 4 + .../importing/pandapower/pandapower_import.py | 13 +- .../importing/powsybl/powsybl_import.py | 38 +++++- src/dc_plus/interfaces/network_information.py | 36 ++++-- src/dc_plus/numpy/__init__.py | 4 + src/dc_plus/preprocess/create_network_data.py | 1 + .../preprocess/preprocess_jacobian_bsdf.py | 7 + .../importing/powsybl/test_powsybl_import.py | 15 ++- tests/importing/test_import_helpers.py | 3 + .../test_powsybl_vs_pandapower_dynamic.py | 2 + tests/numpy/test_injection_outage.py | 122 ++++++++++++++++++ 11 files changed, 217 insertions(+), 28 deletions(-) create mode 100644 tests/numpy/test_injection_outage.py diff --git a/src/dc_plus/importing/import_schema.py b/src/dc_plus/importing/import_schema.py index 5a4db84..81d10e0 100644 --- a/src/dc_plus/importing/import_schema.py +++ b/src/dc_plus/importing/import_schema.py @@ -114,6 +114,10 @@ class BusParamSchema(pa.DataFrameModel): voltage_magnitude: pat.Series[float] = pa.Field(coerce=True, nullable=True) voltage_angle: pat.Series[float] = pa.Field(coerce=True, nullable=True) bus_type: pat.Series[int] = pa.Field(coerce=True, description="0:Slack, 1:PV, 2:PQ") + is_angle_reference: pat.Series[bool] = pa.Field( + coerce=True, + description="True for the bus used as the angle reference, independent of its voltage-control mode.", + ) grid_island_id: pat.Series[int] = pa.Field( coerce=True, description="ID of the grid island the bus belongs to. 0 indicates the main grid." ) diff --git a/src/dc_plus/importing/pandapower/pandapower_import.py b/src/dc_plus/importing/pandapower/pandapower_import.py index 3337a70..ddf2c7c 100644 --- a/src/dc_plus/importing/pandapower/pandapower_import.py +++ b/src/dc_plus/importing/pandapower/pandapower_import.py @@ -272,11 +272,10 @@ def _get_buses_pandapower(net: pandapowerNet, slack_id: NotImplementedError) -> continue res_bus = net.res_bus.loc[idx] + is_angle_reference = idx == slack_id # Determine bus type - if idx == slack_id: - bus_type = BusType.SLACK - elif idx in net.gen["bus"].values or idx in net.sgen["bus"].values: + if idx in net.gen["bus"].values or idx in net.sgen["bus"].values: # Check if voltage control is enabled is_voltage_controlled = False if idx in net.gen["bus"].values: @@ -285,8 +284,13 @@ def _get_buses_pandapower(net: pandapowerNet, slack_id: NotImplementedError) -> if not is_voltage_controlled and idx in net.sgen["bus"].values: sgen_at_bus = net.sgen[net.sgen["bus"] == idx] is_voltage_controlled = sgen_at_bus.get("controllable", pd.Series([False])).any() + else: + is_voltage_controlled = is_angle_reference - bus_type = BusType.PV if is_voltage_controlled else BusType.PQ + if is_angle_reference and is_voltage_controlled: + bus_type = BusType.SLACK + elif is_voltage_controlled: + bus_type = BusType.PV else: bus_type = BusType.PQ @@ -302,6 +306,7 @@ def _get_buses_pandapower(net: pandapowerNet, slack_id: NotImplementedError) -> "voltage_magnitude": res_bus["vm_pu"] if pd.notna(res_bus["vm_pu"]) else np.nan, "voltage_angle": (np.deg2rad(res_bus["va_degree"]) if pd.notna(res_bus["va_degree"]) else np.nan), "bus_type": int(bus_type), + "is_angle_reference": is_angle_reference, "grid_island_id": grid_island_id, } ) diff --git a/src/dc_plus/importing/powsybl/powsybl_import.py b/src/dc_plus/importing/powsybl/powsybl_import.py index 1d0b45e..b2b1dc9 100644 --- a/src/dc_plus/importing/powsybl/powsybl_import.py +++ b/src/dc_plus/importing/powsybl/powsybl_import.py @@ -20,6 +20,7 @@ LimitParamSchema, ShuntParamSchema, ) +from dc_plus.interfaces.network_information import BusType logger = logging.getLogger(__name__) @@ -907,11 +908,21 @@ def _get_dangling_buses(net: Network) -> BusParamSchema: # TODO: get values depending on the connected state and the bus_id of the dangling line dangling_buses["grid_island_id"] = 0 dangling_buses["bus_type"] = 2 # PQ bus + dangling_buses["is_angle_reference"] = False bus_order_int_ids = _get_bus_ids_with_dangling_buses(net) dangling_buses = dangling_buses.merge(bus_order_int_ids, how="left", on=["id_str"]) dangling_buses = dangling_buses[ - ["id_str", "id_int", "name", "voltage_magnitude", "voltage_angle", "bus_type", "grid_island_id"] + [ + "id_str", + "id_int", + "name", + "voltage_magnitude", + "voltage_angle", + "bus_type", + "is_angle_reference", + "grid_island_id", + ] ] dangling_buses = BusParamSchema.validate(dangling_buses) return dangling_buses @@ -942,6 +953,7 @@ def _get_buses_powsybl(net: Network, slack_id: str, injections: InjectionParamSc buses.reset_index(drop=False, inplace=True) buses.rename(columns={"id": "id_str"}, inplace=True) buses["bus_type"] = 2 # PQ bus is default + buses["is_angle_reference"] = buses["v_angle"] == 0.0 buses.rename(columns={"v_mag": "voltage_magnitude", "v_angle": "voltage_angle"}, inplace=True) @@ -966,15 +978,29 @@ def _get_buses_powsybl(net: Network, slack_id: str, injections: InjectionParamSc buses = buses.merge(bus_order_int_ids, how="left", on=["id_str"]) # now the main grid "0" has only connected_component == 0 and synchronous_component == 0 - buses = buses[["id_str", "id_int", "name", "voltage_magnitude", "voltage_angle", "bus_type", "grid_island_id"]] + buses = buses[ + [ + "id_str", + "id_int", + "name", + "voltage_magnitude", + "voltage_angle", + "bus_type", + "is_angle_reference", + "grid_island_id", + ] + ] dangling_buses = _get_dangling_buses(net) buses = pd.concat([buses, dangling_buses], ignore_index=True) buses.sort_values(by=["id_int"], inplace=True) - # set bus types - pv_buses = injections[(injections["voltage_regulation"])]["regulated_bus_id_int"].unique() - buses.loc[buses["id_int"].isin(pv_buses), "bus_type"] = 1 # PV bus - buses.loc[buses["id_str"] == slack_id, "bus_type"] = 0 # slack bus + # Bus type tracks voltage control; the reference-angle role is stored separately. + voltage_control_buses = injections[(injections["connected"]) & (injections["voltage_regulation"])][ + "regulated_bus_id_int" + ].unique() + buses.loc[buses["id_int"].isin(voltage_control_buses), "bus_type"] = BusType.PV + slack_mask = buses["id_str"].eq(slack_id) & buses["id_int"].isin(voltage_control_buses) + buses.loc[slack_mask, "bus_type"] = BusType.SLACK buses = BusParamSchema.validate(buses) return buses diff --git a/src/dc_plus/interfaces/network_information.py b/src/dc_plus/interfaces/network_information.py index 0301d4b..0db7d39 100644 --- a/src/dc_plus/interfaces/network_information.py +++ b/src/dc_plus/interfaces/network_information.py @@ -349,6 +349,14 @@ class DynamicNetworkInformation: 2: pq """ + bus_is_angle_reference: Bool[np.ndarray, " n_buses"] + """Indicates which bus is used as the angle reference. + + This is distinct from the bus voltage-control mode. A reference bus can lose + voltage control after an outage and become a PQ bus while still fixing the + global angle reference. + """ + # injection data injection_to_bus: Int[np.ndarray, " n_injections n_timestep"] @@ -402,14 +410,14 @@ def n_branches(self) -> int: @property def slack_indices(self) -> np.ndarray: - """Return the indices of the slack bus. + """Return the indices of the angle-reference bus. Returns ------- slack_index : np.ndarray - The indices of the slack bus. + The indices of the angle-reference bus. """ - slack_index = np.flatnonzero(self.bus_type == BusType.SLACK) + slack_index = np.flatnonzero(self.bus_is_angle_reference) return slack_index def is_pv_bus(self, bus_index: int) -> bool: @@ -488,7 +496,7 @@ def pv_buses_mask(self) -> np.ndarray: pv_buses_mask : np.ndarray A boolean mask indicating which buses are PV buses. """ - pv_buses_mask = self.bus_type == BusType.PV + pv_buses_mask = (self.bus_type == BusType.PV) & ~self.bus_is_angle_reference return pv_buses_mask @property @@ -505,14 +513,14 @@ def pq_buses_mask(self) -> np.ndarray: @property def pvpq_buses_mask(self) -> np.ndarray: - """Return a boolean mask indicating which buses are PV or PQ buses. + """Return a boolean mask indicating which buses participate in angle equations. Returns ------- pvpq_buses_mask : np.ndarray A boolean mask indicating which buses are PV or PQ buses. """ - pvpq_buses_mask = (self.bus_type == BusType.PV) | (self.bus_type == BusType.PQ) + pvpq_buses_mask = ~self.bus_is_angle_reference return pvpq_buses_mask @property @@ -524,7 +532,7 @@ def pv_buses_indices(self) -> np.ndarray: pv_buses_indices : np.ndarray The indices of the PV buses. """ - pv_buses_indices = np.flatnonzero(self.bus_type == BusType.PV) + pv_buses_indices = np.flatnonzero(self.pv_buses_mask) return pv_buses_indices @property @@ -536,31 +544,32 @@ def pq_buses_indices(self) -> np.ndarray: pq_buses_indices : np.ndarray The indices of the PQ buses. """ - pq_buses_indices = np.flatnonzero(self.bus_type == BusType.PQ) + pq_buses_indices = np.flatnonzero(self.pq_buses_mask) return pq_buses_indices @property def pvpq_buses_indices(self) -> np.ndarray: - """Return the indices of the PV and PQ buses. + """Return the indices of buses participating in angle equations. Returns ------- pvpq_buses_indices : np.ndarray The indices of the PV and PQ buses. """ - pvpq_buses_indices = np.flatnonzero((self.bus_type == BusType.PV) | (self.bus_type == BusType.PQ)) + pvpq_buses_indices = np.flatnonzero(self.pvpq_buses_mask) return pvpq_buses_indices @property def pvpq_buses_indices_pvpq_order(self) -> np.ndarray: - """Return the indices of the PV and PQ buses. + """Return the indices of buses participating in angle equations. Returns ------- pvpq_buses_indices : np.ndarray The indices of the PV and PQ buses. """ - pvpq_buses_indices = np.concatenate((self.pv_buses_indices, self.pq_buses_indices)) + pq_angle_indices = self.pq_buses_indices[~self.bus_is_angle_reference[self.pq_buses_indices]] + pvpq_buses_indices = np.concatenate((self.pv_buses_indices, pq_angle_indices)) return pvpq_buses_indices # injections @@ -691,6 +700,9 @@ def _check_network_data_consistency( assert dynamic_network_data.bus_voltage_magnitudes.shape[0] == dynamic_network_data.bus_type.shape[0], ( "Inconsistent number of buses between bus_voltage_magnitudes and bus_type." ) + assert dynamic_network_data.bus_voltage_magnitudes.shape[0] == dynamic_network_data.bus_is_angle_reference.shape[0], ( + "Inconsistent number of buses between bus_voltage_magnitudes and bus_is_angle_reference." + ) # check injection data assert dynamic_network_data.injection_to_bus.shape[0] == dynamic_network_data.injection_active_power.shape[0], ( diff --git a/src/dc_plus/numpy/__init__.py b/src/dc_plus/numpy/__init__.py index 9c761f9..c2900d1 100644 --- a/src/dc_plus/numpy/__init__.py +++ b/src/dc_plus/numpy/__init__.py @@ -10,3 +10,7 @@ This implementation is intended for reference and testing purposes only. It is currently not optimized for performance, but for readability. """ + +from .injection_outage import injection_outage_dx + +__all__ = ["injection_outage_dx"] diff --git a/src/dc_plus/preprocess/create_network_data.py b/src/dc_plus/preprocess/create_network_data.py index d7f4470..2df1dd3 100644 --- a/src/dc_plus/preprocess/create_network_data.py +++ b/src/dc_plus/preprocess/create_network_data.py @@ -140,6 +140,7 @@ def _create_network_data( bus_active_power=bus_active_power, bus_reactive_power=bus_reactive_power, bus_type=buses["bus_type"].values.astype(int), + bus_is_angle_reference=buses["is_angle_reference"].values.astype(bool), injection_to_bus=injections["bus_index"].values.astype(int), injection_active_power=injections["p"].values, injection_reactive_power=injections["q"].values, diff --git a/src/dc_plus/preprocess/preprocess_jacobian_bsdf.py b/src/dc_plus/preprocess/preprocess_jacobian_bsdf.py index dbf7411..5f7c5de 100644 --- a/src/dc_plus/preprocess/preprocess_jacobian_bsdf.py +++ b/src/dc_plus/preprocess/preprocess_jacobian_bsdf.py @@ -302,6 +302,12 @@ def _pad_axis0(array: np.ndarray, value: float) -> np.ndarray: np.full(split_count, BusType.PQ, dtype=dynamic_network_data.bus_type.dtype), ) ) + extended_bus_is_angle_reference = np.concatenate( + ( + dynamic_network_data.bus_is_angle_reference, + np.zeros(split_count, dtype=dynamic_network_data.bus_is_angle_reference.dtype), + ) + ) return replace( dynamic_network_data, @@ -310,6 +316,7 @@ def _pad_axis0(array: np.ndarray, value: float) -> np.ndarray: bus_active_power=extended_bus_active_power, bus_reactive_power=extended_bus_reactive_power, bus_type=extended_bus_type, + bus_is_angle_reference=extended_bus_is_angle_reference, ) diff --git a/tests/importing/powsybl/test_powsybl_import.py b/tests/importing/powsybl/test_powsybl_import.py index 3d7df2c..5c66694 100644 --- a/tests/importing/powsybl/test_powsybl_import.py +++ b/tests/importing/powsybl/test_powsybl_import.py @@ -39,6 +39,8 @@ _get_trafo_parameter, ) +from dc_plus.interfaces.network_information import BusType + def test_get_tie_line_parameter(): net = pypowsybl.network.create_eurostag_tutorial_example1_with_tie_lines_and_areas() @@ -193,12 +195,13 @@ def test_get_buses(): bus_powsybl = net.get_buses(attributes=[]) dangling_count = len(_get_dangling_bus_ids(net)) assert len(bus_powsybl) + dangling_count == len(buses) - assert 0 in list(buses["bus_type"].values) # slack bus - assert 1 in list(buses["bus_type"].values) # PV bus - assert 2 in list(buses["bus_type"].values) # PQ bus - assert np.sum(buses["bus_type"] == 0) == 1 # only one slack bus - assert np.sum(buses["bus_type"] == 1) == 2 # there are two PV buses - assert np.sum(buses["bus_type"] == 2) >= 1 # at least one PQ bus + assert BusType.SLACK in list(buses["bus_type"].values) # slack bus + assert BusType.PV in list(buses["bus_type"].values) # PV bus + assert BusType.PQ in list(buses["bus_type"].values) # PQ bus + assert np.sum(buses["bus_type"] == BusType.SLACK) == 1 # only one slack bus + assert np.sum(buses["is_angle_reference"]) == 1 # only one reference bus + assert np.sum(buses["bus_type"] == BusType.PV) == 2 # there are two PV buses + assert np.sum(buses["bus_type"] == BusType.PQ) >= 1 # at least one PQ bus def test_no_synchronous_components_returns_connected_ids(): diff --git a/tests/importing/test_import_helpers.py b/tests/importing/test_import_helpers.py index 1e43e76..5633238 100644 --- a/tests/importing/test_import_helpers.py +++ b/tests/importing/test_import_helpers.py @@ -23,6 +23,7 @@ def test_remove_isolated_buses_injections(): "voltage_magnitude": [1.0, 0.98, 1.02], "voltage_angle": [0.0, -0.05, 0.03], "bus_type": [0, 1, 2], + "is_angle_reference": [True, False, False], "grid_island_id": [0, 0, 1], } ) @@ -73,6 +74,7 @@ def test_remove_isolated_branches(): "voltage_magnitude": [1.0, 1.0, 1.0], "voltage_angle": [0.0, 0.0, 0.0], "bus_type": [0, 1, 2], + "is_angle_reference": [True, False, False], "grid_island_id": [0, 0, 1], } ) @@ -125,6 +127,7 @@ def test_remove_isolated_buses(): "voltage_magnitude": [1.0, 0.97, 1.01, 0.95], "voltage_angle": [0.0, -0.03, 0.02, -0.05], "bus_type": [0, 1, 2, 2], + "is_angle_reference": [True, False, False, False], "grid_island_id": [0, 1, 0, 2], } ) diff --git a/tests/importing/test_powsybl_vs_pandapower_dynamic.py b/tests/importing/test_powsybl_vs_pandapower_dynamic.py index e331c58..3bcd137 100644 --- a/tests/importing/test_powsybl_vs_pandapower_dynamic.py +++ b/tests/importing/test_powsybl_vs_pandapower_dynamic.py @@ -296,6 +296,7 @@ def _prepare_powsybl_network( bus_active_power=bus_active_power, bus_reactive_power=bus_reactive_power, bus_type=bus_type, + bus_is_angle_reference=dynamic_psb.bus_is_angle_reference[bus_perm].copy(), injection_to_bus=injection_to_bus, injection_active_power=injection_active_power, injection_reactive_power=injection_reactive_power, @@ -384,6 +385,7 @@ def _prepare_pandapower_network( bus_active_power=dynamic_pp.bus_active_power[bus_perm].copy(), bus_reactive_power=dynamic_pp.bus_reactive_power[bus_perm].copy(), bus_type=dynamic_pp.bus_type[bus_perm].copy(), + bus_is_angle_reference=dynamic_pp.bus_is_angle_reference[bus_perm].copy(), injection_to_bus=dynamic_pp.injection_to_bus[injection_perm].copy(), injection_active_power=dynamic_pp.injection_active_power[injection_perm].copy(), injection_reactive_power=dynamic_pp.injection_reactive_power[injection_perm].copy(), diff --git a/tests/numpy/test_injection_outage.py b/tests/numpy/test_injection_outage.py new file mode 100644 index 0000000..105fb49 --- /dev/null +++ b/tests/numpy/test_injection_outage.py @@ -0,0 +1,122 @@ +# Copyright 2026 50Hertz Transmission GmbH and Elia Transmission Belgium SA/NV +# +# This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0. +# If a copy of the MPL was not distributed with this file, +# you can obtain one at https://mozilla.org/MPL/2.0/. +# Mozilla Public License, version 2.0 + +from copy import deepcopy + +import numpy as np +import pypowsybl +import pytest + +from dc_plus.example_grids.pypowsbl.example_grids import PANDAPOWER_NETWORKS_FOR_POWSYBL, POWSYBL_NETWORKS +from dc_plus.importing.powsybl.powsybl_import import _get_injections_powsybl +from dc_plus.importing.powsybl.powsybl_loadflow_parameter import get_powsybl_loadflow_parameter +from dc_plus.importing.powsybl.powsybl_network_helpers import _load_test_grid +from dc_plus.interfaces.network_information import BusType +from dc_plus.numpy.injection_outage import injection_outage_dx +from dc_plus.preprocess.create_network_data import create_network_data + +powsybl_networks = POWSYBL_NETWORKS +pandapower_networks = PANDAPOWER_NETWORKS_FOR_POWSYBL +TOL = 1e-10 + + +def _connected_injection_indices(injections): + candidates = np.flatnonzero(injections["connected"].to_numpy(dtype=bool)) + if candidates.size == 0: + pytest.skip("No connected injections available for an injection outage test") + return candidates + + +def _expected_bus_delta(dynamic_info, injection_idx, outage_p, outage_q): + delta_p = np.zeros(dynamic_info.n_buses, dtype=float) + delta_q = np.zeros(dynamic_info.n_buses, dtype=float) + bus_idx = int(dynamic_info.injection_to_bus[injection_idx]) + delta_p[bus_idx] -= float(outage_p[injection_idx]) + delta_q[bus_idx] -= float(outage_q[injection_idx]) + return delta_p, delta_q + + +@pytest.mark.parametrize("get_net", powsybl_networks + pandapower_networks) +def test_injection_outage_numpy_compare_powsybl(get_net): + net, _, dynamic_info, string_info, jacobian_data = _load_test_grid(get_net) + injections = _get_injections_powsybl(net).reset_index(drop=True) + injection_powsybl_type = net.get_injections(attributes=["type"]) + outage_p = injections["setpoint_p"].fillna(injections["p"]).to_numpy(dtype=float) + outage_q = injections["setpoint_q"].fillna(injections["q"]).to_numpy(dtype=float) + + assert len(injections) == dynamic_info.n_injections + + compared_cases = 0 + for injection_idx in _connected_injection_indices(injections): + injection_id = injections.loc[injection_idx, "id_str"] + if (injection_powsybl_type.loc[injection_id, "type"]) in (["HVDC_CONVERTER_STATION"]): + # HVDC would need a multi outage + continue + + net_n1 = deepcopy(net) + net_n1.remove_elements(injection_id) + + loadflow_parameter = get_powsybl_loadflow_parameter("one_step") + loadflow_res = pypowsybl.loadflow.run_ac(net_n1, parameters=loadflow_parameter)[0] + assert loadflow_res.iteration_count == 1 + if loadflow_res.status != pypowsybl._pypowsybl.LoadFlowComponentStatus.CONVERGED: + continue + + _, dynamic_info_n1, string_info_n1 = create_network_data(net_n1) + np.testing.assert_array_equal(string_info.bus_ids, string_info_n1.bus_ids) + + same_angle_structure = np.array_equal( + dynamic_info.pvpq_buses_indices_pvpq_order, + dynamic_info_n1.pvpq_buses_indices_pvpq_order, + ) + same_magnitude_structure = np.array_equal( + dynamic_info.pq_buses_indices, + dynamic_info_n1.pq_buses_indices, + ) + if not (same_angle_structure and same_magnitude_structure): + continue + + expected_delta_p, expected_delta_q = _expected_bus_delta(dynamic_info, injection_idx, outage_p, outage_q) + actual_delta_p = dynamic_info_n1.bus_active_power - dynamic_info.bus_active_power + actual_delta_q = dynamic_info_n1.bus_reactive_power - dynamic_info.bus_reactive_power + p_buses = dynamic_info.pvpq_buses_indices_pvpq_order + q_buses = dynamic_info.pq_buses_indices + if not ( + np.allclose(actual_delta_p[p_buses], expected_delta_p[p_buses], rtol=TOL, atol=TOL) + and np.allclose(actual_delta_q[q_buses], expected_delta_q[q_buses], rtol=TOL, atol=TOL) + ): + continue + + dx = injection_outage_dx( + jacobian_inv=jacobian_data.inverse_jacobian, + outage_injection_indices=np.array([injection_idx], dtype=np.int64), + injection_to_bus=dynamic_info.injection_to_bus, + injection_active_power=outage_p, + injection_reactive_power=outage_q, + angle_component_indices=jacobian_data.angle_component_indices, + magnitude_component_indices=jacobian_data.magnitude_component_indices, + ) + + theta_updated = dynamic_info.bus_voltage_angles_rad.copy() + vm_updated = dynamic_info.bus_voltage_magnitudes.copy() + theta_updated[dynamic_info.pvpq_buses_indices_pvpq_order] += dx[jacobian_data.is_angle_component] + vm_updated[dynamic_info.pq_buses_indices] += dx[jacobian_data.is_magnitude_component] + np.testing.assert_allclose( + theta_updated, + dynamic_info_n1.bus_voltage_angles_rad, + rtol=TOL, + atol=TOL, + ) + np.testing.assert_allclose( + vm_updated, + dynamic_info_n1.bus_voltage_magnitudes, + rtol=TOL, + atol=TOL, + ) + compared_cases += 1 + + assert compared_cases > 0, f"No injection outages matched the fixed-injection assumptions for {get_net.__name__}" From 5efee3b9ed3e1c410ab6f8e517c5f60c562ccf1d Mon Sep 17 00:00:00 2001 From: Bpetrick <170433522+BenjPetr@users.noreply.github.com> Date: Fri, 17 Apr 2026 09:39:51 +0000 Subject: [PATCH 2/7] feat: draft injection outage --- src/dc_plus/numpy/injection_outage.py | 79 +++++++++++++++++++++++++++ tests/numpy/test_injection_outage.py | 2 +- 2 files changed, 80 insertions(+), 1 deletion(-) create mode 100644 src/dc_plus/numpy/injection_outage.py diff --git a/src/dc_plus/numpy/injection_outage.py b/src/dc_plus/numpy/injection_outage.py new file mode 100644 index 0000000..2942f54 --- /dev/null +++ b/src/dc_plus/numpy/injection_outage.py @@ -0,0 +1,79 @@ +# Copyright 2026 50Hertz Transmission GmbH and Elia Transmission Belgium SA/NV +# +# This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0. +# If a copy of the MPL was not distributed with this file, +# you can obtain one at https://mozilla.org/MPL/2.0/. +# Mozilla Public License, version 2.0 + +"""One-step state updates for injection outages.""" + +import numpy as np +from jaxtyping import Float, Int + + +def injection_outage_dx( + jacobian_inv: Float[np.ndarray, " n_eq n_eq"], + outage_injection_indices: Int[np.ndarray, " n_outages"], + injection_to_bus: Int[np.ndarray, " n_injections"], + injection_active_power: Float[np.ndarray, " n_injections"], + injection_reactive_power: Float[np.ndarray, " n_injections"], + angle_component_indices: Int[np.ndarray, " n_buses"], + magnitude_component_indices: Int[np.ndarray, " n_buses"], +) -> Float[np.ndarray, " n_eq"]: + """Return the one-step state increment caused by disconnecting injections. + + The Jacobian is unchanged for a pure injection outage because network + topology and admittances remain fixed. The outage only changes the nodal + power mismatch at the buses hosting the disconnected injections. + + This helper assumes the bus type pattern used by the Jacobian remains + unchanged across the outage. That is appropriate for load outages and other + contingencies that do not alter voltage-control structure. + + Parameters + ---------- + jacobian_inv : Float[np.ndarray, " n_eq n_eq"] + Base inverse Jacobian at the hot-start operating point. + outage_injection_indices : Int[np.ndarray, " n_outages"] + Indices of injections to disconnect. + injection_to_bus : Int[np.ndarray, " n_injections"] + Mapping from injection index to bus index. + injection_active_power : Float[np.ndarray, " n_injections"] + Active power removed when each injection is disconnected. + In practice this is typically the injection setpoint, or the solved + injection power when no separate setpoint is available. + injection_reactive_power : Float[np.ndarray, " n_injections"] + Reactive power removed when each injection is disconnected. + In practice this is typically the injection setpoint, or the solved + injection power when no separate setpoint is available. + angle_component_indices : Int[np.ndarray, " n_buses"] + Bus to active-power Jacobian row mapping. + magnitude_component_indices : Int[np.ndarray, " n_buses"] + Bus to reactive-power Jacobian row mapping. + + Returns + ------- + Float[np.ndarray, " n_eq"] + One-step update vector in Jacobian ordering. + """ + outage_indices = outage_injection_indices + mismatch = np.zeros(jacobian_inv.shape[0], dtype=jacobian_inv.dtype) + if outage_indices.size == 0: + return mismatch + + if np.any((outage_indices < 0) | (outage_indices >= injection_to_bus.size)): + raise IndexError("Injection outage index is out of bounds") + + buses = injection_to_bus[outage_indices] + + p_idx = angle_component_indices[buses] + valid_p = p_idx >= 0 + if np.any(valid_p): + np.add.at(mismatch, p_idx[valid_p], -injection_active_power[outage_indices][valid_p]) + + q_idx = magnitude_component_indices[buses] + valid_q = q_idx >= 0 + if np.any(valid_q): + np.add.at(mismatch, q_idx[valid_q], -injection_reactive_power[outage_indices][valid_q]) + + return -jacobian_inv @ mismatch diff --git a/tests/numpy/test_injection_outage.py b/tests/numpy/test_injection_outage.py index 105fb49..e205db0 100644 --- a/tests/numpy/test_injection_outage.py +++ b/tests/numpy/test_injection_outage.py @@ -62,7 +62,7 @@ def test_injection_outage_numpy_compare_powsybl(get_net): loadflow_parameter = get_powsybl_loadflow_parameter("one_step") loadflow_res = pypowsybl.loadflow.run_ac(net_n1, parameters=loadflow_parameter)[0] - assert loadflow_res.iteration_count == 1 + assert loadflow_res.iteration_count <= 1 if loadflow_res.status != pypowsybl._pypowsybl.LoadFlowComponentStatus.CONVERGED: continue From 8ee061724cc9143678b420e640cce0f10a09e68a Mon Sep 17 00:00:00 2001 From: Bpetrick <170433522+BenjPetr@users.noreply.github.com> Date: Sat, 18 Apr 2026 02:57:02 +0000 Subject: [PATCH 3/7] feat: finish injection outage numpy + jax + doc --- docs/dc_plus/injection_outage.md | 32 ++ src/dc_plus/jax/injection_outage.py | 467 +++++++++++++++++++++++++ src/dc_plus/numpy/__init__.py | 4 +- src/dc_plus/numpy/injection_outage.py | 149 ++++++-- tests/jax/test_injection_outage_jax.py | 241 +++++++++++++ tests/numpy/test_injection_outage.py | 103 ++++-- 6 files changed, 935 insertions(+), 61 deletions(-) create mode 100644 docs/dc_plus/injection_outage.md create mode 100644 src/dc_plus/jax/injection_outage.py create mode 100644 tests/jax/test_injection_outage_jax.py diff --git a/docs/dc_plus/injection_outage.md b/docs/dc_plus/injection_outage.md new file mode 100644 index 0000000..0d4e5e4 --- /dev/null +++ b/docs/dc_plus/injection_outage.md @@ -0,0 +1,32 @@ +# Injection Outage + +`non_voltage_regulating_injection_outage_dx` computes one-step AC state updates for a batch of injection outages around the N-0 operating point. + +The Jacobian is kept fixed. Each outage only contributes an active and reactive power mismatch at the affected bus, and the batched update is evaluated with the base inverse Jacobian. + +$$ +\Delta X = - M J^{-T} +$$ + +where `J` is the base-case Jacobian and `M` is the batch of outage mismatch vectors in Jacobian ordering. + +## Batch API + +The function expects a 2D outage array with shape `(n_contingencies, n_outages)`. + +- Each row is one contingency. +- A row may contain multiple outage indices. +- Negative entries are ignored and can be used as padding. +- The return value has shape `(n_contingencies, n_eq)`. + +## Scope + +This update is intended for non-voltage-regulating injections such as loads and other fixed-power injections. + +## Generator Behavior + +The function can also be applied to voltage-regulating generator outages, but this is only an approximation. A voltage-regulating generator outage should normally trigger a PV to PQ bus-type change. The current implementation does not model that switch, so it does not capture the resulting voltage drop at the regulated bus. + +## HVDC + +- HVDC converter outages should be outaged in pairs. DC+ models the HVDC line as a pair of generator + load diff --git a/src/dc_plus/jax/injection_outage.py b/src/dc_plus/jax/injection_outage.py new file mode 100644 index 0000000..d2ae4a6 --- /dev/null +++ b/src/dc_plus/jax/injection_outage.py @@ -0,0 +1,467 @@ +# Copyright 2026 50Hertz Transmission GmbH and Elia Transmission Belgium SA/NV +# +# This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0. +# If a copy of the MPL was not distributed with this file, +# you can obtain one at https://mozilla.org/MPL/2.0/. +# Mozilla Public License, version 2.0 + +"""One-step state updates for injection outages using JAX.""" + +import jax +import jax.numpy as jnp +from jax_dataclasses import pytree_dataclass +from jaxtyping import Array, Complex128, Float, Int + +from .lodf_branches import _compute_monitored_branch_currents, _prepare_monitored_branch_pack +from .lodf_voltages import build_monitor_rows + +# ruff: noqa: PLR0913 + + +@pytree_dataclass +class InjectionOutageMonitoredResults: + """One-step post-contingency results restricted to monitored elements.""" + + n_1_theta: Float[Array, "... n_contingencies n_buses_monitored"] + n_1_voltage: Float[Array, "... n_contingencies n_buses_monitored"] + + n_1_p_from: Float[Array, "... n_contingencies n_branches_monitored"] + n_1_p_to: Float[Array, "... n_contingencies n_branches_monitored"] + n_1_q_from: Float[Array, "... n_contingencies n_branches_monitored"] + n_1_q_to: Float[Array, "... n_contingencies n_branches_monitored"] + n_1_i_from: Complex128[Array, "... n_contingencies n_branches_monitored"] + n_1_i_to: Complex128[Array, "... n_contingencies n_branches_monitored"] + + +def _collect_single_injection_outage_mismatch( + outage_injection_indices: Int[jnp.ndarray, " n_outages"], + jacobian_size: int, + injection_to_bus: Int[jnp.ndarray, " n_injections"], + injection_active_power: Float[jnp.ndarray, " n_injections"], + injection_reactive_power: Float[jnp.ndarray, " n_injections"], + angle_component_indices: Int[jnp.ndarray, " n_buses"], + magnitude_component_indices: Int[jnp.ndarray, " n_buses"], + dtype: jnp.dtype, +) -> Float[jnp.ndarray, " n_eq"]: + """Assemble the outage-induced mismatch vector for one contingency. + + Parameters + ---------- + outage_injection_indices : Int[jnp.ndarray, " n_outages"] + Indices of injections to disconnect for one contingency. Negative entries + are ignored and can be used as padding. + jacobian_size : int + Number of equations in the Jacobian system. + injection_to_bus : Int[jnp.ndarray, " n_injections"] + Map from injection index to hosting bus index. + injection_active_power : Float[jnp.ndarray, " n_injections"] + Active power associated with each injection. + injection_reactive_power : Float[jnp.ndarray, " n_injections"] + Reactive power associated with each injection. + angle_component_indices : Int[jnp.ndarray, " n_buses"] + Map from bus index to active-power mismatch component index. + magnitude_component_indices : Int[jnp.ndarray, " n_buses"] + Map from bus index to reactive-power mismatch component index. + dtype : jnp.dtype + Target dtype of the assembled mismatch vector. + + Returns + ------- + Float[jnp.ndarray, " n_eq"] + Mismatch vector in Jacobian component ordering. + """ + valid_outage_mask = outage_injection_indices >= 0 + safe_outages = jnp.where(valid_outage_mask, outage_injection_indices, 0) + buses = injection_to_bus[safe_outages] + + mismatch = jnp.zeros((jacobian_size,), dtype=dtype) + + p_idx = angle_component_indices[buses] + valid_p = valid_outage_mask & (p_idx >= 0) + safe_p_idx = jnp.where(valid_p, p_idx, 0) + p_values = jnp.where(valid_p, -injection_active_power[safe_outages], 0.0) + mismatch = mismatch.at[safe_p_idx].add(p_values) + + q_idx = magnitude_component_indices[buses] + valid_q = valid_outage_mask & (q_idx >= 0) + safe_q_idx = jnp.where(valid_q, q_idx, 0) + q_values = jnp.where(valid_q, -injection_reactive_power[safe_outages], 0.0) + mismatch = mismatch.at[safe_q_idx].add(q_values) + + return mismatch + + +def _collect_injection_outage_mismatch( + jacobian_size: int, + outage_injection_indices: Int[jnp.ndarray, " n_contingencies n_outages"], + injection_to_bus: Int[jnp.ndarray, " n_injections"], + injection_active_power: Float[jnp.ndarray, " n_injections"], + injection_reactive_power: Float[jnp.ndarray, " n_injections"], + angle_component_indices: Int[jnp.ndarray, " n_buses"], + magnitude_component_indices: Int[jnp.ndarray, " n_buses"], + dtype: jnp.dtype, +) -> Float[jnp.ndarray, " n_contingencies n_eq"]: + """Assemble outage-induced mismatch vectors in Jacobian ordering. + + Parameters + ---------- + jacobian_size : int + Number of equations in the Jacobian system. + outage_injection_indices : Int[jnp.ndarray, " n_contingencies n_outages"] + Batched injection outage indices. Each row represents one contingency and + may include negative padding entries. + injection_to_bus : Int[jnp.ndarray, " n_injections"] + Map from injection index to hosting bus index. + injection_active_power : Float[jnp.ndarray, " n_injections"] + Active power associated with each injection. + injection_reactive_power : Float[jnp.ndarray, " n_injections"] + Reactive power associated with each injection. + angle_component_indices : Int[jnp.ndarray, " n_buses"] + Map from bus index to active-power mismatch component index. + magnitude_component_indices : Int[jnp.ndarray, " n_buses"] + Map from bus index to reactive-power mismatch component index. + dtype : jnp.dtype + Target dtype of the assembled mismatch matrix. + + Returns + ------- + Float[jnp.ndarray, " n_contingencies n_eq"] + Batched mismatch vectors in Jacobian component ordering. + """ + + def collect_single(outage_row: Int[jnp.ndarray, " n_outages"]) -> Float[jnp.ndarray, " n_eq"]: + return _collect_single_injection_outage_mismatch( + outage_injection_indices=outage_row, + jacobian_size=jacobian_size, + injection_to_bus=injection_to_bus, + injection_active_power=injection_active_power, + injection_reactive_power=injection_reactive_power, + angle_component_indices=angle_component_indices, + magnitude_component_indices=magnitude_component_indices, + dtype=dtype, + ) + + return jax.vmap(collect_single)(outage_injection_indices) + + +@jax.jit +def _calculate_injection_outage_dx( + jacobian_inv_transposed: Float[jnp.ndarray, " n_eq n_eq"], + mismatch: Float[jnp.ndarray, " n_contingencies n_eq"], +) -> Float[jnp.ndarray, " n_contingencies n_eq"]: + """Map outage mismatch vectors to one-step state increments. + + Parameters + ---------- + jacobian_inv_transposed : Float[jnp.ndarray, " n_eq n_eq"] + Transposed inverse Jacobian at the base operating point. + mismatch : Float[jnp.ndarray, " n_contingencies n_eq"] + Batched outage-induced mismatch vectors. + + Returns + ------- + Float[jnp.ndarray, " n_contingencies n_eq"] + Batched one-step state increments in Jacobian ordering. + """ + return -(mismatch @ jacobian_inv_transposed) + + +@jax.jit +def _calculate_monitor_bus_state_updates( + jacobian_inv_transposed: Float[jnp.ndarray, " n_eq n_eq"], + mismatch: Float[jnp.ndarray, " n_contingencies n_eq"], + theta_rows: Int[jnp.ndarray, " n_mon_bus"], + vm_rows: Int[jnp.ndarray, " n_mon_bus"], + theta_mask: Float[jnp.ndarray, " n_mon_bus"], + vm_mask: Float[jnp.ndarray, " n_mon_bus"], +) -> tuple[ + Float[jnp.ndarray, " n_contingencies n_mon_bus"], + Float[jnp.ndarray, " n_contingencies n_mon_bus"], +]: + """Map outage mismatch vectors to monitored bus state increments only. + + Parameters + ---------- + jacobian_inv_transposed : Float[jnp.ndarray, " n_eq n_eq"] + Transposed inverse Jacobian at the base operating point. + mismatch : Float[jnp.ndarray, " n_contingencies n_eq"] + Batched outage-induced mismatch vectors. + theta_rows : Int[jnp.ndarray, " n_mon_bus"] + Safe Jacobian row indices for monitored bus angles. + vm_rows : Int[jnp.ndarray, " n_mon_bus"] + Safe Jacobian row indices for monitored bus magnitudes. + theta_mask : Float[jnp.ndarray, " n_mon_bus"] + Mask for monitored buses with angle states in the Jacobian. + vm_mask : Float[jnp.ndarray, " n_mon_bus"] + Mask for monitored buses with magnitude states in the Jacobian. + + Returns + ------- + tuple[Float[jnp.ndarray, " n_contingencies n_mon_bus"], Float[jnp.ndarray, " n_contingencies n_mon_bus"]] + Monitored angle and magnitude increments for each contingency. + """ + dtype = jacobian_inv_transposed.dtype + theta_mask_d = theta_mask.astype(dtype) + vm_mask_d = vm_mask.astype(dtype) + theta_dx = -(mismatch @ jacobian_inv_transposed[:, theta_rows]) * theta_mask_d[None, :] + vm_dx = -(mismatch @ jacobian_inv_transposed[:, vm_rows]) * vm_mask_d[None, :] + return theta_dx, vm_dx + + +def non_voltage_regulating_injection_outage_monitor_buses( + jacobian_inv_transposed: Float[jnp.ndarray, " n_eq n_eq"], + outage_injection_indices: Int[jnp.ndarray, " n_contingencies n_outages"], + injection_to_bus: Int[jnp.ndarray, " n_injections"], + injection_active_power: Float[jnp.ndarray, " n_injections"], + injection_reactive_power: Float[jnp.ndarray, " n_injections"], + angle_component_indices: Int[jnp.ndarray, " n_buses"], + magnitude_component_indices: Int[jnp.ndarray, " n_buses"], + monitor_bus_indices: Int[jnp.ndarray, " n_mon_bus"], + v_mag_hat: Float[jnp.ndarray, " n_buses"], + theta_hat: Float[jnp.ndarray, " n_buses"], +) -> tuple[ + Float[jnp.ndarray, " n_contingencies n_mon_bus"], + Float[jnp.ndarray, " n_contingencies n_mon_bus"], +]: + """Compute monitored post-contingency bus voltages for injection outages only. + + Parameters + ---------- + jacobian_inv_transposed : Float[jnp.ndarray, " n_eq n_eq"] + Base transposed inverse Jacobian at the hot-start operating point. + outage_injection_indices : Int[jnp.ndarray, " n_contingencies n_outages"] + Batched indices of injections to disconnect. Each row is one contingency. + injection_to_bus : Int[jnp.ndarray, " n_injections"] + Map from injection index to hosting bus index. + injection_active_power : Float[jnp.ndarray, " n_injections"] + Active power of each injection at the hot-start operating point. + injection_reactive_power : Float[jnp.ndarray, " n_injections"] + Reactive power of each injection at the hot-start operating point. + angle_component_indices : Int[jnp.ndarray, " n_buses"] + Map from bus index to angle component index in the Jacobian. -1 if not present. + magnitude_component_indices : Int[jnp.ndarray, " n_buses"] + Map from bus index to magnitude component index in the Jacobian. -1 if not present. + monitor_bus_indices : Int[jnp.ndarray, " n_mon_bus"] + Indices of buses whose post-contingency voltages should be returned. + v_mag_hat : Float[jnp.ndarray, " n_buses"] + Base-case voltage magnitudes. + theta_hat : Float[jnp.ndarray, " n_buses"] + Base-case voltage angles. + + Returns + ------- + tuple[Float[jnp.ndarray, " n_contingencies n_mon_bus"], Float[jnp.ndarray, " n_contingencies n_mon_bus"]] + Post-contingency monitored voltage angles and magnitudes. + """ + dtype = jacobian_inv_transposed.dtype + mismatch = _collect_injection_outage_mismatch( + jacobian_size=jacobian_inv_transposed.shape[0], + outage_injection_indices=outage_injection_indices, + injection_to_bus=injection_to_bus, + injection_active_power=injection_active_power, + injection_reactive_power=injection_reactive_power, + angle_component_indices=angle_component_indices, + magnitude_component_indices=magnitude_component_indices, + dtype=dtype, + ) + theta_rows, vm_rows, theta_mask, vm_mask = build_monitor_rows( + angle_component_indices=angle_component_indices, + magnitude_component_indices=magnitude_component_indices, + monitor_bus_indices=monitor_bus_indices, + ) + theta_dx, vm_dx = _calculate_monitor_bus_state_updates( + jacobian_inv_transposed=jacobian_inv_transposed, + mismatch=mismatch, + theta_rows=theta_rows, + vm_rows=vm_rows, + theta_mask=theta_mask, + vm_mask=vm_mask, + ) + base_theta0 = jnp.take(theta_hat, monitor_bus_indices, axis=0).astype(dtype) + base_vm0 = jnp.take(v_mag_hat, monitor_bus_indices, axis=0).astype(dtype) + return base_theta0[None, :] + theta_dx, base_vm0[None, :] + vm_dx + + +def non_voltage_regulating_injection_outage_monitored( + jacobian_inv_transposed: Float[jnp.ndarray, " n_eq n_eq"], + outage_injection_indices: Int[jnp.ndarray, " n_contingencies n_outages"], + injection_to_bus: Int[jnp.ndarray, " n_injections"], + injection_active_power: Float[jnp.ndarray, " n_injections"], + injection_reactive_power: Float[jnp.ndarray, " n_injections"], + angle_component_indices: Int[jnp.ndarray, " n_buses"], + magnitude_component_indices: Int[jnp.ndarray, " n_buses"], + monitor_bus_indices: Int[jnp.ndarray, " n_mon_bus"], + v_mag_hat: Float[jnp.ndarray, " n_buses"], + theta_hat: Float[jnp.ndarray, " n_buses"], + branch_from: Int[jnp.ndarray, " n_branches"], + branch_to: Int[jnp.ndarray, " n_branches"], + y_ff: Complex128[jnp.ndarray, " n_branches"], + y_ft: Complex128[jnp.ndarray, " n_branches"], + y_tf: Complex128[jnp.ndarray, " n_branches"], + y_tt: Complex128[jnp.ndarray, " n_branches"], + monitor_branch_indices: Int[jnp.ndarray, " n_mon_br"], + bus_to_mon_index: Int[jnp.ndarray, " n_buses"], +) -> InjectionOutageMonitoredResults: + """Compute monitored bus states and branch flows for injection outages. + + Parameters + ---------- + jacobian_inv_transposed : Float[jnp.ndarray, " n_eq n_eq"] + Base transposed inverse Jacobian at the hot-start operating point. + outage_injection_indices : Int[jnp.ndarray, " n_contingencies n_outages"] + Batched indices of injections to disconnect. Each row is one contingency. + injection_to_bus : Int[jnp.ndarray, " n_injections"] + Map from injection index to hosting bus index. + injection_active_power : Float[jnp.ndarray, " n_injections"] + Active power of each injection at the hot-start operating point. + injection_reactive_power : Float[jnp.ndarray, " n_injections"] + Reactive power of each injection at the hot-start operating point. + angle_component_indices : Int[jnp.ndarray, " n_buses"] + Map from bus index to angle component index in the Jacobian. -1 if not present. + magnitude_component_indices : Int[jnp.ndarray, " n_buses"] + Map from bus index to magnitude component index in the Jacobian. -1 if not present. + monitor_bus_indices : Int[jnp.ndarray, " n_mon_bus"] + Indices of buses whose post-contingency voltages should be returned. + v_mag_hat : Float[jnp.ndarray, " n_buses"] + Base-case voltage magnitudes. + theta_hat : Float[jnp.ndarray, " n_buses"] + Base-case voltage angles. + branch_from : Int[jnp.ndarray, " n_branches"] + From bus index of each branch. + branch_to : Int[jnp.ndarray, " n_branches"] + To bus index of each branch. + y_ff : Complex128[jnp.ndarray, " n_branches"] + From-from branch admittances. + y_ft : Complex128[jnp.ndarray, " n_branches"] + From-to branch admittances. + y_tf : Complex128[jnp.ndarray, " n_branches"] + To-from branch admittances. + y_tt : Complex128[jnp.ndarray, " n_branches"] + To-to branch admittances. + monitor_branch_indices : Int[jnp.ndarray, " n_mon_br"] + Indices of monitored branches. + bus_to_mon_index : Int[jnp.ndarray, " n_buses"] + Mapping from global bus index to monitored-bus position. Use -1 for unmonitored buses. + + Returns + ------- + InjectionOutageMonitoredResults + Post-contingency monitored bus states and monitored branch flows. + """ + dtype = jacobian_inv_transposed.dtype + theta_all, vm_all = non_voltage_regulating_injection_outage_monitor_buses( + jacobian_inv_transposed=jacobian_inv_transposed, + outage_injection_indices=outage_injection_indices, + injection_to_bus=injection_to_bus, + injection_active_power=injection_active_power, + injection_reactive_power=injection_reactive_power, + angle_component_indices=angle_component_indices, + magnitude_component_indices=magnitude_component_indices, + monitor_bus_indices=monitor_bus_indices, + v_mag_hat=v_mag_hat, + theta_hat=theta_hat, + ) + ( + _, + y_ff_mon, + y_ft_mon, + y_tf_mon, + y_tt_mon, + f_pos_safe, + t_pos_safe, + end_mask, + ) = _prepare_monitored_branch_pack( + branch_from=branch_from, + branch_to=branch_to, + y_ff=y_ff, + y_ft=y_ft, + y_tf=y_tf, + y_tt=y_tt, + monitor_branch_indices=monitor_branch_indices, + bus_to_mon_index=bus_to_mon_index, + dtype=dtype, + ) + v_from_all, v_to_all, i_from_all, i_to_all = _compute_monitored_branch_currents( + theta_all=theta_all, + vm_all=vm_all, + f_pos_safe=f_pos_safe, + t_pos_safe=t_pos_safe, + y_ff_mon=y_ff_mon, + y_ft_mon=y_ft_mon, + y_tf_mon=y_tf_mon, + y_tt_mon=y_tt_mon, + dtype=dtype, + ) + complex_dtype = y_ff_mon.dtype + end_mask_complex = end_mask.astype(complex_dtype)[None, :] + end_mask_real = end_mask.astype(dtype)[None, :] + s_from_all = v_from_all * jnp.conj(i_from_all) * end_mask_complex + s_to_all = v_to_all * jnp.conj(i_to_all) * end_mask_complex + return InjectionOutageMonitoredResults( + n_1_theta=theta_all, + n_1_voltage=vm_all, + n_1_p_from=s_from_all.real.astype(dtype) * end_mask_real, + n_1_p_to=s_to_all.real.astype(dtype) * end_mask_real, + n_1_q_from=s_from_all.imag.astype(dtype) * end_mask_real, + n_1_q_to=s_to_all.imag.astype(dtype) * end_mask_real, + n_1_i_from=i_from_all * end_mask_complex, + n_1_i_to=i_to_all * end_mask_complex, + ) + + +def non_voltage_regulating_injection_outage_dx( + jacobian_inv_transposed: Float[jnp.ndarray, " n_eq n_eq"], + outage_injection_indices: Int[jnp.ndarray, " n_contingencies n_outages"], + injection_to_bus: Int[jnp.ndarray, " n_injections"], + injection_active_power: Float[jnp.ndarray, " n_injections"], + injection_reactive_power: Float[jnp.ndarray, " n_injections"], + angle_component_indices: Int[jnp.ndarray, " n_buses"], + magnitude_component_indices: Int[jnp.ndarray, " n_buses"], +) -> Float[jnp.ndarray, " n_contingencies n_eq"]: + """One-step state increment caused by disconnecting non voltage regulating injections. + + The Jacobian is unchanged for a pure injection outage because network + topology and admittances remain fixed. The outage only changes the nodal + power mismatch at the buses hosting the disconnected injections. + + This helper assumes the bus type pattern used by the Jacobian remains + unchanged across the outage. That is appropriate for load outages and other + contingencies that do not alter voltage-control structure. + + Note: this function still works for voltage regulating injection outages, + but the PV->PQ bus type change is not accounted for. + + Parameters + ---------- + jacobian_inv_transposed : Float[jnp.ndarray, " n_eq n_eq"] + Base transposed inverse Jacobian at the hot-start operating point. + outage_injection_indices : Int[jnp.ndarray, " n_contingencies n_outages"] + Batched indices of injections to disconnect. Each row is one contingency. + injection_to_bus : Int[jnp.ndarray, " n_injections"] + Map from injection index to hosting bus index. + injection_active_power : Float[jnp.ndarray, " n_injections"] + Active power of each injection at the hot-start operating point. + injection_reactive_power : Float[jnp.ndarray, " n_injections"] + Reactive power of each injection at the hot-start operating point. + angle_component_indices : Int[jnp.ndarray, " n_buses"] + Map from bus index to angle component index in the Jacobian. -1 if not present. + magnitude_component_indices : Int[jnp.ndarray, " n_buses"] + Map from bus index to magnitude component index in the Jacobian. -1 if not present. + + Returns + ------- + Float[jnp.ndarray, " n_contingencies n_eq"] + One-step state increment for each outage in Jacobian ordering. + """ + dtype = jacobian_inv_transposed.dtype + mismatch = _collect_injection_outage_mismatch( + jacobian_size=jacobian_inv_transposed.shape[0], + outage_injection_indices=outage_injection_indices, + injection_to_bus=injection_to_bus, + injection_active_power=injection_active_power, + injection_reactive_power=injection_reactive_power, + angle_component_indices=angle_component_indices, + magnitude_component_indices=magnitude_component_indices, + dtype=dtype, + ) + return _calculate_injection_outage_dx(jacobian_inv_transposed=jacobian_inv_transposed, mismatch=mismatch) diff --git a/src/dc_plus/numpy/__init__.py b/src/dc_plus/numpy/__init__.py index c2900d1..dc4b87d 100644 --- a/src/dc_plus/numpy/__init__.py +++ b/src/dc_plus/numpy/__init__.py @@ -11,6 +11,6 @@ It is currently not optimized for performance, but for readability. """ -from .injection_outage import injection_outage_dx +from .injection_outage import non_voltage_regulating_injection_outage_dx -__all__ = ["injection_outage_dx"] +__all__ = ["non_voltage_regulating_injection_outage_dx"] diff --git a/src/dc_plus/numpy/injection_outage.py b/src/dc_plus/numpy/injection_outage.py index 2942f54..666eaa2 100644 --- a/src/dc_plus/numpy/injection_outage.py +++ b/src/dc_plus/numpy/injection_outage.py @@ -11,16 +11,111 @@ from jaxtyping import Float, Int -def injection_outage_dx( +def _collect_injection_outage_mismatch( + jacobian_size: int, + outage_injection_indices: Int[np.ndarray, " n_contingencies n_outages"], + injection_to_bus: Int[np.ndarray, " n_injections"], + injection_active_power: Float[np.ndarray, " n_injections"], + injection_reactive_power: Float[np.ndarray, " n_injections"], + angle_component_indices: Int[np.ndarray, " n_buses"], + magnitude_component_indices: Int[np.ndarray, " n_buses"], + dtype: np.dtype, +) -> Float[np.ndarray, " n_contingencies n_eq"]: + """Assemble outage-induced mismatch vectors in Jacobian ordering. + + Parameters + ---------- + jacobian_size : int + Number of equations in the Jacobian system. + outage_injection_indices : Int[np.ndarray, " n_contingencies n_outages"] + Batched injection outage indices. Each row represents one contingency and + may include negative padding entries. + injection_to_bus : Int[np.ndarray, " n_injections"] + Map from injection index to hosting bus index. + injection_active_power : Float[np.ndarray, " n_injections"] + Active power associated with each injection. + injection_reactive_power : Float[np.ndarray, " n_injections"] + Reactive power associated with each injection. + angle_component_indices : Int[np.ndarray, " n_buses"] + Map from bus index to active-power mismatch component index. + magnitude_component_indices : Int[np.ndarray, " n_buses"] + Map from bus index to reactive-power mismatch component index. + dtype : np.dtype + Target dtype of the assembled mismatch matrix. + + Returns + ------- + Float[np.ndarray, " n_contingencies n_eq"] + Batched mismatch vectors in Jacobian component ordering. + """ + outage_indices = outage_injection_indices + mismatch = np.zeros((outage_indices.shape[0], jacobian_size), dtype=dtype) + if outage_indices.size == 0: + return mismatch + + valid_outage_mask = outage_indices >= 0 + if not np.any(valid_outage_mask): + return mismatch + + valid_rows, valid_cols = np.nonzero(valid_outage_mask) + valid_outages = outage_indices[valid_rows, valid_cols] + if np.any(valid_outages >= injection_to_bus.size): + raise IndexError("Injection outage index is out of bounds") + + buses = injection_to_bus[valid_outages] + + p_idx = angle_component_indices[buses] + valid_p = p_idx >= 0 + if np.any(valid_p): + np.add.at( + mismatch, + (valid_rows[valid_p], p_idx[valid_p]), + -injection_active_power[valid_outages][valid_p], + ) + + q_idx = magnitude_component_indices[buses] + valid_q = q_idx >= 0 + if np.any(valid_q): + np.add.at( + mismatch, + (valid_rows[valid_q], q_idx[valid_q]), + -injection_reactive_power[valid_outages][valid_q], + ) + + return mismatch + + +def _calculate_injection_outage_dx( + jacobian_inv: Float[np.ndarray, " n_eq n_eq"], + mismatch: Float[np.ndarray, " n_contingencies n_eq"], +) -> Float[np.ndarray, " n_contingencies n_eq"]: + """Map outage mismatch vectors to one-step state increments. + + Parameters + ---------- + jacobian_inv : Float[np.ndarray, " n_eq n_eq"] + Inverse Jacobian at the base operating point. + mismatch : Float[np.ndarray, " n_contingencies n_eq"] + Batched outage-induced mismatch vectors. + + Returns + ------- + Float[np.ndarray, " n_contingencies n_eq"] + Batched one-step state increments in Jacobian ordering. + """ + return -(mismatch @ jacobian_inv.T) + + +def non_voltage_regulating_injection_outage_dx( jacobian_inv: Float[np.ndarray, " n_eq n_eq"], - outage_injection_indices: Int[np.ndarray, " n_outages"], + outage_injection_indices: Int[np.ndarray, " n_contingencies n_outages"], injection_to_bus: Int[np.ndarray, " n_injections"], injection_active_power: Float[np.ndarray, " n_injections"], injection_reactive_power: Float[np.ndarray, " n_injections"], angle_component_indices: Int[np.ndarray, " n_buses"], magnitude_component_indices: Int[np.ndarray, " n_buses"], -) -> Float[np.ndarray, " n_eq"]: - """Return the one-step state increment caused by disconnecting injections. +) -> Float[np.ndarray, " n_contingencies n_eq"]: + """One-step state increment caused by disconnecting non voltage regulating injections. The Jacobian is unchanged for a pure injection outage because network topology and admittances remain fixed. The outage only changes the nodal @@ -30,12 +125,18 @@ def injection_outage_dx( unchanged across the outage. That is appropriate for load outages and other contingencies that do not alter voltage-control structure. + Note: this function still works for voltage regulating injection outages, + but the PV->PQ bus type change is not accounted for. + Parameters ---------- jacobian_inv : Float[np.ndarray, " n_eq n_eq"] Base inverse Jacobian at the hot-start operating point. - outage_injection_indices : Int[np.ndarray, " n_outages"] - Indices of injections to disconnect. + outage_injection_indices : Int[np.ndarray, " n_contingencies n_outages"] + Batched indices of injections to disconnect. Each row is one + contingency and may contain multiple outages. Negative entries are + ignored and can be used as padding for rows with different outage + counts. injection_to_bus : Int[np.ndarray, " n_injections"] Mapping from injection index to bus index. injection_active_power : Float[np.ndarray, " n_injections"] @@ -53,27 +154,17 @@ def injection_outage_dx( Returns ------- - Float[np.ndarray, " n_eq"] - One-step update vector in Jacobian ordering. + Float[np.ndarray, " n_contingencies n_eq"] + Batched one-step update vectors in Jacobian ordering. """ - outage_indices = outage_injection_indices - mismatch = np.zeros(jacobian_inv.shape[0], dtype=jacobian_inv.dtype) - if outage_indices.size == 0: - return mismatch - - if np.any((outage_indices < 0) | (outage_indices >= injection_to_bus.size)): - raise IndexError("Injection outage index is out of bounds") - - buses = injection_to_bus[outage_indices] - - p_idx = angle_component_indices[buses] - valid_p = p_idx >= 0 - if np.any(valid_p): - np.add.at(mismatch, p_idx[valid_p], -injection_active_power[outage_indices][valid_p]) - - q_idx = magnitude_component_indices[buses] - valid_q = q_idx >= 0 - if np.any(valid_q): - np.add.at(mismatch, q_idx[valid_q], -injection_reactive_power[outage_indices][valid_q]) - - return -jacobian_inv @ mismatch + mismatch = _collect_injection_outage_mismatch( + jacobian_size=jacobian_inv.shape[0], + outage_injection_indices=outage_injection_indices, + injection_to_bus=injection_to_bus, + injection_active_power=injection_active_power, + injection_reactive_power=injection_reactive_power, + angle_component_indices=angle_component_indices, + magnitude_component_indices=magnitude_component_indices, + dtype=jacobian_inv.dtype, + ) + return _calculate_injection_outage_dx(jacobian_inv=jacobian_inv, mismatch=mismatch) diff --git a/tests/jax/test_injection_outage_jax.py b/tests/jax/test_injection_outage_jax.py new file mode 100644 index 0000000..cc0a021 --- /dev/null +++ b/tests/jax/test_injection_outage_jax.py @@ -0,0 +1,241 @@ +# Copyright 2026 50Hertz Transmission GmbH and Elia Transmission Belgium SA/NV +# +# This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0. +# If a copy of the MPL was not distributed with this file, +# you can obtain one at https://mozilla.org/MPL/2.0/. +# Mozilla Public License, version 2.0 + +import jax +import jax.numpy as jnp +import numpy as np +import pytest + +jax.config.update("jax_enable_x64", True) +jax.config.update("jax_platform_name", "cpu") + +from dc_plus.example_grids.pypowsbl.example_grids import PANDAPOWER_NETWORKS_FOR_POWSYBL, POWSYBL_NETWORKS +from dc_plus.importing.powsybl.powsybl_import import _get_injections_powsybl +from dc_plus.importing.powsybl.powsybl_network_helpers import _load_test_grid +from dc_plus.jax.injection_outage import ( + non_voltage_regulating_injection_outage_monitor_buses, + non_voltage_regulating_injection_outage_monitored, +) +from dc_plus.jax.injection_outage import non_voltage_regulating_injection_outage_dx as jax_injection_outage_dx +from dc_plus.numpy.injection_outage import non_voltage_regulating_injection_outage_dx as numpy_injection_outage_dx + +powsybl_networks = POWSYBL_NETWORKS +pandapower_networks = PANDAPOWER_NETWORKS_FOR_POWSYBL +TOL = 1e-10 + + +def _connected_injection_indices(injections): + candidates = np.flatnonzero(injections["connected"].to_numpy(dtype=bool)) + if candidates.size == 0: + pytest.skip("No connected injections available for an injection outage test") + return candidates + + +def _supported_injection_outage_rows(net, injections, max_rows: int | None = None) -> list[np.ndarray]: + injection_powsybl_type = net.get_injections(attributes=["type"]) + outage_rows: list[np.ndarray] = [] + for injection_idx in _connected_injection_indices(injections): + injection_id = injections.loc[injection_idx, "id_str"] + if (injection_powsybl_type.loc[injection_id, "type"]) in (["HVDC_CONVERTER_STATION", "DANGLING_LINE"]): + continue + outage_rows.append(np.array([injection_idx], dtype=np.int64)) + if max_rows is not None and len(outage_rows) >= max_rows: + break + return outage_rows + + +def _pad_outage_rows(outage_rows: list[np.ndarray]) -> np.ndarray: + if not outage_rows: + return np.empty((0, 0), dtype=np.int64) + + max_outages = max(outage_row.size for outage_row in outage_rows) + padded = np.full((len(outage_rows), max_outages), -1, dtype=np.int64) + for row_idx, outage_row in enumerate(outage_rows): + padded[row_idx, : outage_row.size] = outage_row + return padded + + +def test_non_voltage_regulating_injection_outage_jax_batch_multi_outage() -> None: + jacobian_inv = np.eye(4, dtype=float) + outage_batch = np.array([[0, 2], [1, -1], [-1, -1]], dtype=np.int64) + injection_to_bus = np.array([0, 1, 0], dtype=np.int64) + injection_active_power = np.array([1.0, 2.0, 3.0], dtype=float) + injection_reactive_power = np.array([10.0, 20.0, 30.0], dtype=float) + angle_component_indices = np.array([0, 1], dtype=np.int64) + magnitude_component_indices = np.array([2, 3], dtype=np.int64) + + dx_batch_jax = np.asarray( + jax_injection_outage_dx( + jacobian_inv_transposed=jnp.asarray(jacobian_inv.T), + outage_injection_indices=jnp.asarray(outage_batch), + injection_to_bus=jnp.asarray(injection_to_bus), + injection_active_power=jnp.asarray(injection_active_power), + injection_reactive_power=jnp.asarray(injection_reactive_power), + angle_component_indices=jnp.asarray(angle_component_indices), + magnitude_component_indices=jnp.asarray(magnitude_component_indices), + ) + ) + dx_batch_numpy = numpy_injection_outage_dx( + jacobian_inv=jacobian_inv, + outage_injection_indices=outage_batch, + injection_to_bus=injection_to_bus, + injection_active_power=injection_active_power, + injection_reactive_power=injection_reactive_power, + angle_component_indices=angle_component_indices, + magnitude_component_indices=magnitude_component_indices, + ) + + np.testing.assert_allclose(dx_batch_jax, dx_batch_numpy, rtol=TOL, atol=TOL) + + +@pytest.mark.parametrize("get_net", powsybl_networks + pandapower_networks) +def test_non_voltage_regulating_injection_outage_jax(get_net) -> None: + net, _, dynamic_info, _, jacobian_data = _load_test_grid(get_net) + injections = _get_injections_powsybl(net).reset_index(drop=True) + outage_p = injections["setpoint_p"].fillna(injections["p"]).to_numpy(dtype=float) + outage_q = injections["setpoint_q"].fillna(injections["q"]).to_numpy(dtype=float) + + assert len(injections) == dynamic_info.n_injections + + outage_rows = _supported_injection_outage_rows(net, injections) + assert outage_rows, f"No supported injection outages found for {get_net.__name__}" + + outage_batch = _pad_outage_rows(outage_rows) + dx_batch_jax = np.asarray( + jax_injection_outage_dx( + jacobian_inv_transposed=jnp.asarray(jacobian_data.inverse_jacobian.T), + outage_injection_indices=jnp.asarray(outage_batch), + injection_to_bus=jnp.asarray(dynamic_info.injection_to_bus), + injection_active_power=jnp.asarray(outage_p), + injection_reactive_power=jnp.asarray(outage_q), + angle_component_indices=jnp.asarray(jacobian_data.angle_component_indices), + magnitude_component_indices=jnp.asarray(jacobian_data.magnitude_component_indices), + ) + ) + dx_batch_numpy = numpy_injection_outage_dx( + jacobian_inv=jacobian_data.inverse_jacobian, + outage_injection_indices=outage_batch, + injection_to_bus=dynamic_info.injection_to_bus, + injection_active_power=outage_p, + injection_reactive_power=outage_q, + angle_component_indices=jacobian_data.angle_component_indices, + magnitude_component_indices=jacobian_data.magnitude_component_indices, + ) + + assert dx_batch_jax.shape == (len(outage_rows), jacobian_data.inverse_jacobian.shape[0]) + # the numpy code is tested against powsybl -> we do not need to repeat the test against powsybl + np.testing.assert_allclose(dx_batch_jax, dx_batch_numpy, rtol=TOL, atol=TOL) + + +@pytest.mark.parametrize("get_net", powsybl_networks + pandapower_networks) +def test_non_voltage_regulating_injection_outage_monitored_jax(get_net) -> None: + net, _, dynamic_info, _, jacobian_data = _load_test_grid(get_net) + injections = _get_injections_powsybl(net).reset_index(drop=True) + outage_p = injections["setpoint_p"].fillna(injections["p"]).to_numpy(dtype=float) + outage_q = injections["setpoint_q"].fillna(injections["q"]).to_numpy(dtype=float) + outage_rows = _supported_injection_outage_rows(net, injections, max_rows=5) + assert outage_rows, f"No supported injection outages found for {get_net.__name__}" + outage_batch = _pad_outage_rows(outage_rows) + + branch_connected = np.asarray(dynamic_info.branch_connected, dtype=bool) + monitor_branch_indices = np.flatnonzero(branch_connected)[: min(5, int(branch_connected.sum()))] + if monitor_branch_indices.size == 0: + pytest.skip("No connected monitored branches available") + + branch_from = np.asarray(dynamic_info.branch_from_bus, dtype=np.int64) + branch_to = np.asarray(dynamic_info.branch_to_bus, dtype=np.int64) + monitor_bus_indices = np.unique( + np.concatenate([branch_from[monitor_branch_indices], branch_to[monitor_branch_indices]]) + ).astype(np.int64) + bus_to_mon_index = np.full(dynamic_info.n_buses, -1, dtype=np.int64) + bus_to_mon_index[monitor_bus_indices] = np.arange(monitor_bus_indices.size, dtype=np.int64) + + monitor_theta_jax, monitor_vm_jax = non_voltage_regulating_injection_outage_monitor_buses( + jacobian_inv_transposed=jnp.asarray(jacobian_data.inverse_jacobian.T), + outage_injection_indices=jnp.asarray(outage_batch), + injection_to_bus=jnp.asarray(dynamic_info.injection_to_bus), + injection_active_power=jnp.asarray(outage_p), + injection_reactive_power=jnp.asarray(outage_q), + angle_component_indices=jnp.asarray(jacobian_data.angle_component_indices), + magnitude_component_indices=jnp.asarray(jacobian_data.magnitude_component_indices), + monitor_bus_indices=jnp.asarray(monitor_bus_indices), + v_mag_hat=jnp.asarray(dynamic_info.bus_voltage_magnitudes), + theta_hat=jnp.asarray(dynamic_info.bus_voltage_angles_rad), + ) + monitored_results = non_voltage_regulating_injection_outage_monitored( + jacobian_inv_transposed=jnp.asarray(jacobian_data.inverse_jacobian.T), + outage_injection_indices=jnp.asarray(outage_batch), + injection_to_bus=jnp.asarray(dynamic_info.injection_to_bus), + injection_active_power=jnp.asarray(outage_p), + injection_reactive_power=jnp.asarray(outage_q), + angle_component_indices=jnp.asarray(jacobian_data.angle_component_indices), + magnitude_component_indices=jnp.asarray(jacobian_data.magnitude_component_indices), + monitor_bus_indices=jnp.asarray(monitor_bus_indices), + v_mag_hat=jnp.asarray(dynamic_info.bus_voltage_magnitudes), + theta_hat=jnp.asarray(dynamic_info.bus_voltage_angles_rad), + branch_from=jnp.asarray(branch_from), + branch_to=jnp.asarray(branch_to), + y_ff=jnp.asarray(dynamic_info.branch_effective_admittance_from_from), + y_ft=jnp.asarray(dynamic_info.branch_effective_admittance_from_to), + y_tf=jnp.asarray(dynamic_info.branch_effective_admittance_to_from), + y_tt=jnp.asarray(dynamic_info.branch_effective_admittance_to_to), + monitor_branch_indices=jnp.asarray(monitor_branch_indices), + bus_to_mon_index=jnp.asarray(bus_to_mon_index), + ) + dx_batch_numpy = numpy_injection_outage_dx( + jacobian_inv=jacobian_data.inverse_jacobian, + outage_injection_indices=outage_batch, + injection_to_bus=dynamic_info.injection_to_bus, + injection_active_power=outage_p, + injection_reactive_power=outage_q, + angle_component_indices=jacobian_data.angle_component_indices, + magnitude_component_indices=jacobian_data.magnitude_component_indices, + ) + + expected_theta = np.broadcast_to( + np.asarray(dynamic_info.bus_voltage_angles_rad, dtype=float)[monitor_bus_indices], + (len(outage_rows), monitor_bus_indices.size), + ).copy() + expected_vm = np.broadcast_to( + np.asarray(dynamic_info.bus_voltage_magnitudes, dtype=float)[monitor_bus_indices], + (len(outage_rows), monitor_bus_indices.size), + ).copy() + theta_component_indices = np.asarray(jacobian_data.angle_component_indices, dtype=np.int64)[monitor_bus_indices] + vm_component_indices = np.asarray(jacobian_data.magnitude_component_indices, dtype=np.int64)[monitor_bus_indices] + theta_mask = theta_component_indices >= 0 + vm_mask = vm_component_indices >= 0 + expected_theta[:, theta_mask] += dx_batch_numpy[:, theta_component_indices[theta_mask]] + expected_vm[:, vm_mask] += dx_batch_numpy[:, vm_component_indices[vm_mask]] + + np.testing.assert_allclose(np.asarray(monitor_theta_jax), expected_theta, rtol=TOL, atol=TOL) + np.testing.assert_allclose(np.asarray(monitor_vm_jax), expected_vm, rtol=TOL, atol=TOL) + np.testing.assert_allclose(np.asarray(monitored_results.n_1_theta), expected_theta, rtol=TOL, atol=TOL) + np.testing.assert_allclose(np.asarray(monitored_results.n_1_voltage), expected_vm, rtol=TOL, atol=TOL) + + f_pos = bus_to_mon_index[branch_from[monitor_branch_indices]] + t_pos = bus_to_mon_index[branch_to[monitor_branch_indices]] + assert np.all(f_pos >= 0) + assert np.all(t_pos >= 0) + + v_post = expected_vm * np.exp(1j * expected_theta) + y_ff_mon = np.asarray(dynamic_info.branch_effective_admittance_from_from)[monitor_branch_indices] + y_ft_mon = np.asarray(dynamic_info.branch_effective_admittance_from_to)[monitor_branch_indices] + y_tf_mon = np.asarray(dynamic_info.branch_effective_admittance_to_from)[monitor_branch_indices] + y_tt_mon = np.asarray(dynamic_info.branch_effective_admittance_to_to)[monitor_branch_indices] + v_from = v_post[:, f_pos] + v_to = v_post[:, t_pos] + i_from = y_ff_mon[None, :] * v_from + y_ft_mon[None, :] * v_to + i_to = y_tf_mon[None, :] * v_from + y_tt_mon[None, :] * v_to + s_from = v_from * np.conj(i_from) + s_to = v_to * np.conj(i_to) + + np.testing.assert_allclose(np.asarray(monitored_results.n_1_i_from), i_from, rtol=1e-9, atol=1e-9) + np.testing.assert_allclose(np.asarray(monitored_results.n_1_i_to), i_to, rtol=1e-9, atol=1e-9) + np.testing.assert_allclose(np.asarray(monitored_results.n_1_p_from), s_from.real, rtol=1e-9, atol=1e-9) + np.testing.assert_allclose(np.asarray(monitored_results.n_1_p_to), s_to.real, rtol=1e-9, atol=1e-9) + np.testing.assert_allclose(np.asarray(monitored_results.n_1_q_from), s_from.imag, rtol=1e-9, atol=1e-9) + np.testing.assert_allclose(np.asarray(monitored_results.n_1_q_to), s_to.imag, rtol=1e-9, atol=1e-9) diff --git a/tests/numpy/test_injection_outage.py b/tests/numpy/test_injection_outage.py index e205db0..7c758dc 100644 --- a/tests/numpy/test_injection_outage.py +++ b/tests/numpy/test_injection_outage.py @@ -12,16 +12,15 @@ import pytest from dc_plus.example_grids.pypowsbl.example_grids import PANDAPOWER_NETWORKS_FOR_POWSYBL, POWSYBL_NETWORKS -from dc_plus.importing.powsybl.powsybl_import import _get_injections_powsybl +from dc_plus.importing.powsybl.powsybl_import import DANGLING_BUS_STRING_SUFFIX, _get_injections_powsybl from dc_plus.importing.powsybl.powsybl_loadflow_parameter import get_powsybl_loadflow_parameter from dc_plus.importing.powsybl.powsybl_network_helpers import _load_test_grid -from dc_plus.interfaces.network_information import BusType -from dc_plus.numpy.injection_outage import injection_outage_dx +from dc_plus.numpy.injection_outage import non_voltage_regulating_injection_outage_dx from dc_plus.preprocess.create_network_data import create_network_data powsybl_networks = POWSYBL_NETWORKS pandapower_networks = PANDAPOWER_NETWORKS_FOR_POWSYBL -TOL = 1e-10 +TOL = 5e-10 def _connected_injection_indices(injections): @@ -31,30 +30,66 @@ def _connected_injection_indices(injections): return candidates -def _expected_bus_delta(dynamic_info, injection_idx, outage_p, outage_q): +def _expected_bus_delta(dynamic_info, outage_injection_indices, outage_p, outage_q): delta_p = np.zeros(dynamic_info.n_buses, dtype=float) delta_q = np.zeros(dynamic_info.n_buses, dtype=float) - bus_idx = int(dynamic_info.injection_to_bus[injection_idx]) - delta_p[bus_idx] -= float(outage_p[injection_idx]) - delta_q[bus_idx] -= float(outage_q[injection_idx]) + outage_indices = np.asarray(outage_injection_indices, dtype=np.int64) + valid_mask = outage_indices >= 0 + if not np.any(valid_mask): + return delta_p, delta_q + + valid_outages = outage_indices[valid_mask] + buses = dynamic_info.injection_to_bus[valid_outages] + np.add.at(delta_p, buses, -outage_p[valid_outages]) + np.add.at(delta_q, buses, -outage_q[valid_outages]) return delta_p, delta_q +def _pad_outage_rows(outage_rows: list[np.ndarray]) -> np.ndarray: + if not outage_rows: + return np.empty((0, 0), dtype=np.int64) + + max_outages = max(outage_row.size for outage_row in outage_rows) + padded = np.full((len(outage_rows), max_outages), -1, dtype=np.int64) + for row_idx, outage_row in enumerate(outage_rows): + padded[row_idx, : outage_row.size] = outage_row + return padded + + @pytest.mark.parametrize("get_net", powsybl_networks + pandapower_networks) -def test_injection_outage_numpy_compare_powsybl(get_net): +def test_non_voltage_regulating_injection_outage_numpy(get_net) -> None: + """Test non voltage regulating injection outages with a one-step update against powsybl load flow results. + + This does not test HVDC converter and dangling line outages. + It does test voltage regulatingGenerator outages, as the powsybl one step does also not switch from PV to PQ + node. This therefore does not account for the voltage drop of the disconnected generator. + + In DC+ dangling lines are modelled as an additional bus with a fixed injection p0 and q0 + a one sided line is added between the original bus and the new bus to account for the line losses + of the dangling line. Here numerical errors occur. + The test still verifies the voltages of the bus where the dangling line is attached to with the + full precision. + + Parameters + ---------- + get_net : Callable[[], pypowsybl.network.Network] + A fixture that returns a powsybl network to test on. + """ net, _, dynamic_info, string_info, jacobian_data = _load_test_grid(get_net) injections = _get_injections_powsybl(net).reset_index(drop=True) injection_powsybl_type = net.get_injections(attributes=["type"]) outage_p = injections["setpoint_p"].fillna(injections["p"]).to_numpy(dtype=float) outage_q = injections["setpoint_q"].fillna(injections["q"]).to_numpy(dtype=float) + dangling_bus_mask = np.char.endswith(string_info.bus_ids.astype(str), DANGLING_BUS_STRING_SUFFIX) assert len(injections) == dynamic_info.n_injections - compared_cases = 0 + accepted_cases: list[tuple[str, np.ndarray, object]] = [] for injection_idx in _connected_injection_indices(injections): injection_id = injections.loc[injection_idx, "id_str"] - if (injection_powsybl_type.loc[injection_id, "type"]) in (["HVDC_CONVERTER_STATION"]): - # HVDC would need a multi outage + if (injection_powsybl_type.loc[injection_id, "type"]) in (["HVDC_CONVERTER_STATION", "DANGLING_LINE"]): + # - HVDC would need a multi outage + # - Dangling line removal would change the network for DC+ continue net_n1 = deepcopy(net) @@ -67,7 +102,8 @@ def test_injection_outage_numpy_compare_powsybl(get_net): continue _, dynamic_info_n1, string_info_n1 = create_network_data(net_n1) - np.testing.assert_array_equal(string_info.bus_ids, string_info_n1.bus_ids) + if not np.array_equal(string_info.bus_ids, string_info_n1.bus_ids): + continue same_angle_structure = np.array_equal( dynamic_info.pvpq_buses_indices_pvpq_order, @@ -80,7 +116,8 @@ def test_injection_outage_numpy_compare_powsybl(get_net): if not (same_angle_structure and same_magnitude_structure): continue - expected_delta_p, expected_delta_q = _expected_bus_delta(dynamic_info, injection_idx, outage_p, outage_q) + outage_row = np.array([injection_idx], dtype=np.int64) + expected_delta_p, expected_delta_q = _expected_bus_delta(dynamic_info, outage_row, outage_p, outage_q) actual_delta_p = dynamic_info_n1.bus_active_power - dynamic_info.bus_active_power actual_delta_q = dynamic_info_n1.bus_reactive_power - dynamic_info.bus_reactive_power p_buses = dynamic_info.pvpq_buses_indices_pvpq_order @@ -91,32 +128,38 @@ def test_injection_outage_numpy_compare_powsybl(get_net): ): continue - dx = injection_outage_dx( - jacobian_inv=jacobian_data.inverse_jacobian, - outage_injection_indices=np.array([injection_idx], dtype=np.int64), - injection_to_bus=dynamic_info.injection_to_bus, - injection_active_power=outage_p, - injection_reactive_power=outage_q, - angle_component_indices=jacobian_data.angle_component_indices, - magnitude_component_indices=jacobian_data.magnitude_component_indices, - ) + accepted_cases.append((injection_id, outage_row, dynamic_info_n1)) + assert accepted_cases, f"No injection outages matched the fixed-injection assumptions for {get_net.__name__}" + + outage_batch = _pad_outage_rows([outage_row for _, outage_row, _ in accepted_cases]) + dx_batch = non_voltage_regulating_injection_outage_dx( + jacobian_inv=jacobian_data.inverse_jacobian, + outage_injection_indices=outage_batch, + injection_to_bus=dynamic_info.injection_to_bus, + injection_active_power=outage_p, + injection_reactive_power=outage_q, + angle_component_indices=jacobian_data.angle_component_indices, + magnitude_component_indices=jacobian_data.magnitude_component_indices, + ) + + assert dx_batch.shape == (len(accepted_cases), jacobian_data.inverse_jacobian.shape[0]) + + for case_idx, (_, _, dynamic_info_n1) in enumerate(accepted_cases): + dx = dx_batch[case_idx] theta_updated = dynamic_info.bus_voltage_angles_rad.copy() vm_updated = dynamic_info.bus_voltage_magnitudes.copy() theta_updated[dynamic_info.pvpq_buses_indices_pvpq_order] += dx[jacobian_data.is_angle_component] vm_updated[dynamic_info.pq_buses_indices] += dx[jacobian_data.is_magnitude_component] np.testing.assert_allclose( - theta_updated, - dynamic_info_n1.bus_voltage_angles_rad, + theta_updated[~dangling_bus_mask], + dynamic_info_n1.bus_voltage_angles_rad[~dangling_bus_mask], rtol=TOL, atol=TOL, ) np.testing.assert_allclose( - vm_updated, - dynamic_info_n1.bus_voltage_magnitudes, + vm_updated[~dangling_bus_mask], + dynamic_info_n1.bus_voltage_magnitudes[~dangling_bus_mask], rtol=TOL, atol=TOL, ) - compared_cases += 1 - - assert compared_cases > 0, f"No injection outages matched the fixed-injection assumptions for {get_net.__name__}" From 214938aa971d1845eb07d58ba8d0a32fe638adf0 Mon Sep 17 00:00:00 2001 From: Bpetrick <170433522+BenjPetr@users.noreply.github.com> Date: Sat, 18 Apr 2026 03:10:52 +0000 Subject: [PATCH 4/7] doc: refine text --- docs/dc_plus/injection_outage.md | 2 +- tests/numpy/test_injection_outage.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/dc_plus/injection_outage.md b/docs/dc_plus/injection_outage.md index 0d4e5e4..344ebbb 100644 --- a/docs/dc_plus/injection_outage.md +++ b/docs/dc_plus/injection_outage.md @@ -29,4 +29,4 @@ The function can also be applied to voltage-regulating generator outages, but th ## HVDC -- HVDC converter outages should be outaged in pairs. DC+ models the HVDC line as a pair of generator + load +- HVDC converter outages should be outaged in pairs. DC+ models the HVDC line as a pair of generator + load injections, one at each terminal. diff --git a/tests/numpy/test_injection_outage.py b/tests/numpy/test_injection_outage.py index 7c758dc..4899ba3 100644 --- a/tests/numpy/test_injection_outage.py +++ b/tests/numpy/test_injection_outage.py @@ -61,7 +61,7 @@ def test_non_voltage_regulating_injection_outage_numpy(get_net) -> None: """Test non voltage regulating injection outages with a one-step update against powsybl load flow results. This does not test HVDC converter and dangling line outages. - It does test voltage regulatingGenerator outages, as the powsybl one step does also not switch from PV to PQ + It does test voltage regulating generator outages, as the powsybl one step does also not switch from PV to PQ node. This therefore does not account for the voltage drop of the disconnected generator. In DC+ dangling lines are modelled as an additional bus with a fixed injection p0 and q0 From af4dc4fd4bcd35d0f1bafdc5c6515d5dce6dd8d0 Mon Sep 17 00:00:00 2001 From: Bpetrick <170433522+BenjPetr@users.noreply.github.com> Date: Sat, 18 Apr 2026 03:56:58 +0000 Subject: [PATCH 5/7] refactor: harden reference bus powsybl --- src/dc_plus/importing/powsybl/powsybl_import.py | 6 ++++-- src/dc_plus/interfaces/jacobian_network_data.py | 2 +- src/dc_plus/preprocess/create_network_data.py | 7 ++++++- tests/importing/powsybl/test_powsybl_import.py | 2 +- 4 files changed, 12 insertions(+), 5 deletions(-) diff --git a/src/dc_plus/importing/powsybl/powsybl_import.py b/src/dc_plus/importing/powsybl/powsybl_import.py index b2b1dc9..7768b65 100644 --- a/src/dc_plus/importing/powsybl/powsybl_import.py +++ b/src/dc_plus/importing/powsybl/powsybl_import.py @@ -928,7 +928,7 @@ def _get_dangling_buses(net: Network) -> BusParamSchema: return dangling_buses -def _get_buses_powsybl(net: Network, slack_id: str, injections: InjectionParamSchema) -> BusParamSchema: +def _get_buses_powsybl(net: Network, slack_id: str, injections: InjectionParamSchema, reference_id: str) -> BusParamSchema: """Get the bus parameters of the network. Gets the bus parameters from the network. @@ -941,6 +941,8 @@ def _get_buses_powsybl(net: Network, slack_id: str, injections: InjectionParamSc The id of the slack bus. injections : InjectionParamSchema The injections of the network. + reference_id : str + The id of the reference bus for angle reference. Returns ------- @@ -953,7 +955,7 @@ def _get_buses_powsybl(net: Network, slack_id: str, injections: InjectionParamSc buses.reset_index(drop=False, inplace=True) buses.rename(columns={"id": "id_str"}, inplace=True) buses["bus_type"] = 2 # PQ bus is default - buses["is_angle_reference"] = buses["v_angle"] == 0.0 + buses["is_angle_reference"] = buses["id_str"].eq(reference_id) buses.rename(columns={"v_mag": "voltage_magnitude", "v_angle": "voltage_angle"}, inplace=True) diff --git a/src/dc_plus/interfaces/jacobian_network_data.py b/src/dc_plus/interfaces/jacobian_network_data.py index db6faf5..c784196 100644 --- a/src/dc_plus/interfaces/jacobian_network_data.py +++ b/src/dc_plus/interfaces/jacobian_network_data.py @@ -86,7 +86,7 @@ def get_jacobian_from_network_data( pqpv_indices = dynamic_network_data.pvpq_buses_indices_pvpq_order voltage = voltage_magnitudes * np.exp(1.0j * voltage_angles) - if np.any(voltage == 0.0): + if np.any(np.isclose(voltage, 0.0, atol=1e-10)): raise ValueError("Voltage magnitudes must be strictly positive to construct the Jacobian.") voltage_norm = voltage / abs(voltage) diff --git a/src/dc_plus/preprocess/create_network_data.py b/src/dc_plus/preprocess/create_network_data.py index 2df1dd3..9efafd0 100644 --- a/src/dc_plus/preprocess/create_network_data.py +++ b/src/dc_plus/preprocess/create_network_data.py @@ -188,7 +188,12 @@ def create_network_data( injections = _get_injections_powsybl(network) shunts = _get_shunts_powsybl(network) slack_id = network.get_extensions("slackTerminal")["bus_id"].values[0] - buses = _get_buses_powsybl(net=network, slack_id=slack_id, injections=injections) + references = network.get_extensions("referencePriorities") + if len(references) > 0: + reference_id = network.get_extensions("referencePriorities").index[0] + else: + reference_id = slack_id + buses = _get_buses_powsybl(net=network, slack_id=slack_id, injections=injections, reference_id=reference_id) limits = _get_limits_parameter_powsybl(network) return _create_network_data( diff --git a/tests/importing/powsybl/test_powsybl_import.py b/tests/importing/powsybl/test_powsybl_import.py index 5c66694..f8cfa51 100644 --- a/tests/importing/powsybl/test_powsybl_import.py +++ b/tests/importing/powsybl/test_powsybl_import.py @@ -190,7 +190,7 @@ def test_get_buses(): net = create_complex_grid_battery_hvdc_svc_3w_trafo() injections = _get_injections_powsybl(net) slack_id = net.get_extensions("slackTerminal")["bus_id"].values[0] - buses = _get_buses_powsybl(net=net, slack_id=slack_id, injections=injections) + buses = _get_buses_powsybl(net=net, slack_id=slack_id, injections=injections, reference_id=slack_id) BusParamSchema.validate(buses) bus_powsybl = net.get_buses(attributes=[]) dangling_count = len(_get_dangling_bus_ids(net)) From f19cb90e23d8f6d776b592f77c57668cd82ae23c Mon Sep 17 00:00:00 2001 From: Bpetrick <170433522+BenjPetr@users.noreply.github.com> Date: Sat, 18 Apr 2026 06:28:50 +0000 Subject: [PATCH 6/7] test: refactor injection tests to reduce runtime --- src/dc_plus/example_grids/pypowsbl/example_grids.py | 9 +++++++++ tests/numpy/test_injection_outage.py | 9 ++++----- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/src/dc_plus/example_grids/pypowsbl/example_grids.py b/src/dc_plus/example_grids/pypowsbl/example_grids.py index 2aae34b..5cbd517 100644 --- a/src/dc_plus/example_grids/pypowsbl/example_grids.py +++ b/src/dc_plus/example_grids/pypowsbl/example_grids.py @@ -983,6 +983,15 @@ def basic_node_breaker_network_powsybl() -> pypowsybl.network.Network: pypowsybl.network.create_metrix_tutorial_six_buses_network, # HVDC ] +POWSYBL_NETWORKS_SHORT_LIST = [ + pypowsybl.network.create_ieee9, + pypowsybl.network.create_ieee14, + pypowsybl.network.create_ieee57, + pypowsybl.network.create_micro_grid_nl_network, + basic_node_breaker_network_powsybl, + create_complex_grid_battery_hvdc_svc_3w_trafo, +] + # TODO: check this network and add to POWSYBL_NETWORKS if possible (currently fails) POWSYBL_NETWORKS_NOT_IMPLEMENTED = [ pypowsybl.network.create_micro_grid_be_network, diff --git a/tests/numpy/test_injection_outage.py b/tests/numpy/test_injection_outage.py index 4899ba3..054bb21 100644 --- a/tests/numpy/test_injection_outage.py +++ b/tests/numpy/test_injection_outage.py @@ -11,16 +11,15 @@ import pypowsybl import pytest -from dc_plus.example_grids.pypowsbl.example_grids import PANDAPOWER_NETWORKS_FOR_POWSYBL, POWSYBL_NETWORKS +from dc_plus.example_grids.pypowsbl.example_grids import POWSYBL_NETWORKS_SHORT_LIST from dc_plus.importing.powsybl.powsybl_import import DANGLING_BUS_STRING_SUFFIX, _get_injections_powsybl from dc_plus.importing.powsybl.powsybl_loadflow_parameter import get_powsybl_loadflow_parameter from dc_plus.importing.powsybl.powsybl_network_helpers import _load_test_grid from dc_plus.numpy.injection_outage import non_voltage_regulating_injection_outage_dx from dc_plus.preprocess.create_network_data import create_network_data -powsybl_networks = POWSYBL_NETWORKS -pandapower_networks = PANDAPOWER_NETWORKS_FOR_POWSYBL -TOL = 5e-10 +powsybl_networks = POWSYBL_NETWORKS_SHORT_LIST +TOL = 1e-10 def _connected_injection_indices(injections): @@ -56,7 +55,7 @@ def _pad_outage_rows(outage_rows: list[np.ndarray]) -> np.ndarray: return padded -@pytest.mark.parametrize("get_net", powsybl_networks + pandapower_networks) +@pytest.mark.parametrize("get_net", powsybl_networks) def test_non_voltage_regulating_injection_outage_numpy(get_net) -> None: """Test non voltage regulating injection outages with a one-step update against powsybl load flow results. From 727a9f9b934e1c9b49daee50c620ace4735f1db3 Mon Sep 17 00:00:00 2001 From: Bpetrick <170433522+BenjPetr@users.noreply.github.com> Date: Sat, 18 Apr 2026 07:18:37 +0000 Subject: [PATCH 7/7] fix: fixes warnings and hardens code --- .../importing/powsybl/powsybl_import_helpers.py | 4 ++-- .../importing/powsybl/powsybl_loadflow_parameter.py | 4 ++-- src/dc_plus/interfaces/jacobian_network_data.py | 2 +- src/dc_plus/jax/bsdf.py | 11 ++++++----- src/dc_plus/jax/low_rank_helper.py | 6 ++++-- src/dc_plus/numpy/low_rank_helper.py | 2 +- 6 files changed, 16 insertions(+), 13 deletions(-) diff --git a/src/dc_plus/importing/powsybl/powsybl_import_helpers.py b/src/dc_plus/importing/powsybl/powsybl_import_helpers.py index 8e166bf..e3498b0 100644 --- a/src/dc_plus/importing/powsybl/powsybl_import_helpers.py +++ b/src/dc_plus/importing/powsybl/powsybl_import_helpers.py @@ -9,7 +9,7 @@ import pypowsybl import pypowsybl.loadflow -from pypowsybl.loadflow import ConnectedComponentMode, VoltageInitMode +from pypowsybl.loadflow import ComponentMode, VoltageInitMode from pypowsybl.loadflow import Parameters as LoadFlowParameters from pypowsybl.network import Network @@ -56,7 +56,7 @@ def select_a_generator_as_slack_and_run_loadflow(network: Network) -> None: read_slack_bus=True, distributed_slack=True, use_reactive_limits=True, - connected_component_mode=ConnectedComponentMode.MAIN, # ConnectedComponentMode + component_mode=ComponentMode.MAIN_CONNECTED, ) loadflow_res = pypowsybl.loadflow.run_ac(network, parameters=powsybl_loadflow_parameter)[0] diff --git a/src/dc_plus/importing/powsybl/powsybl_loadflow_parameter.py b/src/dc_plus/importing/powsybl/powsybl_loadflow_parameter.py index f447bff..e2bc69f 100644 --- a/src/dc_plus/importing/powsybl/powsybl_loadflow_parameter.py +++ b/src/dc_plus/importing/powsybl/powsybl_loadflow_parameter.py @@ -10,7 +10,7 @@ import logging from typing import Literal -from pypowsybl.loadflow import ConnectedComponentMode, VoltageInitMode +from pypowsybl.loadflow import ComponentMode, VoltageInitMode from pypowsybl.loadflow import Parameters as LoadFlowParameters logger = logging.getLogger(__name__) @@ -59,7 +59,7 @@ def get_powsybl_loadflow_parameter( # balance_type=BalanceType.PROPORTIONAL_TO_GENERATION_REMAINING_MARGIN, # BalanceType dc_use_transformer_ratio=True, countries_to_balance=None, # Sequence[str] - connected_component_mode=ConnectedComponentMode.MAIN, # ConnectedComponentMode + component_mode=ComponentMode.MAIN_CONNECTED, dc_power_factor=None, provider_parameters=provider_param, # Dict[str, str] ) diff --git a/src/dc_plus/interfaces/jacobian_network_data.py b/src/dc_plus/interfaces/jacobian_network_data.py index c784196..c64fd19 100644 --- a/src/dc_plus/interfaces/jacobian_network_data.py +++ b/src/dc_plus/interfaces/jacobian_network_data.py @@ -51,7 +51,7 @@ def _get_jacobian_data_from_network_data(dynamic_network_data: DynamicNetworkInf jacobian_index_in_use=jacobian_index_in_use, pointer_to_original_bus=np.arange(dynamic_network_data.n_buses, dtype=np.int32), jacobian=jacobian, - inverse_jacobian=sparse_inv(jacobian).toarray(), + inverse_jacobian=sparse_inv(jacobian.tocsc()).toarray(), is_angle_component=is_angle_component, is_magnitude_component=is_magnitude_component, pvpq_indices=pvpq_indices, diff --git a/src/dc_plus/jax/bsdf.py b/src/dc_plus/jax/bsdf.py index a84f597..12ba4b7 100644 --- a/src/dc_plus/jax/bsdf.py +++ b/src/dc_plus/jax/bsdf.py @@ -80,8 +80,9 @@ def _compute_bsdf_update_impl( angle_component_indices: Int[jnp.ndarray, " n_eq_jacobian"], magnitude_component_indices: Int[jnp.ndarray, " n_eq_jacobian"], ) -> Float[jnp.ndarray, " n_eq n_eq"]: - """Legacy implementation: materialize updated inverse transpose.""" - dtype = jacobian_inv_transposed.dtype + """Materialize updated inverse transpose.""" + # get dtype -> either float32 or float64 depending on input + dtype_input = v_mag_hat.dtype n_eq = jacobian_inv_transposed.shape[0] branch_from_old = jnp.take(branch_from, branches_connected_to_bus_b, axis=0) @@ -202,7 +203,7 @@ def _gather_positions(indices: jnp.ndarray) -> Tuple[jnp.ndarray, jnp.ndarray]: ) k_shape = (k_max, k_max) - delta_block = jnp.zeros(k_shape, dtype=dtype) + delta_block = jnp.zeros(k_shape, dtype=dtype_input) def _accumulate( delta: jnp.ndarray, pos: jnp.ndarray, valid: jnp.ndarray, weight: float, target: jnp.ndarray @@ -236,8 +237,8 @@ def _accumulate( theta_row = jnp.where(theta_mask, theta_new_pos, 0) mag_row = jnp.where(mag_mask, mag_new_pos, 0) - minus_one = jnp.asarray(-1.0, dtype=dtype) - zero_val = jnp.asarray(0.0, dtype=dtype) + minus_one = jnp.asarray(-1.0, dtype=dtype_input) + zero_val = jnp.asarray(0.0, dtype=dtype_input) theta_update = jnp.where(theta_mask, minus_one, zero_val) mag_update = jnp.where(mag_mask, minus_one, zero_val) delta_block = delta_block.at[theta_row, theta_row].add(theta_update) diff --git a/src/dc_plus/jax/low_rank_helper.py b/src/dc_plus/jax/low_rank_helper.py index 7089ffe..85f1af5 100644 --- a/src/dc_plus/jax/low_rank_helper.py +++ b/src/dc_plus/jax/low_rank_helper.py @@ -128,7 +128,9 @@ def _compute_branch_delta_submatrix_from_admittance( dqt_dvt = -2.0 * v_t * b_tt + v_f * (g_tf * sin_tf - b_tf * cos_tf) dqt_dvf = v_t * (g_tf * sin_tf - b_tf * cos_tf) - dtype = jnp.result_type(v_mag_from, v_mag_to, theta_from, theta_to, y_ff) + # Jacobian block entries are real-valued; avoid propagating complex dtypes + # from admittance inputs into downstream scatter-add updates. + dtype_input = jnp.result_type(v_mag_from, v_mag_to, theta_from, theta_to, g_ff) delta = jnp.array( [ @@ -137,7 +139,7 @@ def _compute_branch_delta_submatrix_from_admittance( [-dqf_dthf, -dqf_dtht, -dqf_dvf, -dqf_dvt], [-dqt_dthf, -dqt_dtht, -dqt_dvf, -dqt_dvt], ], - dtype=dtype, + dtype=dtype_input, ) return delta diff --git a/src/dc_plus/numpy/low_rank_helper.py b/src/dc_plus/numpy/low_rank_helper.py index 9299b4d..03addd7 100644 --- a/src/dc_plus/numpy/low_rank_helper.py +++ b/src/dc_plus/numpy/low_rank_helper.py @@ -209,7 +209,7 @@ def _compute_branch_delta_submatrix_from_admittance( dqt_dvt = -2.0 * v_t * b_tt + v_f * (g_tf * sin_tf - b_tf * cos_tf) dqt_dvf = v_t * (g_tf * sin_tf - b_tf * cos_tf) - dtype = np.result_type(v_mag_from, v_mag_to, theta_from, theta_to, y_ff) + dtype = np.result_type(v_mag_from, v_mag_to, theta_from, theta_to, g_ff) delta = np.array( [