mambo_power.numerics

Network matrices over scipy.sparse. See the manual page for the formulas and conventions.

mambo_power.numerics

Network matrices over scipy.sparse: the positional pu view, Ybus, Bbus, PTDF and LODF.

The only package module that holds positional indices. NetworkArrays is the single per-unit conversion site; every builder here takes that view, never a Network directly. effective_roles derives the roles a solver must use from the declared ones (W3).

NetworkArrays dataclass

NetworkArrays(
    base_mva: float,
    bus_ids: list[str],
    bus_index: dict[str, int],
    n_bus: int,
    slack: int,
    bus_type: IntArray,
    branch_ids: list[str],
    branch_index: dict[str, int],
    n_branch: int,
    f: IntArray,
    t: IntArray,
    r: FloatArray,
    x: FloatArray,
    b: FloatArray,
    tap: FloatArray,
    shift_rad: FloatArray,
    rating_pu: FloatArray,
    p_load_pu: FloatArray,
    q_load_pu: FloatArray,
    g_shunt_pu: FloatArray,
    b_shunt_pu: FloatArray,
    p_gen_pu: FloatArray,
    q_gen_pu: FloatArray,
    p_min_pu: FloatArray,
    p_max_pu: FloatArray,
    q_min_pu: FloatArray,
    q_max_pu: FloatArray,
    v_set: FloatArray,
    gen_ids: list[str] = list(),
    gen_bus: IntArray = (lambda: zeros(0, dtype=int64))(),
    gen_p_pu: FloatArray = (lambda: zeros(0))(),
    gen_q_pu: FloatArray = (lambda: zeros(0))(),
    gen_p_min_pu: FloatArray = (lambda: zeros(0))(),
    gen_p_max_pu: FloatArray = (lambda: zeros(0))(),
    gen_q_min_pu: FloatArray = (lambda: zeros(0))(),
    gen_q_max_pu: FloatArray = (lambda: zeros(0))(),
    gen_v_set: FloatArray = (lambda: zeros(0))(),
    load_ids: list[str] = list(),
    load_bus: IntArray = (lambda: zeros(0, dtype=int64))(),
    load_p_min_pu: FloatArray = (lambda: zeros(0))(),
    load_p_max_pu: FloatArray = (lambda: zeros(0))(),
    storage_ids: list[str] = list(),
    storage_bus: IntArray = (
        lambda: zeros(0, dtype=int64)
    )(),
    storage_p_max_pu: FloatArray = (lambda: zeros(0))(),
    storage_energy_pu: FloatArray = (lambda: zeros(0))(),
    storage_soc_initial: FloatArray = (lambda: zeros(0))(),
    storage_efficiency_charge: FloatArray = (
        lambda: zeros(0)
    )(),
    storage_efficiency_discharge: FloatArray = (
        lambda: zeros(0)
    )(),
)

Frozen positional arrays over the in-service subset of a network, in per unit.

Positions are 0-based. bus_ids[i] is the id at bus position i; branch_ids[k] and gen_ids[g] likewise. Order follows the network's collection order with the excluded elements removed.

bus_type instance-attribute

bus_type: IntArray

1 = pq, 2 = pv, 3 = slack, as declared on the bus.

b instance-attribute

b: FloatArray

Total line charging susceptance per branch (pu); builders apply b / 2 per end.

tap instance-attribute

tap: FloatArray

Tap ratio magnitude on the from side; 1.0 where the branch has none.

rating_pu instance-attribute

rating_pu: FloatArray

Thermal rating in pu of base_mva; inf where the branch has none.

g_shunt_pu instance-attribute

g_shunt_pu: FloatArray

Shunt conductance per bus in pu (MATPOWER GS sign: positive consumes).

b_shunt_pu instance-attribute

b_shunt_pu: FloatArray

Shunt susceptance per bus in pu (MATPOWER BS sign: positive injects).

v_set instance-attribute

v_set: FloatArray

First in-service generator's v_set_pu at each bus; 1.0 where there is none.

load_p_max_pu class-attribute instance-attribute

load_p_max_pu: FloatArray = field(
    default_factory=lambda: np.zeros(0)
)

Per-load [0, p_mw] bound in pu: the natural bound for a bid-load's demand dispatch is zero up to its own fixed historical p_mw. Built uniformly for every in-service load regardless of Load.bid — Load carries no p_min_mw/p_max_mw fields to mirror the generator-side source data, and this bound formula does not depend on bid presence; whether/how a given load's bound is actually used by opf.dc_opf is decided per-load elsewhere, not here. p_load_pu/q_load_pu (the bus aggregate) are untouched by this addition.

storage_p_max_pu class-attribute instance-attribute

storage_p_max_pu: FloatArray = field(
    default_factory=lambda: np.zeros(0)
)

Charge/discharge power limit (Storage.p_max_mw), pu of base_mva.

storage_energy_pu class-attribute instance-attribute

storage_energy_pu: FloatArray = field(
    default_factory=lambda: np.zeros(0)
)

Energy capacity (Storage.energy_mwh), pu of base_mva (pu-hours) — the same base_mva-division convention every other physical field in this class already uses (ADR-005: physical units in the model, pu in numerics).

storage_soc_initial class-attribute instance-attribute

storage_soc_initial: FloatArray = field(
    default_factory=lambda: np.zeros(0)
)

Initial state of charge, already a fraction of energy_mwh in [0, 1] on the entity — dimensionless, so unlike storage_p_max_pu/storage_energy_pu it carries through unconverted.

storage_efficiency_discharge class-attribute instance-attribute

storage_efficiency_discharge: FloatArray = field(
    default_factory=lambda: np.zeros(0)
)

Charge/discharge efficiency, already dimensionless ratios in (0, 1] — unconverted.

from_network classmethod

from_network(net: Network) -> NetworkArrays

Build the in-service positional view; the one pu-conversion site.

Source code in src/mambo_power/numerics/arrays.py

