From ac08671d472019c7a5132317015bb8b9028dcb60 Mon Sep 17 00:00:00 2001 From: Bpetrick <170433522+BenjPetr@users.noreply.github.com> Date: Thu, 16 Apr 2026 08:42:16 +0000 Subject: [PATCH 1/4] feat: support shunt and non PV injection reassignments for numpy code --- src/dc_plus/interfaces/network_information.py | 3 + src/dc_plus/numpy/bsdf_full_rank.py | 723 ++++++++++++++---- src/dc_plus/preprocess/create_network_data.py | 1 + tests/numpy/test_bsdf.py | 38 +- tests/test_helper/bsdf_helper.py | 332 +++++++- 5 files changed, 918 insertions(+), 179 deletions(-) diff --git a/src/dc_plus/interfaces/network_information.py b/src/dc_plus/interfaces/network_information.py index 0301d4b..6ee5bf6 100644 --- a/src/dc_plus/interfaces/network_information.py +++ b/src/dc_plus/interfaces/network_information.py @@ -618,6 +618,9 @@ class StringNetworkInformation: E.g., load, generator, etc. """ + injection_ids: StringArray + """ids of the injections, shape (n_injections,)""" + def _check_network_data_consistency( dynamic_network_data: DynamicNetworkInformation, diff --git a/src/dc_plus/numpy/bsdf_full_rank.py b/src/dc_plus/numpy/bsdf_full_rank.py index 2adf2df..46efc43 100644 --- a/src/dc_plus/numpy/bsdf_full_rank.py +++ b/src/dc_plus/numpy/bsdf_full_rank.py @@ -17,7 +17,7 @@ _compute_branch_delta_submatrix_from_admittance, ) -# ruff: noqa: ARG001, C901, PLR0913, PLR0915 +# ruff: noqa: PLR0913 def _apply_full_rank_update( @@ -53,7 +53,541 @@ def _apply_full_rank_update( return jacobian_inv - delta_inv_jacobian -# ruff: noqa: C901, PLR0915 +def _compute_shunt_delta_submatrix_from_admittance( + v_mag: Float[np.ndarray, ""], + y_shunt: Complex128[np.ndarray, ""], +) -> Float[np.ndarray, "2 2"]: + """Return the 2x2 Jacobian contribution for a shunt admittance at one bus.""" + conductance = np.real(y_shunt) + susceptance = np.imag(y_shunt) + dtype = np.result_type(v_mag, y_shunt) + + delta = np.array( + [ + [0.0, 2.0 * v_mag * conductance], + [0.0, -2.0 * v_mag * susceptance], + ], + dtype=dtype, + ) + return delta * -1 + + +def _add_bus_component_indices( + container: set[int], + bus_idx: int, + angle_idx_map: Int[np.ndarray, " n_eq_jacobian"], + magnitude_idx_map: Int[np.ndarray, " n_eq_jacobian"], +) -> None: + """Add valid Jacobian component indices for one bus into a set. + + Parameters + ---------- + container : set[int] + Mutable set collecting Jacobian indices affected by the split update. + bus_idx : int + Bus whose angle and magnitude component indices should be added. + angle_idx_map : Int[np.ndarray, " n_eq_jacobian"] + Mapping from bus indices to angle-equation indices. + magnitude_idx_map : Int[np.ndarray, " n_eq_jacobian"] + Mapping from bus indices to magnitude-equation indices. + """ + if 0 <= bus_idx < angle_idx_map.size: + theta_idx = int(angle_idx_map[bus_idx]) + if theta_idx >= 0: + container.add(theta_idx) + if 0 <= bus_idx < magnitude_idx_map.size: + mag_idx = int(magnitude_idx_map[bus_idx]) + if mag_idx >= 0: + container.add(mag_idx) + + +def _collect_targeted_indices( + bus_to_split: int, + new_bus_b_index: int, + branches_connected_to_bus_b: Int[np.ndarray, " n_branches_B"], + shunt_connected_to_bus_b: Int[np.ndarray, " n_shunts_B"], + branch_from: Int[np.ndarray, " n_branches"], + branch_to: Int[np.ndarray, " n_branches"], + shunt_to_bus: Int[np.ndarray, " n_shunts"], + angle_idx_map: Int[np.ndarray, " n_eq_jacobian"], + magnitude_idx_map: Int[np.ndarray, " n_eq_jacobian"], +) -> Int[np.ndarray, " k"]: + """Collect all Jacobian indices affected by branch and shunt reassignment. + + Parameters + ---------- + bus_to_split : int + Index of the original bus being split. + new_bus_b_index : int + Index of the new bus created by the split. + branches_connected_to_bus_b : Int[np.ndarray, " n_branches_B"] + Indices of branches reassigned to the new bus. + shunt_connected_to_bus_b : Int[np.ndarray, " n_shunts_B"] + Indices of shunts reassigned to the new bus. + branch_from : Int[np.ndarray, " n_branches"] + Branch ``from`` bus indices. + branch_to : Int[np.ndarray, " n_branches"] + Branch ``to`` bus indices. + shunt_to_bus : Int[np.ndarray, " n_shunts"] + Shunt bus indices. + angle_idx_map : Int[np.ndarray, " n_eq_jacobian"] + Mapping from bus indices to Jacobian angle-equation indices. + magnitude_idx_map : Int[np.ndarray, " n_eq_jacobian"] + Mapping from bus indices to Jacobian magnitude-equation indices. + + Returns + ------- + Int[np.ndarray, " k"] + Sorted Jacobian indices that appear in the local BSDF update block. + + Raises + ------ + IndexError + Raised when a reassigned branch or shunt index is outside the available arrays. + """ + targeted_indices: set[int] = set() + _add_bus_component_indices(targeted_indices, bus_to_split, angle_idx_map, magnitude_idx_map) + _add_bus_component_indices(targeted_indices, new_bus_b_index, angle_idx_map, magnitude_idx_map) + + for branch_idx in branches_connected_to_bus_b: + if branch_idx < 0 or branch_idx >= branch_from.size: + raise IndexError("Branch index assigned to bus B is out of bounds") + + from_bus_old = int(branch_from[branch_idx]) + to_bus_old = int(branch_to[branch_idx]) + from_bus_new = new_bus_b_index if from_bus_old == bus_to_split else from_bus_old + to_bus_new = new_bus_b_index if to_bus_old == bus_to_split else to_bus_old + + _add_bus_component_indices(targeted_indices, from_bus_old, angle_idx_map, magnitude_idx_map) + _add_bus_component_indices(targeted_indices, to_bus_old, angle_idx_map, magnitude_idx_map) + _add_bus_component_indices(targeted_indices, from_bus_new, angle_idx_map, magnitude_idx_map) + _add_bus_component_indices(targeted_indices, to_bus_new, angle_idx_map, magnitude_idx_map) + + for shunt_idx in shunt_connected_to_bus_b: + if shunt_idx < 0 or shunt_idx >= shunt_to_bus.size: + raise IndexError("Shunt index assigned to bus B is out of bounds") + + shunt_bus_old = int(shunt_to_bus[shunt_idx]) + shunt_bus_new = new_bus_b_index if shunt_bus_old == bus_to_split else shunt_bus_old + + _add_bus_component_indices(targeted_indices, shunt_bus_old, angle_idx_map, magnitude_idx_map) + _add_bus_component_indices(targeted_indices, shunt_bus_new, angle_idx_map, magnitude_idx_map) + + return np.array(sorted(targeted_indices), dtype=int) + + +def _get_branch_component_indices( + bus_from: int, + bus_to: int, + angle_idx_map: Int[np.ndarray, " n_eq_jacobian"], + magnitude_idx_map: Int[np.ndarray, " n_eq_jacobian"], +) -> Int[np.ndarray, "4"]: + """Return the Jacobian component indices touched by one branch contribution. + + Parameters + ---------- + bus_from : int + ``From`` bus of the branch contribution. + bus_to : int + ``To`` bus of the branch contribution. + angle_idx_map : Int[np.ndarray, " n_eq_jacobian"] + Mapping from bus indices to angle-equation indices. + magnitude_idx_map : Int[np.ndarray, " n_eq_jacobian"] + Mapping from bus indices to magnitude-equation indices. + + Returns + ------- + Int[np.ndarray, "4"] + Component indices ordered as angle-from, angle-to, magnitude-from, + magnitude-to. + """ + return np.array( + [ + angle_idx_map[bus_from], + angle_idx_map[bus_to], + magnitude_idx_map[bus_from], + magnitude_idx_map[bus_to], + ], + dtype=int, + ) + + +def _get_shunt_component_indices( + bus_idx: int, + angle_idx_map: Int[np.ndarray, " n_eq_jacobian"], + magnitude_idx_map: Int[np.ndarray, " n_eq_jacobian"], +) -> Int[np.ndarray, "2"]: + """Return the Jacobian component indices touched by one shunt contribution. + + Parameters + ---------- + bus_idx : int + Bus receiving the shunt contribution. + angle_idx_map : Int[np.ndarray, " n_eq_jacobian"] + Mapping from bus indices to angle-equation indices. + magnitude_idx_map : Int[np.ndarray, " n_eq_jacobian"] + Mapping from bus indices to magnitude-equation indices. + + Returns + ------- + Int[np.ndarray, "2"] + Component indices ordered as angle, magnitude. + """ + return np.array([angle_idx_map[bus_idx], magnitude_idx_map[bus_idx]], dtype=int) + + +def _accumulate_sub_delta( + delta_block: Float[np.ndarray, " k k"], + delta_matrix: Float[np.ndarray, " m m"], + component_indices: Int[np.ndarray, " m"], + position_lookup: dict[int, int], + weight: float, +) -> None: + """Add a local branch or shunt Jacobian delta into the global update block. + + Parameters + ---------- + delta_block : Float[np.ndarray, " k k"] + Global Jacobian update block being assembled for the Woodbury update. + delta_matrix : Float[np.ndarray, " m m"] + Local Jacobian contribution for one branch or one shunt. + component_indices : Int[np.ndarray, " m"] + Jacobian indices touched by the local contribution. + position_lookup : dict[int, int] + Mapping from global Jacobian index to local row/column position in + ``delta_block``. + weight : float + Scaling factor used to add or subtract the contribution. + """ + valid_positions = np.flatnonzero(component_indices >= 0) + if valid_positions.size == 0: + return + + local_indices = component_indices[valid_positions] + mapped_positions = [position_lookup[int(idx)] for idx in local_indices] + sub_delta = delta_matrix[np.ix_(valid_positions, valid_positions)] + + for row_offset, pos_row in enumerate(mapped_positions): + for col_offset, pos_col in enumerate(mapped_positions): + delta_block[pos_row, pos_col] += weight * float(sub_delta[row_offset, col_offset]) + + +def _accumulate_branch_reassignment_delta( + delta_block: Float[np.ndarray, " k k"], + position_lookup: dict[int, int], + branches_connected_to_bus_b: Int[np.ndarray, " n_branches_B"], + bus_to_split: int, + new_bus_b_index: int, + branch_from: Int[np.ndarray, " n_branches"], + branch_to: Int[np.ndarray, " n_branches"], + v_mag_hat: Float[np.ndarray, " n_buses"], + theta_hat: Float[np.ndarray, " n_buses"], + y_ff: Complex128[np.ndarray, " n_branches"], + y_ft: Complex128[np.ndarray, " n_branches"], + y_tf: Complex128[np.ndarray, " n_branches"], + y_tt: Complex128[np.ndarray, " n_branches"], + angle_idx_map: Int[np.ndarray, " n_eq_jacobian"], + magnitude_idx_map: Int[np.ndarray, " n_eq_jacobian"], +) -> None: + """Assemble the branch part of the BSDF Jacobian update block. + + Parameters + ---------- + delta_block : Float[np.ndarray, " k k"] + Global Jacobian update block being assembled. + position_lookup : dict[int, int] + Mapping from global Jacobian index to local row/column position. + branches_connected_to_bus_b : Int[np.ndarray, " n_branches_B"] + Branch indices reassigned to the new bus. + bus_to_split : int + Index of the original split bus. + new_bus_b_index : int + Index of the new bus. + branch_from : Int[np.ndarray, " n_branches"] + Branch ``from`` bus indices. + branch_to : Int[np.ndarray, " n_branches"] + Branch ``to`` bus indices. + v_mag_hat : Float[np.ndarray, " n_buses"] + Voltage magnitudes used for the local Jacobian contributions. + theta_hat : Float[np.ndarray, " n_buses"] + Voltage angles used for the local Jacobian contributions. + y_ff : Complex128[np.ndarray, " n_branches"] + Branch self-admittance at the ``from`` side. + y_ft : Complex128[np.ndarray, " n_branches"] + Branch mutual admittance from ``from`` to ``to``. + y_tf : Complex128[np.ndarray, " n_branches"] + Branch mutual admittance from ``to`` to ``from``. + y_tt : Complex128[np.ndarray, " n_branches"] + Branch self-admittance at the ``to`` side. + angle_idx_map : Int[np.ndarray, " n_eq_jacobian"] + Mapping from bus indices to angle-equation indices. + magnitude_idx_map : Int[np.ndarray, " n_eq_jacobian"] + Mapping from bus indices to magnitude-equation indices. + """ + for branch_idx in branches_connected_to_bus_b: + from_bus_old = int(branch_from[branch_idx]) + to_bus_old = int(branch_to[branch_idx]) + + delta_old = _compute_branch_delta_submatrix_from_admittance( + v_mag_from=v_mag_hat[from_bus_old], + v_mag_to=v_mag_hat[to_bus_old], + theta_from=theta_hat[from_bus_old], + theta_to=theta_hat[to_bus_old], + y_ff=y_ff[branch_idx], + y_ft=y_ft[branch_idx], + y_tf=y_tf[branch_idx], + y_tt=y_tt[branch_idx], + ) + _accumulate_sub_delta( + delta_block=delta_block, + delta_matrix=delta_old, + component_indices=_get_branch_component_indices( + from_bus_old, + to_bus_old, + angle_idx_map, + magnitude_idx_map, + ), + position_lookup=position_lookup, + weight=1.0, + ) + + from_bus_new = new_bus_b_index if from_bus_old == bus_to_split else from_bus_old + to_bus_new = new_bus_b_index if to_bus_old == bus_to_split else to_bus_old + delta_new = _compute_branch_delta_submatrix_from_admittance( + v_mag_from=v_mag_hat[from_bus_new], + v_mag_to=v_mag_hat[to_bus_new], + theta_from=theta_hat[from_bus_new], + theta_to=theta_hat[to_bus_new], + y_ff=y_ff[branch_idx], + y_ft=y_ft[branch_idx], + y_tf=y_tf[branch_idx], + y_tt=y_tt[branch_idx], + ) + _accumulate_sub_delta( + delta_block=delta_block, + delta_matrix=delta_new, + component_indices=_get_branch_component_indices( + from_bus_new, + to_bus_new, + angle_idx_map, + magnitude_idx_map, + ), + position_lookup=position_lookup, + weight=-1.0, + ) + + +def _accumulate_shunt_reassignment_delta( + delta_block: Float[np.ndarray, " k k"], + position_lookup: dict[int, int], + shunt_connected_to_bus_b: Int[np.ndarray, " n_shunts_B"], + bus_to_split: int, + new_bus_b_index: int, + shunt_to_bus: Int[np.ndarray, " n_shunts"], + v_mag_hat: Float[np.ndarray, " n_buses"], + y_shunt: Complex128[np.ndarray, " n_shunts"], + angle_idx_map: Int[np.ndarray, " n_eq_jacobian"], + magnitude_idx_map: Int[np.ndarray, " n_eq_jacobian"], +) -> None: + """Assemble the shunt part of the BSDF Jacobian update block. + + Parameters + ---------- + delta_block : Float[np.ndarray, " k k"] + Global Jacobian update block being assembled. + position_lookup : dict[int, int] + Mapping from global Jacobian index to local row/column position. + shunt_connected_to_bus_b : Int[np.ndarray, " n_shunts_B"] + Shunt indices reassigned to the new bus. + bus_to_split : int + Index of the original split bus. + new_bus_b_index : int + Index of the new bus. + shunt_to_bus : Int[np.ndarray, " n_shunts"] + Shunt bus indices. + v_mag_hat : Float[np.ndarray, " n_buses"] + Voltage magnitudes used for the local Jacobian contributions. + y_shunt : Complex128[np.ndarray, " n_shunts"] + Effective shunt admittances. + angle_idx_map : Int[np.ndarray, " n_eq_jacobian"] + Mapping from bus indices to angle-equation indices. + magnitude_idx_map : Int[np.ndarray, " n_eq_jacobian"] + Mapping from bus indices to magnitude-equation indices. + """ + for shunt_idx in shunt_connected_to_bus_b: + shunt_bus_old = int(shunt_to_bus[shunt_idx]) + shunt_bus_new = new_bus_b_index if shunt_bus_old == bus_to_split else shunt_bus_old + + delta_old = _compute_shunt_delta_submatrix_from_admittance( + v_mag=v_mag_hat[shunt_bus_old], + y_shunt=y_shunt[shunt_idx], + ) + _accumulate_sub_delta( + delta_block=delta_block, + delta_matrix=delta_old, + component_indices=_get_shunt_component_indices( + shunt_bus_old, + angle_idx_map, + magnitude_idx_map, + ), + position_lookup=position_lookup, + weight=1.0, + ) + + delta_new = _compute_shunt_delta_submatrix_from_admittance( + v_mag=v_mag_hat[shunt_bus_new], + y_shunt=y_shunt[shunt_idx], + ) + _accumulate_sub_delta( + delta_block=delta_block, + delta_matrix=delta_new, + component_indices=_get_shunt_component_indices( + shunt_bus_new, + angle_idx_map, + magnitude_idx_map, + ), + position_lookup=position_lookup, + weight=-1.0, + ) + + +def _apply_new_bus_diagonal_adjustment( + delta_block: Float[np.ndarray, " k k"], + position_lookup: dict[int, int], + new_bus_b_index: int, + angle_idx_map: Int[np.ndarray, " n_eq_jacobian"], + magnitude_idx_map: Int[np.ndarray, " n_eq_jacobian"], +) -> None: + """Apply the diagonal terms introduced by the additional PQ bus equations. + + Parameters + ---------- + delta_block : Float[np.ndarray, " k k"] + Global Jacobian update block being assembled. + position_lookup : dict[int, int] + Mapping from global Jacobian index to local row/column position. + new_bus_b_index : int + Index of the new bus. + angle_idx_map : Int[np.ndarray, " n_eq_jacobian"] + Mapping from bus indices to angle-equation indices. + magnitude_idx_map : Int[np.ndarray, " n_eq_jacobian"] + Mapping from bus indices to magnitude-equation indices. + """ + theta_idx_new = int(angle_idx_map[new_bus_b_index]) + if theta_idx_new >= 0 and theta_idx_new in position_lookup: + delta_block[position_lookup[theta_idx_new], position_lookup[theta_idx_new]] -= 1.0 + + mag_idx_new = int(magnitude_idx_map[new_bus_b_index]) + if mag_idx_new >= 0 and mag_idx_new in position_lookup: + delta_block[position_lookup[mag_idx_new], position_lookup[mag_idx_new]] -= 1.0 + + +def _build_delta_block( + idx_list: Int[np.ndarray, " k"], + dtype: np.dtype, + branches_connected_to_bus_b: Int[np.ndarray, " n_branches_B"], + shunt_connected_to_bus_b: Int[np.ndarray, " n_shunts_B"], + bus_to_split: int, + new_bus_b_index: int, + branch_from: Int[np.ndarray, " n_branches"], + branch_to: Int[np.ndarray, " n_branches"], + shunt_to_bus: Int[np.ndarray, " n_shunts"], + v_mag_hat: Float[np.ndarray, " n_buses"], + theta_hat: Float[np.ndarray, " n_buses"], + y_ff: Complex128[np.ndarray, " n_branches"], + y_ft: Complex128[np.ndarray, " n_branches"], + y_tf: Complex128[np.ndarray, " n_branches"], + y_tt: Complex128[np.ndarray, " n_branches"], + y_shunt: Complex128[np.ndarray, " n_shunts"], + angle_idx_map: Int[np.ndarray, " n_eq_jacobian"], + magnitude_idx_map: Int[np.ndarray, " n_eq_jacobian"], +) -> Float[np.ndarray, " k k"]: + """Build the local Jacobian update block used in the Woodbury correction. + + Parameters + ---------- + idx_list : Int[np.ndarray, " k"] + Sorted Jacobian indices present in the update block. + dtype : np.dtype + Data type used to allocate the update block. + branches_connected_to_bus_b : Int[np.ndarray, " n_branches_B"] + Branch indices reassigned to the new bus. + shunt_connected_to_bus_b : Int[np.ndarray, " n_shunts_B"] + Shunt indices reassigned to the new bus. + bus_to_split : int + Index of the original split bus. + new_bus_b_index : int + Index of the new bus. + branch_from : Int[np.ndarray, " n_branches"] + Branch ``from`` bus indices. + branch_to : Int[np.ndarray, " n_branches"] + Branch ``to`` bus indices. + shunt_to_bus : Int[np.ndarray, " n_shunts"] + Shunt bus indices. + v_mag_hat : Float[np.ndarray, " n_buses"] + Voltage magnitudes of the base state. + theta_hat : Float[np.ndarray, " n_buses"] + Voltage angles of the base state. + y_ff : Complex128[np.ndarray, " n_branches"] + Branch self-admittance at the ``from`` side. + y_ft : Complex128[np.ndarray, " n_branches"] + Branch mutual admittance from ``from`` to ``to``. + y_tf : Complex128[np.ndarray, " n_branches"] + Branch mutual admittance from ``to`` to ``from``. + y_tt : Complex128[np.ndarray, " n_branches"] + Branch self-admittance at the ``to`` side. + y_shunt : Complex128[np.ndarray, " n_shunts"] + Effective shunt admittances. + angle_idx_map : Int[np.ndarray, " n_eq_jacobian"] + Mapping from bus indices to angle-equation indices. + magnitude_idx_map : Int[np.ndarray, " n_eq_jacobian"] + Mapping from bus indices to magnitude-equation indices. + + Returns + ------- + Float[np.ndarray, " k k"] + Dense local Jacobian update block corresponding to ``idx_list``. + """ + position_lookup = {int(idx): pos for pos, idx in enumerate(idx_list.tolist())} + delta_block = np.zeros((idx_list.size, idx_list.size), dtype=dtype) + + _accumulate_branch_reassignment_delta( + delta_block=delta_block, + position_lookup=position_lookup, + branches_connected_to_bus_b=branches_connected_to_bus_b, + bus_to_split=bus_to_split, + new_bus_b_index=new_bus_b_index, + branch_from=branch_from, + branch_to=branch_to, + v_mag_hat=v_mag_hat, + theta_hat=theta_hat, + y_ff=y_ff, + y_ft=y_ft, + y_tf=y_tf, + y_tt=y_tt, + angle_idx_map=angle_idx_map, + magnitude_idx_map=magnitude_idx_map, + ) + _accumulate_shunt_reassignment_delta( + delta_block=delta_block, + position_lookup=position_lookup, + shunt_connected_to_bus_b=shunt_connected_to_bus_b, + bus_to_split=bus_to_split, + new_bus_b_index=new_bus_b_index, + shunt_to_bus=shunt_to_bus, + v_mag_hat=v_mag_hat, + y_shunt=y_shunt, + angle_idx_map=angle_idx_map, + magnitude_idx_map=magnitude_idx_map, + ) + _apply_new_bus_diagonal_adjustment( + delta_block=delta_block, + position_lookup=position_lookup, + new_bus_b_index=new_bus_b_index, + angle_idx_map=angle_idx_map, + magnitude_idx_map=magnitude_idx_map, + ) + return delta_block def compute_bsdf_update( @@ -78,8 +612,8 @@ def compute_bsdf_update( ) -> Float[np.ndarray, " n_eq n_eq"]: """Compute the BSDF update for a bus split with the full-rank update approach. - Note: Injection and shunt changes are currently not supported, and the function assumes a PQ split. - Only branch reassignments are handled. + Note: injection reassignments do not affect the Jacobian and therefore do not enter this update. + The function currently assumes a PQ split and supports branch and shunt reassignments. Parameters ---------- @@ -113,8 +647,8 @@ def compute_bsdf_update( The "to-from" admittance values for all branches in the original system. y_tt : Complex128[np.ndarray, " n_branches"] The "to-to" admittance values for all branches in the original system. - y_shunt : Complex128[np.ndarray, " n_buses"] - The shunt admittance values for all buses in the original system. + y_shunt : Complex128[np.ndarray, " n_shunts"] + The shunt admittance values for all shunts in the original system. angle_component_indices : Int[np.ndarray, " n_eq_jacobian"] The mapping from bus indices to angle component indices in the Jacobian. magnitude_component_indices : Int[np.ndarray, " n_eq_jacobian"] @@ -128,148 +662,49 @@ def compute_bsdf_update( if new_bus_type != 2: raise NotImplementedError("Only PQ splits are supported") - if np.asarray(shunt_connected_to_bus_b).size: - raise NotImplementedError("Shunt reassignment is not supported") - - if branches_connected_to_bus_b.size == 0: - return np.asarray(jacobian_inv, dtype=float).copy() - - base_inverse = np.asarray(jacobian_inv, dtype=float) - updated_inverse = base_inverse.copy() - - branch_from_arr = np.asarray(branch_from, dtype=int).reshape(-1) - branch_to_arr = np.asarray(branch_to, dtype=int).reshape(-1) - - angle_idx_map = np.asarray(angle_component_indices, dtype=int).reshape(-1) - magnitude_idx_map = np.asarray(magnitude_component_indices, dtype=int).reshape(-1) - - if new_bus_b_index >= angle_idx_map.size or new_bus_b_index >= magnitude_idx_map.size: + if branches_connected_to_bus_b.size == 0 and shunt_connected_to_bus_b.size == 0: + return jacobian_inv.copy() + if new_bus_b_index >= angle_component_indices.size or new_bus_b_index >= magnitude_component_indices.size: raise IndexError("New bus index is out of bounds for component index arrays") + idx_list = _collect_targeted_indices( + bus_to_split=bus_to_split, + new_bus_b_index=new_bus_b_index, + branches_connected_to_bus_b=branches_connected_to_bus_b, + shunt_connected_to_bus_b=shunt_connected_to_bus_b, + branch_from=branch_from, + branch_to=branch_to, + shunt_to_bus=shunt_to_bus, + angle_idx_map=angle_component_indices, + magnitude_idx_map=magnitude_component_indices, + ) + if idx_list.size == 0: + return jacobian_inv.copy() - v_mag_vec = np.asarray(v_mag_hat, dtype=float).reshape(-1) - theta_vec = np.asarray(theta_hat, dtype=float).reshape(-1) - - y_ff_vec = np.asarray(y_ff, dtype=np.complex128).reshape(-1) - y_ft_vec = np.asarray(y_ft, dtype=np.complex128).reshape(-1) - y_tf_vec = np.asarray(y_tf, dtype=np.complex128).reshape(-1) - y_tt_vec = np.asarray(y_tt, dtype=np.complex128).reshape(-1) - - def _add_component_indices(container: set[int], bus_idx: int) -> None: - if 0 <= bus_idx < angle_idx_map.size: - theta_idx = int(angle_idx_map[bus_idx]) - if theta_idx >= 0: - container.add(theta_idx) - if 0 <= bus_idx < magnitude_idx_map.size: - mag_idx = int(magnitude_idx_map[bus_idx]) - if mag_idx >= 0: - container.add(mag_idx) - - targeted_indices: set[int] = set() - - _add_component_indices(targeted_indices, bus_to_split) - _add_component_indices(targeted_indices, new_bus_b_index) - - # remove components from bus_to_split - for branch_idx in branches_connected_to_bus_b: - if branch_idx < 0 or branch_idx >= branch_from_arr.size: - raise IndexError("Branch index assigned to bus B is out of bounds") - - from_bus_old = int(branch_from_arr[branch_idx]) - to_bus_old = int(branch_to_arr[branch_idx]) - - _add_component_indices(targeted_indices, from_bus_old) - _add_component_indices(targeted_indices, to_bus_old) - - from_bus_new = new_bus_b_index if from_bus_old == bus_to_split else from_bus_old - to_bus_new = new_bus_b_index if to_bus_old == bus_to_split else to_bus_old - - _add_component_indices(targeted_indices, from_bus_new) - _add_component_indices(targeted_indices, to_bus_new) - - if not targeted_indices: - return updated_inverse - - idx_list = np.array(sorted(targeted_indices), dtype=int) - position_lookup = {idx: pos for pos, idx in enumerate(idx_list.tolist())} - - delta_block = np.zeros((idx_list.size, idx_list.size), dtype=updated_inverse.dtype) - - def _accumulate_branch_delta( - delta_matrix: Float[np.ndarray, "4 4"], - bus_from: int, - bus_to: int, - weight: float, - ) -> None: - component_indices = np.array( - [ - angle_idx_map[bus_from], - angle_idx_map[bus_to], - magnitude_idx_map[bus_from], - magnitude_idx_map[bus_to], - ], - dtype=int, - ) - valid_positions = np.flatnonzero(component_indices >= 0) - if valid_positions.size == 0: - return - - local_indices = component_indices[valid_positions] - mapped_positions = [position_lookup[idx] for idx in local_indices] - sub_delta = delta_matrix[np.ix_(valid_positions, valid_positions)] - - for row_offset, pos_row in enumerate(mapped_positions): - for col_offset, pos_col in enumerate(mapped_positions): - delta_block[pos_row, pos_col] += weight * float(sub_delta[row_offset, col_offset]) - - # Accumulate the contributions from the branches that are reassigned to bus B - for branch_idx in branches_connected_to_bus_b: - from_bus_old = int(branch_from_arr[branch_idx]) - to_bus_old = int(branch_to_arr[branch_idx]) - - delta_old = _compute_branch_delta_submatrix_from_admittance( - v_mag_from=v_mag_vec[from_bus_old], - v_mag_to=v_mag_vec[to_bus_old], - theta_from=theta_vec[from_bus_old], - theta_to=theta_vec[to_bus_old], - y_ff=y_ff_vec[branch_idx], - y_ft=y_ft_vec[branch_idx], - y_tf=y_tf_vec[branch_idx], - y_tt=y_tt_vec[branch_idx], - ) - _accumulate_branch_delta(delta_old, from_bus_old, to_bus_old, weight=1.0) - - from_bus_new = new_bus_b_index if from_bus_old == bus_to_split else from_bus_old - to_bus_new = new_bus_b_index if to_bus_old == bus_to_split else to_bus_old - - delta_new = _compute_branch_delta_submatrix_from_admittance( - v_mag_from=v_mag_vec[from_bus_new], - v_mag_to=v_mag_vec[to_bus_new], - theta_from=theta_vec[from_bus_new], - theta_to=theta_vec[to_bus_new], - y_ff=y_ff_vec[branch_idx], - y_ft=y_ft_vec[branch_idx], - y_tf=y_tf_vec[branch_idx], - y_tt=y_tt_vec[branch_idx], - ) - _accumulate_branch_delta(delta_new, from_bus_new, to_bus_new, weight=-1.0) - - theta_idx_new = int(angle_idx_map[new_bus_b_index]) - if theta_idx_new >= 0 and theta_idx_new in position_lookup: - pos = position_lookup[theta_idx_new] - delta_block[pos, pos] -= 1.0 - - mag_idx_new = int(magnitude_idx_map[new_bus_b_index]) - if mag_idx_new >= 0 and mag_idx_new in position_lookup: - pos = position_lookup[mag_idx_new] - delta_block[pos, pos] -= 1.0 - + delta_block = _build_delta_block( + idx_list=idx_list, + dtype=jacobian_inv.dtype, + branches_connected_to_bus_b=branches_connected_to_bus_b, + shunt_connected_to_bus_b=shunt_connected_to_bus_b, + bus_to_split=bus_to_split, + new_bus_b_index=new_bus_b_index, + branch_from=branch_from, + branch_to=branch_to, + shunt_to_bus=shunt_to_bus, + v_mag_hat=v_mag_hat, + theta_hat=theta_hat, + y_ff=y_ff, + y_ft=y_ft, + y_tf=y_tf, + y_tt=y_tt, + y_shunt=y_shunt, + angle_idx_map=angle_component_indices, + magnitude_idx_map=magnitude_component_indices, + ) if not np.any(delta_block): - return updated_inverse + return jacobian_inv.copy() - updated_inverse = _apply_full_rank_update( - jacobian_inv=base_inverse, + return _apply_full_rank_update( + jacobian_inv=jacobian_inv, jacobian_delta_submatrix=delta_block, idx_list=idx_list, ) - - return updated_inverse diff --git a/src/dc_plus/preprocess/create_network_data.py b/src/dc_plus/preprocess/create_network_data.py index d7f4470..c89828f 100644 --- a/src/dc_plus/preprocess/create_network_data.py +++ b/src/dc_plus/preprocess/create_network_data.py @@ -158,6 +158,7 @@ def _create_network_data( branch_ids=branches["id_str"].values, limit_names=limits["name"].values, injection_types=injections["injection_type"].values, + injection_ids=injections["id_str"].values, ) _check_network_data_consistency(dynamic_network_data=dynamic_info, string_network_data=string_info) diff --git a/tests/numpy/test_bsdf.py b/tests/numpy/test_bsdf.py index 737669e..93426cb 100644 --- a/tests/numpy/test_bsdf.py +++ b/tests/numpy/test_bsdf.py @@ -8,8 +8,10 @@ import numpy as np import pytest +from dc_plus.importing.powsybl.powsybl_import import DANGLING_BUS_STRING_SUFFIX from dc_plus.numpy.bsdf_full_rank import compute_bsdf_update from tests.test_helper.bsdf_helper import ( + derive_bus_order, get_bsdf_cases, prepare_bsdf_test_context, run_reference_one_step, @@ -27,8 +29,8 @@ def test_bsdf_full_rank(bsdf_test_case): bus_to_split=setup.bus_to_split, new_bus_b_index=setup.new_bus_index, new_bus_type=2, # force select PQ node - branches_connected_to_bus_b=setup.branches_to_move, - shunt_connected_to_bus_b=np.array([], dtype=np.int32), + branches_connected_to_bus_b=setup.branches_connected_to_bus_b, + shunt_connected_to_bus_b=setup.shunt_connected_to_bus_b, branch_from=setup.branch_from_original, branch_to=setup.branch_to_original, shunt_to_bus=setup.dynamic_info.shunt_bus_indices, @@ -54,8 +56,9 @@ def test_bsdf_full_rank(bsdf_test_case): ) # test against powsybl - dynamic_info_one_step = run_reference_one_step(setup.net, bsdf_test_case=bsdf_test_case) - dx = -jacobian_inv_bsdf @ setup.mismatch_n1 + dynamic_info_one_step, string_info_one_step = run_reference_one_step(setup.net, bsdf_test_case=bsdf_test_case) + bus_order = derive_bus_order(setup.split_bus_ids, string_info_one_step.bus_ids) + dx = -jacobian_inv_bsdf @ setup.mismatch_bsdf_reference # Map Jacobian increments back to bus ordering using the Jacobian mapping @@ -71,16 +74,33 @@ def test_bsdf_full_rank(bsdf_test_case): vm_actual[setup.pq_indices] + dx[setup.jacobian_data_with_extra_buses.is_magnitude_component] ) - # reorder bus voltages to match new ordering of manual split + dangling_bus_mask = np.char.endswith(np.asarray(string_info_one_step.bus_ids, dtype=str), DANGLING_BUS_STRING_SUFFIX) + np.testing.assert_allclose( - dynamic_info_one_step.bus_voltage_magnitudes[bsdf_test_case.bus_order], - vm_updated_J, + dynamic_info_one_step.bus_voltage_magnitudes[~dangling_bus_mask], + vm_updated_J[bus_order][~dangling_bus_mask], rtol=1e-10, atol=1e-10, ) np.testing.assert_allclose( - dynamic_info_one_step.bus_voltage_angles_rad[bsdf_test_case.bus_order], - theta_updated_J, + dynamic_info_one_step.bus_voltage_angles_rad[~dangling_bus_mask], + theta_updated_J[bus_order][~dangling_bus_mask], rtol=1e-10, atol=1e-10, ) + + # FIXME + # Powsybl one-step results on dangling buses are not bitwise identical to the + # imported Jacobian linearization + np.testing.assert_allclose( + dynamic_info_one_step.bus_voltage_magnitudes[dangling_bus_mask], + vm_updated_J[bus_order][dangling_bus_mask], + rtol=1e-9, + atol=1e-9, + ) + np.testing.assert_allclose( + dynamic_info_one_step.bus_voltage_angles_rad[dangling_bus_mask], + theta_updated_J[bus_order][dangling_bus_mask], + rtol=1e-8, + atol=1e-8, + ) diff --git a/tests/test_helper/bsdf_helper.py b/tests/test_helper/bsdf_helper.py index 859f4b4..0fe6e0d 100644 --- a/tests/test_helper/bsdf_helper.py +++ b/tests/test_helper/bsdf_helper.py @@ -11,7 +11,10 @@ import numpy as np import pypowsybl -from dc_plus.example_grids.pypowsbl.example_grids import basic_node_breaker_network_powsybl +from dc_plus.example_grids.pypowsbl.example_grids import ( + basic_node_breaker_network_powsybl, + create_complex_grid_battery_hvdc_svc_3w_trafo, +) from dc_plus.importing.powsybl.powsybl_network_helpers import _load_test_grid from dc_plus.importing.powsybl.powsybl_loadflow_parameter import get_powsybl_loadflow_parameter from dc_plus.interfaces.jacobian_network_data import ( @@ -19,19 +22,36 @@ _get_jacobian_data_from_network_data, calculate_nodal_mismatch_network_data, ) -from dc_plus.interfaces.network_information import DynamicNetworkInformation +from dc_plus.interfaces.network_information import DynamicNetworkInformation, StringNetworkInformation from dc_plus.preprocess.create_network_data import create_network_data_pypowsbl from dc_plus.preprocess.preprocess_jacobian_bsdf import preprocess_jacobian_bsdf +_NEW_BUS_PLACEHOLDER_ID = "__new_bus_b__" + + @dataclass class BsdfTestCase: get_net: Callable[..., Any] + """Callable that returns a pypowsybl network instance for the test case.""" bus_to_split: int - branches_to_move: np.ndarray + """Index of the bus to split in the BSDF test case. + Note: this index number refers to the DC+ internal ordering of buses, which differs + from the pypowsybl bus numbering. (string index is sorted by name -> internal index) + """ + branches_connected_to_bus_b_string: list[str] + """Powsybl/DC+ branch IDs connected to the new bus B after the split.""" + + shunt_connected_to_bus_b_string: list[str] + """Powsybl/DC+ shunt IDs connected to the new bus B after the split.""" + injections_connected_to_bus_b_string: list[str] + """Powsybl/DC+ injection IDs connected to the new bus B after the split. + Note: you may only reassign injections which have no PV regulation. + If you want to reassign a PV injection, you need to change the bus b from default type PQ to PV.""" open_switches: tuple[str, ...] + """pypowsybl IDs of switches to open in the reference switch pattern for the test case.""" close_switches: tuple[str, ...] - bus_order: np.ndarray + """pypowsybl IDs of switches to close in the reference switch pattern for the test case.""" @dataclass @@ -43,8 +63,11 @@ class BsdfTestContext: jacobian_data_with_extra_buses: Any jacobian_data_split_manual: Any new_bus_index: int + split_bus_ids: np.ndarray bus_to_split: int - branches_to_move: np.ndarray + branches_connected_to_bus_b: np.ndarray + shunt_connected_to_bus_b: np.ndarray + injections_connected_to_bus_b: np.ndarray y_ff: np.ndarray y_ft: np.ndarray y_tf: np.ndarray @@ -53,7 +76,7 @@ class BsdfTestContext: branch_to_original: np.ndarray v_mag_hat: np.ndarray theta_hat: np.ndarray - mismatch_n1: np.ndarray + mismatch_bsdf_reference: np.ndarray theta_base: np.ndarray vm_base: np.ndarray pvpq_indices: np.ndarray @@ -62,9 +85,13 @@ class BsdfTestContext: def get_bsdf_cases() -> List[BsdfTestCase]: """Return test cases for BSDF tests.""" + + # simple_bsdf_test_case + # based on basic_node_breaker_network_powsybl + # simple reasignment and get_net = basic_node_breaker_network_powsybl bus_to_split = 2 - branches_to_move = np.array([2, 5], dtype=np.int32) + branches_connected_to_bus_b_string = ["L3", "L6"] open_switches = ( "VL3_BREAKER", "L32_DISCONNECTOR_3_0", @@ -78,29 +105,246 @@ def get_bsdf_cases() -> List[BsdfTestCase]: "L72_DISCONNECTOR_7_0", "L62_DISCONNECTOR_5_1", ) - bus_order = np.asarray([0, 1, 2, 4, 5, 3], dtype=np.int32) - return [ - BsdfTestCase( - get_net=get_net, - bus_to_split=bus_to_split, - branches_to_move=branches_to_move, - open_switches=open_switches, - close_switches=close_switches, - bus_order=bus_order, + simple_bsdf_test_case = BsdfTestCase( + get_net=get_net, + bus_to_split=bus_to_split, + branches_connected_to_bus_b_string=branches_connected_to_bus_b_string, + open_switches=open_switches, + close_switches=close_switches, + shunt_connected_to_bus_b_string=[], + injections_connected_to_bus_b_string=[], + ) + + get_net = create_complex_grid_battery_hvdc_svc_3w_trafo + bus_to_split = 13 + branches_connected_to_bus_b_string = ["L12", "L6", "L7"] + open_switches = ( + "VL_MV_BREAKER", + "SHUNT_MV_DISCONNECTOR_21_0", + ) + close_switches = ("SHUNT_MV_DISCONNECTOR_21_1",) + shunt_injection_test_case = BsdfTestCase( + get_net=get_net, + bus_to_split=bus_to_split, + branches_connected_to_bus_b_string=branches_connected_to_bus_b_string, + open_switches=open_switches, + close_switches=close_switches, + shunt_connected_to_bus_b_string=["SHUNT_MV"], + injections_connected_to_bus_b_string=["BAT_MV"], + ) + + return [simple_bsdf_test_case, shunt_injection_test_case] + + +def _reassign_selected_bus_indices( + asset_bus_indices: np.ndarray, + selected_asset_indices: np.ndarray, + bus_to_split: int, + new_bus_index: int, +) -> np.ndarray: + """Reassign the selected assets from the split bus to the new bus.""" + reassigned_bus_indices = np.asarray(asset_bus_indices).copy() + for asset_idx in np.asarray(selected_asset_indices, dtype=np.int32).reshape(-1): + reassigned_bus_indices[asset_idx] = np.where( + reassigned_bus_indices[asset_idx] == bus_to_split, + new_bus_index, + reassigned_bus_indices[asset_idx], + ) + return reassigned_bus_indices + + +def _build_split_bus_ids(original_bus_ids: np.ndarray) -> np.ndarray: + """Append a placeholder bus ID for the newly created split bus.""" + bus_ids = np.asarray(original_bus_ids, dtype=object) + return np.concatenate((bus_ids, np.asarray([_NEW_BUS_PLACEHOLDER_ID], dtype=object))) + + +def derive_bus_order(candidate_bus_ids: np.ndarray, reference_bus_ids: np.ndarray) -> np.ndarray: + """Map candidate bus ordering to the powsybl reference bus ordering. + + The candidate ordering contains all original bus IDs plus one placeholder for the + newly created split bus. The reference ordering contains the actual new bus ID. + """ + candidate_ids = np.asarray(candidate_bus_ids, dtype=str) + reference_ids = np.asarray(reference_bus_ids, dtype=str) + + if candidate_ids.size != reference_ids.size: + raise ValueError( + "Candidate and reference bus ID arrays must have the same length to derive bus ordering. " + f"Got candidate={candidate_ids.size}, reference={reference_ids.size}." + ) + + candidate_lookup = { + bus_id: idx for idx, bus_id in enumerate(candidate_ids.tolist()) if bus_id != _NEW_BUS_PLACEHOLDER_ID + } + bus_order = np.full(reference_ids.size, -1, dtype=np.int32) + used_candidate_indices: set[int] = set() + unresolved_reference_positions: list[int] = [] + + for ref_pos, ref_bus_id in enumerate(reference_ids.tolist()): + candidate_idx = candidate_lookup.get(ref_bus_id) + if candidate_idx is None: + unresolved_reference_positions.append(ref_pos) + continue + + bus_order[ref_pos] = candidate_idx + used_candidate_indices.add(candidate_idx) + + remaining_candidate_indices = [ + idx + for idx in range(candidate_ids.size) + if idx not in used_candidate_indices and candidate_ids[idx] == _NEW_BUS_PLACEHOLDER_ID + ] + + if len(unresolved_reference_positions) != len(remaining_candidate_indices): + raise ValueError( + "Could not derive a unique bus order from bus IDs. " + f"Unresolved reference positions={unresolved_reference_positions}, " + f"remaining candidate indices={remaining_candidate_indices}, " + f"reference bus ids={reference_ids.tolist()}, candidate bus ids={candidate_ids.tolist()}." + ) + + for ref_pos, candidate_idx in zip(unresolved_reference_positions, remaining_candidate_indices): + bus_order[ref_pos] = candidate_idx + + return bus_order + + +def _aggregate_bus_injections( + injection_to_bus: np.ndarray, + injection_active_power: np.ndarray, + injection_reactive_power: np.ndarray, + injection_connected: np.ndarray, + n_buses: int, +) -> tuple[np.ndarray, np.ndarray]: + """Aggregate per-injection powers back to bus totals.""" + injection_to_bus_arr = np.asarray(injection_to_bus, dtype=int) + injection_active_power_arr = np.asarray(injection_active_power, dtype=float) + injection_reactive_power_arr = np.asarray(injection_reactive_power, dtype=float) + injection_connected_arr = np.asarray(injection_connected, dtype=bool) + + if injection_to_bus_arr.ndim == 1: + bus_active_power = np.zeros(n_buses, dtype=injection_active_power_arr.dtype) + bus_reactive_power = np.zeros(n_buses, dtype=injection_reactive_power_arr.dtype) + connected_mask = injection_connected_arr.astype(bool) + np.add.at(bus_active_power, injection_to_bus_arr[connected_mask], injection_active_power_arr[connected_mask]) + np.add.at( + bus_reactive_power, + injection_to_bus_arr[connected_mask], + injection_reactive_power_arr[connected_mask], + ) + return bus_active_power, bus_reactive_power + + n_timesteps = injection_to_bus_arr.shape[1] + bus_active_power = np.zeros((n_buses, n_timesteps), dtype=injection_active_power_arr.dtype) + bus_reactive_power = np.zeros((n_buses, n_timesteps), dtype=injection_reactive_power_arr.dtype) + + for timestep in range(n_timesteps): + connected_mask = injection_connected_arr[:, timestep].astype(bool) + np.add.at( + bus_active_power[:, timestep], + injection_to_bus_arr[connected_mask, timestep], + injection_active_power_arr[connected_mask, timestep], + ) + np.add.at( + bus_reactive_power[:, timestep], + injection_to_bus_arr[connected_mask, timestep], + injection_reactive_power_arr[connected_mask, timestep], ) + + return bus_active_power, bus_reactive_power + + +def _resolve_asset_indices_by_id( + available_ids: np.ndarray, + selected_ids: list[str], + asset_kind: str, +) -> np.ndarray: + """Resolve configured asset IDs to their DC+ internal indices.""" + available_ids_arr = np.asarray(available_ids, dtype=str) + selected_ids_list = list(selected_ids) + + id_to_index = {asset_id: idx for idx, asset_id in enumerate(available_ids_arr.tolist())} + missing_ids = [asset_id for asset_id in selected_ids_list if asset_id not in id_to_index] + if missing_ids: + raise ValueError( + f"Unknown {asset_kind} IDs in BSDF test case: {missing_ids}. Available IDs are: {available_ids_arr.tolist()}" + ) + + return np.asarray([id_to_index[asset_id] for asset_id in selected_ids_list], dtype=np.int32) + + +def _validate_reassignment_targets( + bus_to_split: int, + branches_connected_to_bus_b: np.ndarray, + shunt_connected_to_bus_b: np.ndarray, + injections_connected_to_bus_b: np.ndarray, + branch_from_bus: np.ndarray, + branch_to_bus: np.ndarray, + shunt_bus_indices: np.ndarray, + injection_to_bus: np.ndarray, +) -> None: + """Validate that all selected assets are currently attached to the split bus.""" + invalid_branch_indices = [ + int(branch_idx) + for branch_idx in np.asarray(branches_connected_to_bus_b, dtype=np.int32).reshape(-1) + if not np.any(np.asarray(branch_from_bus[branch_idx]) == bus_to_split) + and not np.any(np.asarray(branch_to_bus[branch_idx]) == bus_to_split) ] + invalid_shunt_indices = [ + int(shunt_idx) + for shunt_idx in np.asarray(shunt_connected_to_bus_b, dtype=np.int32).reshape(-1) + if not np.any(np.asarray(shunt_bus_indices[shunt_idx]) == bus_to_split) + ] + invalid_injection_indices = [ + int(injection_idx) + for injection_idx in np.asarray(injections_connected_to_bus_b, dtype=np.int32).reshape(-1) + if not np.any(np.asarray(injection_to_bus[injection_idx]) == bus_to_split) + ] + + if invalid_branch_indices or invalid_shunt_indices or invalid_injection_indices: + raise ValueError( + "Invalid BSDF test case selection for bus split " + f"{bus_to_split}: branches={invalid_branch_indices}, " + f"shunts={invalid_shunt_indices}, injections={invalid_injection_indices}" + ) def prepare_bsdf_test_context(bsdf_test_case: BsdfTestCase) -> BsdfTestContext: """Build split-bus context shared by BSDF tests to avoid duplication.""" - branches_to_move = np.asarray(bsdf_test_case.branches_to_move, dtype=np.int32) - net, _, dynamic_info, _, jacobian_data = _load_test_grid(bsdf_test_case.get_net) + net, _static_info, dynamic_info, string_info, jacobian_data = _load_test_grid(bsdf_test_case.get_net) + branches_connected_to_bus_b = _resolve_asset_indices_by_id( + available_ids=string_info.branch_ids, + selected_ids=bsdf_test_case.branches_connected_to_bus_b_string, + asset_kind="branch", + ) + shunt_connected_to_bus_b = _resolve_asset_indices_by_id( + available_ids=string_info.shunt_ids, + selected_ids=bsdf_test_case.shunt_connected_to_bus_b_string, + asset_kind="shunt", + ) + injections_connected_to_bus_b = _resolve_asset_indices_by_id( + available_ids=string_info.injection_ids, + selected_ids=bsdf_test_case.injections_connected_to_bus_b_string, + asset_kind="injection", + ) jacobian_data_with_extra_buses, dynamic_info_with_placeholders = preprocess_jacobian_bsdf( jacobian_data=jacobian_data, max_bus_splits=1, dynamic_network_data=dynamic_info, ) + _validate_reassignment_targets( + bus_to_split=bsdf_test_case.bus_to_split, + branches_connected_to_bus_b=branches_connected_to_bus_b, + shunt_connected_to_bus_b=shunt_connected_to_bus_b, + injections_connected_to_bus_b=injections_connected_to_bus_b, + branch_from_bus=dynamic_info.branch_from_bus, + branch_to_bus=dynamic_info.branch_to_bus, + shunt_bus_indices=dynamic_info.shunt_bus_indices, + injection_to_bus=dynamic_info.injection_to_bus, + ) new_bus_index = dynamic_info_with_placeholders.n_buses - 1 + split_bus_ids = _build_split_bus_ids(string_info.bus_ids) v_mag_placeholder = dynamic_info_with_placeholders.bus_voltage_magnitudes.copy() theta_placeholder = dynamic_info_with_placeholders.bus_voltage_angles_rad.copy() v_mag_placeholder[new_bus_index] = dynamic_info.bus_voltage_magnitudes[bsdf_test_case.bus_to_split] @@ -112,7 +356,7 @@ def prepare_bsdf_test_context(bsdf_test_case: BsdfTestCase) -> BsdfTestContext: ) branch_from_split = dynamic_info.branch_from_bus.copy() branch_to_split = dynamic_info.branch_to_bus.copy() - for branch_idx in branches_to_move: + for branch_idx in branches_connected_to_bus_b: branch_from_split[branch_idx] = np.where( branch_from_split[branch_idx] == bsdf_test_case.bus_to_split, new_bus_index, @@ -123,10 +367,35 @@ def prepare_bsdf_test_context(bsdf_test_case: BsdfTestCase) -> BsdfTestContext: new_bus_index, branch_to_split[branch_idx], ) + + shunt_bus_split = _reassign_selected_bus_indices( + asset_bus_indices=dynamic_info_with_placeholders.shunt_bus_indices, + selected_asset_indices=shunt_connected_to_bus_b, + bus_to_split=bsdf_test_case.bus_to_split, + new_bus_index=new_bus_index, + ) + injection_to_bus_split = _reassign_selected_bus_indices( + asset_bus_indices=dynamic_info_with_placeholders.injection_to_bus, + selected_asset_indices=injections_connected_to_bus_b, + bus_to_split=bsdf_test_case.bus_to_split, + new_bus_index=new_bus_index, + ) + bus_active_power_split, bus_reactive_power_split = _aggregate_bus_injections( + injection_to_bus=injection_to_bus_split, + injection_active_power=dynamic_info_with_placeholders.injection_active_power, + injection_reactive_power=dynamic_info_with_placeholders.injection_reactive_power, + injection_connected=dynamic_info_with_placeholders.injection_connected, + n_buses=dynamic_info_with_placeholders.n_buses, + ) + dynamic_info_split_manual = replace( dynamic_info_with_placeholders, branch_from_bus=branch_from_split, branch_to_bus=branch_to_split, + shunt_bus_indices=shunt_bus_split, + injection_to_bus=injection_to_bus_split, + bus_active_power=bus_active_power_split, + bus_reactive_power=bus_reactive_power_split, ) y_ff = np.asarray(dynamic_info.branch_effective_admittance_from_from, dtype=np.complex128) y_ft = np.asarray(dynamic_info.branch_effective_admittance_from_to, dtype=np.complex128) @@ -137,7 +406,9 @@ def prepare_bsdf_test_context(bsdf_test_case: BsdfTestCase) -> BsdfTestContext: v_mag_hat = np.asarray(dynamic_info_with_placeholders.bus_voltage_magnitudes, dtype=float).flatten() theta_hat = np.asarray(dynamic_info_with_placeholders.bus_voltage_angles_rad, dtype=float).flatten() y_matrix_n1 = _get_admittance_matrix_from_network_data(dynamic_info_split_manual) - mismatch_n1 = calculate_nodal_mismatch_network_data(dynamic_network_data=dynamic_info_split_manual, y_matrix=y_matrix_n1) + mismatch_bsdf_reference = calculate_nodal_mismatch_network_data( + dynamic_network_data=dynamic_info_split_manual, y_matrix=y_matrix_n1 + ) theta_base = np.asarray(dynamic_info_split_manual.bus_voltage_angles_rad, dtype=float).flatten() vm_base = np.asarray(dynamic_info_split_manual.bus_voltage_magnitudes, dtype=float).flatten() jacobian_data_split_manual = _get_jacobian_data_from_network_data(dynamic_info_split_manual) @@ -151,8 +422,11 @@ def prepare_bsdf_test_context(bsdf_test_case: BsdfTestCase) -> BsdfTestContext: jacobian_data_with_extra_buses=jacobian_data_with_extra_buses, jacobian_data_split_manual=jacobian_data_split_manual, new_bus_index=new_bus_index, + split_bus_ids=split_bus_ids, bus_to_split=bsdf_test_case.bus_to_split, - branches_to_move=branches_to_move, + branches_connected_to_bus_b=branches_connected_to_bus_b, + shunt_connected_to_bus_b=shunt_connected_to_bus_b, + injections_connected_to_bus_b=injections_connected_to_bus_b, y_ff=y_ff, y_ft=y_ft, y_tf=y_tf, @@ -161,7 +435,7 @@ def prepare_bsdf_test_context(bsdf_test_case: BsdfTestCase) -> BsdfTestContext: branch_to_original=branch_to_original, v_mag_hat=v_mag_hat, theta_hat=theta_hat, - mismatch_n1=mismatch_n1, + mismatch_bsdf_reference=mismatch_bsdf_reference, theta_base=theta_base, vm_base=vm_base, pvpq_indices=pvpq_indices, @@ -169,7 +443,10 @@ def prepare_bsdf_test_context(bsdf_test_case: BsdfTestCase) -> BsdfTestContext: ) -def run_reference_one_step(net: Any, bsdf_test_case: BsdfTestCase) -> DynamicNetworkInformation: +def run_reference_one_step( + net: Any, + bsdf_test_case: BsdfTestCase, +) -> tuple[DynamicNetworkInformation, StringNetworkInformation]: """Apply the reference switch pattern and run one-step AC load-flow for the given case.""" for switch_id in bsdf_test_case.open_switches: net.open_switch(switch_id) @@ -177,6 +454,9 @@ def run_reference_one_step(net: Any, bsdf_test_case: BsdfTestCase) -> DynamicNet net.close_switch(switch_id) loadflow_parameter = get_powsybl_loadflow_parameter("one_step") - pypowsybl.loadflow.run_ac(net, parameters=loadflow_parameter)[0] - _static_info, dynamic_info_one_step, _string_info = create_network_data_pypowsbl(net) - return dynamic_info_one_step + lf_res = pypowsybl.loadflow.run_ac(net, parameters=loadflow_parameter)[0] + assert lf_res.iteration_count == 1, ( + f"Expected the reference load-flow to converge in one step, but got {lf_res.iteration_count} iterations." + ) + _static_info, dynamic_info_one_step, string_info_one_step = create_network_data_pypowsbl(net) + return dynamic_info_one_step, string_info_one_step From 6595c1ea7c86b6c82270cc161625118b53017a24 Mon Sep 17 00:00:00 2001 From: Bpetrick <170433522+BenjPetr@users.noreply.github.com> Date: Thu, 16 Apr 2026 08:45:18 +0000 Subject: [PATCH 2/4] feat: support shunt and non PV injection reassignments for jax code --- src/dc_plus/jax/bsdf.py | 70 ++++++++++++++++++++++++++++++++++++-- tests/jax/test_bsdf_jax.py | 35 ++++++++++++++----- 2 files changed, 93 insertions(+), 12 deletions(-) diff --git a/src/dc_plus/jax/bsdf.py b/src/dc_plus/jax/bsdf.py index a84f597..e13f66a 100644 --- a/src/dc_plus/jax/bsdf.py +++ b/src/dc_plus/jax/bsdf.py @@ -63,20 +63,42 @@ def _full_rank_lodf( return lodf_transposed +def _compute_shunt_delta_submatrix_from_admittance( + v_mag: Float[jnp.ndarray, ""], + y_shunt: Complex128[jnp.ndarray, ""], +) -> Float[jnp.ndarray, "2 2"]: + """Return the 2x2 Jacobian contribution for a shunt admittance at one bus.""" + conductance = jnp.real(y_shunt) + susceptance = jnp.imag(y_shunt) + dtype = jnp.result_type(v_mag, y_shunt) + + delta = jnp.array( + [ + [0.0, 2.0 * v_mag * conductance], + [0.0, -2.0 * v_mag * susceptance], + ], + dtype=dtype, + ) + return delta * -1 + + @jax.jit def _compute_bsdf_update_impl( jacobian_inv_transposed: Float[jnp.ndarray, " n_eq n_eq"], bus_to_split: int, new_bus_b_index: int, branches_connected_to_bus_b: Int[jnp.ndarray, " n_branches_B"], + shunt_connected_to_bus_b: Int[jnp.ndarray, " n_shunts_B"], branch_from: Int[jnp.ndarray, " n_branches"], branch_to: Int[jnp.ndarray, " n_branches"], + shunt_to_bus: Int[jnp.ndarray, " n_shunts"], v_mag_hat: Float[jnp.ndarray, " n_buses"], theta_hat: Float[jnp.ndarray, " n_buses"], 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"], + y_shunt: Complex128[jnp.ndarray, " n_shunts"], 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"]: @@ -86,9 +108,11 @@ def _compute_bsdf_update_impl( branch_from_old = jnp.take(branch_from, branches_connected_to_bus_b, axis=0) branch_to_old = jnp.take(branch_to, branches_connected_to_bus_b, axis=0) + shunt_bus_old = jnp.take(shunt_to_bus, shunt_connected_to_bus_b, axis=0) branch_from_new = jnp.where(branch_from_old == bus_to_split, new_bus_b_index, branch_from_old) branch_to_new = jnp.where(branch_to_old == bus_to_split, new_bus_b_index, branch_to_old) + shunt_bus_new = jnp.where(shunt_bus_old == bus_to_split, new_bus_b_index, shunt_bus_old) base_bus_indices = jnp.array([bus_to_split, new_bus_b_index], dtype=jnp.int32) bus_candidates = jnp.concatenate( @@ -98,6 +122,8 @@ def _compute_bsdf_update_impl( branch_to_old, branch_from_new, branch_to_new, + shunt_bus_old, + shunt_bus_new, ], axis=0, ) @@ -141,6 +167,11 @@ def _gather_positions(indices: jnp.ndarray) -> Tuple[jnp.ndarray, jnp.ndarray]: mag_from_new_idx = jnp.take(magnitude_component_indices, branch_from_new, axis=0) mag_to_new_idx = jnp.take(magnitude_component_indices, branch_to_new, axis=0) + theta_shunt_old_idx = jnp.take(angle_component_indices, shunt_bus_old, axis=0) + mag_shunt_old_idx = jnp.take(magnitude_component_indices, shunt_bus_old, axis=0) + theta_shunt_new_idx = jnp.take(angle_component_indices, shunt_bus_new, axis=0) + mag_shunt_new_idx = jnp.take(magnitude_component_indices, shunt_bus_new, axis=0) + old_indices = jnp.stack( [ theta_from_old_idx, @@ -161,8 +192,26 @@ def _gather_positions(indices: jnp.ndarray) -> Tuple[jnp.ndarray, jnp.ndarray]: axis=1, ) + shunt_old_indices = jnp.stack( + [ + theta_shunt_old_idx, + mag_shunt_old_idx, + ], + axis=1, + ) + + shunt_new_indices = jnp.stack( + [ + theta_shunt_new_idx, + mag_shunt_new_idx, + ], + axis=1, + ) + old_pos, old_valid = _gather_positions(old_indices) new_pos, new_valid = _gather_positions(new_indices) + shunt_old_pos, shunt_old_valid = _gather_positions(shunt_old_indices) + shunt_new_pos, shunt_new_valid = _gather_positions(shunt_new_indices) vm_from_old = jnp.take(v_mag_hat, branch_from_old, axis=0) vm_to_old = jnp.take(v_mag_hat, branch_to_old, axis=0) @@ -178,6 +227,7 @@ def _gather_positions(indices: jnp.ndarray) -> Tuple[jnp.ndarray, jnp.ndarray]: y_ft_sel = jnp.take(y_ft, branches_connected_to_bus_b, axis=0) y_tf_sel = jnp.take(y_tf, branches_connected_to_bus_b, axis=0) y_tt_sel = jnp.take(y_tt, branches_connected_to_bus_b, axis=0) + y_shunt_sel = jnp.take(y_shunt, shunt_connected_to_bus_b, axis=0) compute_delta = jax.vmap(_compute_branch_delta_submatrix_from_admittance) delta_old = compute_delta( @@ -201,6 +251,16 @@ def _gather_positions(indices: jnp.ndarray) -> Tuple[jnp.ndarray, jnp.ndarray]: y_tt_sel, ) + compute_shunt_delta = jax.vmap(_compute_shunt_delta_submatrix_from_admittance) + delta_shunt_old = compute_shunt_delta( + jnp.take(v_mag_hat, shunt_bus_old, axis=0), + y_shunt_sel, + ) + delta_shunt_new = compute_shunt_delta( + jnp.take(v_mag_hat, shunt_bus_new, axis=0), + y_shunt_sel, + ) + k_shape = (k_max, k_max) delta_block = jnp.zeros(k_shape, dtype=dtype) @@ -223,6 +283,8 @@ def _accumulate( delta_block = _accumulate(delta_old, old_pos, old_valid, 1.0, delta_block) delta_block = _accumulate(delta_new, new_pos, new_valid, -1.0, delta_block) + delta_block = _accumulate(delta_shunt_old, shunt_old_pos, shunt_old_valid, 1.0, delta_block) + delta_block = _accumulate(delta_shunt_new, shunt_new_pos, shunt_new_valid, -1.0, delta_block) theta_new_idx = angle_component_indices[new_bus_b_index] mag_new_idx = magnitude_component_indices[new_bus_b_index] @@ -274,18 +336,17 @@ def compute_bsdf_update( 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 dense inverse transpose update for BSDF. + """Inverse transpose update for BSDF. This function computes the updated Jacobian inverse transpose after a bus split with branch re-attachments, using a full-rank update approach. It is intended for reference and testing purposes, and is not optimized for performance. Note: currently not supported: - - Shunt reassignments - Changes in branch parameters (e.g., series admittance, taps, phase shifts) (This would involve computing the different delta_blocks) - Changes in the type of the new bus (e.g., PQ, PV, slack) - - injection reassignments + - Injection reassignments of regulating elements, resulting in a new PV bus Parameters ---------- @@ -339,14 +400,17 @@ def compute_bsdf_update( bus_to_split=int(bus_to_split), new_bus_b_index=int(new_bus_b_index), branches_connected_to_bus_b=jnp.asarray(branches_connected_to_bus_b, dtype=jnp.int32), + shunt_connected_to_bus_b=jnp.asarray(shunt_connected_to_bus_b, dtype=jnp.int32), branch_from=jnp.asarray(branch_from, dtype=jnp.int32), branch_to=jnp.asarray(branch_to, dtype=jnp.int32), + shunt_to_bus=jnp.asarray(shunt_to_bus, dtype=jnp.int32), v_mag_hat=jnp.asarray(v_mag_hat, dtype=real_dtype), theta_hat=jnp.asarray(theta_hat, dtype=real_dtype), y_ff=jnp.asarray(y_ff, dtype=jnp.complex128), y_ft=jnp.asarray(y_ft, dtype=jnp.complex128), y_tf=jnp.asarray(y_tf, dtype=jnp.complex128), y_tt=jnp.asarray(y_tt, dtype=jnp.complex128), + y_shunt=jnp.asarray(y_shunt, dtype=jnp.complex128), angle_component_indices=jnp.asarray(angle_component_indices, dtype=jnp.int32), magnitude_component_indices=jnp.asarray(magnitude_component_indices, dtype=jnp.int32), ) diff --git a/tests/jax/test_bsdf_jax.py b/tests/jax/test_bsdf_jax.py index 6e0853c..8ed5703 100644 --- a/tests/jax/test_bsdf_jax.py +++ b/tests/jax/test_bsdf_jax.py @@ -11,8 +11,10 @@ import numpy as np import pytest +from dc_plus.importing.powsybl.powsybl_import import DANGLING_BUS_STRING_SUFFIX from dc_plus.jax.bsdf import compute_bsdf_update as compute_bsdf_update_jax from tests.test_helper.bsdf_helper import ( + derive_bus_order, get_bsdf_cases, prepare_bsdf_test_context, run_reference_one_step, @@ -33,8 +35,8 @@ def test_bsdf_full_rank_jax(bsdf_test_case): bus_to_split=setup.bus_to_split, new_bus_b_index=setup.new_bus_index, new_bus_type=2, - branches_connected_to_bus_b=jnp.asarray(setup.branches_to_move, dtype=jnp.int32), - shunt_connected_to_bus_b=jnp.asarray([], dtype=jnp.int32), + branches_connected_to_bus_b=jnp.asarray(setup.branches_connected_to_bus_b, dtype=jnp.int32), + shunt_connected_to_bus_b=jnp.asarray(setup.shunt_connected_to_bus_b, dtype=jnp.int32), branch_from=setup.dynamic_info.branch_from_bus, branch_to=setup.dynamic_info.branch_to_bus, shunt_to_bus=setup.dynamic_info.shunt_bus_indices, @@ -67,9 +69,9 @@ def test_bsdf_full_rank_jax(bsdf_test_case): ) dynamic_info_split_manual = setup.dynamic_info_split_manual - dynamic_info_one_step = run_reference_one_step(setup.net, bsdf_test_case=bsdf_test_case) - bus_order = bsdf_test_case.bus_order - dx = -jacobian_inv_jax @ setup.mismatch_n1 + dynamic_info_one_step, string_info_one_step = run_reference_one_step(setup.net, bsdf_test_case=bsdf_test_case) + bus_order = derive_bus_order(setup.split_bus_ids, string_info_one_step.bus_ids) + dx = -jacobian_inv_jax @ setup.mismatch_bsdf_reference theta_actual = dynamic_info_split_manual.bus_voltage_angles_rad vm_actual = dynamic_info_split_manual.bus_voltage_magnitudes @@ -80,16 +82,31 @@ def test_bsdf_full_rank_jax(bsdf_test_case): theta_updated_J[pvpq] = theta_actual[pvpq] + dx[setup.jacobian_data_with_extra_buses.is_angle_component] vm_updated_J[pq] = vm_actual[pq] + dx[setup.jacobian_data_with_extra_buses.is_magnitude_component] + dangling_bus_mask = np.char.endswith(np.asarray(string_info_one_step.bus_ids, dtype=str), DANGLING_BUS_STRING_SUFFIX) + # reorder bus voltages to match new ordering of manual split np.testing.assert_allclose( - dynamic_info_one_step.bus_voltage_magnitudes[bus_order], - vm_updated_J, + dynamic_info_one_step.bus_voltage_magnitudes[~dangling_bus_mask], + vm_updated_J[bus_order][~dangling_bus_mask], rtol=1e-10, atol=1e-10, ) np.testing.assert_allclose( - dynamic_info_one_step.bus_voltage_angles_rad[bus_order], - theta_updated_J, + dynamic_info_one_step.bus_voltage_angles_rad[~dangling_bus_mask], + theta_updated_J[bus_order][~dangling_bus_mask], rtol=1e-10, atol=1e-10, ) + + np.testing.assert_allclose( + dynamic_info_one_step.bus_voltage_magnitudes[dangling_bus_mask], + vm_updated_J[bus_order][dangling_bus_mask], + rtol=1e-8, + atol=1e-8, + ) + np.testing.assert_allclose( + dynamic_info_one_step.bus_voltage_angles_rad[dangling_bus_mask], + theta_updated_J[bus_order][dangling_bus_mask], + rtol=1e-8, + atol=1e-8, + ) From 6daf5f2101d0ef2957dcc67bfd77bb6ee28387dc Mon Sep 17 00:00:00 2001 From: Bpetrick <170433522+BenjPetr@users.noreply.github.com> Date: Thu, 16 Apr 2026 09:02:45 +0000 Subject: [PATCH 3/4] chore: fix typehint and docstr for shunts --- src/dc_plus/jax/bsdf.py | 4 ++-- src/dc_plus/numpy/bsdf_full_rank.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/dc_plus/jax/bsdf.py b/src/dc_plus/jax/bsdf.py index e13f66a..b761743 100644 --- a/src/dc_plus/jax/bsdf.py +++ b/src/dc_plus/jax/bsdf.py @@ -332,7 +332,7 @@ def compute_bsdf_update( y_ft: Complex128[jnp.ndarray, " n_branches"], y_tf: Complex128[jnp.ndarray, " n_branches"], y_tt: Complex128[jnp.ndarray, " n_branches"], - y_shunt: Complex128[jnp.ndarray, " n_buses"], + y_shunt: Complex128[jnp.ndarray, " n_shunts"], 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"]: @@ -380,7 +380,7 @@ def compute_bsdf_update( "To-From" admittance for all branches. y_tt : Complex128[jnp.ndarray, " n_branches"] "To-To" admittance for all branches. - y_shunt : Complex128[jnp.ndarray, " n_buses"] + y_shunt : Complex128[jnp.ndarray, " n_shunts"] Shunt admittance for all buses. angle_component_indices : Int[jnp.ndarray, " n_eq_jacobian"] Mapping from bus indices to angle component indices in the Jacobian. diff --git a/src/dc_plus/numpy/bsdf_full_rank.py b/src/dc_plus/numpy/bsdf_full_rank.py index 46efc43..d4ed824 100644 --- a/src/dc_plus/numpy/bsdf_full_rank.py +++ b/src/dc_plus/numpy/bsdf_full_rank.py @@ -606,7 +606,7 @@ def compute_bsdf_update( y_ft: Complex128[np.ndarray, " n_branches"], y_tf: Complex128[np.ndarray, " n_branches"], y_tt: Complex128[np.ndarray, " n_branches"], - y_shunt: Complex128[np.ndarray, " n_buses"], + y_shunt: Complex128[np.ndarray, " n_shunts"], angle_component_indices: Int[np.ndarray, " n_eq_jacobian"], magnitude_component_indices: Int[np.ndarray, " n_eq_jacobian"], ) -> Float[np.ndarray, " n_eq n_eq"]: From ba94d7f3cd2411c1c4d23445e3ab435b349133fe Mon Sep 17 00:00:00 2001 From: Bpetrick <170433522+BenjPetr@users.noreply.github.com> Date: Thu, 16 Apr 2026 09:05:13 +0000 Subject: [PATCH 4/4] fix: add missing consistency check --- src/dc_plus/interfaces/network_information.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/dc_plus/interfaces/network_information.py b/src/dc_plus/interfaces/network_information.py index 6ee5bf6..371ca31 100644 --- a/src/dc_plus/interfaces/network_information.py +++ b/src/dc_plus/interfaces/network_information.py @@ -739,3 +739,6 @@ def _check_network_data_consistency( assert string_network_data.injection_types.shape[0] == dynamic_network_data.injection_to_bus.shape[0], ( "Inconsistent number of injections between injection_types and injection_to_bus." ) + assert string_network_data.injection_ids.shape[0] == dynamic_network_data.injection_to_bus.shape[0], ( + "Inconsistent number of injections between injection_ids and injection_to_bus." + )