Handle negative net demand in DC load shedding

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}
@@ -29,7 +29,8 @@ 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
 
 ## 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 load shed upper bound constraints through the clipped demand
+``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 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}``)
 
 Switching affects the Laplacian ``B``, weight matrix ``W``, and flow definition through:

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

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
 

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
 

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)

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)

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