@classmethod
def from_network(cls, net: Network) -> NetworkArrays:
    """Build the in-service positional view; the one pu-conversion site."""
    base = float(net.base_mva)

    live_buses = [bus for bus in net.buses if bus.in_service]
    bus_ids = [bus.id for bus in live_buses]
    bus_index = {bus_id: i for i, bus_id in enumerate(bus_ids)}
    n_bus = len(bus_ids)
    slack_positions = [i for i, bus in enumerate(live_buses) if bus.type == "slack"]
    if len(slack_positions) != 1:
        raise ValueError(
            f"expected exactly one in-service slack bus, found {len(slack_positions)}"
        )
    bus_type = np.fromiter(
        (BUS_TYPE_CODE[bus.type] for bus in live_buses), dtype=np.int64, count=n_bus
    )

    branches = [
        br
        for br in net.branches
        if br.in_service and br.from_bus in bus_index and br.to_bus in bus_index
    ]
    n_branch = len(branches)
    branch_ids = [br.id for br in branches]
    branch_index = {br_id: k for k, br_id in enumerate(branch_ids)}
    f = np.fromiter((bus_index[br.from_bus] for br in branches), np.int64, n_branch)
    t = np.fromiter((bus_index[br.to_bus] for br in branches), np.int64, n_branch)
    r = np.fromiter((br.r for br in branches), np.float64, n_branch)
    x = np.fromiter((br.x for br in branches), np.float64, n_branch)
    b = np.fromiter((br.b for br in branches), np.float64, n_branch)
    tap = np.fromiter(
        (1.0 if br.tap_ratio is None else br.tap_ratio for br in branches),
        np.float64,
        n_branch,
    )
    shift_rad = np.fromiter(
        (0.0 if br.shift_deg is None else math.radians(br.shift_deg) for br in branches),
        np.float64,
        n_branch,
    )
    rating_pu = np.fromiter(
        (math.inf if br.rating_mva is None else br.rating_mva / base for br in branches),
        np.float64,
        n_branch,
    )

    def per_bus(pairs: list[tuple[int, float]]) -> FloatArray:
        positions = np.fromiter((p for p, _ in pairs), np.int64, len(pairs))
        values = np.fromiter((v for _, v in pairs), np.float64, len(pairs))
        summed = np.asarray(
            np.bincount(positions, weights=values, minlength=n_bus), dtype=np.float64
        )
        return summed / base

    loads = [ld for ld in net.loads if ld.in_service and ld.bus in bus_index]
    p_load_pu = per_bus([(bus_index[ld.bus], ld.p_mw) for ld in loads])
    q_load_pu = per_bus([(bus_index[ld.bus], ld.q_mvar) for ld in loads])

    n_load = len(loads)
    load_bus = np.fromiter((bus_index[ld.bus] for ld in loads), np.int64, n_load)

    def per_load(values: list[float]) -> FloatArray:
        return np.fromiter(values, np.float64, n_load) / base

    shunts = [sh for sh in net.shunts if sh.in_service and sh.bus in bus_index]
    g_shunt_pu = per_bus([(bus_index[sh.bus], sh.g_mw) for sh in shunts])
    b_shunt_pu = per_bus([(bus_index[sh.bus], sh.b_mvar) for sh in shunts])

    gens = [g for g in net.generators if g.in_service and g.bus in bus_index]
    n_gen = len(gens)
    gen_bus = np.fromiter((bus_index[g.bus] for g in gens), np.int64, n_gen)
    p_gen_pu = per_bus([(bus_index[g.bus], g.p_mw) for g in gens])
    q_gen_pu = per_bus([(bus_index[g.bus], g.q_mvar) for g in gens])
    p_min_pu = per_bus([(bus_index[g.bus], g.p_min_mw) for g in gens])
    p_max_pu = per_bus([(bus_index[g.bus], g.p_max_mw) for g in gens])
    q_min_pu = per_bus([(bus_index[g.bus], g.q_min_mvar) for g in gens])
    q_max_pu = per_bus([(bus_index[g.bus], g.q_max_mvar) for g in gens])
    v_set = np.ones(n_bus)
    seen: set[int] = set()
    for g in gens:
        position = bus_index[g.bus]
        if position not in seen:
            seen.add(position)
            v_set[position] = g.v_set_pu

    def per_gen(values: list[float]) -> FloatArray:
        return np.fromiter(values, np.float64, n_gen) / base

    storage_units = [s for s in net.storage if s.in_service and s.bus in bus_index]
    n_storage = len(storage_units)
    storage_bus = np.fromiter((bus_index[s.bus] for s in storage_units), np.int64, n_storage)

    def per_storage(values: list[float]) -> FloatArray:
        return np.fromiter(values, np.float64, n_storage)

    return cls(
        base_mva=base,
        bus_ids=bus_ids,
        bus_index=bus_index,
        n_bus=n_bus,
        slack=slack_positions[0],
        bus_type=bus_type,
        branch_ids=branch_ids,
        branch_index=branch_index,
        n_branch=n_branch,
        f=f,
        t=t,
        r=r,
        x=x,
        b=b,
        tap=tap,
        shift_rad=shift_rad,
        rating_pu=rating_pu,
        p_load_pu=p_load_pu,
        q_load_pu=q_load_pu,
        g_shunt_pu=g_shunt_pu,
        b_shunt_pu=b_shunt_pu,
        p_gen_pu=p_gen_pu,
        q_gen_pu=q_gen_pu,
        p_min_pu=p_min_pu,
        p_max_pu=p_max_pu,
        q_min_pu=q_min_pu,
        q_max_pu=q_max_pu,
        v_set=v_set,
        gen_ids=[g.id for g in gens],
        gen_bus=gen_bus,
        gen_p_pu=per_gen([g.p_mw for g in gens]),
        gen_q_pu=per_gen([g.q_mvar for g in gens]),
        gen_p_min_pu=per_gen([g.p_min_mw for g in gens]),
        gen_p_max_pu=per_gen([g.p_max_mw for g in gens]),
        gen_q_min_pu=per_gen([g.q_min_mvar for g in gens]),
        gen_q_max_pu=per_gen([g.q_max_mvar for g in gens]),
        gen_v_set=np.fromiter((g.v_set_pu for g in gens), np.float64, n_gen),
        load_ids=[ld.id for ld in loads],
        load_bus=load_bus,
        load_p_min_pu=np.zeros(n_load),
        load_p_max_pu=per_load([ld.p_mw for ld in loads]),
        storage_ids=[s.id for s in storage_units],
        storage_bus=storage_bus,
        storage_p_max_pu=per_storage([s.p_max_mw / base for s in storage_units]),
        storage_energy_pu=per_storage([s.energy_mwh / base for s in storage_units]),
        storage_soc_initial=per_storage([s.soc_initial for s in storage_units]),
        storage_efficiency_charge=per_storage([s.efficiency_charge for s in storage_units]),
        storage_efficiency_discharge=per_storage(
            [s.efficiency_discharge for s in storage_units]
        ),
    )

NoSlackGeneratorError

NoSlackGeneratorError(bus_id: str, position: int)

Bases: Exception

The slack bus has no in-service generator, so nothing can close the power balance.

Raised by mambo_power.numerics.effective_roles. MATPOWER's bustypes would silently hand the reference role to the first PV bus; the M2 spec rejects that re-slacking (Not Doing: "MATPOWER slack-limiting re-slack") and names the condition instead. bus_id is the slack bus id, position its index in the arrays.

Source code in src/mambo_power/numerics/errors.py

def __init__(self, bus_id: str, position: int) -> None:
    self.bus_id = bus_id
    self.position = position
    super().__init__(
        f'slack bus "{bus_id}" (position {position}) has no in-service generator; '
        "a power flow cannot close the balance"
    )

SetpointConflictWarning

Bases: UserWarning

Several in-service generators at one bus carry different voltage setpoints.

Emitted by mambo_power.numerics.effective_roles (via warnings.warn). The last generator's setpoint is used, following MATPOWER; pandapower raises a UserWarning and aborts in the same situation, which is why this is surfaced rather than resolved silently.

UnsolvableNetworkError

Bases: Exception

