Memoized AC-Jacobian sparse structure, stored in its OWN PowerFlowData field (not the shared solver_cache slot): the NR/TR AC Jacobian and a SolverCache can both be live in one solve — e.g. a FastDecoupled solve that hands off to NR uses a FastDecoupledCacheand this structure — so the two must not contend for a single slot. Cache key is the network-matrix identity + slack nonzero pattern (nzind); see _get_or_build_jacobian_structure.

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.

An abstract supertype for the persistent per-solve caches stored in PowerFlowData.solver_cache[]. Concrete subtypes (DCSolverCache for the DC/PTDF path, FastDecoupledCache for the polar fast-decoupled solver) are type-disjoint, so the slot's type discriminates which path populated it — no sentinel tag is needed and a cross-use is a plain MethodError rather than a silent reuse.

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 PowerSystems.Systemsys. 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 PowerSystems.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 PowerSystems.Systemsys. 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 PowerSystems.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 PowerSystems.Systemsys. 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 PowerSystems.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 PowerSystems.Systemsys. 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 PowerSystems.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 one time step's arc angle differences using precomputed from/to bus index vectors. Used by the AC solver where fb_ix/tb_ix are already available from the branch flow calculation.

get_bus_lookup(pfd::PowerFlowData)

Bus number → row index lookup for matrices stored in pfd (via the metadata PowerNetworkMatrices matrix). Use this when mapping device buses from a PowerSystems.System onto PowerFlowData injection and withdrawal arrays.

get_time_steps(pfd::PowerFlowData)

Number of time steps configured on the embedded PowerFlowEvaluationModel.

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 PowerSystems.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.

_ac_power_flow_with_area_relax!(data, pf, time_step; kwargs...) -> Bool

Wraps _ac_power_flow with greedy-relax handling for embedded area net-interchange control: on non-convergence with areas still enrolled, de-enroll the worst-|r_a| area and re-solve, warm-started with surviving areas' ΔP_a re-seeded from the delta_p mirror; repeat until convergence or exhaustion. Exhaustion while still failing is genuine network non-convergence plus a terminal diagnostic (_report_area_interchange_failure). Relaxation is never silent: an @error at each de-enrollment and a solve-end summary; converging after a relax still returns true.

Resets to the full pristine enrollment before each time step's attempt (_ensure_pristine_area_set!) so a previous step's relax never carries over. The never-enrolled short-circuit deliberately tests the PRISTINE set, not the WORKING one — a previous time step's relax may have emptied the working set, and short-circuiting on it would permanently disable area control for the rest of data's lifetime.

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), and calculates that time step's branch power flows. If the power flow does not converge, it sets the corresponding entries in data to NaN.

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 before that time step's branch flows are computed.

Examples

solve_power_flow!(data)

DCSolverCache{M, B, C, S} <: SolverCache

DC/PTDF persistent cache stored in data.solver_cache[]. The factored network matrix matrix and the backend form the invalidation key (rebuild when either changes); cache is the actual factorization and scratch the per-solve buffers. Subtypes SolverCache so it is type-disjoint from the AC FastDecoupledCache that shares the same slot (no sentinel tag needed).

_adjust_dc_slack_injections!(
    data::Union{PTDFPowerFlowData, vPTDFPowerFlowData, ABAPowerFlowData},
    power_injections::Matrix{Float64},
)

After a DC power-flow solve, the reference bus of each subnetwork is excluded from the linear system and its injection is left at the original generator setpoint. This function closes the active-power balance per subnetwork by subtracting the total subnetwork imbalance from the reference-bus injection. The imbalance is computed from bus_active_power_injections - bus_active_power_withdrawals + bus_hvdc_net_power over all buses in the subnetwork (including the reference bus). The optional power_injections matrix is updated in lockstep so that downstream optional post-processing, such as DC loss factors, uses the balanced injection vector.

_distribute_dc_slack!(data)

