An abstract supertype for AC power flow evaluation models, parametrized by the solver type S <: ACPowerFlowSolverType. Concrete subtypes select the formulation: ACPolarPowerFlow uses the polar voltage state; a rectangular current-injection formulation is provided separately. The solver and the formulation are orthogonal.

A type alias for a PowerFlowData struct whose type parameters are configured for the DCPowerFlow method.

A type alias for a PowerFlowData struct whose type parameters are configured for an AC power flow method (ACPolarPowerFlow or ACRectangularPowerFlow, i.e. any AbstractACPowerFlow).

A type alias for a PowerFlowData struct whose type parameters are configured for the PTDFDCPowerFlow method .

PowerFlowData{M <: PNM.PowerNetworkMatrix, N <: Union{PNM.PowerNetworkMatrix, Nothing}}

Structure containing all the data required for the evaluation of the power flows and angles, as well as these ones.

All fields starting with bus_ are ordered according to bus_lookup, and all fields starting with arc_ are ordered according to arc_lookup: one row per bus/arc, one column per time period. Here, buses should be understood as "buses remaining, after the network reduction." Similarly, we use "arcs" instead of "branches" to distinguish between network elements (post-reduction) and system objects (pre-reduction).

Generally, do not construct this directly. Instead, use one of the later constructors to pass in a PowerFlowEvaluationModel and a PowerSystems.System. aux\_network\_matrix and power\_network\_matrix will then be set to the appropriate matrices that are needed for computing that type of power flow. See also ACPowerFlowData, ABAPowerFlowData, PTDFPowerFlowData, and vPTDFPowerFlowData: these are all aliases for PowerFlowData{N, M} with specific N,M, that are used for the respective type of power flow evaluations.

Fields:

• bus_active_power_injections::Matrix{Float64}: matrix containing the bus active power injections.
• bus_reactive_power_injections::Matrix{Float64}: matrix containing the bus reactive power injections.
• bus_active_power_withdrawals::Matrix{Float64}: matrix containing the bus reactive power withdrawals.
• bus_reactive_power_withdrawals::Matrix{Float64}: matrix containing the bus reactive power withdrawals.
• bus_active_power_constant_current_withdrawals::Matrix{Float64}: matrix containing the bus active power constant current withdrawals.
• bus_reactive_power_constant_current_withdrawals::Matrix{Float64}: matrix containing the bus reactive power constant current withdrawals.
• bus_active_power_constant_impedance_withdrawals::Matrix{Float64}: matrix containing the bus active power constant impedance withdrawals.
• bus_reactive_power_constant_impedance_withdrawals::Matrix{Float64}: matrix containing the bus reactive power constant impedance withdrawals.
• bus_reactive_power_bounds::Matrix{Float64}: matrix containing upper and lower bounds for the reactive supply at each bus at each time period.
• bus_type::Matrix{PSY.ACBusTypes}: matrix containing type of buses present in the system.
• bus_magnitude::Matrix{Float64}: matrix containing the bus voltage magnitudes.
• bus_angles::Matrix{Float64}: matrix containing the bus voltage angles.
• arc_active_power_flow_from_to::Matrix{Float64}: matrix containing the active power flows measured at the from bus.
• arc_reactive_power_flow_from_to::Matrix{Float64}: matrix containing the reactive power flows measured at the from bus.
• arc_active_power_flow_to_from::Matrix{Float64}: matrix containing the active power flows measured at the to bus.
• arc_reactive_power_flow_to_from::Matrix{Float64}: matrix containing the reactive power flows measured at the to bus.
• arc_angle_differences::Matrix{Float64}: matrix containing the voltage angle difference (θfrom − θto) across each arc.
• generic_hvdc_flows::Dict{Tuple{Int, Int}, Tuple{Float64, Float64}}: dictionary mapping each generic HVDC line (represented as a tuple of the from and to bus numbers) to a tuple of (P_from_to, P_to_from) active power flows.
• bus_hvdc_net_power::Matrix{Float64}: "(b, t)" matrix containing the net power injections from all HVDC lines at each bus. b: number of buses, t: number of time period. Only contains HVDCs handled as separate injection/withdrawal pairs: LCCs and generic for DC, or just generic for AC.
• time_step_map::Dict{Int, S}: dictionary mapping the number of the time periods (corresponding to the column number of the previously mentioned matrices) and their names.
• power_network_matrix::M: matrix used for the evaluation of either the power flows or bus angles, depending on the method considered.
• aux_network_matrix::N: matrix used for the evaluation of either the power flows or bus angles, depending on the method considered.
• neighbors::Vector{Set{Int}}: Vector with the sets of adjacent buses.

PowerFlowData(
    pf::DCPowerFlow,
    sys::PSY.System
) -> ABAPowerFlowData

Creates a PowerFlowData structure configured for a standard DC power flow calculation, given the System sys. Configuration options like time_steps, time_step_names, network_reductions, and correct_bustypes are taken from the DCPowerFlow object.

Calling this function will not evaluate the power flows and angles. Note that first input is of type DCPowerFlow: this version is used to solve DC power flows, and returns an ABAPowerFlowData object.

Arguments:

• pf::DCPowerFlow: Run a DC power flow: internally, store the ABA matrix as power_network_matrix and the BA matrix as aux_network_matrix. Configuration options are taken from this object.
• sys::PSY.System: A System object that represents the power grid under consideration.

PowerFlowData(
    pf::PTDFDCPowerFlow,
    sys::PSY.System
) -> PTDFPowerFlowData

Creates a PowerFlowData structure configured for a Partial Transfer Distribution Factor Matrix DC power flow calculation, given the System sys. Configuration options like time_steps, time_step_names, network_reductions, and correct_bustypes are taken from the PTDFDCPowerFlow object.

Calling this function will not evaluate the power flows and angles. Note that first input is of type PTDFDCPowerFlow: this version is used to solve DC power flows via the Power Transfer Distribution Factor (PTDF) matrix. This function returns a PTDFPowerFlowData object.

Arguments:

• pf::PTDFDCPowerFlow: Run a DC power flow with PTDF matrix: internally, store the PTDF matrix as power_network_matrix and the ABA matrix as aux_network_matrix. Configuration options are taken from this object.
• sys::PSY.System: A System object that represents the power grid under consideration.

PowerFlowData(
    pf::AbstractACPowerFlow{<:ACPowerFlowSolverType},
    sys::PSY.System
) -> ACPowerFlowData{<:ACPowerFlowSolverType}

Creates the structure for an AC power flow calculation, given the System sys. Configuration options like time_steps, time_step_names, network_reductions, and correct_bustypes are taken from the AbstractACPowerFlow object (either ACPolarPowerFlow or ACRectangularPowerFlow).

Calling this function will not evaluate the power flows and angles. This version is used to solve AC power flows and returns an ACPowerFlowData object.

Arguments:

• pf::AbstractACPowerFlow: the settings for the AC power flow solver, including time_steps, time_step_names, network_reductions, and correct_bustypes.
• sys::PSY.System: A System object that represents the power grid under consideration.

WARNING: functions for the evaluation of the multi-period AC PF still to be implemented.

PowerFlowData(
    pf::vPTDFDCPowerFlow,
    sys::PSY.System
) -> vPTDFPowerFlowData

Creates a PowerFlowData structure configured for a virtual Partial Transfer Distribution Factor Matrix DC power flow calculation, given the System sys. Configuration options like time_steps, time_step_names, network_reductions, and correct_bustypes are taken from the vPTDFDCPowerFlow object.

Calling this function will not evaluate the power flows and angles. Note that first input is of type vPTDFDCPowerFlow: this version is used to solve DC power flows using a virtual Power Transfer Distribution Factor (PTDF) matrix. This function returns a vPTDFPowerFlowData object.

