diff --git a/docs/src/advanced.md b/docs/src/advanced.md index 3545400..04fdef2 100644 --- a/docs/src/advanced.md +++ b/docs/src/advanced.md @@ -36,14 +36,31 @@ 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 | -| `ref_bus` | `Int` | Reference bus index (sequential) | +| `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. 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 +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. @@ -73,7 +90,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/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/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 084e484..3170f49 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,13 @@ 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 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 +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 @@ -51,13 +58,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 \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`` @@ -66,7 +73,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 @@ -78,7 +85,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 +95,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 nonsmooth; the implementation uses the fixed zero shedding convention already required by the collapsed bound ``0 \leq \text{psh} \leq 0``. @@ -102,6 +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 +or closing a bridge splits or merges islands and changes the reference set, so +the derivative is nonsmooth at that topology boundary. ### Cost Coefficients (``c_q``, ``c_l``) @@ -132,13 +142,13 @@ 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: ```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 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..ce2ce13 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: @@ -23,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]``: @@ -38,6 +45,11 @@ 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. 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 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 +66,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..d5c200e 100644 --- a/ext/PowerDiffAPFExt.jl +++ b/ext/PowerDiffAPFExt.jl @@ -53,24 +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 networks lack generators, costs, and limits, so this is one way. +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`. +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) +# 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 + +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) @@ -84,6 +126,19 @@ 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]) 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..68a1ec0 100644 --- a/src/prob/kkt_dc_opf.jl +++ b/src/prob/kkt_dc_opf.jl @@ -281,19 +281,25 @@ 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 +function kkt_dims(prob::DCOPFProblem) + net = getfield(prob, :network) + return _dc_kkt_dims(getfield(net, :n), getfield(net, :m), getfield(net, :k), + 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 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 +309,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 +332,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,8 +344,30 @@ function kkt_indices(n::Int, m::Int, k::Int) ) end -kkt_indices(net::DCNetwork) = kkt_indices(net.n, net.m, net.k) -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), + getfield(prob, :_n_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 = getfield(prob, :_n_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 @@ -351,8 +381,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 +400,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 +451,13 @@ 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' * 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 @@ -460,13 +491,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 +542,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 +613,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 +626,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 +645,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 +692,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 +836,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 +859,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(n, m, k) - idx = kkt_indices(n, m, k) + dim, idx = _dc_kkt_layout(net) mu_ub = getfield(sol, :mu_ub) colptr = Vector{Int}(undef, n + 1) @@ -853,8 +887,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 @@ -888,7 +922,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 @@ -899,9 +933,6 @@ function calc_kkt_jacobian_switching(prob::DCOPFProblem, sol::DCOPFSolution) J_s = spzeros(dim, m) - # Use centralized index calculation - idx = kkt_indices(n, m, k) - # Precompute A * θ once (invariant across branches) Aθ = A * θ @@ -922,7 +953,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 @@ -955,8 +986,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) @@ -970,11 +1001,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 @@ -995,7 +1026,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/cost.jl b/src/sens/cost.jl index 571d295..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(n, m, k) - idx = kkt_indices(n, m, k) + 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(n, m, k) - idx = kkt_indices(n, m, k) + 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 87438c2..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(n, m, k) - idx = kkt_indices(n, m, k) + 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/lmp.jl b/src/sens/lmp.jl index 3162e63..2d447d5 100644 --- a/src/sens/lmp.jl +++ b/src/sens/lmp.jl @@ -25,8 +25,8 @@ # 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) -# energy_component = ν_bal - congestion_component (uniform for connected network) +# 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): # - Our LMPs are positive (cost increases when demand increases) @@ -76,14 +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 the 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. @@ -91,10 +95,16 @@ 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) + # Angle difference constraints are gated by `sw` in the model, so their + # 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) result = zeros(net.n) result[non_ref] = B_r_factor \ rhs_full[non_ref] @@ -106,8 +116,9 @@ 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 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. @@ -176,4 +187,3 @@ function calc_qlmp(prob::ACOPFProblem) sol = _ensure_ac_solved!(prob) return calc_qlmp(sol, prob) end - diff --git a/src/sens/susceptance.jl b/src/sens/susceptance.jl index 37d57c1..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 @@ -50,8 +51,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(n, m, k) + dim, idx = _dc_kkt_layout(prob) theta = sol.va nu_bal = sol.nu_bal @@ -91,8 +91,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/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..b7e5bfa 100644 --- a/src/types/dc_network.jl +++ b/src/types/dc_network.jl @@ -19,6 +19,31 @@ # 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. + +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 + # 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 # 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 + +_DCTopologyCache() = _DCTopologyCache(Int[], Int[], BitVector(), Int[], Int[], false) + """ DCNetwork <: AbstractPowerNetwork @@ -35,12 +60,15 @@ 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) -- `ref_bus`: Reference bus index (phase angle = 0) +- `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 - `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 @@ -63,6 +91,7 @@ struct DCNetwork <: AbstractPowerNetwork ref_bus::Int tau::Float64 id_map::IDMapping + topology_cache::_DCTopologyCache end # ============================================================================= @@ -85,7 +114,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 +133,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 +145,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 @@ -164,7 +193,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 @@ -274,7 +303,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 @@ -404,7 +433,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) @@ -508,7 +537,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. @@ -535,8 +564,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) + 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 """ @@ -566,7 +598,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), @@ -576,8 +608,11 @@ function DCNetwork( Float64.(c_shed), Float64.(demand), Float64.(pg_init), ref_bus, tau, - IDMapping(n, m, k) + IDMapping(n, m, k), + _DCTopologyCache() ) + _refresh_topology_cache!(net) + return net end # ============================================================================= @@ -624,24 +659,140 @@ function calc_susceptance_matrix(network::DCNetwork) return sparse(network.A' * W * network.A) end +@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 + +# 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) + 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 + + # 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] # +1 entry -> from bus + else + 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) + 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!(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) + cache.initialized = true + return cache +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 + +_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) = copy(_topology_cache(net).non_ref) + """ _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 +821,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 +852,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..3de0be0 100644 --- a/src/types/dc_opf_problem.jl +++ b/src/types/dc_opf_problem.jl @@ -50,8 +50,8 @@ 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) -- `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) +- `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 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 @@ -122,6 +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 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) """ @@ -135,6 +139,7 @@ mutable struct DCOPFProblem{O} <: AbstractOPFProblem d::Vector{Float64} cons::NamedTuple cache::DCSensitivityCache + _n_ref::Int _optimizer::O _silent::Bool end @@ -175,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(), 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 @@ -192,6 +207,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 @@ -244,8 +265,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) @@ -256,6 +278,9 @@ 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, flow_def = flow_def, diff --git a/src/types/show.jl b/src/types/show.jl index ad1bd71..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 # ============================================================================= @@ -76,11 +85,11 @@ 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))]") - 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/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_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 diff --git a/test/test_apf_integration.jl b/test/test_apf_integration.jl index d00c75d..8816fc9 100644 --- a/test/test_apf_integration.jl +++ b/test/test_apf_integration.jl @@ -171,6 +171,42 @@ 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 + +@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) + # Default: fractional sw is rejected (APF cannot represent partial switching). + dc_net.sw[1] = 0.3 + @test_throws ArgumentError to_apf_network(dc_net) + # 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 # ========================================================================= diff --git a/test/test_dc_islands.jl b/test/test_dc_islands.jl new file mode 100644 index 0000000..e17483b --- /dev/null +++ b/test/test_dc_islands.jl @@ -0,0 +1,191 @@ +# 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] + @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 + 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 "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] + + # 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 ["case6.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 + + 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(parsed) + d = calc_demand_vector(parsed) + @test getfield(net, :topology_cache).initialized + @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 "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