Units migration: adopt unitful PowerSystems getter API

Project.toml

@@ -20,6 +20,10 @@ Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 [extensions]
 MKLPardisoExt = "Pardiso"
 
+[sources]
+InfrastructureSystems = {url = "https://github.com/NREL-Sienna/InfrastructureSystems.jl", rev = "IS4"}
+PowerSystems = {url = "https://github.com/NREL-Sienna/PowerSystems.jl", rev = "psy6"}
+
 [compat]
 DataStructures = "^0.19"
 DocStringExtensions = "~0.8, ~0.9"

src/BA_ABA_matrices.jl

@@ -110,7 +110,7 @@ function BA_Matrix(ybus::Ybus)
         # Get series susceptance from components, not the equivalent ybus for reductions of degree two nodes
         # This results in reduced error relative to the DC power flow result without reductions
         if is_arc_in_series_map(nr_data, arc)
-            b = get_series_susceptance(get_mapped_series_branch(nr_data, arc))
+            b = get_series_susceptance(get_mapped_series_branch(nr_data, arc), PSY.SU)
         else
             Yt = -1 * ybus.data[ix_from_bus, ix_to_bus]
             Zt = 1 / Yt
@@ -136,8 +136,8 @@ function BA_Matrix(ybus::Ybus)
     return BA_Matrix(data, axes, lookup, subnetwork_axes, ybus.network_reduction_data)
 end
 
-function get_series_susceptance(segment::PSY.ACTransmission)
-    return PSY.get_series_susceptance(segment)
+function get_series_susceptance(segment::PSY.ACTransmission, units::IS.AbstractUnitSystem)
+    return PSY.get_series_susceptance(segment, units)
 end
 
 """

src/BranchesParallel.jl

@@ -73,15 +73,18 @@ function compute_parallel_multiplier(
     b_branch = 0.0
     for br in parallel_branch_set
         if PSY.get_name(br) == branch_name
-            b_branch += PSY.get_series_susceptance(br)
+            b_branch += PSY.get_series_susceptance(br, PSY.SU)
         end
-        b_total += PSY.get_series_susceptance(br)
+        b_total += PSY.get_series_susceptance(br, PSY.SU)
     end
     return b_branch / b_total
 end
 
-function get_series_susceptance(segment::AbstractBranchesParallel)
-    return sum(get_series_susceptance(branch) for branch in segment.branches)
+function get_series_susceptance(
+    segment::AbstractBranchesParallel,
+    units::IS.AbstractUnitSystem,
+)
+    return sum(get_series_susceptance(branch, units) for branch in segment.branches)
 end
 
 function get_equivalent_physical_branch_parameters(bp::AbstractBranchesParallel)
@@ -129,7 +132,12 @@ physically splits across a parallel group. Throws `ArgumentError` if the total
 series susceptance is zero or non-finite.
 """
 function get_impedance_averaged_rating(bp::AbstractBranchesParallel)