Arguments:

• pf::vPTDFDCPowerFlow: Run a virtual PTDF power flow: internally, store the virtual PTDF matrix as power_network_matrix and the ABA matrix as aux_network_matrix. Configuration options are taken from this object.
• sys::PSY.System: A System object that represents the power grid under consideration.

Sets the two PowerNetworkMatrix fields and a few others (time_steps, time_step_map), then creates arrays of default values (usually zeros) for the rest.

A PowerFlowContainer that represents its data as a PSY.System.

A type alias for a PowerFlowData struct whose type parameters are configured for the vPTDFDCPowerFlow method.

Compute arc angle differences for all arcs and all time steps, looking up bus indices from the arc axis and bus lookup stored in data. Used by DC solvers.

Compute arc angle differences using precomputed from/to bus index vectors over specified time steps. Used by the AC solver where fb_ix/tb_ix are already available from the branch flow calculation.

Create an appropriate PowerFlowContainer for the given PowerFlowEvaluationModel and initialize it from the given PSY.System.

Configuration options like time_steps, time_step_names, network_reductions, and correct_bustypes are taken from the PowerFlowEvaluationModel object.

Arguments:

• pfem::PowerFlowEvaluationModel: power flow model to construct a container for (e.g., DCPowerFlow())
• sys::PSY.System: the System from which to initialize the power flow container

Trait signifying whether the PowerFlowContainer can represent multi-period data. Must be implemented for all concrete subtypes.

solve_power_flow!(data::ACPowerFlowData; kwargs...)

Solve the multiperiod AC power flow problem for the given power flow data.

The bus types can be changed from PV to PQ if the reactive power limits are violated. The power flow solver settings are taken from the ACPowerFlow object stored in data.

Arguments

• data::ACPowerFlowData: The power flow data containing the grid information and initial conditions.
• kwargs...: Additional keyword arguments. If these overlap with those in the solver_settings of the ACPowerFlow object, the values in kwargs take precedence.

Keyword Arguments

• time_steps: Specifies the time steps to solve. Defaults to sorting and collecting the keys of get_time_step_map(data).

Description

This function solves the AC power flow problem for each time step specified in data. It preallocates memory for the results and iterates over the sorted time steps. For each time step, it calls the _ac_power_flow function to solve the power flow equations and updates the data object with the results. If the power flow converges, it updates the active and reactive power injections, as well as the voltage magnitudes and angles for different bus types (REF, PV, PQ). If the power flow does not converge, it sets the corresponding entries in data to NaN. Finally, it calculates the branch power flows and updates the data object.

Notes

• If the grid topology changes (e.g., tap positions of transformers or in-service status of branches), the admittance matrices Yft and Ytf must be updated.
• If Yft and Ytf change between time steps, the branch flow calculations must be moved inside the loop.

Examples

solve_power_flow!(data)

_get_arc_resistances(data::Union{PTDFPowerFlowData, vPTDFPowerFlowData, ABAPowerFlowData}) -> Vector{Float64}

Look up the equivalent resistance of each arc from the network reduction data. Delegates to _get_arc_branch_params and returns only the resistance vector.

Adjust the power injections vector to account for the power flows through LCCs.

Relies on the fact that we calculate those flows during initialization and save them to the active_power_flow_from_to and active_power_flow_to_from fields of the LCCParameters struct.

dc_loss_factors(
    data::Union{PTDFPowerFlowData, vPTDFPowerFlowData},
    P::Matrix{Float64},
) -> Matrix{Float64}

Compute the gradient of total system active power losses with respect to bus injections using the DC power flow approximation:

∂Loss/∂P = 2 · PTDFᵀ · diag(R) · PTDF · P

This is equivalent to the per-element form:

∂Loss/∂Pᵢ = Σₖ 2·Rₖ·PTDFₖᵢ·Σⱼ PTDFₖⱼ·Pⱼ

Arguments

• data::Union{PTDFPowerFlowData, vPTDFPowerFlowData}: solved power flow data containing the PTDF matrix and network reduction data for looking up branch resistances.
• P::Matrix{Float64}: bus injection matrix of size (num_buses, num_timesteps).

Returns

• Matrix{Float64}: loss factor matrix of size (num_buses, num_timesteps), where each entry [i, t] is the marginal change in total system losses per unit injection at bus i in time step t.

solve_power_flow!(data::PTDFPowerFlowData)

Evaluates the PTDF power flow and writes the result to the fields of the PTDFPowerFlowData structure.

This function modifies the following fields of data, setting them to the computed values:

• data.bus_angles: the bus angles for each bus in the system.
• data.branch_active_power_flow_from_to: the active power flow from the "from" bus to the "to" bus of each branch
• data.branch_active_power_flow_to_from: the active power flow from the "to" bus to the "from" bus of each branch

Additionally, it sets data.converged to true, indicating that the power flow calculation was successful.

solve_power_flow!(data::vPTDFPowerFlowData)

Evaluates the virtual PTDF power flow and writes the results to the fields of the vPTDFPowerFlowData structure.

This function modifies the following fields of data, setting them to the computed values:

• data.bus_angles: the bus angles for each bus in the system.
• data.branch_active_power_flow_from_to: the active power flow from the "from" bus to the "to" bus of each branch
• data.branch_active_power_flow_to_from: the active power flow from the "to" bus to the "from" bus of each branch

Additionally, it sets data.converged to true, indicating that the power flow calculation was successful.

_get_arc_branch_params(data) -> (rs, xs, taps, shifts)

Return vectors of series resistance, reactance, tap ratio, and phase shift for each arc, looking up every branch map in the NetworkReductionData.

This is the single source of truth for per-arc electrical parameters; _get_arc_resistances delegates to it.

_populate_loss_injections!(data::ABAPowerFlowData, sys::PSY.System)

Compute DCLF-style loss injections from the AC voltage profile stored in sys and add them as withdrawals in data.bus_active_power_withdrawals.

For each in-service branch k with series admittance y_k = g_k + j·b_k, the branch losses are computed using complex voltages:

P_loss_k = g_k · |V_i / tap_k − V_j|²

Losses are withdrawn at the sending-end bus (determined by power flow direction). This is a single-pass, non-iterative computation.

When the system is flat-start (V=1, θ=0 everywhere), all losses are zero and the method degenerates to standard lossless DCLF.

Partitions the state vector's variables based on what physical quantity each represents. Returns a NamedTuple, with the 4 keys Va, Vm, P, and Q. The 4 values are vectors of length equal to the number of buses, with NaNs in the positions where that physical quantity is not part of the state vector for that bus. (Currently not intended for use in spots where performance is critical.)

Update the fields of data based on the values of the state vector.

Update P_net and Q_net based on the values of the state vector.

Update state vector based on values of fields of data.

Sets the fields of a PowerFlowData struct to match the given System.

When solving AC power flows, if the initial guess has large residual, we run a DC power flow as a fallback. This runs a DC power flow on data::ACPowerFlowData for the given time_step, and writes the solution to data.bus_angles.

Rectangular/MCPB analog of _enhanced_flat_start: per subnetwork, set PV/PQ bus angles to the mean REF-bus angle and PQ magnitudes to the mean PV setpoint magnitude, written back as (e, f) = (Vm·cosθ, Vm·sinθ). PV buses keep their setpoint magnitude (only the angle changes); REF blocks and the PV Q / REF (P,Q) slots are left as in x0. Uses residual.subnetworks (ref-bus index → member bus indices) for the partition. Identical for the rectangular and MCPB layouts (both use 2-slot (e, f) PV/PQ blocks and never touch a PV Q slot).

Use state variables from a previous converged time step (prev) as a candidate starting point.

Calculate x0 from data.

If initial residual is large, run a DC power flow and see if that gives a better starting point for angles. If so, then overwrite x0 with the result of the DC power flow. If not, keep the original x0.