A Network that passes validate_network but cannot be solved by the numerics it was handed to — e.g. DC susceptance is undefined when a branch carries x == 0 with r != 0 (legal under the model, since BAD_RANGE only rejects r == x == 0). This is user data, not a solver bug: mambo_power.jobs.run maps it to the structured UNSOLVABLE_NETWORK failure code rather than INTERNAL.

EffectiveRoles dataclass

EffectiveRoles(
    bus_type: IntArray,
    v_set: FloatArray,
    demoted_pv: IntArray,
    setpoint_conflicts: list[
        tuple[str, list[str], list[float]]
    ],
)

The roles and setpoints a solver must use, positional over the arrays' bus order.

bus_type instance-attribute

bus_type: IntArray

Effective role per bus: 1 = pq, 2 = pv, 3 = slack (same codes as the arrays).

v_set instance-attribute

v_set: FloatArray

Effective voltage setpoint per bus (pu): last in-service generator's; 1.0 if none.

demoted_pv instance-attribute

demoted_pv: IntArray

Positions of buses declared PV but solved as PQ (no in-service generator), ascending.

setpoint_conflicts instance-attribute

setpoint_conflicts: list[tuple[str, list[str], list[float]]]

(bus_id, gen_ids, setpoints) for every bus whose in-service generators disagree.

bf

bf(arr: NetworkArrays) -> Any

Bf: n_branch × n_bus CSC matrix; Bf @ θ is the from-side DC flow (pu).

Source code in src/mambo_power/numerics/bbus.py

def bf(arr: NetworkArrays) -> Any:
    """``Bf``: ``n_branch × n_bus`` CSC matrix; ``Bf @ θ`` is the from-side DC flow (pu)."""
    b = branch_susceptance(arr)
    rows = np.concatenate([np.arange(arr.n_branch), np.arange(arr.n_branch)])
    cols = np.concatenate([arr.f, arr.t])
    data = np.concatenate([b, -b])
    return sparse.csc_matrix((data, (rows, cols)), shape=(arr.n_branch, arr.n_bus))

flow_from_ptdf

flow_from_ptdf(
    ptdf: FloatArray,
    injection_mw: FloatArray,
    arr: NetworkArrays,
) -> FloatArray

Branch flow, MW, from a PTDF matrix and a full bus net-injection vector, MW.

flow = ptdf @ (injection_mw − p_shift·base_mva) + pf_shift·base_mva — the phase-shifter injection is subtracted out of the bus injection before the PTDF product, then each branch's own from-side shift flow is added back on. This is exactly :func:mambo_power.pf.dc. solve's construction (its module docstring: rhs = P − p_shift feeds the angle solve, and p_from = Bf·θ + pf_shift; combined with θ = B'⁻¹(P − p_shift) and PTDF = Bf·B'⁻¹ on the reduced system, p_from = PTDF·(P − p_shift) + pf_shift) — the model every DC PTDF-based flow in this package must match. Omitting the − p_shift term (all of opf.dc_opf, opf.solve_dc_opf and market._clearing did until M8 finding F1 / A19) reproduces pf.solve_dc's flow only when no branch has a shift, since p_shift(arr) == 0 identically in that case.

injection_mw must be the full net injection per bus (generation minus load minus shunt, MW) — callers that instead fold some of that into a decision-variable-relative LP constant (opf.dc_opf's own flow-limit rows) derive the identical correction by hand rather than calling this helper, since their injection is not one vector (see that module's own derivation).

Source code in src/mambo_power/numerics/bbus.py

def flow_from_ptdf(ptdf: FloatArray, injection_mw: FloatArray, arr: NetworkArrays) -> FloatArray:
    """Branch flow, MW, from a PTDF matrix and a full bus net-injection vector, MW.

    ``flow = ptdf @ (injection_mw − p_shift·base_mva) + pf_shift·base_mva`` — the phase-shifter
    injection is subtracted out of the bus injection *before* the PTDF product, then each
    branch's own from-side shift flow is added back on. This is exactly :func:`mambo_power.pf.dc.
    solve`'s construction (its module docstring: ``rhs = P − p_shift`` feeds the angle solve, and
    ``p_from = Bf·θ + pf_shift``; combined with ``θ = B'⁻¹(P − p_shift)`` and ``PTDF = Bf·B'⁻¹``
    on the reduced system, ``p_from = PTDF·(P − p_shift) + pf_shift``) — the model every DC
    PTDF-based flow in this package must match. Omitting the ``− p_shift`` term (all of
    ``opf.dc_opf``, ``opf.solve_dc_opf`` and ``market._clearing`` did until M8 finding F1 / A19)
    reproduces ``pf.solve_dc``'s flow only when no branch has a shift, since
    ``p_shift(arr) == 0`` identically in that case.

    ``injection_mw`` must be the *full* net injection per bus (generation minus load minus
    shunt, MW) — callers that instead fold some of that into a decision-variable-relative LP
    constant (``opf.dc_opf``'s own flow-limit rows) derive the identical correction by hand
    rather than calling this helper, since their ``injection`` is not one vector (see that
    module's own derivation).
    """
    result: FloatArray = (
        ptdf @ (injection_mw - p_shift(arr) * arr.base_mva) + pf_shift(arr) * arr.base_mva
    )
    return result

p_shift

p_shift(arr: NetworkArrays) -> FloatArray

Per-bus phase-shifter injection Cftᵀ · pf_shift (pu); P = Bbus·θ + p_shift.

Source code in src/mambo_power/numerics/bbus.py

def p_shift(arr: NetworkArrays) -> FloatArray:
    """Per-bus phase-shifter injection ``Cftᵀ · pf_shift`` (pu); ``P = Bbus·θ + p_shift``."""
    result: FloatArray = np.asarray(incidence(arr).T @ pf_shift(arr), dtype=np.float64).ravel()
    return result

bridges

bridges(arr: NetworkArrays) -> list[int]

Positions of branches whose removal disconnects the in-service graph (sorted).

Iterative Tarjan lowpoint search over the multigraph; parallel branches between the same pair of buses are never bridges because the search skips only the edge it arrived by.

Source code in src/mambo_power/numerics/lodf.py

def bridges(arr: NetworkArrays) -> list[int]:
    """Positions of branches whose removal disconnects the in-service graph (sorted).

    Iterative Tarjan lowpoint search over the multigraph; parallel branches between the same
    pair of buses are never bridges because the search skips only the *edge* it arrived by.
    """
    n_bus, n_branch = arr.n_bus, arr.n_branch
    adjacency: list[list[tuple[int, int]]] = [[] for _ in range(n_bus)]
    for k in range(n_branch):
        u, v = int(arr.f[k]), int(arr.t[k])
        adjacency[u].append((v, k))
        adjacency[v].append((u, k))

    disc = [-1] * n_bus
    low = [0] * n_bus
    found: list[int] = []
    clock = 0
    for root in range(n_bus):
        if disc[root] != -1:
            continue
        disc[root] = low[root] = clock
        clock += 1
        # stack entries: (node, edge used to enter it, next adjacency cursor)
        stack: list[tuple[int, int, int]] = [(root, -1, 0)]
        while stack:
            node, via, cursor = stack[-1]
            if cursor < len(adjacency[node]):
                stack[-1] = (node, via, cursor + 1)
                nxt, edge = adjacency[node][cursor]
                if edge == via:
                    continue
                if disc[nxt] == -1:
                    disc[nxt] = low[nxt] = clock
                    clock += 1
                    stack.append((nxt, edge, 0))
                else:
                    low[node] = min(low[node], disc[nxt])
            else:
                stack.pop()
                if stack:
                    parent = stack[-1][0]
                    low[parent] = min(low[parent], low[node])
                    if low[node] > disc[parent]:
                        found.append(via)
    found.sort()
    return found

