Handle negative net demand in DC load shedding

[cameronkhanpour]

Summary

Addresses #53.

This PR fixes the IEEE300 DC OPF failure caused by negative net demand interacting with load-shedding bounds.

Root Cause

DC OPF previously imposed 0 <= psh <= d. IEEE300 contains buses with negative net demand, so those buses received contradictory shedding bounds: psh >= 0 and psh <= d < 0.

Fix

• Define the curtailable load at each bus as the positive part of signed net demand: d_plus = max(d, 0).
• Keep the original signed d in nodal power balance, so negative net demand remains an injection.
• Fix psh = 0 where d_plus = 0, because injections cannot be shed.
• Apply the same rule in DC OPF construction, update_demand!, solve post-processing, KKT residuals, analytical Jacobians, and allocation-light JVP/VJP paths.
• Drop negative-demand reporting from @warn to @debug now that embedded generation / negative net demand is supported input.
• Document the positive-part shedding-cap design choice and its non-smooth point at d = 0.
• Add a regression test covering negative net demand solves, KKT residuals, finite-difference sensitivity behavior, JVP/VJP consistency, and update_demand!.

AC OPF does not need the same code change: it has no load-shedding variable or shedding upper bound and already keeps signed load directly in nodal balance.

Follow-Up

The broader disconnected DC topology / multi-island work that was previously bundled here has been moved to draft follow-up #55. That PR is stacked on this one so the issue #53 fix can be reviewed and merged independently.

Validation

• Clean temporary Julia environment: Pkg.test("PowerDiff")
• Clean temporary docs environment: SITE_BUILD=true include("docs/make.jl")
• GitHub checks on the rewritten PR are green: CI test, Documentation build, and Benchmark benchmark

[github-actions]

Benchmark Results (Julia v1)

Time benchmarks

	main	fb62986...	main / fb62986...
ac_opf/kkt_jacobian/case30.m	2.44 ± 0.031 ms	2.48 ± 0.035 ms	0.982 ± 0.019
ac_opf/kkt_param/case30.m/switching	0.0839 ± 0.0058 ms	0.0856 ± 0.0054 ms	0.98 ± 0.092
dc_opf/kkt_jacobian/case30.m/cost_linear	0.131 ± 0.01 μs	0.121 ± 0.011 μs	1.08 ± 0.13
dc_opf/kkt_jacobian/case30.m/cost_quadratic	0.121 ± 0.01 μs	0.14 ± 0.02 μs	0.864 ± 0.14
dc_opf/kkt_jacobian/case30.m/demand	0.22 ± 0.03 μs	0.24 ± 0.099 μs	0.917 ± 0.4
dc_opf/kkt_jacobian/case30.m/flowlimit	0.311 ± 0.049 μs	0.331 ± 0.13 μs	0.94 ± 0.39
dc_opf/kkt_jacobian/case30.m/full	12.8 ± 4.7 μs	12.1 ± 6.4 μs	1.06 ± 0.68
dc_opf/kkt_jacobian/case30.m/susceptance	0.0723 ± 0.0043 ms	0.0768 ± 0.0038 ms	0.941 ± 0.073
time_to_load	1.42 ± 0.0023 s	1.43 ± 0.0063 s	0.995 ± 0.0046

Memory benchmarks

	main	fb62986...	main / fb62986...