MCPB analog of the rectangular improve_x0: base flat start (via mixed_initial_state!) → previous-converged-timestep warm start (via _mixed_fill_state!) → enhanced flat start (gated on get_enhanced_flat_start(pf)) → large-residual warning. Mirrors the rectangular path verbatim, swapping rect_initial_state!→mixed_initial_state!, _rect_fill_state!→_mixed_fill_state!, and the rectangular _enhanced_flat_start→the MCPB ACMixedCPBResidual overload. No DC robust fallback: ACMixedPowerFlow has no robust_power_flow field (get_robust_power_flow(::AbstractACPowerFlow) == false).

Rectangular analog of the polar improve_x0: base flat start → previous-converged-timestep warm start → enhanced flat start (gated on get_enhanced_flat_start(pf)) → large-residual warning. No DC robust fallback: ACRectangularPowerFlow has no robust_power_flow field and a CI-aware DC fallback is out of scope (see the formulation/solver-split spec).

_calculate_dP_dV_lcc(t, I_dc, x_t, Vm, ϕ) -> Float64

True-ϕ derivative of P_lcc = Vm · t · √6/π · I_dc · cos(ϕ(Vm, t, α)) with respect to Vm. Two regimes:

• Interior (ϕ unclamped): ∂ϕ/∂Vm = -∂raw/∂Vm / sin(ϕ) is nonzero; the sin(ϕ) factor from the chain rule cancels the -sin(ϕ) from differentiating cos(ϕ), giving the second (chain) term below.
• Clamp (sin(ϕ) ≈ 0, i.e. ϕ ∈ {0, π}): ϕ is locally pinned (∂ϕ/∂x = 0) and the residual sees only the leading Vm · cos(ϕ) dependence on Vm. The chain term must be dropped — otherwise the analytic Jacobian disagrees with the residual at the clamp, exactly analogous to the sin(ϕ)→0 guard in _calculate_dQ_dV_lcc.

Caller passes I_dc > 0 and the side-specific ϕ. Rectifier: phi_r. Inverter: phi_i (already encodes the sign convention via _calculate_ϕ_lcc(-I_dc, …); positive I_dc is still passed here).

_calculate_dP_dt_lcc(t, I_dc, x_t, Vm, ϕ) -> Float64

True-ϕ derivative of P_lcc with respect to the transformer tap t. Same two-regime structure as _calculate_dP_dV_lcc: chain term only when unclamped (sin(ϕ) ≥ LCC_sinϕ_TOLERANCE); leading Vm · cos(ϕ) term always present.

_calculate_dP_dα_lcc(t, I_dc, Vm, α, ϕ) -> Float64

True-ϕ derivative of P_lcc with respect to the firing/extinction angle α, rectifier sign convention. In the interior, ∂ϕ/∂α = sin(α)/sin(ϕ) and combines with the -sin(ϕ) from differentiating cos(ϕ) to give the closed form below (no sin(ϕ) in the result). At the clamp, ∂ϕ/∂α = 0 and the true derivative is zero — same boundary handling as the dQ helpers. Inverter callers must negate the helper output (the inverter ϕ convention flips ∂ϕ_i/∂α_i).

_calculate_dQ_dV_lcc(t::Float64, I_dc::Float64, x_t::Float64, Vm::Float64, ϕ::Float64) -> Float64

Compute the derivative of reactive power Q with respect to voltage magnitude Vm for LCC converter calculations.

_calculate_dQ_dt_lcc(t::Float64, I_dc::Float64, x_t::Float64, Vm::Float64, ϕ::Float64) -> Float64

Compute the derivative of reactive power Q with respect to transformer tap t for LCC converter calculations.

_calculate_dQ_dα_lcc(t::Float64, I_dc::Float64, x_t::Float64, Vm::Float64, ϕ::Float64, α::Float64) -> Float64

Compute the derivative of reactive power Q with respect to firing/extinction angle α for LCC converter calculations.

_calculate_y_lcc(t::Float64, I_dc::Float64, Vm::Float64, ϕ::Float64) -> ComplexF64

Compute the admittance value Y for LCC converter calculations.

_calculate_ϕ_lcc(α::Float64, I_dc::Float64, x_t::Float64, Vm::Float64) -> Float64

Compute the phase angle ϕ for LCC converter calculations.

_dphi_dV_lcc(x_t, I_dc, V, t, ϕ) -> Float64

∂ϕ/∂V with sin(ϕ) → 0 clamp guard returning 0. In the interior, ∂ϕ/∂V = -∂raw/∂V / sin(ϕ) = -x_t·I_dc / (√2·V²·t·sin(ϕ)). At the clamp ϕ is pinned (∂ϕ/∂V = 0). Same form on both sides — the inverter passes the same positive I_dc, only its ϕ differs.

_dphi_dt_lcc(x_t, I_dc, V, t, ϕ) -> Float64

∂ϕ/∂t (tap) with clamp guard. -x_t·I_dc / (√2·V·t²·sin(ϕ)) in the interior, 0 at the clamp.

_dphi_dα_lcc(α, ϕ) -> Float64

∂ϕ/∂α (rectifier sign) with clamp guard. sin(α)/sin(ϕ) in the interior, 0 at the clamp. Inverter convention flips the sign — callers on the inverter side negate the helper output.

_lcc_i_dc_from_p_set(r, p) -> Float64

DC current set point from the active-power set point p and total DC-side resistance r: the positive root of r·I² + I − p = 0. Have to special-case r == 0, where the quadratic formula gives I = 0/0 = NaN instead of I = p.

_lcc_jacobian_scalars(data, i, time_step, Vm_fb, Vm_tb)

Precompute the scalar coefficients used by both the polar and the rectangular LCC Jacobian assembly for LCC i at time_step. Vm_fb / Vm_tb are the AC-side voltage magnitudes — polar reads them from data.bus_magnitude; rectangular computes sqrt(e² + f²) from state.

The returned NamedTuple includes the tail-row × {tail-column, bus-V} entries that are shared between formulations. These are computed via the true-ϕ helpers (_calculate_dP_dt_lcc, _calculate_dP_dα_lcc, _calculate_dP_dV_lcc), which apply the sin(ϕ) → 0 boundary guard: in the interior the algebraic identity makes the result equal to the α-approximation form, and at the clamp the guard correctly drops the chain term so the Jacobian matches the residual (which sees ∂ϕ/∂x = 0 at the clamp).

The P-setpoint row F_t_fb depends on the rectifier-side state (V_fb, tap_r, α_r) when data.lcc.setpoint_at_rectifier[i] and on the inverter-side state (V_tb, tap_i, α_i) otherwise; the helper branches on that flag and zeroes the inactive side so each assembly writes its F_t_fb slots unconditionally.

_set_lcc_tail_residuals!(F, data, base_offset, time_step) [polar]
_set_lcc_tail_residuals!(F, data, base_offset, time_step, e_state, f_state) [rect]

Write the 4 LCC tail residual rows (P-setpoint, DC-line balance, two α limit constraints) for each LCC into F, starting at slot base_offset + 1. The i-th LCC occupies slots base_offset + 4(i-1) + 1 .. base_offset + 4i. The polar method reads |V| from data.bus_magnitude; the rectangular method reads sqrt(e² + f²) from the (e, f) state (since data.bus_magnitude holds V_set at PV buses, not the actual state magnitude). Mirrors the two-method layout of _update_ybus_lcc!.

_update_ybus_lcc!(data, time_step, e_state, f_state)

Rectangular variant: reads |V| at each AC terminal from sqrt(e_state[i]^2 + f_state[i]^2) so the LCC math stays consistent with the rectangular CI residual / Jacobian (which operate on (e, f) instead of (|V|, θ)).