Distribute each subnetwork's active-power imbalance across its participating buses (including the reference bus) in proportion to bus_slack_participation_factors, writing the result into bus_active_power_injections. No-op unless the user configured distributed slack (explicit factors or headroom mode), detected via a non-empty computed_generator_slack_participation_factors.

_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.

_shift_angles_to_stored_reference!(data, time_steps = :)

The DC solve computes angles relative to 0 at each subnetwork's ref bus; if the subnetwork's ref bus has nonzero angle, shift the angles accordingly.

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.

1-based x-index offset for the area-interchange tail in the polar state/residual layout: ΔP_a for tail slot tail_ix lives at x[area_tail_offset(data, dcn) + tail_ix]; residual row r_a at the same index.

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

Length of the shared non-bus state tail (LCC + VSC + area interchange), identical across polar/rectangular/mixed formulations. The ONLY place the tail composition is defined — every state/residual sizing site calls this.

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

_initialize_residual_x0(pf::ACPolarPowerFlow, data, time_step; kwargs...)
    -> (residual, x0_computed)

Build the polar residual and the (warm-started, validated) initial state vector WITHOUT constructing the formulation Jacobian. Shared by initialize_power_flow_variables (which adds the Jacobian) and by the fast-decoupled :decoupled driver, whose B′/B″ half-steps never use the formulation Jacobian — that driver only materializes J when a handoff solver or loss/voltage- stability factors are requested, so building it eagerly here would waste a full sparse-Jacobian allocation + evaluation per solve (and per time step in multi-period runs).

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.

_d2P_lcc(V, t, α, I_dc, ϕ, σ) -> NamedTuple

Second partials of the LCC active-power contribution P_s = K·I·V·t·cos ϕ_s with respect to (V_s, t_s, α_s). σ = +1 rectifier, σ = -1 inverter; ϕ is the side-specific phase angle (its cos carries the side sign, as in _d2Q_lcc).

Interior regime: cos ϕ_r = cos α − β/(Vt) (inverter: −cos α − β/(Vt)), so P_s is linear in V·t with an α-dependent coefficient. Clamp regime (sin ϕ < LCC_sinϕ_TOLERANCE): cos ϕ is pinned at ±1, locally constant, so the mixed V–t partial K·I·cos ϕ is the only survivor and all α-partials vanish — matching the clamp-guarded first-derivative helpers.

_d2Q_lcc(V, t, α, x_t, I_dc, ϕ, σ) -> NamedTuple

Second partials of the LCC reactive-power contribution Q_s = V t K I sin ϕ_s with respect to (V_s, t_s, α_s). σ = +1 for the rectifier, σ = -1 for the inverter. I_dc > 0, x_t > 0, and ϕ is the side-specific phase angle (so cos ϕ_s already carries the side sign).

Uses the sin ϕ → 0 clamp guard: when sin ϕ falls below LCC_sinϕ_TOLERANCE, returns all-zero second partials, mirroring the existing _calculate_dQ_*_lcc helpers. In that regime Q_s ≈ 0 and the analytic formulas (which contain 1/sin³ϕ factors) are singular but the residual sees no Q_s dependence locally.

_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.

_fd_converter_substep!(data, time_step; max_iter=20, tol=1e-10) -> Float64

Sequential AC–DC converter sub-solve for the fast-decoupled (FDDecoupled) path: holding the AC bus voltages in data.bus_magnitude FIXED, Newton-iterate each LCC's four tail states [tap_r, tap_i, α_r, α_i] until its four control residuals (P-setpoint, DC-line balance, α limits; _write_lcc_tail!) fall below tol, refreshing data.lcc.branch_admittances via _update_ybus_lcc! so the converter's equivalent injection at the AC terminals is current. Each LCC is an independent 4×4 solve. Returns the final ‖tail residual‖∞ across all LCCs. This is the per-AC-iteration converter solve of the sequential AC–DC method: given V it produces the DC boundary conditions the AC half-steps then balance. The caller must copy the updated data.lcc.* states back into the state vector's trailing slots before the next residual evaluation (the residual reads those states from x).

