From 61bd88337139b3d64e9d41e358d566b925a3533c Mon Sep 17 00:00:00 2001 From: Cameron Khanpour <99142483+cameronkhanpour@users.noreply.github.com> Date: Mon, 1 Jun 2026 23:37:18 -0400 Subject: [PATCH 1/2] Handle negative net demand in DC load shedding --- docs/src/math/dc-opf.md | 14 +++++++--- src/prob/dc_opf.jl | 28 ++++++++++---------- src/prob/kkt_dc_opf.jl | 24 +++++++++-------- src/sens/vjp_jvp.jl | 17 ++++++------ src/types/dc_opf_problem.jl | 15 ++++++++--- test/runtests.jl | 6 ++--- test/test_psh.jl | 53 +++++++++++++++++++++++++++++++++---- 7 files changed, 110 insertions(+), 47 deletions(-) diff --git a/docs/src/math/dc-opf.md b/docs/src/math/dc-opf.md index b1a958e..8940718 100644 --- a/docs/src/math/dc-opf.md +++ b/docs/src/math/dc-opf.md @@ -16,7 +16,7 @@ G_{\text{inc}} g + \text{psh} - d &= B \theta & (\nu_{\text{bal}}) \\ f &= W A \theta & (\nu_{\text{flow}}) \\ -f_{\max} \leq f &\leq f_{\max} & (\lambda_{\text{lb}}, \lambda_{\text{ub}}) \\ 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}}) \\ +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}}) \end{aligned} @@ -30,6 +30,7 @@ where: - ``C_q = \operatorname{diag}(c_q)`` contains quadratic cost coefficients - ``c_l`` contains linear cost coefficients - ``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 ## KKT System for Implicit Differentiation @@ -77,13 +78,20 @@ For each parameter type, we compute ``\partial K / \partial p`` and then the ful ### Demand (``d``) -Demand enters the power balance and the upper shedding bound: +Demand enters the power balance and the upper shedding bound through +``d_+ = \max(d, 0)``: ```math \frac{\partial K_{\nu_{\text{bal}}}}{\partial d} = -I, \qquad -\frac{\partial K_{\mu_{\text{ub}}}}{\partial d} = \operatorname{diag}(\mu_{\text{ub}}) +\frac{\partial K_{\mu_{\text{ub}}}}{\partial d} = +\operatorname{diag}\left(\mu_{\text{ub}} \circ \frac{\partial d_+}{\partial d}\right) ``` +For strictly positive demand, ``\partial d_+ / \partial d = 1``. For negative +demand, it is ``0``. At zero demand, the positive-part function is non-smooth; +the implementation uses the fixed-zero shedding convention already required by +the collapsed bound ``0 \leq \text{psh} \leq 0``. + ### Switching (``\mathrm{sw}``) Switching affects the Laplacian ``B``, weight matrix ``W``, and flow definition through: diff --git a/src/prob/dc_opf.jl b/src/prob/dc_opf.jl index e6dac22..525bb50 100644 --- a/src/prob/dc_opf.jl +++ b/src/prob/dc_opf.jl @@ -47,14 +47,15 @@ function _check_solve_status(model, label::String) end """ - _warn_negative_demand(d) + _debug_negative_demand(d) -Warn if any demand entries are negative, which makes shedding bounds infeasible. +Log if any demand entries are negative. Negative net demand is retained in +power balance, while load shedding is disabled at those buses. """ -function _warn_negative_demand(d) +function _debug_negative_demand(d) neg_buses = findall(d .< 0) if !isempty(neg_buses) - _SILENCE_WARNINGS[] || @warn "Negative demand at buses $neg_buses; shedding bounds 0 ≤ psh ≤ d will be infeasible at those buses" + @debug "Negative demand treated as net injection with psh fixed at zero" buses=neg_buses end end @@ -119,10 +120,8 @@ function solve!(prob::DCOPFProblem) # Snap to strict complementarity for clean KKT sensitivity computation. d = prob.d for i in eachindex(psh_val) - if d[i] < -TOL - _SILENCE_WARNINGS[] || @warn "Negative demand at bus $i (d=$(d[i])); psh snap may be unreliable" - end - if abs(d[i]) < TOL + shed_capacity = _shed_capacity(d[i]) + if shed_capacity < TOL # Degenerate: both bounds collapse to 0 ≤ psh ≤ 0 psh_val[i] = 0.0 # Canonicalize the duplicate dual pair so the sensitivity KKT sees @@ -133,9 +132,9 @@ function solve!(prob::DCOPFProblem) # Lower bound active: psh = 0 psh_val[i] = 0.0 μ_ub[i] = 0.0 - elseif d[i] - psh_val[i] < TOL - # Upper bound active: psh = d - psh_val[i] = d[i] + elseif shed_capacity - psh_val[i] < TOL + # Upper bound active: psh = max(d, 0) + psh_val[i] = shed_capacity μ_lb[i] = 0.0 else # Strictly interior: both bounds inactive @@ -165,7 +164,8 @@ end Rewrite the RHS of the power balance and shedding upper bound constraints to the new demand `d`, mutate `prob.d` in place, and invalidate the sensitivity -cache. Avoids rebuilding the JuMP model. +cache. The shedding upper bound is `max(d, 0)` because negative net demand is +an injection rather than curtailable load. Avoids rebuilding the JuMP model. To keep the JuMP constraints and the sensitivity cache consistent with the stored demand, always go through this function — do not mutate `prob.d` in @@ -175,13 +175,13 @@ function update_demand!(prob::DCOPFProblem, d::AbstractVector) n = prob.network.n length(d) == n || throw(DimensionMismatch("Demand vector length $(length(d)) must match number of buses $n")) - _warn_negative_demand(d) + _debug_negative_demand(d) invalidate!(prob.cache) prob.d .= d for i in 1:n set_normalized_rhs(prob.cons.power_bal[i], d[i]) - set_normalized_rhs(prob.cons.shed_ub[i], d[i]) + set_normalized_rhs(prob.cons.shed_ub[i], _shed_capacity(d[i])) end return prob diff --git a/src/prob/kkt_dc_opf.jl b/src/prob/kkt_dc_opf.jl index c8cd940..aaf07df 100644 --- a/src/prob/kkt_dc_opf.jl +++ b/src/prob/kkt_dc_opf.jl @@ -72,7 +72,8 @@ function _ensure_kkt_factor!(prob::DCOPFProblem) end # Must match solve!'s snap tolerance so the canonicalized duals agree with the KKT system. -@inline _is_fixed_zero_shed(d::Real) = abs(d) < COMPLEMENTARITY_SNAP_TOL +@inline _is_fixed_zero_shed(d::Real) = _shed_capacity(d) < COMPLEMENTARITY_SNAP_TOL +@inline _shed_capacity_derivative(d::Real) = _is_fixed_zero_shed(d) ? zero(d) : one(d) # ============================================================================= # Cached Derivative Computation Functions @@ -421,7 +422,7 @@ s.t. G_inc * g + psh - d = B * θ (ν_bal) g ≥ gmin (ρ_lb) g ≤ gmax (ρ_ub) 0 ≤ psh (μ_lb) - psh ≤ d (μ_ub) + psh ≤ max(d, 0) (μ_ub) θ[ref] = 0 (η_ref) ``` @@ -492,9 +493,9 @@ function kkt(z::AbstractVector, net::DCNetwork, d::AbstractVector) K_ρ_ub = ρ_ub .* (net.gmax - g) # 8. Load shedding bounds - # If d[i] == 0, psh[i] is structurally fixed to zero and the upper-bound - # dual is redundant. Replace the degenerate complementarity pair by - # psh[i] = 0 and μ_ub[i] = 0 to avoid an exactly singular KKT block. + # If max(d[i], 0) == 0, psh[i] is structurally fixed to zero and the + # upper-bound dual is redundant. Replace the degenerate complementarity + # pair by psh[i] = 0 and μ_ub[i] = 0 to avoid an exactly singular KKT block. K_μ_lb = similar(psh) K_μ_ub = similar(psh) @inbounds for i in 1:n @@ -503,7 +504,7 @@ function kkt(z::AbstractVector, net::DCNetwork, d::AbstractVector) K_μ_ub[i] = μ_ub[i] else K_μ_lb[i] = μ_lb[i] * psh[i] - K_μ_ub[i] = μ_ub[i] * (d[i] - psh[i]) + K_μ_ub[i] = μ_ub[i] * (_shed_capacity(d[i]) - psh[i]) end end @@ -762,7 +763,7 @@ function calc_kkt_jacobian(net::DCNetwork, d::AbstractVector, prob::DCOPFProblem # ∂K_psh/∂μ_lb = -I # ∂K_μ_lb/∂μ_lb = Diag(psh), except fixed zero-shed buses omit this term # ∂K_psh/∂μ_ub = I - # ∂K_μ_ub/∂μ_ub = Diag(d - psh), except fixed zero-shed buses use I + # ∂K_μ_ub/∂μ_ub = Diag(max(d, 0) - psh), except fixed zero-shed buses use I @inbounds for i in 1:n start_col!(idx.mu_lb[i]) _push_csc_entry!(rowval, nzval, idx.psh[i], -1.0) @@ -774,7 +775,7 @@ function calc_kkt_jacobian(net::DCNetwork, d::AbstractVector, prob::DCOPFProblem if _is_fixed_zero_shed(d[i]) _push_csc_entry!(rowval, nzval, idx.mu_ub[i], 1.0) else - _push_csc_entry!(rowval, nzval, idx.mu_ub[i], d[i] - psh[i]) + _push_csc_entry!(rowval, nzval, idx.mu_ub[i], _shed_capacity(d[i]) - psh[i]) end end @@ -835,10 +836,11 @@ function calc_kkt_jacobian_demand(net::DCNetwork, d::AbstractVector, sol::DCOPFS sizehint!(rowval, 2n) sizehint!(nzval, 2n) - # ∂K_power_bal/∂d = -I and ∂K_μ_ub/∂d = Diag(μ_ub) for positive-demand buses. + # ∂K_power_bal/∂d = -I and ∂K_μ_ub/∂d = Diag(μ_ub .* ∂max(d, 0)/∂d). @inbounds for i in 1:n colptr[i] = length(rowval) + 1 - _is_fixed_zero_shed(d[i]) || _push_csc_entry!(rowval, nzval, idx.mu_ub[i], mu_ub[i]) + shed_capacity_derivative = _shed_capacity_derivative(d[i]) + _push_csc_entry!(rowval, nzval, idx.mu_ub[i], mu_ub[i] * shed_capacity_derivative) _push_csc_entry!(rowval, nzval, idx.nu_bal[i], -1.0) end colptr[n + 1] = length(rowval) + 1 @@ -856,7 +858,7 @@ function calc_kkt_jacobian_demand_column(net::DCNetwork, d::AbstractVector, sol: col = zeros(kkt_dims(net)) idx = kkt_indices(net) col[idx.nu_bal[j]] = -1.0 - _is_fixed_zero_shed(d[j]) || (col[idx.mu_ub[j]] = sol.mu_ub[j]) + col[idx.mu_ub[j]] = sol.mu_ub[j] * _shed_capacity_derivative(d[j]) return col end diff --git a/src/sens/vjp_jvp.jl b/src/sens/vjp_jvp.jl index c0cc8c4..aa04fa0 100644 --- a/src/sens/vjp_jvp.jl +++ b/src/sens/vjp_jvp.jl @@ -55,7 +55,7 @@ function _dc_param_vjp!(out::AbstractVector, prob::DCOPFProblem, sol::DCOPFSolution, param::Symbol, u::AbstractVector, idx::NamedTuple) if param === :d - _dc_demand_vjp!(out, sol, u, idx, prob.network.n) + _dc_demand_vjp!(out, prob.d, sol, u, idx, prob.network.n) elseif param === :cl _dc_cost_linear_vjp!(out, u, idx, prob.network.k) elseif param === :cq @@ -70,10 +70,11 @@ function _dc_param_vjp!(out::AbstractVector, prob::DCOPFProblem, return out end -# `:d` — 2 nonzeros/col at nu_bal[j]=-1, mu_ub[j]=sol.mu_ub[j] -function _dc_demand_vjp!(out, sol, u, idx, n) +# `:d` — nu_bal[j]=-1, plus mu_ub[j]=sol.mu_ub[j] where max(d[j], 0) is differentiable. +function _dc_demand_vjp!(out, d, sol, u, idx, n) @inbounds for j in 1:n - out[j] = u[idx.nu_bal[j]] - sol.mu_ub[j] * u[idx.mu_ub[j]] + out[j] = u[idx.nu_bal[j]] - + sol.mu_ub[j] * _shed_capacity_derivative(d[j]) * u[idx.mu_ub[j]] end end @@ -122,7 +123,7 @@ function _dc_param_jvp!(v::AbstractVector, prob::DCOPFProblem, sol::DCOPFSolution, param::Symbol, tang::AbstractVector, idx::NamedTuple) if param === :d - _dc_demand_jvp!(v, sol, tang, idx, prob.network.n) + _dc_demand_jvp!(v, prob.d, sol, tang, idx, prob.network.n) elseif param === :cl _dc_cost_linear_jvp!(v, tang, idx, prob.network.k) elseif param === :cq @@ -137,11 +138,11 @@ function _dc_param_jvp!(v::AbstractVector, prob::DCOPFProblem, return v end -# `:d` — scatter -tang[j] into nu_bal, mu_ub[j]*tang[j] into mu_ub -function _dc_demand_jvp!(v, sol, tang, idx, n) +# `:d` — scatter -tang[j] into nu_bal and the positive-part bound derivative into mu_ub. +function _dc_demand_jvp!(v, d, sol, tang, idx, n) @inbounds for j in 1:n v[idx.nu_bal[j]] += -tang[j] - v[idx.mu_ub[j]] += sol.mu_ub[j] * tang[j] + v[idx.mu_ub[j]] += sol.mu_ub[j] * _shed_capacity_derivative(d[j]) * tang[j] end end diff --git a/src/types/dc_opf_problem.jl b/src/types/dc_opf_problem.jl index 52cd1fe..a03602c 100644 --- a/src/types/dc_opf_problem.jl +++ b/src/types/dc_opf_problem.jl @@ -143,6 +143,14 @@ end # DCOPFProblem Constructors # ============================================================================= +""" + _shed_capacity(d::Real) + +Return the curtailable portion of bus demand. Negative net demand represents an +injection, so it remains in power balance but cannot be shed. +""" +@inline _shed_capacity(d::Real) = max(d, zero(d)) + """ DCOPFProblem(network::DCNetwork, d::AbstractVector; optimizer=Ipopt.Optimizer, silent=true) @@ -165,7 +173,7 @@ solve!(prob) function DCOPFProblem(network::DCNetwork, d::AbstractVector; optimizer=Ipopt.Optimizer, silent::Bool=true) length(d) == network.n || throw(DimensionMismatch("Demand vector length $(length(d)) must match number of buses $(network.n)")) - _warn_negative_demand(d) + _debug_negative_demand(d) prob = DCOPFProblem( JuMP.Model(), network, VariableRef[], VariableRef[], VariableRef[], VariableRef[], @@ -232,9 +240,10 @@ function _rebuild_jump_model!(prob::DCOPFProblem) gen_lb = @constraint(model, pg .>= network.gmin) gen_ub = @constraint(model, pg .<= network.gmax) - # Load shedding bounds: 0 ≤ psh ≤ d + # Load shedding bounds: 0 ≤ psh ≤ max(d, 0). Signed net demand remains in + # power balance, but negative net demand is an injection and cannot be shed. shed_lb = @constraint(model, psh .>= 0) - shed_ub = @constraint(model, psh .<= d) + shed_ub = @constraint(model, psh .<= _shed_capacity.(d)) # Reference bus ref_con = @constraint(model, va[network.ref_bus] == 0.0) diff --git a/test/runtests.jl b/test/runtests.jl index 9c11c2b..7141e47 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -870,16 +870,16 @@ include("test_apf_integration.jl") @testset "silence()" begin PowerDiff._SILENCE_WARNINGS[] = false - # Before silence: negative demand warning should fire + # Negative net demand is supported and should not warn. net5 = load_test_case("case5.m") if !isnothing(net5) dc_net = DCNetwork(net5) d = calc_demand_vector(dc_net) d_neg = copy(d) d_neg[1] = -1.0 - @test_warn "Negative demand" DCOPFProblem(dc_net, d_neg) + @test_nowarn DCOPFProblem(dc_net, d_neg) - # After silence: same call should be quiet + # silence() still toggles package warning suppression. PowerDiff.silence() @test PowerDiff._SILENCE_WARNINGS[] == true @test_nowarn DCOPFProblem(dc_net, d_neg) diff --git a/test/test_psh.jl b/test/test_psh.jl index 0111532..ea69069 100644 --- a/test/test_psh.jl +++ b/test/test_psh.jl @@ -420,12 +420,55 @@ end end end - @testset "Negative demand warning" begin - dc_net = _make_2bus_psh(gmax=10.0) - @test_warn "Negative demand" DCOPFProblem(dc_net, [0.0, -0.5]) + @testset "Negative net demand is a non-curtailable injection" begin + dc_net = _make_2bus_psh(gmax=10.0, cq=1.0, tau=0.01) + d = [-0.5, 1.0] + prob = DCOPFProblem(dc_net, d) + sol = solve!(prob) + + # Bus 1 injects 0.5 MW, so the generator only needs to provide the + # remaining 0.5 MW. Negative net demand cannot be shed. + @test sol.pg[1] ≈ 0.5 atol=1e-4 + @test all(abs.(sol.psh) .< 1e-6) - prob = DCOPFProblem(dc_net, [0.0, 1.0]) - @test_warn "Negative demand" update_demand!(prob, [0.0, -0.5]) + B_mat = PowerDiff.calc_susceptance_matrix(dc_net) + residual = dc_net.G_inc * sol.pg + sol.psh - d - B_mat * sol.va + @test norm(residual) < 1e-4 + + # The fixed-zero shedding convention must remain consistent with the + # KKT system and with demand derivatives away from the d=0 kink. + z = flatten_variables(sol, prob) + @test norm(kkt(z, prob, d)) < 1e-4 + + dpg_dd = calc_sensitivity(prob, :pg, :d) + dpsh_dd = calc_sensitivity(prob, :psh, :d) + @test all(isfinite, Matrix(dpg_dd)) + @test all(isfinite, Matrix(dpsh_dd)) + + dlmp_dd = calc_sensitivity(prob, :lmp, :d) + adj = [0.3, -0.2] + tang = [0.1, -0.4] + prob_vjp = DCOPFProblem(dc_net, d) + solve!(prob_vjp) + @test vjp(prob_vjp, :lmp, :d, adj) ≈ dlmp_dd.matrix' * adj atol=1e-8 + prob_jvp = DCOPFProblem(dc_net, d) + solve!(prob_jvp) + @test jvp(prob_jvp, :lmp, :d, tang) ≈ dlmp_dd.matrix * tang atol=1e-8 + + delta = 1e-5 + d_pert = copy(d) + d_pert[1] += delta + sol_pert = solve!(DCOPFProblem(dc_net, d_pert)) + fd = (sol_pert.pg - sol.pg) / delta + @test Matrix(dpg_dd)[:, 1] ≈ fd atol=1e-3 + + # Updating an existing problem must rewrite the shedding cap to max(d, 0). + prob_updated = DCOPFProblem(dc_net, [0.0, 1.0]) + solve!(prob_updated) + update_demand!(prob_updated, d) + sol_updated = solve!(prob_updated) + @test sol_updated.pg ≈ sol.pg atol=1e-5 + @test sol_updated.psh ≈ sol.psh atol=1e-6 end @testset "Degenerate complementarity (d=0 everywhere)" begin From fb629864aa4f90a14691e322a22d67fd953ce2e7 Mon Sep 17 00:00:00 2001 From: Samuel Talkington <10187005+samtalki@users.noreply.github.com> Date: Tue, 9 Jun 2026 00:25:44 -0400 Subject: [PATCH 2/2] Fix wording and clarify demand constraints in dc-opf.md Corrected wording for load shedding cost vector and clarified demand constraints. --- docs/src/math/dc-opf.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/math/dc-opf.md b/docs/src/math/dc-opf.md index 8940718..084e484 100644 --- a/docs/src/math/dc-opf.md +++ b/docs/src/math/dc-opf.md @@ -29,7 +29,7 @@ where: - ``G_{\text{inc}}`` is the ``n \times k`` generator-bus incidence matrix - ``C_q = \operatorname{diag}(c_q)`` contains quadratic cost coefficients - ``c_l`` contains linear cost coefficients -- ``c_{\text{shed}}`` is the load-shedding cost vector +- ``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 @@ -78,7 +78,7 @@ For each parameter type, we compute ``\partial K / \partial p`` and then the ful ### Demand (``d``) -Demand enters the power balance and the upper shedding bound through +Demand enters the power balance and the load shed upper bound constraints through the clipped demand ``d_+ = \max(d, 0)``: ```math @@ -88,8 +88,8 @@ Demand enters the power balance and the upper shedding bound through ``` For strictly positive demand, ``\partial d_+ / \partial d = 1``. For negative -demand, it is ``0``. At zero demand, the positive-part function is non-smooth; -the implementation uses the fixed-zero shedding convention already required by +demand, it is ``0``. At zero demand, the clipping function is non-smooth; +the implementation uses the fixed zero shedding convention already required by the collapsed bound ``0 \leq \text{psh} \leq 0``. ### Switching (``\mathrm{sw}``)