_update_ybus_lcc!(data, time_step)

Recompute data.lcc.rectifier.phi, data.lcc.inverter.phi, and data.lcc.branch_admittances for each LCC at time_step. Reads |V| at each AC terminal from data.bus_magnitude (the polar convention). The (e_state, f_state) method below covers the rectangular CI case where |V_state| = sqrt(e² + f²) must be used instead — at PV buses, data.bus_magnitude holds V_set rather than the actual state magnitude.

Adjust the power injections/withdrawal vectors to account for all HVDC lines of a given type, modeling those HVDC lines as a simple fixed injection/withdrawal at each terminal.

Initialize the arcs and bus_indices fields of the LCCParameters structure in the PowerFlowData.

struct ACPowerFlowResidual

A struct to keep track of the residuals in the Newton-Raphson AC power flow calculation.

Fields

• data::ACPowerFlowData: The grid model data.
• Rf!::Function: A function that updates the residuals based on the latest values stored in the grid at the given iteration.
• Rv::Vector{Float64}: A vector of the values of the residuals.
• P_net::Vector{Float64}: A vector of net active power injections.
• Q_net::Vector{Float64}: A vector of net reactive power injections.
• P_net_set::Vector{Float64}: A vector of the set-points for active power injections (their initial values before power flow calculation).
• bus_slack_participation_factors::SparseVector{Float64, Int}: A sparse vector of the slack participation factors aggregated at the bus level.
• subnetworks::Dict{Int64, Vector{Int64}}: The dictionary that identifies subnetworks (connected components), with the key defining the REF bus, values defining the corresponding buses in the subnetwork.
• P_slack_buf::Vector{Float64}: Scratch buffer of length n_buses used by _update_residual_values! to write the per-subnetwork slack distribution in place, avoiding a per-iteration allocation when indexing bus_slack_participation_factors by subnetwork_buses.
• validate_indices::Vector{Int}: precomputed x-indices of PQ-bus |V| entries for the per-iteration voltage-magnitude diagnostic.

ACPowerFlowResidual(data::ACPowerFlowData, time_step::Int64)

Create an instance of ACPowerFlowResidual for a given time step.

Arguments

• data::ACPowerFlowData: The power flow data representing the power system model.
• time_step::Int64: The time step for which the power flow calculation is executed.

Returns

• ACPowerFlowResidual: An instance containing the residual values, net bus active power injections, and net bus reactive power injections.

(Residual::ACPowerFlowResidual)(x::Vector{Float64}, time_step::Int64)

Update the AC power flow residuals inplace and store the result in the attribute Rv of the struct. The inputs are the values of state vector x and the current time step time_step. This function implements the functor approach for the ACPowerFlowResidual struct. This makes the struct callable. Calling the ACPowerFlowResidual will also update the values of P, Q, V, Θ in the data struct.

Arguments

• x::Vector{Float64}: The state vector values.
• time_step::Int64: The current time step.

(Residual::ACPowerFlowResidual)(Rv::Vector{Float64}, x::Vector{Float64}, time_step::Int64)

Evaluate the AC power flow residuals and store the result in Rv using the provided state vector x and the current time step time_step. The residuals are updated inplace in the struct and additionally copied to the provided array. This function implements the functor approach for the ACPowerFlowResidual struct. This makes the struct callable. Calling the ACPowerFlowResidual will also update the values of P, Q, V, Θ in the data struct.

Arguments

• Rv::Vector{Float64}: The vector to store the calculated residuals.
• x::Vector{Float64}: The state vector.
• time_step::Int64: The current time step.

_build_bus_slack_participation_factors(data, bus_type, subnetworks, time_step)

Collect the per-bus generator-slack-participation factors (REF and PV buses only), validate that the sum is positive and no value is negative, and normalize so that each subnetwork's participating buses sum to 1. Returns a SparseVector{Float64, Int} of length n_buses.

Shared between the polar ACPowerFlowResidual and the rectangular current-injection (CI) residual (ACRectangularCIResidual) constructors — both need identical slack-distribution semantics.

_update_residual_values!(
    F::Vector{Float64},
    x::Vector{Float64},
    P_net::Vector{Float64},
    Q_net::Vector{Float64},
    data::ACPowerFlowData,
    time_step::Int64,
)

Update the residual values for the Newton-Raphson AC power flow calculation. This function is used internally in the ACPowerFlowResidual struct. This function also updates the values of P, Q, V, Θ in the data struct.

Arguments

• F::Vector{Float64}: Vector of the values of the residuals.
• x::Vector{Float64}: State vector values.
• P_net::Vector{Float64}: Vector of net active power injections at each bus.
• Q_net::Vector{Float64}: Vector of net reactive power injections at each bus.
• P_net_set::Vector{Float64}: Vector of the set-points for active power injections (their initial values before power flow calculation).
• bus_slack_participation_factors::SparseVector{Float64, Int}: Sparse vector of the slack participation factors aggregated at the bus level.
• ref_bus::Int: The index of the reference bus to be used for the total slack power.
• data::ACPowerFlowData: Data structure representing the grid model for the AC power flow calculation.
• time_step::Int64: The current time step for which the residual values are being updated.

struct ACPowerFlowJacobian

A struct that represents the Jacobian matrix for AC power flow calculations.

This struct uses the functor pattern, meaning instances of ACPowerFlowJacobian store the data (Jacobian matrix) internally and can be called as a function at the same time. Calling the instance as a function updates the stored Jacobian matrix.

Fields

• data::ACPowerFlowData: The grid model data used for power flow calculations.
• Jf!::Function: A function that calculates the Jacobian matrix inplace.
• Jv::SparseArrays.SparseMatrixCSC{Float64, Int32}: The Jacobian matrix, which is updated by the function Jf!.
• diag_elements::MVector{4, Float64}: Temporary storage for diagonal elements during Jacobian update.
• bus_slack_participation_factors::SparseVector{Float64, Int}: Normalized per-bus slack participation factors for the current time step (from the ACPowerFlowResidual). Used for the distributed slack Jacobian entries.
• subnetworks::Dict{Int64, Vector{Int64}}: Subnetwork mapping from REF bus to bus list (from the ACPowerFlowResidual). Used for the distributed slack Jacobian entries.

(J::ACPowerFlowJacobian)(time_step::Int64)

Update the Jacobian matrix Jv using the function Jf! and the provided data and time step.

Defining this method allows an instance of ACPowerFlowJacobian to be called as a function, following the functor pattern.

Arguments

• time_step::Int64: The time step for the calculations.

Example

residual = ACPowerFlowResidual(data, time_step)
J = ACPowerFlowJacobian(residual, time_step)
J(time_step)  # Updates the Jacobian matrix Jv

ACPowerFlowJacobian(residual::ACPowerFlowResidual, time_step::Int64) -> ACPowerFlowJacobian

Constructor for ACPowerFlowJacobian. The returned instance has its sparsity pattern initialized and shares the residual's slack-participation, subnetwork, and ZIP-coefficient caches — the residual must be constructed first against the same data and time_step.

Arguments

• residual::ACPowerFlowResidual: The companion residual; supplies data, bus_slack_participation_factors, subnetworks, and the per-bus ZIP load coefficient vectors.
• time_step::Int64: The time step for the calculations.

Example

residual = ACPowerFlowResidual(data, time_step)
J = ACPowerFlowJacobian(residual, time_step)
J(time_step)  # Updates the Jacobian matrix stored internally in J.
J.Jv  # Access the Jacobian matrix stored internally in J.

(J::ACPowerFlowJacobian)(J::SparseArrays.SparseMatrixCSC{Float64, Int32}, time_step::Int64)

Use the ACPowerFlowJacobian to update the provided Jacobian matrix J inplace.