-    b_total = sum(PSY.get_series_susceptance, bp.branches)
+    # The susceptance weights must share a consistent impedance base across the
+    # group, so use system base (SU) like the sibling `compute_parallel_multiplier`.
+    # Within a parallel group (a single bus pair) this equals the natural-units
+    # weighting; device base would mix bases when the branches differ in base power.
+    # Requires the branches to be attached to a system.
+    b_total = sum(PSY.get_series_susceptance(br, PSY.SU) for br in bp.branches)
     if !isfinite(b_total) || iszero(b_total)
         throw(
             ArgumentError(
@@ -138,7 +146,7 @@ function get_impedance_averaged_rating(bp::AbstractBranchesParallel)
         )
     end
     return sum(
-        PSY.get_series_susceptance(br) / b_total * get_equivalent_rating(br)
+        PSY.get_series_susceptance(br, PSY.SU) / b_total * get_equivalent_rating(br)
         for br in bp.branches
     )
 end

src/BranchesSeries.jl

@@ -104,8 +104,9 @@ Base.length(bs::BranchesSeries) =
 
 Base.eltype(::Type{BranchesSeries}) = PSY.ACTransmission
 
-function get_series_susceptance(series_chain::BranchesSeries)
-    series_susceptances_sum = sum(inv(get_series_susceptance(x)) for x in series_chain)
+function get_series_susceptance(series_chain::BranchesSeries, units::IS.AbstractUnitSystem)
+    series_susceptances_sum =
+        sum(inv(get_series_susceptance(x, units)) for x in series_chain)
     total_susceptance = 1 / series_susceptances_sum
     return total_susceptance
 end
@@ -145,7 +146,7 @@ _series_member_rating(branch::PSY.ACTransmission) = get_equivalent_rating(branch
 Return the rating for PSY.ACTransmission branches.
 """
 function get_equivalent_rating(bs::PSY.ACTransmission)
-    return PSY.get_rating(bs)
+    return PSY.get_rating(bs, PSY.DU)
 end
 
 """
@@ -154,7 +155,8 @@ end
 Rating is assumed to be max_flow for GenericArcImpedance.
 """
 function get_equivalent_rating(bs::PSY.GenericArcImpedance)
-    return PSY.get_max_flow(bs)
+    # Detached synthetic ward equivalent: read the stored value with device base.
+    return PSY.get_max_flow(bs, PSY.DU)
 end
 
 """
@@ -174,12 +176,12 @@ end
 Return the emergency rating for PSY.ACTransmission branches.
 """
 function get_equivalent_emergency_rating(branch::PSY.ACTransmission)
-    if isnothing(PSY.get_rating_b(branch))
+    if isnothing(PSY.get_rating_b(branch, PSY.DU))
         @debug "Branch $(get_name(branch)) has no 'rating_b' defined. Post-contingency limit is going to be set using normal-operation rating.
             \n Consider including post-contingency limits using set_rating_b!()."
-        return PSY.get_rating(branch)
+        return PSY.get_rating(branch, PSY.DU)
     end
-    return PSY.get_rating_b(branch)
+    return PSY.get_rating_b(branch, PSY.DU)
 end
 
 """
@@ -189,7 +191,7 @@ Return the emergency rating for PSY.GenericArcImpedance.
 """
 function get_equivalent_emergency_rating(branch::PSY.GenericArcImpedance)
     @debug "GenericArcImpedance $(get_name(branch)) has no emergency rating. Using max_flow as a proxy instead."
-    return PSY.get_max_flow(branch)
+    return PSY.get_max_flow(branch, PSY.DU)
 end
 
 """

src/ThreeWindingTransformerWinding.jl

@@ -31,10 +31,13 @@ function get_name(three_wt_winding::ThreeWindingTransformerWinding)
     return PSY.get_name(transformer) * "_winding_$winding"
 end
 
-function get_series_susceptance(segment::ThreeWindingTransformerWinding)
+function get_series_susceptance(
+    segment::ThreeWindingTransformerWinding,
+    units::IS.AbstractUnitSystem,
+)
     tfw = get_transformer(segment)
     winding_int = get_winding_number(segment)
-    return PSY.get_series_susceptances(tfw)[winding_int]
+    return PSY.get_series_susceptances(tfw, units)[winding_int]
 end
 
 """
@@ -48,11 +51,11 @@ function get_equivalent_r(tw::ThreeWindingTransformerWinding)
     winding_num = get_winding_number(tw)
 
     if winding_num == 1
-        return PSY.get_r_primary(tfw)
+        return PSY.get_r_primary(tfw, PSY.SU)
     elseif winding_num == 2
-        return PSY.get_r_secondary(tfw)
+        return PSY.get_r_secondary(tfw, PSY.SU)
     elseif winding_num == 3
-        return PSY.get_r_tertiary(tfw)
+        return PSY.get_r_tertiary(tfw, PSY.SU)
     else
         throw(ArgumentError("Invalid winding number: $winding_num"))
     end
@@ -69,11 +72,11 @@ function get_equivalent_x(tw::ThreeWindingTransformerWinding)
     winding_num = get_winding_number(tw)
 
     if winding_num == 1
-        return PSY.get_x_primary(tfw)
+        return PSY.get_x_primary(tfw, PSY.SU)
     elseif winding_num == 2
-        return PSY.get_x_secondary(tfw)
+        return PSY.get_x_secondary(tfw, PSY.SU)
     elseif winding_num == 3
-        return PSY.get_x_tertiary(tfw)
+        return PSY.get_x_tertiary(tfw, PSY.SU)
     else
         throw(ArgumentError("Invalid winding number: $winding_num"))
     end
@@ -92,7 +95,7 @@ function get_equivalent_b(tw::ThreeWindingTransformerWinding)
 
     if winding_num == 1
         # Only the primary winding has the shunt susceptance
-        return (from = PSY.get_b(tfw), to = 0.0)
+        return (from = PSY.get_b(tfw, PSY.SU), to = 0.0)
     elseif winding_num == 2 || winding_num == 3
         # Secondary and tertiary windings don't have shunt susceptance
         return (from = 0.0, to = 0.0)
@@ -112,20 +115,20 @@ function get_equivalent_rating(tw::ThreeWindingTransformerWinding)
     winding_num = get_winding_number(tw)
 
     winding_rating = if winding_num == 1
-        PSY.get_rating_primary(tfw)
+        PSY.get_rating_primary(tfw, PSY.DU)
     elseif winding_num == 2
-        PSY.get_rating_secondary(tfw)
+        PSY.get_rating_secondary(tfw, PSY.DU)
     elseif winding_num == 3
-        PSY.get_rating_tertiary(tfw)
+        PSY.get_rating_tertiary(tfw, PSY.DU)
     else
         throw(ArgumentError("Invalid winding number: $winding_num"))
     end
     if winding_rating != 0.0
         return winding_rating
-    elseif isnothing(PSY.get_rating(tfw))
+    elseif isnothing(PSY.get_rating(tfw, PSY.DU))
         return 0.0
     else
-        return PSY.get_rating(tfw)
+        return PSY.get_rating(tfw, PSY.DU)
     end
 end
 

src/Ybus.jl

@@ -423,7 +423,7 @@ function add_branch_entries_to_indexing_maps!(
 end
 
 _get_shunt(br::PSY.ACTransmission, node::Symbol) =
-    PSY.get_g(br)[node] + 1im * PSY.get_b(br)[node]
+    PSY.get_g(br, PSY.SU)[node] + 1im * PSY.get_b(br, PSY.SU)[node]
 _get_shunt(::PSY.DiscreteControlledACBranch, ::Symbol) = zero(YBUS_ELTYPE)
 
 """Ybus entries for a `Line` or `DiscreteControlledACBranch`. `min_x_eps` substitutes
@@ -432,8 +432,8 @@ function ybus_branch_entries(
     br::PSY.ACTransmission;
     min_x_eps::Float64 = ZERO_IMPEDANCE_X_EPSILON,
 )
-    r = PSY.get_r(br)
-    x = PSY.get_x(br)
+    r = PSY.get_r(br, PSY.SU)
+    x = PSY.get_x(br, PSY.SU)
     if r == 0.0 && x == 0.0
         @warn "Branch $(PSY.get_name(br)) has r=0.0 and x=0.0; substituting x=$(min_x_eps) to avoid division by zero. This branch will be reduced by ZeroImpedanceBranchReduction unless its endpoints are irreducible."
         x = min_x_eps
@@ -442,7 +442,7 @@ function ybus_branch_entries(
     Y11 = Y_l + _get_shunt(br, :from)
     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))",
+            "Data in $(PSY.get_name(br)) is incorrect. r = $(PSY.get_r(br, PSY.SU)), x = $(PSY.get_x(br, PSY.SU))",
         )
     end
     Y12 = -Y_l