effective_roles

effective_roles(arr: NetworkArrays) -> EffectiveRoles

Derive the effective roles and setpoints from arr (see the module docstring).

Raises NoSlackGeneratorError when the slack bus has no in-service generator. Emits one SetpointConflictWarning per bus whose generators disagree.

Source code in src/mambo_power/numerics/roles.py

def effective_roles(arr: NetworkArrays) -> EffectiveRoles:
    """Derive the effective roles and setpoints from ``arr`` (see the module docstring).

    Raises :class:`NoSlackGeneratorError` when the slack bus has no in-service generator.
    Emits one :class:`SetpointConflictWarning` per bus whose generators disagree.
    """
    pq, pv, slack = BUS_TYPE_CODE["pq"], BUS_TYPE_CODE["pv"], BUS_TYPE_CODE["slack"]
    n_gen = len(arr.gen_ids)
    gens_at: dict[int, list[int]] = {}
    for g in range(n_gen):
        gens_at.setdefault(int(arr.gen_bus[g]), []).append(g)

    bus_type = arr.bus_type.copy()
    has_gen = np.zeros(arr.n_bus, dtype=bool)
    has_gen[list(gens_at)] = True

    if bus_type[arr.slack] == slack and not has_gen[arr.slack]:
        raise NoSlackGeneratorError(arr.bus_ids[arr.slack], int(arr.slack))

    demoted = np.flatnonzero((bus_type == pv) & ~has_gen).astype(np.int64)
    bus_type[demoted] = pq

    v_set = np.ones(arr.n_bus)
    conflicts: list[tuple[str, list[str], list[float]]] = []
    for position in sorted(gens_at):
        rows = gens_at[position]
        setpoints = [float(arr.gen_v_set[g]) for g in rows]
        v_set[position] = setpoints[-1]
        if len(rows) > 1 and any(
            not math.isclose(s, setpoints[-1], rel_tol=0.0, abs_tol=SETPOINT_TOL) for s in setpoints
        ):
            bus_id = arr.bus_ids[position]
            gen_ids = [arr.gen_ids[g] for g in rows]
            conflicts.append((bus_id, gen_ids, setpoints))
            pairs = ", ".join(f"{gid}={s:g}" for gid, s in zip(gen_ids, setpoints, strict=True))
            warnings.warn(
                SetpointConflictWarning(
                    f'bus "{bus_id}": in-service generators disagree on the voltage setpoint '
                    f"({pairs}); using the last one, {setpoints[-1]:g} pu (MATPOWER rule)"
                ),
                stacklevel=2,
            )

    return EffectiveRoles(
        bus_type=bus_type, v_set=v_set, demoted_pv=demoted, setpoint_conflicts=conflicts
    )

yf_yt

yf_yt(arr: NetworkArrays) -> tuple[Any, Any]

(Yf, Yt): n_branch × n_bus complex CSC matrices giving from/to branch currents.

Yf @ V is the current injected into each branch at its from bus, Yt @ V at its to bus — the inputs M2 needs for branch flows.

Source code in src/mambo_power/numerics/ybus.py

def yf_yt(arr: NetworkArrays) -> tuple[Any, Any]:
    """``(Yf, Yt)``: ``n_branch × n_bus`` complex CSC matrices giving from/to branch currents.

    ``Yf @ V`` is the current injected into each branch at its from bus, ``Yt @ V`` at its to
    bus — the inputs M2 needs for branch flows.
    """
    yff, yft, ytf, ytt = branch_admittances(arr)
    rows = np.concatenate([np.arange(arr.n_branch), np.arange(arr.n_branch)])
    cols = np.concatenate([arr.f, arr.t])
    shape = (arr.n_branch, arr.n_bus)
    yf = sparse.csc_matrix((np.concatenate([yff, yft]), (rows, cols)), shape=shape)
    yt = sparse.csc_matrix((np.concatenate([ytf, ytt]), (rows, cols)), shape=shape)
    return yf, yt

NetworkArrays

mambo_power.numerics.arrays

NetworkArrays: the positional, per-unit view of a Network.

This is the only place in the package that holds positional indices and the single site where physical units (MW, MVAr, MVA) are divided by base_mva (wave M1 design items 1 and 7). Every matrix builder in mambo_power.numerics consumes this view; nothing else in the package divides by base_mva.

Scope: the in-service subset. Out-of-service buses are dropped, and so is every branch, generator, load or shunt that is itself out of service or attached to a dropped bus. The network's own validation guarantees the surviving buses form one connected component with exactly one slack.

BUS_TYPE_CODE module-attribute

BUS_TYPE_CODE = {'pq': 1, 'pv': 2, 'slack': 3}

MATPOWER bus type codes used in NetworkArrays.bus_type.

NetworkArrays dataclass

NetworkArrays(
    base_mva: float,
    bus_ids: list[str],
    bus_index: dict[str, int],
    n_bus: int,
    slack: int,
    bus_type: IntArray,
    branch_ids: list[str],
    branch_index: dict[str, int],
    n_branch: int,
    f: IntArray,
    t: IntArray,
    r: FloatArray,
    x: FloatArray,
    b: FloatArray,
    tap: FloatArray,
    shift_rad: FloatArray,
    rating_pu: FloatArray,
    p_load_pu: FloatArray,
    q_load_pu: FloatArray,
    g_shunt_pu: FloatArray,
    b_shunt_pu: FloatArray,
    p_gen_pu: FloatArray,
    q_gen_pu: FloatArray,
    p_min_pu: FloatArray,
    p_max_pu: FloatArray,
    q_min_pu: FloatArray,
    q_max_pu: FloatArray,
    v_set: FloatArray,
    gen_ids: list[str] = list(),
    gen_bus: IntArray = (lambda: zeros(0, dtype=int64))(),
    gen_p_pu: FloatArray = (lambda: zeros(0))(),
    gen_q_pu: FloatArray = (lambda: zeros(0))(),
    gen_p_min_pu: FloatArray = (lambda: zeros(0))(),
    gen_p_max_pu: FloatArray = (lambda: zeros(0))(),
    gen_q_min_pu: FloatArray = (lambda: zeros(0))(),
    gen_q_max_pu: FloatArray = (lambda: zeros(0))(),
    gen_v_set: FloatArray = (lambda: zeros(0))(),
    load_ids: list[str] = list(),
    load_bus: IntArray = (lambda: zeros(0, dtype=int64))(),
    load_p_min_pu: FloatArray = (lambda: zeros(0))(),
    load_p_max_pu: FloatArray = (lambda: zeros(0))(),
    storage_ids: list[str] = list(),
    storage_bus: IntArray = (
        lambda: zeros(0, dtype=int64)
    )(),
    storage_p_max_pu: FloatArray = (lambda: zeros(0))(),
    storage_energy_pu: FloatArray = (lambda: zeros(0))(),
    storage_soc_initial: FloatArray = (lambda: zeros(0))(),
    storage_efficiency_charge: FloatArray = (
        lambda: zeros(0)
    )(),
    storage_efficiency_discharge: FloatArray = (
        lambda: zeros(0)
    )(),
)