Update the internally stored Jacobian matrix Jv using the function Jf! and the provided data and time step, and write the updated Jacobian values to J.

This method allows an instance of ACPowerFlowJacobian to be called as a function, following the functor pattern.

Arguments

• J::SparseArrays.SparseMatrixCSC{Float64, Int32}: A sparse matrix to be updated with new values of the Jacobian matrix.
• time_step::Int64: The time step for the calculations.

Example

residual = ACPowerFlowResidual(data, time_step)
J = ACPowerFlowJacobian(residual, time_step)
Jv = SparseArrays.sparse(Float64[], J_INDEX_TYPE[], J_INDEX_TYPE[])
J(Jv, time_step)  # Updates the Jacobian matrix Jv and writes it to J

_block_J_indices(data::ACPowerFlowData, time_step::Int) -> (Vector{Int32}, Vector{Int32})

Get the indices to reindex the Jacobian matrix from the interleaved form to the block form:

\[\begin{bmatrix} \frac{\partial P}{\partial \theta} & \frac{\partial P}{\partial V} \\ \frac{\partial Q}{\partial \theta} & \frac{\partial Q}{\partial V} \end{bmatrix}\]

Arguments

• pvpq::Vector{Int32}: Indices of the buses that are PV or PQ buses.
• pq::Vector{Int32}: Indices of the buses that are PQ buses.

Returns

• rows::Vector{Int32}: Row indices for the block Jacobian matrix.
• cols::Vector{Int32}: Column indices for the block Jacobian matrix.

calculate_loss_factors(data::ACPowerFlowData, Jv::SparseMatrixCSC{Float64, Int32}, time_step::Int)

Calculate and store the active power loss factors in the loss_factors matrix of the ACPowerFlowData structure for a given time step.

The loss factors are computed using the Jacobian matrix Jv and the vector dSbus_dV_ref, which contains the partial derivatives of slack power with respect to bus voltages. The function interprets changes in slack active power injections as indicative of changes in grid active power losses. KLU is used to factorize the sparse Jacobian matrix to solve for the loss factors.

Arguments

• data::ACPowerFlowData: The data structure containing power flow information, including the loss_factors matrix.
• Jv::SparseMatrixCSC{Float64, Int32}: The sparse Jacobian matrix of the power flow system.
• time_step::Int: The time step index for which the loss factors are calculated.

calculate_voltage_stability_factors(data::ACPowerFlowData, J::ACPowerFlowJacobian, time_step::Integer)

Calculate and store the voltage stability factors in the voltage_stability_factors matrix of the ACPowerFlowData structure for a given time step. The voltage stability factors are computed using the Jacobian matrix J in block format after a converged power flow calculation. The results are stored in the voltage_stability_factors matrix in the data instance. The factor for the grid as a whole (σ) is stored in the position of the REF bus. The values of the singular vector v indicate the sensitivity of the buses and are stored in the positions of the PQ buses. The values of v for PV buses are set to zero. The function uses the method described in "Fast calculation of a voltage stability index" by PA Lof et. al.

Arguments

• data::ACPowerFlowData: The instance containing the grid model data.
• J::ACPowerFlowJacobian: The Jacobian matrix cache.
• time_step::Integer: The calculated time step.

_create_jacobian_matrix_structure(data::ACPowerFlowData, time_step::Int64) -> SparseMatrixCSC{Float64, Int32}

Create the structure of the Jacobian matrix for an AC power flow problem.

Arguments

• data::ACPowerFlowData: The power flow model.
• time_step::Int64: The specific time step for which the Jacobian matrix structure is created.

Returns

• SparseMatrixCSC{Float64, Int32}: A sparse matrix with structural zeros representing the structure of the Jacobian matrix.

Description

This function initializes the structure of the Jacobian matrix for an AC power flow problem. The Jacobian matrix is used in power flow analysis to represent the partial derivatives of bus active and reactive power injections with respect to bus voltage magnitudes and angles.

Unlike some commonly used approaches where the Jacobian matrix is constructed as four submatrices, each grouping values for the four types of partial derivatives, this function groups the partial derivatives by bus. The structure is organized as groups of 4 values per bus.

This approach is more memory-efficient. Furthermore, this structure results in a more efficient factorization because the values are more likely to be grouped close to the diagonal. Refer to Electric Energy Systems: Analysis and Operation by Antonio Gomez-Exposito and Fernando L. Alvarado for more details.

The function initializes three arrays (rows, columns, and values) to store the row indices, column indices, and values of the non-zero elements of the Jacobian matrix, respectively.

For each bus in the system, the function iterates over its neighboring buses and determines the type of each neighboring bus (REF, PV, or PQ). Depending on the bus type, the function adds the appropriate entries to the Jacobian matrix structure.

• For REF buses, entries are added for local active and reactive power.
• For PV buses, entries are added for active and reactive power with respect to angle, and for local reactive power.
• For PQ buses, entries are added for active and reactive power with respect to voltage magnitude and angle.

Example Structure

For a system with 3 buses where bus 1 is REF, bus 2 is PV, and bus 3 is PQ:

Let $\Delta P_j$, $\Delta Q_j$ be the active, reactive power balance at the $j$th bus. Let $P_j$ and $Q_j$ be the active and reactive power generated at the $j$th bus (REF and PV only). The state vector is $x = [P_1, Q_1, Q_2, \theta_2, V_3, \theta_3]$, and the residual vector is $F(x) = [\Delta P_1, \Delta Q_1, \Delta P_2, \Delta Q_2, \Delta P_3, \Delta Q_3]$.

The Jacobian matrix $J = \nabla F(x)$ has the structure:

\[J = \begin{bmatrix} \frac{\partial \vec{F}}{\partial P_1} & \frac{\partial \vec{F}}{\partial Q_1} & \frac{\partial \vec{F}}{\partial Q_2} & \frac{\partial \vec{F}}{\partial \theta_2} & \frac{\partial \vec{F}}{\partial V_3} & \frac{\partial \vec{F}}{\partial \theta_3} \end{bmatrix}\]

In reality, for large networks, this matrix would be sparse, and each 2×2 block would only be nonzero when there's a line between the respective buses.

Finally, the function constructs a sparse matrix from the collected indices and values and returns it.

Create the Jacobian matrix structure for a PV bus. Currently unused: we fill all four values even for PV buses with structiural zeros using the same function as for PQ buses.

Create the Jacobian matrix structure for a reference bus (REF). Currently unused: we fill all four values even for PV buses with structiural zeros using the same function as for PQ buses.

Create the Jacobian matrix structure for a PQ bus. Using this for all buses because a) for REF buses it doesn't matter if there are 2 values or 4 values - there are not many of them in the grid b) for PV buses we fill all four values because we can have a PV -> PQ transition and then we need to fill all four values

_create_jacobian_matrix_structure_lcc(
    data::ACPowerFlowData,
    rows::Vector{Int32},
    columns::Vector{Int32},
    values::Vector{Float64},
    num_buses::Int
)

Create the Jacobian matrix structure for LCC HVDC systems.

Description

The function iterates over each LCC system and adds the non-zero entries to the Jacobian matrix structure. The state vector for every LCC contains 4 variables: tap position and thyristor angle for both the rectifier and inverter sides. The indices of non-zero entries correspond to the positions of these variables in the extended state vector.

For an LCC system connecting bus $i$ (rectifier side) and bus $j$ (inverter side), the state variables are:

• $t_i$: tap position at rectifier
• $t_j$: tap position at inverter
• $\alpha_i$: thyristor angle at rectifier
• $\alpha_j$: thyristor angle at inverter

The residuals include:

• $F_{t_i}$: Active power balance at rectifier (controls $P_i$ to match setpoint)
• $F_{t_j}$: Total active power balance across LCC system
• $F_{\alpha_i}$: Rectifier thyristor angle constraint (maintains $\alpha_i$ at minimum)
• $F_{\alpha_j}$: Inverter thyristor angle constraint (maintains $\alpha_j$ at minimum)