_lcc_ac_active_powers(data, i, time_step, Vm_fb, Vm_tb) -> (P_lcc_from, P_lcc_to)

AC-side active power at LCC i's two converter terminals, each signed as power flowing from that AC bus INTO the DC link (rectifier > 0 draws, inverter < 0 receives), so the DC-line balance reads P_lcc_from + P_lcc_to = R·i_dc². Reads the tap/ϕ/idc state off data.lcc (refreshed from the iterate earlier in `updateresidualvalues!`). Shared by the LCC tail residual and the area-interchange DC-tie residual so both evaluate the exact same expression.

_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.

_lcc_tail_jacobian_block(data, i, time_step, Vm_fb, Vm_tb) -> Matrix{Float64}

The 4×4 Jacobian of LCC i's four tail residual rows (P-setpoint, DC-line balance, rectifier α-limit, inverter α-limit; see _write_lcc_tail!) with respect to its four tail states [tap_r, tap_i, α_r, α_i], at fixed AC terminal magnitudes Vm_fb/Vm_tb. Used by the sequential fast-decoupled LCC sub-solve (_fd_converter_substep!) to drive the converter control equations to zero given the AC voltages — it is the tail×tail diagonal block of the unified Jacobian. Requires data.lcc.*.phi to be current (call _update_ybus_lcc! first). The zero-DC-current case (i_dc == 0) returns the identity, matching the tap-pinning branch of _write_lcc_tail!.

_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|, θ)).

_write_lcc_state_to_x!(x, data, time_step)

Copy the per-LCC converter states (tap_r, tap_i, α_r, α_i) from data.lcc.* into the trailing 4·n_lcc slots of the state vector x — the inverse of the read the residual functor performs. Used by the sequential converter sub-solve so x stays consistent with the freshly-solved converter state before the next residual evaluation.

_write_lcc_tail!(F, data, base_offset, time_step, i, fb, tb, Vm_fb, Vm_tb)

Write LCC i's four tail residual rows (P-setpoint, DC-line balance, rectifier α-limit, inverter α-limit) into F.

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.

Internal lowering of all PSY DC components into a single DC network solved jointly with the AC buses. Fields ordered: converters (length n_conv), DC nodes (length n_node), DC branches (length n_branch). Built once from the System in initialize_DCNetwork! and stored on PowerFlowData behind a Ref (the struct is immutable; the network is assigned after the system is scanned, exactly like solver_cache).

initialize_DCNetwork!(data, sys, bus_lookup, reverse_bus_search_map, removed_buses)

Scan sys for DC components, lower them into a single DCNetwork, and store it on data. Handles point-to-point TwoTerminalVSCLine (each → 2 converters + 2 DC nodes + 1 DC branch). A system with no DC components leaves the empty DCNetwork() placeholder in place. The joint AC↔DC tail model is AC-only; for DC power flow the VSC stays a fixed-power injection (see lcc_vsc_fixed_injections!), so this is a no-op on non-AC data.

AreaInterchangeData

Holds the enrolled controlled areas and their ties for embedded net-interchange control. Always present on PowerFlowData — empty vectors when control is off, mirroring data.lcc with zero LCCs (no Union{Nothing} sentinel).

Fields

