Support ThreeWindingTransformer outages in NetworkModification

src/network_modification.jl

@@ -113,16 +113,74 @@ 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)
     mods = _classify_branch_modification(nr, arc_lookup, arc_sus, branch)
     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)
 
@@ -197,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
 
@@ -238,11 +310,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 +338,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 +407,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 +450,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 +482,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 +509,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

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