From da322eced2caebf60a5f942f9527eb93552ad187 Mon Sep 17 00:00:00 2001 From: Cameron Khanpour <99142483+cameronkhanpour@users.noreply.github.com> Date: Mon, 1 Jun 2026 23:37:18 -0400 Subject: [PATCH 01/13] Handle negative net demand in DC load shedding --- docs/src/math/dc-opf.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/math/dc-opf.md b/docs/src/math/dc-opf.md index 084e484..43b40c3 100644 --- a/docs/src/math/dc-opf.md +++ b/docs/src/math/dc-opf.md @@ -78,7 +78,7 @@ For each parameter type, we compute ``\partial K / \partial p`` and then the ful ### Demand (``d``) -Demand enters the power balance and the load shed upper bound constraints through the clipped demand +Demand enters the power balance and the load shedding upper bound through ``d_+ = \max(d, 0)``: ```math @@ -88,7 +88,7 @@ Demand enters the power balance and the load shed upper bound constraints throug ``` For strictly positive demand, ``\partial d_+ / \partial d = 1``. For negative -demand, it is ``0``. At zero demand, the clipping function is non-smooth; +demand, it is ``0``. At zero demand, the positive part function is non-smooth; the implementation uses the fixed zero shedding convention already required by the collapsed bound ``0 \leq \text{psh} \leq 0``. From 93e43479c0791aaea5ac209e984cccc1d3a39db4 Mon Sep 17 00:00:00 2001 From: Cameron Khanpour <99142483+cameronkhanpour@users.noreply.github.com> Date: Tue, 2 Jun 2026 01:19:49 -0400 Subject: [PATCH 02/13] Support disconnected DC network islands --- docs/src/advanced.md | 6 +- docs/src/getting-started.md | 2 +- docs/src/math/dc-opf.md | 16 ++-- docs/src/math/dc-power-flow.md | 15 ++-- ext/PowerDiffAPFExt.jl | 7 +- src/PowerDiff.jl | 2 +- src/prob/dc_opf.jl | 3 +- src/prob/kkt_dc_opf.jl | 74 +++++++++------- src/sens/cost.jl | 8 +- src/sens/flowlimit.jl | 4 +- src/sens/lmp.jl | 13 +-- src/sens/susceptance.jl | 4 +- src/sens/topology.jl | 2 +- src/types/dc_network.jl | 93 ++++++++++++++++---- src/types/dc_opf_problem.jl | 10 ++- src/types/show.jl | 2 +- test/runtests.jl | 1 + test/test_apf_integration.jl | 5 ++ test/test_dc_islands.jl | 156 +++++++++++++++++++++++++++++++++ 19 files changed, 339 insertions(+), 84 deletions(-) create mode 100644 test/test_dc_islands.jl diff --git a/docs/src/advanced.md b/docs/src/advanced.md index 3545400..84ad36c 100644 --- a/docs/src/advanced.md +++ b/docs/src/advanced.md @@ -37,13 +37,17 @@ Stores the DC network topology and parameters. | `angmax`, `angmin` | `Vector{Float64}` | Phase angle difference limits | | `cq`, `cl` | `Vector{Float64}` | Cost coefficients (quadratic, linear) | | `c_shed` | `Vector{Float64}` | Load-shedding cost per bus | -| `ref_bus` | `Int` | Reference bus index (sequential) | +| `ref_bus` | `Int` | Preferred reference bus index (sequential) | | `tau` | `Float64` | Regularization parameter | | `id_map` | `IDMapping` | Bidirectional element ID mapping (original ↔ sequential) | Construct from a parsed MATPOWER network with `DCNetwork(parse_file("case14.m"))`, or with explicit parameters: `DCNetwork(n, m, k, A, G_inc, b; ...)`. +Use `reference_buses(net)` to obtain the effective reference set. It preserves +`ref_bus` for its energized island and deterministically adds one reference for +each additional island, including isolated buses. + ### ACNetwork Stores the AC network with vectorized admittance representation. diff --git a/docs/src/getting-started.md b/docs/src/getting-started.md index 8a85bf8..c44fffc 100644 --- a/docs/src/getting-started.md +++ b/docs/src/getting-started.md @@ -29,7 +29,7 @@ suggested `calc_sensitivity` commands to try. ## DC Power Flow -DC power flow computes voltage angles from the reduced system ``\theta_r = B_r^{-1} p_r``, where ``B_r`` is the susceptance-weighted Laplacian with the reference bus row and column deleted. +DC power flow computes voltage angles from the reduced system ``\theta_r = B_r^{-1} p_r``, where ``B_r`` is the susceptance-weighted Laplacian with one reference row and column deleted per energized island. ```julia dc_net = DCNetwork(net) diff --git a/docs/src/math/dc-opf.md b/docs/src/math/dc-opf.md index 43b40c3..40962e8 100644 --- a/docs/src/math/dc-opf.md +++ b/docs/src/math/dc-opf.md @@ -18,7 +18,7 @@ f &= W A \theta & (\nu_{\text{flow}}) \\ g_{\min} \leq g &\leq g_{\max} & (\rho_{\text{lb}}, \rho_{\text{ub}}) \\ 0 \leq \text{psh} &\leq d_+ & (\mu_{\text{lb}}, \mu_{\text{ub}}) \\ \alpha_{\min} \leq A\theta &\leq \alpha_{\max} & (\gamma_{\text{lb}}, \gamma_{\text{ub}}) \\ -\theta_{\text{ref}} &= 0 & (\eta_{\text{ref}}) +\theta_{\text{refs}} &= 0 & (\eta_{\text{ref}}) \end{aligned} ``` @@ -32,6 +32,7 @@ where: - ``c_{\text{shed}}`` is the load shedding cost vector - ``d_+ = \max(d, 0)`` is the curtailable portion of signed net demand; negative net demand remains in power balance as an injection - ``\tau`` is a small regularization parameter for numerical conditioning +- ``\text{refs}`` contains one deterministic reference bus per energized island, including isolated buses; the configured ``\text{ref_bus}`` remains the reference for its island ## KKT System for Implicit Differentiation @@ -51,13 +52,13 @@ The KKT variable vector ``z`` is structured as: z = [\theta, g, f, \text{psh}, \lambda_{\text{lb}}, \lambda_{\text{ub}}, \gamma_{\text{lb}}, \gamma_{\text{ub}}, \rho_{\text{lb}}, \rho_{\text{ub}}, \mu_{\text{lb}}, \mu_{\text{ub}}, \nu_{\text{bal}}, \nu_{\text{flow}}, \eta_{\text{ref}}] ``` -with total dimension ``5n + 6m + 3k + 1``. +with total dimension ``5n + 6m + 3k + n_{\text{ref}}``. ### KKT Conditions The KKT residual ``K(z, p)`` consists of: -1. **Stationarity w.r.t. ``\theta``**: ``B^\top \nu_{\text{bal}} + (WA)^\top \nu_{\text{flow}} + e_{\text{ref}} \eta_{\text{ref}} + A^\top (\gamma_{\text{ub}} - \gamma_{\text{lb}}) = 0`` +1. **Stationarity w.r.t. ``\theta``**: ``B^\top \nu_{\text{bal}} + (WA)^\top \nu_{\text{flow}} + E_{\text{ref}} \eta_{\text{ref}} + A^\top (\gamma_{\text{ub}} - \gamma_{\text{lb}}) = 0`` 2. **Stationarity w.r.t. ``g``**: ``2 C_q g + c_l - G_{\text{inc}}^\top \nu_{\text{bal}} - \rho_{\text{lb}} + \rho_{\text{ub}} = 0`` 3. **Stationarity w.r.t. ``f``**: ``\tau^2 f - \nu_{\text{flow}} - \lambda_{\text{lb}} + \lambda_{\text{ub}} = 0`` 4. **Stationarity w.r.t. psh**: ``c_{\text{shed}} - \nu_{\text{bal}} - \mu_{\text{lb}} + \mu_{\text{ub}} = 0`` @@ -66,7 +67,7 @@ The KKT residual ``K(z, p)`` consists of: 5c. **Complementary slackness (generation/shedding bounds)**: ``\rho \circ (\cdot) = 0``, ``\mu \circ (\cdot) = 0`` 6. **Primal feasibility**: ``G_{\text{inc}} g + \text{psh} - d - B\theta = 0`` 7. **Flow definition**: ``f - WA\theta = 0`` -8. **Reference bus**: ``\theta_{\text{ref}} = 0`` +8. **Reference buses**: ``\theta_{\text{refs}} = 0`` ### Analytical Sparse KKT Jacobian @@ -102,6 +103,9 @@ Switching affects the Laplacian ``B``, weight matrix ``W``, and flow definition ``` This propagates into the stationarity, power balance, and flow definition blocks of the KKT system. +Sensitivities are defined while the energized-island partition is fixed. Opening +or closing a bridge splits or merges islands and changes the reference set, so +the derivative is non-smooth at that topology boundary. ### Cost Coefficients (``c_q``, ``c_l``) @@ -132,7 +136,7 @@ Susceptances affect the same blocks as switching (through ``B`` and ``W``), but Locational marginal prices are the power balance duals ``\nu_{\text{bal}}``, decomposed as: ```math -\text{LMP} = \underbrace{\text{energy}}_{\text{uniform component}} + \underbrace{\text{congestion}}_{\text{flow-limit component}} +\text{LMP} = \underbrace{\text{energy}}_{\text{per-island uniform component}} + \underbrace{\text{congestion}}_{\text{flow-limit component}} ``` The congestion component is extracted by solving: @@ -141,4 +145,4 @@ The congestion component is extracted by solving: \text{congestion}[\text{non-ref}] = B_r^{-1} \left(A_r^\top W (\lambda_{\text{ub}}^{\text{std}} - \lambda_{\text{lb}}^{\text{std}}) + A_r^\top (\gamma_{\text{ub}}^{\text{std}} - \gamma_{\text{lb}}^{\text{std}})\right) ``` -where ``\lambda^{\text{std}}``, ``\gamma^{\text{std}}`` use the standard sign convention (non-negative for binding constraints). The ``\gamma`` terms capture congestion from binding phase angle difference limits. The energy component is uniform across all buses in a connected network and reflects the marginal cost of generation. +where ``\lambda^{\text{std}}``, ``\gamma^{\text{std}}`` use the standard sign convention (non-negative for binding constraints). The ``\gamma`` terms capture congestion from binding phase angle difference limits. The energy component is uniform within each energized island and reflects that island's marginal cost of generation. diff --git a/docs/src/math/dc-power-flow.md b/docs/src/math/dc-power-flow.md index 62c63ef..9cab6e6 100644 --- a/docs/src/math/dc-power-flow.md +++ b/docs/src/math/dc-power-flow.md @@ -4,16 +4,17 @@ For non-OPF power flow with fixed generation, the DC approximation linearizes the power flow equations using the susceptance-weighted Laplacian. The voltage -angles satisfy the reduced system obtained by eliminating the reference (slack) bus: +angles satisfy the reduced system obtained by eliminating one reference bus per +energized island: ```math \theta_r = B_r^{-1} \, p_r ``` where: -- ``B_r`` is the susceptance-weighted Laplacian with the reference bus row and column deleted (invertible for a connected network) -- ``p_r = g_r - d_r`` is the net injection vector with the reference entry removed -- ``\theta_{\text{ref}} = 0`` by convention +- ``B_r`` is the susceptance-weighted Laplacian with one reference row and column deleted per energized island, including isolated buses +- ``p_r = g_r - d_r`` is the net injection vector with the reference entries removed +- ``\theta_{\text{refs}} = 0`` by convention The susceptance-weighted Laplacian is: @@ -38,6 +39,10 @@ where the perturbation is a rank-1 update from the incidence column of branch `` \frac{\partial B_r}{\partial \mathrm{sw}_e} = -b_e \, a_{e,r} \, a_{e,r}^\top ``` +These derivatives apply while the energized-island partition is fixed. A +bridge opening or closing changes the reference set, so sensitivities are +non-smooth at the split or merge boundary. + ### Flow Sensitivity to Switching Branch flows are ``f = W A \theta`` where ``W = \operatorname{diag}(-b \circ \mathrm{sw})``. The flow sensitivity has both indirect (via angle changes) and direct (via the switching coefficient) contributions: @@ -54,7 +59,7 @@ Since ``p = g - d`` and generation is fixed, ``\partial p / \partial d = -I``. T \frac{\partial \theta}{\partial d} = -B_r^{-1} ``` -embedded in the non-reference block (with zero rows/columns for the reference bus). The flow sensitivity follows as: +embedded in the non-reference block (with zero rows/columns for the reference buses). The flow sensitivity follows as: ```math \frac{\partial f}{\partial d} = W A \frac{\partial \theta}{\partial d} diff --git a/ext/PowerDiffAPFExt.jl b/ext/PowerDiffAPFExt.jl index 77f6742..c837ba8 100644 --- a/ext/PowerDiffAPFExt.jl +++ b/ext/PowerDiffAPFExt.jl @@ -58,6 +58,7 @@ end Convert a `DCNetwork` to an `APF.Network`. APF networks lack generators, costs, and limits, so this is one-way. +APF exposes one slack bus, so disconnected `DCNetwork`s are rejected. Bus demand is set to zero (PD separates demand from network topology). Branch `status` is derived from switching state: `sw[e] > 0.5`. @@ -67,6 +68,10 @@ use `APF.from_power_models(pm_data)` directly instead. """ function PowerDiff.to_apf_network(net::PowerDiff.DCNetwork) n, m = net.n, net.m + refs = PowerDiff.reference_buses(net) + length(refs) == 1 || throw(ArgumentError( + "AcceleratedDCPowerFlows conversion requires one energized island; found $(length(refs))" + )) # All buses are active: DCNetwork is built from PM.build_ref() which filters # out inactive buses, so every bus in net.n is active by construction. @@ -89,7 +94,7 @@ function PowerDiff.to_apf_network(net::PowerDiff.DCNetwork) branches[e] = APF.Branch(e, net.sw[e] > 0.5, net.b[e], net.fmax[e], from_bus[e], to_bus[e]) end - return APF.Network("PowerDiff", buses, net.ref_bus, branches) + return APF.Network("PowerDiff", buses, only(refs), branches) end # ----------------------------------------------------------------------------- diff --git a/src/PowerDiff.jl b/src/PowerDiff.jl index 703d1aa..f8b7a89 100644 --- a/src/PowerDiff.jl +++ b/src/PowerDiff.jl @@ -111,7 +111,7 @@ export DCNetwork, DCPowerFlowState export DCOPFProblem, DCOPFSolution export DCSensitivityCache, invalidate! export solve!, update_demand!, update_fmax!, update_switching! -export calc_demand_vector, calc_susceptance_matrix +export calc_demand_vector, calc_susceptance_matrix, reference_buses # DC Sensitivity Functions (convenience wrappers) export calc_generation_participation_factors, calc_ptdf_from_sensitivity diff --git a/src/prob/dc_opf.jl b/src/prob/dc_opf.jl index 50d31c9..4abd08a 100644 --- a/src/prob/dc_opf.jl +++ b/src/prob/dc_opf.jl @@ -96,7 +96,8 @@ function solve!(prob::DCOPFProblem) μ_ub = -dual.(prob.cons.shed_ub) γ_lb = dual.(prob.cons.phase_diff_lb) γ_ub = -dual.(prob.cons.phase_diff_ub) - η_ref = dual(prob.cons.ref) + refs = reference_buses(prob.network) + η_ref = [dual(prob.cons.ref[i]) for i in refs] # Post-process phase angle difference duals for strict complementarity. # Interior point solvers leave gamma ≈ 1e-8 for non binding constraints. diff --git a/src/prob/kkt_dc_opf.jl b/src/prob/kkt_dc_opf.jl index 53f97b2..2da407a 100644 --- a/src/prob/kkt_dc_opf.jl +++ b/src/prob/kkt_dc_opf.jl @@ -281,19 +281,21 @@ The KKT system includes: - Primal: va (n), pg (k), f (m), psh (n) - Dual (inequality): lam_lb (m), lam_ub (m), gamma_lb (m), gamma_ub (m), rho_lb (k), rho_ub (k), mu_lb (n), mu_ub (n) - Dual (equality): nu_bal (n), nu_flow (m) -- Reference bus constraint: 1 +- Reference bus constraints: n_ref -Total: 5n + 6m + 3k + 1 +Total: 5n + 6m + 3k + n_ref """ kkt_dims(prob::DCOPFProblem) = kkt_dims(prob.network) -kkt_dims(n::Int, m::Int, k::Int) = 5n + 6m + 3k + 1 +kkt_dims(n::Int, m::Int, k::Int) = _dc_kkt_dims(n, m, k, 1) +_dc_kkt_dims(n::Int, m::Int, k::Int, n_ref::Int) = 5n + 6m + 3k + n_ref function kkt_dims(net::DCNetwork) n = getfield(net, :n) m = getfield(net, :m) k = getfield(net, :k) - # va(n) + pg(k) + f(m) + psh(n) + lam_lb(m) + lam_ub(m) + gamma_lb(m) + gamma_ub(m) + rho_lb(k) + rho_ub(k) + mu_lb(n) + mu_ub(n) + nu_bal(n) + nu_flow(m) + ref(1) - return kkt_dims(n, m, k) + n_ref = length(reference_buses(net)) + # va(n) + pg(k) + f(m) + psh(n) + lam_lb(m) + lam_ub(m) + gamma_lb(m) + gamma_ub(m) + rho_lb(k) + rho_ub(k) + mu_lb(n) + mu_ub(n) + nu_bal(n) + nu_flow(m) + ref(n_ref) + return _dc_kkt_dims(n, m, k, n_ref) end """ @@ -303,12 +305,14 @@ Compute all KKT variable indices from network dimensions. Single source of truth for index calculations. # Variable ordering -[va(n), pg(k), f(m), psh(n), lam_lb(m), lam_ub(m), gamma_lb(m), gamma_ub(m), rho_lb(k), rho_ub(k), mu_lb(n), mu_ub(n), nu_bal(n), nu_flow(m), eta(1)] +[va(n), pg(k), f(m), psh(n), lam_lb(m), lam_ub(m), gamma_lb(m), gamma_ub(m), rho_lb(k), rho_ub(k), mu_lb(n), mu_ub(n), nu_bal(n), nu_flow(m), eta(n_ref)] # Returns NamedTuple with index ranges for each variable block. """ -function kkt_indices(n::Int, m::Int, k::Int) +kkt_indices(n::Int, m::Int, k::Int) = _dc_kkt_indices(n, m, k, 1) + +function _dc_kkt_indices(n::Int, m::Int, k::Int, n_ref::Int) i = 0 idx_θ = (i+1):(i+n); i += n idx_g = (i+1):(i+k); i += k @@ -324,7 +328,7 @@ function kkt_indices(n::Int, m::Int, k::Int) idx_μ_ub = (i+1):(i+n); i += n idx_ν_bal = (i+1):(i+n); i += n idx_ν_flow = (i+1):(i+m); i += m - idx_η = i + 1 + idx_η = (i+1):(i+n_ref); i += n_ref return ( va = idx_θ, pg = idx_g, f = idx_f, psh = idx_psh, @@ -336,7 +340,7 @@ function kkt_indices(n::Int, m::Int, k::Int) ) end -kkt_indices(net::DCNetwork) = kkt_indices(net.n, net.m, net.k) +kkt_indices(net::DCNetwork) = _dc_kkt_indices(net.n, net.m, net.k, length(reference_buses(net))) kkt_indices(prob::DCOPFProblem) = kkt_indices(prob.network) # ============================================================================= @@ -351,8 +355,8 @@ Flatten solution primal and dual variables into a single vector for KKT evaluati # Variable ordering [va; pg; f; psh; lam_lb; lam_ub; gamma_lb; gamma_ub; rho_lb; rho_ub; mu_lb; mu_ub; nu_bal; nu_flow; eta] -where eta is the dual for the reference bus constraint `va[ref_bus] == 0`, packed -from the recovered `sol.eta_ref` (the KKT operator includes this entry). +where eta contains one dual per energized island reference constraint, packed +from the recovered `sol.eta_ref` vector. """ function flatten_variables(sol::DCOPFSolution, prob::DCOPFProblem) return vcat( @@ -370,7 +374,7 @@ function flatten_variables(sol::DCOPFSolution, prob::DCOPFProblem) sol.mu_ub, sol.nu_bal, sol.nu_flow, - [sol.eta_ref] + sol.eta_ref ) end @@ -421,12 +425,12 @@ s.t. G_inc * g + psh - d = B * θ (ν_bal) g ≤ gmax (ρ_ub) 0 ≤ psh (μ_lb) psh ≤ max(d, 0) (μ_ub) - θ[ref] = 0 (η_ref) + θ[refs] = 0 (η_ref) ``` # Returns Vector of KKT residuals (should be zero at optimum): -1. Stationarity w.r.t. θ: B' * ν_bal + (W*A)' * ν_flow + e_ref * η_ref + A'*(γ_ub - γ_lb) = 0 +1. Stationarity w.r.t. θ: B' * ν_bal + (W*A)' * ν_flow + E_ref * η_ref + A'*(γ_ub - γ_lb) = 0 2. Stationarity w.r.t. g: 2*Cq * g + cl - G_inc' * ν_bal - ρ_lb + ρ_ub = 0 3. Stationarity w.r.t. f: τ² * f - ν_flow - λ_lb + λ_ub = 0 4. Stationarity w.r.t. psh: c_shed - ν_bal - μ_lb + μ_ub = 0 @@ -460,13 +464,12 @@ function kkt(z::AbstractVector, net::DCNetwork, d::AbstractVector) B_mat = net.A' * W * net.A WA = W * net.A - # Reference bus indicator - e_ref = zeros(n) - e_ref[net.ref_bus] = 1.0 + refs = reference_buses(net) # KKT conditions # 1. Stationarity w.r.t. θ - K_θ = B_mat' * ν_bal + WA' * ν_flow + e_ref * η_ref + net.A' * (net.sw .* (γ_ub - γ_lb)) + K_θ = B_mat' * ν_bal + WA' * ν_flow + net.A' * (net.sw .* (γ_ub - γ_lb)) + K_θ[refs] .+= η_ref # 2. Stationarity w.r.t. g K_g = 2 * Diagonal(net.cq) * g + net.cl - net.G_inc' * ν_bal - ρ_lb + ρ_ub @@ -512,10 +515,10 @@ function kkt(z::AbstractVector, net::DCNetwork, d::AbstractVector) # 10. Primal feasibility: flow definition K_flow_def = f - WA * θ - # 11. Reference bus - K_ref = θ[net.ref_bus] + # 11. Reference buses + K_ref = θ[refs] - return vcat(K_θ, K_g, K_f, K_psh, K_λ_lb, K_λ_ub, K_γ_lb, K_γ_ub, K_ρ_lb, K_ρ_ub, K_μ_lb, K_μ_ub, K_power_bal, K_flow_def, [K_ref]) + return vcat(K_θ, K_g, K_f, K_psh, K_λ_lb, K_λ_ub, K_γ_lb, K_γ_ub, K_ρ_lb, K_ρ_ub, K_μ_lb, K_μ_ub, K_power_bal, K_flow_def, K_ref) end # ============================================================================= @@ -583,7 +586,9 @@ function calc_kkt_jacobian(net::DCNetwork, d::AbstractVector, prob::DCOPFProblem n = getfield(net, :n) m = getfield(net, :m) k = getfield(net, :k) - dim = kkt_dims(n, m, k) + refs = reference_buses(net) + n_ref = length(refs) + dim = _dc_kkt_dims(n, m, k, n_ref) A = getfield(net, :A) G_inc = getfield(net, :G_inc) b = getfield(net, :b) @@ -594,7 +599,6 @@ function calc_kkt_jacobian(net::DCNetwork, d::AbstractVector, prob::DCOPFProblem cq = getfield(net, :cq) angmin = getfield(net, :angmin) angmax = getfield(net, :angmax) - ref_bus = getfield(net, :ref_bus) tau = getfield(net, :tau) va = getfield(sol, :va) pg = getfield(sol, :pg) @@ -614,16 +618,18 @@ function calc_kkt_jacobian(net::DCNetwork, d::AbstractVector, prob::DCOPFProblem B_mat = sparse(A' * W * A) # Build Jacobian blocks using centralized index calculation - idx = kkt_indices(n, m, k) + idx = _dc_kkt_indices(n, m, k, n_ref) A_rowptr, A_rowcols, A_rowvals = _sparse_row_storage(A) G_rowptr, G_rowcols, G_rowvals = _sparse_row_storage(G_inc) rowval = Int[] nzval = Float64[] - entry_hint = 4 * nnz(A) + 2 * nnz(B_mat) + 2 * nnz(G_inc) + 5n + 10m + 7k + 1 + entry_hint = 4 * nnz(A) + 2 * nnz(B_mat) + 2 * nnz(G_inc) + 5n + 10m + 7k + 2n_ref sizehint!(rowval, entry_hint) sizehint!(nzval, entry_hint) colptr = Vector{Int}(undef, dim + 1) Aθ = A * va + ref_pos = zeros(Int, n) + ref_pos[refs] .= eachindex(refs) @inline function start_col!(col::Int) colptr[col] = length(rowval) + 1 @@ -659,7 +665,7 @@ function calc_kkt_jacobian(net::DCNetwork, d::AbstractVector, prob::DCOPFProblem e = rowvals(A)[p] _push_csc_entry!(rowval, nzval, idx.nu_flow[e], b[e] * sw[e] * nonzeros(A)[p]) end - j == ref_bus && _push_csc_entry!(rowval, nzval, idx.η, 1.0) + ref_pos[j] > 0 && _push_csc_entry!(rowval, nzval, idx.η[ref_pos[j]], 1.0) end # pg columns: @@ -803,10 +809,12 @@ function calc_kkt_jacobian(net::DCNetwork, d::AbstractVector, prob::DCOPFProblem _push_csc_entry!(rowval, nzval, idx.f[e], -1.0) end - # eta column: - # ∂K_θ/∂η = e_ref - start_col!(idx.η) - _push_csc_entry!(rowval, nzval, idx.va[ref_bus], 1.0) + # eta columns: + # ∂K_θ/∂η = E_ref + @inbounds for (i, ref_bus) in enumerate(refs) + start_col!(idx.η[i]) + _push_csc_entry!(rowval, nzval, idx.va[ref_bus], 1.0) + end colptr[dim + 1] = length(rowval) + 1 return SparseMatrixCSC(dim, dim, colptr, rowval, nzval) @@ -824,8 +832,8 @@ function calc_kkt_jacobian_demand(net::DCNetwork, d::AbstractVector, sol::DCOPFS n = getfield(net, :n) m = getfield(net, :m) k = getfield(net, :k) - dim = kkt_dims(n, m, k) - idx = kkt_indices(n, m, k) + dim = kkt_dims(net) + idx = kkt_indices(net) mu_ub = getfield(sol, :mu_ub) colptr = Vector{Int}(undef, n + 1) @@ -900,7 +908,7 @@ function calc_kkt_jacobian_switching(prob::DCOPFProblem, sol::DCOPFSolution) J_s = spzeros(dim, m) # Use centralized index calculation - idx = kkt_indices(n, m, k) + idx = kkt_indices(net) # Precompute A * θ once (invariant across branches) Aθ = A * θ diff --git a/src/sens/cost.jl b/src/sens/cost.jl index 571d295..980c8a6 100644 --- a/src/sens/cost.jl +++ b/src/sens/cost.jl @@ -37,8 +37,8 @@ function calc_kkt_jacobian_cost_linear(net::DCNetwork) n = getfield(net, :n) m = getfield(net, :m) k = getfield(net, :k) - dim = kkt_dims(n, m, k) - idx = kkt_indices(n, m, k) + dim = kkt_dims(net) + idx = kkt_indices(net) colptr = Vector{Int}(undef, k + 1) rowval = Int[] @@ -93,8 +93,8 @@ function calc_kkt_jacobian_cost_quadratic(prob::DCOPFProblem, sol::DCOPFSolution n = getfield(net, :n) m = getfield(net, :m) k = getfield(net, :k) - dim = kkt_dims(n, m, k) - idx = kkt_indices(n, m, k) + dim = kkt_dims(net) + idx = kkt_indices(net) g = getfield(sol, :pg) diff --git a/src/sens/flowlimit.jl b/src/sens/flowlimit.jl index 87438c2..fe0571a 100644 --- a/src/sens/flowlimit.jl +++ b/src/sens/flowlimit.jl @@ -53,8 +53,8 @@ function calc_kkt_jacobian_flowlimit(prob::DCOPFProblem, sol::DCOPFSolution) n = getfield(net, :n) m = getfield(net, :m) k = getfield(net, :k) - dim = kkt_dims(n, m, k) - idx = kkt_indices(n, m, k) + dim = kkt_dims(net) + idx = kkt_indices(net) lambda_lb = getfield(sol, :lam_lb) lambda_ub = getfield(sol, :lam_ub) diff --git a/src/sens/lmp.jl b/src/sens/lmp.jl index 3162e63..c558023 100644 --- a/src/sens/lmp.jl +++ b/src/sens/lmp.jl @@ -26,7 +26,7 @@ # LMP = ν_bal = energy_component + congestion_component # where: # congestion_component = B_r⁻¹ [A_r' Diag(-b .* sw) (λ_ub - λ_lb) + A_r'(γ_ub - γ_lb)] (non-ref block) -# energy_component = ν_bal - congestion_component (uniform for connected network) +# energy_component = ν_bal - congestion_component (uniform within each island) # # Sign conventions (DC OPF): # - Our LMPs are positive (cost increases when demand increases) @@ -77,13 +77,14 @@ end Extract the congestion component of LMPs for analysis. From the θ-stationarity KKT condition: - B' * ν_bal + (WA)' * ν_flow + e_ref * η_ref + A'*(γ_ub - γ_lb) = 0 + B' * ν_bal + (WA)' * ν_flow + E_ref * η_ref + A'*(γ_ub - γ_lb) = 0 The congestion RHS includes both flow limit duals and angle difference duals: congestion[non_ref] = B_r \\ (A' W (λ_ub - λ_lb) + A'(γ_ub - γ_lb))[non_ref] The congestion component captures price differentiation due to binding flow and angle -constraints, with the reference bus congestion component equal to zero. +constraints, with each energized-island reference bus congestion component equal +to zero. # Returns Vector (length n) of congestion contributions to each bus's LMP. @@ -91,7 +92,7 @@ Vector (length n) of congestion contributions to each bus's LMP. function calc_congestion_component(sol::DCOPFSolution, net::DCNetwork; B_r_factor=sol.B_r_factor) w = -net.b # positive weights (b < 0 for inductive lines) - non_ref = setdiff(1:net.n, net.ref_bus) + non_ref = _non_reference_buses(net) At = net.A' rhs_full = At * Diagonal(w .* net.sw) * (sol.lam_ub - sol.lam_lb) + At * (sol.gamma_ub - sol.gamma_lb) @@ -106,8 +107,8 @@ end Extract the energy (non-congestion) component of LMPs for analysis. -This is the uniform price component: energy = ν_bal - congestion -For a connected network, this should be approximately constant across all buses. +This is the uniform price component: energy = ν_bal - congestion. +It is approximately constant within each energized island. # Returns Vector (length n) of energy contributions to each bus's LMP. diff --git a/src/sens/susceptance.jl b/src/sens/susceptance.jl index 37d57c1..4729e10 100644 --- a/src/sens/susceptance.jl +++ b/src/sens/susceptance.jl @@ -38,7 +38,7 @@ Susceptance b affects: - Flow definition: f = W * A * theta The affected KKT conditions are: -- K_theta = B' * nu_bal + (WA)' * nu_flow + e_ref * eta_ref + A'*(gamma_ub - gamma_lb) +- K_theta = B' * nu_bal + (WA)' * nu_flow + E_ref * eta_ref + A'*(gamma_ub - gamma_lb) (gamma term has no b-dependence, so ∂K_theta/∂b only comes from B and WA terms) - K_power_bal = G_inc * g + psh - d - B * theta - K_flow_def = f - W * A * theta @@ -51,7 +51,7 @@ function calc_kkt_jacobian_susceptance(prob::DCOPFProblem, sol::DCOPFSolution) net = prob.network n, m, k = net.n, net.m, net.k dim = kkt_dims(net) - idx = kkt_indices(n, m, k) + idx = kkt_indices(net) theta = sol.va nu_bal = sol.nu_bal diff --git a/src/sens/topology.jl b/src/sens/topology.jl index 89d458a..04cb5c8 100644 --- a/src/sens/topology.jl +++ b/src/sens/topology.jl @@ -27,7 +27,7 @@ For DC power flow `θ_r = B_r⁻¹ p_r`, the sensitivity of angles w.r.t. switch where `∂B_r/∂swₑ = -bₑ · a_{e,r} · a_{e,r}'` is a rank-1 update from the incidence column of branch `e` restricted to non-reference buses, and `B_r` is the susceptance -matrix with the reference bus row and column deleted. +matrix with one reference row and column deleted per energized island. # Arguments - `state`: DCPowerFlowState containing the solved power flow diff --git a/src/types/dc_network.jl b/src/types/dc_network.jl index ed98782..2d32ff9 100644 --- a/src/types/dc_network.jl +++ b/src/types/dc_network.jl @@ -36,7 +36,7 @@ topology sensitivity analysis. - `angmax`, `angmin`: Phase angle difference limits - `cq`, `cl`: Quadratic and linear generation cost coefficients - `c_shed`: Load-shedding cost per bus (penalty for involuntary load curtailment) -- `ref_bus`: Reference bus index (phase angle = 0) +- `ref_bus`: Preferred reference bus index (phase angle = 0) - `tau`: Regularization parameter for strong convexity - `id_map`: Bidirectional mapping between original and sequential element IDs - `demand`: Real power demand aggregated per bus @@ -85,7 +85,7 @@ Solution container for DC OPF problem, storing both primal and dual variables. - `rho_ub`, `rho_lb`: Generator upper/lower bound duals - `mu_lb`, `mu_ub`: Load shedding lower/upper bound duals - `gamma_lb`, `gamma_ub`: Phase angle difference lower/upper bound duals -- `eta_ref`: Reference bus constraint dual (`va[ref_bus] == 0`) +- `eta_ref`: Reference bus constraint duals (`va[reference_buses(net)] == 0`) - `objective`: Optimal objective value - `B_r_factor`: Cached factorization of reduced susceptance matrix `B[non_ref, non_ref]` """ @@ -104,7 +104,7 @@ struct DCOPFSolution{F<:Factorization{Float64}} <: AbstractOPFSolution mu_ub::Vector{Float64} gamma_lb::Vector{Float64} gamma_ub::Vector{Float64} - eta_ref::Float64 + eta_ref::Vector{Float64} objective::Float64 B_r_factor::F end @@ -116,19 +116,19 @@ DC power flow solution (phase angles from reduced-Laplacian solve, no optimizati Supports both generation and demand for flexible sensitivity analysis. Unlike DCOPFSolution, this represents a simple power flow solution -`θ_r = B_r \\ p_r` where `B_r` is the susceptance matrix with the reference bus row and -column deleted (invertible for a connected network), without optimal dispatch or +`θ_r = B_r \\ p_r` where `B_r` is the susceptance matrix with one reference row and +column deleted per energized island, without optimal dispatch or constraint handling. # Fields - `net`: DCNetwork data -- `va`: Phase angles (rad), with `va[ref_bus] = 0` +- `va`: Phase angles (rad), with `va[reference_buses(net)] = 0` - `p`: Net injection vector (p = pg - d) - `pg`: Generation vector - `d`: Demand vector - `f`: Branch flows (computed from va) - `B_r_factor`: Factorization of `B[non_ref, non_ref]` (Cholesky for inductive networks, LU fallback) -- `non_ref`: Indices of non-reference buses +- `non_ref`: Indices excluding one reference bus per energized island """ struct DCPowerFlowState{F<:Factorization{Float64}} <: AbstractPowerFlowState net::DCNetwork @@ -624,24 +624,87 @@ function calc_susceptance_matrix(network::DCNetwork) return sparse(network.A' * W * network.A) end +""" + reference_buses(net::DCNetwork) → Vector{Int} + +Return one deterministic reference bus for each energized island. + +The configured `net.ref_bus` is preserved for its island. Every other island, +including an isolated bus, uses its lowest sequential bus index. A branch is +energized when `b[e] * sw[e] != 0`. +""" +function reference_buses(net::DCNetwork) + parent = collect(1:net.n) + + function find_root(i::Int) + while parent[i] != i + parent[i] = parent[parent[i]] + i = parent[i] + end + return i + end + + function union_roots!(i::Int, j::Int) + root_i = find_root(i) + root_j = find_root(j) + root_i == root_j && return nothing + if root_i < root_j + parent[root_j] = root_i + else + parent[root_i] = root_j + end + return nothing + end + + from_bus = zeros(Int, net.m) + to_bus = zeros(Int, net.m) + rows, cols, nz_values = findnz(net.A) + @inbounds for p in eachindex(rows) + if nz_values[p] > 0 + from_bus[rows[p]] = cols[p] + else + to_bus[rows[p]] = cols[p] + end + end + + @inbounds for e in 1:net.m + iszero(net.b[e] * net.sw[e]) && continue + from_bus[e] > 0 && to_bus[e] > 0 || + error("Incidence matrix row $e must have exactly two nonzero entries") + union_roots!(from_bus[e], to_bus[e]) + end + + refs_by_root = Dict{Int,Int}() + @inbounds for bus in 1:net.n + root = find_root(bus) + refs_by_root[root] = min(get(refs_by_root, root, bus), bus) + end + refs_by_root[find_root(net.ref_bus)] = net.ref_bus + return sort!(collect(values(refs_by_root))) +end + +"""Return bus indices after removing one reference bus per energized island.""" +_non_reference_buses(net::DCNetwork) = setdiff(1:net.n, reference_buses(net)) + """ _factorize_B_r(net::DCNetwork) → (factor, non_ref) Factorize the reduced susceptance matrix `B[non_ref, non_ref]`. Uses Cholesky for standard inductive networks (~2x faster), with LU fallback -for edge cases (capacitive branches or disconnected networks) where B_r is not -positive definite. Follows the approach of AcceleratedDCPowerFlows.jl. +for edge cases such as capacitive branches where B_r is not positive definite. +One reference row and column is removed per energized island, so disconnected +networks and isolated buses remain well-defined. """ function _factorize_B_r(net::DCNetwork) B = calc_susceptance_matrix(net) - non_ref = setdiff(1:net.n, net.ref_bus) + non_ref = _non_reference_buses(net) B_r = B[non_ref, non_ref] factor = try cholesky(Symmetric(B_r)) catch e e isa PosDefException || rethrow() - _SILENCE_WARNINGS[] || @warn "Reduced susceptance matrix B_r is not positive definite (e.g., capacitive branches or disconnected subnetwork); falling back to LU factorization. Results remain correct." + _SILENCE_WARNINGS[] || @warn "Reduced susceptance matrix B_r is not positive definite (e.g., capacitive branches); falling back to LU factorization. Results remain correct." lu(B_r) end return factor, non_ref @@ -670,9 +733,9 @@ Solve DC power flow for given generation and demand. Computes phase angles θ by solving the reduced system: B_r * θ_r = p_r -where B_r is the susceptance-weighted Laplacian with the reference bus row and -column deleted (invertible for a connected network), and p_r is the net injection -with the reference entry removed. The reference bus angle is zero by construction. +where B_r is the susceptance-weighted Laplacian with one reference bus row and +column deleted per energized island, and p_r is the net injection with those +reference entries removed. Reference bus angles are zero by construction. # Arguments - `net`: DCNetwork containing topology and parameters @@ -701,7 +764,7 @@ function DCPowerFlowState(net::DCNetwork, g::AbstractVector{<:Real}, d::Abstract # Factorize reduced susceptance matrix (Cholesky with LU fallback) F, non_ref = _factorize_B_r(net) - # Solve reduced system: θ[non_ref] = B_r \ p[non_ref], θ[ref] = 0 + # Solve reduced system: θ[non_ref] = B_r \ p[non_ref], θ[refs] = 0 θ = zeros(n) θ[non_ref] = F \ p[non_ref] diff --git a/src/types/dc_opf_problem.jl b/src/types/dc_opf_problem.jl index e35ee87..30873f3 100644 --- a/src/types/dc_opf_problem.jl +++ b/src/types/dc_opf_problem.jl @@ -51,7 +51,7 @@ precomputed at construction time). - `dz_dfmax`: Full KKT derivative w.r.t. flow limits (or nothing) - `dz_db`: Full KKT derivative w.r.t. susceptances (or nothing) - `b_r_factor`: Cached reduced susceptance factorization (topology-dependent, survives demand changes) -- `work`: Scratch workspace for VJP/JVP KKT solves (lazily allocated on first call, survives invalidation since its size depends only on `(n, m, k)` which are fixed at construction) +- `work`: Scratch workspace for VJP/JVP KKT solves (lazily allocated on first call) """ mutable struct DCSensitivityCache solution::Union{Nothing,DCOPFSolution} @@ -97,12 +97,13 @@ end """ invalidate_topology!(cache::DCSensitivityCache) -Clear all cached data including topology-dependent `b_r_factor`. +Clear all cached data including topology-dependent `b_r_factor` and `work`. Called when network topology changes (switching, susceptances). """ function invalidate_topology!(cache::DCSensitivityCache) invalidate!(cache) cache.b_r_factor = nothing + cache.work = nothing return nothing end @@ -244,8 +245,9 @@ function _rebuild_jump_model!(prob::DCOPFProblem) shed_lb = @constraint(model, psh .>= 0) shed_ub = @constraint(model, psh .<= _shed_capacity.(d)) - # Reference bus - ref_con = @constraint(model, va[network.ref_bus] == 0.0) + # One angle reference per energized island removes the Laplacian nullspace. + refs = reference_buses(network) + ref_con = @constraint(model, [i in refs], va[i] == 0.0) # Open lines should not constrain angle differences. phase_diff_lb = @constraint(model, network.sw .* (network.A * va) .>= network.sw .* network.angmin) diff --git a/src/types/show.jl b/src/types/show.jl index ad1bd71..ccb7d98 100644 --- a/src/types/show.jl +++ b/src/types/show.jl @@ -76,7 +76,7 @@ end function Base.show(io::IO, ::MIME"text/plain", net::DCNetwork) println(io, "DCNetwork ($(net.n) buses, $(net.m) branches, $(net.k) gens)") - println(io, " Reference bus: $(net.ref_bus)") + println(io, " Reference buses: $(reference_buses(net))") n_open = count(x -> x < 1.0, net.sw) println(io, " Open branches: $n_open/$(net.m)") println(io, " Flow limits: [$(round(minimum(net.fmax); digits=2)), $(round(maximum(net.fmax); digits=2))]") diff --git a/test/runtests.jl b/test/runtests.jl index 4040409..bc3df27 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -863,6 +863,7 @@ include("test_ac_pf_verification.jl") include("test_update_switching.jl") include("test_update_fmax.jl") include("test_psh.jl") +include("test_dc_islands.jl") include("test_nonbasic.jl") include("test_jvp_vjp.jl") include("test_kkt_vjp_jvp.jl") diff --git a/test/test_apf_integration.jl b/test/test_apf_integration.jl index d00c75d..4631b98 100644 --- a/test/test_apf_integration.jl +++ b/test/test_apf_integration.jl @@ -171,6 +171,11 @@ APF = AcceleratedDCPowerFlows end end +@testset "to_apf_network rejects multiple islands" begin + net = DCNetwork(2, 0, 0, spzeros(0, 2), spzeros(2, 0), Float64[]) + @test_throws ArgumentError to_apf_network(net) +end + # ========================================================================= # Index alignment # ========================================================================= diff --git a/test/test_dc_islands.jl b/test/test_dc_islands.jl new file mode 100644 index 0000000..3506443 --- /dev/null +++ b/test/test_dc_islands.jl @@ -0,0 +1,156 @@ +# Copyright 2026 Samuel Talkington and contributors +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 + +function _make_bridge_network() + A = sparse([ + 1.0 -1.0 0.0 0.0 + 0.0 1.0 -1.0 0.0 + 0.0 0.0 1.0 -1.0 + ]) + G_inc = sparse([ + 1.0 0.0 + 0.0 0.0 + 0.0 1.0 + 0.0 0.0 + ]) + return DCNetwork(4, 3, 2, A, G_inc, fill(-10.0, 3); + fmax=fill(10.0, 3), gmax=fill(10.0, 2), gmin=zeros(2), + cq=fill(1.0, 2), cl=[10.0, 20.0], ref_bus=1, tau=1e-3) +end + +function _make_isolated_load_network() + A = sparse([1.0 -1.0 0.0]) + G_inc = sparse(reshape([1.0, 0.0, 0.0], 3, 1)) + return DCNetwork(3, 1, 1, A, G_inc, [-10.0]; + fmax=[10.0], gmax=[10.0], gmin=[0.0], + cq=[1.0], cl=[10.0], ref_bus=1, tau=1e-3) +end + +function _make_fully_isolated_network() + return DCNetwork(2, 0, 2, spzeros(0, 2), sparse(Matrix(1.0I, 2, 2)), Float64[]; + fmax=Float64[], gmax=fill(10.0, 2), gmin=zeros(2), + cq=fill(1.0, 2), cl=[10.0, 20.0], ref_bus=1, tau=1e-3) +end + +@testset "Multi-island DC support" begin + @testset "isolated load shedding and LMP decomposition" begin + net = _make_isolated_load_network() + d = [0.0, 0.4, 0.3] + @test reference_buses(net) == [1, 3] + + state = DCPowerFlowState(net, [0.4, 0.0, 0.0], d) + @test state.non_ref == [2] + @test state.va[reference_buses(net)] == zeros(2) + @test all(isfinite, state.f) + + prob = DCOPFProblem(net, d) + sol = solve!(prob) + @test sol.psh[3] ≈ d[3] atol=1e-5 + @test sol.va[reference_buses(net)] ≈ zeros(2) atol=1e-8 + @test norm(kkt(flatten_variables(sol, prob), prob, d), Inf) < 1e-4 + + congestion = calc_congestion_component(sol, net) + energy = calc_energy_component(sol, net) + @test calc_lmp(sol, net) ≈ energy + congestion atol=1e-10 + @test congestion[reference_buses(net)] == zeros(2) + end + + @testset "fully isolated buses" begin + net = _make_fully_isolated_network() + d = [0.2, 0.3] + @test reference_buses(net) == [1, 2] + + state = DCPowerFlowState(net, d, d) + @test isempty(state.non_ref) + @test isempty(state.f) + @test state.va == zeros(2) + @test Matrix(calc_sensitivity(state, :va, :d)) == zeros(2, 2) + + prob = DCOPFProblem(net, d) + sol = solve!(prob) + @test sol.pg ≈ d atol=1e-5 + @test sol.va == zeros(2) + @test norm(kkt(flatten_variables(sol, prob), prob, d), Inf) < 1e-4 + end + + @testset "bridge opening resizes DC topology workspaces" begin + net = _make_bridge_network() + d = [0.0, 0.5, 0.0, 0.5] + prob = DCOPFProblem(net, d) + solve!(prob) + + @test reference_buses(net) == [1] + dim_connected = kkt_dims(prob) + out = zeros(net.n) + jvp!(out, prob, :lmp, :d, ones(net.n)) + @test length(prob.cache.work) == dim_connected + + sw = copy(net.sw) + sw[2] = 0.0 + update_switching!(prob, sw) + @test reference_buses(net) == [1, 3] + @test isnothing(prob.cache.work) + @test isnothing(prob.cache.b_r_factor) + @test length(prob.cons.ref) == 2 + + sol = solve!(prob) + @test kkt_dims(prob) == dim_connected + 1 + @test sol.va[reference_buses(net)] ≈ zeros(2) atol=1e-8 + @test norm(kkt(flatten_variables(sol, prob), prob, d), Inf) < 1e-4 + + energy = calc_energy_component(sol, net) + @test energy[1] ≈ energy[2] atol=1e-5 + @test energy[3] ≈ energy[4] atol=1e-5 + + sensitivity = calc_sensitivity(prob, :pg, :d) + @test all(isfinite, Matrix(sensitivity)) + jvp!(out, prob, :lmp, :d, ones(net.n)) + @test all(isfinite, out) + @test length(prob.cache.work) == kkt_dims(prob) + @test all(isfinite, vjp(prob, :lmp, :d, ones(net.n))) + + state = DCPowerFlowState(net, [0.5, 0.0, 0.5, 0.0], d) + @test state.va[reference_buses(net)] == zeros(2) + @test all(isfinite, Matrix(calc_sensitivity(state, :f, :d))) + end + + @testset "bundled PowerModels disconnected cases" begin + for case_name in ["case5_db.m", "case6.m"] + @testset "$case_name" begin + raw = load_raw_case(case_name) + if isnothing(raw) + @test_skip false + continue + end + + pm_result = PowerModels.solve_dc_opf(raw, + optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)) + @test pm_result["termination_status"] in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + + net = DCNetwork(raw) + d = calc_demand_vector(raw) + @test length(reference_buses(net)) > 1 + sol = solve!(DCOPFProblem(net, d)) + @test all(isfinite, sol.va) + @test all(isfinite, sol.pg) + @test sol.objective ≈ pm_result["objective"] rtol=1e-4 atol=1e-2 + end + end + end + + @testset "case7_tplgy specialized-component classification" begin + raw = load_raw_case("case7_tplgy.m") + if isnothing(raw) + @test_skip false + else + specialized = sum(length(get(raw, key, Dict())) for key in ("dcline", "storage", "switch")) + @test specialized > 0 + @test length(reference_buses(DCNetwork(raw))) > 1 + end + end +end From f68b1c682f5d40a8d536fe4af355615693df407c Mon Sep 17 00:00:00 2001 From: Cameron Khanpour <99142483+cameronkhanpour@users.noreply.github.com> Date: Tue, 2 Jun 2026 02:04:20 -0400 Subject: [PATCH 03/13] Cache DC island topology for KKT builders --- docs/src/api.md | 1 + src/prob/kkt_dc_opf.jl | 54 ++++++++++++++------ src/sens/cost.jl | 14 +++--- src/sens/flowlimit.jl | 7 ++- src/sens/susceptance.jl | 7 ++- src/types/dc_network.jl | 98 ++++++++++++++++++++++++++++--------- src/types/dc_opf_problem.jl | 2 +- test/test_dc_islands.jl | 23 +++++++++ 8 files changed, 149 insertions(+), 57 deletions(-) diff --git a/docs/src/api.md b/docs/src/api.md index 6ad29c5..c4ffa4a 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -61,6 +61,7 @@ update_switching! update_fmax! calc_demand_vector calc_susceptance_matrix +reference_buses calc_lmp calc_qlmp calc_congestion_component diff --git a/src/prob/kkt_dc_opf.jl b/src/prob/kkt_dc_opf.jl index 2da407a..4b2f243 100644 --- a/src/prob/kkt_dc_opf.jl +++ b/src/prob/kkt_dc_opf.jl @@ -285,7 +285,11 @@ The KKT system includes: Total: 5n + 6m + 3k + n_ref """ -kkt_dims(prob::DCOPFProblem) = kkt_dims(prob.network) +function kkt_dims(prob::DCOPFProblem) + net = getfield(prob, :network) + return _dc_kkt_dims(getfield(net, :n), getfield(net, :m), getfield(net, :k), + length(getfield(prob, :cons).ref)) +end kkt_dims(n::Int, m::Int, k::Int) = _dc_kkt_dims(n, m, k, 1) _dc_kkt_dims(n::Int, m::Int, k::Int, n_ref::Int) = 5n + 6m + 3k + n_ref @@ -293,7 +297,7 @@ function kkt_dims(net::DCNetwork) n = getfield(net, :n) m = getfield(net, :m) k = getfield(net, :k) - n_ref = length(reference_buses(net)) + n_ref = length(_reference_buses(net)) # va(n) + pg(k) + f(m) + psh(n) + lam_lb(m) + lam_ub(m) + gamma_lb(m) + gamma_ub(m) + rho_lb(k) + rho_ub(k) + mu_lb(n) + mu_ub(n) + nu_bal(n) + nu_flow(m) + ref(n_ref) return _dc_kkt_dims(n, m, k, n_ref) end @@ -340,8 +344,30 @@ function _dc_kkt_indices(n::Int, m::Int, k::Int, n_ref::Int) ) end -kkt_indices(net::DCNetwork) = _dc_kkt_indices(net.n, net.m, net.k, length(reference_buses(net))) -kkt_indices(prob::DCOPFProblem) = kkt_indices(prob.network) +kkt_indices(net::DCNetwork) = _dc_kkt_indices(net.n, net.m, net.k, length(_reference_buses(net))) + +function kkt_indices(prob::DCOPFProblem) + net = getfield(prob, :network) + return _dc_kkt_indices(getfield(net, :n), getfield(net, :m), getfield(net, :k), + length(getfield(prob, :cons).ref)) +end + +function _dc_kkt_layout(net::DCNetwork) + n_ref = length(_reference_buses(net)) + n = getfield(net, :n) + m = getfield(net, :m) + k = getfield(net, :k) + return _dc_kkt_dims(n, m, k, n_ref), _dc_kkt_indices(n, m, k, n_ref) +end + +function _dc_kkt_layout(prob::DCOPFProblem) + net = getfield(prob, :network) + n_ref = length(getfield(prob, :cons).ref) + n = getfield(net, :n) + m = getfield(net, :m) + k = getfield(net, :k) + return _dc_kkt_dims(n, m, k, n_ref), _dc_kkt_indices(n, m, k, n_ref) +end # ============================================================================= # Variable Flattening/Unflattening @@ -464,7 +490,7 @@ function kkt(z::AbstractVector, net::DCNetwork, d::AbstractVector) B_mat = net.A' * W * net.A WA = W * net.A - refs = reference_buses(net) + refs = _reference_buses(net) # KKT conditions # 1. Stationarity w.r.t. θ @@ -586,7 +612,7 @@ function calc_kkt_jacobian(net::DCNetwork, d::AbstractVector, prob::DCOPFProblem n = getfield(net, :n) m = getfield(net, :m) k = getfield(net, :k) - refs = reference_buses(net) + refs = _reference_buses(net) n_ref = length(refs) dim = _dc_kkt_dims(n, m, k, n_ref) A = getfield(net, :A) @@ -832,8 +858,7 @@ function calc_kkt_jacobian_demand(net::DCNetwork, d::AbstractVector, sol::DCOPFS n = getfield(net, :n) m = getfield(net, :m) k = getfield(net, :k) - dim = kkt_dims(net) - idx = kkt_indices(net) + dim, idx = _dc_kkt_layout(net) mu_ub = getfield(sol, :mu_ub) colptr = Vector{Int}(undef, n + 1) @@ -861,8 +886,8 @@ Compute column `j` of the KKT parameter Jacobian ∂K/∂d. Only 1-2 nonzeros: always `nu_bal[j]`, plus `mu_ub[j]` if `d[j] > 0`. """ function calc_kkt_jacobian_demand_column(net::DCNetwork, d::AbstractVector, sol::DCOPFSolution, j::Int) - col = zeros(kkt_dims(net)) - idx = kkt_indices(net) + dim, idx = _dc_kkt_layout(net) + col = zeros(dim) col[idx.nu_bal[j]] = -1.0 col[idx.mu_ub[j]] = sol.mu_ub[j] * _shed_capacity_derivative(d[j]) return col @@ -896,7 +921,7 @@ enabling gradient-based optimization for topology control. function calc_kkt_jacobian_switching(prob::DCOPFProblem, sol::DCOPFSolution) net = prob.network n, m, k = net.n, net.m, net.k - dim = kkt_dims(net) + dim, idx = _dc_kkt_layout(prob) θ = sol.va @@ -907,9 +932,6 @@ function calc_kkt_jacobian_switching(prob::DCOPFProblem, sol::DCOPFSolution) J_s = spzeros(dim, m) - # Use centralized index calculation - idx = kkt_indices(net) - # Precompute A * θ once (invariant across branches) Aθ = A * θ @@ -963,8 +985,8 @@ Compute column `e` of the KKT parameter Jacobian ∂K/∂sw. """ function calc_kkt_jacobian_switching_column(prob::DCOPFProblem, sol::DCOPFSolution, e::Int) net = prob.network - col = zeros(kkt_dims(net)) - idx = kkt_indices(net) + dim, idx = _dc_kkt_layout(prob) + col = zeros(dim) A = net.A; b = net.b; θ = sol.va Aθ_e = dot(A[e, :], θ) f_bus, t_bus = _branch_bus_indices(A, e) diff --git a/src/sens/cost.jl b/src/sens/cost.jl index 980c8a6..72aa176 100644 --- a/src/sens/cost.jl +++ b/src/sens/cost.jl @@ -37,8 +37,7 @@ function calc_kkt_jacobian_cost_linear(net::DCNetwork) n = getfield(net, :n) m = getfield(net, :m) k = getfield(net, :k) - dim = kkt_dims(net) - idx = kkt_indices(net) + dim, idx = _dc_kkt_layout(net) colptr = Vector{Int}(undef, k + 1) rowval = Int[] @@ -62,8 +61,8 @@ end Compute column `j` of ∂K/∂cl. Only 1 nonzero: `pg[j] = 1.0`. """ function calc_kkt_jacobian_cost_linear_column(net::DCNetwork, j::Int) - col = zeros(kkt_dims(net)) - idx = kkt_indices(net) + dim, idx = _dc_kkt_layout(net) + col = zeros(dim) col[idx.pg[j]] = 1.0 return col end @@ -93,8 +92,7 @@ function calc_kkt_jacobian_cost_quadratic(prob::DCOPFProblem, sol::DCOPFSolution n = getfield(net, :n) m = getfield(net, :m) k = getfield(net, :k) - dim = kkt_dims(net) - idx = kkt_indices(net) + dim, idx = _dc_kkt_layout(prob) g = getfield(sol, :pg) @@ -122,8 +120,8 @@ end Compute column `j` of ∂K/∂cq. Only 1 nonzero: `pg[j] = 2*g[j]`. """ function calc_kkt_jacobian_cost_quadratic_column(net::DCNetwork, sol::DCOPFSolution, j::Int) - col = zeros(kkt_dims(net)) - idx = kkt_indices(net) + dim, idx = _dc_kkt_layout(net) + col = zeros(dim) col[idx.pg[j]] = 2.0 * sol.pg[j] return col end diff --git a/src/sens/flowlimit.jl b/src/sens/flowlimit.jl index fe0571a..e820c84 100644 --- a/src/sens/flowlimit.jl +++ b/src/sens/flowlimit.jl @@ -53,8 +53,7 @@ function calc_kkt_jacobian_flowlimit(prob::DCOPFProblem, sol::DCOPFSolution) n = getfield(net, :n) m = getfield(net, :m) k = getfield(net, :k) - dim = kkt_dims(net) - idx = kkt_indices(net) + dim, idx = _dc_kkt_layout(prob) lambda_lb = getfield(sol, :lam_lb) lambda_ub = getfield(sol, :lam_ub) @@ -84,8 +83,8 @@ end Compute column `e` of ∂K/∂fmax. Only 2 nonzeros: `lam_lb[e]` and `lam_ub[e]`. """ function calc_kkt_jacobian_flowlimit_column(net::DCNetwork, sol::DCOPFSolution, e::Int) - col = zeros(kkt_dims(net)) - idx = kkt_indices(net) + dim, idx = _dc_kkt_layout(net) + col = zeros(dim) col[idx.lam_lb[e]] = sol.lam_lb[e] col[idx.lam_ub[e]] = sol.lam_ub[e] return col diff --git a/src/sens/susceptance.jl b/src/sens/susceptance.jl index 4729e10..37ba257 100644 --- a/src/sens/susceptance.jl +++ b/src/sens/susceptance.jl @@ -50,8 +50,7 @@ Derivatives: function calc_kkt_jacobian_susceptance(prob::DCOPFProblem, sol::DCOPFSolution) net = prob.network n, m, k = net.n, net.m, net.k - dim = kkt_dims(net) - idx = kkt_indices(net) + dim, idx = _dc_kkt_layout(prob) theta = sol.va nu_bal = sol.nu_bal @@ -91,8 +90,8 @@ Compute column `e` of ∂K/∂b. ~6 nonzeros from the incidence structure of bra """ function calc_kkt_jacobian_susceptance_column(prob::DCOPFProblem, sol::DCOPFSolution, e::Int) net = prob.network - col = zeros(kkt_dims(net)) - idx = kkt_indices(net) + dim, idx = _dc_kkt_layout(prob) + col = zeros(dim) A = net.A; sw = net.sw; θ = sol.va Aθ_e = dot(A[e, :], θ) f_bus, t_bus = _branch_bus_indices(A, e) diff --git a/src/types/dc_network.jl b/src/types/dc_network.jl index 2d32ff9..b526e8e 100644 --- a/src/types/dc_network.jl +++ b/src/types/dc_network.jl @@ -19,6 +19,21 @@ # DC network representation for B-theta OPF formulation with susceptance-weighted # Laplacian B = A' * Diag(-b .* sw) * A. +""" +Internal cache for the energized DC topology. The incidence matrix structure is +fixed after construction, while `b` and `sw` may change in place. +""" +mutable struct _DCTopologyCache + from_bus::Vector{Int} + to_bus::Vector{Int} + energized::BitVector + refs::Vector{Int} + non_ref::Vector{Int} + initialized::Bool +end + +_DCTopologyCache() = _DCTopologyCache(Int[], Int[], BitVector(), Int[], Int[], false) + """ DCNetwork <: AbstractPowerNetwork @@ -63,6 +78,7 @@ struct DCNetwork <: AbstractPowerNetwork ref_bus::Int tau::Float64 id_map::IDMapping + topology_cache::_DCTopologyCache end # ============================================================================= @@ -536,7 +552,8 @@ function DCNetwork(data::NamedTuple; tau::Float64=DEFAULT_TAU, ref_bus::Union{No end return DCNetwork(n, m, k, A, G_inc, b, sw, fmax, gmax, gmin, angmax, angmin, - cq, cl, c_shed, demand, pg_init, ref_bus, tau, id_map) + cq, cl, c_shed, demand, pg_init, ref_bus, tau, id_map, + _DCTopologyCache()) end """ @@ -576,7 +593,8 @@ function DCNetwork( Float64.(c_shed), Float64.(demand), Float64.(pg_init), ref_bus, tau, - IDMapping(n, m, k) + IDMapping(n, m, k), + _DCTopologyCache() ) end @@ -624,16 +642,23 @@ function calc_susceptance_matrix(network::DCNetwork) return sparse(network.A' * W * network.A) end -""" - reference_buses(net::DCNetwork) → Vector{Int} - -Return one deterministic reference bus for each energized island. +@inline _is_energized(net::DCNetwork, e::Int) = + !iszero(getfield(net, :b)[e] * getfield(net, :sw)[e]) + +function _topology_cache_valid(net::DCNetwork) + cache = getfield(net, :topology_cache) + cache.initialized || return false + energized = cache.energized + m = getfield(net, :m) + length(energized) == m || return false + @inbounds for e in 1:m + energized[e] == _is_energized(net, e) || return false + end + return true +end -The configured `net.ref_bus` is preserved for its island. Every other island, -including an isolated bus, uses its lowest sequential bus index. A branch is -energized when `b[e] * sw[e] != 0`. -""" -function reference_buses(net::DCNetwork) +function _refresh_topology_cache!(net::DCNetwork) + cache = net.topology_cache parent = collect(1:net.n) function find_root(i::Int) @@ -656,22 +681,26 @@ function reference_buses(net::DCNetwork) return nothing end - from_bus = zeros(Int, net.m) - to_bus = zeros(Int, net.m) - rows, cols, nz_values = findnz(net.A) - @inbounds for p in eachindex(rows) - if nz_values[p] > 0 - from_bus[rows[p]] = cols[p] - else - to_bus[rows[p]] = cols[p] + if length(cache.from_bus) != net.m + cache.from_bus = zeros(Int, net.m) + cache.to_bus = zeros(Int, net.m) + rows, cols, nz_values = findnz(net.A) + @inbounds for p in eachindex(rows) + if nz_values[p] > 0 + cache.from_bus[rows[p]] = cols[p] + else + cache.to_bus[rows[p]] = cols[p] + end end end + resize!(cache.energized, net.m) @inbounds for e in 1:net.m - iszero(net.b[e] * net.sw[e]) && continue - from_bus[e] > 0 && to_bus[e] > 0 || + cache.energized[e] = _is_energized(net, e) + cache.energized[e] || continue + cache.from_bus[e] > 0 && cache.to_bus[e] > 0 || error("Incidence matrix row $e must have exactly two nonzero entries") - union_roots!(from_bus[e], to_bus[e]) + union_roots!(cache.from_bus[e], cache.to_bus[e]) end refs_by_root = Dict{Int,Int}() @@ -680,11 +709,32 @@ function reference_buses(net::DCNetwork) refs_by_root[root] = min(get(refs_by_root, root, bus), bus) end refs_by_root[find_root(net.ref_bus)] = net.ref_bus - return sort!(collect(values(refs_by_root))) + cache.refs = sort!(collect(values(refs_by_root))) + cache.non_ref = setdiff(1:net.n, cache.refs) + cache.initialized = true + return cache +end + +function _topology_cache(net::DCNetwork) + _topology_cache_valid(net) || _refresh_topology_cache!(net) + return getfield(net, :topology_cache) end +_reference_buses(net::DCNetwork) = _topology_cache(net).refs + +""" + reference_buses(net::DCNetwork) → Vector{Int} + +Return one deterministic reference bus for each energized island. + +The configured `net.ref_bus` is preserved for its island. Every other island, +including an isolated bus, uses its lowest sequential bus index. A branch is +energized when `b[e] * sw[e] != 0`. +""" +reference_buses(net::DCNetwork) = copy(_reference_buses(net)) + """Return bus indices after removing one reference bus per energized island.""" -_non_reference_buses(net::DCNetwork) = setdiff(1:net.n, reference_buses(net)) +_non_reference_buses(net::DCNetwork) = copy(_topology_cache(net).non_ref) """ _factorize_B_r(net::DCNetwork) → (factor, non_ref) diff --git a/src/types/dc_opf_problem.jl b/src/types/dc_opf_problem.jl index 30873f3..11bc075 100644 --- a/src/types/dc_opf_problem.jl +++ b/src/types/dc_opf_problem.jl @@ -246,7 +246,7 @@ function _rebuild_jump_model!(prob::DCOPFProblem) shed_ub = @constraint(model, psh .<= _shed_capacity.(d)) # One angle reference per energized island removes the Laplacian nullspace. - refs = reference_buses(network) + refs = _reference_buses(network) ref_con = @constraint(model, [i in refs], va[i] == 0.0) # Open lines should not constrain angle differences. diff --git a/test/test_dc_islands.jl b/test/test_dc_islands.jl index 3506443..68a50e1 100644 --- a/test/test_dc_islands.jl +++ b/test/test_dc_islands.jl @@ -119,6 +119,29 @@ end @test all(isfinite, Matrix(calc_sensitivity(state, :f, :d))) end + @testset "reference bus cache follows direct topology mutation" begin + net = _make_bridge_network() + dim_connected = kkt_dims(net) + @test reference_buses(net) == [1] + + # Returned references are defensive copies of the internal cache. + refs = reference_buses(net) + refs[1] = 4 + @test reference_buses(net) == [1] + + net.sw[2] = 0.0 + @test reference_buses(net) == [1, 3] + @test kkt_dims(net) == dim_connected + 1 + + net.sw[2] = 1.0 + @test reference_buses(net) == [1] + @test kkt_dims(net) == dim_connected + + net.b[2] = 0.0 + @test reference_buses(net) == [1, 3] + @test kkt_dims(net) == dim_connected + 1 + end + @testset "bundled PowerModels disconnected cases" begin for case_name in ["case5_db.m", "case6.m"] @testset "$case_name" begin From 87687e578e4004383adc33de13defe39f6dcf23c Mon Sep 17 00:00:00 2001 From: Cameron Khanpour <99142483+cameronkhanpour@users.noreply.github.com> Date: Tue, 2 Jun 2026 04:33:24 -0400 Subject: [PATCH 04/13] Cache DC problem reference count --- src/prob/kkt_dc_opf.jl | 6 +++--- src/types/dc_opf_problem.jl | 5 ++++- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/prob/kkt_dc_opf.jl b/src/prob/kkt_dc_opf.jl index 4b2f243..d61dca0 100644 --- a/src/prob/kkt_dc_opf.jl +++ b/src/prob/kkt_dc_opf.jl @@ -288,7 +288,7 @@ Total: 5n + 6m + 3k + n_ref function kkt_dims(prob::DCOPFProblem) net = getfield(prob, :network) return _dc_kkt_dims(getfield(net, :n), getfield(net, :m), getfield(net, :k), - length(getfield(prob, :cons).ref)) + getfield(prob, :_n_ref)) end kkt_dims(n::Int, m::Int, k::Int) = _dc_kkt_dims(n, m, k, 1) _dc_kkt_dims(n::Int, m::Int, k::Int, n_ref::Int) = 5n + 6m + 3k + n_ref @@ -349,7 +349,7 @@ kkt_indices(net::DCNetwork) = _dc_kkt_indices(net.n, net.m, net.k, length(_refer function kkt_indices(prob::DCOPFProblem) net = getfield(prob, :network) return _dc_kkt_indices(getfield(net, :n), getfield(net, :m), getfield(net, :k), - length(getfield(prob, :cons).ref)) + getfield(prob, :_n_ref)) end function _dc_kkt_layout(net::DCNetwork) @@ -362,7 +362,7 @@ end function _dc_kkt_layout(prob::DCOPFProblem) net = getfield(prob, :network) - n_ref = length(getfield(prob, :cons).ref) + n_ref = getfield(prob, :_n_ref) n = getfield(net, :n) m = getfield(net, :m) k = getfield(net, :k) diff --git a/src/types/dc_opf_problem.jl b/src/types/dc_opf_problem.jl index 11bc075..a91deb9 100644 --- a/src/types/dc_opf_problem.jl +++ b/src/types/dc_opf_problem.jl @@ -123,6 +123,7 @@ B-θ formulation of DC OPF wrapped around a JuMP model. - `d`: Demand parameter (can be updated for sensitivity analysis) - `cons`: Named tuple of constraint references - `cache`: Mutable sensitivity cache for avoiding redundant KKT solves +- `_n_ref`: Number of energized-island reference constraints (internal) - `_optimizer`: Optimizer factory for model rebuilds (internal) - `_silent`: Whether to suppress solver output (internal) """ @@ -136,6 +137,7 @@ mutable struct DCOPFProblem{O} <: AbstractOPFProblem d::Vector{Float64} cons::NamedTuple cache::DCSensitivityCache + _n_ref::Int _optimizer::O _silent::Bool end @@ -178,7 +180,7 @@ function DCOPFProblem(network::DCNetwork, d::AbstractVector; optimizer=Ipopt.Opt prob = DCOPFProblem( JuMP.Model(), network, VariableRef[], VariableRef[], VariableRef[], VariableRef[], - Float64.(d), (;), DCSensitivityCache(), optimizer, silent + Float64.(d), (;), DCSensitivityCache(), 0, optimizer, silent ) _rebuild_jump_model!(prob) return prob @@ -258,6 +260,7 @@ function _rebuild_jump_model!(prob::DCOPFProblem) prob.pg = pg prob.f = f prob.psh = psh + prob._n_ref = length(refs) prob.cons = ( power_bal = power_bal, flow_def = flow_def, From 57a446c0ce3fde700ad09699f91ca56b18a8f651 Mon Sep 17 00:00:00 2001 From: ckhanpour3 Date: Sun, 7 Jun 2026 03:05:02 -0400 Subject: [PATCH 05/13] Document DC topology cache contracts --- docs/src/advanced.md | 11 ++++++++++- docs/src/math/dc-opf.md | 6 ++++++ docs/src/math/dc-power-flow.md | 6 ++++++ src/types/dc_network.jl | 24 ++++++++++++++++++++---- src/types/dc_opf_problem.jl | 12 +++++++++++- test/test_dc_islands.jl | 2 ++ 6 files changed, 55 insertions(+), 6 deletions(-) diff --git a/docs/src/advanced.md b/docs/src/advanced.md index 84ad36c..cf6289d 100644 --- a/docs/src/advanced.md +++ b/docs/src/advanced.md @@ -48,6 +48,15 @@ Use `reference_buses(net)` to obtain the effective reference set. It preserves `ref_bus` for its energized island and deterministically adds one reference for each additional island, including isolated buses. +`DCNetwork` precomputes an internal energized-topology cache and refreshes it +when topology readers observe a direct `b` or `sw` change. This cache is not a +thread-safety mechanism. Sharing a `DCNetwork` across threads is supported only +when topology fields are treated as read-only; callers that mutate `b` or `sw` +directly must serialize the mutation and the next topology-dependent read. For +`DCOPFProblem`, switch changes should go through [`update_switching!`](@ref), and +topology-changing susceptance edits require rebuilding the problem so the JuMP +model and KKT layout keep the same reference constraints. + ### ACNetwork Stores the AC network with vectorized admittance representation. @@ -77,7 +86,7 @@ The [`DCOPFProblem`](@ref) maintains a `DCSensitivityCache` that avoids redundan Calling `calc_sensitivity` with different operands for the same parameter reuses the cached KKT solve. For example, computing both `:va` and `:pg` w.r.t. `:d` only solves the KKT system once. -Cache invalidation happens automatically when `solve!`, `update_demand!`, `update_switching!`, or `update_fmax!` is called. +Cache invalidation happens automatically when `solve!`, `update_demand!`, `update_switching!`, or `update_fmax!` is called. Direct mutation of fields inside `prob.network` bypasses this contract. ### ACSensitivityCache diff --git a/docs/src/math/dc-opf.md b/docs/src/math/dc-opf.md index 40962e8..ae8d860 100644 --- a/docs/src/math/dc-opf.md +++ b/docs/src/math/dc-opf.md @@ -34,6 +34,12 @@ where: - ``\tau`` is a small regularization parameter for numerical conditioning - ``\text{refs}`` contains one deterministic reference bus per energized island, including isolated buses; the configured ``\text{ref_bus}`` remains the reference for its island +The built OPF model stores one reference constraint for each entry of +``\text{refs}``. Consequently, topology changes that can alter the energized +island partition must rebuild the model: use [`update_switching!`](@ref) for +switch changes, and rebuild the `DCOPFProblem` after direct susceptance edits +that move a branch across zero. + ## KKT System for Implicit Differentiation OPF sensitivities are computed via the implicit function theorem applied to the KKT conditions. At an optimal solution ``z^*``, the KKT residual satisfies ``K(z^*, p) = 0`` where ``z`` collects all primal and dual variables and ``p`` is a parameter. diff --git a/docs/src/math/dc-power-flow.md b/docs/src/math/dc-power-flow.md index 9cab6e6..902d443 100644 --- a/docs/src/math/dc-power-flow.md +++ b/docs/src/math/dc-power-flow.md @@ -24,6 +24,12 @@ B = A^\top \operatorname{diag}(-b \circ \mathrm{sw}) \, A where ``A`` is the ``m \times n`` incidence matrix and ``b`` stores the imaginary part of the inverse impedance (``b_e = \operatorname{Im}(1/z_e) < 0`` for inductive branches, so ``-b > 0``). +`DCNetwork` caches the energized-island partition used to choose reference +buses. Constructors initialize this cache, and topology-dependent readers +refresh it if direct `b` or `sw` edits change which branches are energized. The +cache assumes serialized topology mutation: a shared `DCNetwork` may be read +from multiple threads only while its topology fields are not being mutated. + ## Switching Sensitivity Switching sensitivity follows from matrix perturbation theory. For a branch ``e`` with switching state ``\mathrm{sw}_e \in [0,1]``: diff --git a/src/types/dc_network.jl b/src/types/dc_network.jl index b526e8e..e54c4dd 100644 --- a/src/types/dc_network.jl +++ b/src/types/dc_network.jl @@ -22,6 +22,11 @@ """ Internal cache for the energized DC topology. The incidence matrix structure is fixed after construction, while `b` and `sw` may change in place. + +The cache is prewarmed by constructors and refreshed by topology readers when +`b` or `sw` changes. It is not a synchronization primitive: callers sharing a +`DCNetwork` across threads must treat topology fields as read-only, or serialize +mutations and the first topology read after each mutation. """ mutable struct _DCTopologyCache from_bus::Vector{Int} @@ -56,6 +61,9 @@ topology sensitivity analysis. - `id_map`: Bidirectional mapping between original and sequential element IDs - `demand`: Real power demand aggregated per bus - `pg_init`: Initial real generation aggregated per bus +- `topology_cache`: Internal energized island cache. This cache is mutable even + though `DCNetwork` is an immutable struct; concurrent direct mutations of + `b`/`sw` are unsupported. """ struct DCNetwork <: AbstractPowerNetwork n::Int @@ -551,9 +559,11 @@ function DCNetwork(data::NamedTuple; tau::Float64=DEFAULT_TAU, ref_bus::Union{No end end - return DCNetwork(n, m, k, A, G_inc, b, sw, fmax, gmax, gmin, angmax, angmin, - cq, cl, c_shed, demand, pg_init, ref_bus, tau, id_map, - _DCTopologyCache()) + net = DCNetwork(n, m, k, A, G_inc, b, sw, fmax, gmax, gmin, angmax, angmin, + cq, cl, c_shed, demand, pg_init, ref_bus, tau, id_map, + _DCTopologyCache()) + _refresh_topology_cache!(net) + return net end """ @@ -583,7 +593,7 @@ function DCNetwork( length(demand) == n || throw(DimensionMismatch("demand length $(length(demand)) must match number of buses $n")) length(pg_init) == n || throw(DimensionMismatch("pg_init length $(length(pg_init)) must match number of buses $n")) all(c_shed .> 0) || throw(ArgumentError("c_shed must be strictly positive at all buses")) - return DCNetwork( + net = DCNetwork( n, m, k, sparse(Float64.(A)), sparse(Float64.(G_inc)), Float64.(b), Float64.(sw), @@ -596,6 +606,8 @@ function DCNetwork( IDMapping(n, m, k), _DCTopologyCache() ) + _refresh_topology_cache!(net) + return net end # ============================================================================= @@ -716,6 +728,10 @@ function _refresh_topology_cache!(net::DCNetwork) end function _topology_cache(net::DCNetwork) + # Refresh mutates `net.topology_cache`. Constructors prewarm this cache, so + # normal read-only sharing across threads does not first-touch it. If callers + # mutate `b` or `sw` directly, they must serialize that mutation and the next + # topology read; the exposed vectors themselves are not thread-safe. _topology_cache_valid(net) || _refresh_topology_cache!(net) return getfield(net, :topology_cache) end diff --git a/src/types/dc_opf_problem.jl b/src/types/dc_opf_problem.jl index a91deb9..f3ebf20 100644 --- a/src/types/dc_opf_problem.jl +++ b/src/types/dc_opf_problem.jl @@ -123,7 +123,9 @@ B-θ formulation of DC OPF wrapped around a JuMP model. - `d`: Demand parameter (can be updated for sensitivity analysis) - `cons`: Named tuple of constraint references - `cache`: Mutable sensitivity cache for avoiding redundant KKT solves -- `_n_ref`: Number of energized-island reference constraints (internal) +- `_n_ref`: Number of energized-island reference constraints in the currently + built JuMP model (internal). Invariant after `_rebuild_jump_model!`: + `_n_ref == length(cons.ref)`. - `_optimizer`: Optimizer factory for model rebuilds (internal) - `_silent`: Whether to suppress solver output (internal) """ @@ -195,6 +197,12 @@ _is_ipopt_optimizer(opt::MOI.OptimizerWithAttributes) = opt.optimizer_constructo Build (or rebuild) the JuMP model from current network parameters. Called by the constructor and by `update_switching!` after mutating `network.sw`. + +This function owns the `_n_ref == length(prob.cons.ref)` invariant. Directly +mutating `prob.network.sw` or moving `prob.network.b` across zero changes the +energized-island topology without rebuilding the JuMP model, so callers must use +`update_switching!` for switch changes and rebuild the problem after +topology-changing susceptance edits. """ function _rebuild_jump_model!(prob::DCOPFProblem) network = prob.network @@ -260,6 +268,8 @@ function _rebuild_jump_model!(prob::DCOPFProblem) prob.pg = pg prob.f = f prob.psh = psh + # Cache the reference-constraint count for KKT layout hot paths. This value + # belongs to the built JuMP model and must stay equal to length(ref_con). prob._n_ref = length(refs) prob.cons = ( power_bal = power_bal, diff --git a/test/test_dc_islands.jl b/test/test_dc_islands.jl index 68a50e1..9c07d62 100644 --- a/test/test_dc_islands.jl +++ b/test/test_dc_islands.jl @@ -121,6 +121,7 @@ end @testset "reference bus cache follows direct topology mutation" begin net = _make_bridge_network() + @test getfield(net, :topology_cache).initialized dim_connected = kkt_dims(net) @test reference_buses(net) == [1] @@ -157,6 +158,7 @@ end net = DCNetwork(raw) d = calc_demand_vector(raw) + @test getfield(net, :topology_cache).initialized @test length(reference_buses(net)) > 1 sol = solve!(DCOPFProblem(net, d)) @test all(isfinite, sol.va) From 7ee4b3cb0439cb260573d7367dfeaeaa30b6081e Mon Sep 17 00:00:00 2001 From: samtalki <10187005+samtalki@users.noreply.github.com> Date: Mon, 15 Jun 2026 15:02:09 -0400 Subject: [PATCH 06/13] Fix DC island tests for PowerIO parser path --- test/test_dc_islands.jl | 31 +++++++++++++++++++------------ 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/test/test_dc_islands.jl b/test/test_dc_islands.jl index 9c07d62..3b968f6 100644 --- a/test/test_dc_islands.jl +++ b/test/test_dc_islands.jl @@ -144,10 +144,11 @@ end end @testset "bundled PowerModels disconnected cases" begin - for case_name in ["case5_db.m", "case6.m"] + for case_name in ["case6.m"] @testset "$case_name" begin raw = load_raw_case(case_name) - if isnothing(raw) + parsed = load_test_case(case_name) + if isnothing(raw) || isnothing(parsed) @test_skip false continue end @@ -156,8 +157,8 @@ end optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)) @test pm_result["termination_status"] in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) - net = DCNetwork(raw) - d = calc_demand_vector(raw) + net = DCNetwork(parsed) + d = calc_demand_vector(parsed) @test getfield(net, :topology_cache).initialized @test length(reference_buses(net)) > 1 sol = solve!(DCOPFProblem(net, d)) @@ -168,14 +169,20 @@ end end end - @testset "case7_tplgy specialized-component classification" begin - raw = load_raw_case("case7_tplgy.m") - if isnothing(raw) - @test_skip false - else - specialized = sum(length(get(raw, key, Dict())) for key in ("dcline", "storage", "switch")) - @test specialized > 0 - @test length(reference_buses(DCNetwork(raw))) > 1 + @testset "specialized PowerModels components stay unsupported" begin + for case_name in ["case5_db.m", "case7_tplgy.m"] + @testset "$case_name" begin + raw = load_raw_case(case_name) + parsed = load_test_case(case_name) + if isnothing(raw) || isnothing(parsed) + @test_skip false + continue + end + + specialized = sum(length(get(raw, key, Dict())) for key in ("dcline", "storage", "switch")) + @test specialized > 0 + @test_throws ArgumentError DCNetwork(parsed) + end end end end From f7c628e58fce38aceb4b8608b9f6956cb2b0f9b5 Mon Sep 17 00:00:00 2001 From: samtalki <10187005+samtalki@users.noreply.github.com> Date: Mon, 15 Jun 2026 15:07:17 -0400 Subject: [PATCH 07/13] Handle empty DC display ranges --- docs/src/advanced.md | 12 ++++++------ docs/src/math/dc-opf.md | 8 ++++---- docs/src/math/dc-power-flow.md | 8 ++++---- ext/PowerDiffAPFExt.jl | 2 +- src/prob/kkt_dc_opf.jl | 8 ++++---- src/sens/lmp.jl | 3 +-- src/types/dc_network.jl | 16 ++++++++-------- src/types/dc_opf_problem.jl | 10 +++++----- src/types/show.jl | 19 ++++++++++++++----- test/test_dc_islands.jl | 3 +++ 10 files changed, 50 insertions(+), 39 deletions(-) diff --git a/docs/src/advanced.md b/docs/src/advanced.md index cf6289d..57da28a 100644 --- a/docs/src/advanced.md +++ b/docs/src/advanced.md @@ -36,7 +36,7 @@ Stores the DC network topology and parameters. | `gmax`, `gmin` | `Vector{Float64}` | Generator limits | | `angmax`, `angmin` | `Vector{Float64}` | Phase angle difference limits | | `cq`, `cl` | `Vector{Float64}` | Cost coefficients (quadratic, linear) | -| `c_shed` | `Vector{Float64}` | Load-shedding cost per bus | +| `c_shed` | `Vector{Float64}` | Load shedding cost per bus | | `ref_bus` | `Int` | Preferred reference bus index (sequential) | | `tau` | `Float64` | Regularization parameter | | `id_map` | `IDMapping` | Bidirectional element ID mapping (original ↔ sequential) | @@ -48,13 +48,13 @@ Use `reference_buses(net)` to obtain the effective reference set. It preserves `ref_bus` for its energized island and deterministically adds one reference for each additional island, including isolated buses. -`DCNetwork` precomputes an internal energized-topology cache and refreshes it +`DCNetwork` precomputes an internal energized topology cache and refreshes it when topology readers observe a direct `b` or `sw` change. This cache is not a -thread-safety mechanism. Sharing a `DCNetwork` across threads is supported only -when topology fields are treated as read-only; callers that mutate `b` or `sw` -directly must serialize the mutation and the next topology-dependent read. For +thread safety mechanism. Sharing a `DCNetwork` across threads is supported only +when topology fields are treated as read only; callers that mutate `b` or `sw` +directly must serialize the mutation and the next topology dependent read. For `DCOPFProblem`, switch changes should go through [`update_switching!`](@ref), and -topology-changing susceptance edits require rebuilding the problem so the JuMP +topology changing susceptance edits require rebuilding the problem so the JuMP model and KKT layout keep the same reference constraints. ### ACNetwork diff --git a/docs/src/math/dc-opf.md b/docs/src/math/dc-opf.md index ae8d860..ac43451 100644 --- a/docs/src/math/dc-opf.md +++ b/docs/src/math/dc-opf.md @@ -95,7 +95,7 @@ Demand enters the power balance and the load shedding upper bound through ``` For strictly positive demand, ``\partial d_+ / \partial d = 1``. For negative -demand, it is ``0``. At zero demand, the positive part function is non-smooth; +demand, it is ``0``. At zero demand, the positive part function is nonsmooth; the implementation uses the fixed zero shedding convention already required by the collapsed bound ``0 \leq \text{psh} \leq 0``. @@ -109,9 +109,9 @@ Switching affects the Laplacian ``B``, weight matrix ``W``, and flow definition ``` This propagates into the stationarity, power balance, and flow definition blocks of the KKT system. -Sensitivities are defined while the energized-island partition is fixed. Opening +Sensitivities are defined while the energized island partition is fixed. Opening or closing a bridge splits or merges islands and changes the reference set, so -the derivative is non-smooth at that topology boundary. +the derivative is nonsmooth at that topology boundary. ### Cost Coefficients (``c_q``, ``c_l``) @@ -142,7 +142,7 @@ Susceptances affect the same blocks as switching (through ``B`` and ``W``), but Locational marginal prices are the power balance duals ``\nu_{\text{bal}}``, decomposed as: ```math -\text{LMP} = \underbrace{\text{energy}}_{\text{per-island uniform component}} + \underbrace{\text{congestion}}_{\text{flow-limit component}} +\text{LMP} = \underbrace{\text{energy}}_{\text{per island uniform component}} + \underbrace{\text{congestion}}_{\text{flow limit component}} ``` The congestion component is extracted by solving: diff --git a/docs/src/math/dc-power-flow.md b/docs/src/math/dc-power-flow.md index 902d443..c48f72f 100644 --- a/docs/src/math/dc-power-flow.md +++ b/docs/src/math/dc-power-flow.md @@ -24,8 +24,8 @@ B = A^\top \operatorname{diag}(-b \circ \mathrm{sw}) \, A where ``A`` is the ``m \times n`` incidence matrix and ``b`` stores the imaginary part of the inverse impedance (``b_e = \operatorname{Im}(1/z_e) < 0`` for inductive branches, so ``-b > 0``). -`DCNetwork` caches the energized-island partition used to choose reference -buses. Constructors initialize this cache, and topology-dependent readers +`DCNetwork` caches the energized island partition used to choose reference +buses. Constructors initialize this cache, and topology dependent readers refresh it if direct `b` or `sw` edits change which branches are energized. The cache assumes serialized topology mutation: a shared `DCNetwork` may be read from multiple threads only while its topology fields are not being mutated. @@ -45,9 +45,9 @@ where the perturbation is a rank-1 update from the incidence column of branch `` \frac{\partial B_r}{\partial \mathrm{sw}_e} = -b_e \, a_{e,r} \, a_{e,r}^\top ``` -These derivatives apply while the energized-island partition is fixed. A +These derivatives apply while the energized island partition is fixed. A bridge opening or closing changes the reference set, so sensitivities are -non-smooth at the split or merge boundary. +nonsmooth at the split or merge boundary. ### Flow Sensitivity to Switching diff --git a/ext/PowerDiffAPFExt.jl b/ext/PowerDiffAPFExt.jl index c837ba8..de87d1b 100644 --- a/ext/PowerDiffAPFExt.jl +++ b/ext/PowerDiffAPFExt.jl @@ -57,7 +57,7 @@ end Convert a `DCNetwork` to an `APF.Network`. -APF networks lack generators, costs, and limits, so this is one-way. +APF networks lack generators, costs, and limits, so this is one way. APF exposes one slack bus, so disconnected `DCNetwork`s are rejected. Bus demand is set to zero (PD separates demand from network topology). Branch `status` is derived from switching state: `sw[e] > 0.5`. diff --git a/src/prob/kkt_dc_opf.jl b/src/prob/kkt_dc_opf.jl index d61dca0..7d6968d 100644 --- a/src/prob/kkt_dc_opf.jl +++ b/src/prob/kkt_dc_opf.jl @@ -952,7 +952,7 @@ function calc_kkt_jacobian_switching(prob::DCOPFProblem, sol::DCOPFSolution) Ae_dot_ν = dot(A_e_vec, ν_bal) J_s[idx.va, e] = -b[e] * A_e_vec * (Ae_dot_ν + ν_flow[e]) - # ∂K_θ/∂s_e from gated angle-difference bounds + # ∂K_θ/∂s_e from gated angle difference bounds J_s[idx.va, e] += A_e_vec * (sol.gamma_ub[e] - sol.gamma_lb[e]) # ∂K_γ/∂s_e from gated complementary slackness @@ -1000,11 +1000,11 @@ function calc_kkt_jacobian_switching_column(prob::DCOPFProblem, sol::DCOPFSoluti coeff = -b[e] * (sol.nu_bal[f_bus] - sol.nu_bal[t_bus] + sol.nu_flow[e]) col[idx.va[f_bus]] += coeff col[idx.va[t_bus]] -= coeff - # gated angle-difference stationarity contribution + # gated angle difference stationarity contribution coeff_ang = sol.gamma_ub[e] - sol.gamma_lb[e] col[idx.va[f_bus]] += coeff_ang col[idx.va[t_bus]] -= coeff_ang - # gated angle-difference complementary slackness + # gated angle difference complementary slackness col[idx.gamma_lb[e]] = sol.gamma_lb[e] * (Aθ_e - net.angmin[e]) col[idx.gamma_ub[e]] = sol.gamma_ub[e] * (net.angmax[e] - Aθ_e) return col @@ -1025,7 +1025,7 @@ function update_switching!(prob::DCOPFProblem, s::AbstractVector) length(s) == m || throw(DimensionMismatch("Switching vector length $(length(s)) must match number of branches $m")) all(0 .<= s .<= 1) || throw(ArgumentError("Switching values must be in [0,1]")) - # Invalidate all cached data including topology-dependent b_r_factor + # Invalidate all cached data including topology dependent b_r_factor invalidate_topology!(prob.cache) # Update network switching state diff --git a/src/sens/lmp.jl b/src/sens/lmp.jl index c558023..65bce40 100644 --- a/src/sens/lmp.jl +++ b/src/sens/lmp.jl @@ -83,7 +83,7 @@ The congestion RHS includes both flow limit duals and angle difference duals: congestion[non_ref] = B_r \\ (A' W (λ_ub - λ_lb) + A'(γ_ub - γ_lb))[non_ref] The congestion component captures price differentiation due to binding flow and angle -constraints, with each energized-island reference bus congestion component equal +constraints, with each energized island reference bus congestion component equal to zero. # Returns @@ -177,4 +177,3 @@ function calc_qlmp(prob::ACOPFProblem) sol = _ensure_ac_solved!(prob) return calc_qlmp(sol, prob) end - diff --git a/src/types/dc_network.jl b/src/types/dc_network.jl index e54c4dd..494dbaa 100644 --- a/src/types/dc_network.jl +++ b/src/types/dc_network.jl @@ -25,7 +25,7 @@ fixed after construction, while `b` and `sw` may change in place. The cache is prewarmed by constructors and refreshed by topology readers when `b` or `sw` changes. It is not a synchronization primitive: callers sharing a -`DCNetwork` across threads must treat topology fields as read-only, or serialize +`DCNetwork` across threads must treat topology fields as read only, or serialize mutations and the first topology read after each mutation. """ mutable struct _DCTopologyCache @@ -55,7 +55,7 @@ topology sensitivity analysis. - `fmax`, `gmax`, `gmin`: Flow and generation limits - `angmax`, `angmin`: Phase angle difference limits - `cq`, `cl`: Quadratic and linear generation cost coefficients -- `c_shed`: Load-shedding cost per bus (penalty for involuntary load curtailment) +- `c_shed`: Load shedding cost per bus (penalty for involuntary load curtailment) - `ref_bus`: Preferred reference bus index (phase angle = 0) - `tau`: Regularization parameter for strong convexity - `id_map`: Bidirectional mapping between original and sequential element IDs @@ -188,7 +188,7 @@ const DEFAULT_SHED_COST_MULTIPLIER = 10 # the network tables the DCNetwork and ACNetwork constructors consume. The only # logic beyond re-keying to source bus ids is the OPF solver modeling PowerIO leaves # to the consumer: polynomial cost interpretation, finite flow limits, default -# angle-difference bounds, and rejection of records PowerDiff does not model. +# angle difference bounds, and rejection of records PowerDiff does not model. """ parse_file(path::String; library=nothing, from=nothing, filetype=nothing) -> PowerIO.Network @@ -298,7 +298,7 @@ load/shunt aggregation, reference-bus inference (`type == 3`), source bus ids on This adapter keys bus references back to source bus ids (so [`IDMapping`](@ref)'s sorted ordering is preserved) and applies the OPF modeling PowerIO leaves to the consumer: polynomial cost interpretation (rejecting PWL and higher-than-quadratic), -a finite flow-limit fallback when `rate_a == 0`, default angle-difference bounds, +a finite flow limit fallback when `rate_a == 0`, default angle difference bounds, and rejection of storage / HVDC records that PowerDiff does not model. The returned `bus`/`gen`/`branch` rows mirror the field names the network @@ -428,7 +428,7 @@ function _fallback_rate_a(r::Float64, x::Float64, angmin::Float64, angmax::Float return ymag * max(fr_vmax, to_vmax) * cmax end -# Default angle-difference bounds (radians in, radians out). MATPOWER angmin == angmax +# Default angle difference bounds (radians in, radians out). MATPOWER angmin == angmax # == 0 means unbounded; treat ±90° or wider and the zero case as a ±60° window, the # MATPOWER/PowerModels convention. PowerIO's `to_powerdata` already converts to radians. function _normalize_angle_bounds(angmin::Float64, angmax::Float64) @@ -532,7 +532,7 @@ function DCNetwork(data::NamedTuple; tau::Float64=DEFAULT_TAU, ref_bus::Union{No demand = calc_demand_vector(data, id_map) pg_init = _calc_generation_vector(data, id_map) - # Load-shedding cost: high penalty to discourage shedding when feasible. + # Load shedding cost: high penalty to discourage shedding when feasible. # Guard the reduction so a generator-free network (valid for pure DC power flow # built via the NamedTuple constructor) falls back to a unit marginal cost # instead of `maximum` throwing on an empty collection. @@ -729,9 +729,9 @@ end function _topology_cache(net::DCNetwork) # Refresh mutates `net.topology_cache`. Constructors prewarm this cache, so - # normal read-only sharing across threads does not first-touch it. If callers + # normal read only sharing across threads does not first touch it. If callers # mutate `b` or `sw` directly, they must serialize that mutation and the next - # topology read; the exposed vectors themselves are not thread-safe. + # topology read; the exposed vectors themselves are not thread safe. _topology_cache_valid(net) || _refresh_topology_cache!(net) return getfield(net, :topology_cache) end diff --git a/src/types/dc_opf_problem.jl b/src/types/dc_opf_problem.jl index f3ebf20..c671b4c 100644 --- a/src/types/dc_opf_problem.jl +++ b/src/types/dc_opf_problem.jl @@ -50,7 +50,7 @@ precomputed at construction time). - `dz_dsw`: Full KKT derivative w.r.t. switching (or nothing) - `dz_dfmax`: Full KKT derivative w.r.t. flow limits (or nothing) - `dz_db`: Full KKT derivative w.r.t. susceptances (or nothing) -- `b_r_factor`: Cached reduced susceptance factorization (topology-dependent, survives demand changes) +- `b_r_factor`: Cached reduced susceptance factorization (topology dependent, survives demand changes) - `work`: Scratch workspace for VJP/JVP KKT solves (lazily allocated on first call) """ mutable struct DCSensitivityCache @@ -97,7 +97,7 @@ end """ invalidate_topology!(cache::DCSensitivityCache) -Clear all cached data including topology-dependent `b_r_factor` and `work`. +Clear all cached data including topology dependent `b_r_factor` and `work`. Called when network topology changes (switching, susceptances). """ function invalidate_topology!(cache::DCSensitivityCache) @@ -123,7 +123,7 @@ B-θ formulation of DC OPF wrapped around a JuMP model. - `d`: Demand parameter (can be updated for sensitivity analysis) - `cons`: Named tuple of constraint references - `cache`: Mutable sensitivity cache for avoiding redundant KKT solves -- `_n_ref`: Number of energized-island reference constraints in the currently +- `_n_ref`: Number of energized island reference constraints in the currently built JuMP model (internal). Invariant after `_rebuild_jump_model!`: `_n_ref == length(cons.ref)`. - `_optimizer`: Optimizer factory for model rebuilds (internal) @@ -200,9 +200,9 @@ Called by the constructor and by `update_switching!` after mutating `network.sw` This function owns the `_n_ref == length(prob.cons.ref)` invariant. Directly mutating `prob.network.sw` or moving `prob.network.b` across zero changes the -energized-island topology without rebuilding the JuMP model, so callers must use +energized island topology without rebuilding the JuMP model, so callers must use `update_switching!` for switch changes and rebuild the problem after -topology-changing susceptance edits. +topology changing susceptance edits. """ function _rebuild_jump_model!(prob::DCOPFProblem) network = prob.network diff --git a/src/types/show.jl b/src/types/show.jl index ccb7d98..c6d6eec 100644 --- a/src/types/show.jl +++ b/src/types/show.jl @@ -66,6 +66,15 @@ Base.getproperty(sol::DCOPFSolution, s::Symbol) = _alias_getproperty(sol, _DC_OP Base.getproperty(prob::DCOPFProblem, s::Symbol) = _alias_getproperty(prob, _DC_OPF_PROBLEM_ALIASES, s) Base.getproperty(sol::ACOPFSolution, s::Symbol) = _alias_getproperty(sol, _AC_OPF_SOLUTION_ALIASES, s) +_range_str(v; digits=2) = + isempty(v) ? "empty" : "[$(round(minimum(v); digits=digits)), $(round(maximum(v); digits=digits))]" + +_range_str(lo, hi; digits=2) = + isempty(lo) || isempty(hi) ? "empty" : "[$(round(minimum(lo); digits=digits)), $(round(maximum(hi); digits=digits))]" + +_max_abs_str(v; digits=2) = + isempty(v) ? "empty" : string(round(maximum(abs, v); digits=digits)) + # ============================================================================= # DCNetwork # ============================================================================= @@ -79,8 +88,8 @@ function Base.show(io::IO, ::MIME"text/plain", net::DCNetwork) println(io, " Reference buses: $(reference_buses(net))") n_open = count(x -> x < 1.0, net.sw) println(io, " Open branches: $n_open/$(net.m)") - println(io, " Flow limits: [$(round(minimum(net.fmax); digits=2)), $(round(maximum(net.fmax); digits=2))]") - println(io, " Gen capacity: [$(round(minimum(net.gmin); digits=2)), $(round(maximum(net.gmax); digits=2))]") + println(io, " Flow limits: $(_range_str(net.fmax))") + println(io, " Gen capacity: $(_range_str(net.gmin, net.gmax))") print(io, " tau = $(net.tau)") end @@ -111,7 +120,7 @@ end function Base.show(io::IO, ::MIME"text/plain", state::DCPowerFlowState) println(io, "DCPowerFlowState ($(state.net.n) buses, $(state.net.m) branches)") println(io, " max|va|: $(round(maximum(abs, state.va); digits=4)) rad") - println(io, " max|f|: $(round(maximum(abs, state.f); digits=4)) p.u.") + println(io, " max|f|: $(_max_abs_str(state.f; digits=4)) p.u.") print(io, " Total demand: $(round(sum(state.d); digits=2)) p.u.") end @@ -142,8 +151,8 @@ end function Base.show(io::IO, ::MIME"text/plain", sol::DCOPFSolution) println(io, "DCOPFSolution (objective = $(round(sol.objective; digits=2)))") - println(io, " Generators: $(length(sol.pg)) (range: [$(round(minimum(sol.pg); digits=2)), $(round(maximum(sol.pg); digits=2))])") - println(io, " Flows: $(length(sol.f)) (max |f| = $(round(maximum(abs, sol.f); digits=2)))") + println(io, " Generators: $(length(sol.pg)) (range: $(_range_str(sol.pg)))") + println(io, " Flows: $(length(sol.f)) (max |f| = $(_max_abs_str(sol.f)))") print(io, " Shedding: $(round(sum(sol.psh); digits=2)) p.u. total") end diff --git a/test/test_dc_islands.jl b/test/test_dc_islands.jl index 3b968f6..e17483b 100644 --- a/test/test_dc_islands.jl +++ b/test/test_dc_islands.jl @@ -64,18 +64,21 @@ end net = _make_fully_isolated_network() d = [0.2, 0.3] @test reference_buses(net) == [1, 2] + @test !isempty(sprint(show, MIME"text/plain"(), net)) state = DCPowerFlowState(net, d, d) @test isempty(state.non_ref) @test isempty(state.f) @test state.va == zeros(2) @test Matrix(calc_sensitivity(state, :va, :d)) == zeros(2, 2) + @test !isempty(sprint(show, MIME"text/plain"(), state)) prob = DCOPFProblem(net, d) sol = solve!(prob) @test sol.pg ≈ d atol=1e-5 @test sol.va == zeros(2) @test norm(kkt(flatten_variables(sol, prob), prob, d), Inf) < 1e-4 + @test !isempty(sprint(show, MIME"text/plain"(), sol)) end @testset "bridge opening resizes DC topology workspaces" begin From 5bafb47ef48a56c70c83d9e9ef55f6b2659d617f Mon Sep 17 00:00:00 2001 From: ckhanpour3 Date: Mon, 15 Jun 2026 16:20:12 -0400 Subject: [PATCH 08/13] Address review feedback on DC island support - Reject non-binary switching states in to_apf_network so PowerDiff's energized-island count (b*sw != 0) cannot disagree with APF's on/off branch binarization (sw > 0.5); add a regression test. - Annotate the DCOPFProblem placeholder constructor call so the positional arguments (empty variable/constraint fields, normalized demand, initial _n_ref, optimizer, silent) are self-explanatory. - Make the documented theta-stationarity consistent with the implementation: define E_ref as the n x n_ref reference-selection matrix and include the Diag(sw) gate on the angle-difference dual term across lmp.jl, susceptance.jl, kkt_dc_opf.jl, and dc-opf.md; carry the same gate in calc_congestion_component (a no-op for energized branches where sw == 1). - Clarify that the per-island reference choice is deterministic (configured ref_bus for its island, lowest sequential index otherwise) and that the energy component is exactly uniform per island only in the unregularized limit, in lmp.jl and the math/advanced docs. - Complete the DCNetwork field table in advanced.md (demand, pg_init, topology_cache). --- docs/src/advanced.md | 10 +++++++--- docs/src/math/dc-opf.md | 4 ++-- docs/src/math/dc-power-flow.md | 7 ++++--- ext/PowerDiffAPFExt.jl | 16 +++++++++++++++- src/prob/kkt_dc_opf.jl | 3 ++- src/sens/lmp.jl | 25 +++++++++++++++++-------- src/sens/susceptance.jl | 7 ++++--- src/types/dc_opf_problem.jl | 14 ++++++++++++-- test/test_apf_integration.jl | 15 +++++++++++++++ 9 files changed, 78 insertions(+), 23 deletions(-) diff --git a/docs/src/advanced.md b/docs/src/advanced.md index 57da28a..531f6b5 100644 --- a/docs/src/advanced.md +++ b/docs/src/advanced.md @@ -37,16 +37,20 @@ Stores the DC network topology and parameters. | `angmax`, `angmin` | `Vector{Float64}` | Phase angle difference limits | | `cq`, `cl` | `Vector{Float64}` | Cost coefficients (quadratic, linear) | | `c_shed` | `Vector{Float64}` | Load shedding cost per bus | +| `demand` | `Vector{Float64}` | Real power demand aggregated per bus | +| `pg_init` | `Vector{Float64}` | Initial real generation aggregated per bus | | `ref_bus` | `Int` | Preferred reference bus index (sequential) | | `tau` | `Float64` | Regularization parameter | | `id_map` | `IDMapping` | Bidirectional element ID mapping (original ↔ sequential) | +| `topology_cache` | `_DCTopologyCache` | Internal energized-island cache (not part of the public API) | Construct from a parsed MATPOWER network with `DCNetwork(parse_file("case14.m"))`, or with explicit parameters: `DCNetwork(n, m, k, A, G_inc, b; ...)`. -Use `reference_buses(net)` to obtain the effective reference set. It preserves -`ref_bus` for its energized island and deterministically adds one reference for -each additional island, including isolated buses. +Use `reference_buses(net)` to obtain the effective reference set. The choice is +deterministic: `ref_bus` is kept as the reference for its energized island, and +every other island (including an isolated bus) uses its lowest sequential bus +index. `DCNetwork` precomputes an internal energized topology cache and refreshes it when topology readers observe a direct `b` or `sw` change. This cache is not a diff --git a/docs/src/math/dc-opf.md b/docs/src/math/dc-opf.md index ac43451..145e0ca 100644 --- a/docs/src/math/dc-opf.md +++ b/docs/src/math/dc-opf.md @@ -32,7 +32,7 @@ where: - ``c_{\text{shed}}`` is the load shedding cost vector - ``d_+ = \max(d, 0)`` is the curtailable portion of signed net demand; negative net demand remains in power balance as an injection - ``\tau`` is a small regularization parameter for numerical conditioning -- ``\text{refs}`` contains one deterministic reference bus per energized island, including isolated buses; the configured ``\text{ref_bus}`` remains the reference for its island +- ``\text{refs}`` contains one reference bus per energized island, including isolated buses. The choice is deterministic: the configured ``\text{ref_bus}`` is the reference for its island, and every other island uses its lowest sequential bus index The built OPF model stores one reference constraint for each entry of ``\text{refs}``. Consequently, topology changes that can alter the energized @@ -64,7 +64,7 @@ with total dimension ``5n + 6m + 3k + n_{\text{ref}}``. The KKT residual ``K(z, p)`` consists of: -1. **Stationarity w.r.t. ``\theta``**: ``B^\top \nu_{\text{bal}} + (WA)^\top \nu_{\text{flow}} + E_{\text{ref}} \eta_{\text{ref}} + A^\top (\gamma_{\text{ub}} - \gamma_{\text{lb}}) = 0`` +1. **Stationarity w.r.t. ``\theta``**: ``B^\top \nu_{\text{bal}} + (WA)^\top \nu_{\text{flow}} + E_{\text{ref}} \eta_{\text{ref}} + A^\top \operatorname{diag}(\mathrm{sw}) (\gamma_{\text{ub}} - \gamma_{\text{lb}}) = 0``, where ``E_{\text{ref}} \in \mathbb{R}^{n \times n_{\text{ref}}}`` selects the per-island reference buses 2. **Stationarity w.r.t. ``g``**: ``2 C_q g + c_l - G_{\text{inc}}^\top \nu_{\text{bal}} - \rho_{\text{lb}} + \rho_{\text{ub}} = 0`` 3. **Stationarity w.r.t. ``f``**: ``\tau^2 f - \nu_{\text{flow}} - \lambda_{\text{lb}} + \lambda_{\text{ub}} = 0`` 4. **Stationarity w.r.t. psh**: ``c_{\text{shed}} - \nu_{\text{bal}} - \mu_{\text{lb}} + \mu_{\text{ub}} = 0`` diff --git a/docs/src/math/dc-power-flow.md b/docs/src/math/dc-power-flow.md index c48f72f..ce2ce13 100644 --- a/docs/src/math/dc-power-flow.md +++ b/docs/src/math/dc-power-flow.md @@ -45,9 +45,10 @@ where the perturbation is a rank-1 update from the incidence column of branch `` \frac{\partial B_r}{\partial \mathrm{sw}_e} = -b_e \, a_{e,r} \, a_{e,r}^\top ``` -These derivatives apply while the energized island partition is fixed. A -bridge opening or closing changes the reference set, so sensitivities are -nonsmooth at the split or merge boundary. +These derivatives apply while the energized island partition is fixed. Toggling +a bridge (a branch whose energization changes the island partition) adds or +removes a reference bus and changes the dimension of the reduced system ``B_r``, +so sensitivities are nonsmooth at that split or merge boundary. ### Flow Sensitivity to Switching diff --git a/ext/PowerDiffAPFExt.jl b/ext/PowerDiffAPFExt.jl index de87d1b..a888842 100644 --- a/ext/PowerDiffAPFExt.jl +++ b/ext/PowerDiffAPFExt.jl @@ -60,7 +60,10 @@ Convert a `DCNetwork` to an `APF.Network`. APF networks lack generators, costs, and limits, so this is one way. APF exposes one slack bus, so disconnected `DCNetwork`s are rejected. Bus demand is set to zero (PD separates demand from network topology). -Branch `status` is derived from switching state: `sw[e] > 0.5`. +Branch `status` is derived from switching state: `sw[e] > 0.5`. Because APF +models branches as strictly on/off, switching states must be binary (`sw[e] ∈ +{0, 1}`); fractional values are rejected so that PD's energized-island count +(`b[e] * sw[e] != 0`) cannot disagree with the converted APF topology. Note: `to_apf_network` sets bus demand to zero because PD separates demand from topology. For APF workflows that need demand data (e.g., `compute_flow!`), @@ -68,6 +71,17 @@ use `APF.from_power_models(pm_data)` directly instead. """ function PowerDiff.to_apf_network(net::PowerDiff.DCNetwork) n, m = net.n, net.m + + # APF binarizes branch status via sw[e] > 0.5, but PD treats a branch as + # energized whenever b[e] * sw[e] != 0. These two views only agree when every + # switching state is binary, so reject fractional sw to keep the island count + # used for the slack-bus check consistent with the topology APF actually builds. + all(s -> s == 0 || s == 1, net.sw) || throw(ArgumentError( + "AcceleratedDCPowerFlows conversion requires binary switching states (sw[e] ∈ {0, 1}); " * + "APF models branches as on/off via sw[e] > 0.5, so fractional switching would make the " * + "energized-island count disagree with the converted APF topology" + )) + refs = PowerDiff.reference_buses(net) length(refs) == 1 || throw(ArgumentError( "AcceleratedDCPowerFlows conversion requires one energized island; found $(length(refs))" diff --git a/src/prob/kkt_dc_opf.jl b/src/prob/kkt_dc_opf.jl index 7d6968d..68a1ec0 100644 --- a/src/prob/kkt_dc_opf.jl +++ b/src/prob/kkt_dc_opf.jl @@ -456,7 +456,8 @@ s.t. G_inc * g + psh - d = B * θ (ν_bal) # Returns Vector of KKT residuals (should be zero at optimum): -1. Stationarity w.r.t. θ: B' * ν_bal + (W*A)' * ν_flow + E_ref * η_ref + A'*(γ_ub - γ_lb) = 0 +1. Stationarity w.r.t. θ: B' * ν_bal + (W*A)' * ν_flow + E_ref * η_ref + A' * Diag(sw) * (γ_ub - γ_lb) = 0 + (E_ref is the n × n_ref selection matrix for the per-island reference buses) 2. Stationarity w.r.t. g: 2*Cq * g + cl - G_inc' * ν_bal - ρ_lb + ρ_ub = 0 3. Stationarity w.r.t. f: τ² * f - ν_flow - λ_lb + λ_ub = 0 4. Stationarity w.r.t. psh: c_shed - ν_bal - μ_lb + μ_ub = 0 diff --git a/src/sens/lmp.jl b/src/sens/lmp.jl index 65bce40..8b51248 100644 --- a/src/sens/lmp.jl +++ b/src/sens/lmp.jl @@ -76,15 +76,18 @@ end Extract the congestion component of LMPs for analysis. -From the θ-stationarity KKT condition: - B' * ν_bal + (WA)' * ν_flow + E_ref * η_ref + A'*(γ_ub - γ_lb) = 0 +From the θ-stationarity KKT condition, where `E_ref` is the `n × n_ref` matrix +that selects the one reference bus per energized island (so `E_ref * η_ref` is a +length-`n` vector aligned with the bus stationarity rows): + B' * ν_bal + (WA)' * ν_flow + E_ref * η_ref + A' Diag(sw) (γ_ub - γ_lb) = 0 -The congestion RHS includes both flow limit duals and angle difference duals: - congestion[non_ref] = B_r \\ (A' W (λ_ub - λ_lb) + A'(γ_ub - γ_lb))[non_ref] +Neglecting the O(τ²) flow regularization (so ν_flow ≈ λ_ub - λ_lb), the congestion +RHS gathers the flow limit and angle difference dual contributions: + congestion[non_ref] = B_r \\ (A' W (λ_ub - λ_lb) + A' Diag(sw) (γ_ub - γ_lb))[non_ref] The congestion component captures price differentiation due to binding flow and angle -constraints, with each energized island reference bus congestion component equal -to zero. +constraints. Only the non-reference rows are populated, so every energized island +reference bus has a congestion component of exactly zero. # Returns Vector (length n) of congestion contributions to each bus's LMP. @@ -95,7 +98,12 @@ function calc_congestion_component(sol::DCOPFSolution, net::DCNetwork; non_ref = _non_reference_buses(net) At = net.A' - rhs_full = At * Diagonal(w .* net.sw) * (sol.lam_ub - sol.lam_lb) + At * (sol.gamma_ub - sol.gamma_lb) + # Angle difference constraints are gated by `sw` in the model, so their + # stationarity contribution carries the same `Diag(sw)` factor. For energized + # branches `sw == 1` (a no-op); on de-energized branches the gate matches the + # zero angle dual the solver returns there. + rhs_full = At * Diagonal(w .* net.sw) * (sol.lam_ub - sol.lam_lb) + + At * Diagonal(net.sw) * (sol.gamma_ub - sol.gamma_lb) result = zeros(net.n) result[non_ref] = B_r_factor \ rhs_full[non_ref] @@ -108,7 +116,8 @@ end Extract the energy (non-congestion) component of LMPs for analysis. This is the uniform price component: energy = ν_bal - congestion. -It is approximately constant within each energized island. +It is exactly uniform within each energized island in the unregularized limit; +the O(τ²) flow regularization perturbs that uniformity slightly. # Returns Vector (length n) of energy contributions to each bus's LMP. diff --git a/src/sens/susceptance.jl b/src/sens/susceptance.jl index 37ba257..1fe599e 100644 --- a/src/sens/susceptance.jl +++ b/src/sens/susceptance.jl @@ -37,9 +37,10 @@ Susceptance b affects: - Weight matrix: W = Diag(-b .* sw) - Flow definition: f = W * A * theta -The affected KKT conditions are: -- K_theta = B' * nu_bal + (WA)' * nu_flow + E_ref * eta_ref + A'*(gamma_ub - gamma_lb) - (gamma term has no b-dependence, so ∂K_theta/∂b only comes from B and WA terms) +The affected KKT conditions are (E_ref is the n × n_ref selection matrix for the +per-island reference buses, so E_ref * eta_ref is a length-n vector): +- K_theta = B' * nu_bal + (WA)' * nu_flow + E_ref * eta_ref + A' * Diag(sw) * (gamma_ub - gamma_lb) + (the eta_ref and gamma terms have no b-dependence, so ∂K_theta/∂b only comes from B and WA terms) - K_power_bal = G_inc * g + psh - d - B * theta - K_flow_def = f - W * A * theta diff --git a/src/types/dc_opf_problem.jl b/src/types/dc_opf_problem.jl index c671b4c..3de0be0 100644 --- a/src/types/dc_opf_problem.jl +++ b/src/types/dc_opf_problem.jl @@ -180,9 +180,19 @@ function DCOPFProblem(network::DCNetwork, d::AbstractVector; optimizer=Ipopt.Opt _debug_negative_demand(d) + # The model/variable/constraint fields start empty; `_rebuild_jump_model!` + # below builds the JuMP model and fills va/pg/f/psh, cons, and _n_ref. + # Arguments are positional to match the DCOPFProblem field order. prob = DCOPFProblem( - JuMP.Model(), network, VariableRef[], VariableRef[], VariableRef[], VariableRef[], - Float64.(d), (;), DCSensitivityCache(), 0, optimizer, silent + JuMP.Model(), # model: replaced on rebuild + network, + VariableRef[], VariableRef[], VariableRef[], VariableRef[], # va, pg, f, psh: filled on rebuild + Float64.(d), # d: demand parameter, normalized to Float64 + (;), # cons: filled on rebuild + DCSensitivityCache(), # cache: empty sensitivity cache + 0, # _n_ref: set to length(refs) on rebuild + optimizer, # _optimizer: factory for model (re)builds + silent, # _silent: suppress solver output ) _rebuild_jump_model!(prob) return prob diff --git a/test/test_apf_integration.jl b/test/test_apf_integration.jl index 4631b98..3198e32 100644 --- a/test/test_apf_integration.jl +++ b/test/test_apf_integration.jl @@ -176,6 +176,21 @@ end @test_throws ArgumentError to_apf_network(net) end +@testset "to_apf_network rejects fractional switching" begin + net = load_test_case("case14.m") + if isnothing(net) + @test_skip false + else + dc_net = DCNetwork(net) + # Fractional sw cannot be represented by APF's on/off branch status and + # would make PD's energized-island count disagree with the APF topology. + dc_net.sw[1] = 0.3 + @test_throws ArgumentError to_apf_network(dc_net) + dc_net.sw[1] = 1.0 # restore; binary states convert fine + @test to_apf_network(dc_net) isa APF.Network + end +end + # ========================================================================= # Index alignment # ========================================================================= From 661908a6ac86e0997ca160a49c025f427bd975af Mon Sep 17 00:00:00 2001 From: ckhanpour3 Date: Mon, 15 Jun 2026 17:03:29 -0400 Subject: [PATCH 09/13] Add regression test for sw-gated angle duals in congestion split Locks in the fix from the prior commit: calc_congestion_component must gate the angle-difference dual term by sw to match the KKT theta-stationarity A' Diag(sw) (gamma_ub - gamma_lb). The test uses a fractionally switched branch with a binding angle limit and asserts the gated result differs from the old ungated form, which it does not for fully closed networks (sw == 1). --- test/test_angle_diff_duals.jl | 37 +++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/test/test_angle_diff_duals.jl b/test/test_angle_diff_duals.jl index 2ad5619..28d185b 100644 --- a/test/test_angle_diff_duals.jl +++ b/test/test_angle_diff_duals.jl @@ -384,4 +384,41 @@ import PowerDiff: kkt, kkt_indices, flatten_variables @test rel_err < 0.05 end end + + @testset "Congestion component gates angle duals by sw" begin + # Regression (PR #55): calc_congestion_component must gate the angle + # difference dual term by sw, matching the KKT theta-stationarity term + # A' * Diag(sw) * (gamma_ub - gamma_lb). That term used to be ungated, + # which gave a wrong energy/congestion split on any branch with sw != 1 + # and a nonzero angle dual. With sw == 1 everywhere the two forms are + # identical, so this only shows up on partially switched branches. + net = DCNetwork(n, m, k, A, G_inc, b; sw=[1.0, 0.5, 1.0], + fmax=fill(100.0, m), gmax=[5.0, 5.0], gmin=[0.0, 0.0], + angmax=[Float64(pi), 0.025, Float64(pi)], angmin=fill(-Float64(pi), m), + cl=[10.0, 50.0], cq=[1.0, 1.0], ref_bus=1, tau=0.01) + prob = DCOPFProblem(net, d) + sol = solve!(prob) + + # The fractionally switched branch (2) carries an active angle dual. + @test abs(sol.gamma_ub[2] - sol.gamma_lb[2]) > 1e-3 + + cong = calc_congestion_component(sol, net) + w = -net.b + non_ref = PowerDiff._non_reference_buses(net) + At = net.A' + + # Gated RHS: matches the implementation and the documented stationarity. + rhs_gated = At * Diagonal(w .* net.sw) * (sol.lam_ub - sol.lam_lb) + + At * Diagonal(net.sw) * (sol.gamma_ub - sol.gamma_lb) + expected = zeros(net.n) + expected[non_ref] = sol.B_r_factor \ rhs_gated[non_ref] + @test isapprox(cong, expected; atol=1e-8) + + # Ungated RHS (the pre-fix behavior) disagrees, so the gate is load bearing. + rhs_ungated = At * Diagonal(w .* net.sw) * (sol.lam_ub - sol.lam_lb) + + At * (sol.gamma_ub - sol.gamma_lb) + ungated = zeros(net.n) + ungated[non_ref] = sol.B_r_factor \ rhs_ungated[non_ref] + @test !isapprox(cong, ungated; atol=1e-3) + end end From 744680a37907f7857c18840a9638b5131c1652cd Mon Sep 17 00:00:00 2001 From: ckhanpour3 Date: Mon, 15 Jun 2026 17:05:25 -0400 Subject: [PATCH 10/13] Document DC topology cache and refresh Explain that the cache's from_bus/to_bus are the branch endpoints already encoded by the incidence matrix A, materialized as dense Int vectors so the connectivity refresh reads them in O(1). Annotate _refresh_topology_cache!: union-find over energized branches builds the island partition, then one reference bus is chosen per island (lowest index, except the configured ref_bus which is forced for its own island). --- src/types/dc_network.jl | 32 +++++++++++++++++++++++++++----- 1 file changed, 27 insertions(+), 5 deletions(-) diff --git a/src/types/dc_network.jl b/src/types/dc_network.jl index 494dbaa..b7e5bfa 100644 --- a/src/types/dc_network.jl +++ b/src/types/dc_network.jl @@ -29,11 +29,16 @@ The cache is prewarmed by constructors and refreshed by topology readers when mutations and the first topology read after each mutation. """ mutable struct _DCTopologyCache + # Branch endpoints as sequential bus indices, cached from the incidence matrix: + # from_bus[e]/to_bus[e] are the columns of the +1/-1 entries in row `e` of `A`. + # This is the same information `A` already encodes; it is materialized as dense + # Int vectors (once, since `A`'s structure is fixed) so the connectivity refresh + # can read each branch's endpoints in O(1) instead of rescanning sparse rows. from_bus::Vector{Int} to_bus::Vector{Int} - energized::BitVector - refs::Vector{Int} - non_ref::Vector{Int} + energized::BitVector # energized[e] == (b[e] * sw[e] != 0) + refs::Vector{Int} # one reference bus per energized island (sorted) + non_ref::Vector{Int} # all buses except `refs` initialized::Bool end @@ -669,8 +674,19 @@ function _topology_cache_valid(net::DCNetwork) return true end +# Recompute the energized-island partition and its reference buses. +# +# Energized branches (b[e]*sw[e] != 0) define the connectivity graph; buses joined +# by them form one island. We run union-find over those branches, then pick one +# reference bus per island: the lowest bus index in each, except the configured +# `ref_bus`, which is forced to be the reference for its own island. Buses with no +# energized branch are singleton islands that reference themselves. function _refresh_topology_cache!(net::DCNetwork) cache = net.topology_cache + + # Union-find over buses. `parent[i]` points toward i's component root; a bus is + # its own parent until merged. find_root uses path halving, union_roots! keeps + # the lower index as root so the partition (and thus refs) is deterministic. parent = collect(1:net.n) function find_root(i::Int) @@ -693,19 +709,23 @@ function _refresh_topology_cache!(net::DCNetwork) return nothing end + # Cache the branch endpoints from `A` once. `A`'s structure is fixed after + # construction, so this only runs the first time (or if `m` somehow changed). if length(cache.from_bus) != net.m cache.from_bus = zeros(Int, net.m) cache.to_bus = zeros(Int, net.m) rows, cols, nz_values = findnz(net.A) @inbounds for p in eachindex(rows) if nz_values[p] > 0 - cache.from_bus[rows[p]] = cols[p] + cache.from_bus[rows[p]] = cols[p] # +1 entry -> from bus else - cache.to_bus[rows[p]] = cols[p] + cache.to_bus[rows[p]] = cols[p] # -1 entry -> to bus end end end + # Recompute the energized flags (these track `b`/`sw`) and union the endpoints + # of every energized branch so each island collapses to a single root. resize!(cache.energized, net.m) @inbounds for e in 1:net.m cache.energized[e] = _is_energized(net, e) @@ -715,11 +735,13 @@ function _refresh_topology_cache!(net::DCNetwork) union_roots!(cache.from_bus[e], cache.to_bus[e]) end + # One reference per island: start with the lowest bus index in each component... refs_by_root = Dict{Int,Int}() @inbounds for bus in 1:net.n root = find_root(bus) refs_by_root[root] = min(get(refs_by_root, root, bus), bus) end + # ...then force the configured ref_bus to be the reference for its island. refs_by_root[find_root(net.ref_bus)] = net.ref_bus cache.refs = sort!(collect(values(refs_by_root))) cache.non_ref = setdiff(1:net.n, cache.refs) From ffcefaeca39f5c6d62824fb9fada64f9c7e806dd Mon Sep 17 00:00:00 2001 From: ckhanpour3 Date: Mon, 15 Jun 2026 17:21:07 -0400 Subject: [PATCH 11/13] Allow opt-in fractional switching in to_apf_network Replace the hard rejection of non-binary sw with an allow_fractional kwarg (default false). Fractional switching states are the basis of switching sensitivity, so rejecting them outright was too strict; the only constraint is that APF cannot represent partial switching. With allow_fractional=true the network is converted under APF's own sw > 0.5 binarization. Base the single-island requirement on that same binarization rule (closed and nonzero susceptance) rather than PD's b*sw != 0 energized check, so the slack-bus assumption matches the topology APF builds even for fractional sw. For binary sw the two rules coincide, so the default path is unchanged. Add a union-find connectivity helper and tests for both the opt-in conversion and the still-rejected disconnected-under-binarization case. --- ext/PowerDiffAPFExt.jl | 88 +++++++++++++++++++++++++----------- test/test_apf_integration.jl | 24 ++++++++-- 2 files changed, 82 insertions(+), 30 deletions(-) diff --git a/ext/PowerDiffAPFExt.jl b/ext/PowerDiffAPFExt.jl index a888842..d5c200e 100644 --- a/ext/PowerDiffAPFExt.jl +++ b/ext/PowerDiffAPFExt.jl @@ -53,43 +53,66 @@ end # ----------------------------------------------------------------------------- """ - to_apf_network(net::DCNetwork) → APF.Network + to_apf_network(net::DCNetwork; allow_fractional=false) → APF.Network Convert a `DCNetwork` to an `APF.Network`. APF networks lack generators, costs, and limits, so this is one way. -APF exposes one slack bus, so disconnected `DCNetwork`s are rejected. +APF exposes one slack bus, so the converted topology must be a single island. Bus demand is set to zero (PD separates demand from network topology). -Branch `status` is derived from switching state: `sw[e] > 0.5`. Because APF -models branches as strictly on/off, switching states must be binary (`sw[e] ∈ -{0, 1}`); fractional values are rejected so that PD's energized-island count -(`b[e] * sw[e] != 0`) cannot disagree with the converted APF topology. +Branch `status` is derived from switching state: `sw[e] > 0.5`, and the +susceptance `b[e]` is passed unscaled. + +Because APF models branches as strictly on/off, a fractional `sw[e]` has no +faithful APF representation. By default (`allow_fractional=false`) such switching +states are rejected. Switching sensitivity, however, treats `sw ∈ [0, 1]` as +continuous, so pass `allow_fractional=true` to convert anyway under APF's own +`sw[e] > 0.5` binarization. In either case the single-island requirement is +checked with that same binarization rule, so the slack-bus assumption matches +the topology APF actually builds. Note: `to_apf_network` sets bus demand to zero because PD separates demand from topology. For APF workflows that need demand data (e.g., `compute_flow!`), use `APF.from_power_models(pm_data)` directly instead. """ -function PowerDiff.to_apf_network(net::PowerDiff.DCNetwork) - n, m = net.n, net.m - - # APF binarizes branch status via sw[e] > 0.5, but PD treats a branch as - # energized whenever b[e] * sw[e] != 0. These two views only agree when every - # switching state is binary, so reject fractional sw to keep the island count - # used for the slack-bus check consistent with the topology APF actually builds. - all(s -> s == 0 || s == 1, net.sw) || throw(ArgumentError( - "AcceleratedDCPowerFlows conversion requires binary switching states (sw[e] ∈ {0, 1}); " * - "APF models branches as on/off via sw[e] > 0.5, so fractional switching would make the " * - "energized-island count disagree with the converted APF topology" - )) +# Single-island check under APF's branch-status rule. APF includes a branch in +# its susceptance matrix iff it is closed (`sw[e] > 0.5`) and carries nonzero +# susceptance, so islands must be counted with that exact rule (not PD's +# `b[e]*sw[e] != 0`) to match the topology APF builds. For binary sw the two +# rules coincide. Returns true iff all buses fall in one connected component. +function _apf_single_island(n::Int, from_bus::Vector{Int}, to_bus::Vector{Int}, active::AbstractVector{Bool}) + parent = collect(1:n) + function find(i::Int) + while parent[i] != i + parent[i] = parent[parent[i]] + i = parent[i] + end + return i + end + @inbounds for e in eachindex(active) + active[e] || continue + ri, rj = find(from_bus[e]), find(to_bus[e]) + ri == rj && continue + ri < rj ? (parent[rj] = ri) : (parent[ri] = rj) + end + n <= 1 && return true + r0 = find(1) + return all(find(i) == r0 for i in 2:n) +end - refs = PowerDiff.reference_buses(net) - length(refs) == 1 || throw(ArgumentError( - "AcceleratedDCPowerFlows conversion requires one energized island; found $(length(refs))" - )) +function PowerDiff.to_apf_network(net::PowerDiff.DCNetwork; allow_fractional::Bool=false) + n, m = net.n, net.m - # All buses are active: DCNetwork is built from PM.build_ref() which filters - # out inactive buses, so every bus in net.n is active by construction. - buses = [APF.Bus(i, true, 0.0) for i in 1:n] + # APF models branches as strictly on/off and passes b[e] unscaled, so a + # fractional sw has no faithful APF representation. Reject it unless the caller + # opts in, in which case we convert under APF's own sw[e] > 0.5 binarization. + if !allow_fractional + all(s -> s == 0 || s == 1, net.sw) || throw(ArgumentError( + "AcceleratedDCPowerFlows conversion requires binary switching states (sw[e] ∈ {0, 1}); " * + "APF binarizes status via sw[e] > 0.5 and uses b[e] unscaled, so fractional switching has " * + "no faithful APF representation. Pass allow_fractional=true to convert under the sw > 0.5 rule." + )) + end from_bus = zeros(Int, m) to_bus = zeros(Int, m) @@ -103,12 +126,25 @@ function PowerDiff.to_apf_network(net::PowerDiff.DCNetwork) end @assert all(>(0), from_bus) && all(>(0), to_bus) "Incidence matrix A must have exactly one +1 and one -1 per row" + # Count islands with the rule APF will apply (closed and nonzero susceptance), + # so the single-slack requirement matches the converted topology even when sw + # is fractional. For binary sw this matches PD's energized check exactly. + active = [net.sw[e] > 0.5 && !iszero(net.b[e]) for e in 1:m] + _apf_single_island(n, from_bus, to_bus, active) || throw(ArgumentError( + "AcceleratedDCPowerFlows conversion requires one energized island under APF's sw[e] > 0.5 " * + "binarization; the converted topology would be disconnected" + )) + + # All buses are active: DCNetwork is built from PM.build_ref() which filters + # out inactive buses, so every bus in net.n is active by construction. + buses = [APF.Bus(i, true, 0.0) for i in 1:n] + branches = Vector{APF.Branch}(undef, m) for e in 1:m branches[e] = APF.Branch(e, net.sw[e] > 0.5, net.b[e], net.fmax[e], from_bus[e], to_bus[e]) end - return APF.Network("PowerDiff", buses, only(refs), branches) + return APF.Network("PowerDiff", buses, net.ref_bus, branches) end # ----------------------------------------------------------------------------- diff --git a/test/test_apf_integration.jl b/test/test_apf_integration.jl index 3198e32..8816fc9 100644 --- a/test/test_apf_integration.jl +++ b/test/test_apf_integration.jl @@ -176,21 +176,37 @@ end @test_throws ArgumentError to_apf_network(net) end -@testset "to_apf_network rejects fractional switching" begin +@testset "to_apf_network fractional switching handling" begin net = load_test_case("case14.m") if isnothing(net) @test_skip false else dc_net = DCNetwork(net) - # Fractional sw cannot be represented by APF's on/off branch status and - # would make PD's energized-island count disagree with the APF topology. + # Default: fractional sw is rejected (APF cannot represent partial switching). dc_net.sw[1] = 0.3 @test_throws ArgumentError to_apf_network(dc_net) - dc_net.sw[1] = 1.0 # restore; binary states convert fine + # Opt-in: allow_fractional converts under APF's sw > 0.5 binarization. + # case14 stays connected with branch 1 binarized open, so it is one island. + apf_frac = to_apf_network(dc_net; allow_fractional=true) + @test apf_frac isa APF.Network + @test !apf_frac.branches[1].status # sw = 0.3 binarizes to open + dc_net.sw[1] = 1.0 # restore; binary states convert under the default path @test to_apf_network(dc_net) isa APF.Network end end +@testset "to_apf_network rejects disconnected binarized topology" begin + # Single branch between two buses. allow_fractional still checks connectivity + # under sw > 0.5: at sw = 0.3 the branch binarizes open, leaving two islands. + A = sparse([1.0 -1.0]) + G = sparse(reshape([1.0, 0.0], 2, 1)) + net_open = DCNetwork(2, 1, 1, A, G, [-10.0]; sw=[0.3], ref_bus=1) + @test_throws ArgumentError to_apf_network(net_open; allow_fractional=true) + # At sw = 0.7 the branch binarizes closed, so the two buses form one island. + net_closed = DCNetwork(2, 1, 1, A, G, [-10.0]; sw=[0.7], ref_bus=1) + @test to_apf_network(net_closed; allow_fractional=true) isa APF.Network +end + # ========================================================================= # Index alignment # ========================================================================= From f5229bf3680fcb7de7054719ad5457316e90e8aa Mon Sep 17 00:00:00 2001 From: Samuel Talkington <10187005+samtalki@users.noreply.github.com> Date: Mon, 15 Jun 2026 19:27:02 -0400 Subject: [PATCH 12/13] Fix wording in advanced.md for topology_cache Corrected wording in the description of 'topology_cache'. --- docs/src/advanced.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/advanced.md b/docs/src/advanced.md index 531f6b5..04fdef2 100644 --- a/docs/src/advanced.md +++ b/docs/src/advanced.md @@ -42,7 +42,7 @@ Stores the DC network topology and parameters. | `ref_bus` | `Int` | Preferred reference bus index (sequential) | | `tau` | `Float64` | Regularization parameter | | `id_map` | `IDMapping` | Bidirectional element ID mapping (original ↔ sequential) | -| `topology_cache` | `_DCTopologyCache` | Internal energized-island cache (not part of the public API) | +| `topology_cache` | `_DCTopologyCache` | Internal energized island cache (not part of the public API) | Construct from a parsed MATPOWER network with `DCNetwork(parse_file("case14.m"))`, or with explicit parameters: `DCNetwork(n, m, k, A, G_inc, b; ...)`. From 915caaac535ff03694888a6873925cbd221e528b Mon Sep 17 00:00:00 2001 From: ckhanpour3 Date: Mon, 15 Jun 2026 19:43:22 -0400 Subject: [PATCH 13/13] Gate angle duals consistently in LMP decomposition docs Address Copilot review: the file-header decomposition comment in lmp.jl, the inline comment in calc_congestion_component, and the LMP decomposition equation in dc-opf.md still showed the angle-difference dual term ungated. Add the Diag(sw) gate so they match the implementation and the theta-stationarity, and reword the inline note to cover fully closed, fractional, and open branches rather than implying energized always means sw == 1. --- docs/src/math/dc-opf.md | 2 +- src/sens/lmp.jl | 9 +++++---- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/docs/src/math/dc-opf.md b/docs/src/math/dc-opf.md index 145e0ca..3170f49 100644 --- a/docs/src/math/dc-opf.md +++ b/docs/src/math/dc-opf.md @@ -148,7 +148,7 @@ Locational marginal prices are the power balance duals ``\nu_{\text{bal}}``, dec The congestion component is extracted by solving: ```math -\text{congestion}[\text{non-ref}] = B_r^{-1} \left(A_r^\top W (\lambda_{\text{ub}}^{\text{std}} - \lambda_{\text{lb}}^{\text{std}}) + A_r^\top (\gamma_{\text{ub}}^{\text{std}} - \gamma_{\text{lb}}^{\text{std}})\right) +\text{congestion}[\text{non-ref}] = B_r^{-1} \left(A_r^\top W (\lambda_{\text{ub}}^{\text{std}} - \lambda_{\text{lb}}^{\text{std}}) + A_r^\top \operatorname{diag}(\mathrm{sw}) (\gamma_{\text{ub}}^{\text{std}} - \gamma_{\text{lb}}^{\text{std}})\right) ``` where ``\lambda^{\text{std}}``, ``\gamma^{\text{std}}`` use the standard sign convention (non-negative for binding constraints). The ``\gamma`` terms capture congestion from binding phase angle difference limits. The energy component is uniform within each energized island and reflects that island's marginal cost of generation. diff --git a/src/sens/lmp.jl b/src/sens/lmp.jl index 8b51248..2d447d5 100644 --- a/src/sens/lmp.jl +++ b/src/sens/lmp.jl @@ -25,7 +25,7 @@ # LMP Decomposition (for analysis): # LMP = ν_bal = energy_component + congestion_component # where: -# congestion_component = B_r⁻¹ [A_r' Diag(-b .* sw) (λ_ub - λ_lb) + A_r'(γ_ub - γ_lb)] (non-ref block) +# congestion_component = B_r⁻¹ [A_r' Diag(-b .* sw) (λ_ub - λ_lb) + A_r' Diag(sw) (γ_ub - γ_lb)] (non-ref block) # energy_component = ν_bal - congestion_component (uniform within each island) # # Sign conventions (DC OPF): @@ -99,9 +99,10 @@ function calc_congestion_component(sol::DCOPFSolution, net::DCNetwork; At = net.A' # Angle difference constraints are gated by `sw` in the model, so their - # stationarity contribution carries the same `Diag(sw)` factor. For energized - # branches `sw == 1` (a no-op); on de-energized branches the gate matches the - # zero angle dual the solver returns there. + # stationarity contribution carries the same `Diag(sw)` factor: the gate is the + # identity on fully closed branches (sw == 1), scales the contribution on + # fractional branches (0 < sw < 1), and zeroes the term on open branches + # (sw == 0), matching the gated angle dual the solver returns. rhs_full = At * Diagonal(w .* net.sw) * (sol.lam_ub - sol.lam_lb) + At * Diagonal(net.sw) * (sol.gamma_ub - sol.gamma_lb)