From c0791432d0693e454eb77c595de5c8df6cfa0d43 Mon Sep 17 00:00:00 2001 From: ckhanpour3 Date: Sun, 9 Aug 2026 20:55:12 -0400 Subject: [PATCH] Gate DC angle limits by energized branches --- docs/src/math/dc-opf.md | 26 ++++--- src/prob/dc_opf.jl | 7 ++ src/prob/kkt_dc_opf.jl | 60 ++++++++------- src/sens/lmp.jl | 16 ++-- src/sens/susceptance.jl | 10 ++- src/types/dc_network.jl | 21 ++++++ src/types/dc_opf_problem.jl | 8 +- test/test_angle_diff_duals.jl | 137 ++++++++++++++++++++++++++++++++++ 8 files changed, 233 insertions(+), 52 deletions(-) diff --git a/docs/src/math/dc-opf.md b/docs/src/math/dc-opf.md index b305573..05e4a13 100644 --- a/docs/src/math/dc-opf.md +++ b/docs/src/math/dc-opf.md @@ -17,7 +17,7 @@ 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}}) \\ -\mathrm{sw} \circ \alpha_{\min} \leq \mathrm{sw} \circ A\theta &\leq \mathrm{sw} \circ \alpha_{\max} & (\gamma_{\text{lb}}, \gamma_{\text{ub}}) \\ +\mathrm{sw}_{\mathrm{eff}} \circ \alpha_{\min} \leq \mathrm{sw}_{\mathrm{eff}} \circ A\theta &\leq \mathrm{sw}_{\mathrm{eff}} \circ \alpha_{\max} & (\gamma_{\text{lb}}, \gamma_{\text{ub}}) \\ \theta_{\text{refs}} &= 0 & (\eta_{\text{ref}}) \end{aligned} ``` @@ -32,7 +32,7 @@ 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 -- the angle difference limits are gated by ``\mathrm{sw}`` so an open branch (``\mathrm{sw}_e = 0``) imposes no limit; the factor cancels for a branch in service (``\mathrm{sw}_e = 1``) +- ``\mathrm{sw}_{\mathrm{eff},e} = \mathrm{sw}_e`` when ``b_e \mathrm{sw}_e \ne 0`` and is zero otherwise, so angle difference limits use the same energized predicate as the island partition while retaining fractional switching on energized branches - ``\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 @@ -65,12 +65,12 @@ with total dimension ``5n + 6m + 3k + n_{\text{ref}}``. 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 \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 +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}_{\mathrm{eff}}) (\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`` 5. **Complementary slackness (flow bounds)**: ``\lambda_{\text{lb}} \circ (f + f_{\max}) = 0``, ``\lambda_{\text{ub}} \circ (f_{\max} - f) = 0`` -5b. **Complementary slackness (angle differences)**: ``\gamma_{\text{lb}} \circ \mathrm{sw} \circ (A\theta - \alpha_{\min}) = 0``, ``\gamma_{\text{ub}} \circ \mathrm{sw} \circ (\alpha_{\max} - A\theta) = 0`` +5b. **Complementary slackness (angle differences)**: ``\gamma_{\text{lb}} \circ \mathrm{sw}_{\mathrm{eff}} \circ (A\theta - \alpha_{\min}) = 0``, ``\gamma_{\text{ub}} \circ \mathrm{sw}_{\mathrm{eff}} \circ (\alpha_{\max} - A\theta) = 0`` 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`` @@ -102,7 +102,8 @@ the collapsed bound ``0 \leq \text{psh} \leq 0``. ### Switching (``\mathrm{sw}``) -Switching enters the Laplacian ``B``, the weight matrix ``W``, and the angle difference limits gated by ``\mathrm{sw}``, with elementary perturbations: +Switching enters the Laplacian ``B``, the weight matrix ``W``, and the angle difference limits gated by ``\mathrm{sw}_{\mathrm{eff}}``. Within a fixed energized regime, +``\partial \mathrm{sw}_{\mathrm{eff},e}/\partial \mathrm{sw}_e = 1`` when ``b_e \ne 0`` and is zero when ``b_e = 0``. The elementary network perturbations are: ```math \frac{\partial B}{\partial \mathrm{sw}_e} = -b_e \, a_e a_e^\top, \qquad @@ -121,17 +122,20 @@ where ``a_e^\top`` is row ``e`` of ``A`` (so ``a_e^\top \theta = (A\theta)_e``) && \text{(flow definition)} \\ \frac{\partial K_{\theta}}{\partial \mathrm{sw}_e} &= -b_e\,\bigl(a_e^\top \nu_{\text{bal}} + (\nu_{\text{flow}})_e\bigr)\, a_e - + (\gamma_{\text{ub},e} - \gamma_{\text{lb},e})\, a_e + + \frac{\partial \mathrm{sw}_{\mathrm{eff},e}}{\partial \mathrm{sw}_e} + (\gamma_{\text{ub},e} - \gamma_{\text{lb},e})\, a_e && (\theta\text{ stationarity}) \\ \frac{\partial K_{\gamma_{\text{lb}}}}{\partial \mathrm{sw}_e} - &= \gamma_{\text{lb},e}\,\bigl((A\theta)_e - \alpha_{\min,e}\bigr)\, e_e, \quad + &= \frac{\partial \mathrm{sw}_{\mathrm{eff},e}}{\partial \mathrm{sw}_e} + \gamma_{\text{lb},e}\,\bigl((A\theta)_e - \alpha_{\min,e}\bigr)\, e_e, \quad \frac{\partial K_{\gamma_{\text{ub}}}}{\partial \mathrm{sw}_e} - = \gamma_{\text{ub},e}\,\bigl(\alpha_{\max,e} - (A\theta)_e\bigr)\, e_e + = \frac{\partial \mathrm{sw}_{\mathrm{eff},e}}{\partial \mathrm{sw}_e} + \gamma_{\text{ub},e}\,\bigl(\alpha_{\max,e} - (A\theta)_e\bigr)\, e_e && \text{(angle limits)} \end{aligned} ``` -Each block is rank-1 in the incidence row ``a_e`` or supported on branch ``e`` alone, so every column of ``\partial K / \partial \mathrm{sw}`` has only a handful of nonzeros. The stationarity block combines the Laplacian term (``\partial B / \partial \mathrm{sw}_e``), the flow coupling term (``\partial W / \partial \mathrm{sw}_e``), and the gated angle limit term (``\partial \operatorname{diag}(\mathrm{sw}) / \partial \mathrm{sw}_e = e_e e_e^\top``). +Each block is rank-1 in the incidence row ``a_e`` or supported on branch ``e`` alone, so every column of ``\partial K / \partial \mathrm{sw}`` has only a handful of nonzeros. The stationarity block combines the Laplacian term (``\partial B / \partial \mathrm{sw}_e``), the flow coupling term (``\partial W / \partial \mathrm{sw}_e``), and the effective angle gate derivative above. 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 @@ -159,7 +163,7 @@ Flow limits enter the complementary slackness conditions: ### Susceptances (``b``) -Susceptances affect the same blocks as switching (through ``B`` and ``W``), but with different partial derivatives since ``B = A^\top \operatorname{diag}(-b \circ \mathrm{sw}) A``. +Susceptances affect the same blocks as switching through ``B`` and ``W``, but with different partial derivatives since ``B = A^\top \operatorname{diag}(-b \circ \mathrm{sw}) A``. The effective angle gate is locally constant in ``b`` on either side of zero, so its ``b`` derivative is zero within a fixed energized regime. At ``b_e = 0`` the gate and, potentially, the island partition change nonsmoothly; the analytical parameter Jacobian uses a zero gate derivative at that boundary. ## LMP Decomposition @@ -172,7 +176,7 @@ Locational marginal prices are the power balance duals ``\nu_{\text{bal}}``, dec 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 \operatorname{diag}(\mathrm{sw}) (\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}_{\mathrm{eff}}) (\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 within each energized island and reflects that island's marginal cost of generation. diff --git a/src/prob/dc_opf.jl b/src/prob/dc_opf.jl index 4abd08a..5c7a57f 100644 --- a/src/prob/dc_opf.jl +++ b/src/prob/dc_opf.jl @@ -105,6 +105,13 @@ function solve!(prob::DCOPFProblem) Atheta = net.A * θ_val TOL = COMPLEMENTARITY_SNAP_TOL for e in 1:net.m + if iszero(_angle_difference_gate(net, e)) + # A de-energized branch has a constant 0 <= 0 pair of constraints. + # Canonicalize its otherwise arbitrary duals to zero. + γ_lb[e] = 0.0 + γ_ub[e] = 0.0 + continue + end if net.angmax[e] - Atheta[e] > TOL # upper angle not binding γ_ub[e] = 0.0 else diff --git a/src/prob/kkt_dc_opf.jl b/src/prob/kkt_dc_opf.jl index 68a1ec0..79bfba2 100644 --- a/src/prob/kkt_dc_opf.jl +++ b/src/prob/kkt_dc_opf.jl @@ -445,8 +445,8 @@ s.t. G_inc * g + psh - d = B * θ (ν_bal) f = W * A * θ (ν_flow) f ≥ -fmax (λ_lb) f ≤ fmax (λ_ub) - A * θ ≥ angmin (γ_lb) - A * θ ≤ angmax (γ_ub) + sw_eff .* (A * θ) ≥ sw_eff .* angmin (γ_lb) + sw_eff .* (A * θ) ≤ sw_eff .* angmax (γ_ub) g ≥ gmin (ρ_lb) g ≤ gmax (ρ_ub) 0 ≤ psh (μ_lb) @@ -454,9 +454,11 @@ s.t. G_inc * g + psh - d = B * θ (ν_bal) θ[refs] = 0 (η_ref) ``` +Here `sw_eff[e] = sw[e]` when `b[e] * sw[e] != 0` and zero otherwise. + # Returns Vector of KKT residuals (should be zero at optimum): -1. Stationarity w.r.t. θ: B' * ν_bal + (W*A)' * ν_flow + E_ref * η_ref + A' * Diag(sw) * (γ_ub - γ_lb) = 0 +1. Stationarity w.r.t. θ: B' * ν_bal + (W*A)' * ν_flow + E_ref * η_ref + A' * Diag(sw_eff) * (γ_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 @@ -490,12 +492,13 @@ function kkt(z::AbstractVector, net::DCNetwork, d::AbstractVector) W = Diagonal(-net.b .* net.sw) B_mat = net.A' * W * net.A WA = W * net.A + angle_gate = _angle_difference_gates(net) refs = _reference_buses(net) # KKT conditions # 1. Stationarity w.r.t. θ - K_θ = B_mat' * ν_bal + WA' * ν_flow + net.A' * (net.sw .* (γ_ub - γ_lb)) + K_θ = B_mat' * ν_bal + WA' * ν_flow + net.A' * (angle_gate .* (γ_ub - γ_lb)) K_θ[refs] .+= η_ref # 2. Stationarity w.r.t. g @@ -513,8 +516,8 @@ function kkt(z::AbstractVector, net::DCNetwork, d::AbstractVector) # 6. Complementary slackness: phase angle difference bounds Aθ = net.A * θ - K_γ_lb = γ_lb .* net.sw .* (Aθ - net.angmin) - K_γ_ub = γ_ub .* net.sw .* (net.angmax - Aθ) + K_γ_lb = γ_lb .* angle_gate .* (Aθ - net.angmin) + K_γ_ub = γ_ub .* angle_gate .* (net.angmax - Aθ) # 7. Complementary slackness: generation bounds K_ρ_lb = ρ_lb .* (g - net.gmin) @@ -620,6 +623,7 @@ function calc_kkt_jacobian(net::DCNetwork, d::AbstractVector, prob::DCOPFProblem G_inc = getfield(net, :G_inc) b = getfield(net, :b) sw = getfield(net, :sw) + angle_gate = _angle_difference_gates(net) fmax = getfield(net, :fmax) gmin = getfield(net, :gmin) gmax = getfield(net, :gmax) @@ -668,8 +672,8 @@ function calc_kkt_jacobian(net::DCNetwork, d::AbstractVector, prob::DCOPFProblem # SparseMatrixCSC diagonals and structural setindex! insertions. # va columns: - # ∂K_γ_lb/∂θ = Diag(γ_lb .* sw) * A - # ∂K_γ_ub/∂θ = -Diag(γ_ub .* sw) * A + # ∂K_γ_lb/∂θ = Diag(γ_lb .* sw_eff) * A + # ∂K_γ_ub/∂θ = -Diag(γ_ub .* sw_eff) * A # ∂K_power_bal/∂θ = -B # ∂K_flow_def/∂θ = -W*A, where W = Diag(-b .* sw) # ∂K_ref/∂θ_ref = 1 @@ -678,12 +682,12 @@ function calc_kkt_jacobian(net::DCNetwork, d::AbstractVector, prob::DCOPFProblem for p in nzrange(A, j) e = rowvals(A)[p] aej = nonzeros(A)[p] - _push_csc_entry!(rowval, nzval, idx.gamma_lb[e], gamma_lb[e] * sw[e] * aej) + _push_csc_entry!(rowval, nzval, idx.gamma_lb[e], gamma_lb[e] * angle_gate[e] * aej) end for p in nzrange(A, j) e = rowvals(A)[p] aej = nonzeros(A)[p] - _push_csc_entry!(rowval, nzval, idx.gamma_ub[e], -gamma_ub[e] * sw[e] * aej) + _push_csc_entry!(rowval, nzval, idx.gamma_ub[e], -gamma_ub[e] * angle_gate[e] * aej) end for p in nzrange(B_mat, j) _push_csc_entry!(rowval, nzval, idx.nu_bal[rowvals(B_mat)[p]], -nonzeros(B_mat)[p]) @@ -755,23 +759,23 @@ function calc_kkt_jacobian(net::DCNetwork, d::AbstractVector, prob::DCOPFProblem end # gamma columns: - # ∂K_θ/∂γ_lb = -A' * Diag(sw) - # ∂K_γ_lb/∂γ_lb = Diag(sw .* (A*θ - angmin)) - # ∂K_θ/∂γ_ub = A' * Diag(sw) - # ∂K_γ_ub/∂γ_ub = Diag(sw .* (angmax - A*θ)) + # ∂K_θ/∂γ_lb = -A' * Diag(sw_eff) + # ∂K_γ_lb/∂γ_lb = Diag(sw_eff .* (A*θ - angmin)) + # ∂K_θ/∂γ_ub = A' * Diag(sw_eff) + # ∂K_γ_ub/∂γ_ub = Diag(sw_eff .* (angmax - A*θ)) @inbounds for e in 1:m start_col!(idx.gamma_lb[e]) for p in A_rowptr[e]:(A_rowptr[e + 1] - 1) - _push_csc_entry!(rowval, nzval, idx.va[A_rowcols[p]], -sw[e] * A_rowvals[p]) + _push_csc_entry!(rowval, nzval, idx.va[A_rowcols[p]], -angle_gate[e] * A_rowvals[p]) end - _push_csc_entry!(rowval, nzval, idx.gamma_lb[e], sw[e] * (Aθ[e] - angmin[e])) + _push_csc_entry!(rowval, nzval, idx.gamma_lb[e], angle_gate[e] * (Aθ[e] - angmin[e])) end @inbounds for e in 1:m start_col!(idx.gamma_ub[e]) for p in A_rowptr[e]:(A_rowptr[e + 1] - 1) - _push_csc_entry!(rowval, nzval, idx.va[A_rowcols[p]], sw[e] * A_rowvals[p]) + _push_csc_entry!(rowval, nzval, idx.va[A_rowcols[p]], angle_gate[e] * A_rowvals[p]) end - _push_csc_entry!(rowval, nzval, idx.gamma_ub[e], sw[e] * (angmax[e] - Aθ[e])) + _push_csc_entry!(rowval, nzval, idx.gamma_ub[e], angle_gate[e] * (angmax[e] - Aθ[e])) end # rho columns: @@ -926,8 +930,6 @@ function calc_kkt_jacobian_switching(prob::DCOPFProblem, sol::DCOPFSolution) θ = sol.va - # Current switching state - s = net.sw b = net.b A = net.A @@ -953,12 +955,15 @@ 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 - J_s[idx.va, e] += A_e_vec * (sol.gamma_ub[e] - sol.gamma_lb[e]) + # ∂K_θ/∂s_e from gated angle difference bounds. When b[e] == 0, + # changing switching alone cannot energize the branch, so this derivative + # vanishes with the gate. + dgate_dsw = _angle_difference_gate_dsw(net, e) + J_s[idx.va, e] += dgate_dsw * A_e_vec * (sol.gamma_ub[e] - sol.gamma_lb[e]) # ∂K_γ/∂s_e from gated complementary slackness - J_s[idx.gamma_lb[e], e] = sol.gamma_lb[e] * (Aθ_e - net.angmin[e]) - J_s[idx.gamma_ub[e], e] = sol.gamma_ub[e] * (net.angmax[e] - Aθ_e) + J_s[idx.gamma_lb[e], e] = dgate_dsw * sol.gamma_lb[e] * (Aθ_e - net.angmin[e]) + J_s[idx.gamma_ub[e], e] = dgate_dsw * sol.gamma_ub[e] * (net.angmax[e] - Aθ_e) end return J_s @@ -1002,12 +1007,13 @@ function calc_kkt_jacobian_switching_column(prob::DCOPFProblem, sol::DCOPFSoluti col[idx.va[f_bus]] += coeff col[idx.va[t_bus]] -= coeff # gated angle difference stationarity contribution - coeff_ang = sol.gamma_ub[e] - sol.gamma_lb[e] + dgate_dsw = _angle_difference_gate_dsw(net, e) + coeff_ang = dgate_dsw * (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 - 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) + col[idx.gamma_lb[e]] = dgate_dsw * sol.gamma_lb[e] * (Aθ_e - net.angmin[e]) + col[idx.gamma_ub[e]] = dgate_dsw * sol.gamma_ub[e] * (net.angmax[e] - Aθ_e) return col end diff --git a/src/sens/lmp.jl b/src/sens/lmp.jl index 2d447d5..c3dfce8 100644 --- a/src/sens/lmp.jl +++ b/src/sens/lmp.jl @@ -25,7 +25,7 @@ # LMP Decomposition (for analysis): # LMP = ν_bal = energy_component + congestion_component # where: -# congestion_component = B_r⁻¹ [A_r' Diag(-b .* sw) (λ_ub - λ_lb) + A_r' Diag(sw) (γ_ub - γ_lb)] (non-ref block) +# congestion_component = B_r⁻¹ [A_r' Diag(-b .* sw) (λ_ub - λ_lb) + A_r' Diag(sw_eff) (γ_ub - γ_lb)] (non-ref block) # energy_component = ν_bal - congestion_component (uniform within each island) # # Sign conventions (DC OPF): @@ -79,11 +79,11 @@ Extract the congestion component of LMPs for analysis. 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 + B' * ν_bal + (WA)' * ν_flow + E_ref * η_ref + A' Diag(sw_eff) (γ_ub - γ_lb) = 0 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] + congestion[non_ref] = B_r \\ (A' W (λ_ub - λ_lb) + A' Diag(sw_eff) (γ_ub - γ_lb))[non_ref] The congestion component captures price differentiation due to binding flow and angle constraints. Only the non-reference rows are populated, so every energized island @@ -96,15 +96,13 @@ 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 = _non_reference_buses(net) + angle_gate = _angle_difference_gates(net) At = net.A' - # 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. + # Match the model's effective angle gate: preserve fractional switching on + # energized branches and remove the contribution whenever b[e] * sw[e] == 0. rhs_full = At * Diagonal(w .* net.sw) * (sol.lam_ub - sol.lam_lb) + - At * Diagonal(net.sw) * (sol.gamma_ub - sol.gamma_lb) + At * Diagonal(angle_gate) * (sol.gamma_ub - sol.gamma_lb) result = zeros(net.n) result[non_ref] = B_r_factor \ rhs_full[non_ref] diff --git a/src/sens/susceptance.jl b/src/sens/susceptance.jl index 1fe599e..ee04c36 100644 --- a/src/sens/susceptance.jl +++ b/src/sens/susceptance.jl @@ -39,8 +39,8 @@ Susceptance b affects: 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_theta = B' * nu_bal + (WA)' * nu_flow + E_ref * eta_ref + A' * Diag(sw_eff) * (gamma_ub - gamma_lb) + (within a fixed energized regime, sw_eff is locally constant in b, 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 @@ -66,6 +66,9 @@ function calc_kkt_jacobian_susceptance(prob::DCOPFProblem, sol::DCOPFSolution) # Precompute A * theta once (invariant across branches) Atheta = A * theta + # The effective angle gate is locally constant in b on either side of zero. + # Its derivative is also defined as zero at the nonsmooth b == 0 boundary, + # so the gamma blocks make no contribution to this Jacobian. for e in 1:m A_e_vec = Vector(A[e, :]) Atheta_e = Atheta[e] @@ -79,6 +82,7 @@ function calc_kkt_jacobian_susceptance(prob::DCOPFProblem, sol::DCOPFSolution) # dK_flow_def/db_e: only row e is nonzero J_b[idx.nu_flow[e], e] = sw[e] * Atheta_e + end return J_b @@ -106,5 +110,7 @@ function calc_kkt_jacobian_susceptance_column(prob::DCOPFProblem, sol::DCOPFSolu coeff = -sw[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 + # The effective angle gate has zero b derivative within a fixed energized + # regime and by convention at the nonsmooth b[e] == 0 boundary. return col end diff --git a/src/types/dc_network.jl b/src/types/dc_network.jl index 090b93d..9e9dc63 100644 --- a/src/types/dc_network.jl +++ b/src/types/dc_network.jl @@ -662,6 +662,27 @@ end @inline _is_energized(net::DCNetwork, e::Int) = !iszero(getfield(net, :b)[e] * getfield(net, :sw)[e]) +# Phase-angle bounds use the switching value on energized branches and vanish +# everywhere else. This preserves the existing continuous switching scaling while +# making the bound topology agree with the energized-island predicate above. +@inline _angle_difference_gate(net::DCNetwork, e::Int) = + _is_energized(net, e) ? getfield(net, :sw)[e] : zero(getfield(net, :sw)[e]) + +function _angle_difference_gates(net::DCNetwork) + gates = similar(getfield(net, :sw)) + @inbounds for e in eachindex(gates) + gates[e] = _angle_difference_gate(net, e) + end + return gates +end + +# Within a fixed energized regime, the gate is `sw` when b != 0 and zero when +# b == 0. Crossing b == 0 is a nonsmooth topology boundary, so its local +# susceptance derivative follows the same zero convention used at other topology +# boundaries. +@inline _angle_difference_gate_dsw(net::DCNetwork, e::Int) = + iszero(getfield(net, :b)[e]) ? zero(getfield(net, :sw)[e]) : one(getfield(net, :sw)[e]) + function _topology_cache_valid(net::DCNetwork) cache = getfield(net, :topology_cache) cache.initialized || return false diff --git a/src/types/dc_opf_problem.jl b/src/types/dc_opf_problem.jl index f6f4637..3894014 100644 --- a/src/types/dc_opf_problem.jl +++ b/src/types/dc_opf_problem.jl @@ -269,9 +269,11 @@ function _rebuild_jump_model!(prob::DCOPFProblem) 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) - phase_diff_ub = @constraint(model, network.sw .* (network.A * va) .<= network.sw .* network.angmax) + # De-energized lines should not constrain angle differences. Use the same + # b[e] * sw[e] != 0 predicate as the energized-island partition. + angle_gate = _angle_difference_gates(network) + phase_diff_lb = @constraint(model, angle_gate .* (network.A * va) .>= angle_gate .* network.angmin) + phase_diff_ub = @constraint(model, angle_gate .* (network.A * va) .<= angle_gate .* network.angmax) prob.model = model prob.va = va diff --git a/test/test_angle_diff_duals.jl b/test/test_angle_diff_duals.jl index 28d185b..b553e9c 100644 --- a/test/test_angle_diff_duals.jl +++ b/test/test_angle_diff_duals.jl @@ -421,4 +421,141 @@ import PowerDiff: kkt, kkt_indices, flatten_variables ungated[non_ref] = sol.B_r_factor \ rhs_ungated[non_ref] @test !isapprox(cong, ungated; atol=1e-3) end + + @testset "Zero-susceptance angle gating" begin + A_gate = sparse([ + 1.0 -1.0 0.0 + 0.0 1.0 -1.0 + 1.0 0.0 -1.0 + ]) + G_gate = sparse(Matrix(1.0I, 3, 3)) + d_gate = [0.1, 0.1, 0.1] + + function make_gate_network(b3; exclude_raw_angle=false) + angle_min = exclude_raw_angle ? [-Float64(pi), -Float64(pi), 0.1] : fill(-Float64(pi), 3) + angle_max = exclude_raw_angle ? [Float64(pi), Float64(pi), 0.2] : fill(Float64(pi), 3) + return DCNetwork(3, 3, 3, A_gate, G_gate, [-10.0, -10.0, b3]; + sw=[1.0, 0.8, 0.6], fmax=fill(10.0, 3), + gmax=fill(10.0, 3), gmin=zeros(3), + angmin=angle_min, angmax=angle_max, + cq=ones(3), cl=[10.0, 20.0, 30.0], ref_bus=1, tau=1e-3) + end + + function with_test_gamma(sol, e) + gamma_lb = copy(sol.gamma_lb) + gamma_ub = copy(sol.gamma_ub) + gamma_lb[e] = 1.25 + gamma_ub[e] = 3.5 + return DCOPFSolution( + copy(sol.va), copy(sol.pg), copy(sol.f), copy(sol.psh), + copy(sol.nu_bal), copy(sol.nu_flow), + copy(sol.lam_ub), copy(sol.lam_lb), + copy(sol.rho_ub), copy(sol.rho_lb), + copy(sol.mu_lb), copy(sol.mu_ub), + gamma_lb, gamma_ub, copy(sol.eta_ref), sol.objective, sol.B_r_factor, + ) + end + + function finite_difference_kkt_column(net, z, field, e, h) + values = getfield(net, field) + base = values[e] + try + values[e] = base + h + plus = kkt(z, net, d_gate) + values[e] = base - h + minus = kkt(z, net, d_gate) + return (plus - minus) / (2h) + finally + values[e] = base + reference_buses(net) # refresh the topology cache after restoration + end + end + + # A zero-b bridge splits the buses into separate energized islands. Its + # raw angle interval excludes zero, so the old sw-only gate made this + # otherwise independent two-island problem infeasible. + A_bridge = sparse(reshape([1.0, -1.0], 1, 2)) + net_bridge = DCNetwork(2, 1, 2, A_bridge, sparse(Matrix(1.0I, 2, 2)), [0.0]; + sw=[1.0], fmax=[10.0], gmax=[10.0, 10.0], gmin=zeros(2), + angmin=[0.1], angmax=[0.2], cq=ones(2), cl=[10.0, 20.0], + ref_bus=1, tau=1e-3) + @test reference_buses(net_bridge) == [1, 2] + prob_bridge = DCOPFProblem(net_bridge, [0.1, 0.1]) + @test length(prob_bridge.cons.ref) == 2 + sol_bridge = solve!(prob_bridge) + @test (net_bridge.A * sol_bridge.va)[1] < net_bridge.angmin[1] + @test iszero(sol_bridge.gamma_lb[1]) + @test iszero(sol_bridge.gamma_ub[1]) + + # The zero-b branch is redundant for connectivity, but its raw angle + # interval deliberately excludes the solved angle. It must impose no model + # constraint because b[3] * sw[3] == 0. + net_zero = make_gate_network(0.0; exclude_raw_angle=true) + @test PowerDiff._angle_difference_gates(net_zero) == [1.0, 0.8, 0.0] + @test reference_buses(net_zero) == [1] + prob_zero = DCOPFProblem(net_zero, d_gate) + sol_zero = solve!(prob_zero) + @test (net_zero.A * sol_zero.va)[3] < net_zero.angmin[3] + @test iszero(sol_zero.gamma_lb[3]) + @test iszero(sol_zero.gamma_ub[3]) + + z_zero = flatten_variables(sol_zero, prob_zero) + idx_zero = kkt_indices(prob_zero) + @test norm(kkt(z_zero, prob_zero, d_gate), Inf) < 1e-5 + + # Nonzero test multipliers make every gamma path load bearing without + # relying on an arbitrary solver dual for the constant 0 <= 0 constraints. + sol_zero_gamma = with_test_gamma(sol_zero, 3) + z_zero_gamma = flatten_variables(sol_zero_gamma, prob_zero) + @test kkt(z_zero_gamma, prob_zero, d_gate) ≈ kkt(z_zero, prob_zero, d_gate) atol=1e-12 + + Jz_zero = Matrix(PowerDiff.calc_kkt_jacobian(prob_zero; sol=sol_zero_gamma)) + Jz_zero_fd = ForwardDiff.jacobian(z -> kkt(z, prob_zero, d_gate), z_zero_gamma) + @test Jz_zero ≈ Jz_zero_fd atol=1e-10 rtol=1e-10 + @test all(iszero, Jz_zero[:, idx_zero.gamma_lb[3]]) + @test all(iszero, Jz_zero[:, idx_zero.gamma_ub[3]]) + + Jsw_zero = Matrix(PowerDiff.calc_kkt_jacobian_switching(prob_zero, sol_zero_gamma)) + Jsw_zero_col = PowerDiff.calc_kkt_jacobian_switching_column(prob_zero, sol_zero_gamma, 3) + Jsw_zero_fd = finite_difference_kkt_column(net_zero, z_zero_gamma, :sw, 3, 1e-6) + @test all(iszero, Jsw_zero[:, 3]) + @test Jsw_zero_col == Jsw_zero[:, 3] + @test Jsw_zero[:, 3] ≈ Jsw_zero_fd atol=1e-10 + + Jb_zero = Matrix(PowerDiff.calc_kkt_jacobian_susceptance(prob_zero, sol_zero_gamma)) + Jb_zero_plain = Matrix(PowerDiff.calc_kkt_jacobian_susceptance(prob_zero, sol_zero)) + @test Jb_zero ≈ Jb_zero_plain atol=1e-12 + @test PowerDiff.calc_kkt_jacobian_susceptance_column(prob_zero, sol_zero_gamma, 3) ≈ Jb_zero[:, 3] atol=1e-12 + + # A synthetic gamma on the zero-b branch must not leak into congestion. + @test calc_congestion_component(sol_zero_gamma, net_zero) ≈ + calc_congestion_component(sol_zero, net_zero) atol=1e-12 + + # On either side of zero, sw_eff == sw and is locally constant in b. + # Verify both analytic parameter Jacobians without crossing the nonsmooth + # b == 0 boundary. + for b3 in (-1e-3, 1e-3) + net_side = make_gate_network(b3) + prob_side = DCOPFProblem(net_side, d_gate) + sol_side = solve!(prob_side) + sol_side_gamma = with_test_gamma(sol_side, 3) + z_side_gamma = flatten_variables(sol_side_gamma, prob_side) + + Jsw = Matrix(PowerDiff.calc_kkt_jacobian_switching(prob_side, sol_side_gamma)) + Jsw_plain = Matrix(PowerDiff.calc_kkt_jacobian_switching(prob_side, sol_side)) + Jsw_col = PowerDiff.calc_kkt_jacobian_switching_column(prob_side, sol_side_gamma, 3) + Jsw_fd = finite_difference_kkt_column(net_side, z_side_gamma, :sw, 3, 1e-6) + @test Jsw_col ≈ Jsw[:, 3] atol=1e-12 + @test Jsw[:, 3] ≈ Jsw_fd atol=1e-8 rtol=1e-7 + @test norm(Jsw[:, 3] - Jsw_plain[:, 3]) > 1.0 + + Jb = Matrix(PowerDiff.calc_kkt_jacobian_susceptance(prob_side, sol_side_gamma)) + Jb_plain = Matrix(PowerDiff.calc_kkt_jacobian_susceptance(prob_side, sol_side)) + Jb_col = PowerDiff.calc_kkt_jacobian_susceptance_column(prob_side, sol_side_gamma, 3) + Jb_fd = finite_difference_kkt_column(net_side, z_side_gamma, :b, 3, 1e-5) + @test Jb_col ≈ Jb[:, 3] atol=1e-12 + @test Jb[:, 3] ≈ Jb_fd atol=1e-8 rtol=1e-7 + @test Jb ≈ Jb_plain atol=1e-12 + end + end end