ac_opf/kkt_jacobian/case30.m	0.0339 M allocs: 1.18 MB	0.0339 M allocs: 1.18 MB	1
ac_opf/kkt_param/case30.m/switching	1.49 k allocs: 0.602 MB	1.49 k allocs: 0.602 MB	1
dc_opf/kkt_jacobian/case30.m/cost_linear	6 allocs: 0.328 kB	6 allocs: 0.328 kB	1
dc_opf/kkt_jacobian/case30.m/cost_quadratic	6 allocs: 0.328 kB	6 allocs: 0.328 kB	1
dc_opf/kkt_jacobian/case30.m/demand	6 allocs: 1.42 kB	6 allocs: 1.42 kB	1
dc_opf/kkt_jacobian/case30.m/flowlimit	6 allocs: 1.89 kB	6 allocs: 1.89 kB	1
dc_opf/kkt_jacobian/case30.m/full	0.079 k allocs: 0.082 MB	0.079 k allocs: 0.082 MB	1
dc_opf/kkt_jacobian/case30.m/susceptance	2.28 k allocs: 0.292 MB	2.28 k allocs: 0.292 MB	1
time_to_load	0.149 k allocs: 11.1 kB	0.149 k allocs: 11.1 kB	1

[cameronkhanpour]

Follow-up DC island support is pushed in 18f81f5.

Multi-island root cause

The native DC stack assumed a single ref_bus. For networks with multiple energized islands or isolated buses, removing one Laplacian row and column left additional island nullspaces in OPF, power flow, decomposition, and sensitivity paths.

Fix

• Add reference_buses(net), preserving configured net.ref_bus for its island and deterministically adding the lowest-index bus for every other energized island, including isolated buses.
• Use the component-aware reduction in DCPF, DC OPF, PTDF-related sensitivity paths, LMP decomposition, KKT residuals, analytical Jacobians, JVPs, and VJPs.
• Generalize DC KKT reference duals from one scalar to one value per energized island.
• Clear topology-dependent JVP/VJP workspaces when switching changes can resize the KKT system.
• Document per-island references, per-island energy components, and the non-smooth sensitivity boundary when a switch splits or merges islands.
• Reject multi-island AcceleratedDCPowerFlows conversion clearly because APF exposes one slack bus.

AC OPF is intentionally unchanged.

Committed regressions

• Programmatic two-island, isolated-load-shedding, fully isolated bus, and bridge-opening cases.
• OPF, DCPF, KKT residual, LMP decomposition, cache invalidation, and finite-sensitivity checks.
• Bundled PowerModels regressions for case5_db.m and case6.m.
• Separate classification for case7_tplgy.m, which includes unsupported specialized components.

Validation

• Pkg.test(PowerDiff)
• Focused multi-island suite: 46/46 passing
• Optional APF extension runtime suite could not run locally: AcceleratedDCPowerFlows is not installed and is not resolvable from the configured Julia registry.

Disposable PGLib-OPF v23.07 sanity sweep

Ran an external, uncommitted sweep over all 198 MATPOWER files: 66 base, 66 api, and 66 sad. The checkout, runner, temp environments, and reports were deleted afterward.

• 112 supported cases solved with PowerDiff and matched PowerModels generation cost, including MATPOWER generator constant terms.
• 8 cases exercised intentional PowerDiff load-shedding behavior.
• 24 cases were classified separately for active shunts outside the native DC formulation.
• 42 cases did not solve in PowerModels under the sweep limits: 38 locally infeasible and 4 time-limited.
• 12 large PowerDiff runs reached the package's 30-second Ipopt cap. These were solver-limited rows, not disconnected-topology factorization failures; a smaller example remained time-limited after a 180-second probe.
• 0 supported solved cases had an unexplained PowerDiff-vs-PowerModels cost mismatch after classification.

The issue #53 and disconnected-topology failures are addressed in this draft. The sweep leaves separate future work for large Ipopt-limited networks and currently unsupported PowerModels features such as active shunts.

[cameronkhanpour]

Pushed follow-up commit 60d7980 to fix the Documentation failure and the DC benchmark regression introduced by multi-island support.

Root cause

The multi-island change made each lightweight DC KKT parameter-Jacobian builder call reference_buses(net) through both kkt_dims(net) and kkt_indices(net). reference_buses rebuilt graph components from the sparse incidence matrix each time. For tiny builders such as linear-cost and demand Jacobians, that topology work dominated runtime and increased allocations from a few hundred bytes to about 65 KB per call.

Fix