• areas::Vector{ControlledArea}, ties::Vector{AreaTie}: the WORKING (currently enrolled) set — what the residual/Jacobian tail actually sees. The greedy relax loop shrinks these mid-solve on a non-converged time step; see pristine_areas/pristine_ties below for the set this is derived from.
• dc_ties::Vector{DCTie}: DC-line (LCC/VSC) ties, enumerated alongside ties by build_area_ties and consumed by the residual/Jacobian/relax machinery exactly like ties. Empty on any system with no cross-area DC converter (zero behavior change for AC-only systems).
• tolerance::Float64: interchange convergence tolerance.
• ni_scratch::Vector{Float64}: pre-allocated per-area net-interchange accumulator, length length(areas), reused every residual evaluation by _set_area_tail_residuals! (see area_residual.jl) to avoid a per-iteration allocation on the hot path.
• delta_p::Matrix{Float64}: per-area, per-time-step ΔP_a mirror, sized (length(areas), n_time_steps). Unlike ni_scratch this is NOT reset every call — _update_residual_values! writes the last-evaluated ΔP_a into it on every residual evaluation, and update_state! reads it back to seed a warm re-solve's x0; the per-time-step column keeps a multi-period data's warm starts independent.
• pristine_areas::Vector{ControlledArea}, pristine_ties::Vector{AreaTie}, pristine_dc_ties::Vector{DCTie}: the FULL enrolled set as originally built by build_area_interchange_data, never mutated. Greedy relax shrinks only the WORKING set, for the rest of the CURRENT time step's attempts; _ensure_pristine_area_set! resets the working set from these before the next time step. Also the only source for a relaxed area's achieved net interchange in the results table, since its tail is translated to 0 out of the WORKING ties the moment it is de-enrolled.
• pristine_delta_p::Matrix{Float64}: persistent, PRISTINE-tail_ix-indexed mirror of delta_p — its row layout never changes under de-enrollment renumbering, so it survives across time steps. _ensure_pristine_area_set! reseeds the working delta_p from it; _sync_pristine_delta_p! writes back into it once a time step's solve converges.
• relaxed::Dict{Int, Vector{RelaxedAreaRecord}}: per-time-step record of areas relaxed away that time step (absent/empty = none), keyed by time_step so a multi-period data's relax decisions stay independent.

An AC branch whose endpoints straddle a controlled-area boundary (or touch one).

Fields

• from_bus_ix::Int, to_bus_ix::Int: reduced-network bus indices.
• nz_offsets::NTuple{4, Int}: Ybus nzval offsets Y11,Y12,Y21,Y22.
• metered_from::Bool: metered end is the from bus.
• from_area_tail::Int: tail_ix of the controlled area owning the from bus; 0 = uncontrolled.
• to_area_tail::Int: tail_ix of the controlled area owning the to bus; 0 = uncontrolled.
• diag_pollution::NTuple{2, ComplexF64}: (from, to) correction so the kernel can recover the corridor's OWN self-admittance from the aggregate Y-bus diagonal, which sums EVERY branch/shunt incident at that bus, not just this tie's members. Cached once at tie-build time as ybus_diag − Σ_{corridor members} y11/y22_primitive; live evaluation then does ybus_nzval[o[1]] − diag_pollution[1] (and the to-side analog). This keeps flows exact under a controlled tap ON the corridor (the tap's own delta cancels against the live diagonal read, so live−constant stays exact) while still being exact for any static topology; a non-member device mutating an endpoint diagonal is out of scope and fenced by an enrollment-time @warn (see enrollment.jl), not corrected here.

A PSY Area enrolled for embedded net-interchange control.

Fields

• name::String: PSY area name.
• slack_bus_ix::Int: reduced-network bus index of the area slack (PSS/E ISW).
• pdes::Float64: net interchange target, pu system base.
• tail_ix::Int: 1-based slot in the area tail.

A converter whose AC bus sits in one controlled area while its DC counterpart terminal sits in another — the DC-line analogue of AreaTie, carrying converter identity instead of a Y-bus block: a DC converter's AC-side active-power injection (from _lcc_ac_active_powers, or a VSC converter's P_c tail state) is what enters NI_a, not a branch flow read off Ybus.

