Skip to content

mambo_power.pf

Power-flow solvers. See the manual page for the formulation, the slack convention and the AC solver's contract.

mambo_power.pf

Power-flow solvers (epic Design §2 pf/): AC Newton-Raphson (W1) and DC (W2).

Public entry points take and return pydantic models (a Network in, a typed result out) and stamp provenance; the array-level solvers (mambo_power.pf.ac_newton.newton, mambo_power.pf.dc.solve) work on NetworkArrays only. Both entry points derive the bus roles through mambo_power.numerics.effective_roles (W3) — a PV bus without an in-service generator solves as PQ, a slack without one raises NoSlackGeneratorError, and a SetpointConflictWarning propagates to the caller.

AcOptions

Bases: BaseModel

Options of the AC Newton-Raphson solve (spec design item 1).

tol is compared against the infinity norm of the per-unit power mismatch (MATPOWER pf.tol semantics; pandapower's tolerance_mva is the same pu quantity despite its name). init="auto" warm-starts from the buses' stored vm_pu/va_deg when every in-service bus carries both, else flat; PV and slack magnitudes are always the setpoint.

tol class-attribute instance-attribute 

tol: float = 1e-08

Mismatch ∞-norm tolerance, pu.

max_iter class-attribute instance-attribute 

max_iter: int = 20

Newton iterations per solve.

q_limits class-attribute instance-attribute 

q_limits: bool = True

Enforce generator reactive limits.

max_q_rounds class-attribute instance-attribute 

max_q_rounds: int = 10

Maximum Q-limit re-solves.

init class-attribute instance-attribute 

init: Literal['auto', 'flat'] = 'auto'

Starting point rule.

AcSolution dataclass 

AcSolution(
    v: ComplexArray,
    converged: bool,
    iterations: int,
    max_mismatch_pu: float,
    q_limit_rounds: int,
    q_limited: IntArray,
    bus_type: IntArray,
    s_bus_pu: ComplexArray,
    gen_p_pu: FloatArray,
    gen_q_pu: FloatArray,
    message: str | None = None,
)

Positional AC solution in per unit, in NetworkArrays order.

v instance-attribute 

v: ComplexArray

Complex bus voltages (pu, radians inside).

converged instance-attribute 

converged: bool

Final Newton solve met tol and no Q-limit violation remained.

iterations instance-attribute 

iterations: int

Newton iterations summed over every Q-limit round.

max_mismatch_pu instance-attribute 

max_mismatch_pu: float

Infinity norm of the final mismatch vector, pu.

q_limit_rounds instance-attribute 

q_limit_rounds: int

Number of re-solves triggered by pinning (0 when nothing was pinned).

q_limited instance-attribute 

q_limited: IntArray

Per bus: 0 free, +1 pinned at ΣQmax, -1 pinned at ΣQmin.

bus_type instance-attribute 

bus_type: IntArray

Effective bus types after pinning (1 = pq, 2 = pv, 3 = slack).

s_bus_pu instance-attribute 

s_bus_pu: ComplexArray

Realised net complex injection per bus V·conj(Y V).

gen_p_pu instance-attribute 

gen_p_pu: FloatArray

Per-generator active output; the first slack-bus generator absorbs the balance.

gen_q_pu instance-attribute 

gen_q_pu: FloatArray

Per-generator reactive output, split by the MATPOWER pfsoln rule.

message class-attribute instance-attribute 

message: str | None = None

Diagnostic when converged is False; None otherwise.

DcSolution dataclass 

DcSolution(
    theta_rad: FloatArray,
    p_from_pu: FloatArray,
    p_inj_pu: FloatArray,
    gen_p_pu: FloatArray,
)

Positional DC solution in per unit, in NetworkArrays order.

theta_rad instance-attribute 

theta_rad: FloatArray

Bus angles, radians; theta_rad[arr.slack] == 0.

p_from_pu instance-attribute 

p_from_pu: FloatArray

From-side branch flow Bf·θ + pf_shift; the to-side flow is its negative.

p_inj_pu instance-attribute 

p_inj_pu: FloatArray

Realised net injection per bus B'·θ + p_shift (slack closes the balance).

gen_p_pu instance-attribute 

gen_p_pu: FloatArray

Per-generator output; the first in-service slack-bus generator absorbs the balance.

initial_voltage 

initial_voltage(
    net: Network,
    arr: NetworkArrays,
    roles: EffectiveRoles,
    options: AcOptions,
) -> ComplexArray

Starting voltages for solve_ac under options.init.

"flat": mambo_power.pf.ac_newton.flat_start. "auto": when every in-service bus carries both vm_pu and va_deg the stored state is the start (angles in radians, the slack keeping its stored angle), with PV and slack magnitudes replaced by the effective setpoint; otherwise flat.

Source code in src/mambo_power/pf/__init__.py 
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
def initial_voltage(
    net: Network, arr: NetworkArrays, roles: EffectiveRoles, options: AcOptions
) -> ComplexArray:
    """Starting voltages for :func:`solve_ac` under ``options.init``.

    ``"flat"``: :func:`mambo_power.pf.ac_newton.flat_start`. ``"auto"``: when every in-service
    bus carries both ``vm_pu`` and ``va_deg`` the stored state is the start (angles in radians,
    the slack keeping its stored angle), with PV and slack magnitudes replaced by the effective
    setpoint; otherwise flat.
    """
    if options.init == "auto":
        stored = {b.id: b for b in net.buses if b.in_service}
        if all(b.vm_pu is not None and b.va_deg is not None for b in stored.values()):
            vm = np.array([float(stored[i].vm_pu or 0.0) for i in arr.bus_ids])
            va = np.array([math.radians(float(stored[i].va_deg or 0.0)) for i in arr.bus_ids])
            held = roles.bus_type != BUS_TYPE_CODE["pq"]
            vm[held] = roles.v_set[held]
            return np.asarray(vm * np.exp(1j * va), dtype=np.complex128)
    return ac_newton.flat_start(arr, roles)

solve_ac 

solve_ac(
    net: Network, *, options: AcOptions | None = None
) -> AcPowerFlowResult

AC power flow of net by Newton-Raphson (mambo_power.pf.ac_newton).

Builds the in-service NetworkArrays and the effective roles, solves with Q-limit enforcement per options, computes branch flows S_from = V_f · conj(Yf V) and S_to = V_t · conj(Yt V), and returns an AcPowerFlowResult in MW/MVAr keyed by ids with provenance (kind = "pf.ac", solver = scipy.sparse.linalg.splu, the options as run). A solve that does not converge is reported through converged = False — never raised. The network is not modified.

Source code in src/mambo_power/pf/__init__.py 
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
def solve_ac(net: Network, *, options: AcOptions | None = None) -> AcPowerFlowResult:
    """AC power flow of ``net`` by Newton-Raphson (:mod:`mambo_power.pf.ac_newton`).

    Builds the in-service :class:`NetworkArrays` and the effective roles, solves with Q-limit
    enforcement per ``options``, computes branch flows ``S_from = V_f · conj(Yf V)`` and
    ``S_to = V_t · conj(Yt V)``, and returns an :class:`~mambo_power.results.AcPowerFlowResult`
    in MW/MVAr keyed by ids with provenance (``kind = "pf.ac"``, ``solver =
    scipy.sparse.linalg.splu``, the options as run). A solve that does not converge is
    reported through ``converged = False`` — never raised. The network is not modified.
    """
    opts = options if options is not None else AcOptions()
    started_at = datetime.now(UTC)
    clock = time.perf_counter()
    arr = NetworkArrays.from_network(net)
    roles = effective_roles(arr)
    v0 = initial_voltage(net, arr, roles, opts)
    sol = ac_newton.newton(arr, roles, opts, v0=v0)
    yf, yt = yf_yt(arr)
    s_from = np.asarray(sol.v[arr.f] * np.conj(yf @ sol.v), dtype=np.complex128)
    s_to = np.asarray(sol.v[arr.t] * np.conj(yt @ sol.v), dtype=np.complex128)
    # reported bus injection = generation − load − shunt: take the shunt's |V|²·conj(g + jb)
    # out of V·conj(Y V) (Y carries the shunts), matching the DC result and pandapower's res_bus
    vm2 = np.abs(sol.v) ** 2
    s_inj = np.asarray(
        sol.s_bus_pu - vm2 * (arr.g_shunt_pu - 1j * arr.b_shunt_pu), dtype=np.complex128
    )
    elapsed_s = time.perf_counter() - clock
    provenance = ResultProvenance(
        engine="mambo-power",
        version=mambo_power.__version__,
        kind="pf.ac",
        solver=ac_newton.SOLVER,
        started_at=started_at,
        elapsed_s=elapsed_s,
        options=opts.model_dump(),
    )
    return ac_result_from_arrays(
        arr,
        v=sol.v,
        s_bus_pu=s_inj,
        s_from_pu=s_from,
        s_to_pu=s_to,
        gen_p_pu=sol.gen_p_pu,
        gen_q_pu=sol.gen_q_pu,
        bus_type=sol.bus_type,
        q_limited=sol.q_limited,
        converged=sol.converged,
        iterations=sol.iterations,
        max_mismatch_pu=sol.max_mismatch_pu,
        q_limit_rounds=sol.q_limit_rounds,
        provenance=provenance,
        message=sol.message,
    )

solve_dc 

solve_dc(net: Network) -> DcPowerFlowResult

DC power flow of net: lossless B'θ = P with phase shifts, flows via Bf.

Builds the in-service NetworkArrays, runs mambo_power.pf.dc.solve, and returns a DcPowerFlowResult in MW keyed by ids, with provenance (version = mambo_power.__version__, solver = scipy.sparse.linalg.splu, UTC start time, wall-clock duration). role_effective reports the effective roles (W3): the solve itself needs no setpoints, but a slack without an in-service generator is still an error and a gen-less PV bus is reported as PQ. The network is not modified.

Source code in src/mambo_power/pf/__init__.py 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
def solve_dc(net: Network) -> DcPowerFlowResult:
    """DC power flow of ``net``: lossless ``B'θ = P`` with phase shifts, flows via ``Bf``.

    Builds the in-service :class:`NetworkArrays`, runs :func:`mambo_power.pf.dc.solve`, and
    returns a :class:`~mambo_power.results.DcPowerFlowResult` in MW keyed by ids, with
    provenance (``version = mambo_power.__version__``, ``solver = scipy.sparse.linalg.splu``,
    UTC start time, wall-clock duration). ``role_effective`` reports the effective roles (W3):
    the solve itself needs no setpoints, but a slack without an in-service generator is still
    an error and a gen-less PV bus is reported as PQ. The network is not modified.
    """
    started_at = datetime.now(UTC)
    clock = time.perf_counter()
    arr = NetworkArrays.from_network(net)
    roles = effective_roles(arr)
    sol = dc.solve(arr)
    elapsed_s = time.perf_counter() - clock
    provenance = ResultProvenance(
        engine="mambo-power",
        version=mambo_power.__version__,
        kind="pf.dc",
        solver=dc.SOLVER,
        started_at=started_at,
        elapsed_s=elapsed_s,
        options={},
    )
    return dc_result_from_arrays(
        arr,
        theta_rad=sol.theta_rad,
        p_from_pu=sol.p_from_pu,
        p_inj_pu=sol.p_inj_pu,
        gen_p_pu=sol.gen_p_pu,
        provenance=provenance,
        bus_type=roles.bus_type,
    )

DC solver over arrays

mambo_power.pf.dc

DC power flow over NetworkArrays (MATPOWER rundcpf).

Formulation. With the DC susceptance matrix B' (mambo_power.numerics.bbus), the from-side flow matrix Bf (mambo_power.numerics.bf), the phase-shifter injections p_shift = Cftᵀ · pf_shift and pf_shift = -b · shift (mambo_power.numerics.bbus), and the declared net injection per bus in per unit::

P_bus = P_gen − P_load − G_shunt          (G_shunt: conductance consumption at 1.0 pu)

the angles solve the linear system with the slack row and column removed and θ_slack = 0::

B'[keep, keep] · θ[keep] = (P_bus − p_shift)[keep]          θ[slack] = 0

and the flows and realised injections follow::

p_from = Bf · θ + pf_shift          p_to = −p_from          p_inj = B' · θ + p_shift

p_inj equals P_bus on every non-slack bus; at the slack it is whatever closes the balance (lossless, so Σ p_inj = 0). Those are exactly MATPOWER's rundcpf steps — Pbus = real(makeSbus) − Pbusinj − GS/baseMVA, Va = dcpf(B, Pbus, Va0, ref, pv, pq), PF = (Bf·Va + Pfinj)·baseMVA, PT = −PF — which pandapower's rundcpp copies verbatim (record/m2-research.md §2).

Slack generation. The slack-bus balance goes entirely to the first in-service generator at the slack bus, every other generator keeping its dispatch (MATPOWER rundcpf: gen(on(refgen(1)), PG) += (B(ref,:)·Va − Pbus(ref))·baseMVA; pandapower reports the same number on res_ext_grid). Bus-level generation is therefore engine-independent; the per-generator split is a documented convention. If the slack bus carries no in-service generator the balance is still visible on the bus injection; naming that situation is W3's effective_roles.

The reduced system is factorised with scipy.sparse.linalg.splu (the same backend as the PTDF builder and the AC Newton solve). Bus roles are the declared roles from the arrays — DC needs no generator setpoints, so effective-role derivation (W3) does not enter here.

SOLVER module-attribute 

SOLVER = 'scipy.sparse.linalg.splu'

Linear-algebra backend name stamped into the result provenance.

DcSolution dataclass 

DcSolution(
    theta_rad: FloatArray,
    p_from_pu: FloatArray,
    p_inj_pu: FloatArray,
    gen_p_pu: FloatArray,
)

Positional DC solution in per unit, in NetworkArrays order.

theta_rad instance-attribute 

theta_rad: FloatArray

Bus angles, radians; theta_rad[arr.slack] == 0.

p_from_pu instance-attribute 

p_from_pu: FloatArray

From-side branch flow Bf·θ + pf_shift; the to-side flow is its negative.

p_inj_pu instance-attribute 

p_inj_pu: FloatArray

Realised net injection per bus B'·θ + p_shift (slack closes the balance).

gen_p_pu instance-attribute 

gen_p_pu: FloatArray

Per-generator output; the first in-service slack-bus generator absorbs the balance.

declared_injection 

declared_injection(arr: NetworkArrays) -> FloatArray

P_gen − P_load − G_shunt per bus in pu — the right-hand side before phase shifts.

Source code in src/mambo_power/pf/dc.py 
69
70
71
72
def declared_injection(arr: NetworkArrays) -> FloatArray:
    """``P_gen − P_load − G_shunt`` per bus in pu — the right-hand side before phase shifts."""
    result: FloatArray = arr.p_gen_pu - arr.p_load_pu - arr.g_shunt_pu
    return result

solve 

solve(arr: NetworkArrays) -> DcSolution

Solve B'θ = P − p_shift with the slack at 0 and return angles, flows and injections.

Raises UnsolvableNetworkError when a branch has x == 0 (susceptance undefined; user data DC cannot solve, distinct from a malformed-input ValueError) or ValueError when the reduced B' is singular / yields non-finite angles (an islanded bus set).

Source code in src/mambo_power/pf/dc.py 
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
def solve(arr: NetworkArrays) -> DcSolution:
    """Solve ``B'θ = P − p_shift`` with the slack at 0 and return angles, flows and injections.

    Raises :class:`~mambo_power.numerics.UnsolvableNetworkError` when a branch has ``x == 0``
    (susceptance undefined; user data DC cannot solve, distinct from a malformed-input
    ``ValueError``) or ``ValueError`` when the reduced ``B'`` is singular / yields non-finite
    angles (an islanded bus set).
    """
    b_matrix = bbus(arr)
    p_declared = declared_injection(arr)
    rhs = p_declared - p_shift(arr)

    theta = np.zeros(arr.n_bus)
    keep = np.array([i for i in range(arr.n_bus) if i != arr.slack], dtype=np.int64)
    if keep.size:
        reduced = b_matrix[keep][:, keep].tocsc()
        try:
            theta[keep] = splu(reduced).solve(rhs[keep])
        except RuntimeError as exc:  # SuperLU: "Factor is exactly singular"
            raise ValueError(f"DC power flow: reduced B' is singular ({exc})") from exc
    if not np.all(np.isfinite(theta)):
        raise ValueError("DC power flow: non-finite angles (reduced B' is singular)")

    p_from: FloatArray = np.asarray(bf(arr) @ theta, dtype=np.float64).ravel() + pf_shift(arr)
    p_inj: FloatArray = np.asarray(b_matrix @ theta, dtype=np.float64).ravel() + p_shift(arr)

    # realised gross generation the slack bus supplies = realised net injection plus what
    # declared_injection subtracted to declare it (load, shunt); absorb_slack_p undoes that
    # subtraction on the other side (arr.p_gen_pu[arr.slack]) — see its docstring.
    p_bus = float(p_inj[arr.slack] + arr.p_load_pu[arr.slack] + arr.g_shunt_pu[arr.slack])
    gen_p = absorb_slack_p(arr, p_bus)

    return DcSolution(theta_rad=theta, p_from_pu=p_from, p_inj_pu=p_inj, gen_p_pu=gen_p)

AC solver over arrays

mambo_power.pf.ac_newton

AC power flow by polar Newton-Raphson over NetworkArrays (W1).

Formulation (MATPOWER newtonpf). With the bus admittance matrix Y (mambo_power.numerics.ybus), the complex voltages V = Vm·e^{jVa} and the specified net injections S_spec = (P_gen − P_load) + j(Q_gen − Q_load) in per unit (shunts live in Y), the mismatch is::

ΔS = V · conj(Y V) − S_spec
F  = [ real(ΔS)[pv ∪ pq] ; imag(ΔS)[pq] ]

The state is x = [Va[pv ∪ pq]; Vm[pq]] and each iteration solves J·Δx = −F with the Jacobian assembled from the sparse partial derivatives (MATPOWER dSbus_dV, polar)::

∂S/∂Vm = diag(V) · conj(Y · diag(V/|V|)) + conj(diag(Y V)) · diag(V/|V|)
∂S/∂Va = j · diag(V) · conj(diag(Y V) − Y · diag(V))

J = [ real(∂S/∂Va)[pvpq, pvpq]   real(∂S/∂Vm)[pvpq, pq] ]
    [ imag(∂S/∂Va)[pq,   pvpq]   imag(∂S/∂Vm)[pq,   pq] ]

factorised with scipy.sparse.linalg.splu. The mismatch is tested before each step, so a start that already satisfies ‖F‖∞ ≤ tol reports zero iterations; the loop stops with converged = False after max_iter updates, on a singular Jacobian, or when an update produces a non-finite voltage (the last finite iterate is returned).

Start. Flat: Vm = 1, Va = 0 at PQ buses, Vm = v_set (the effective setpoint from mambo_power.numerics.effective_roles) and Va = 0 at PV and slack buses. Warm: a caller-supplied v0 (solve_ac builds one from the buses' stored vm_pu/va_deg under init="auto"); PV and slack magnitudes are always overridden by the setpoint. The slack angle is whatever the start carries (0 for flat).

Q-limit enforcement (pandapower semantics, spec design item 3; pandapower 3.3.0 pf/run_newton_raphson_pf.py:182-249 _run_ac_pf_with_qlims_enforced, itself MATPOWER runpf.m:366-440 with pf.enforce_q_lims = 1). After every converged Newton solve the reactive generation per bus is Qg = imag(V·conj(Y V)) + Q_load; every bus still PV whose Qg > ΣQmax or Qg < ΣQmin (aggregate over its in-service generators, strict comparison — pandapower :199-200; MATPOWER adds a 5e-6 opf.violation slack, pandapower does not) is converted to PQ with Q_spec = Q_limit − Q_load (:224-242: the generator's QG is pinned at the limit and folded into the bus load). All violators of a round are converted together (enforce_q_lims=1, simultaneous), the next solve warm-starts from the current voltages, and pins accumulate — a pinned bus is never restored to PV (limited = r_[limited, mx], :235; the spec rejects the restore). The slack bus is never converted (setdiff1d(changed_gens, ref), :227). The loop ends when a converged solve shows no new violation; if violations persist after max_q_rounds re-solves the result carries converged = False and a diagnostic (pandapower would raise LoadflowNotConverged). A Newton solve that fails to converge ends the loop immediately without pinning.

Generator allocation (MATPOWER pfsoln; pandapower pypower/pfsoln.py:109-141 is a verbatim copy). Active power: every generator keeps its dispatch except the first in-service generator at the slack bus, which absorbs the slack-bus balance (the rule pf.dc already applies). Reactive power: the bus total Qg_bus = imag(S) + Q_load is split among the bus's in-service generators — equally when every generator's range is zero, otherwise Qg_i = Qmin_i + (Qg_bus − ΣQmin) / (ΣQmax − ΣQmin) · (Qmax_i − Qmin_i) (proportional to each generator's reactive range). A pinned bus's generators therefore sit exactly at their individual limits.

SOLVER module-attribute 

SOLVER = 'scipy.sparse.linalg.splu'

Linear-algebra backend name stamped into the result provenance.

AcOptions

Bases: BaseModel

Options of the AC Newton-Raphson solve (spec design item 1).

tol is compared against the infinity norm of the per-unit power mismatch (MATPOWER pf.tol semantics; pandapower's tolerance_mva is the same pu quantity despite its name). init="auto" warm-starts from the buses' stored vm_pu/va_deg when every in-service bus carries both, else flat; PV and slack magnitudes are always the setpoint.

tol class-attribute instance-attribute 

tol: float = 1e-08

Mismatch ∞-norm tolerance, pu.

max_iter class-attribute instance-attribute 

max_iter: int = 20

Newton iterations per solve.

q_limits class-attribute instance-attribute 

q_limits: bool = True

Enforce generator reactive limits.

max_q_rounds class-attribute instance-attribute 

max_q_rounds: int = 10

Maximum Q-limit re-solves.

init class-attribute instance-attribute 

init: Literal['auto', 'flat'] = 'auto'

Starting point rule.

AcSolution dataclass 

AcSolution(
    v: ComplexArray,
    converged: bool,
    iterations: int,
    max_mismatch_pu: float,
    q_limit_rounds: int,
    q_limited: IntArray,
    bus_type: IntArray,
    s_bus_pu: ComplexArray,
    gen_p_pu: FloatArray,
    gen_q_pu: FloatArray,
    message: str | None = None,
)

Positional AC solution in per unit, in NetworkArrays order.

v instance-attribute 

v: ComplexArray

Complex bus voltages (pu, radians inside).

converged instance-attribute 

converged: bool

Final Newton solve met tol and no Q-limit violation remained.

iterations instance-attribute 

iterations: int

Newton iterations summed over every Q-limit round.

max_mismatch_pu instance-attribute 

max_mismatch_pu: float

Infinity norm of the final mismatch vector, pu.

q_limit_rounds instance-attribute 

q_limit_rounds: int

Number of re-solves triggered by pinning (0 when nothing was pinned).

q_limited instance-attribute 

q_limited: IntArray

Per bus: 0 free, +1 pinned at ΣQmax, -1 pinned at ΣQmin.

bus_type instance-attribute 

bus_type: IntArray

Effective bus types after pinning (1 = pq, 2 = pv, 3 = slack).

s_bus_pu instance-attribute 

s_bus_pu: ComplexArray

Realised net complex injection per bus V·conj(Y V).

gen_p_pu instance-attribute 

gen_p_pu: FloatArray

Per-generator active output; the first slack-bus generator absorbs the balance.

gen_q_pu instance-attribute 

gen_q_pu: FloatArray

Per-generator reactive output, split by the MATPOWER pfsoln rule.

message class-attribute instance-attribute 

message: str | None = None

Diagnostic when converged is False; None otherwise.

flat_start 

flat_start(
    arr: NetworkArrays, roles: EffectiveRoles
) -> ComplexArray

1∠0 at PQ buses, v_set∠0 at PV and slack buses (effective roles).

Source code in src/mambo_power/pf/ac_newton.py 
137
138
139
140
141
142
def flat_start(arr: NetworkArrays, roles: EffectiveRoles) -> ComplexArray:
    """``1∠0`` at PQ buses, ``v_set∠0`` at PV and slack buses (effective roles)."""
    vm = np.ones(arr.n_bus)
    held = roles.bus_type != _PQ
    vm[held] = roles.v_set[held]
    return np.asarray(vm, dtype=np.complex128)

specified_injection 

specified_injection(arr: NetworkArrays) -> ComplexArray

(P_gen − P_load) + j(Q_gen − Q_load) per bus, pu; shunts are in Y, not here.

Source code in src/mambo_power/pf/ac_newton.py 
145
146
147
148
149
150
def specified_injection(arr: NetworkArrays) -> ComplexArray:
    """``(P_gen − P_load) + j(Q_gen − Q_load)`` per bus, pu; shunts are in ``Y``, not here."""
    return np.asarray(
        (arr.p_gen_pu - arr.p_load_pu) + 1j * (arr.q_gen_pu - arr.q_load_pu),
        dtype=np.complex128,
    )

newton_raphson 

newton_raphson(
    y: Any,
    s_spec: ComplexArray,
    v0: ComplexArray,
    pv: IntArray,
    pq: IntArray,
    *,
    tol: float,
    max_iter: int,
) -> tuple[ComplexArray, bool, int, float, str | None]

One Newton solve (MATPOWER newtonpf): (V, converged, iterations, ‖F‖∞, message).

Source code in src/mambo_power/pf/ac_newton.py 
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
def newton_raphson(
    y: Any,
    s_spec: ComplexArray,
    v0: ComplexArray,
    pv: IntArray,
    pq: IntArray,
    *,
    tol: float,
    max_iter: int,
) -> tuple[ComplexArray, bool, int, float, str | None]:
    """One Newton solve (MATPOWER ``newtonpf``): ``(V, converged, iterations, ‖F‖∞, message)``."""
    v = v0.astype(np.complex128, copy=True)
    va = np.angle(v)
    vm = np.abs(v)
    pvpq = np.concatenate([pv, pq])
    n_a = pvpq.size
    n_m = pq.size

    def mismatch(v: ComplexArray) -> FloatArray:
        ds = v * np.conj(y @ v) - s_spec
        return np.asarray(np.concatenate([ds.real[pvpq], ds.imag[pq]]), dtype=np.float64)

    f = mismatch(v)
    norm = float(np.max(np.abs(f))) if f.size else 0.0
    norm0 = norm
    converged = bool(np.isfinite(norm) and norm <= tol)
    iterations = 0
    message: str | None = None
    while not converged and iterations < max_iter:
        iterations += 1
        ds_dvm, ds_dva = _dsbus_dv(y, v)
        j11 = ds_dva[pvpq, :][:, pvpq].real
        j12 = ds_dvm[pvpq, :][:, pq].real
        j21 = ds_dva[pq, :][:, pvpq].imag
        j22 = ds_dvm[pq, :][:, pq].imag
        jac = sparse.bmat([[j11, j12], [j21, j22]], format="csc")
        try:
            dx = -splu(jac).solve(f)
        except RuntimeError as exc:  # SuperLU: "Factor is exactly singular"
            message = f"singular Jacobian at iteration {iterations} ({exc})"
            break
        va_new = va.copy()
        vm_new = vm.copy()
        va_new[pvpq] += dx[:n_a]
        vm_new[pq] += dx[n_a : n_a + n_m]
        v_new = np.asarray(vm_new * np.exp(1j * va_new), dtype=np.complex128)
        if not np.all(np.isfinite(v_new)):
            message = f"non-finite voltage at iteration {iterations}"
            break
        v = v_new
        vm = np.abs(v)  # MATPOWER: re-normalise so a negative Vm step flips the angle
        va = np.angle(v)
        f = mismatch(v)
        norm = float(np.max(np.abs(f))) if f.size else 0.0
        if not np.isfinite(norm):
            message = f"non-finite mismatch at iteration {iterations}"
            break
        if norm0 > 0.0 and norm > _DIVERGENCE_FACTOR * norm0:
            message = (
                f"diverging: ‖F‖∞ = {norm:.3e} pu exceeds "
                f"{_DIVERGENCE_FACTOR:.0e}× the starting mismatch ({norm0:.3e} pu) "
                f"at iteration {iterations}"
            )
            break
        converged = norm <= tol
    if not converged and message is None:
        message = f"did not converge in {max_iter} iterations (‖F‖∞ = {norm:.3e} pu)"
    return v, converged, iterations, norm, message

allocate_generation 

allocate_generation(
    arr: NetworkArrays,
    s_bus: ComplexArray,
    q_limited: IntArray,
) -> tuple[FloatArray, FloatArray]

Per-generator (P, Q) from the bus totals: MATPOWER pfsoln rules (module docstring).

A pinned bus (q_limited ±1) reports exactly its aggregate limit rather than the solved imag(S) + Q_load, which differs from it by the convergence tolerance (pandapower restores fixedQg the same way, run_newton_raphson_pf.py:246).

Source code in src/mambo_power/pf/ac_newton.py 
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
def allocate_generation(
    arr: NetworkArrays, s_bus: ComplexArray, q_limited: IntArray
) -> tuple[FloatArray, FloatArray]:
    """Per-generator ``(P, Q)`` from the bus totals: MATPOWER ``pfsoln`` rules (module docstring).

    A pinned bus (``q_limited`` ±1) reports exactly its aggregate limit rather than the solved
    ``imag(S) + Q_load``, which differs from it by the convergence tolerance (pandapower restores
    ``fixedQg`` the same way, ``run_newton_raphson_pf.py:246``).
    """
    n_gen = len(arr.gen_ids)
    gen_q = np.zeros(n_gen)
    if n_gen == 0:
        return arr.gen_p_pu.copy(), gen_q
    p_bus = float(s_bus[arr.slack].real + arr.p_load_pu[arr.slack])
    gen_p = absorb_slack_p(arr, p_bus)
    q_bus = s_bus.imag + arr.q_load_pu
    q_bus[q_limited == 1] = arr.q_max_pu[q_limited == 1]
    q_bus[q_limited == -1] = arr.q_min_pu[q_limited == -1]
    counts = np.bincount(arr.gen_bus, minlength=arr.n_bus)
    for position in np.flatnonzero(counts):
        rows = np.flatnonzero(arr.gen_bus == position)
        q_min = arr.gen_q_min_pu[rows]
        q_max = arr.gen_q_max_pu[rows]
        total_range = float(np.sum(q_max - q_min))
        if total_range > 0.0:
            share = (q_bus[position] - float(np.sum(q_min))) / total_range
            gen_q[rows] = q_min + share * (q_max - q_min)
        else:
            gen_q[rows] = q_bus[position] / rows.size
    return gen_p, gen_q

newton 

newton(
    arr: NetworkArrays,
    roles: EffectiveRoles,
    opts: AcOptions,
    v0: ComplexArray | None = None,
) -> AcSolution

Solve the AC power flow of arr with the effective roles (module docstring).

v0 is the starting voltage (flat start when None); PV and slack magnitudes in it are replaced by the effective setpoints. Never raises for a non-converged solve — the result carries converged = False and a message.

Source code in src/mambo_power/pf/ac_newton.py 
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
def newton(
    arr: NetworkArrays,
    roles: EffectiveRoles,
    opts: AcOptions,
    v0: ComplexArray | None = None,
) -> AcSolution:
    """Solve the AC power flow of ``arr`` with the effective ``roles`` (module docstring).

    ``v0`` is the starting voltage (flat start when ``None``); PV and slack magnitudes in it are
    replaced by the effective setpoints. Never raises for a non-converged solve — the result
    carries ``converged = False`` and a ``message``.
    """
    y = ybus(arr)
    s_spec = specified_injection(arr)
    bus_type = roles.bus_type.copy()
    q_limited = np.zeros(arr.n_bus, dtype=np.int64)
    if v0 is None:
        v = flat_start(arr, roles)
    else:
        v = np.asarray(v0, dtype=np.complex128).copy()
        held = bus_type != _PQ
        v[held] = roles.v_set[held] * np.exp(1j * np.angle(v[held]))

    total_iterations = 0
    rounds = 0
    message: str | None = None
    while True:
        pv = np.flatnonzero(bus_type == _PV).astype(np.int64)
        pq = np.flatnonzero(bus_type == _PQ).astype(np.int64)
        v, converged, iterations, norm, message = newton_raphson(
            y, s_spec, v, pv, pq, tol=opts.tol, max_iter=opts.max_iter
        )
        total_iterations += iterations
        if not converged or not opts.q_limits:
            break
        s_calc = v * np.conj(y @ v)
        q_gen = s_calc.imag + arr.q_load_pu
        over = pv[q_gen[pv] > arr.q_max_pu[pv]]
        under = pv[q_gen[pv] < arr.q_min_pu[pv]]
        if over.size == 0 and under.size == 0:
            break
        if rounds >= opts.max_q_rounds:
            converged = False
            violators = [arr.bus_ids[i] for i in np.concatenate([over, under])]
            message = (
                f"Q-limit enforcement did not settle within max_q_rounds={opts.max_q_rounds}; "
                f"still violating: {violators}"
            )
            break
        rounds += 1
        bus_type[over] = _PQ
        bus_type[under] = _PQ
        q_limited[over] = 1
        q_limited[under] = -1
        s_spec[over] = s_spec[over].real + 1j * (arr.q_max_pu[over] - arr.q_load_pu[over])
        s_spec[under] = s_spec[under].real + 1j * (arr.q_min_pu[under] - arr.q_load_pu[under])

    s_bus = np.asarray(v * np.conj(y @ v), dtype=np.complex128)
    gen_p, gen_q = allocate_generation(arr, s_bus, q_limited)
    return AcSolution(
        v=v,
        converged=converged,
        iterations=total_iterations,
        max_mismatch_pu=norm,
        q_limit_rounds=rounds,
        q_limited=q_limited,
        bus_type=bus_type,
        s_bus_pu=s_bus,
        gen_p_pu=gen_p,
        gen_q_pu=gen_q,
        message=None if converged else message,
    )