Example Structure

For a system with 2 buses connected by one LCC where bus 1 is the rectifier side and bus 2 is the inverter side, the Jacobian matrix would have non-zero entries at positions like:

\[\begin{array}{c|cccccccc} & V_1 & \delta_1 & V_2 & \delta_2 & t_1 & t_2 & \alpha_1 & \alpha_2 \\ \hline P_1 & \frac{\partial P_1}{\partial V_1} & & & & \frac{\partial P_1}{\partial t_1} & & \frac{\partial P_1}{\partial \alpha_1} & \\ Q_1 & \frac{\partial Q_1}{\partial V_1} & & & & \frac{\partial Q_1}{\partial t_1} & & \frac{\partial Q_1}{\partial \alpha_1} & \\ P_2 & & & & & & & & \\ Q_2 & & & & & & & & \\ F_{t_1} & \frac{\partial F_{t_1}}{\partial V_1} & & & & \frac{\partial F_{t_1}}{\partial t_1} & & \frac{\partial F_{t_1}}{\partial \alpha_1} & \\ F_{t_2} & \frac{\partial F_{t_2}}{\partial V_1} & & \frac{\partial F_{t_2}}{\partial V_2} & & \frac{\partial F_{t_2}}{\partial t_1} & \frac{\partial F_{t_2}}{\partial t_2} & \frac{\partial F_{t_2}}{\partial \alpha_1} & \frac{\partial F_{t_2}}{\partial \alpha_2} \\ F_{\alpha_1} & & & & & & & \frac{\partial F_{\alpha_1}}{\partial \alpha_1} & \\ F_{\alpha_2} & & & & & & & & \frac{\partial F_{\alpha_2}}{\partial \alpha_2} \end{array}\]

This function sets up the indices of these non-zero entries in the sparse Jacobian matrix structure.

Arguments

• data::ACPowerFlowData: The power flow data containing LCC system information.
• rows::Vector{Int32}: Vector to store row indices of non-zero Jacobian entries.
• columns::Vector{Int32}: Vector to store column indices of non-zero Jacobian entries.
• values::Vector{Float64}: Vector to store initial values of non-zero Jacobian entries.
• num_buses::Int: Total number of buses in the system.

_singular_value_decomposition(J::SparseMatrixCSC{Float64, Int32}, npvpq::Integer; tol::Float64 = 1e-9, max_iter::Integer = 100,)

Estimate the smallest singular value σ and corresponding left and right singular vectors u and v of a sparse matrix G_s (a sub-matrix of J). This function uses an iterative method involving LU factorization of the Jacobian matrix to estimate the smallest singular value of G_s. The algorithm alternates between updating u and v, normalizing, and checking for convergence based on the change in the estimated singular value σ. The function uses the method described in Algorithm 3 of "Fast calculation of a voltage stability index" by PA Lof et. al.

Arguments

• J::SparseMatrixCSC{Float64, Int32}: The sparse block-form Jacobian matrix.
• npvpq::Integer: Number of PV and PQ buses in J.

Keyword Arguments

• tol::Float64=1e-9: Convergence tolerance for the iterative algorithm.
• max_iter::Integer=100: Maximum number of iterations.

Returns

• σ::Float64: The estimated smallest singular value.
• left::Vector{Float64}: The estimated left singular vector (referred to as u in the cited paper).
• right::Vector{Float64}: The estimated right singular vector (referred to as v in the cited paper).

Used to update Jv based on the bus voltages, angles, etc. in data.

_rect_fill_state!(x, data, bus_state_offset, type_time_step, value_time_step)

Fill the rectangular state vector x. The state-block layout (which buses are REF/PV/PQ and therefore the 2- vs 3-slot blocks) is taken from data.bus_type[:, type_time_step]; the values (voltages, injections, LCC taps/angles) are read from *[:, value_time_step].

With type_time_step == value_time_step this is the plain flat start. With value_time_step pointing at a previously converged step it produces the previous-solution warm-start candidate while keeping the offsets valid for the current step (the rectangular analog of polar _previous_solution_start / update_state!).

compute_bus_state_offsets(bus_type)

Compute per-bus state-vector offsets and block sizes for the augmented current-injection (rectangular) formulation. PQ and REF buses occupy 2 entries each (e, f) or (P_gen, Q_gen); PV buses occupy 3 entries (e, f, Q).

Returns (offsets, block_sizes, total_bus_state) where

• offsets[i] is the 1-based start index of bus i's block in the state vector
• offsets[end] is the start of the LCC tail (== total_bus_state + 1)
• block_sizes[i] ∈ {2, 3}
• total_bus_state is the total count of bus-state slots (excluding LCC tail)

fold_zip_constant_z!(Y_bus_eff, data, time_step)

Add the constant-impedance ZIP load components into the Y_bus_eff diagonal as fixed shunt admittances. Sienna's ZIP load model parameterizes the load at |V| = 1.0 pu directly: a load with constant_impedance_active_power = β_P and constant_impedance_reactive_power = β_Q draws S = (β_P + jβ_Q)·|V|². As a shunt admittance this is Y_sh = (β_P − jβ_Q) because |V|²·conj(Y_sh) = (β_P + jβ_Q)·|V|².

rect_finalize_bus_injections!(data, x, bus_state_offset, P_net_set,
                              bus_slack_participation_factors, subnetworks, time_step)

Distribute the converged subnetwork slack across participating buses and write bus_active_power_injections and bus_reactive_power_injections accordingly.

Mirrors polar _set_state_variables_at_bus! semantics: at every participating REF and PV bus, P_gen = P_net_set[i] + c_i · P_slack_total, where P_slack_total = x[ref_off] - P_net_set[ref_bus]. At REF, Qgen is taken from x[off + 1]; at PV, Qgen is taken from x[off + 2]. At PQ buses no slack attribution is needed: bus_active_power_injections / bus_reactive_power_injections already hold the load setpoint from PowerFlowData construction.

Called once per time step after the NR loop converges (not on every iteration), because the slack distribution is only meaningful at the converged x.

rect_initial_state!(x, data, bus_state_offset, bus_block_size, time_step)

Initialize the state vector x from data.bus_magnitude, data.bus_angles, and the bus power-injection fields, plus the LCC tap/angle fields. Counterpart of rect_update_data!. At REF buses, the first two slots hold (P_gen, Q_gen) (including any distributed-slack increment); elsewhere the first two slots hold (e, f), and PV buses' third slot holds Q_gen.

rect_update_data!(data, x, bus_state_offset, bus_block_size, time_step)

Write the state-derived voltage fields (bus_magnitude, bus_angles) and LCC taps/angles from x back into data. Per-iteration helper invoked by the rectangular CI residual.

Does NOT write bus_active_power_injections / bus_reactive_power_injections: those are finalized once after convergence with the correct distributed-slack share by rect_finalize_bus_injections!. At REF buses x[off] carries the entire subnetwork slack and would over-attribute it; at PV buses the gen Q has not yet been combined with the per-bus slack share.

Counterpart of rect_initial_state!.

struct ACRectangularCIResidual

Residual functor for the augmented current-injection (rectangular) AC power flow. Mirrors ACPowerFlowResidual but operates on the per-bus variable-block state representation: PQ/REF blocks are 2 entries (e,f) or (P_gen, Q_gen); PV blocks are 3 entries (e, f, Q).

Fields