Scope: PSY.TwoTerminalLCCLine and point-to-point PSY.TwoTerminalVSCLine only. PSY.InterconnectingConverter (multi-terminal DC) is NOT enumerated — pairing a converter with "the" other-area counterpart terminal is ambiguous on an N-terminal DC subnet (unlike a point-to-point line's unambiguous from/to) and needs its own union-find-style traversal.

Fields

• kind::DCTieKind: DC_TIE_LCC or DC_TIE_VSC.
• lcc_ix::Int: 1-based index into data.lcc.* (data.lcc.bus_indices[lcc_ix] == (from_bus_ix, to_bus_ix), rectifier="from"/inverter="to" always, independent of setpoint_at_rectifier); 0 when kind == DC_TIE_VSC.
• from_conv_ix::Int, to_conv_ix::Int: DCNetwork converter indices (dcn.p_c/q_c row) for the VSC line's from/to sides; both 0 when kind == DC_TIE_LCC.
• from_bus_ix::Int, to_bus_ix::Int: reduced-network AC bus indices of the two terminals.
• metered_from::Bool: metered end is the from terminal (mirrors AreaTie/_metered_from — same ext["metered_end"] convention, defaulting to from).
• from_area_tail::Int, to_area_tail::Int: owning-area tail indices; 0 = uncontrolled, exactly like AreaTie. A DC tie whose in-area terminal is that area's own SLACK bus is still emitted — NI_a and the bus's own power-balance row are independent residual rows, so there is no double-count to guard against.

Which converter model backs a DCTie. Int8-backed, mirroring VSCControlMode (vsc_parameters.jl): use-site code dispatches on it by equality (tie.kind == DC_TIE_LCC), never isa/<:, keeping DCTie a single concrete struct.

A controlled area that the greedy relax loop de-enrolled mid-solve because its interchange schedule proved unenforceable given the network/tie capacity at a failed Newton iterate. Recorded per time step on AreaInterchangeData.relaxed so the results table (post_processing.jl) can still report it after it drops out of the WORKING areas (pdes is otherwise unavailable once an area has been de-enrolled).

Fields

• name::String: PSY area name — matches a pristine_areas entry by name.
• pdes::Float64: the area's original net-interchange target, pu system base — unaffected by relaxation or by any tail renumbering a de-enrollment does to OTHER areas.

Length of the area-interchange state/residual tail: one ΔP_a per enrolled area.

Deduplicate tie candidates to one AreaTie per unordered reduced bus pair. _ybus_block_offsets already returns the AGGREGATE Y-bus block for a bus pair (parallel admittances are summed into one sparse entry — that IS the corridor flow), so pushing one AreaTie per parallel branch would double-count the corridor when downstream code sums per-tie flows. The first-seen candidate for a pair fixes the tie's direction and metered end; later parallel members are folded in (their own primitive diag accumulated into the corridor's diag_pollution sum), and a disagreement on the (direction-adjusted) metered end is a data inconsistency worth a @warn, not a silent pick.

diag_pollution is finalized only after every candidate for a key has been folded in: it is aggregate_ybus_diag − Σ_corridor_members primitive_diag, so the kernel can later recover the corridor's own self-admittance from the aggregate Y-bus diagonal (which sums EVERY branch/shunt incident at that bus, not just this corridor's members).

Rebuild the SAME filtered, ordered TwoTerminalLCCLine list initialize_LCCParameters! built, so index i here lines up with lcc.bus_indices[i] (rectifier="from"/inverter="to" always; see DCTie docstring). A reduction merging both LCC terminals into one bus shows up as fix == tix: the tie vanishes (no boundary left to cross).

Whether ext["metered_end"] marks the from end as metered. Absent key silently defaults to from-metered (the legitimate default); a present-but-unrecognized value (typo, wrong case, stray whitespace) warns and also defaults to from-metered rather than silently coercing it — a "to"-only special case would mask bad data.

Decompose a three-winding transformer into its three star-node windings — the same per-winding Arcs and availability gating PNM uses when stamping the Y-bus. The star bus is a real network node assigned the primary terminal's area, so each winding is a normal two-terminal tie candidate; KCL at the star bus makes the boundary-crossing windings sum to the transformer's true net export from each area.