• Added an internal energized-topology cache to DCNetwork with cached branch endpoints, energized mask, effective reference buses, and non-reference buses.
• Cache hits validate only whether an edge crossed the energized/non-energized boundary; graph components rebuild only after a real split or merge.
• Direct in-place sw mutations and zero-crossing b mutations are still detected.
• Public reference_buses(net) returns a defensive copy, while native hot paths reuse the cached layout.
• KKT builders now recover dimensions and indices from one layout lookup. Builders with a DCOPFProblem use the already-built reference-constraint count directly.
• Added reference_buses to the canonical API docs block, fixing Documenter's checkdocs = :exports failure.
• Added a regression covering defensive copies plus direct sw split/merge and b zero-crossing cache refreshes.

Local verification

• Pkg.test(PowerDiff) passes.
• Local Documenter build passes.
• APF extension package was not installed locally, so its optional extension-specific path could not be run; the regular APF integration fallback tests passed.
• Quick case30.m allocation checks are back near baseline: cost_linear 640 B, cost_quadratic 432 B, demand 1.8 KB, and flowlimit 2.0 KB per call, instead of roughly 65 KB each in the regressed benchmark.

[cameronkhanpour]

CI follow-up for 60d7980: all relevant checks are green.

• CI: passed
• Documentation: passed
• Benchmark: passed

The GitHub benchmark confirms the accidental topology-scan regression is removed. Representative case30.m results:

DC KKT builder	Regressed 18f81f5	Fixed 60d7980	main
cost_linear	39.1 μs, 477 allocs / 64.8 KB	0.261 μs, 6 allocs / 328 B	0.170 μs, 6 allocs / 328 B
cost_quadratic	39.6 μs, 477 allocs / 64.8 KB	0.260 μs, 7 allocs / 375 B	0.151 μs, 6 allocs / 328 B
demand	40.1 μs, 477 allocs / 65.9 KB	0.421 μs, 6 allocs / 1.42 KB	0.290 μs, 6 allocs / 1.42 KB
flowlimit	15.7 μs, 477 allocs / 66.3 KB	0.501 μs, 7 allocs / 1.94 KB	0.421 μs, 6 allocs / 1.89 KB
susceptance	106 μs, 2.75k allocs / 357 KB	100 μs, 2.28k allocs / 292 KB	98.8 μs, 2.28k allocs / 292 KB

The small residual sub-microsecond overhead is the energized-mask validation that preserves correctness for direct in-place topology mutations.

[cameronkhanpour]

Pushed 60a7400 as a small benchmark refinement.

The remaining table mixed two cases:

• cost_linear and demand receive a bare DCNetwork, so they intentionally validate the energized-edge mask to remain correct after direct in-place sw or zero-crossing b mutations.
• cost_quadratic and flowlimit receive a built DCOPFProblem, but were still reading the reference count through the dynamically typed prob.cons tuple.

The follow-up stores the built model's effective reference count in a concrete internal prob._n_ref field and refreshes it whenever _rebuild_jump_model! runs. This removes the avoidable problem-backed lookup overhead while preserving the topology semantics.

Local verification:

• Pkg.test(PowerDiff) passes.
• CI-style local BenchmarkTools smoke run returns cost_quadratic and flowlimit to 6 allocations, matching the original allocation shape.

[cameronkhanpour]

CI benchmark follow-up for 60a7400: the avoidable problem-backed overhead is gone.

DC KKT builder	main	60a7400	Interpretation
cost_quadratic	0.121 μs	0.130 μs	matches within noise
flowlimit	0.331 μs	0.321 μs	matches within noise
susceptance	80.7 μs	80.1 μs	matches within noise
cost_linear	0.121 μs	0.201 μs	expected bare-network validation cost
demand	0.221 μs	0.310 μs	expected bare-network validation cost