• data::ACPowerFlowData
• Rf!::Function — inplace residual update
• Rv::Vector{Float64} — current residual values, length total_bus_state + 4·n_LCC
• Y_bus_eff::SparseMatrixCSC{ComplexF64, Int} — Y_bus with ZIP constant-Z folded in
• P_net_const::Vector{Float64} — constant-power net injection (no |V| dependence)
• Q_net_const::Vector{Float64} — constant-power net reactive injection
• const_I_P::Vector{Float64} — constant-current P-withdrawal coefficient per bus
• const_I_Q::Vector{Float64} — constant-current Q-withdrawal coefficient per bus
• P_net_set::Vector{Float64} — initial P_net for distributed-slack delta computation
• bus_slack_participation_factors::SparseVector{Float64, Int}
• subnetworks::Dict{Int64, Vector{Int64}}
• bus_state_offset::Vector{REC_INDEX_TYPE}
• bus_block_size::Vector{Int8}
• total_bus_state::Int
• validate_offsets::Vector{Int} — precomputed x-offsets of PQ/PV buses for the per-iteration voltage-magnitude diagnostic

Update residual values F for the augmented current-injection formulation.

Strategy: walk Ybuseff once to accumulate Y·V into the F slots, then add the specified-current contribution per bus type. PV buses add the ΔV² row. The 4 LCC tail residuals are appended. ZIP constant-Z is already folded into Ybuseff at setup, so it contributes via Y·V. ZIP constant-current is subtracted from the effective P/Q before computing I_spec.

struct ACRectangularCIJacobian

Jacobian functor for the augmented current-injection (rectangular) AC power flow. Mirrors ACPowerFlowJacobian but operates on the per-bus variable-block state representation.

Per-iteration updates write directly into nonzeros(Jv) via nzval-index caches built once at construction time, so the hot-path cost is O(N + n_LCC) rather than O((N + n_LCC) · log(nnz_per_col)) of Jv[r, c] = v setindex.

Fields

• data::ACPowerFlowData
• Jf!::Function — inplace Jacobian update
• Jv::SparseMatrixCSC{Float64, J_INDEX_TYPE} — Jacobian values
• Y_bus_eff::SparseMatrixCSC{ComplexF64, Int} — Y_bus with ZIP-Z folded in
• Y_diag::Vector{ComplexF64} — cached Ybuseff diagonal
• Bus-diagonal nzval caches diag_base_nz (4×nbuses), `pvextranz` (4×nbuses with sentinel 0 for non-PV buses)
• Slack cross-term nzval caches slack_nz_idx_e, slack_nz_idx_f, plus slack_bus_k / slack_c_k for the corresponding per-iteration data
• LCC tail nzval cache lcc_nz (24×n_lccs; the last 2 identity diagonals stay 1.0)

_build_diag_nz_cache(Jv, bus_state_offset, bus_types)

Return (diag_base_nz::Matrix{Int}, pv_extra_nz::Matrix{Int}), each 4 × n_buses, containing the nzval indices for the per-bus diagonal block entries.

Row layout of diag_base_nz (always populated, all bus types): 1: Jv[off, off], 2: Jv[off, off+1], 3: Jv[off+1, off], 4: Jv[off+1, off+1]

Row layout of pv_extra_nz (only populated for PV; 0 for non-PV): 1: Jv[off, off+2], 2: Jv[off+1, off+2], 3: Jv[off+2, off], 4: Jv[off+2, off+1]

_build_lcc_nz_cache(Jv, data, bus_state_offset, total_bus_state, n_lccs)

Return a 24 × n_lccs matrix of nzval indices for the per-LCC tail entries that get updated each iteration. The two identity diagonals (Jv[idx_alpha_r, idx_alpha_r] and Jv[idx_alpha_i, idx_alpha_i]) are not included — they are set to 1.0 at structure-build time and never updated. The 8 FB/TB-side diagonal-block overlay entries are NOT included either — they share nzval slots with diag_base_nz for buses fb and tb and are addressed through that cache.

Rows 9, 10, 15, 16, 21–24 belong to the P-setpoint row F_t_fb (idx_tap_r): rows 9, 10, 15, 16 hold its rectifier-side dependence and rows 21–24 its inverter-side dependence; _lcc_jacobian_scalars zeroes whichever side the set point is not on.

Row layout (matches order pushed by [_create_rect_ci_lcc_structure!]): 1: Jv[colefb, idxtapr], 2: Jv[colefb, idxalphar], 3: Jv[colffb, idxtapr], 4: Jv[colffb, idxalphar], 5: Jv[coletb, idxtapi], 6: Jv[coletb, idxalphai], 7: Jv[colftb, idxtapi], 8: Jv[colftb, idxalphai], 9: Jv[idxtapr, colefb], 10: Jv[idxtapr, colffb], 11: Jv[idxtapi, colefb], 12: Jv[idxtapi, colffb], 13: Jv[idxtapi, coletb], 14: Jv[idxtapi, colftb], 15: Jv[idxtapr, idxtapr], 16: Jv[idxtapr, idxalphar], 17: Jv[idxtapi, idxtapr], 18: Jv[idxtapi, idxtapi], 19: Jv[idxtapi, idxalphar], 20: Jv[idxtapi, idxalphai], 21: Jv[idxtapr, coletb], 22: Jv[idxtapr, colftb], 23: Jv[idxtapr, idxtapi], 24: Jv[idxtapr, idxalphai],

_build_slack_nz_cache(Jv, bus_state_offset, subnetworks, bus_slack_participation_factors)

Return (slack_nz_idx_e, slack_nz_idx_f, slack_bus_k, slack_c_k). Each entry corresponds to one (busk != refbus, ck != 0) slack cross-term. The nzval indices point at `Jv[koff, refoff]andJv[koff+1, ref_off]`.

