Skip to content

Gate DC angle limits by energized branches #71

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
cameronkhanpour wants to merge 1 commit into from
Open

Gate DC angle limits by energized branches #71

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 15 additions & 11 deletions docs/src/math/dc-opf.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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}
```
Expand All @@ -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
Expand Down Expand Up @@ -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``
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -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.
7 changes: 7 additions & 0 deletions src/prob/dc_opf.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
60 changes: 33 additions & 27 deletions src/prob/kkt_dc_opf.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -445,18 +445,20 @@ 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)
psh ≤ max(d, 0) (μ_ub)
θ[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
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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])
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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

Expand Down
16 changes: 7 additions & 9 deletions src/sens/lmp.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down Expand Up @@ -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
Expand All @@ -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]
Expand Down
Loading
Loading