Rebuild the SAME _available_vsc_lines list _lower_vsc_lines! iterated when lowering dcn, so converter pair (2i-1, 2i) here matches line i there (from-side, to-side; see DCNetwork/_lower_vsc_lines!). dcn.converter_ac_bus_ix is already the reduced-network AC bus index. iszero(n_vsc_converters(dcn)) with nonempty VSC lines means DC-network joint modeling was turned off (solver_settings[:model_dc_network] = false) — no P_c state exists to feed NI_a then, so no DC ties are enumerated (silent, mirrors "DC network off" already meaning "VSC ignored by the AC solve" elsewhere).

build_area_ties(sys, bus_lookup, ybus, nrd, bus_area_map, lcc, dcn) -> (Vector{AreaTie}, Vector{DCTie})

Enumerate AC-branch ties AND DC-line ties whose post-reduction endpoints straddle an area boundary.

bus_lookup maps PSY bus number to reduced-network index; ybus is the assembled AC_Ybus_Matrix; nrd is the PNM.NetworkReductionData used to resolve merged buses. bus_area_map maps a reduced-network bus index to its enrolled area tail_ix (0 for an uncontrolled or area-less bus); ties where both endpoints are 0 are not stored, since no residual row references that boundary. lcc/dcn are data.lcc/get_dc_network(data).

nz_offsets are cached via _ybus_block_offsets — offsets only, never admittance copies — so a controlled tap that is also a tie keeps feeding correct flows after apply_parameter! mutates the Y-bus in place. Out-of-service branches are excluded (PSY.get_available_components already filters on available; _tie_in_service filters the additional discrete-control CLOSED requirement). Parallel branches between the same reduced bus pair are deduplicated to one AreaTie (see _dedup_ties); DC ties are never paralleled (one TwoTerminalLCCLine/TwoTerminalVSCLine per DC-tie candidate) so need no analogous dedup pass.

build_dc_ties(sys, lcc, dcn, nrd, bus_area_map) -> Vector{DCTie}

Enumerate DC-line ties: PSY.TwoTerminalLCCLine and point-to-point PSY.TwoTerminalVSCLine whose two AC terminals resolve to different areas on the reduced network. PSY.InterconnectingConverter (multi-terminal DC) is not enumerated — see DCTie docstring. Empty on any system with no cross-area DC converter, including every AC-only system (zero overhead / no behavior change).

build_area_interchange_data(pf, sys, data) -> AreaInterchangeData

Derive the enrolled ControlledArea/AreaTie set for embedded PSS/E-style net-interchange control from sys, applying the enrollment guards in order. Called once, at PowerFlowData construction, gated on get_area_interchange_control(pf); see initialize_power_flow_data! for how the result populates data.area_interchange.

_area_net_interchange(ties, dc_ties, tail_ix, data, time_step) -> Float64

Net interchange for the area at tail_ix, computed directly from the tie-flow kernels at data's current bus state — independent of whether that area is still in data.area_interchange.areas. Used by the results table to report an achieved NI for a relaxed area (pass aid.pristine_ties/pristine_dc_ties with its PRISTINE tail_ix). ties and dc_ties must come from the SAME snapshot (both pristine or both working).

_area_residual_gaps(data, time_step) -> Vector{Float64}

r_a = NI_a - PDES_a for every CURRENTLY enrolled (working) controlled area — the SAME kernel _set_area_tail_residuals! uses inside the Newton residual, run against a scratch vector sized to just the area tail (offset 0) so it can be called standalone after a solve returns, converged or not.

_dc_tie_metered_active_power(data, dcn, tie, Vm, time_step) -> Float64

Active power flowing out of tie's METERED end over the DC link — the DC-line analogue of _tie_metered_active_power. For an LCC the metered terminal's _lcc_ac_active_powers value is already "power into the DC link" signed (rectifier from> 0, inverter to< 0), so it is returned as-is. For a VSC the converter's P_c state is bus-injection signed (into the AC bus), so power leaving the area is its negation.

_deenroll_area!(data, drop_tail_ix) -> ControlledArea

