feat: draft injection outage

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 injections, one at each terminal.

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,

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."
     )

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,
             }
         )

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,17 +908,27 @@ 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
 
 
-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.
@@ -930,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
     -------
@@ -942,6 +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["id_str"].eq(reference_id)
 
     buses.rename(columns={"v_mag": "voltage_magnitude", "v_angle": "voltage_angle"}, inplace=True)
 
@@ -966,15 +980,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

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]

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]
     )

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,
@@ -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)

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], (

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)