Frozen positional arrays over the in-service subset of a network, in per unit.

Positions are 0-based. bus_ids[i] is the id at bus position i; branch_ids[k] and gen_ids[g] likewise. Order follows the network's collection order with the excluded elements removed.

bus_type instance-attribute

bus_type: IntArray

1 = pq, 2 = pv, 3 = slack, as declared on the bus.

b instance-attribute

b: FloatArray

Total line charging susceptance per branch (pu); builders apply b / 2 per end.

tap instance-attribute

tap: FloatArray

Tap ratio magnitude on the from side; 1.0 where the branch has none.

rating_pu instance-attribute

rating_pu: FloatArray

Thermal rating in pu of base_mva; inf where the branch has none.

g_shunt_pu instance-attribute

g_shunt_pu: FloatArray

Shunt conductance per bus in pu (MATPOWER GS sign: positive consumes).

b_shunt_pu instance-attribute

b_shunt_pu: FloatArray

Shunt susceptance per bus in pu (MATPOWER BS sign: positive injects).

v_set instance-attribute

v_set: FloatArray

First in-service generator's v_set_pu at each bus; 1.0 where there is none.

load_p_max_pu class-attribute instance-attribute

load_p_max_pu: FloatArray = field(
    default_factory=lambda: np.zeros(0)
)

Per-load [0, p_mw] bound in pu: the natural bound for a bid-load's demand dispatch is zero up to its own fixed historical p_mw. Built uniformly for every in-service load regardless of Load.bid — Load carries no p_min_mw/p_max_mw fields to mirror the generator-side source data, and this bound formula does not depend on bid presence; whether/how a given load's bound is actually used by opf.dc_opf is decided per-load elsewhere, not here. p_load_pu/q_load_pu (the bus aggregate) are untouched by this addition.

storage_p_max_pu class-attribute instance-attribute

storage_p_max_pu: FloatArray = field(
    default_factory=lambda: np.zeros(0)
)

Charge/discharge power limit (Storage.p_max_mw), pu of base_mva.

storage_energy_pu class-attribute instance-attribute

storage_energy_pu: FloatArray = field(
    default_factory=lambda: np.zeros(0)
)

Energy capacity (Storage.energy_mwh), pu of base_mva (pu-hours) — the same base_mva-division convention every other physical field in this class already uses (ADR-005: physical units in the model, pu in numerics).

storage_soc_initial class-attribute instance-attribute

storage_soc_initial: FloatArray = field(
    default_factory=lambda: np.zeros(0)
)

Initial state of charge, already a fraction of energy_mwh in [0, 1] on the entity — dimensionless, so unlike storage_p_max_pu/storage_energy_pu it carries through unconverted.

storage_efficiency_discharge class-attribute instance-attribute

storage_efficiency_discharge: FloatArray = field(
    default_factory=lambda: np.zeros(0)
)

Charge/discharge efficiency, already dimensionless ratios in (0, 1] — unconverted.

from_network classmethod

from_network(net: Network) -> NetworkArrays

Build the in-service positional view; the one pu-conversion site.

Source code in src/mambo_power/numerics/arrays.py