Remove the WORKING area at drop_tail_ix from data.area_interchange (greedy relax): the surviving areas' tail_ix are renumbered contiguously (the state vector has no room for gaps), ties/dc_ties referencing the dropped tail are translated to 0 (uncontrolled — the kernel already skips a zero tail), and delta_p is rebuilt at the new (smaller) size with each survivor's row carried over from its OLD tail_ix, for every time step column (so OTHER time steps' already-converged mirrors are not corrupted by a de-enrollment that only THIS time step's relax loop decided on).

Invalidates the Jacobian-structure/NR-factorization caches for the reason documented on _ensure_pristine_area_set!. The dropped ControlledArea's ORIGINAL pdes survives on it, unaffected by the renumbering.

_ensure_pristine_area_set!(data, time_step)

Reset the WORKING areas/ties/dc_ties/ni_scratch/delta_p to the full PRISTINE enrollment before a time step's greedy relax loop runs — relax decisions are per time step, never permanent.

Explicitly invalidates the Jacobian-structure and NR symbolic-factorization caches (data.ac_jacobian_structure_cache/data.polar_nr_cache): neither cache key sees a tail-size change made by mutating the SAME AreaInterchangeData object in place (the structure cache holds area_data by IDENTITY; the NR cache keys only on the network matrix, slack pattern, and backend type), so without this reset a wrong-sized cached structure/factorization would be silently reused. In-place mutation is the only option: PowerFlowData is immutable, so data.area_interchange itself can't be reassigned.

_set_area_tail_residuals!(F, x, data, area_off, time_step)

Write the area-interchange residual rows F[area_off + a] = NI_a - PDES_a for every enrolled controlled area. Sign convention: the metered side's controlled area (if any) gets +P_m, the other side's (if controlled) gets -P_m. x is unused: the polar formulation reads the current iterate off data.bus_magnitude/data.bus_angles, already updated in place earlier in _update_residual_values!.

_sync_pristine_delta_p!(data, time_step)

Write the WORKING delta_p[:, time_step] (one row per CURRENTLY enrolled area, keyed by its renumbered working tail_ix) back into pristine_delta_p[:, time_step] (keyed by the area's never-changing PRISTINE tail_ix, looked up by name), once a time step's solve converges. The persistent mirror is what the NEXT time step's _ensure_pristine_area_set! reseeds from, so a later time step's warm start recovers THIS time step's converged ΔP_a for every area that survived. A relaxed area's row is left untouched — it has no converged ΔP_a to record; the results table reports delta_p = 0.0 for it directly, not from this mirror.

_tie_admittances(tie, ybus_nzval) -> (g11, b11, g12, b12, g21, b21, g22, b22)

Read tie's corridor admittances from the aggregate Y-bus 2×2 block through tie.nz_offsets (never cached — a controlled tap mutates ybus_nzval in place, so this must be re-read every evaluation). The DIAGONAL reads (o[1]/o[4]) subtract the cached tie.diag_pollution correction to recover the corridor's own self-term — see AreaTie's docstring for the derivation. One shared read so the residual kernel and the Jacobian partials always differentiate identical values.

_tie_metered_active_power(tie, Vm_f, θ_f, Vm_t, θ_t, ybus_nzval) -> Float64

Active power flowing out of tie's METERED end. Admittances are read via _tie_admittances (diagonal-pollution-corrected — see that function's docstring).

Metered at from: P_m = Vf² g11 + Vf·Vt·(g12·cos(θf−θt) + b12·sin(θf−θt)). Metered at to swaps roles: P_m = Vt² g22 + Vt·Vf·(g21·cos(θt−θf) + b21·sin(θt−θf)).

_warn_area_violations(data, time_step)

Solve-end diagnostic: @warn for any CURRENTLY enrolled (working) area whose achieved NI still misses its target by more than data.area_interchange.tolerance after a converged solve. Defensive — the area residual row is driven to ~0 by the same Newton tolerance as every other row, so this should not normally fire.

