diff --git a/docs/src/tutorials/tutorial_DFAX.md b/docs/src/tutorials/tutorial_DFAX.md index f1779acd..2f37f64a 100644 --- a/docs/src/tutorials/tutorial_DFAX.md +++ b/docs/src/tutorials/tutorial_DFAX.md @@ -263,7 +263,7 @@ mod_ab3 = NetworkModification(vmodf, (123, 217)); # AB3 outage ctg_n2 = NetworkModification( "AB2_and_AB3_out", vcat(collect(mod_ab2.arc_modifications), - collect(mod_ab3.arc_modifications)), + collect(mod_ab3.arc_modifications)), ); row_n2 = vmodf[m, ctg_n2]; @@ -303,7 +303,7 @@ ctg_ab3 = NetworkModification(vmodf, (123, 217)); ctg_ab2_ab3 = NetworkModification( "AB2_and_AB3_out", vcat(collect(ctg_ab2.arc_modifications), - collect(ctg_ab3.arc_modifications)), + collect(ctg_ab3.arc_modifications)), ); contingencies = [ diff --git a/src/NetworkReductionData.jl b/src/NetworkReductionData.jl index f7c82563..15081137 100644 --- a/src/NetworkReductionData.jl +++ b/src/NetworkReductionData.jl @@ -134,6 +134,7 @@ network reduction algorithms. - `component_to_reduction_name_map::Dict{Type, Dict{String, String}}`: Lazily filled with the call to [`populate_branch_maps_by_type!`](@ref), maps component names to the names of the reduction entries used in name_to_arc_map. - `filters_applied::Dict{Type, Function}`: Filters applied when populating branch maps by type - `direct_branch_name_map::Dict{String, Tuple{Int, Int}}`: Lazily filled, maps branch names to their corresponding arc tuples for direct branches +- `icd_map::Union{Nothing, ICDCorrectionMap}`: Precomputed impedance-correction-data factors for 2W and 3W transformers, built once via [`_build_icd_correction_map`](@ref) """ @kwdef mutable struct NetworkReductionData irreducible_buses::Set{Int} = Set{Int}() # Buses that are not reduced in the network reduction @@ -181,6 +182,7 @@ network reduction algorithms. filters_applied = Dict{Type, Function}() #Filters applied when populating branch maps by type direct_branch_name_map::Dict{String, Tuple{Int, Int}} = Dict{String, Tuple{Int, Int}}() + icd_map::Union{Nothing, ICDCorrectionMap} = nothing end function add_to_map(device::T, filters::Dict) where {T <: PSY.ACTransmission} @@ -485,6 +487,9 @@ has_filtered_branches(rb::NetworkReductionData) = !isempty(rb.filters_applied) function Base.isempty(rb::NetworkReductionData) for field in fieldnames(NetworkReductionData) + # icd_map is a precomputed cache (Union{Nothing, ICDCorrectionMap}), not a + # reduction-state collection; isempty is undefined for both possible types. + field == :icd_map && continue if !isempty(getfield(rb, field)) return false end @@ -494,6 +499,7 @@ end function Base.empty!(rb::NetworkReductionData) for field in fieldnames(NetworkReductionData) + field == :icd_map && continue empty!(getfield(rb, field)) end return @@ -557,6 +563,9 @@ function isequal( # direct_branch_name_map is populated when indexing into matrices with branch names # this should not prevent using matrices for downstream computations (e.g. LODF(A, BA, ABA)) field == :direct_branch_name_map && continue + # icd_map has no `==` override, so two structurally-identical maps from separate + # _build_icd_correction_map calls (or deepcopy) would compare unequal via `===`. + field == :icd_map && continue if getfield(rb1, field) != getfield(rb2, field) return false end diff --git a/src/PowerNetworkMatrices.jl b/src/PowerNetworkMatrices.jl index 6abb526b..62b3c9a2 100644 --- a/src/PowerNetworkMatrices.jl +++ b/src/PowerNetworkMatrices.jl @@ -121,9 +121,9 @@ include("ward_reduction.jl") # `apply_zero_impedance_reduction.jl`, included after `Ybus.jl` below. include("zero_impedance_branch_reduction.jl") include("ReductionContainer.jl") +include("YbusACBranches.jl") include("NetworkReductionData.jl") include("ArcAdmittanceMatrix.jl") -include("YbusACBranches.jl") include("Ybus.jl") include("apply_zero_impedance_reduction.jl") include("IncidenceMatrix.jl") diff --git a/src/Ybus.jl b/src/Ybus.jl index 032b8493..acbe6141 100644 --- a/src/Ybus.jl +++ b/src/Ybus.jl @@ -464,6 +464,30 @@ function ybus_branch_entries( return ybus_branch_entries(br; min_x_eps = min_x_eps) end +function ybus_branch_entries( + br::PSY.TwoWindingTransformer, + nr::NetworkReductionData; + min_x_eps::Float64 = ZERO_IMPEDANCE_X_EPSILON, +) + return ybus_branch_entries( + br, + _get_impedance_correction_factor(br, nr); + min_x_eps = min_x_eps, + ) +end + +function ybus_branch_entries( + tp::ThreeWindingTransformerWinding, + nr::NetworkReductionData; + min_x_eps::Float64 = ZERO_IMPEDANCE_X_EPSILON, +) + return ybus_branch_entries( + tp, + _get_impedance_correction_factor(tp, nr); + min_x_eps = min_x_eps, + ) +end + function ybus_branch_entries( br::PSY.GenericArcImpedance; min_x_eps::Float64 = ZERO_IMPEDANCE_X_EPSILON, @@ -481,6 +505,16 @@ function ybus_branch_entries( return (Y11, Y12, Y21, Y22) end +# A single branch has unambiguous orientation; this `nr` overload lets callers iterating +# heterogeneous segments (single branches and parallel groups) pass `nr` uniformly. +function ybus_branch_entries( + br::PSY.GenericArcImpedance, + ::NetworkReductionData; + min_x_eps::Float64 = ZERO_IMPEDANCE_X_EPSILON, +) + return ybus_branch_entries(br; min_x_eps = min_x_eps) +end + function ybus_branch_entries( parallel_br::AbstractBranchesParallel, nr::NetworkReductionData; @@ -493,7 +527,7 @@ function ybus_branch_entries( # full-precision; callers narrow to YBUS_ELTYPE when storing. Y11 = Y12 = Y21 = Y22 = zero(ComplexF64) for br in parallel_br - (y11, y12, y21, y22) = ybus_branch_entries(br) + (y11, y12, y21, y22) = ybus_branch_entries(br, nr) if get_arc_tuple(br, nr) != reference Y11 += y22 Y12 += y21 @@ -522,12 +556,109 @@ end _get_tap(::PSY.Transformer2W) = one(YBUS_ELTYPE) _get_tap(br::PSY.TwoWindingTransformer) = PSY.get_tap(br) +function _interpolate_correction_factor(curve::IS.PiecewiseLinearData, x::Real) + points = IS.get_points(curve) + x = clamp(x, points[1].x, points[end].x) + for i in 1:(length(points) - 1) + if x <= points[i + 1].x + dx = points[i + 1].x - points[i].x + iszero(dx) && return points[i].y + t = (x - points[i].x) / dx + return points[i].y + t * (points[i + 1].y - points[i].y) + end + end + return points[end].y +end + +function _compute_icd_factor( + br::PSY.TwoWindingTransformer, + ict::PSY.ImpedanceCorrectionData, +) + curve = PSY.get_impedance_correction_curve(ict) + mode = PSY.get_transformer_control_mode(ict) + x = if mode == PSY.ImpedanceCorrectionTransformerControlMode.TAP_RATIO + abs(_get_tap(br)) + else # PHASE_SHIFT_ANGLE — table x-values are in degrees; α is stored in radians + rad2deg(PSY.get_α(br)) + end + return _interpolate_correction_factor(curve, x) +end + +function _compute_icd_factor( + br::PSY.ThreeWindingTransformer, + ict::PSY.ImpedanceCorrectionData, + winding_num::Int, +) + curve = PSY.get_impedance_correction_curve(ict) + mode = PSY.get_transformer_control_mode(ict) + x = if mode == PSY.ImpedanceCorrectionTransformerControlMode.TAP_RATIO + abs( + ( + PSY.get_primary_turns_ratio(br), + PSY.get_secondary_turns_ratio(br), + PSY.get_tertiary_turns_ratio(br), + )[winding_num], + ) + else # PHASE_SHIFT_ANGLE — table x-values are in degrees; angles stored in radians + rad2deg( + ( + PSY.get_α_primary(br), + PSY.get_α_secondary(br), + PSY.get_α_tertiary(br), + )[winding_num], + ) + end + return _interpolate_correction_factor(curve, x) +end + +function _winding_category_to_number(cat::PSY.WindingCategoryModule.WindingCategory)::Int + cat == PSY.WindingCategory.PRIMARY_WINDING && return 1 + cat == PSY.WindingCategory.SECONDARY_WINDING && return 2 + return 3 +end + +function _build_icd_correction_map(sys::PSY.System)::Union{Nothing, ICDCorrectionMap} + pairs_2w = PSY.get_component_supplemental_attribute_pairs( + PSY.TwoWindingTransformer, + PSY.ImpedanceCorrectionData, + sys, + ) + pairs_3w = PSY.get_component_supplemental_attribute_pairs( + PSY.ThreeWindingTransformer, + PSY.ImpedanceCorrectionData, + sys, + ) + isempty(pairs_2w) && isempty(pairs_3w) && return nothing + map_2w = Dict{Base.UUID, Float64}() + for pair in pairs_2w + map_2w[IS.get_uuid(pair.component)] = + _compute_icd_factor(pair.component, pair.supplemental_attribute) + end + map_3w = Dict{Tuple{Base.UUID, Int}, Float64}() + for pair in pairs_3w + wnum = _winding_category_to_number( + PSY.get_transformer_winding(pair.supplemental_attribute), + ) + map_3w[(IS.get_uuid(pair.component), wnum)] = + _compute_icd_factor(pair.component, pair.supplemental_attribute, wnum) + end + return ICDCorrectionMap(map_2w, map_3w) +end + +function _get_impedance_correction_factor( + br::PSY.TwoWindingTransformer, + nr::NetworkReductionData, +) + return _get_correction(nr.icd_map, br) +end + """Ybus entries for a `Transformer2W`, `TapTransformer`, or `PhaseShiftingTransformer`.""" function ybus_branch_entries( - br::PSY.TwoWindingTransformer; + br::PSY.TwoWindingTransformer, + correction::Float64; min_x_eps::Float64 = ZERO_IMPEDANCE_X_EPSILON, ) - Y_t = 1 / (PSY.get_r(br) + PSY.get_x(br) * 1im) + Y_t = 1 / ((PSY.get_r(br) + PSY.get_x(br) * 1im) * correction) tap = _get_tap(br) * exp(PSY.get_α(br) * 1im) c_tap = _get_tap(br) * exp(-1 * PSY.get_α(br) * 1im) y_shunt = PSY.get_primary_shunt(br) @@ -543,15 +674,26 @@ function ybus_branch_entries( return (Y11 + y_shunt, Y12, Y21, Y22) end +function _get_impedance_correction_factor( + tp::ThreeWindingTransformerWinding, + nr::NetworkReductionData, +) + isnothing(nr.icd_map) && return one(Float64) + br = get_transformer(tp) + winding_num = get_winding_number(tp) + return get(nr.icd_map.map_3w, (IS.get_uuid(br), winding_num), 1.0) +end + """Ybus branch entries for an arc in the wye model of a `ThreeWindingTransformer`.""" function ybus_branch_entries( - tp::ThreeWindingTransformerWinding; + tp::ThreeWindingTransformerWinding, + correction::Float64; min_x_eps::Float64 = ZERO_IMPEDANCE_X_EPSILON, ) br = get_transformer(tp) winding_number = get_winding_number(tp) if winding_number == 1 - Y_t = 1 / (PSY.get_r_primary(br) + PSY.get_x_primary(br) * 1im) + Y_t = 1 / ((PSY.get_r_primary(br) + PSY.get_x_primary(br) * 1im) * correction) tap = PSY.get_primary_turns_ratio(br) * exp(PSY.get_α_primary(br) * 1im) c_tap = PSY.get_primary_turns_ratio(br) * exp(-1 * PSY.get_α_primary(br) * 1im) Y11 = Y_t / abs2(tap) @@ -565,9 +707,10 @@ function ybus_branch_entries( # primary bus alone gets the shunt term Y11 += y_shunt elseif winding_number == 2 - Y_t = 1 / (PSY.get_r_secondary(br) + PSY.get_x_secondary(br) * 1im) + Y_t = 1 / ((PSY.get_r_secondary(br) + PSY.get_x_secondary(br) * 1im) * correction) tap = PSY.get_secondary_turns_ratio(br) * exp(PSY.get_α_secondary(br) * 1im) - c_tap = PSY.get_secondary_turns_ratio(br) * exp(-1 * PSY.get_α_secondary(br) * 1im) + c_tap = + PSY.get_secondary_turns_ratio(br) * exp(-1 * PSY.get_α_secondary(br) * 1im) Y11 = Y_t / abs2(tap) if !isfinite(Y11) error( @@ -576,7 +719,7 @@ function ybus_branch_entries( ) end elseif winding_number == 3 - Y_t = 1 / (PSY.get_r_tertiary(br) + PSY.get_x_tertiary(br) * 1im) + Y_t = 1 / ((PSY.get_r_tertiary(br) + PSY.get_x_tertiary(br) * 1im) * correction) tap = PSY.get_tertiary_turns_ratio(br) * exp(PSY.get_α_tertiary(br) * 1im) c_tap = PSY.get_tertiary_turns_ratio(br) * exp(-1 * PSY.get_α_tertiary(br) * 1im) Y11 = Y_t / abs2(tap) @@ -602,6 +745,7 @@ function _ybus!( y21::Vector{YBUS_ELTYPE}, y22::Vector{YBUS_ELTYPE}, br::PSY.ACTransmission, + correction::Float64, num_bus::Dict{Int, Int}, branch_ix::Int, fb::Vector{Int}, @@ -619,30 +763,42 @@ function _ybus!( return end +function _ybus!( + y11::Vector{YBUS_ELTYPE}, + y12::Vector{YBUS_ELTYPE}, + y21::Vector{YBUS_ELTYPE}, + y22::Vector{YBUS_ELTYPE}, + br::PSY.TwoWindingTransformer, + correction::Float64, + num_bus::Dict{Int, Int}, + branch_ix::Int, + fb::Vector{Int}, + tb::Vector{Int}, + nr::NetworkReductionData, +) + add_branch_entries_to_indexing_maps!(num_bus, branch_ix, nr, fb, tb, br) + Y11, Y12, Y21, Y22 = ybus_branch_entries(br, correction) + y11[branch_ix] = Y11 + y12[branch_ix] = Y12 + y21[branch_ix] = Y21 + y22[branch_ix] = Y22 + return +end + function _ybus!( y11::Vector{YBUS_ELTYPE}, y12::Vector{YBUS_ELTYPE}, y21::Vector{YBUS_ELTYPE}, y22::Vector{YBUS_ELTYPE}, br::PSY.DynamicBranch, + correction::Float64, num_bus::Dict{Int, Int}, branch_ix::Int, fb::Vector{Int}, tb::Vector{Int}, nr::NetworkReductionData, ) - _ybus!( - y11, - y12, - y21, - y22, - br.branch, - num_bus, - branch_ix, - fb, - tb, - nr, - ) + _ybus!(y11, y12, y21, y22, br.branch, correction, num_bus, branch_ix, fb, tb, nr) return end @@ -658,6 +814,9 @@ function _ybus!( tb::Vector{Int}, ix::Int, nr::NetworkReductionData, + c1::Float64, + c2::Float64, + c3::Float64, ) primary_star_arc = PSY.get_primary_star_arc(br) secondary_star_arc = PSY.get_secondary_star_arc(br) @@ -671,7 +830,8 @@ function _ybus!( primary_ix, star_ix = get_bus_indices(primary_star_arc, num_bus, nr) fb[offset_ix + ix + n_entries] = primary_ix tb[offset_ix + ix + n_entries] = star_ix - (Y11, Y12, Y21, Y22) = ybus_branch_entries(ThreeWindingTransformerWinding(br, 1)) + (Y11, Y12, Y21, Y22) = + ybus_branch_entries(ThreeWindingTransformerWinding(br, 1), c1) y11[offset_ix + ix + n_entries] = Y11 y12[offset_ix + ix + n_entries] = Y12 y21[offset_ix + ix + n_entries] = Y21 @@ -682,7 +842,8 @@ function _ybus!( secondary_ix, star_ix = get_bus_indices(secondary_star_arc, num_bus, nr) fb[offset_ix + ix + n_entries] = secondary_ix tb[offset_ix + ix + n_entries] = star_ix - (Y11, Y12, Y21, Y22) = ybus_branch_entries(ThreeWindingTransformerWinding(br, 2)) + (Y11, Y12, Y21, Y22) = + ybus_branch_entries(ThreeWindingTransformerWinding(br, 2), c2) y11[offset_ix + ix + n_entries] = Y11 y12[offset_ix + ix + n_entries] = Y12 y21[offset_ix + ix + n_entries] = Y21 @@ -693,7 +854,8 @@ function _ybus!( tertiary_ix, star_ix = get_bus_indices(tertiary_star_arc, num_bus, nr) fb[offset_ix + ix + n_entries] = tertiary_ix tb[offset_ix + ix + n_entries] = star_ix - (Y11, Y12, Y21, Y22) = ybus_branch_entries(ThreeWindingTransformerWinding(br, 3)) + (Y11, Y12, Y21, Y22) = + ybus_branch_entries(ThreeWindingTransformerWinding(br, 3), c3) y11[offset_ix + ix + n_entries] = Y11 y12[offset_ix + ix + n_entries] = Y12 y21[offset_ix + ix + n_entries] = Y21 @@ -766,6 +928,7 @@ function _buildybus!( fixed_admittances::Vector{PSY.FixedAdmittance}, switched_admittances::Vector{PSY.SwitchedAdmittance}, standard_loads::Vector{PSY.StandardLoad}, + icd_map::Union{Nothing, ICDCorrectionMap} = nothing, ) branch_entries_transformer_3w = 0 for br in transformer_3w @@ -789,11 +952,23 @@ function _buildybus!( y22 = zeros(YBUS_ELTYPE, branchcount) ysh = zeros(YBUS_ELTYPE, fa_count + sa_count + sl_count) - _foreach_ybus_branch(branches) do b, ix + _foreach_ybus_branch(branches, icd_map) do b, ix, correction if PSY.get_name(b) == "init" throw(DataFormatError("The data in Branch is invalid")) end - _ybus!(y11, y12, y21, y22, b, num_bus, ix, fb, tb, network_reduction_data) + _ybus!( + y11, + y12, + y21, + y22, + b, + correction, + num_bus, + ix, + fb, + tb, + network_reduction_data, + ) end ix = 1 @@ -801,6 +976,14 @@ function _buildybus!( if PSY.get_name(b) == "init" throw(DataFormatError("The data in Transformer3W is invalid")) end + uuid = IS.get_uuid(b) + c1, c2, c3 = if isnothing(icd_map) + 1.0, 1.0, 1.0 + else + get(icd_map.map_3w, (uuid, 1), 1.0), + get(icd_map.map_3w, (uuid, 2), 1.0), + get(icd_map.map_3w, (uuid, 3), 1.0) + end n_entries = _ybus!( y11, y12, @@ -813,6 +996,9 @@ function _buildybus!( tb, ix, network_reduction_data, + c1, + c2, + c3, ) ix += n_entries end @@ -1005,6 +1191,9 @@ function Ybus( for (ix, b) in enumerate(bus_ax) bus_lookup[b] = ix end + adj = SparseArrays.spdiagm(ones(Int8, busnumber)) + icd_map = _build_icd_correction_map(sys) + nr.icd_map = icd_map branches = _get_ybus_two_terminal_ac_branches(sys) transformer_3W = _get_filtered_components(PSY.ThreeWindingTransformer, sys, PSY.get_available) @@ -1026,6 +1215,7 @@ function Ybus( fixed_admittances, switched_admittances, standard_loads, + icd_map, ) # Build adjacency matrix from COO triplets in a single sparse() call to avoid # ~2×branchcount structural insertions into a growing CSC matrix. @@ -2141,7 +2331,7 @@ function add_segment_to_ybus!( segment_orientation::Symbol, nr::NetworkReductionData, ) - (Y11, Y12, Y21, Y22) = ybus_branch_entries(segment) + (Y11, Y12, Y21, Y22) = ybus_branch_entries(segment, nr) push!(fb, ix) push!(tb, ix + 1) if segment_orientation == :FromTo diff --git a/src/YbusACBranches.jl b/src/YbusACBranches.jl index 7d221dff..f91619c9 100644 --- a/src/YbusACBranches.jl +++ b/src/YbusACBranches.jl @@ -1,3 +1,8 @@ +struct ICDCorrectionMap + map_2w::Dict{Base.UUID, Float64} + map_3w::Dict{Tuple{Base.UUID, Int}, Float64} +end + struct YbusACBranches lines::Vector{PSY.Line} monitored_lines::Vector{PSY.MonitoredLine} @@ -32,6 +37,13 @@ function _populate_ybus_branch_vector!( return end +_get_correction(::Nothing, ::PSY.ACTransmission) = 1.0 +_get_correction(::Nothing, ::PSY.DynamicBranch) = 1.0 +_get_correction(::ICDCorrectionMap, ::PSY.ACTransmission) = 1.0 +_get_correction(map::ICDCorrectionMap, br::PSY.TwoWindingTransformer) = + get(map.map_2w, IS.get_uuid(br), 1.0) +_get_correction(map::ICDCorrectionMap, br::PSY.DynamicBranch) = + get(map.map_2w, IS.get_uuid(br.branch), 1.0) function _populate_ybus_branch_vector!( vec::Vector{T}, sys::PSY.System, @@ -71,25 +83,27 @@ end function _foreach_ybus_branch( f::F, branches::YbusACBranches, + icd_map::Union{Nothing, ICDCorrectionMap}, ) where {F <: Function} - ix = _foreach_typed_branches(f, branches.lines, 0) - ix = _foreach_typed_branches(f, branches.monitored_lines, ix) - ix = _foreach_typed_branches(f, branches.generic_arc_impedances, ix) - ix = _foreach_typed_branches(f, branches.tap_transformers, ix) - ix = _foreach_typed_branches(f, branches.phase_shifting_transformers, ix) - ix = _foreach_typed_branches(f, branches.transformer2w, ix) - ix = _foreach_typed_branches(f, branches.dynamic_branches, ix) - ix = _foreach_typed_branches(f, branches.breaker_switches, ix) + ix = _foreach_typed_branches(f, branches.lines, icd_map, 0) + ix = _foreach_typed_branches(f, branches.monitored_lines, icd_map, ix) + ix = _foreach_typed_branches(f, branches.generic_arc_impedances, icd_map, ix) + ix = _foreach_typed_branches(f, branches.tap_transformers, icd_map, ix) + ix = _foreach_typed_branches(f, branches.phase_shifting_transformers, icd_map, ix) + ix = _foreach_typed_branches(f, branches.transformer2w, icd_map, ix) + ix = _foreach_typed_branches(f, branches.dynamic_branches, icd_map, ix) + ix = _foreach_typed_branches(f, branches.breaker_switches, icd_map, ix) return ix end function _foreach_typed_branches( f::F, vec::Vector{T}, + icd_map::Union{Nothing, ICDCorrectionMap}, offset::Int, ) where {F <: Function, T <: PSY.ACTransmission} for (i, br) in enumerate(vec) - f(br, offset + i) + f(br, offset + i, _get_correction(icd_map, br)) end return offset + length(vec) end diff --git a/src/network_modification.jl b/src/network_modification.jl index f45af29c..047584cd 100644 --- a/src/network_modification.jl +++ b/src/network_modification.jl @@ -26,18 +26,19 @@ function _scaled_pi_model( end """ - _compute_parallel_partial_ybus_delta(bp, delta_b) -> NTuple{4, YBUS_ELTYPE} + _compute_parallel_partial_ybus_delta(bp, nr, delta_b) -> NTuple{4, YBUS_ELTYPE} Compute the Pi-model Ybus delta for a partial outage on a parallel branch group. Finds the circuit whose susceptance matches `delta_b` and returns its negated Pi-model entries. """ function _compute_parallel_partial_ybus_delta( bp::AbstractBranchesParallel, + nr::NetworkReductionData, delta_b::Float64, )::NTuple{4, YBUS_ELTYPE} for br in bp.branches if _is_full_outage(delta_b, get_series_susceptance(br)) - return _negated_pi_model(ybus_branch_entries(br)) + return _negated_pi_model(ybus_branch_entries(br, nr)) end end error( @@ -49,10 +50,11 @@ end # Direct arc: full outage negates the Pi-model; otherwise scale it by `delta_b / b_arc`. function _direct_arc_ybus_delta( br::PSY.ACTransmission, + nr::NetworkReductionData, delta_b::Float64, )::NTuple{4, YBUS_ELTYPE} b_arc = get_series_susceptance(br) - entries = ybus_branch_entries(br) + entries = ybus_branch_entries(br, nr) if _is_full_outage(delta_b, b_arc) return _negated_pi_model(entries) end @@ -68,7 +70,7 @@ function _parallel_arc_ybus_delta( if _is_full_outage(delta_b, get_series_susceptance(bp)) return _negated_pi_model(ybus_branch_entries(bp, nr)) end - return _compute_parallel_partial_ybus_delta(bp, delta_b) + return _compute_parallel_partial_ybus_delta(bp, nr, delta_b) end # Series-reduced arc: only full outage of the equivalent chain is supported. @@ -90,6 +92,7 @@ end # 3-winding transformer arc: only full outage of the winding is supported. function _transformer3W_arc_ybus_delta( tr::ThreeWindingTransformerWinding, + nr::NetworkReductionData, arc_tuple::Tuple{Int, Int}, delta_b::Float64, )::NTuple{4, YBUS_ELTYPE} @@ -99,7 +102,7 @@ function _transformer3W_arc_ybus_delta( "Arc $(arc_tuple), Δb=$(delta_b).", ) end - return _negated_pi_model(ybus_branch_entries(tr)) + return _negated_pi_model(ybus_branch_entries(tr, nr)) end """ @@ -114,7 +117,7 @@ function _compute_arc_ybus_delta( delta_b::Float64, )::NTuple{4, YBUS_ELTYPE} if haskey(nr.direct_branch_map, arc_tuple) - return _direct_arc_ybus_delta(nr.direct_branch_map[arc_tuple], delta_b) + return _direct_arc_ybus_delta(nr.direct_branch_map[arc_tuple], nr, delta_b) elseif haskey(nr.parallel_branch_map, arc_tuple) return _parallel_arc_ybus_delta(nr.parallel_branch_map[arc_tuple], nr, delta_b) elseif haskey(nr.series_branch_map, arc_tuple) @@ -127,6 +130,7 @@ function _compute_arc_ybus_delta( elseif haskey(nr.transformer3W_map, arc_tuple) return _transformer3W_arc_ybus_delta( nr.transformer3W_map[arc_tuple], + nr, arc_tuple, delta_b, ) diff --git a/test/performance/performance_test.jl b/test/performance/performance_test.jl index daa41111..ff13062a 100644 --- a/test/performance/performance_test.jl +++ b/test/performance/performance_test.jl @@ -79,7 +79,10 @@ for (group, name) in systems end else try - record_samples("$(name)-Build VirtualPTDF", sample_times(() -> VirtualPTDF(sys))) + record_samples( + "$(name)-Build VirtualPTDF", + sample_times(() -> VirtualPTDF(sys)), + ) vptdf = VirtualPTDF(sys) n_arcs = length(vptdf.axes[1]) n_query = min(10, n_arcs) diff --git a/test/test_arc_types_and_reductions.jl b/test/test_arc_types_and_reductions.jl index 1129ca94..38bd7942 100644 --- a/test/test_arc_types_and_reductions.jl +++ b/test/test_arc_types_and_reductions.jl @@ -155,7 +155,7 @@ end # ybus_branch_entries on the mixed group should equal the sum of the parts. Y11_l, Y12_l, Y21_l, Y22_l = PNM.ybus_branch_entries(line) - Y11_t, Y12_t, Y21_t, Y22_t = PNM.ybus_branch_entries(tap) + Y11_t, Y12_t, Y21_t, Y22_t = PNM.ybus_branch_entries(tap, PNM.NetworkReductionData()) Y11_m, Y12_m, Y21_m, Y22_m = PNM.ybus_branch_entries(mbp, PNM.NetworkReductionData()) @test Y11_m ≈ Y11_l + Y11_t diff --git a/test/test_data/modified_14bus_system_no_icd.txt b/test/test_data/modified_14bus_system_no_icd.txt new file mode 100644 index 00000000..84de40b1 --- /dev/null +++ b/test/test_data/modified_14bus_system_no_icd.txt @@ -0,0 +1,73 @@ + 107, 107, 3.10060620307922 , -10019.4659313190 + 107, 1100002, -0.00000000000000 , 10000.0000000000 + 107, 108, -3.10060620307922 , 19.4662113189697 + 1100002, 1100002, 0.00000000000000 , -30000.0000000000 + 1100002, 107, -0.00000000000000 , 10000.0000000000 + 1100002, 104, -0.00000000000000 , 10000.0000000000 + 1100002, 109, -0.00000000000000 , 10000.0000000000 + 103, 103, 132.904750823975 , -168.593001680914 + 103, 102, -19.0339088439941 , 120.448959350586 + 103, 104, -113.870841979980 , 48.1467323303223 + 102, 102, 5412.94931411743 , -4402.66704886226 + 102, 103, -19.0339088439941 , 120.448959350586 + 102, 104, -515.366577148438 , 376.782653808594 + 102, 101, -4878.04882812500 , 3902.43872070312 + 104, 104, 1028.59844207764 , -30440.8992626672 + 104, 1100002, -0.00000000000000 , 10000.0000000000 + 104, 103, -113.870841979980 , 48.1467323303223 + 104, 102, -515.366577148438 , 376.782653808594 + 104, 109, -0.00000000000000 , 10000.0000000000 + 104, 106, -0.00000000000000 , 10000.0000000000 + 104, 101, -399.361022949219 , 15.9744415283203 + 601, 601, 255.515943527222 , -266.450859436765 + 601, 301, -18.8712596893311 , 2.16857123374939 + 601, 109, -236.644683837891 , 264.337158203125 + 501, 501, 170.927175521851 , -68.8906958228908 + 501, 201, -157.686218261719 , 68.1742706298828 + 501, 109, -13.2409572601318 , 0.720955193042755 + 401, 401, 174.032964706421 , -239.601025610114 + 401, 200, -150.343994140625 , 102.884414672852 + 401, 109, -23.6889705657959 , 136.741210937500 + 301, 301, 24.4906101226807 , -51.2181056383997 + 301, 601, -18.8712596893311 , 2.16857123374939 + 301, 114, -5.61935043334961 , 49.1082344055176 + 201, 201, 253.106613159180 , -143.700978166540 + 201, 501, -157.686218261719 , 68.1742706298828 + 201, 114, -95.4203948974609 , 75.5306625366211 + 200, 200, 951.584594726562 , -167.455695092678 + 200, 401, -150.343994140625 , 102.884414672852 + 200, 114, -801.240600585938 , 64.6161804199219 + 114, 114, 902.280345916748 , -10189.1999023607 + 114, 301, -5.61935043334961 , 49.1082344055176 + 114, 201, -95.4203948974609 , 75.5306625366211 + 114, 200, -801.240600585938 , 64.6161804199219 + 114, 1100001, -0.00000000000000 , 10000.0000000000 + 112, 112, 76.9152069091797 , -10216.5844952832 + 112, 106, -76.9152069091797 , 216.588165283203 + 112, 1100001, -0.00000000000000 , 10000.0000000000 + 109, 109, 273.574611663818 , -30401.7677043328 + 109, 1100002, -0.00000000000000 , 10000.0000000000 + 109, 104, -0.00000000000000 , 10000.0000000000 + 109, 601, -236.644683837891 , 264.337158203125 + 109, 501, -13.2409572601318 , 0.720955193042755 + 109, 401, -23.6889705657959 , 136.741210937500 + 109, 110, -0.00000000000000 , 10000.0000000000 + 111, 111, 802.165511131287 , -149.270100883659 + 111, 110, -5.99848079681396 , 40.6076927185059 + 111, 106, -795.167030334473 , 110.667713165283 + 110, 110, 6.19848079979420 , -20036.6051827185 + 110, 109, -0.00000000000000 , 10000.0000000000 + 110, 111, -5.99848079681396 , 40.6076927185059 + 110, 1100001, -0.00000000000000 , 10000.0000000000 + 106, 106, 872.882237255573 , -10326.9494134365 + 106, 104, -0.00000000000000 , 10000.0000000000 + 106, 112, -76.9152069091797 , 216.588165283203 + 106, 111, -795.167030334473 , 110.667713165283 + 1100001, 1100001, 0.00000000000000 , -30000.0000000000 + 1100001, 114, -0.00000000000000 , 10000.0000000000 + 1100001, 112, -0.00000000000000 , 10000.0000000000 + 1100001, 110, -0.00000000000000 , 10000.0000000000 + + ZERO IMPEDANCE LINE CONNECTED BUSES (BUS IN MATRIX LISTED FIRST): + 104, 105 + 112, 113 diff --git a/test/test_data/modified_14bus_system_off_nominal_icd_ymatrix.txt b/test/test_data/modified_14bus_system_off_nominal_icd_ymatrix.txt new file mode 100644 index 00000000..60337383 --- /dev/null +++ b/test/test_data/modified_14bus_system_off_nominal_icd_ymatrix.txt @@ -0,0 +1,73 @@ + 107, 107, 3.10060620307922 , -10019.4659313190 + 107, 1100002, -0.00000000000000 , 10000.0000000000 + 107, 108, -3.10060620307922 , 19.4662113189697 + 1100002, 1100002, 0.00000000000000 , -30000.0000000000 + 1100002, 107, -0.00000000000000 , 10000.0000000000 + 1100002, 104, -0.00000000000000 , 10000.0000000000 + 1100002, 109, -0.00000000000000 , 10000.0000000000 + 103, 103, 132.904750823975 , -168.593001680914 + 103, 102, -19.0339088439941 , 120.448959350586 + 103, 104, -113.870841979980 , 48.1467323303223 + 102, 102, 5412.94931411743 , -4402.66704886226 + 102, 103, -19.0339088439941 , 120.448959350586 + 102, 104, -515.366577148438 , 376.782653808594 + 102, 101, -4878.04882812500 , 3902.43872070312 + 104, 104, 1028.59844207764 , -30149.6375439172 + 104, 1100002, -0.00000000000000 , 10000.0000000000 + 104, 103, -113.870841979980 , 48.1467323303223 + 104, 102, -515.366577148438 , 376.782653808594 + 104, 109, 0.00000000000000 , 10219.7246093750 + 104, 106, -0.00000000000000 , 10000.0000000000 + 104, 101, -399.361022949219 , 15.9744415283203 + 601, 601, 255.515943527222 , -266.450859436765 + 601, 301, -18.8712596893311 , 2.16857123374939 + 601, 109, -236.644683837891 , 264.337158203125 + 501, 501, 170.927175521851 , -68.8906958228908 + 501, 201, -157.686218261719 , 68.1742706298828 + 501, 109, -13.2409572601318 , 0.720955193042755 + 401, 401, 174.032964706421 , -239.601025610114 + 401, 200, -150.343994140625 , 102.884414672852 + 401, 109, -23.6889705657959 , 136.741210937500 + 301, 301, 24.4906101226807 , -51.2181056383997 + 301, 601, -18.8712596893311 , 2.16857123374939 + 301, 114, -5.61935043334961 , 49.1082344055176 + 201, 201, 253.106613159180 , -143.700978166540 + 201, 501, -157.686218261719 , 68.1742706298828 + 201, 114, -95.4203948974609 , 75.5306625366211 + 200, 200, 951.584594726562 , -167.455695092678 + 200, 401, -150.343994140625 , 102.884414672852 + 200, 114, -801.240600585938 , 64.6161804199219 + 114, 114, 902.280345916748 , -10189.1999023607 + 114, 301, -5.61935043334961 , 49.1082344055176 + 114, 201, -95.4203948974609 , 75.5306625366211 + 114, 200, -801.240600585938 , 64.6161804199219 + 114, 1100001, -0.00000000000000 , 10000.0000000000 + 112, 112, 76.9152069091797 , -9307.49367497070 + 112, 106, -76.9152069091797 , 216.588165283203 + 112, 1100001, -1578.61914062500 , 8952.79785156250 + 109, 109, 273.574611663818 , -30401.5167277703 + 109, 1100002, -0.00000000000000 , 10000.0000000000 + 109, 104, 0.00000000000000 , 10219.7246093750 + 109, 601, -236.644683837891 , 264.337158203125 + 109, 501, -13.2409572601318 , 0.720955193042755 + 109, 401, -23.6889705657959 , 136.741210937500 + 109, 110, 1953.09338378906 , 9033.41796875000 + 111, 111, 802.165511131287 , -149.270100883659 + 111, 110, -5.99848079681396 , 40.6076927185059 + 111, 106, -795.167030334473 , 110.667713165283 + 110, 110, 6.19848079979420 , -19278.7487374060 + 110, 109, -1953.09338378906 , 9033.41796875000 + 110, 111, -5.99848079681396 , 40.6076927185059 + 110, 1100001, -0.00000000000000 , 10000.0000000000 + 106, 106, 872.882237255573 , -10326.9494134365 + 106, 104, -0.00000000000000 , 10000.0000000000 + 106, 112, -76.9152069091797 , 216.588165283203 + 106, 111, -795.167030334473 , 110.667713165283 + 1100001, 1100001, 0.00000000000000 , -29090.9091796875 + 1100001, 114, -0.00000000000000 , 10000.0000000000 + 1100001, 112, 1578.61914062500 , 8952.79785156250 + 1100001, 110, -0.00000000000000 , 10000.0000000000 + + ZERO IMPEDANCE LINE CONNECTED BUSES (BUS IN MATRIX LISTED FIRST): + 104, 105 + 112, 113 diff --git a/test/test_ybus_contingencies.jl b/test/test_ybus_contingencies.jl index a41f3e36..c75f8132 100644 --- a/test/test_ybus_contingencies.jl +++ b/test/test_ybus_contingencies.jl @@ -230,7 +230,7 @@ end arc_tuple = PNM.get_arc_axis(nr)[arc_mod.arc_index] if haskey(nr.direct_branch_map, arc_tuple) br = nr.direct_branch_map[arc_tuple] - _, Y12, _, _ = PNM.ybus_branch_entries(br) + _, Y12, _, _ = PNM.ybus_branch_entries(br, nr) @test isapprox(arc_mod.delta_y12, PNM.YBUS_ELTYPE(-Y12); atol = 1e-6) end end diff --git a/test/test_ybus_psse.jl b/test/test_ybus_psse.jl index d9f1138a..960ce1b5 100644 --- a/test/test_ybus_psse.jl +++ b/test/test_ybus_psse.jl @@ -85,7 +85,7 @@ end if col_bus ∈ keys(nr.reverse_bus_search_map) col_bus = nr.reverse_bus_search_map[col_bus] end - @test isapprox(Ybus_pnm[row_bus, col_bus], val, atol = 1e-3) + @test isapprox(Ybus_pnm[row_bus, col_bus], val, rtol = 2 * eps(Float32), atol = 0.0) end end @@ -118,7 +118,7 @@ end if col_bus ∈ keys(nr.reverse_bus_search_map) col_bus = nr.reverse_bus_search_map[col_bus] end - @test isapprox(Ybus_pnm[row_bus, col_bus], val, atol = 2e-3) + @test isapprox(Ybus_pnm[row_bus, col_bus], val, rtol = 4 * eps(Float32), atol = 0.0) end end @@ -202,7 +202,7 @@ end if col_bus ∈ keys(nr.reverse_bus_search_map) col_bus = nr.reverse_bus_search_map[col_bus] end - @test isapprox(Ybus_pnm[row_bus, col_bus], val; atol = 1e-1) + @test isapprox(Ybus_pnm[row_bus, col_bus], val; rtol = 9 * eps(Float32), atol = 0.0) end end @@ -238,6 +238,77 @@ end if col_bus ∈ keys(nr.reverse_bus_search_map) col_bus = nr.reverse_bus_search_map[col_bus] end - @test isapprox(Ybus_pnm[row_bus, col_bus], val; atol = 1e-5) + @test isapprox(Ybus_pnm[row_bus, col_bus], val; rtol = 2 * eps(Float32), atol = 0.0) + end +end + +@testset "psse_modified_14bus_off_nominal_icd" begin + # Three transformers have off-nominal settings that activate their ICTs: + # BUS 109-BUS 104-i_1 (TapTransformer): WINDV1=0.95, table 4 → F=1.030 + # BUS 110-BUS 109-i_1 (PhaseShiftingTransformer): ANG1=12.2°, table 3 → F=1.082 + # BUS 113-BUS 110-BUS 114-i_1 winding 1 (PST3W): ANG1=10.0°, table 9 → F=1.100 + sys = build_system(PSSEParsingTestSystems, "psse_modified_14bus_off_nominal_icd") + Ybus_pnm = Ybus(sys) + nr = Ybus_pnm.network_reduction_data + + Ybus_psse, b_ix_psse, row_buses, col_buses, y_values, reduced_bus_pairs_psse = + parse_psse_ybus( + joinpath(TEST_DATA_DIR, "modified_14bus_system_off_nominal_icd_ymatrix.txt"), + ) + for x in reduced_bus_pairs_psse + _test_psse_reduction_row(x, nr.reverse_bus_search_map) + end + + psse_ref_bus_numbers = [101, 108] + psse_skip_indices = indexin(psse_ref_bus_numbers, Ybus_pnm.axes[1]) + n_psse_excluded_elements = nnz(Ybus_pnm.data[psse_skip_indices, :]) + @test nnz(Ybus_pnm.data) == length(filter(!iszero, y_values)) + n_psse_excluded_elements + + psse_to_pnm_star_bus = Dict(1100001 => 1002, 1100002 => 1001) + for (row_bus, col_bus, val) in zip(row_buses, col_buses, y_values) + row_bus = get(psse_to_pnm_star_bus, row_bus, row_bus) + col_bus = get(psse_to_pnm_star_bus, col_bus, col_bus) + if row_bus ∈ keys(nr.reverse_bus_search_map) + row_bus = nr.reverse_bus_search_map[row_bus] + end + if col_bus ∈ keys(nr.reverse_bus_search_map) + col_bus = nr.reverse_bus_search_map[col_bus] + end + @test isapprox(Ybus_pnm[row_bus, col_bus], val; rtol = 2 * eps(Float32), atol = 0.0) + end +end + +@testset "psse_modified_14bus_nominal_icd" begin + sys = build_system(PSSEParsingTestSystems, "psse_modified_14bus_nominal_icd") + Ybus_pnm = Ybus(sys) + nr = Ybus_pnm.network_reduction_data + + Ybus_psse, b_ix_psse, row_buses, col_buses, y_values, reduced_bus_pairs_psse = + parse_psse_ybus(joinpath(TEST_DATA_DIR, "modified_14bus_system_no_icd.txt")) + for x in reduced_bus_pairs_psse + _test_psse_reduction_row(x, nr.reverse_bus_search_map) + end + + # PSS/E excluded rows for both type-3 buses it saw (101 and 108); count those + # PNM entries to reconcile the NNZ check. + psse_ref_bus_numbers = [101, 108] + psse_skip_indices = indexin(psse_ref_bus_numbers, Ybus_pnm.axes[1]) + n_psse_excluded_elements = nnz(Ybus_pnm.data[psse_skip_indices, :]) + @test nnz(Ybus_pnm.data) == length(filter(!iszero, y_values)) + n_psse_excluded_elements + + # PSS/E numbers its switch star-buses as 1100001/1100002; PNM assigns 1002/1001. + # Adjacency-derived mapping: PSS/E 1100001 ↔ {114,112,110} = PNM 1002; + # PSS/E 1100002 ↔ {107,104,109} = PNM 1001. + psse_to_pnm_star_bus = Dict(1100001 => 1002, 1100002 => 1001) + for (row_bus, col_bus, val) in zip(row_buses, col_buses, y_values) + row_bus = get(psse_to_pnm_star_bus, row_bus, row_bus) + col_bus = get(psse_to_pnm_star_bus, col_bus, col_bus) + if row_bus ∈ keys(nr.reverse_bus_search_map) + row_bus = nr.reverse_bus_search_map[row_bus] + end + if col_bus ∈ keys(nr.reverse_bus_search_map) + col_bus = nr.reverse_bus_search_map[col_bus] + end + @test isapprox(Ybus_pnm[row_bus, col_bus], val; rtol = 2 * eps(Float32), atol = 0.0) end end diff --git a/test/test_ybus_reductions.jl b/test/test_ybus_reductions.jl index fb4c093f..2cdf5328 100644 --- a/test/test_ybus_reductions.jl +++ b/test/test_ybus_reductions.jl @@ -774,10 +774,81 @@ end blind11 = zero(eltype(ybus.data)) blind12 = zero(eltype(ybus.data)) for br in bp.branches - entries = PNM.ybus_branch_entries(br) + entries = PNM.ybus_branch_entries(br, nr) blind11 += entries[1] blind12 += entries[2] end @test !isapprox(blind11, ybus.data[ip, ip]) @test !isapprox(blind12, ybus.data[ip, iq]) end + +@testset "GenericArcImpedance in series (degree-two) reduction" begin + sys = PSB.build_system(PSITestSystems, "c_sys5") + template = first(get_components(ACBus, sys)) + grid_bus = first( + b for b in get_components(ACBus, sys) if + get_bustype(b) != PSY.ACBusTypes.REF && get_bustype(b) != PSY.ACBusTypes.ISOLATED + ) + grid_bus_num = get_number(grid_bus) + function _mk_bus(num, name) + b = deepcopy(template) + b.internal = IS.InfrastructureSystemsInternal() + set_number!(b, num) + set_name!(b, name) + set_bustype!(b, PSY.ACBusTypes.PQ) + return b + end + junction = _mk_bus(901, "GAI_JUNCTION") + stub = _mk_bus(902, "GAI_STUB") + foreach(b -> add_component!(sys, b), (junction, stub)) + + # grid_bus -- junction (Line) + arc1 = Arc(grid_bus, junction) + add_component!(sys, arc1) + add_component!( + sys, + Line( + "GAI_grid_to_junction", + true, + 0.0, + 0.0, + arc1, + 0.01, + 0.10, + (from = 0.0, to = 0.0), + 100.0, + (-1.5, 1.5), + ), + ) + + # junction -- stub (GenericArcImpedance), forming the second segment of the chain. + arc2 = Arc(junction, stub) + add_component!(sys, arc2) + add_component!( + sys, + PSY.GenericArcImpedance(; + name = "GAI_junction_to_stub", + available = true, + active_power_flow = 0.0, + reactive_power_flow = 0.0, + max_flow = 100.0, + arc = arc2, + r = 0.02, + x = 0.08, + ), + ) + + ybus_full = Ybus(sys) + ybus = Ybus(sys; network_reductions = NetworkReduction[DegreeTwoReduction()]) + nrd = PNM.get_network_reduction_data(ybus) + + @test 901 ∉ PNM.get_bus_axis(ybus) + @test 902 ∈ PNM.get_bus_axis(ybus) + @test PSY.GenericArcImpedance in PNM.get_ac_transmission_types(nrd) + + @test isapprox( + ybus[grid_bus_num, 902]^-1, + ybus_full[grid_bus_num, 901]^-1 + ybus_full[901, 902]^-1; + rtol = sqrt(eps(real(YBUS_ELTYPE))), + ) +end