@@ -455,11 +455,14 @@ function ybus_branch_entries(
     br::PSY.GenericArcImpedance;
     min_x_eps::Float64 = ZERO_IMPEDANCE_X_EPSILON,
 )
-    Y_l = (1 / (PSY.get_r(br) + PSY.get_x(br) * 1im))
+    # GenericArcImpedance is a detached ward equivalent whose r/x are already the
+    # system-base values; device base (DU) reads them back as identity (system base
+    # would need the system base power, which a detached component cannot resolve).
+    Y_l = (1 / (PSY.get_r(br, PSY.DU) + PSY.get_x(br, PSY.DU) * 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))",
+            "Data in $(PSY.get_name(br)) is incorrect. r = $(PSY.get_r(br, PSY.DU)), x = $(PSY.get_x(br, PSY.DU))",
         )
     end
     Y12 = -Y_l
@@ -502,10 +505,10 @@ function ybus_branch_entries(
     br::PSY.TwoWindingTransformer;
     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.SU) + PSY.get_x(br, PSY.SU) * 1im)
     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)
+    y_shunt = PSY.get_primary_shunt(br, PSY.SU)
     Y11 = Y_t / abs2(tap)
     if !isfinite(Y11) || !isfinite(Y_t) || !isfinite(y_shunt * c_tap)
         error(
@@ -526,39 +529,39 @@ function ybus_branch_entries(
     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.SU) + PSY.get_x_primary(br, PSY.SU) * 1im)
         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)