Create the Jacobian matrix structure for the area-interchange tail (polar) — the bordered block. Column ΔPa (`areaoff + area.tailix): one structural entry at the area's slack-bus P-mismatch row (2*slackbusix - 1`), bus-type-invariant (survives a PV<->PQ Q-limit flip of that bus — the row position never depends on bus type). Row ra (area_off + area.tail_ix): the union pattern of both state columns at every non-REF tie endpoint bus of every tie incident to area a, via _push_area_row_endpoint_cols!. A tie with BOTH endpoints controlled (different areas) contributes to TWO area rows. ∂r_a/∂ΔP_a is structurally ABSENT (zero diagonal border — KLU's full pivoting handles it; do not stamp).

Fill the area-interchange tail Jacobian entries (polar). Column ΔPa is the constant -1.0 at each area's slack-bus P-mismatch row. Row ra is filled by a zero-then-accumulate two-pass sweep over every tie, since a boundary bus of degree > 1 feeds the same area row from multiple ties.

_tie_metered_active_power_partials(tie, Vm_f, θ_f, Vm_t, θ_t, ybus_nzval)
    -> (dPm_dVf, dPm_dθf, dPm_dVt, dPm_dθt)

∂P_m/∂(Vf,θf,Vt,θt) for tie's metered-end active power, differentiated term-by-term from _tie_metered_active_power's exact expression so both always agree — see _tie_admittances (area_residual.jl) for the shared, diagonal-pollution-corrected admittance read.

Sign convention: Δθ = θf − θt throughout (a single shared sincos), including the metered-at-to case — its natural angle argument θt − θf = −Δθ is folded into the g21/b21 terms below instead of needing a second sincos call.

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.
• 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.
• Jv::SparseArrays.SparseMatrixCSC{Float64, Int32}: The Jacobian matrix, which is updated by _update_jacobian_matrix_values!.
• 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.
• independent_ref::Set{Int}: Multi-swing REF bus indices, from _multi_swing_ref_indices. Computed once at construction because the Q-limit loop only flips PV↔PQ, never REF.

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

Update the Jacobian matrix Jv using _update_jacobian_matrix_values! 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

(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 _update_jacobian_matrix_values! 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.

Bus indices of REF buses sharing an island with another REF (multi-swing). Each self-balances its own P-slot (∂F_P/∂x[2i−1] = −1) instead of the distributed island scalar; single-swing islands are excluded and keep the distributed-slack path.

_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,
                              independent_ref, 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.

A multi-swing island (REF bus in independent_ref) holds each swing at its own fixed voltage and self-balances its own P-slot instead of sharing the island's distributed scalar (see ACRectangularCIResidual's REF branch): such a REF's P_gen = x[off] directly, bypassing c_k/P_slack_total entirely.

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
• 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}}
• independent_ref::Set{Int} — REF buses that share an island with another REF (multi-swing); precomputed once here (bus REF-status is fixed across a solve) so the hot per-iteration path never allocates a Set.
• 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
• 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, independent_ref)

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, refoff]. Islands keyed by a REF bus inindependentref` (multi-swing) are skipped entirely — those islands have no distributed-slack cross-terms in the structural pattern.

Pre-compute the nonzeros(Jv) indices for the VSC tail (layout-generic over rectangular CI and MCPB — both share _create_rect_ci_vsc_structure! and _set_entries_for_vsc_rect_mcpb!). The conv row order matches the slot push order in _create_rect_ci_vsc_structure!:

1-4   bus current-injection coupling: (off,pc) (off,qc) (off+1,pc) (off+1,qc)
5-6   control row r1:                 (pc,pc) (pc,vk)
7-9   control row r2:                 (qc,qc) (qc,off) (qc,off+1)
10-13 DC-KCL converter coupling:      (vk,pc) (vk,qc) (vk,off) (vk,off+1)

The (vk,vk) node diagonal is shared by every converter on a node, so it lives in node (set to G_dc[k,k] then accumulated) rather than per-converter.

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; `independentrefis the set of REF buses that share an island with another REF (multi-swing), precomputed once here (bus REF-status is fixed across a solve) so the hot per-iteration path never allocates aSet;validate_offsetsare 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.