@classmethod
def from_network(cls, net: Network) -> NetworkArrays:
    """Build the in-service positional view; the one pu-conversion site."""
    base = float(net.base_mva)

    live_buses = [bus for bus in net.buses if bus.in_service]
    bus_ids = [bus.id for bus in live_buses]
    bus_index = {bus_id: i for i, bus_id in enumerate(bus_ids)}
    n_bus = len(bus_ids)
    slack_positions = [i for i, bus in enumerate(live_buses) if bus.type == "slack"]
    if len(slack_positions) != 1:
        raise ValueError(
            f"expected exactly one in-service slack bus, found {len(slack_positions)}"
        )
    bus_type = np.fromiter(
        (BUS_TYPE_CODE[bus.type] for bus in live_buses), dtype=np.int64, count=n_bus
    )

    branches = [
        br
        for br in net.branches
        if br.in_service and br.from_bus in bus_index and br.to_bus in bus_index
    ]
    n_branch = len(branches)
    branch_ids = [br.id for br in branches]
    branch_index = {br_id: k for k, br_id in enumerate(branch_ids)}
    f = np.fromiter((bus_index[br.from_bus] for br in branches), np.int64, n_branch)
    t = np.fromiter((bus_index[br.to_bus] for br in branches), np.int64, n_branch)
    r = np.fromiter((br.r for br in branches), np.float64, n_branch)
    x = np.fromiter((br.x for br in branches), np.float64, n_branch)
    b = np.fromiter((br.b for br in branches), np.float64, n_branch)
    tap = np.fromiter(
        (1.0 if br.tap_ratio is None else br.tap_ratio for br in branches),
        np.float64,
        n_branch,
    )
    shift_rad = np.fromiter(
        (0.0 if br.shift_deg is None else math.radians(br.shift_deg) for br in branches),
        np.float64,
        n_branch,
    )
    rating_pu = np.fromiter(
        (math.inf if br.rating_mva is None else br.rating_mva / base for br in branches),
        np.float64,
        n_branch,
    )

    def per_bus(pairs: list[tuple[int, float]]) -> FloatArray:
        positions = np.fromiter((p for p, _ in pairs), np.int64, len(pairs))
        values = np.fromiter((v for _, v in pairs), np.float64, len(pairs))
        summed = np.asarray(
            np.bincount(positions, weights=values, minlength=n_bus), dtype=np.float64
        )
        return summed / base

    loads = [ld for ld in net.loads if ld.in_service and ld.bus in bus_index]
    p_load_pu = per_bus([(bus_index[ld.bus], ld.p_mw) for ld in loads])
    q_load_pu = per_bus([(bus_index[ld.bus], ld.q_mvar) for ld in loads])

    n_load = len(loads)
    load_bus = np.fromiter((bus_index[ld.bus] for ld in loads), np.int64, n_load)

    def per_load(values: list[float]) -> FloatArray:
        return np.fromiter(values, np.float64, n_load) / base

    shunts = [sh for sh in net.shunts if sh.in_service and sh.bus in bus_index]
    g_shunt_pu = per_bus([(bus_index[sh.bus], sh.g_mw) for sh in shunts])
    b_shunt_pu = per_bus([(bus_index[sh.bus], sh.b_mvar) for sh in shunts])

    gens = [g for g in net.generators if g.in_service and g.bus in bus_index]
    n_gen = len(gens)
    gen_bus = np.fromiter((bus_index[g.bus] for g in gens), np.int64, n_gen)
    p_gen_pu = per_bus([(bus_index[g.bus], g.p_mw) for g in gens])
    q_gen_pu = per_bus([(bus_index[g.bus], g.q_mvar) for g in gens])
    p_min_pu = per_bus([(bus_index[g.bus], g.p_min_mw) for g in gens])
    p_max_pu = per_bus([(bus_index[g.bus], g.p_max_mw) for g in gens])
    q_min_pu = per_bus([(bus_index[g.bus], g.q_min_mvar) for g in gens])
    q_max_pu = per_bus([(bus_index[g.bus], g.q_max_mvar) for g in gens])
    v_set = np.ones(n_bus)
    seen: set[int] = set()
    for g in gens:
        position = bus_index[g.bus]
        if position not in seen:
            seen.add(position)
            v_set[position] = g.v_set_pu

    def per_gen(values: list[float]) -> FloatArray:
        return np.fromiter(values, np.float64, n_gen) / base

    storage_units = [s for s in net.storage if s.in_service and s.bus in bus_index]
    n_storage = len(storage_units)
    storage_bus = np.fromiter((bus_index[s.bus] for s in storage_units), np.int64, n_storage)

    def per_storage(values: list[float]) -> FloatArray:
        return np.fromiter(values, np.float64, n_storage)

    return cls(
        base_mva=base,
        bus_ids=bus_ids,
        bus_index=bus_index,
        n_bus=n_bus,
        slack=slack_positions[0],
        bus_type=bus_type,
        branch_ids=branch_ids,
        branch_index=branch_index,
        n_branch=n_branch,
        f=f,
        t=t,
        r=r,
        x=x,
        b=b,
        tap=tap,
        shift_rad=shift_rad,
        rating_pu=rating_pu,
        p_load_pu=p_load_pu,
        q_load_pu=q_load_pu,
        g_shunt_pu=g_shunt_pu,
        b_shunt_pu=b_shunt_pu,
        p_gen_pu=p_gen_pu,
        q_gen_pu=q_gen_pu,
        p_min_pu=p_min_pu,
        p_max_pu=p_max_pu,
        q_min_pu=q_min_pu,
        q_max_pu=q_max_pu,
        v_set=v_set,
        gen_ids=[g.id for g in gens],
        gen_bus=gen_bus,
        gen_p_pu=per_gen([g.p_mw for g in gens]),
        gen_q_pu=per_gen([g.q_mvar for g in gens]),
        gen_p_min_pu=per_gen([g.p_min_mw for g in gens]),
        gen_p_max_pu=per_gen([g.p_max_mw for g in gens]),
        gen_q_min_pu=per_gen([g.q_min_mvar for g in gens]),
        gen_q_max_pu=per_gen([g.q_max_mvar for g in gens]),
        gen_v_set=np.fromiter((g.v_set_pu for g in gens), np.float64, n_gen),
        load_ids=[ld.id for ld in loads],
        load_bus=load_bus,
        load_p_min_pu=np.zeros(n_load),
        load_p_max_pu=per_load([ld.p_mw for ld in loads]),
        storage_ids=[s.id for s in storage_units],
        storage_bus=storage_bus,
        storage_p_max_pu=per_storage([s.p_max_mw / base for s in storage_units]),
        storage_energy_pu=per_storage([s.energy_mwh / base for s in storage_units]),
        storage_soc_initial=per_storage([s.soc_initial for s in storage_units]),
        storage_efficiency_charge=per_storage([s.efficiency_charge for s in storage_units]),
        storage_efficiency_discharge=per_storage(
            [s.efficiency_discharge for s in storage_units]
        ),
    )

Bus admittance matrix

mambo_power.numerics.ybus

Bus admittance matrix and branch admittance matrices (MATPOWER makeYbus conventions).

Per branch with series admittance y = 1 / (r + jx), total charging b and from-side complex tap a = tap · e^{j·shift}::

Yff = (y + j·b/2) / |a|²      Yft = -y / conj(a)
Ytf = -y / a                  Ytt =  y + j·b/2

Ybus = Cfᵀ·Yf + Ctᵀ·Yt + diag(g_shunt + j·b_shunt) with the shunt admittance already in pu (NetworkArrays divided by base_mva).

branch_admittances

branch_admittances(
    arr: NetworkArrays,
) -> tuple[
    ComplexArray, ComplexArray, ComplexArray, ComplexArray
]

Per-branch (Yff, Yft, Ytf, Ytt) vectors.

Source code in src/mambo_power/numerics/ybus.py

def branch_admittances(
    arr: NetworkArrays,
) -> tuple[ComplexArray, ComplexArray, ComplexArray, ComplexArray]:
    """Per-branch ``(Yff, Yft, Ytf, Ytt)`` vectors."""
    shorted = (arr.r == 0.0) & (arr.x == 0.0)
    if np.any(shorted):
        zero = [arr.branch_ids[k] for k in np.flatnonzero(shorted)]
        raise ValueError(
            f"series admittance undefined: r == x == 0 on in-service branch(es) {zero}"
        )
    ys = np.asarray(1.0 / (arr.r + 1j * arr.x), dtype=np.complex128)
    bc = np.asarray(1j * arr.b / 2.0, dtype=np.complex128)
    a = np.asarray(arr.tap * np.exp(1j * arr.shift_rad), dtype=np.complex128)
    yff: ComplexArray = (ys + bc) / (a * np.conj(a))
    yft: ComplexArray = -ys / np.conj(a)
    ytf: ComplexArray = -ys / a
    ytt: ComplexArray = ys + bc
    return yff, yft, ytf, ytt

yf_yt

yf_yt(arr: NetworkArrays) -> tuple[Any, Any]

(Yf, Yt): n_branch × n_bus complex CSC matrices giving from/to branch currents.

Yf @ V is the current injected into each branch at its from bus, Yt @ V at its to bus — the inputs M2 needs for branch flows.

Source code in src/mambo_power/numerics/ybus.py

def yf_yt(arr: NetworkArrays) -> tuple[Any, Any]:
    """``(Yf, Yt)``: ``n_branch × n_bus`` complex CSC matrices giving from/to branch currents.

    ``Yf @ V`` is the current injected into each branch at its from bus, ``Yt @ V`` at its to
    bus — the inputs M2 needs for branch flows.
    """
    yff, yft, ytf, ytt = branch_admittances(arr)
    rows = np.concatenate([np.arange(arr.n_branch), np.arange(arr.n_branch)])
    cols = np.concatenate([arr.f, arr.t])
    shape = (arr.n_branch, arr.n_bus)
    yf = sparse.csc_matrix((np.concatenate([yff, yft]), (rows, cols)), shape=shape)
    yt = sparse.csc_matrix((np.concatenate([ytf, ytt]), (rows, cols)), shape=shape)
    return yf, yt

ybus

ybus(arr: NetworkArrays) -> Any

The n_bus × n_bus complex CSC bus admittance matrix over the in-service subset.

