Support disconnected DC network islands

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
 

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

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)

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.

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}

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

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

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.