All lightweight builders now match main allocations exactly. The remaining cost_linear and demand delta is about 0.08–0.09 μs per call and comes from allocation-free validation of the energized-edge mask. Those functions receive only a mutable DCNetwork, so this check is what keeps KKT dimensions correct if callers mutate net.sw directly or move net.b across zero.

Removing that last delta would require a broader API decision: tracked mutation-aware vectors, or requiring all topology changes to go through an invalidating update function. I do not think that complexity or behavioral restriction is justified for a sub-microsecond bare-network helper cost.

[samtalki]

Solid fix, and the test coverage here is genuinely good. A few things before it lands.

Scope: I'd split this. The actual #53 fix (psh ≤ max(d,0), keep signed d in power balance, fix psh=0 where d₊=0) is small, self-contained, and exactly what the issue asked for. The multi-island reference_buses work plus _DCTopologyCache and the _n_ref caching are a separate feature. They came out of the sweep that "surfaced separate edge cases," and they carry their own design decisions and two rounds of perf fixes for regressions the feature introduced. Bundling them means a one-line bug fix is now gated on reviewing a much larger, riskier change. I'd land the negative net demand fix on its own and put multi-island in a separate PR. If you'd rather keep them together that's fine, but say so in the title, since "Handle negative net demand" undersells what's actually in here.

Topology cache and threads. DCNetwork is an immutable struct, but topology_cache is a mutable struct that now gets written lazily on read. reference_buses, kkt_dims(net), kkt_indices(net), kkt(z, net, d), _factorize_B_r, all of them trigger _refresh_topology_cache! on first touch (and after any energized-edge change). No lock. So a DCNetwork that used to be safe to share across threads (truly immutable, read only) is now a data race if two threads first-touch it at the same time, or one reads while another flips sw/b. Either guard the refresh or document that a DCNetwork can no longer be shared across threads.

_n_ref is a second source of truth. kkt_dims(prob)/kkt_indices(prob) read the cached prob._n_ref, while kkt(z, net, d) and flatten_variables recompute _reference_buses(net) live. They agree as long as topology only changes through update_switching! (which rebuilds and refreshes _n_ref). Mutate prob.network.sw/b directly and they diverge: kkt_dims(prob) goes stale while the residual uses the fresh ref set, which is a confusing path to a dimension mismatch. It matches the existing "go through update_switching!" contract so it isn't wrong, but worth a comment pinning the invariant _n_ref == length(prob.cons.ref).

Minor:

• _warn_negative_demand still @warns now that negative net demand is a supported case rather than a failure. Anyone with legit negative net injection (embedded gen) gets warned on every build. Consider dropping it to @debug.
• kkt_indices(...).η changed from an Int to a UnitRange (length 1 in the single reference case). Fine inside the repo, just noting it for any external callers.

The negative demand path and the island math both have real KKT residual plus finite difference plus JVP/VJP coverage, which is the right way to do it. Nice work there.

[cameronkhanpour]

Implemented the requested scope split.

• Force-pushed Handle negative net demand in DC load shedding #54 down to the focused issue Issue with DC OPF solve (PowerDiff vs PowerModels) #53 fix: negative net demand now stays signed in power balance, load shedding is capped by max(d, 0), and psh is fixed at zero where there is no positive curtailable load.
• Removed the disconnected-topology / multi-island work from this PR, including reference_buses, _DCTopologyCache, and _n_ref changes, so the topology-cache/thread-safety and cached-invariant design questions no longer gate Handle negative net demand in DC load shedding #54.
• Lowered negative-demand reporting from @warn to @debug and updated the tests so supported negative net demand is quiet by default.
• Opened draft follow-up Support disconnected DC network islands #55 for the separated multi-island work, stacked on this branch.

Validation after the rewrite:

• clean temporary Julia environment: Pkg.test("PowerDiff")
• clean temporary docs environment: SITE_BUILD=true include("docs/make.jl")
• GitHub checks on Handle negative net demand in DC load shedding #54 are green: CI, Documentation, and Benchmark