From 76e49718f23a01eae9861e2455aa9716f2a20abc Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Thu, 19 Feb 2026 14:04:27 -0700 Subject: [PATCH 1/2] implement split ybus --- src/PowerNetworkMatrices.jl | 2 + src/Ybus_split.jl | 572 ++++++++++++++++++++++++++ test/test_BA_ABA_matrix.jl | 2 +- test/test_arc_types_and_reductions.jl | 4 +- test/test_ybus_split.jl | 63 +++ 5 files changed, 641 insertions(+), 2 deletions(-) create mode 100644 src/Ybus_split.jl create mode 100644 test/test_ybus_split.jl diff --git a/src/PowerNetworkMatrices.jl b/src/PowerNetworkMatrices.jl index 8454bc643..c271cc590 100644 --- a/src/PowerNetworkMatrices.jl +++ b/src/PowerNetworkMatrices.jl @@ -29,6 +29,7 @@ export validate_connectivity export VirtualLODF export VirtualPTDF export Ybus +export YbusSplit export ArcAdmittanceMatrix export DC_ABA_Matrix_Factorized export DC_ABA_Matrix_Unfactorized @@ -91,6 +92,7 @@ include("ReductionContainer.jl") include("NetworkReductionData.jl") include("ArcAdmittanceMatrix.jl") include("Ybus.jl") +include("Ybus_split.jl") include("IncidenceMatrix.jl") include("AdjacencyMatrix.jl") include("connectivity_checks.jl") diff --git a/src/Ybus_split.jl b/src/Ybus_split.jl new file mode 100644 index 000000000..89ae0cef2 --- /dev/null +++ b/src/Ybus_split.jl @@ -0,0 +1,572 @@ +const TransformerRecord = NamedTuple{ + (:from, :to, :turns, :angle, :Y, :shunt), + Tuple{Int, Int, Float32, Float32, ComplexF32, ComplexF32}, +} + +struct YbusSplit{Ax, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{ComplexF32} + y_series::SparseArrays.SparseMatrixCSC{ComplexF32, Int} + y_shunt::Vector{ComplexF32} + y_lambda::SparseArrays.SparseMatrixCSC{ComplexF32, Int} + adjacency_data::SparseArrays.SparseMatrixCSC{Int8, Int} + axes::Ax + lookup::L + subnetwork_axes::Dict{Int, Ax} + arc_subnetwork_axis::Dict{Int, Vector{Tuple{Int, Int}}} + network_reduction_data::NetworkReductionData + # need record of transformer data: (from, to, turns, angle, Y, shunt) for each. + # store separately (6 arrays), or together (one array of 6-tuples)? + transformers::Vector{TransformerRecord} +end + +get_y_series(M::YbusSplit) = M.y_series +get_y_shunt(M::YbusSplit) = M.y_shunt +get_axes(M::YbusSplit) = M.axes +get_lookup(M::YbusSplit) = M.lookup +get_ref_bus(M::YbusSplit) = sort!(collect(keys(M.subnetwork_axes))) +get_ref_bus_position(M::YbusSplit) = + [get_bus_lookup(M)[x] for x in keys(M.subnetwork_axes)] +get_network_reduction_data(M::YbusSplit) = M.network_reduction_data +get_bus_axis(M::YbusSplit) = M.axes[1] +get_bus_lookup(M::YbusSplit) = M.lookup[1] + +# --- Branch entry routines (series-only + separate shunts) --- + +"""Series-only Ybus entries for a `Line` or `DiscreteControlledACBranch`.""" +function ybus_split_branch_entries(br::PSY.ACTransmission) + Y_l = (1 / (PSY.get_r(br) + PSY.get_x(br) * 1im)) + Y11 = Y_l + if !isfinite(Y11) || !isfinite(Y_l) + error( + "Data in $(PSY.get_name(br)) is incorrect. r = $(PSY.get_r(br)), x = $(PSY.get_x(br))", + ) + end + Y12 = -Y_l + Y21 = Y12 + Y22 = Y_l + return (Y11, Y12, Y21, Y22, _get_shunt(br, :from), _get_shunt(br, :to)) +end + +# this covers phase shifting, tap transformers, and plain 2W (fixed winding number) +"""Structurally nonzero but numerically zero entries for two-winding transformers.""" +function ybus_split_branch_entries(::PSY.TwoWindingTransformer) + z = zero(ComplexF32) + return (z, z, z, z, z, z) +end + +function ybus_split_branch_entries(parallel_br::BranchesParallel) + Y11 = Y12 = Y21 = Y22 = Yshunt_from = Yshunt_to = zero(ComplexF32) + for br in parallel_br + # All branches in BranchesParallel have the same orientation when constructed in add_to_branch_maps! + (y11, y12, y21, y22, yshunt_from, yshunt_to) = ybus_split_branch_entries(br) + Y11 += y11 + Y12 += y12 + Y21 += y21 + Y22 += y22 + Yshunt_from += yshunt_from + Yshunt_to += yshunt_to + end + return (Y11, Y12, Y21, Y22, Yshunt_from, Yshunt_to) +end + +function ybus_split_branch_entries(::BranchesSeries) + error( + "Equivalent Circuit Formulation techniques aren't compatible with degree 2 reductions.", + ) +end + +# --- Storing branch entries into COO vectors --- + +function add_branch_entries_to_ybus_split!( + y11::Vector{ComplexF32}, + y12::Vector{ComplexF32}, + y21::Vector{ComplexF32}, + y22::Vector{ComplexF32}, + branch_ix::Int, + y_shunt_from::Vector{ComplexF32}, + y_shunt_to::Vector{ComplexF32}, + br::PSY.ACTransmission, +) + Y11, Y12, Y21, Y22, Yshunt_from, Yshunt_to = ybus_split_branch_entries(br) + y11[branch_ix] = Y11 + y12[branch_ix] = Y12 + y21[branch_ix] = Y21 + y22[branch_ix] = Y22 + y_shunt_from[branch_ix] = Yshunt_from + y_shunt_to[branch_ix] = Yshunt_to + return +end + +# --- Per-branch dispatch (indexing + admittance) --- + +function _add_to_transformer_records!( + transformer_records::Vector{TransformerRecord}, + comp::PSY.TwoWindingTransformer, + from::Int, + to::Int, +) + push!( + transformer_records, + ( + from = from, + to = to, + turns = _get_tap(comp), + angle = PSY.get_α(comp), + Y = 1 / (PSY.get_r(comp) + PSY.get_x(comp) * 1im), + shunt = PSY.get_primary_shunt(comp), + ), + ) + return +end + +_add_to_transformer_records!( + ::Vector{TransformerRecord}, + ::PSY.ACTransmission, + ::Int, + ::Int, +) = nothing + +function _ybus_split!( + y11::Vector{ComplexF32}, + y12::Vector{ComplexF32}, + y21::Vector{ComplexF32}, + y22::Vector{ComplexF32}, + y_shunt_from::Vector{ComplexF32}, + y_shunt_to::Vector{ComplexF32}, + br::PSY.ACTransmission, + num_bus::Dict{Int, Int}, + branch_ix::Int, + fb::Vector{Int}, + tb::Vector{Int}, + nr::NetworkReductionData, + adj::SparseArrays.SparseMatrixCSC{Int8, Int}, + transformer_records::Vector{TransformerRecord}, +) + add_branch_entries_to_indexing_maps!(num_bus, branch_ix, nr, fb, tb, adj, br) + add_branch_entries_to_ybus_split!( + y11, y12, y21, y22, branch_ix, y_shunt_from, y_shunt_to, br, + ) + _add_to_transformer_records!(transformer_records, br, fb[branch_ix], tb[branch_ix]) + return +end + +function _ybus_split!( + y11::Vector{ComplexF32}, + y12::Vector{ComplexF32}, + y21::Vector{ComplexF32}, + y22::Vector{ComplexF32}, + y_shunt_from::Vector{ComplexF32}, + y_shunt_to::Vector{ComplexF32}, + br::PSY.DynamicBranch, + num_bus::Dict{Int, Int}, + branch_ix::Int, + fb::Vector{Int}, + tb::Vector{Int}, + nr::NetworkReductionData, + adj::SparseArrays.SparseMatrixCSC{Int8, Int}, + transformer_records::Vector{TransformerRecord}, +) + _ybus_split!( + y11, y12, y21, y22, y_shunt_from, y_shunt_to, + br.branch, num_bus, branch_ix, fb, tb, nr, adj, + transformer_records, + ) + return +end + +function _ybus_split!( + br::PSY.ThreeWindingTransformer, + num_bus::Dict{Int, Int}, + offset_ix::Int, + fb::Vector{Int}, + tb::Vector{Int}, + ix::Int, + nr::NetworkReductionData, + adj::SparseArrays.SparseMatrixCSC{Int8, Int}, + transformer_records::Vector{TransformerRecord}, +) + primary_star_arc = PSY.get_primary_star_arc(br) + secondary_star_arc = PSY.get_secondary_star_arc(br) + tertiary_star_arc = PSY.get_tertiary_star_arc(br) + add_to_branch_maps!(nr, primary_star_arc, secondary_star_arc, tertiary_star_arc, br) + # All admittance values are zero for 3W transformers in the split formulation; + # only the structural entries (fb/tb) and adjacency need to be set. + n_entries = 0 + all_windings_available = + PSY.get_available_primary(br) && PSY.get_available_secondary(br) && + PSY.get_available_tertiary(br) + all_windings_available || error( + "YbusSplit doesn't support 3WT's with unavailable windings. Found: $(PSY.get_name(br))", + ) + # primary + primary_ix, star_ix = get_bus_indices(primary_star_arc, num_bus, nr) + adj[primary_ix, star_ix] = 1 + adj[star_ix, primary_ix] = -1 + entry_ix = offset_ix + ix + n_entries + fb[entry_ix] = primary_ix + tb[entry_ix] = star_ix + n_entries += 1 + push!( + transformer_records, + ( + from = primary_ix, + to = star_ix, + turns = PSY.get_primary_turns_ratio(br), + angle = PSY.get_α_primary(br), + Y = 1 / (PSY.get_r_primary(br) + PSY.get_x_primary(br) * 1im), + shunt = PSY.get_g(br) + im * PSY.get_b(br), + ), + ) + # secondary + secondary_ix, star_ix = get_bus_indices(secondary_star_arc, num_bus, nr) + adj[secondary_ix, star_ix] = 1 + adj[star_ix, secondary_ix] = -1 + entry_ix = offset_ix + ix + n_entries + fb[entry_ix] = secondary_ix + tb[entry_ix] = star_ix + n_entries += 1 + push!( + transformer_records, + ( + from = secondary_ix, + to = star_ix, + turns = PSY.get_secondary_turns_ratio(br), + angle = PSY.get_α_secondary(br), + Y = 1 / (PSY.get_r_secondary(br) + PSY.get_x_secondary(br) * 1im), + shunt = zero(ComplexF32), + ), + ) + # tertiary + tertiary_ix, star_ix = get_bus_indices(tertiary_star_arc, num_bus, nr) + adj[tertiary_ix, star_ix] = 1 + adj[star_ix, tertiary_ix] = -1 + entry_ix = offset_ix + ix + n_entries + fb[entry_ix] = tertiary_ix + tb[entry_ix] = star_ix + n_entries += 1 + push!( + transformer_records, + ( + from = tertiary_ix, + to = star_ix, + turns = PSY.get_tertiary_turns_ratio(br), + angle = PSY.get_α_tertiary(br), + Y = 1 / (PSY.get_r_tertiary(br) + PSY.get_x_tertiary(br) * 1im), + shunt = zero(ComplexF32), + ), + ) + return n_entries +end + +# --- Shunt device handlers --- + +function _ybus_split_shunt!( + ysh::Vector{ComplexF32}, + fa::PSY.FixedAdmittance, + num_bus::Dict{Int, Int}, + fa_ix::Int, + sb::Vector{Int}, + nr::NetworkReductionData, +) + bus_no = get_bus_index(fa, num_bus, nr) + Y = PSY.get_Y(fa) + sb[fa_ix] = bus_no + if !isfinite(Y) + error( + "Data in $(PSY.get_name(fa)) is incorrect. Y = $(Y)", + ) + end + ysh[fa_ix] = Y + return +end + +function _ybus_split_shunt!( + ysh::Vector{ComplexF32}, + fa::PSY.SwitchedAdmittance, + num_bus::Dict{Int, Int}, + fa_ix::Int, + sb::Vector{Int}, + nr::NetworkReductionData, +) + bus_no = get_bus_index(fa, num_bus, nr) + sb[fa_ix] = bus_no + ysh[fa_ix] = 0.0 + return +end + +# --- Core assembly --- + +function _buildybus_split!( + network_reduction_data::NetworkReductionData, + adj::SparseArrays.SparseMatrixCSC{Int8, Int}, + branches, + transformer_3w::Vector{PSY.ThreeWindingTransformer}, + num_bus::Dict{Int, Int}, + fixed_admittances::Vector{PSY.FixedAdmittance}, + switched_admittances::Vector{PSY.SwitchedAdmittance}, + transformer_records::Vector{TransformerRecord}, +) + branch_entries_transformer_3w = 0 + for br in transformer_3w + branch_entries_transformer_3w += count([ + PSY.get_available_primary(br), + PSY.get_available_secondary(br), + PSY.get_available_tertiary(br), + ]) + end + branchcount = length(branches) + branch_entries_transformer_3w + branchcount_no_3w = length(branches) + fa_count = length(fixed_admittances) + sa_count = length(switched_admittances) + fb = zeros(Int, branchcount) + tb = zeros(Int, branchcount) + sb = zeros(Int, fa_count + sa_count) + + y11 = zeros(ComplexF32, branchcount) + y12 = zeros(ComplexF32, branchcount) + y21 = zeros(ComplexF32, branchcount) + y22 = zeros(ComplexF32, branchcount) + y_shunt_from = zeros(ComplexF32, branchcount) + y_shunt_to = zeros(ComplexF32, branchcount) + ysh = zeros(ComplexF32, fa_count + sa_count) + + for (ix, b) in enumerate(branches) + if PSY.get_name(b) == "init" + throw(DataFormatError("The data in Branch is invalid")) + end + _ybus_split!( + y11, y12, y21, y22, y_shunt_from, y_shunt_to, + b, num_bus, ix, fb, tb, network_reduction_data, adj, + transformer_records, + ) + end + + ix = 1 + for b in transformer_3w + if PSY.get_name(b) == "init" + throw(DataFormatError("The data in Transformer3W is invalid")) + end + n_entries = _ybus_split!(b, num_bus, branchcount_no_3w, fb, tb, ix, + network_reduction_data, adj, transformer_records, + ) + ix += n_entries + end + for (ix, fa) in enumerate([fixed_admittances; switched_admittances]) + _ybus_split_shunt!(ysh, fa, num_bus, ix, sb, network_reduction_data) + end + return ( + y11, + y12, + y21, + y22, + y_shunt_from, + y_shunt_to, + ysh, + fb, + tb, + sb, + ) +end + +# --- Constructor --- + +function YbusSplit( + sys::PSY.System; + subnetwork_algorithm = iterative_union_find, + network_reductions = NetworkReduction[], +) + units_base = PSY.get_units_base(sys) + if units_base != "SYSTEM_BASE" + @warn "Setting the system unit base from $units_base to SYSTEM_BASE for matrix construction" + PSY.set_units_base_system!(sys, "SYSTEM_BASE") + end + + if !isempty(network_reductions) + error("YbusSplit doesn't support network reductions.") + end + + # Prototype guards: error on unsupported component types + for b in PSY.get_available_components(PSY.ACBus, sys) + if PSY.get_bustype(b) == ACBusTypes.ISOLATED + error("YbusSplit does not support isolated buses. Found: $(PSY.get_name(b))") + end + end + if !isempty( + collect( + PSY.get_available_components(PSY.DiscreteControlledACBranch, sys), + ), + ) + error("YbusSplit does not support DiscreteControlledACBranch components.") + end + if !isempty( + collect( + PSY.get_available_components(PSY.StandardLoad, sys), + ), + ) + error("YbusSplit does not support StandardLoad components.") + end + + ref_bus_numbers = Set{Int}() + nr = NetworkReductionData() + bus_reduction_map = get_bus_reduction_map(nr) + + for b in PSY.get_available_components(PSY.ACBus, sys) + bus_reduction_map[PSY.get_number(b)] = Set{Int}() + if PSY.get_bustype(b) == ACBusTypes.REF + push!(ref_bus_numbers, PSY.get_number(b)) + end + end + + bus_ax = sort!(collect(keys(bus_reduction_map))) + axes = (bus_ax, bus_ax) + bus_lookup = Dict{Int, Int}() + lookup = (bus_lookup, bus_lookup) + busnumber = length(bus_ax) + for (ix, b) in enumerate(bus_ax) + bus_lookup[b] = ix + end + adj = SparseArrays.spdiagm(ones(Int8, busnumber)) + branches = collect( + PSY.get_components( + x -> + PSY.get_available(x) && + typeof(x) ∉ [ + PSY.Transformer3W, + PSY.PhaseShiftingTransformer3W, + ], + PSY.ACTransmission, + sys, + ), + ) + transformer_3W = + collect( + PSY.get_available_components(PSY.ThreeWindingTransformer, sys), + ) + fixed_admittances = collect( + PSY.get_available_components(PSY.FixedAdmittance, sys), + ) + switched_admittances = + collect( + PSY.get_available_components(PSY.SwitchedAdmittance, sys), + ) + n_2W_transformers = length( + collect( + PSY.get_available_components(PSY.TwoWindingTransformer, sys), + ), + ) + transformer_records = Vector{TransformerRecord}() + sizehint!(transformer_records, n_2W_transformers + 3 * length(transformer_3W)) + y11, y12, y21, y22, y_shunt_from, y_shunt_to, ysh, fb, tb, sb = + _buildybus_split!( + nr, + adj, + branches, + transformer_3W, + bus_lookup, + fixed_admittances, + switched_admittances, + transformer_records, + ) + + # Assemble series-only sparse matrix (no shunts folded into diagonal). + # Transformer entries are structurally present but numerically zero. + y_series = SparseArrays.sparse( + [fb; fb; tb; tb], + [fb; tb; fb; tb], + [y11; y12; y21; y22], + busnumber, + busnumber, + ) + function has_structural_diagonal(A::SparseArrays.SparseMatrixCSC) + rv = SparseArrays.rowvals(A) + for j in 1:size(A, 2) + r = SparseArrays.nzrange(A, j) + idx = searchsortedfirst(rv, j, first(r), last(r), Base.Order.Forward) + (idx > last(r) || rv[idx] != j) && return false + end + return true + end + @assert has_structural_diagonal(y_series) "YbusSplit requires structurally nonzero diagonal entries for connectivity detection. Found zero diagonal entry at bus index $(findfirst(x -> x == 0, diag(y_series)))" + + I, J, V = SparseArrays.findnz(y_series) + y_lambda = SparseArrays.sparse( + I, J, zeros(ComplexF32, length(V)), + ) + + # Accumulate per-bus shunt vector: line charging from branches + device shunts. + y_shunt = zeros(ComplexF32, busnumber) + for i in eachindex(fb) + y_shunt[fb[i]] += y_shunt_from[i] + y_shunt[tb[i]] += y_shunt_to[i] + end + for i in eachindex(sb) + y_shunt[sb[i]] += ysh[i] + end + + if length(bus_lookup) > 1 + # Use the combined matrix for connectivity detection + subnetworks = assign_reference_buses!( + find_subnetworks( + y_series + SparseArrays.spdiagm(y_shunt), + bus_ax; + subnetwork_algorithm = subnetwork_algorithm, + ), + ref_bus_numbers, + ) + if length(subnetworks) > 1 + @warn "More than one island found; Network is not connected" + end + else + subnetworks = Dict{Int, Set{Int}}(bus_lookup[1] => Set(bus_ax)) + end + subnetwork_axes = _make_bus_subnetwork_axes(subnetworks) + arc_subnetwork_axis = _make_arc_subnetwork_axis(subnetworks, nr) + return YbusSplit( + y_series, + y_shunt, + y_lambda, + adj, + axes, + lookup, + subnetwork_axes, + arc_subnetwork_axis, + nr, + transformer_records, + ) +end + +function generic_transfomer_ybus_entries(turns::Number, angle::Number, Y::Complex) + tap = turns * exp(1im * angle) + Y11 = Y / abs2(tap) + Y12 = -Y / conj(tap) + Y21 = -Y / tap + Y22 = Y + return (Y11, Y12, Y21, Y22) +end + +function update_y_lambda!(M::YbusSplit, lambda::Number, gamma::Number) + nonzeros(M.y_lambda) .= 0.0 + # transmission piece: (1 + λ γ) * Y_series + (1 - λ) * Diagonal(Y_shunt) + # λ = 0 gives original, Y_series + Diagonal(Y_shunt) + # λ = 1 gives (1 + γ) * Y_series, shunts removed and all series entries shorted. + @assert rowvals(M.y_lambda) == rowvals(M.y_series) + @assert colptr(M.y_lambda) == colptr(M.y_series) + nonzeros(M.y_lambda) .+= nonzeros(M.y_series) .* (1 + lambda * gamma) + for j in 1:size(M.y_lambda, 2) + M.y_lambda[j, j] += M.y_shunt[j] * (1 - lambda) + end + # transformer piece: + for record in M.transformers + from = record.from + to = record.to + # λ = 0 gives original + # λ = 1 gives turns = 1, angle = 0, admittance = Y (1 + γ) ~ infinity, shunt = 0. + turns = record.turns + lambda * (1 - record.turns) + angle = (1 - lambda) * record.angle + Y = (1 + lambda * gamma) * record.Y + Y11, Y12, Y21, Y22 = generic_transfomer_ybus_entries(turns, angle, Y) + M.y_lambda[from, from] += Y11 + record.shunt * (1 - lambda) + M.y_lambda[from, to] += Y12 + M.y_lambda[to, from] += Y21 + M.y_lambda[to, to] += Y22 + end +end diff --git a/test/test_BA_ABA_matrix.jl b/test/test_BA_ABA_matrix.jl index f5f32a16b..d86841fb7 100644 --- a/test/test_BA_ABA_matrix.jl +++ b/test/test_BA_ABA_matrix.jl @@ -12,7 +12,7 @@ # check the is_factorized function @test is_factorized(ABA_no_lu) == false # factorize the ABA matrix - ABA_no_lu = factorize(ABA_no_lu) + ABA_no_lu = PNM.factorize(ABA_no_lu) @test is_factorized(ABA_no_lu) == true # get the ABA matrix with the current method ABA_lu = ABA_Matrix(sys; factorize = true) diff --git a/test/test_arc_types_and_reductions.jl b/test/test_arc_types_and_reductions.jl index 820135c91..7601f19ba 100644 --- a/test/test_arc_types_and_reductions.jl +++ b/test/test_arc_types_and_reductions.jl @@ -29,7 +29,7 @@ end function test_all_subtypes(sys::System, network_reductions) for T in subtypes(PNM.PowerNetworkMatrix) # arc admittance matrix constructor has different args. - (T == PNM.Ybus || T == PNM.ArcAdmittanceMatrix) && continue + (T == PNM.Ybus || T == PNM.ArcAdmittanceMatrix || T == PNM.YbusSplit) && continue M = T(sys; network_reductions = deepcopy(network_reductions)) # test that it runs without error @test M isa T @@ -37,6 +37,8 @@ function test_all_subtypes(sys::System, network_reductions) # return Ybus so we can inspect the network reduction data. Y = Ybus(sys; network_reductions = deepcopy(network_reductions)) @test Y isa PNM.Ybus + #Y_split = YbusSplit(sys) + #@test Y_split isa PNM.YbusSplit return Y end diff --git a/test/test_ybus_split.jl b/test/test_ybus_split.jl new file mode 100644 index 000000000..333b5fd37 --- /dev/null +++ b/test/test_ybus_split.jl @@ -0,0 +1,63 @@ +using SparseArrays +using LinearAlgebra + +@testset "YbusSplit" begin + + # --- No transformers: y_series + Diagonal(y_shunt) should equal the normal Ybus --- + @testset "No transformers: recombined split equals Ybus" begin + sys = PSB.build_system(PSB.PSITestSystems, "c_sys5") + ybus = Ybus(sys) + ybus_split = YbusSplit(sys) + + y_recombined = ybus_split.y_series + spdiagm(ybus_split.y_shunt) + # dropzeros so the comparison is clean + dropzeros!(y_recombined) + + @test size(ybus.data) == size(y_recombined) + I_ref, J_ref, V_ref = findnz(ybus.data) + I_split, J_split, V_split = findnz(y_recombined) + @test Set(zip(I_ref, J_ref)) == Set(zip(I_split, J_split)) + for (i, j, v) in zip(I_ref, J_ref, V_ref) + @test isapprox(y_recombined[i, j], v; atol = 1e-5) + end + end + + # --- With transformers: same sparsity structure, off-diagonal match away from TX buses --- + @testset "With transformers: structure and non-TX entries" begin + sys = build_system(MatpowerTestSystems, "matpower_case5_transformer") + ybus = Ybus(sys) + ybus_split = YbusSplit(sys) + + y_recombined = ybus_split.y_series + spdiagm(ybus_split.y_shunt) + + # Sparsity pattern of y_series should be a superset of the normal Ybus + # (transformer arcs are structurally present but numerically zero in y_series) + I_ref, J_ref, _ = findnz(ybus.data) + ref_pattern = Set(zip(I_ref, J_ref)) + I_split, J_split, _ = findnz(ybus_split.y_series) + split_pattern = Set(zip(I_split, J_split)) + @test ref_pattern ⊆ split_pattern + + # Find buses connected to transformers so we can exclude them + bus_lookup = PNM.get_bus_lookup(ybus) + tx_bus_indices = Set{Int}() + for br in PSY.get_components( + x -> PSY.get_available(x), + PSY.TwoWindingTransformer, + sys, + ) + arc = PSY.get_arc(br) + from_no = PSY.get_number(PSY.get_from(arc)) + to_no = PSY.get_number(PSY.get_to(arc)) + push!(tx_bus_indices, bus_lookup[from_no]) + push!(tx_bus_indices, bus_lookup[to_no]) + end + + # Off-diagonal entries not touching any transformer bus should match exactly + for (i, j, v) in zip(I_ref, J_ref, findnz(ybus.data)[3]) + i == j && continue + (i ∈ tx_bus_indices || j ∈ tx_bus_indices) && continue + @test isapprox(y_recombined[i, j], v; atol = 1e-5) + end + end +end From 6b713d637319e4f5e9b290ab65ff435a9bbcbfc8 Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Tue, 24 Feb 2026 09:12:59 -0700 Subject: [PATCH 2/2] switch to float64 and other tweaks --- src/PowerNetworkMatrices.jl | 1 + src/Ybus_split.jl | 119 ++++++++++++++++++++---------------- 2 files changed, 67 insertions(+), 53 deletions(-) diff --git a/src/PowerNetworkMatrices.jl b/src/PowerNetworkMatrices.jl index c271cc590..2ccb9966d 100644 --- a/src/PowerNetworkMatrices.jl +++ b/src/PowerNetworkMatrices.jl @@ -30,6 +30,7 @@ export VirtualLODF export VirtualPTDF export Ybus export YbusSplit +export update_y_lambda! export ArcAdmittanceMatrix export DC_ABA_Matrix_Factorized export DC_ABA_Matrix_Unfactorized diff --git a/src/Ybus_split.jl b/src/Ybus_split.jl index 89ae0cef2..042f21346 100644 --- a/src/Ybus_split.jl +++ b/src/Ybus_split.jl @@ -1,12 +1,12 @@ const TransformerRecord = NamedTuple{ (:from, :to, :turns, :angle, :Y, :shunt), - Tuple{Int, Int, Float32, Float32, ComplexF32, ComplexF32}, + Tuple{Int, Int, Float64, Float64, ComplexF64, ComplexF64}, } -struct YbusSplit{Ax, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{ComplexF32} - y_series::SparseArrays.SparseMatrixCSC{ComplexF32, Int} - y_shunt::Vector{ComplexF32} - y_lambda::SparseArrays.SparseMatrixCSC{ComplexF32, Int} +struct YbusSplit{Ax, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{ComplexF64} + y_series::SparseArrays.SparseMatrixCSC{ComplexF64, Int} + y_shunt::Vector{ComplexF64} + y_lambda::SparseArrays.SparseMatrixCSC{ComplexF64, Int} adjacency_data::SparseArrays.SparseMatrixCSC{Int8, Int} axes::Ax lookup::L @@ -49,12 +49,12 @@ end # this covers phase shifting, tap transformers, and plain 2W (fixed winding number) """Structurally nonzero but numerically zero entries for two-winding transformers.""" function ybus_split_branch_entries(::PSY.TwoWindingTransformer) - z = zero(ComplexF32) + z = zero(ComplexF64) return (z, z, z, z, z, z) end function ybus_split_branch_entries(parallel_br::BranchesParallel) - Y11 = Y12 = Y21 = Y22 = Yshunt_from = Yshunt_to = zero(ComplexF32) + Y11 = Y12 = Y21 = Y22 = Yshunt_from = Yshunt_to = zero(ComplexF64) for br in parallel_br # All branches in BranchesParallel have the same orientation when constructed in add_to_branch_maps! (y11, y12, y21, y22, yshunt_from, yshunt_to) = ybus_split_branch_entries(br) @@ -77,13 +77,13 @@ end # --- Storing branch entries into COO vectors --- function add_branch_entries_to_ybus_split!( - y11::Vector{ComplexF32}, - y12::Vector{ComplexF32}, - y21::Vector{ComplexF32}, - y22::Vector{ComplexF32}, + y11::Vector{ComplexF64}, + y12::Vector{ComplexF64}, + y21::Vector{ComplexF64}, + y22::Vector{ComplexF64}, branch_ix::Int, - y_shunt_from::Vector{ComplexF32}, - y_shunt_to::Vector{ComplexF32}, + y_shunt_from::Vector{ComplexF64}, + y_shunt_to::Vector{ComplexF64}, br::PSY.ACTransmission, ) Y11, Y12, Y21, Y22, Yshunt_from, Yshunt_to = ybus_split_branch_entries(br) @@ -126,12 +126,12 @@ _add_to_transformer_records!( ) = nothing function _ybus_split!( - y11::Vector{ComplexF32}, - y12::Vector{ComplexF32}, - y21::Vector{ComplexF32}, - y22::Vector{ComplexF32}, - y_shunt_from::Vector{ComplexF32}, - y_shunt_to::Vector{ComplexF32}, + y11::Vector{ComplexF64}, + y12::Vector{ComplexF64}, + y21::Vector{ComplexF64}, + y22::Vector{ComplexF64}, + y_shunt_from::Vector{ComplexF64}, + y_shunt_to::Vector{ComplexF64}, br::PSY.ACTransmission, num_bus::Dict{Int, Int}, branch_ix::Int, @@ -150,12 +150,12 @@ function _ybus_split!( end function _ybus_split!( - y11::Vector{ComplexF32}, - y12::Vector{ComplexF32}, - y21::Vector{ComplexF32}, - y22::Vector{ComplexF32}, - y_shunt_from::Vector{ComplexF32}, - y_shunt_to::Vector{ComplexF32}, + y11::Vector{ComplexF64}, + y12::Vector{ComplexF64}, + y21::Vector{ComplexF64}, + y22::Vector{ComplexF64}, + y_shunt_from::Vector{ComplexF64}, + y_shunt_to::Vector{ComplexF64}, br::PSY.DynamicBranch, num_bus::Dict{Int, Int}, branch_ix::Int, @@ -232,7 +232,7 @@ function _ybus_split!( turns = PSY.get_secondary_turns_ratio(br), angle = PSY.get_α_secondary(br), Y = 1 / (PSY.get_r_secondary(br) + PSY.get_x_secondary(br) * 1im), - shunt = zero(ComplexF32), + shunt = zero(ComplexF64), ), ) # tertiary @@ -251,7 +251,7 @@ function _ybus_split!( turns = PSY.get_tertiary_turns_ratio(br), angle = PSY.get_α_tertiary(br), Y = 1 / (PSY.get_r_tertiary(br) + PSY.get_x_tertiary(br) * 1im), - shunt = zero(ComplexF32), + shunt = zero(ComplexF64), ), ) return n_entries @@ -260,7 +260,7 @@ end # --- Shunt device handlers --- function _ybus_split_shunt!( - ysh::Vector{ComplexF32}, + ysh::Vector{ComplexF64}, fa::PSY.FixedAdmittance, num_bus::Dict{Int, Int}, fa_ix::Int, @@ -280,7 +280,7 @@ function _ybus_split_shunt!( end function _ybus_split_shunt!( - ysh::Vector{ComplexF32}, + ysh::Vector{ComplexF64}, fa::PSY.SwitchedAdmittance, num_bus::Dict{Int, Int}, fa_ix::Int, @@ -321,13 +321,13 @@ function _buildybus_split!( tb = zeros(Int, branchcount) sb = zeros(Int, fa_count + sa_count) - y11 = zeros(ComplexF32, branchcount) - y12 = zeros(ComplexF32, branchcount) - y21 = zeros(ComplexF32, branchcount) - y22 = zeros(ComplexF32, branchcount) - y_shunt_from = zeros(ComplexF32, branchcount) - y_shunt_to = zeros(ComplexF32, branchcount) - ysh = zeros(ComplexF32, fa_count + sa_count) + y11 = zeros(ComplexF64, branchcount) + y12 = zeros(ComplexF64, branchcount) + y21 = zeros(ComplexF64, branchcount) + y22 = zeros(ComplexF64, branchcount) + y_shunt_from = zeros(ComplexF64, branchcount) + y_shunt_to = zeros(ComplexF64, branchcount) + ysh = zeros(ComplexF64, fa_count + sa_count) for (ix, b) in enumerate(branches) if PSY.get_name(b) == "init" @@ -385,11 +385,6 @@ function YbusSplit( end # Prototype guards: error on unsupported component types - for b in PSY.get_available_components(PSY.ACBus, sys) - if PSY.get_bustype(b) == ACBusTypes.ISOLATED - error("YbusSplit does not support isolated buses. Found: $(PSY.get_name(b))") - end - end if !isempty( collect( PSY.get_available_components(PSY.DiscreteControlledACBranch, sys), @@ -410,6 +405,11 @@ function YbusSplit( bus_reduction_map = get_bus_reduction_map(nr) for b in PSY.get_available_components(PSY.ACBus, sys) + if PSY.get_bustype(b) == ACBusTypes.ISOLATED + @debug "Found available isolated bus $(PSY.get_name(b)) with number $(PSY.get_number(b)). This is excluded from the YbusSplit build." + push!(nr.removed_buses, PSY.get_number(b)) + continue + end bus_reduction_map[PSY.get_number(b)] = Set{Int}() if PSY.get_bustype(b) == ACBusTypes.REF push!(ref_bus_numbers, PSY.get_number(b)) @@ -442,11 +442,23 @@ function YbusSplit( PSY.get_available_components(PSY.ThreeWindingTransformer, sys), ) fixed_admittances = collect( - PSY.get_available_components(PSY.FixedAdmittance, sys), + PSY.get_components( + x -> + PSY.get_available(x) && + PSY.get_bustype(PSY.get_bus(x)) != ACBusTypes.ISOLATED, + PSY.FixedAdmittance, + sys, + ), ) switched_admittances = collect( - PSY.get_available_components(PSY.SwitchedAdmittance, sys), + PSY.get_components( + x -> + PSY.get_available(x) && + PSY.get_bustype(PSY.get_bus(x)) != ACBusTypes.ISOLATED, + PSY.SwitchedAdmittance, + sys, + ), ) n_2W_transformers = length( collect( @@ -467,8 +479,6 @@ function YbusSplit( transformer_records, ) - # Assemble series-only sparse matrix (no shunts folded into diagonal). - # Transformer entries are structurally present but numerically zero. y_series = SparseArrays.sparse( [fb; fb; tb; tb], [fb; tb; fb; tb], @@ -489,11 +499,11 @@ function YbusSplit( I, J, V = SparseArrays.findnz(y_series) y_lambda = SparseArrays.sparse( - I, J, zeros(ComplexF32, length(V)), + I, J, zeros(ComplexF64, length(V)), ) # Accumulate per-bus shunt vector: line charging from branches + device shunts. - y_shunt = zeros(ComplexF32, busnumber) + y_shunt = zeros(ComplexF64, busnumber) for i in eachindex(fb) y_shunt[fb[i]] += y_shunt_from[i] y_shunt[tb[i]] += y_shunt_to[i] @@ -503,10 +513,12 @@ function YbusSplit( end if length(bus_lookup) > 1 - # Use the combined matrix for connectivity detection + # Use the adjacency matrix for connectivity detection. + # y_series + spdiagm(y_shunt) would drop structural zeros (transformer + # entries that are zero at lambda != 0), incorrectly splitting the network. subnetworks = assign_reference_buses!( find_subnetworks( - y_series + SparseArrays.spdiagm(y_shunt), + adj, bus_ax; subnetwork_algorithm = subnetwork_algorithm, ), @@ -544,13 +556,14 @@ function generic_transfomer_ybus_entries(turns::Number, angle::Number, Y::Comple end function update_y_lambda!(M::YbusSplit, lambda::Number, gamma::Number) - nonzeros(M.y_lambda) .= 0.0 + SparseArrays.nonzeros(M.y_lambda) .= 0.0 # transmission piece: (1 + λ γ) * Y_series + (1 - λ) * Diagonal(Y_shunt) # λ = 0 gives original, Y_series + Diagonal(Y_shunt) # λ = 1 gives (1 + γ) * Y_series, shunts removed and all series entries shorted. - @assert rowvals(M.y_lambda) == rowvals(M.y_series) - @assert colptr(M.y_lambda) == colptr(M.y_series) - nonzeros(M.y_lambda) .+= nonzeros(M.y_series) .* (1 + lambda * gamma) + @assert SparseArrays.rowvals(M.y_lambda) == SparseArrays.rowvals(M.y_series) + @assert M.y_lambda.colptr == M.y_series.colptr + SparseArrays.nonzeros(M.y_lambda) .+= + SparseArrays.nonzeros(M.y_series) .* (1 + lambda * gamma) for j in 1:size(M.y_lambda, 2) M.y_lambda[j, j] += M.y_shunt[j] * (1 - lambda) end