Source code in src/mambo_power/numerics/ybus.py

def ybus(arr: NetworkArrays) -> Any:
    """The ``n_bus × n_bus`` complex CSC bus admittance matrix over the in-service subset."""
    yff, yft, ytf, ytt = branch_admittances(arr)
    rows = np.concatenate([arr.f, arr.f, arr.t, arr.t, np.arange(arr.n_bus)])
    cols = np.concatenate([arr.f, arr.t, arr.f, arr.t, np.arange(arr.n_bus)])
    data = np.concatenate([yff, yft, ytf, ytt, arr.g_shunt_pu + 1j * arr.b_shunt_pu])
    matrix = sparse.csc_matrix((data, (rows, cols)), shape=(arr.n_bus, arr.n_bus))
    matrix.sum_duplicates()
    return matrix

DC susceptance matrices

mambo_power.numerics.bbus

DC susceptance matrices and phase-shift injections (MATPOWER makeBdc conventions).

Per branch b = 1 / (x · tap) (tap magnitude only; r and line charging ignored). With Cft the n_branch × n_bus from-minus-to incidence matrix::

Bf   = diag(b) · Cft                 Pf   = Bf · θ + pf_shift
Bbus = Cftᵀ · Bf                     P    = Bbus · θ + p_shift

where pf_shift = -b · shift_rad and p_shift = Cftᵀ · pf_shift. A DC solve therefore reads Bbus · θ = P - p_shift.

branch_susceptance

branch_susceptance(arr: NetworkArrays) -> FloatArray

Per-branch DC susceptance 1 / (x · tap).

Source code in src/mambo_power/numerics/bbus.py

def branch_susceptance(arr: NetworkArrays) -> FloatArray:
    """Per-branch DC susceptance ``1 / (x · tap)``."""
    if np.any(arr.x == 0.0):
        zero = [arr.branch_ids[k] for k in np.flatnonzero(arr.x == 0.0)]
        raise UnsolvableNetworkError(
            f"DC susceptance undefined: x == 0 on in-service branch(es) {zero}"
        )
    result: FloatArray = 1.0 / (arr.x * arr.tap)
    return result

incidence

incidence(arr: NetworkArrays) -> Any

Cft: n_branch × n_bus sparse matrix with +1 at the from bus and -1 at the to bus.

Source code in src/mambo_power/numerics/bbus.py

def incidence(arr: NetworkArrays) -> Any:
    """``Cft``: ``n_branch × n_bus`` sparse matrix with +1 at the from bus and -1 at the to bus."""
    rows = np.concatenate([np.arange(arr.n_branch), np.arange(arr.n_branch)])
    cols = np.concatenate([arr.f, arr.t])
    data = np.concatenate([np.ones(arr.n_branch), -np.ones(arr.n_branch)])
    return sparse.csc_matrix((data, (rows, cols)), shape=(arr.n_branch, arr.n_bus))

bf

bf(arr: NetworkArrays) -> Any

Bf: n_branch × n_bus CSC matrix; Bf @ θ is the from-side DC flow (pu).

Source code in src/mambo_power/numerics/bbus.py

def bf(arr: NetworkArrays) -> Any:
    """``Bf``: ``n_branch × n_bus`` CSC matrix; ``Bf @ θ`` is the from-side DC flow (pu)."""
    b = branch_susceptance(arr)
    rows = np.concatenate([np.arange(arr.n_branch), np.arange(arr.n_branch)])
    cols = np.concatenate([arr.f, arr.t])
    data = np.concatenate([b, -b])
    return sparse.csc_matrix((data, (rows, cols)), shape=(arr.n_branch, arr.n_bus))

bbus

bbus(arr: NetworkArrays) -> Any

Bbus: n_bus × n_bus real CSC DC susceptance matrix (Cftᵀ · Bf).

Source code in src/mambo_power/numerics/bbus.py

def bbus(arr: NetworkArrays) -> Any:
    """``Bbus``: ``n_bus × n_bus`` real CSC DC susceptance matrix (``Cftᵀ · Bf``)."""
    matrix = (incidence(arr).T @ bf(arr)).tocsc()
    matrix.sum_duplicates()
    return matrix

pf_shift

pf_shift(arr: NetworkArrays) -> FloatArray

Per-branch phase-shifter flow injection -b · shift_rad (pu), at the from bus.

Source code in src/mambo_power/numerics/bbus.py

def pf_shift(arr: NetworkArrays) -> FloatArray:
    """Per-branch phase-shifter flow injection ``-b · shift_rad`` (pu), at the from bus."""
    result: FloatArray = branch_susceptance(arr) * (-arr.shift_rad)
    return result

p_shift

p_shift(arr: NetworkArrays) -> FloatArray

Per-bus phase-shifter injection Cftᵀ · pf_shift (pu); P = Bbus·θ + p_shift.

Source code in src/mambo_power/numerics/bbus.py

def p_shift(arr: NetworkArrays) -> FloatArray:
    """Per-bus phase-shifter injection ``Cftᵀ · pf_shift`` (pu); ``P = Bbus·θ + p_shift``."""
    result: FloatArray = np.asarray(incidence(arr).T @ pf_shift(arr), dtype=np.float64).ravel()
    return result

flow_from_ptdf

flow_from_ptdf(
    ptdf: FloatArray,
    injection_mw: FloatArray,
    arr: NetworkArrays,
) -> FloatArray

Branch flow, MW, from a PTDF matrix and a full bus net-injection vector, MW.

flow = ptdf @ (injection_mw − p_shift·base_mva) + pf_shift·base_mva — the phase-shifter injection is subtracted out of the bus injection before the PTDF product, then each branch's own from-side shift flow is added back on. This is exactly :func:mambo_power.pf.dc. solve's construction (its module docstring: rhs = P − p_shift feeds the angle solve, and p_from = Bf·θ + pf_shift; combined with θ = B'⁻¹(P − p_shift) and PTDF = Bf·B'⁻¹ on the reduced system, p_from = PTDF·(P − p_shift) + pf_shift) — the model every DC PTDF-based flow in this package must match. Omitting the − p_shift term (all of opf.dc_opf, opf.solve_dc_opf and market._clearing did until M8 finding F1 / A19) reproduces pf.solve_dc's flow only when no branch has a shift, since p_shift(arr) == 0 identically in that case.

injection_mw must be the full net injection per bus (generation minus load minus shunt, MW) — callers that instead fold some of that into a decision-variable-relative LP constant (opf.dc_opf's own flow-limit rows) derive the identical correction by hand rather than calling this helper, since their injection is not one vector (see that module's own derivation).

Source code in src/mambo_power/numerics/bbus.py