Build the sparsity pattern for the rectangular CI Jacobian. Per-bus blocks have variable size (2 or 3). Off-diagonal blocks have entries only for the (e, f) columns of the neighbor (current injection from neighbor is independent of neighbor's Q variable). Slack cross-terms and LCC tail entries are added with structural zeros.

_jv_nz_index(Jv, row, col)

Return the nzval index for Jv[row, col]. Assumes the entry is structurally present (errors otherwise). Used at construction time to pre-compute indices for the hot-path update functions.

Populate the Ybus off-diagonal blocks (constant across NR iterations) and the REF row off-diagonal Ybus blocks. These entries are filled once and not touched during per-iteration updates.

Off-diagonal Ybus block: F = Ispec − Iinj, so the contribution to the Jacobian from `−Iinj = −Y·Vis the 2×2 real representation of−Y_ij`:

∂(−I_inj_r)/∂e_j = −G_ij,  ∂(−I_inj_r)/∂f_j =  B_ij
∂(−I_inj_i)/∂e_j = −B_ij,  ∂(−I_inj_i)/∂f_j = −G_ij

For PV neighbor columns, only the (e, f) sub-columns are populated; the Q column has structural zeros in off-diagonals (captured by pattern construction).

Write the LCC Jacobian entries. Mirrors polar _set_entries_for_lcc with these changes:

• Polar's idx_p_fb (Vm slot, single column) becomes two rectangular columns (col_e_fb, col_f_fb). Polar partials ∂(·)/∂Vm_fb translate via chain rule: ∂(·)/∂e = ∂(·)/∂Vm · e/|V|, similarly for f.
• The bus residual rows for fb are ΔI_r_fb and ΔI_i_fb (current mismatch), not polar's ΔP_fb / ΔQ_fb. The LCC contribution −Re(Y_lcc·V) = −A·u(ϕ) and −Im(Y_lcc·V) = −A·w(ϕ) (with the true ϕ from _calculate_ϕ_lcc, not the α-approximation) gives the bus-diagonal additions and tail cross-terms below. The inverter sign convention is absorbed into cos(ϕ_i), sin(ϕ_i) — no separate handling needed.

Tail row entries (∂Ft*) use the true-ϕ ∂P/∂V helper from _lcc_jacobian_scalars chain-ruled into (e, f) columns via ∂V/∂e = e/V, ∂V/∂f = f/V. Every chain term that picks up a 1/sin(ϕ) factor (i.e. cos(ϕ)·∂ϕ/∂x chain) goes through _dphi_dV_lcc / _dphi_dt_lcc / _dphi_dα_lcc, which return 0 at the sin(ϕ) → 0 clamp.

Update state-dependent Jacobian entries: per-bus diagonal blocks, slack cross-terms, LCC tail entries. Reads state from the residual's state caches (e_state, f_state, Q_state, P_eff_cache, Q_eff_cache) — these must be up to date (call the residual on x first).

All writes go through the pre-computed nzval-index caches into nonzeros(Jv), so the per-iteration cost scales as O(N + n_LCC) rather than incurring O(log(nnz_per_col)) per assignment.

_mixed_fill_state!(x, data, bus_state_offset, type_time_step, value_time_step)

Fill the MCPB state vector x: block layout from data.bus_type[:, type_time_step], values from *[:, value_time_step]. Equal time steps give the flat start; a converged value_time_step gives the previous-solution warm start while keeping the current step's offsets valid.

compute_mixed_bus_state_offsets(bus_type)

Per-bus state-vector offsets and block sizes for MCPB. Every bus occupies 2 slots (no PV→3 expansion). Returns (offsets, block_sizes, total_bus_state): offsets[i] is bus i's 1-based block start, offsets[end] the LCC-tail start (== total_bus_state + 1), block_sizes[i] == 2.

Distribute the converged subnetwork slack and write the bus power injections. Mirrors rect_finalize_bus_injections! / polar _setpq. Difference from rect: MCPB has no Q slot, so S_net,i = V_i·conj((Y_raw·V)_i) is recovered from a raw-Y matvec (raw Y excludes const-P/I/Z loads, avoiding the const-Z/const-I double-count). Q_net = imag(S_net); P_net is redistributed two-pass (sum total physical slack, then P_net = P_net_set_polar[i] + c_i·P_slack_total). Called once per step after convergence.

mixed_initial_state!(x, data, bus_state_offset, bus_block_size, time_step)

Initialize the MCPB state vector x (flat start) from data. REF slots hold (P_gen, Q_gen); all other buses hold (e, f). Counterpart of mixed_update_data!.

mixed_update_data!(data, x, bus_state_offset, bus_block_size, time_step)

Write state-derived voltages and LCC taps/angles from x back into data (per-iteration). Does NOT write power injections: at REF x[off] carries the whole subnetwork slack and at PV the gen Q is only known post-convergence, so both are finalized by mixed_finalize_bus_injections!. Counterpart of mixed_initial_state!.

struct ACMixedCPBResidual

Residual functor for the mixed current/power-balance (MCPB) AC power flow. Mirrors ACRectangularCIResidual 1:1, but every bus uses a 2-slot block (no PV→3 expansion).

Non-obvious fields: Y_bus_eff folds in ZIP constant-Z; P_net_const/ Q_net_const are the |V|-independent net injections; const_I_P/const_I_Q are the constant-current withdrawal coefficients; P_net_set is the initial Pnet for the distributed-slack delta; `validateoffsetsare the precomputed PQ/PVx`-offsets for the voltage-magnitude diagnostic. Remaining fields are named after their roles.

Update MCPB residual F (paper §IV). Mirrors _update_rect_ci_residual_values! step-for-step. Network current Y_bus_eff·V (+ LCC) is accumulated into per-bus Ir_acc/Ii_acc (not subtracted into F), keeping rect's residual = I_spec − I_network sign so the Jacobian mirrors rect. Per-bus blocks are 2 slots: PQ = divided-current balance, imag-first (so nonzero B_ii lands on the block diagonal); PV = eq.7 power balance then eq.8 |V|² − V_set²; REF = rect's REF branch verbatim.

struct ACMixedCPBJacobian

Jacobian functor for the mixed current/power-balance (MCPB) AC power flow. Mirrors ACRectangularCIJacobian 1:1, but every bus uses a 2-slot block (no PV→3 expansion): rect's pv_extra_nz is dropped, and offdiag_pv_nz is added for the PV power-balance row's off-diagonals, which are nonlinear in MCPB and rewritten each iteration. PQ off-diagonals are constant ±Y (_populate_mixed_constant_yb_blocks!). Per-iteration updates write into nonzeros(Jv) through nzval-index caches built once at construction, so the hot path is O(N + n_LCC). Field roles are in the inline comments below.

_build_mixed_diag_nz_cache(Jv, bus_state_offset)

diag_base_nz::Matrix{Int} (4 × n_buses): nzval indices of each per-bus 2×2 diagonal block, rows ordered (off,off), (off,off+1), (off+1,off), (off+1,off+1). No pv_extra_nz (MCPB PV blocks are 2×2, no Q column).

_build_offdiag_pv_nz_cache(Jv, Y_bus_eff, bus_state_offset, bus_types)

Cache nzval indices for the PV power-balance row's off-diagonal (e_k, f_k) columns. Each 2 × n_pv_pairs column is one ordered (PV bus i, neighbor k≠i) pair (REF neighbors excluded), rows Jv[i_off, k_off] and Jv[i_off, k_off+1]. The PV voltage-constraint row has no off-diagonals.

Build the MCPB Jacobian sparsity pattern (all blocks 2×2). Off-diagonal blocks reserve both (e, f) neighbor columns for PQ and PV rows; PQ entries are later filled constant by _populate_mixed_constant_yb_blocks!, PV entries are structural zeros written per-iteration. Slack cross-terms, LCC tail, and REF handling mirror rect verbatim.

Fill the constant Ybus off-diagonal blocks for PQ rows only. PV off-diagonals are nonlinear (left 0.0 here, written per-iteration via `offdiagpv_nz); REF rows have no(e, f)` off-diagonals.

Only −I_acc carries the off-diagonal dependence (the Ispec terms depend on bus i alone), and `∂Iracc/∂(ek,fk) = (Gik, −Bik),∂Iiacc/∂(ek,fk) = (Bik, G_ik)`. With the MCPB imag-first slot order this gives the 2×2 block

Jv[off,   k_off] = −B_ik   Jv[off,   k_off+1] = −G_ik
Jv[off+1, k_off] = −G_ik   Jv[off+1, k_off+1] = +B_ik

i.e. rect's off-diagonal block with its two rows swapped (rect is real-first; MCPB PQ is imag-first), matching the residual's PQ slot swap.

MCPB LCC tail. PQ/REF terminals match rect's _set_entries_for_lcc_rect!, with PQ using the imag-first slot order. At a PV terminal the LCC current enters via Iracc/Iiacc, so its contribution lands in the eq.7 power-balance row (F[off] += e_i·∂ΔIr_lcc + f_i·∂ΔIi_lcc); the eq.8 |V|² row gets none. Tail-row and tail×tail entries are identical to rect.

Update state-dependent MCPB Jacobian entries (per-bus diagonal blocks, PV off-diagonals, slack cross-terms, LCC tail) through the nzval-index caches. Reads the residual's shared state caches, so the residual must have been evaluated on the current x first.

MCPB PQ diagonal 2×2 (imag-first: slot off = imag balance, off+1 = real balance). Same divided-current + const-I expressions as rect's _update_pq_diag_block! with the two rows swapped — diag rows 1,2,3,4 are rect's rows 3,4,1,2 — matching the residual's PQ slot swap.

MCPB PV diagonal 2×2 (slot off = eq.7 power balance e·Ir + f·Ii − P, off+1 = eq.8 |V|² − V_set²). Ir/Ii include the Y_ii diagonal term; a const-I P_eff = P_net_const − cIP·Vm adds ∂P/∂(e,f) = −cIP·(e,f)/Vm.