diff --git a/docs/src/how_to_guides/compute_network_matrices.md b/docs/src/how_to_guides/compute_network_matrices.md index cea182ed3..b103fcd00 100644 --- a/docs/src/how_to_guides/compute_network_matrices.md +++ b/docs/src/how_to_guides/compute_network_matrices.md @@ -79,6 +79,21 @@ vlodf_matrix = PNM.VirtualLODF(sys) vlodf_matrix[(1, 2), (3, 4)] ``` +## Computing Virtual MODF Matrix + +For post-contingency / post-modification PTDF rows using [`VirtualMODF`](@ref): + +```julia +vmodf_matrix = PNM.VirtualMODF(sys) + +# Inspect contingencies auto-registered from PSY.Outage supplemental attributes +PNM.get_registered_contingencies(vmodf_matrix) + +# Access the post-modification PTDF row for monitored arc (1, 2) under a contingency +contingency = first(values(PNM.get_registered_contingencies(vmodf_matrix))) +vmodf_matrix[(1, 2), contingency] +``` + ## Computing Incidence and BA Matrices For the fundamental network topology matrices: @@ -134,6 +149,7 @@ For matrices involving branches (IncidenceMatrix, PTDF, LODF), branches are repr | `Ybus` | Bus numbers | Bus numbers | | `VirtualPTDF` | Arc tuples | Bus numbers | | `VirtualLODF` | Arc tuples | Arc tuples | +| `VirtualMODF` | Arc tuples | Bus numbers | !!! note diff --git a/docs/src/reference/network_matrices_overview.md b/docs/src/reference/network_matrices_overview.md index b3af6d2ef..b04d560a7 100644 --- a/docs/src/reference/network_matrices_overview.md +++ b/docs/src/reference/network_matrices_overview.md @@ -93,6 +93,19 @@ LODF is derived from PTDF using matrix operations that model the effect of remov See [`VirtualLODF`](@ref) for cases where it is not possible to compute or store the full LODF matrix. +### Modification Distribution Factors ([`VirtualMODF`](@ref)) + +The MODF answers: "Under a registered contingency or network modification (one or more branch outages, partial-susceptance changes, or shunt changes), what is the resulting PTDF row for a given monitored arc?" + +**Key Characteristics:** + + - On-demand computation of post-modification PTDF rows via the Woodbury matrix identity over a base PTDF. + - Caches Woodbury factors per modification and post-modification PTDF rows per `(monitored_arc, modification)` pair. + - Auto-registers `PSY.Outage` supplemental attributes attached to the source system. + - Rows are indexed by monitored arc tuples or arc indices; the second index is a [`ContingencySpec`](@ref), a [`NetworkModification`](@ref), or a `PSY.Outage`. + +**Use it when:** you need post-contingency or post-modification flow sensitivities for one or a few contingencies on a large system, and constructing the full LODF or recomputing PTDF per contingency is too expensive. + ## Arc-Based Indexing All matrices that involve branches use **arc tuples** as identifiers instead of branch name strings. An arc tuple is a `Tuple{Int, Int}` of the form `(from_bus_number, to_bus_number)`, representing the directed connection between two buses. @@ -119,6 +132,7 @@ Each matrix stores `axes` and `lookup` fields: | `Ybus` | Bus numbers | Bus numbers | | `VirtualPTDF` | Arc tuples | Bus numbers | | `VirtualLODF` | Arc tuples | Arc tuples | +| `VirtualMODF` | Arc tuples | Bus numbers | !!! note diff --git a/src/virtual_lodf_calculations.jl b/src/virtual_lodf_calculations.jl index d6f4be258..322c146e2 100644 --- a/src/virtual_lodf_calculations.jl +++ b/src/virtual_lodf_calculations.jl @@ -17,8 +17,9 @@ every libklu solve runs under `_LIBKLU_LOCK` (process-wide) and the per-cache callers can issue requests concurrently; the libklu work runs one at a time. # Arguments -- `K::KLULinSolveCache{Float64}`: - ABA factorization (single cache; solves serialize via the locks above). +- `K`: + ABA factorization cache (`KLULinSolveCache{Float64}` or `AAFactorization`). + Solves serialize via the locks above. - `BA::SparseArrays.SparseMatrixCSC{Float64, Int}`: BA matrix. - `A::SparseArrays.SparseMatrixCSC{Int8, Int}`: @@ -109,8 +110,24 @@ function Base.show(io::IO, ::MIME{Symbol("text/plain")}, array::VirtualLODF) return end +# Map bus index (1:n_buses) to its position in `valid_ix`; ref buses map to 0. +function _build_bus_to_valid_idx(n_buses::Int, valid_ix::Vector{Int})::Vector{Int} + bus_to_valid = zeros(Int, n_buses) + for (vi, bus_ix) in enumerate(valid_ix) + bus_to_valid[bus_ix] = vi + end + return bus_to_valid +end + +""" + _get_PTDF_A_diag(K, BA, A, ref_bus_positions) -> Vector{Float64} + +Compute `diag(PTDF · A)`. Each row of `A` has exactly two nonzeros (+1 at the +from-bus, -1 at the to-bus), so the per-arc dot product reduces to two indexed +reads into the solved PTDF row after a one-time transpose of `A`. +""" function _get_PTDF_A_diag( - K::KLULinSolveCache{Float64}, + K, BA::SparseArrays.SparseMatrixCSC{Float64, Int}, A::SparseArrays.SparseMatrixCSC{Int8, Int}, ref_bus_positions::Set{Int}, @@ -119,35 +136,61 @@ function _get_PTDF_A_diag( n_buses = size(BA, 1) diag_ = zeros(n_branches) - # Pre-compute valid indices (non-reference buses) valid_ix = setdiff(1:n_buses, ref_bus_positions) n_valid = length(valid_ix) + bus_to_valid_idx = _build_bus_to_valid_idx(n_buses, valid_ix) + + # Per-arc (from_valid, to_valid) via one transpose of A; 0 = ref bus. + A_T = SparseArrays.sparse(transpose(A)) + arc_from_valid = Vector{Int}(undef, n_branches) + arc_to_valid = Vector{Int}(undef, n_branches) + at_rv = SparseArrays.rowvals(A_T) + at_nz = SparseArrays.nonzeros(A_T) + for i in 1:n_branches + f_valid = 0 + t_valid = 0 + @inbounds for k in SparseArrays.nzrange(A_T, i) + bus_ix = at_rv[k] + v = at_nz[k] + valid_i = bus_to_valid_idx[bus_ix] + if v > 0 + f_valid = valid_i + elseif v < 0 + t_valid = valid_i + end + end + arc_from_valid[i] = f_valid + arc_to_valid[i] = t_valid + end - # Pre-allocate work arrays for efficiency ba_col = zeros(n_valid) - ptdf_row = zeros(n_buses) + ba_rv = SparseArrays.rowvals(BA) + ba_nz = SparseArrays.nonzeros(BA) - # For each branch, compute PTDF row and dot with incidence column for i in 1:n_branches - # Extract BA column for valid indices fill!(ba_col, 0.0) - for idx in 1:n_valid - bus_idx = valid_ix[idx] - ba_col[idx] = BA[bus_idx, i] + @inbounds for k in SparseArrays.nzrange(BA, i) + valid_i = bus_to_valid_idx[ba_rv[k]] + valid_i > 0 || continue + ba_col[valid_i] = ba_nz[k] end - solve!(K, ba_col) + _solve_factorization(K, ba_col) - # Map back to full bus indices - fill!(ptdf_row, 0.0) - for idx in 1:n_valid - ptdf_row[valid_ix[idx]] = ba_col[idx] + # ba_col is now PTDF row i in valid-index space; H[e,e] = ptdf[from] - ptdf[to]. + f = arc_from_valid[i] + t = arc_to_valid[i] + v_f = if f > 0 + ba_col[f] + else + 0.0 end - - # Compute diagonal element: sum of PTDF[i,j] * A[i,j] for all buses j - for j in 1:n_buses - diag_[i] += ptdf_row[j] * A[i, j] + v_t = if t > 0 + ba_col[t] + else + 0.0 end + @inbounds diag_[i] = v_f - v_t end return diag_ end diff --git a/src/virtual_modf_calculations.jl b/src/virtual_modf_calculations.jl index cbe596087..29c2ba7d5 100644 --- a/src/virtual_modf_calculations.jl +++ b/src/virtual_modf_calculations.jl @@ -30,7 +30,8 @@ cache and skips the recomputation. - `A::SparseArrays.SparseMatrixCSC{Int8, Int}`: Incidence matrix. - `PTDF_A_diag::Vector{Float64}`: - Raw diagonal elements of PTDF·A (H[e,e] values). + Diagonal of `PTDF·A` (H[e,e] values). Lazily populated on the first + read of `vmodf.PTDF_A_diag`; empty until then. - `arc_susceptances::Vector{Float64}`: Effective susceptance for each arc. - `branch_susceptances_by_arc::Vector{Vector{Float64}}`: @@ -111,6 +112,49 @@ get_bus_axis(mat::VirtualMODF) = mat.axes[2] get_tol(mat::VirtualMODF) = mat.tol[] get_system_uuid(M::VirtualMODF) = M.system_uuid +function Base.getproperty(vmodf::VirtualMODF, name::Symbol) + name === :PTDF_A_diag && return _get_ptdf_a_diag_lazy!(vmodf) + return getfield(vmodf, name) +end + +# Use `getfield` for every field read here — `vmodf.PTDF_A_diag` would +# recurse via `getproperty`. +function _get_ptdf_a_diag_lazy!(vmodf::VirtualMODF) + diag = getfield(vmodf, :PTDF_A_diag) + !isempty(diag) && return diag + @lock getfield(vmodf, :solver_lock) begin + diag = getfield(vmodf, :PTDF_A_diag) + !isempty(diag) && return diag + n_arcs = length(getfield(vmodf, :axes)[1]) + @info "Computing VirtualMODF.PTDF_A_diag on first access ($n_arcs arcs)." + t0 = time_ns() + K = getfield(vmodf, :K) + BA = getfield(vmodf, :BA) + A = getfield(vmodf, :A) + ref_bus_set = _ref_bus_positions(vmodf) + new_diag = _get_PTDF_A_diag(K, BA, A, ref_bus_set) + resize!(diag, length(new_diag)) + copyto!(diag, new_diag) + elapsed = (time_ns() - t0) / 1e9 + @info "Computed VirtualMODF.PTDF_A_diag in $(round(elapsed; digits = 2)) s (cached)." + return diag + end +end + +function _ref_bus_positions(vmodf::VirtualMODF) + n_buses = length(getfield(vmodf, :axes)[2]) + valid_ix = getfield(vmodf, :valid_ix) + return Set{Int}(setdiff(1:n_buses, valid_ix)) +end + +""" +$(TYPEDSIGNATURES) + +Return `H[e, e]` for each arc `e`. Same as `vmodf.PTDF_A_diag`; both trigger +the lazy compute on first call and return the cached vector thereafter. +""" +get_PTDF_A_diag(vmodf::VirtualMODF) = vmodf.PTDF_A_diag + # Woodbury kernel accessors. The Woodbury kernel always goes through # `with_solver` so it picks up the `solver_lock` + `_LIBKLU_LOCK` chain. _get_BA(m::VirtualMODF) = m.BA @@ -242,8 +286,8 @@ function VirtualMODF( valid_ix = setdiff(1:length(bus_ax), ref_bus_positions) - # PTDF diagonal precompute runs serially on the dispatcher thread. - PTDF_A_diag = _get_PTDF_A_diag(K, BA.data, A.data, Set(ref_bus_positions)) + # Empty: populated lazily on first read of `vmodf.PTDF_A_diag`. + PTDF_A_diag = Float64[] arc_susceptances = _extract_arc_susceptances(BA.data) branch_susceptances_by_arc = _extract_branch_susceptances_by_arc(BA.data, arc_ax, Ymatrix.network_reduction_data) diff --git a/test/test_modf_lodf_reductions.jl b/test/test_modf_lodf_reductions.jl index de89ab299..68d577df0 100644 --- a/test/test_modf_lodf_reductions.jl +++ b/test/test_modf_lodf_reductions.jl @@ -62,6 +62,90 @@ function verify_modf_lodf_identity( return end +""" +Verify the N-2 simultaneous-outage identity for VirtualMODF. + +For full outage of arcs `e1` and `e2` simultaneously, checks: + + VirtualMODF[m, ctg] ≈ PTDF[m,:] + eff1*PTDF[e1,:] + eff2*PTDF[e2,:] + +for every monitored arc m, where eff1/eff2 come from the 2×2 LODF +correction (Woodbury identity on the pre-contingency LODF matrix): + + denom = 1 - LODF[e1,e2] * LODF[e2,e1] + eff1 = (LODF[m,e1] + LODF[e2,e1] * LODF[m,e2]) / denom + eff2 = (LODF[e1,e2] * LODF[m,e1] + LODF[m,e2]) / denom + +Raises an error if the pair {e1, e2} is an N-2 islanding pair +(denom ≈ 0), since the closed-form formula is undefined in that case. +""" +function verify_modf_n2_lodf_identity( + vmodf::VirtualMODF, + vlodf::VirtualLODF, + ptdf::PTDF, + arc_idx1::Int, + arc_idx2::Int; + atol = 1e-6, +) + n_arcs = size(vlodf, 1) + b_arc1 = vmodf.arc_susceptances[arc_idx1] + b_arc2 = vmodf.arc_susceptances[arc_idx2] + + L12 = vlodf[arc_idx1, arc_idx2] + L21 = vlodf[arc_idx2, arc_idx1] + denom = 1 - L12 * L21 + abs(denom) < 1e-6 && + error( + "Arc pair ($arc_idx1, $arc_idx2) is an N-2 islanding pair (denom=$denom); use a non-islanding pair", + ) + + ctg_uuid = Base.UUID(UInt128(hash((arc_idx1, arc_idx2, :n2)))) + ctg = ContingencySpec( + ctg_uuid, + NetworkModification( + "n2_test_$(arc_idx1)_$(arc_idx2)", + [ArcModification(arc_idx1, -b_arc1), ArcModification(arc_idx2, -b_arc2)], + ), + ) + vmodf.contingency_cache[ctg_uuid] = ctg + + for m in 1:n_arcs + Lm1 = vlodf[m, arc_idx1] + Lm2 = vlodf[m, arc_idx2] + eff1 = (Lm1 + L21 * Lm2) / denom + eff2 = (L12 * Lm1 + Lm2) / denom + + modf_row = vmodf[m, ctg] + expected = ptdf[m, :] .+ eff1 .* ptdf[arc_idx1, :] .+ eff2 .* ptdf[arc_idx2, :] + @test isapprox(modf_row, expected; atol = atol) + end + + empty!(vmodf.woodbury_cache) + return +end + +""" +Return two distinct non-islanding arcs from `branch_map` whose simultaneous +outage does not island the network (L12*L21 ≠ 1 in the pre-contingency LODF). +Raises an error if no such pair exists in the map. +""" +function _find_two_non_islanding_arcs(vmodf, vlodf, branch_map) + candidates = Tuple{Any, Int}[] + for (arc_tuple, _) in branch_map + arc_idx = vmodf.lookup[1][arc_tuple] + if abs(vmodf.PTDF_A_diag[arc_idx]) < 1.0 - 1e-6 + push!(candidates, (arc_tuple, arc_idx)) + end + end + for i in eachindex(candidates), j in (i + 1):length(candidates) + _, idx1 = candidates[i] + _, idx2 = candidates[j] + denom = 1 - vlodf[idx1, idx2] * vlodf[idx2, idx1] + abs(denom) >= 1e-6 && return candidates[i], candidates[j] + end + error("No non-islanding arc pair found in map ($(length(candidates)) candidates)") +end + """ Find a representative non-islanding arc from a branch map. @@ -185,3 +269,16 @@ end verify_modf_lodf_identity(vmodf, vlodf, ptdf, arc_idx, delta_b) end end + +@testset "MODF vs LODF N-2: direct branch pair" begin + sys = PSB.build_system(PSSEParsingTestSystems, "psse_14_network_reduction_test_system") + reductions = NetworkReduction[] + ptdf = PTDF(sys; network_reductions = reductions) + vlodf = VirtualLODF(sys; network_reductions = reductions) + vmodf = VirtualMODF(sys; network_reductions = reductions) + nrd = get_network_reduction_data(vmodf) + + (_, arc_idx1), (_, arc_idx2) = + _find_two_non_islanding_arcs(vmodf, vlodf, nrd.direct_branch_map) + verify_modf_n2_lodf_identity(vmodf, vlodf, ptdf, arc_idx1, arc_idx2) +end diff --git a/test/test_virtual_lodf.jl b/test/test_virtual_lodf.jl index 1c70d05d6..7301f3b70 100644 --- a/test/test_virtual_lodf.jl +++ b/test/test_virtual_lodf.jl @@ -147,3 +147,45 @@ end end @test test_value end + +@testset "_get_PTDF_A_diag: matches reference implementation" begin + # Reference: full-bus dot product. Slow, unambiguously correct. + function _reference_ptdf_a_diag(K, BA, A, ref_bus_positions::Set{Int}) + n_branches = size(BA, 2) + n_buses = size(BA, 1) + valid_ix = setdiff(1:n_buses, ref_bus_positions) + n_valid = length(valid_ix) + diag_ = zeros(n_branches) + ba_col = zeros(n_valid) + ptdf_row = zeros(n_buses) + for i in 1:n_branches + fill!(ba_col, 0.0) + for idx in 1:n_valid + bus_idx = valid_ix[idx] + ba_col[idx] = BA[bus_idx, i] + end + PNM._solve_factorization(K, ba_col) + fill!(ptdf_row, 0.0) + for idx in 1:n_valid + ptdf_row[valid_ix[idx]] = ba_col[idx] + end + for j in 1:n_buses + diag_[i] += ptdf_row[j] * A[i, j] + end + end + return diag_ + end + + for case in ("c_sys5", "c_sys14") + sys = PSB.build_system(PSB.PSITestSystems, case) + A = PNM.IncidenceMatrix(sys) + BA = PNM.BA_Matrix(sys) + ref_pos = Set(PNM.get_ref_bus_position(A)) + ABA = PNM.calculate_ABA_matrix(A.data, BA.data, ref_pos) + K = PNM.klu_factorize(ABA) + + fast = PNM._get_PTDF_A_diag(K, BA.data, A.data, ref_pos) + ref = _reference_ptdf_a_diag(K, BA.data, A.data, ref_pos) + @test fast ≈ ref atol = 1e-12 rtol = 0 + end +end diff --git a/test/test_virtual_modf.jl b/test/test_virtual_modf.jl index 4dc9464c8..40f4dbbfd 100644 --- a/test/test_virtual_modf.jl +++ b/test/test_virtual_modf.jl @@ -468,3 +468,181 @@ end end end end + +@testset "VirtualMODF: N-2 post-contingency PTDF matches N-2 LODF correction" begin + sys5 = PSB.build_system(PSB.PSITestSystems, "c_sys5") + vlodf = VirtualLODF(sys5) + ptdf_ref = PTDF(sys5) + vmodf = VirtualMODF(sys5) + + n_arcs = size(vlodf, 1) + + # N-2 expected row from pre-contingency PTDF + LODF values. + # Derived from Woodbury/Sherman-Morrison on the 2×2 susceptance update: + # denom = 1 - LODF[e1,e2]*LODF[e2,e1] + # eff1 = (LODF[m,e1] + LODF[e2,e1]*LODF[m,e2]) / denom + # eff2 = (LODF[e1,e2]*LODF[m,e1] + LODF[m,e2]) / denom + function n2_lodf_expected(m, e1, e2) + Lm1 = vlodf[m, e1] + Lm2 = vlodf[m, e2] + L12 = vlodf[e1, e2] + L21 = vlodf[e2, e1] + denom = 1 - L12 * L21 + eff1 = (Lm1 + L21 * Lm2) / denom + eff2 = (L12 * Lm1 + Lm2) / denom + return ptdf_ref[m, :] .+ eff1 .* ptdf_ref[e1, :] .+ eff2 .* ptdf_ref[e2, :] + end + + # Bridge arcs (PTDF_A_diag ≈ 1) produce N-1 islanding; skip them. + is_bridge(e) = abs(vmodf.PTDF_A_diag[e]) >= 1.0 - 1e-6 + + # Two non-bridge arcs can still form an N-2 island (bridge pair) where + # removing both disconnects the network. This is detected by + # L12*L21 ≈ 1 (the 2×2 Woodbury denominator collapses to zero). + # In that regime VirtualMODF uses a pseudoinverse while the closed-form + # LODF formula is undefined; skip those pairs. + n2_is_islanding(e1, e2) = abs(1 - vlodf[e1, e2] * vlodf[e2, e1]) < 1e-6 + + uuid_counter = UInt128(10_000) + for e1 in 1:n_arcs + is_bridge(e1) && continue + for e2 in (e1 + 1):n_arcs + is_bridge(e2) && continue + n2_is_islanding(e1, e2) && continue + + b_e1 = vmodf.arc_susceptances[e1] + b_e2 = vmodf.arc_susceptances[e2] + + ctg_uuid = Base.UUID(uuid_counter) + uuid_counter += 1 + ctg = ContingencySpec( + ctg_uuid, + NetworkModification( + "n2_outage_$(e1)_$(e2)", + [ArcModification(e1, -b_e1), ArcModification(e2, -b_e2)], + ), + ) + vmodf.contingency_cache[ctg_uuid] = ctg + + for m in 1:n_arcs + modf_row = PNM._compute_modf_entry(vmodf, m, ctg.modification) + expected = n2_lodf_expected(m, e1, e2) + @test isapprox(modf_row, expected; atol = 1e-6) + end + + # Clear Woodbury cache between contingencies to avoid stale entries. + empty!(vmodf.woodbury_cache) + end + end +end + +@testset "VirtualMODF: PTDF_A_diag is lazy, logs on first access, caches" begin + sys5 = PSB.build_system(PSB.PSITestSystems, "c_sys5") + vmodf = VirtualMODF(sys5) + n_arcs = length(PNM.get_arc_axis(vmodf)) + + # `getfield` bypasses the `getproperty` hook under test. + @test isempty(getfield(vmodf, :PTDF_A_diag)) + + diag1 = @test_logs (:info, r"Computing.*PTDF_A_diag.*first access") ( + :info, r"Computed.*PTDF_A_diag", + ) vmodf.PTDF_A_diag + @test length(diag1) == n_arcs + @test !isempty(getfield(vmodf, :PTDF_A_diag)) + + # Second read is a cache hit: same identity, no further logs. + diag2 = @test_logs min_level = Logging.Info vmodf.PTDF_A_diag + @test diag2 === diag1 + @test PNM.get_PTDF_A_diag(vmodf) === diag1 + + # Cross-check against VirtualLODF (which still computes eagerly). + vlodf = VirtualLODF(sys5) + @test diag1 ≈ vlodf.PTDF_A_diag atol = 1e-10 +end + +@testset "VirtualMODF: N-2 non-bridge islanding pair — connected subnetwork matches rebuilt PTDF" begin + # Real-world flowgate scenario: two individually non-critical lines (neither is + # a bridge on its own) whose simultaneous loss isolates bus 222 in the RTS-GMLC + # network. Lines "B34" (221–222) and "B30" (217–222) are the only connections + # to bus 222; removing both makes it electrically isolated. + # + # The 2×2 Woodbury matrix is singular (L12·L21 ≈ 1), so VirtualMODF falls back + # to its pseudoinverse path. We verify correctness against a freshly built PTDF + # with both lines disabled, skipping the isolated bus column because the + # pseudoinverse does not force that column to zero. + sys = PSB.build_system(PSB.PSITestSystems, "test_RTS_GMLC_sys") + vmodf = VirtualMODF(sys) + + arc_ax = PNM.get_arc_axis(vmodf) + bus_ax = PNM.get_bus_axis(vmodf) + arc_lookup = PNM.get_arc_lookup(vmodf) + n_arcs = length(arc_ax) + + line_b34 = PSY.get_component(PSY.ACBranch, sys, "B34") + line_b30 = PSY.get_component(PSY.ACBranch, sys, "B30") + + function arc_idx_for_line(line) + arc = PSY.get_arc(line) + return arc_lookup[(arc.from.number, arc.to.number)] + end + + e1 = arc_idx_for_line(line_b34) # arc (221, 222) + e2 = arc_idx_for_line(line_b30) # arc (217, 222) + + b_e1 = vmodf.arc_susceptances[e1] + b_e2 = vmodf.arc_susceptances[e2] + ctg_uuid = Base.UUID(UInt128(88_000)) + ctg = ContingencySpec( + ctg_uuid, + NetworkModification( + "n2_island_B34_B30", + [ArcModification(e1, -b_e1), ArcModification(e2, -b_e2)], + ), + ) + vmodf.contingency_cache[ctg_uuid] = ctg + + # Rebuild the PTDF with both lines disabled — gold-standard ground truth. + sys_mod = PSB.build_system(PSB.PSITestSystems, "test_RTS_GMLC_sys") + for name in ("B34", "B30") + PSY.set_available!(PSY.get_component(PSY.ACBranch, sys_mod, name), false) + end + ptdf_rebuilt = PTDF(sys_mod) + rebuilt_arc_lookup = PNM.get_arc_lookup(ptdf_rebuilt) + rebuilt_bus_lookup = PNM.get_bus_lookup(ptdf_rebuilt) + + # Connected buses: every bus appearing in at least one surviving arc. + # Bus 222 is absent because both its lines are outaged. + connected_buses = Set{Int}() + for (fb, tb) in keys(rebuilt_arc_lookup) + push!(connected_buses, fb) + push!(connected_buses, tb) + end + + for m in 1:n_arcs + modf_row = vmodf[m, ctg] + + # Pseudoinverse must produce a finite result even for islanding events. + @test all(isfinite, modf_row) + + monitored_arc = arc_ax[m] + + # Outaged arcs: full loss of flow sensitivity → zero row. + if !haskey(rebuilt_arc_lookup, monitored_arc) + @test all(abs.(modf_row) .< 1e-8) + continue + end + + # Surviving arcs: match rebuilt PTDF for every electrically connected bus. + rebuilt_m = rebuilt_arc_lookup[monitored_arc] + for (b_idx, bus_num) in enumerate(bus_ax) + bus_num ∈ connected_buses || continue + @test isapprox( + modf_row[b_idx], + ptdf_rebuilt[rebuilt_m, rebuilt_bus_lookup[bus_num]]; + atol = 1e-6, + ) + end + end + + PNM.clear_caches!(vmodf) +end