From b0ed0220421fb1349e19d5b0746b1de8b3c07efe Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 14 Apr 2026 18:52:19 +0000 Subject: [PATCH 1/4] Initial plan From df968ebf99955d47d38d09f516eeb43b7f17097c Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 14 Apr 2026 19:23:04 +0000 Subject: [PATCH 2/4] Support ThreeWindingTransformer outages in NetworkModification - Replace error in _classify_outage_component! for PSY.ThreeWindingTransformer with automatic decomposition into available winding arcs - Add NetworkModification(mat, branch::PSY.ThreeWindingTransformer) constructor - Add _classify_branch_modification for PSY.ThreeWindingTransformer - Handle 3WT Ybus delta directly via ybus_branch_entries to avoid BA-matrix vs winding susceptance mismatch - Use get_name instead of PSY.get_name for ThreeWindingTransformerWinding compat Agent-Logs-Url: https://github.com/NREL-Sienna/PowerNetworkMatrices.jl/sessions/7950c45a-eba1-4cd1-a6d6-3e4d6238e999 Co-authored-by: jd-lara <16385323+jd-lara@users.noreply.github.com> --- src/network_modification.jl | 114 +++++++++++++++++++++++++++++------- 1 file changed, 94 insertions(+), 20 deletions(-) diff --git a/src/network_modification.jl b/src/network_modification.jl index a00ea3d6c..e75b350e1 100644 --- a/src/network_modification.jl +++ b/src/network_modification.jl @@ -113,6 +113,28 @@ Construct a `NetworkModification` from a branch component using network reduction reverse maps to classify the branch as direct, parallel, or series. """ function NetworkModification(mat::PowerNetworkMatrix, branch::PSY.ACTransmission) + nr = get_network_reduction_data(mat) + arc_lookup = get_arc_lookup(mat) + arc_sus = _get_arc_susceptances(mat) + mods = _classify_branch_modification(nr, arc_lookup, arc_sus, branch) + return NetworkModification( + get_name(branch), + mods, + ) +end + +""" +$(TYPEDSIGNATURES) + +Construct a `NetworkModification` from a `ThreeWindingTransformer` component. +Automatically decomposes the transformer into its three winding arcs and classifies +each one. For a partial outage (single winding trip), use a +`ThreeWindingTransformerWinding` instead. +""" +function NetworkModification( + mat::PowerNetworkMatrix, + branch::PSY.ThreeWindingTransformer, +) nr = get_network_reduction_data(mat) arc_lookup = get_arc_lookup(mat) arc_sus = _get_arc_susceptances(mat) @@ -238,11 +260,21 @@ function _classify_outage_component!( ) tag, arc_tuple = _resolve_branch_arc(nr, component) - if tag === :direct || tag === :transformer3w + if tag === :direct arc_idx = arc_lookup[arc_tuple] b_arc = arc_susceptances[arc_idx] dy11, dy12, dy21, dy22 = _compute_arc_ybus_delta(nr, arc_tuple, -b_arc) push!(direct_mods, ArcModification(arc_idx, -b_arc, dy11, dy12, dy21, dy22)) + elseif tag === :transformer3w + arc_idx = arc_lookup[arc_tuple] + b_arc = arc_susceptances[arc_idx] + tr = nr.transformer3W_map[arc_tuple] + Y11, Y12, Y21, Y22 = ybus_branch_entries(tr) + push!(direct_mods, ArcModification( + arc_idx, -b_arc, + YBUS_ELTYPE(-Y11), YBUS_ELTYPE(-Y12), + YBUS_ELTYPE(-Y21), YBUS_ELTYPE(-Y22), + )) elseif tag === :parallel arc_idx = arc_lookup[arc_tuple] b_circuit = PSY.get_series_susceptance(component) @@ -256,11 +288,11 @@ function _classify_outage_component!( end push!(series_components_by_arc[arc_idx], component) else - @info "Branch $(PSY.get_name(component)) not found in any reduction map. " * + @info "Branch $(get_name(component)) not found in any reduction map. " * "The component may have been eliminated by a radial reduction." return end - push!(component_names, PSY.get_name(component)) + push!(component_names, get_name(component)) return end @@ -325,23 +357,31 @@ function _classify_outage_component!( end function _classify_outage_component!( - ::NetworkReductionData, - ::Dict, - ::Vector{Float64}, - ::Dict{Int, Int}, + nr::NetworkReductionData, + arc_lookup::Dict, + arc_susceptances::Vector{Float64}, + bus_lookup::Dict{Int, Int}, component::PSY.ThreeWindingTransformer, - ::Vector{ArcModification}, - ::Vector{ArcModification}, - ::Dict{Int, Vector{PSY.ACTransmission}}, - ::Dict{Int, Tuple{Int, Int}}, - ::Vector{ShuntModification}, - ::Vector{String}, + direct_mods::Vector{ArcModification}, + parallel_mods::Vector{ArcModification}, + series_components_by_arc::Dict{Int, Vector{PSY.ACTransmission}}, + series_arc_tuples::Dict{Int, Tuple{Int, Int}}, + shunt_mods::Vector{ShuntModification}, + component_names::Vector{String}, ) - error( - "Outages on ThreeWindingTransformer components are not yet supported. " * - "Component: $(PSY.get_name(component)). " * - "Use individual ThreeWindingTransformerWinding arcs instead.", - ) + for winding_num in 1:3 + winding = ThreeWindingTransformerWinding(component, winding_num) + if !get_equivalent_available(winding) + continue + end + _classify_outage_component!( + nr, arc_lookup, arc_susceptances, bus_lookup, winding, + direct_mods, parallel_mods, + series_components_by_arc, series_arc_tuples, + shunt_mods, component_names, + ) + end + return end """ @@ -360,6 +400,30 @@ function _classify_branch_modification( _assert_not_phase_shifting(branch) end +""" + _classify_branch_modification(nr, arc_lookup, arc_susceptances, branch::PSY.ThreeWindingTransformer) -> Vector{ArcModification} + +Classify a `ThreeWindingTransformer` by decomposing it into its three winding arcs +and classifying each one individually. Returns arc modifications for all windings +present in the network. +""" +function _classify_branch_modification( + nr::NetworkReductionData, + arc_lookup::Dict, + arc_susceptances::Vector{Float64}, + branch::PSY.ThreeWindingTransformer, +)::Vector{ArcModification} + mods = ArcModification[] + for winding_num in 1:3 + winding = ThreeWindingTransformerWinding(branch, winding_num) + if !get_equivalent_available(winding) + continue + end + append!(mods, _classify_branch_modification(nr, arc_lookup, arc_susceptances, winding)) + end + return mods +end + function _classify_branch_modification( nr::NetworkReductionData, arc_lookup::Dict, @@ -368,11 +432,21 @@ function _classify_branch_modification( )::Vector{ArcModification} tag, arc_tuple = _resolve_branch_arc(nr, branch) - if tag === :direct || tag === :transformer3w + if tag === :direct arc_idx = arc_lookup[arc_tuple] b_arc = arc_susceptances[arc_idx] dy11, dy12, dy21, dy22 = _compute_arc_ybus_delta(nr, arc_tuple, -b_arc) return [ArcModification(arc_idx, -b_arc, dy11, dy12, dy21, dy22)] + elseif tag === :transformer3w + arc_idx = arc_lookup[arc_tuple] + b_arc = arc_susceptances[arc_idx] + tr = nr.transformer3W_map[arc_tuple] + Y11, Y12, Y21, Y22 = ybus_branch_entries(tr) + return [ArcModification( + arc_idx, -b_arc, + YBUS_ELTYPE(-Y11), YBUS_ELTYPE(-Y12), + YBUS_ELTYPE(-Y21), YBUS_ELTYPE(-Y22), + )] elseif tag === :parallel arc_idx = arc_lookup[arc_tuple] b_circuit = PSY.get_series_susceptance(branch) @@ -385,7 +459,7 @@ function _classify_branch_modification( dy11, dy12, dy21, dy22 = _compute_arc_ybus_delta(nr, arc_tuple, delta_b) return [ArcModification(arc_idx, delta_b, dy11, dy12, dy21, dy22)] else - @info "Branch $(PSY.get_name(branch)) not found in any reduction map. " * + @info "Branch $(get_name(branch)) not found in any reduction map. " * "The component may have been eliminated by a radial reduction." return ArcModification[] end From bcfc2a4cb307e9a43788c05dfbc607248194226c Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 14 Apr 2026 19:32:34 +0000 Subject: [PATCH 3/4] Add tests for ThreeWindingTransformer outage support in NetworkModification Agent-Logs-Url: https://github.com/NREL-Sienna/PowerNetworkMatrices.jl/sessions/7950c45a-eba1-4cd1-a6d6-3e4d6238e999 Co-authored-by: jd-lara <16385323+jd-lara@users.noreply.github.com> --- test/test_network_modification.jl | 69 +++++++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) diff --git a/test/test_network_modification.jl b/test/test_network_modification.jl index f3f714502..4ce41ef35 100644 --- a/test/test_network_modification.jl +++ b/test/test_network_modification.jl @@ -250,3 +250,72 @@ end end end end + +@testset "NetworkModification: full ThreeWindingTransformer outage" begin + sys = PSB.build_system(PSB.PSITestSystems, "case10_radial_series_reductions") + trf = first(PSY.get_components(PSY.ThreeWindingTransformer, sys)) + vptdf = VirtualPTDF(sys) + + # Full 3WT outage should produce 3 arc modifications (one per winding) + mod = NetworkModification(vptdf, trf) + @test length(mod.arc_modifications) == 3 + @test mod.label == PSY.get_name(trf) + + # Each modification should have negative delta_b (removing susceptance) + for am in mod.arc_modifications + @test am.delta_b < 0 + end + + # Should produce valid PTDF rows + row = get_post_modification_ptdf_row(vptdf, 1, mod) + @test length(row) == length(PNM.get_bus_axis(vptdf)) +end + +@testset "NetworkModification: single ThreeWindingTransformerWinding outage" begin + sys = PSB.build_system(PSB.PSITestSystems, "case10_radial_series_reductions") + trf = first(PSY.get_components(PSY.ThreeWindingTransformer, sys)) + vptdf = VirtualPTDF(sys) + + # Single winding outage + for w in 1:3 + winding = PNM.ThreeWindingTransformerWinding(trf, w) + mod = NetworkModification(vptdf, winding) + @test length(mod.arc_modifications) == 1 + @test mod.arc_modifications[1].delta_b < 0 + end +end + +@testset "NetworkModification: partial ThreeWindingTransformer (disabled winding)" begin + sys = PSB.build_system(PSB.PSITestSystems, "case10_radial_series_reductions") + trf = first(PSY.get_components(PSY.ThreeWindingTransformer, sys)) + + # Disable one winding before building the matrix + PSY.set_available_secondary!(trf, false) + vptdf = VirtualPTDF(sys) + + # Full 3WT outage should only produce 2 mods (secondary is unavailable) + mod = NetworkModification(vptdf, trf) + @test length(mod.arc_modifications) == 2 +end + +@testset "NetworkModification: ThreeWindingTransformer via Outage attribute" begin + sys = PSB.build_system(PSB.PSITestSystems, "case10_radial_series_reductions") + trf = first(PSY.get_components(PSY.ThreeWindingTransformer, sys)) + + # Attach an outage supplemental attribute to the 3WT + outage = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 0.0, + outage_transition_probability = 0.0, + ) + add_supplemental_attribute!(sys, trf, outage) + vptdf = VirtualPTDF(sys) + + # Construct modification via Outage path + mod = NetworkModification(vptdf, sys, outage) + @test length(mod.arc_modifications) == 3 + @test !isempty(mod.label) + + # Should produce valid PTDF rows + row = get_post_modification_ptdf_row(vptdf, 1, mod) + @test length(row) == length(PNM.get_bus_axis(vptdf)) +end From 4b8589277682a76d0ec439126481c11ab1e84017 Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Fri, 29 May 2026 12:21:20 -0600 Subject: [PATCH 4/4] updates for 3W buses --- src/network_modification.jl | 52 +++++++++++++- test/test_3wt_outage.jl | 137 ++++++++++++++++++++++++++++++++++++ 2 files changed, 188 insertions(+), 1 deletion(-) create mode 100644 test/test_3wt_outage.jl diff --git a/src/network_modification.jl b/src/network_modification.jl index 6da00b3d6..fb8496081 100644 --- a/src/network_modification.jl +++ b/src/network_modification.jl @@ -142,9 +142,45 @@ function NetworkModification( return NetworkModification( PSY.get_name(branch), mods, + ShuntModification[], + _3wt_real_bus_islanding(mat, mods), ) end +""" + _3wt_real_bus_islanding(mat, mods) -> Bool + +True if outaging the 3WT winding arcs `mods` disconnects a real (non-star) bus +from its reference. A full 3WT outage always isolates the fictitious star bus, +which is benign (it carries no injection); this flags only the case where a +load/gen-bearing terminal also disconnects. The star bus is kept in the +union-find so it connects its terminals in the baseline, but is excluded from +the component count — a higher post-outage count means a real bus was islanded. +""" +function _3wt_real_bus_islanding( + mat::PowerNetworkMatrix, + mods::Vector{ArcModification}, +) + isempty(mods) && return false + arc_ax = get_arc_axis(mat) + bus_lookup = get_bus_lookup(mat) + nbus = length(get_bus_axis(mat)) + star_idx = bus_lookup[arc_ax[mods[1].arc_index][2]] + removed = Set(m.arc_index for m in mods) + + uf_before = collect(1:nbus) + uf_after = collect(1:nbus) + for (e, arc) in enumerate(arc_ax) + f = bus_lookup[arc[1]] + t = bus_lookup[arc[2]] + union_sets!(uf_before, f, t) + e in removed || union_sets!(uf_after, f, t) + end + + _count(uf) = length(Set(get_representative(uf, b) for b in 1:nbus if b != star_idx)) + return _count(uf_after) > _count(uf_before) +end + """ $(TYPEDSIGNATURES) @@ -219,9 +255,23 @@ function NetworkModification( outage_uuid = IS.get_uuid(outage) ctg_name = isempty(component_names) ? string(outage_uuid) : join(component_names, "+") - return NetworkModification(ctg_name, mods, shunt_mods, false) + + # A fully-outaged ThreeWindingTransformer isolates its star bus and may + # island a real terminal bus; flag that on `is_islanding`. + is_island = false + arc_ax = get_arc_axis(mat) + for component in all_components + _is_three_winding_transformer(component) || continue + star_num = get_arc_tuple(ThreeWindingTransformerWinding(component, 1))[2] + t3w_mods = [m for m in direct_mods if arc_ax[m.arc_index][2] == star_num] + is_island = is_island || _3wt_real_bus_islanding(mat, t3w_mods) + end + return NetworkModification(ctg_name, mods, shunt_mods, is_island) end +_is_three_winding_transformer(::Any) = false +_is_three_winding_transformer(::PSY.ThreeWindingTransformer) = true + """ _classify_outage_component!(nr, arc_lookup, arc_sus, bus_lookup, component, ...) -> nothing diff --git a/test/test_3wt_outage.jl b/test/test_3wt_outage.jl new file mode 100644 index 000000000..f680b74df --- /dev/null +++ b/test/test_3wt_outage.jl @@ -0,0 +1,137 @@ +import LinearAlgebra as LA + +@testset "full 3WT outage matches deflated-ABA reference (case10)" begin + sys = PSB.build_system(PSB.PSITestSystems, "case10_radial_series_reductions") + vptdf = VirtualPTDF(sys) + bus_ax = PNM.get_bus_axis(vptdf) + arc_ax = PNM.get_arc_axis(vptdf) + BA = vptdf.BA + A = vptdf.A + valid = vptdf.valid_ix + nbus = length(bus_ax) + + trf = first(PSY.get_components(PSY.ThreeWindingTransformer, sys)) + mod = NetworkModification(vptdf, trf) + wind = sort([am.arc_index for am in mod.arc_modifications]) + + # Ground truth: post-outage ABA over valid buses (winding arcs removed), + # solved with the Moore-Penrose pseudo-inverse, which is the convention the + # production pinv-based Woodbury path implements for islanding. + BAm = Matrix(BA) + for e in wind + BAm[:, e] .= 0.0 + end + ABAm = (BAm * A)[valid, valid] + mon = first(setdiff(1:length(arc_ax), wind)) + rhs = Vector(BA[:, mon])[valid] + x_ref = LA.pinv(Matrix(ABAm); atol = 1e-8) * rhs + ref = zeros(nbus) + ref[valid] .= x_ref + + row = get_post_modification_ptdf_row(vptdf, mon, mod) + + # Compare only on buses still connected post-outage (non-zero ABA row). + keep = [valid[k] for k in 1:length(valid) if !all(iszero, ABAm[k, :])] + @test maximum(abs.(row[keep] .- ref[keep])) < 1e-7 +end + +@testset "star bus column is ~zero in post-mod row" begin + sys = PSB.build_system(PSB.PSITestSystems, "case10_radial_series_reductions") + vptdf = VirtualPTDF(sys) + bus_ax = PNM.get_bus_axis(vptdf) + arc_ax = PNM.get_arc_axis(vptdf) + trf = first(PSY.get_components(PSY.ThreeWindingTransformer, sys)) + mod = NetworkModification(vptdf, trf) + wind = sort([am.arc_index for am in mod.arc_modifications]) + star_num = unique([arc_ax[e][2] for e in wind])[1] + star_pos = findfirst(==(star_num), bus_ax) + mon = first(setdiff(1:length(arc_ax), wind)) + row = get_post_modification_ptdf_row(vptdf, mon, mod) + # Isolated zero-injection star bus carries no sensitivity. + @test isapprox(row[star_pos], 0.0; atol = 1e-9) +end + +@testset "is_islanding is true for radial-terminal full 3WT outage" begin + sys = PSB.build_system(PSB.PSITestSystems, "case10_radial_series_reductions") + vptdf = VirtualPTDF(sys) + trf = first(PSY.get_components(PSY.ThreeWindingTransformer, sys)) + mod = NetworkModification(vptdf, trf) + @test mod.is_islanding == true +end + +@testset "is_islanding via Outage attribute path" begin + sys = PSB.build_system(PSB.PSITestSystems, "case10_radial_series_reductions") + trf = first(PSY.get_components(PSY.ThreeWindingTransformer, sys)) + outage = PSY.GeometricDistributionForcedOutage(; + mean_time_to_recovery = 0.0, + outage_transition_probability = 0.0, + ) + PSY.add_supplemental_attribute!(sys, trf, outage) + vptdf = VirtualPTDF(sys) + mod = NetworkModification(vptdf, sys, outage) + @test mod.is_islanding == true +end + +@testset "outaged winding arc returns zero row" begin + sys = PSB.build_system(PSB.PSITestSystems, "case10_radial_series_reductions") + vptdf = VirtualPTDF(sys) + trf = first(PSY.get_components(PSY.ThreeWindingTransformer, sys)) + mod = NetworkModification(vptdf, trf) + wind = [am.arc_index for am in mod.arc_modifications] + for e in wind + row = get_post_modification_ptdf_row(vptdf, e, mod) + @test all(iszero, row) + end +end + +@testset "VirtualMODF cache path equals one-shot (full 3WT)" begin + sys = PSB.build_system(PSB.PSITestSystems, "case10_radial_series_reductions") + vptdf = VirtualPTDF(sys) + vmodf = VirtualMODF(sys) + trf = first(PSY.get_components(PSY.ThreeWindingTransformer, sys)) + mod = NetworkModification(vptdf, trf) + arc_ax = PNM.get_arc_axis(vptdf) + wind = Set(am.arc_index for am in mod.arc_modifications) + mon = first(setdiff(1:length(arc_ax), wind)) + oneshot = get_post_modification_ptdf_row(vptdf, mon, mod) + cached = vmodf[mon, mod] + @test isapprox(oneshot, cached; atol = 1e-9) +end + +@testset "Ybus delta matches star-surgery rebuild on surviving buses" begin + sys = PSB.build_system(PSB.PSITestSystems, "case10_radial_series_reductions") + ybus0 = Ybus(sys) + vptdf = VirtualPTDF(sys) + trf = first(PSY.get_components(PSY.ThreeWindingTransformer, sys)) + mod = NetworkModification(vptdf, trf) + arc_ax = PNM.get_arc_axis(vptdf) + wind = [am.arc_index for am in mod.arc_modifications] + star_num = unique([arc_ax[e][2] for e in wind])[1] + + # Raw post-modification Ybus (original-system positions, star row/col zeroed). + ybus_mod = apply_ybus_modification(ybus0, mod) + lookup_mod = PNM.get_bus_lookup(ybus0) + + # Star-surgery rebuild: remove 3WT + orphan star bus. + sys2 = deepcopy(sys) + trf2 = first(PSY.get_components(PSY.ThreeWindingTransformer, sys2)) + PSY.remove_component!(sys2, trf2) + star_bus = first(b for b in PSY.get_components(PSY.ACBus, sys2) + if PSY.get_number(b) == star_num) + PSY.remove_component!(sys2, star_bus) + ybus_ref = Ybus(sys2) + ref_data = ybus_ref.data + lookup_ref = PNM.get_bus_lookup(ybus_ref) + + # Compare entries for every real bus pair present in the surgery system. + bus_ax_ref = PNM.get_bus_axis(ybus_ref) + for i in bus_ax_ref, j in bus_ax_ref + # ComplexF32 entries: compare with float32-appropriate relative tol. + @test isapprox( + ybus_mod[lookup_mod[i], lookup_mod[j]], + ref_data[lookup_ref[i], lookup_ref[j]]; + atol = 1e-4, + rtol = 1e-4, + ) + end +end