-        y_shunt = PSY.get_g(br) + im * PSY.get_b(br)
+        y_shunt = PSY.get_g(br, PSY.SU) + im * PSY.get_b(br, PSY.SU)
         if !isfinite(Y11) || !isfinite(y_shunt)
             error(
                 "Data in $(PSY.get_name(br)) is incorrect.
-                r_p = $(PSY.get_r_primary(br)), x_p = $(PSY.get_x_primary(br)), primary_turns_ratio = $(PSY.get_primary_turns_ratio(br))",
+                r_p = $(PSY.get_r_primary(br, PSY.SU)), x_p = $(PSY.get_x_primary(br, PSY.SU)), primary_turns_ratio = $(PSY.get_primary_turns_ratio(br))",
             )
         end
         # 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.SU) + PSY.get_x_secondary(br, PSY.SU) * 1im)
         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)
         Y11 = Y_t / abs2(tap)
         if !isfinite(Y11)
             error(
                 "Data in $(PSY.get_name(br)) is incorrect.
-                r_s = $(PSY.get_r_secondary(br)), x_s = $(PSY.get_x_secondary(br)), secondary_turns_ratio = $(PSY.get_secondary_turns_ratio(br))",
+                r_s = $(PSY.get_r_secondary(br, PSY.SU)), x_s = $(PSY.get_x_secondary(br, PSY.SU)), secondary_turns_ratio = $(PSY.get_secondary_turns_ratio(br))",
             )
         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.SU) + PSY.get_x_tertiary(br, PSY.SU) * 1im)
         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)
         if !isfinite(Y11)
             error(
                 "Data in $(PSY.get_name(br)) is incorrect.
-                r_t = $(PSY.get_r_tertiary(br)), x_t = $(PSY.get_x_tertiary(br)), tertiary_turns_ratio = $(PSY.get_tertiary_turns_ratio(br))",
+                r_t = $(PSY.get_r_tertiary(br, PSY.SU)), x_t = $(PSY.get_x_tertiary(br, PSY.SU)), tertiary_turns_ratio = $(PSY.get_tertiary_turns_ratio(br))",
             )
         end
     end
@@ -722,7 +725,9 @@ function _ybus!(
     nr::NetworkReductionData,
 )
     bus_no = get_bus_index(fa, num_bus, nr)
-    Y = PSY.get_impedance_active_power(fa) - im * PSY.get_impedance_reactive_power(fa)
+    Y =
+        PSY.get_impedance_active_power(fa, PSY.SU) -
+        im * PSY.get_impedance_reactive_power(fa, PSY.SU)
     if !isfinite(Y)
         error(
             "Data in $(PSY.get_name(fa)) is incorrect. Y = $(Y)",

src/common.jl

@@ -490,7 +490,7 @@ function _segment_susceptance_after_outage(
     segment::PSY.ACTransmission,
     tripped_set::Set{<:PSY.ACTransmission},
 )::Float64
-    return segment ∈ tripped_set ? 0.0 : get_series_susceptance(segment)
+    return segment ∈ tripped_set ? 0.0 : get_series_susceptance(segment, PSY.SU)
 end
 
 function _segment_susceptance_after_outage(
@@ -500,7 +500,7 @@ function _segment_susceptance_after_outage(
     b_remaining = 0.0
     for branch in segment.branches
         if branch ∉ tripped_set
-            b_remaining += get_series_susceptance(branch)
+            b_remaining += get_series_susceptance(branch, PSY.SU)
         end
     end
     return b_remaining
@@ -539,7 +539,7 @@ function _compute_series_outage_delta_b(
     series_chain::BranchesSeries,
     tripped::Vector{<:PSY.ACTransmission},
 )::Float64
-    b_old = get_series_susceptance(series_chain)
+    b_old = get_series_susceptance(series_chain, PSY.SU)
     tripped_set = Set{PSY.ACTransmission}(tripped)
     remaining_inv_sum = 0.0
     for segment in series_chain