def flow_from_ptdf(ptdf: FloatArray, injection_mw: FloatArray, arr: NetworkArrays) -> FloatArray:
    """Branch flow, MW, from a PTDF matrix and a full bus net-injection vector, MW.

    ``flow = ptdf @ (injection_mw − p_shift·base_mva) + pf_shift·base_mva`` — the phase-shifter
    injection is subtracted out of the bus injection *before* the PTDF product, then each
    branch's own from-side shift flow is added back on. This is exactly :func:`mambo_power.pf.dc.
    solve`'s construction (its module docstring: ``rhs = P − p_shift`` feeds the angle solve, and
    ``p_from = Bf·θ + pf_shift``; combined with ``θ = B'⁻¹(P − p_shift)`` and ``PTDF = Bf·B'⁻¹``
    on the reduced system, ``p_from = PTDF·(P − p_shift) + pf_shift``) — the model every DC
    PTDF-based flow in this package must match. Omitting the ``− p_shift`` term (all of
    ``opf.dc_opf``, ``opf.solve_dc_opf`` and ``market._clearing`` did until M8 finding F1 / A19)
    reproduces ``pf.solve_dc``'s flow only when no branch has a shift, since
    ``p_shift(arr) == 0`` identically in that case.

    ``injection_mw`` must be the *full* net injection per bus (generation minus load minus
    shunt, MW) — callers that instead fold some of that into a decision-variable-relative LP
    constant (``opf.dc_opf``'s own flow-limit rows) derive the identical correction by hand
    rather than calling this helper, since their ``injection`` is not one vector (see that
    module's own derivation).
    """
    result: FloatArray = (
        ptdf @ (injection_mw - p_shift(arr) * arr.base_mva) + pf_shift(arr) * arr.base_mva
    )
    return result

PTDF

mambo_power.numerics.ptdf

Power transfer distribution factors from the DC model.

PTDF = Bf · Bbus⁻¹ with the slack row and column removed before the inverse and the slack column of the result set to zero: flows = PTDF @ P for any injection vector P (the slack absorbs the imbalance). The reduced Bbus is factorised once with a sparse LU and solved against the dense transposed Bf; the full matrix is never inverted densely.

ptdf

ptdf(
    arr: NetworkArrays, slack: int | None = None
) -> FloatArray

Dense n_branch × n_bus PTDF with a zero column at slack (default: the network's).

Source code in src/mambo_power/numerics/ptdf.py

def ptdf(arr: NetworkArrays, slack: int | None = None) -> FloatArray:
    """Dense ``n_branch × n_bus`` PTDF with a zero column at ``slack`` (default: the network's)."""
    ref = arr.slack if slack is None else slack
    if not 0 <= ref < arr.n_bus:
        raise ValueError(f"slack position {ref} out of range for {arr.n_bus} buses")
    keep = np.array([i for i in range(arr.n_bus) if i != ref], dtype=np.int64)
    result = np.zeros((arr.n_branch, arr.n_bus))
    if arr.n_branch == 0 or keep.size == 0:
        return result
    b_reduced = bbus(arr)[keep][:, keep].tocsc()
    bf_reduced = bf(arr)[:, keep].toarray()
    # Bbus is symmetric, so solving Bᵀ·X = Bfᵀ gives X = (Bf·B⁻¹)ᵀ.
    lu = splu(b_reduced)
    solved = lu.solve(np.ascontiguousarray(bf_reduced.T))
    result[:, keep] = solved.T
    return result

LODF and bridges

mambo_power.numerics.lodf

Line outage distribution factors and graph-theoretic bridge detection.

LODF[l, k] is the fraction of branch k's pre-outage flow that appears on branch l after k is removed. With h_k = PTDF·(e_f(k) − e_t(k)) the flows caused by a unit transfer across k::

LODF[l, k] = h_k[l] / (1 − h_k[k])     (l ≠ k)         LODF[k, k] = −1

A branch whose removal disconnects the network (a bridge) has h_k[k] = 1 and no finite LODF; its whole column is NaN. The numeric test (|1 − h_kk| < 1e-10) and the graph-theoretic bridges are independent and must agree — the test suite checks that.

BRIDGE_TOL module-attribute

BRIDGE_TOL = 1e-10

|1 − PTDF_kk| < BRIDGE_TOL marks branch k as a bridge in lodf.

lodf

lodf(
    arr: NetworkArrays,
    ptdf_matrix: FloatArray | None = None,
) -> FloatArray

Dense n_branch × n_branch LODF; bridge columns are NaN, diagonal is −1.

Source code in src/mambo_power/numerics/lodf.py

def lodf(arr: NetworkArrays, ptdf_matrix: FloatArray | None = None) -> FloatArray:
    """Dense ``n_branch × n_branch`` LODF; bridge columns are ``NaN``, diagonal is ``−1``."""
    h = ptdf(arr) if ptdf_matrix is None else ptdf_matrix
    if h.shape != (arr.n_branch, arr.n_bus):
        raise ValueError(f"ptdf_matrix has shape {h.shape}, expected {(arr.n_branch, arr.n_bus)}")
    # Column k: flows on every branch for a unit transfer from f[k] to t[k].
    transfer = h[:, arr.f] - h[:, arr.t]
    denominator = 1.0 - np.diag(transfer)
    is_bridge = np.abs(denominator) < BRIDGE_TOL
    safe = np.where(is_bridge, 1.0, denominator)
    result = transfer / safe[np.newaxis, :]
    np.fill_diagonal(result, -1.0)
    result[:, is_bridge] = np.nan
    return result

bridges

bridges(arr: NetworkArrays) -> list[int]

Positions of branches whose removal disconnects the in-service graph (sorted).

Iterative Tarjan lowpoint search over the multigraph; parallel branches between the same pair of buses are never bridges because the search skips only the edge it arrived by.

Source code in src/mambo_power/numerics/lodf.py

def bridges(arr: NetworkArrays) -> list[int]:
    """Positions of branches whose removal disconnects the in-service graph (sorted).

    Iterative Tarjan lowpoint search over the multigraph; parallel branches between the same
    pair of buses are never bridges because the search skips only the *edge* it arrived by.
    """
    n_bus, n_branch = arr.n_bus, arr.n_branch
    adjacency: list[list[tuple[int, int]]] = [[] for _ in range(n_bus)]
    for k in range(n_branch):
        u, v = int(arr.f[k]), int(arr.t[k])
        adjacency[u].append((v, k))
        adjacency[v].append((u, k))

    disc = [-1] * n_bus
    low = [0] * n_bus
    found: list[int] = []
    clock = 0
    for root in range(n_bus):
        if disc[root] != -1:
            continue
        disc[root] = low[root] = clock
        clock += 1
        # stack entries: (node, edge used to enter it, next adjacency cursor)
        stack: list[tuple[int, int, int]] = [(root, -1, 0)]
        while stack:
            node, via, cursor = stack[-1]
            if cursor < len(adjacency[node]):
                stack[-1] = (node, via, cursor + 1)
                nxt, edge = adjacency[node][cursor]
                if edge == via:
                    continue
                if disc[nxt] == -1:
                    disc[nxt] = low[nxt] = clock
                    clock += 1
                    stack.append((nxt, edge, 0))
                else:
                    low[node] = min(low[node], disc[nxt])
            else:
                stack.pop()
                if stack:
                    parent = stack[-1][0]
                    low[parent] = min(low[parent], low[node])
                    if low[node] > disc[parent]:
                        found.append(via)
    found.sort()
    return found