diff --git a/Project.toml b/Project.toml
index 97404bc9..48f1e3f7 100644
--- a/Project.toml
+++ b/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"
diff --git a/src/BA_ABA_matrices.jl b/src/BA_ABA_matrices.jl
index 7902317a..2b53cf3d 100644
--- a/src/BA_ABA_matrices.jl
+++ b/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
 
 """
diff --git a/src/BranchesParallel.jl b/src/BranchesParallel.jl
index cb66e831..909f0e50 100644
--- a/src/BranchesParallel.jl
+++ b/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
diff --git a/src/BranchesSeries.jl b/src/BranchesSeries.jl
index 0b38a7ae..1fa87c78 100644
--- a/src/BranchesSeries.jl
+++ b/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
 
 """
diff --git a/src/ThreeWindingTransformerWinding.jl b/src/ThreeWindingTransformerWinding.jl
index c4a9e290..c6b841ce 100644
--- a/src/ThreeWindingTransformerWinding.jl
+++ b/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
 
diff --git a/src/Ybus.jl b/src/Ybus.jl
index 7ac38aaf..b23ae42e 100644
--- a/src/Ybus.jl
+++ b/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)",
diff --git a/src/common.jl b/src/common.jl
index 660d0779..16c48eac 100644
--- a/src/common.jl
+++ b/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
diff --git a/src/connectivity_checks.jl b/src/connectivity_checks.jl
index 3ffe5e3d..7a8ce254 100644
--- a/src/connectivity_checks.jl
+++ b/src/connectivity_checks.jl
@@ -83,24 +83,30 @@ function find_connected_components(sys::PSY.System)
     return find_connected_components(a.data, a.lookup[1])
 end
 
-# this function extends the PowerModels.jl implementation to accept an adjacency matrix and bus lookup
-function find_connected_components(M, bus_lookup::Dict{Int64, Int64})
-    pm_buses = Dict([i => Dict("bus_type" => 1, "bus_i" => b) for (i, b) in bus_lookup])
-
-    arcs = findall((LinearAlgebra.UpperTriangular(M) - LinearAlgebra.I) .!= 0)
-    pm_branches = Dict([
-        i => Dict("f_bus" => a[1], "t_bus" => a[2], "br_status" => 1) for
-        (i, a) in enumerate(arcs)
-    ],)
-
-    data = Dict("bus" => pm_buses, "branch" => pm_branches)
-    cc = PSY.calc_connected_components(data)
-    bus_decode = Dict(value => key for (key, value) in bus_lookup)
-    connected_components = Vector{Set{Int64}}()
-    for c in cc
-        push!(connected_components, Set([bus_decode[b] for b in c]))
+# Group bus numbers into connected components of the graph whose edges are the
+# off-diagonal nonzeros of `M`. This replaces the former PowerModels.jl
+# `calc_connected_components` round-trip, which PowerSystems no longer re-exports.
+function find_connected_components(
+    M::SparseArrays.SparseMatrixCSC,
+    bus_lookup::Dict{Int64, Int64},
+)
+    n = length(bus_lookup)
+    bus_decode = Dict(index => bus_number for (bus_number, index) in bus_lookup)
+    uf = collect(1:n)
+    rows = SparseArrays.rowvals(M)
+    for col in 1:n
+        for k in SparseArrays.nzrange(M, col)
+            row = rows[k]
+            # Any off-diagonal nonzero is an undirected edge; ignore self-loops.
+            row != col && union_sets!(uf, col, row)
+        end
+    end
+    components = Dict{Int, Set{Int64}}()
+    for index in 1:n
+        root = get_representative(uf, index)
+        push!(get!(() -> Set{Int64}(), components, root), bus_decode[index])
     end
-    return Set(connected_components)
+    return Set(values(components))
 end
 
 function dfs_connectivity(M, ::Vector{PSY.ACBus}, bus_lookup::Dict{Int64, Int64})
diff --git a/src/network_modification.jl b/src/network_modification.jl
index 08b94b07..b11e3309 100644
--- a/src/network_modification.jl
+++ b/src/network_modification.jl
@@ -9,7 +9,7 @@ function _compute_parallel_partial_ybus_delta(
     delta_b::Float64,
 )::NTuple{4, YBUS_ELTYPE}
     for br in bp.branches
-        b_circuit = get_series_susceptance(br)
+        b_circuit = get_series_susceptance(br, PSY.SU)
         if isapprox(-b_circuit, delta_b; atol = YBUS_DELTA_TOL, rtol = 0)
             Y11, Y12, Y21, Y22 = ybus_branch_entries(br)
             return (YBUS_ELTYPE(-Y11), YBUS_ELTYPE(-Y12),
@@ -35,7 +35,7 @@ function _compute_arc_ybus_delta(
 )::NTuple{4, YBUS_ELTYPE}
     if haskey(nr.direct_branch_map, arc_tuple)
         br = nr.direct_branch_map[arc_tuple]
-        b_arc = get_series_susceptance(br)
+        b_arc = get_series_susceptance(br, PSY.SU)
         Y11, Y12, Y21, Y22 = ybus_branch_entries(br)
         if isapprox(delta_b, -b_arc; atol = YBUS_DELTA_TOL, rtol = 0)
             return (YBUS_ELTYPE(-Y11), YBUS_ELTYPE(-Y12),
@@ -47,7 +47,7 @@ function _compute_arc_ybus_delta(
         end
     elseif haskey(nr.parallel_branch_map, arc_tuple)
         bp = nr.parallel_branch_map[arc_tuple]
-        b_arc = get_series_susceptance(bp)
+        b_arc = get_series_susceptance(bp, PSY.SU)
         if isapprox(delta_b, -b_arc; atol = YBUS_DELTA_TOL, rtol = 0)
             Y11, Y12, Y21, Y22 = ybus_branch_entries(bp)
             return (YBUS_ELTYPE(-Y11), YBUS_ELTYPE(-Y12),
@@ -57,7 +57,7 @@ function _compute_arc_ybus_delta(
         end
     elseif haskey(nr.series_branch_map, arc_tuple)
         series_chain = nr.series_branch_map[arc_tuple]
-        b_arc = get_series_susceptance(series_chain)
+        b_arc = get_series_susceptance(series_chain, PSY.SU)
         if isapprox(delta_b, -b_arc; atol = YBUS_DELTA_TOL, rtol = 0)
             Y11, Y12, Y21, Y22 = ybus_branch_entries(series_chain)
             return (YBUS_ELTYPE(-Y11), YBUS_ELTYPE(-Y12),
@@ -70,7 +70,7 @@ function _compute_arc_ybus_delta(
         end
     elseif haskey(nr.transformer3W_map, arc_tuple)
         tr = nr.transformer3W_map[arc_tuple]
-        b_arc = get_series_susceptance(tr)
+        b_arc = get_series_susceptance(tr, PSY.SU)
         if !isapprox(delta_b, -b_arc; atol = YBUS_DELTA_TOL, rtol = 0)
             error(
                 "Partial Ybus delta is not supported on 3-winding transformer arcs. " *
@@ -245,7 +245,7 @@ function _classify_outage_component!(
         push!(direct_mods, ArcModification(arc_idx, -b_arc, dy11, dy12, dy21, dy22))
     elseif tag === :parallel
         arc_idx = arc_lookup[arc_tuple]
-        b_circuit = PSY.get_series_susceptance(component)
+        b_circuit = PSY.get_series_susceptance(component, PSY.SU)
         dy11, dy12, dy21, dy22 = _compute_arc_ybus_delta(nr, arc_tuple, -b_circuit)
         push!(parallel_mods, ArcModification(arc_idx, -b_circuit, dy11, dy12, dy21, dy22))
     elseif tag === :series
@@ -299,8 +299,8 @@ function _classify_outage_component!(
 )
     bus_ix = get_bus_index(component, bus_lookup, nr)
     Y =
-        PSY.get_impedance_active_power(component) -
-        im * PSY.get_impedance_reactive_power(component)
+        PSY.get_impedance_active_power(component, PSY.SU) -
+        im * PSY.get_impedance_reactive_power(component, PSY.SU)
     push!(shunt_mods, ShuntModification(bus_ix, YBUS_ELTYPE(-Y)))
     push!(component_names, PSY.get_name(component))
     return
@@ -375,7 +375,7 @@ function _classify_branch_modification(
         return [ArcModification(arc_idx, -b_arc, dy11, dy12, dy21, dy22)]
     elseif tag === :parallel
         arc_idx = arc_lookup[arc_tuple]
-        b_circuit = PSY.get_series_susceptance(branch)
+        b_circuit = PSY.get_series_susceptance(branch, PSY.SU)
         dy11, dy12, dy21, dy22 = _compute_arc_ybus_delta(nr, arc_tuple, -b_circuit)
         return [ArcModification(arc_idx, -b_circuit, dy11, dy12, dy21, dy22)]
     elseif tag === :series
diff --git a/src/virtual_factor_helpers.jl b/src/virtual_factor_helpers.jl
index 3188845f..bdc941e3 100644
--- a/src/virtual_factor_helpers.jl
+++ b/src/virtual_factor_helpers.jl
@@ -97,12 +97,12 @@ function _extract_branch_susceptances_by_arc(
         if haskey(nr_data.parallel_branch_map, arc)
             bp = nr_data.parallel_branch_map[arc]
             result[j] = Float64[
-                get_series_susceptance(branch) for branch in bp.branches
+                get_series_susceptance(branch, PSY.SU) for branch in bp.branches
             ]
         elseif haskey(nr_data.series_branch_map, arc)
             bs = nr_data.series_branch_map[arc]
             result[j] = Float64[
-                get_series_susceptance(segment) for segment in bs
+                get_series_susceptance(segment, PSY.SU) for segment in bs
             ]
         else
             result[j] = [arc_b]
diff --git a/test/Project.toml b/test/Project.toml
index 8196a905..ecdeb90c 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+AppleAccelerate = "13e28ba4-7ad8-5781-acae-3021b1ed3924"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
@@ -9,6 +10,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+PowerFlowFileParser = "bed98974-b02e-5e2f-9ee0-a103f5c450dd"
 PowerNetworkMatrices = "bed98974-b02a-5e2f-9fe0-a103f5c450dd"
 PowerSystemCaseBuilder = "f00506e0-b84f-492a-93c2-c0a9afc4364e"
 PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
@@ -18,5 +20,11 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TimeSeries = "9e3dc215-6440-5c97-bce1-76c03772f85e"
 
+[sources]
+InfrastructureSystems = {url = "https://github.com/NREL-Sienna/InfrastructureSystems.jl", rev = "IS4"}
+PowerFlowFileParser = {url = "https://github.com/NREL-Sienna/PowerFlowFileParser.jl", rev = "main"}
+PowerSystemCaseBuilder = {url = "https://github.com/NREL-Sienna/PowerSystemCaseBuilder.jl", rev = "psy6"}
+PowerSystems = {url = "https://github.com/NREL-Sienna/PowerSystems.jl", rev = "psy6"}
+
 [compat]
 julia = "^1.10"
diff --git a/test/test_A_BA_ABA_with_radial_lines.jl b/test/test_A_BA_ABA_with_radial_lines.jl
index 3ec32646..39d4c575 100644
--- a/test/test_A_BA_ABA_with_radial_lines.jl
+++ b/test/test_A_BA_ABA_with_radial_lines.jl
@@ -78,7 +78,7 @@ end
             !PSY.get_available(source) && continue
             bus = PSY.get_bus(source)
             bus_ix = bus_lookup[PSY.get_number(bus)]
-            bus_activepower_injection[bus_ix] += PSY.get_active_power(source)
+            bus_activepower_injection[bus_ix] += PSY.get_active_power(source, PSY.SU)
         end
         bus_activepower_withdrawals = zeros(Float64, n_buses)
         loads = PSY.get_components(x -> !isa(x, PSY.FixedAdmittance), PSY.ElectricLoad, sys)
@@ -86,7 +86,7 @@ end
             !PSY.get_available(l) && continue
             bus = PSY.get_bus(l)
             bus_ix = bus_lookup[PSY.get_number(bus)]
-            bus_activepower_withdrawals[bus_ix] += PSY.get_active_power(l)
+            bus_activepower_withdrawals[bus_ix] += PSY.get_active_power(l, PSY.SU)
         end
         power_injection =
             deepcopy(bus_activepower_injection - bus_activepower_withdrawals)
diff --git a/test/test_arc_types_and_reductions.jl b/test/test_arc_types_and_reductions.jl
index c734d29a..059d71f1 100644
--- a/test/test_arc_types_and_reductions.jl
+++ b/test/test_arc_types_and_reductions.jl
@@ -134,6 +134,14 @@ end
         winding_group_number = WindingGroupNumber.GROUP_11,
     )
 
+    # Attach the branches to a system so the unit-aware getters used by
+    # ybus_branch_entries can resolve the system base.
+    sys = System(100.0)
+    add_component!(sys, bus1)
+    add_component!(sys, bus2)
+    add_component!(sys, line)
+    add_component!(sys, tap)
+
     # Homogeneous group dispatches to BranchesParallel{Line}.
     bp_homog = PNM.BranchesParallel([line])
     @test bp_homog isa PNM.BranchesParallel{PSY.Line}
diff --git a/test/test_connectivity.jl b/test/test_connectivity.jl
index 0fe8f01c..92354fdd 100644
--- a/test/test_connectivity.jl
+++ b/test/test_connectivity.jl
@@ -127,17 +127,26 @@ end
         irreducible_buses = Set(collect(1:14)),
     )
 
+    # The three constructions agree to floating-point noise. Compare with an
+    # absolute tolerance: many entries are physically zero, so a relative
+    # tolerance (the `≈`/isapprox default) is meaningless there, and the dense
+    # and on-demand paths drift by a few ULPs run-to-run on the IS4/psy6 stack.
+    atol = 1e-10
     for i in ptdf_1.axes[1], j in ptdf_1.axes[2]
-        @test ptdf_1[j, i] == ptdf_2[j, i] == ptdf_3[j, i]
+        @test isapprox(ptdf_1[j, i], ptdf_2[j, i]; atol = atol) &&
+              isapprox(ptdf_2[j, i], ptdf_3[j, i]; atol = atol)
     end
     for i in lodf_1.axes[1], j in lodf_1.axes[2]
-        @test lodf_1[i, j] == lodf_2[i, j] == lodf_3[i, j]
+        @test isapprox(lodf_1[i, j], lodf_2[i, j]; atol = atol) &&
+              isapprox(lodf_2[i, j], lodf_3[i, j]; atol = atol)
     end
     for i in vptdf_1.axes[1], j in vptdf_1.axes[2]
-        @test vptdf_1[i, j] == vptdf_2[i, j] == vptdf_3[i, j]
+        @test isapprox(vptdf_1[i, j], vptdf_2[i, j]; atol = atol) &&
+              isapprox(vptdf_2[i, j], vptdf_3[i, j]; atol = atol)
     end
     for i in vlodf_1.axes[1], j in vlodf_1.axes[2]
-        @test vlodf_1[i, j] == vlodf_2[i, j] == vlodf_3[i, j]
+        @test isapprox(vlodf_1[i, j], vlodf_2[i, j]; atol = atol) &&
+              isapprox(vlodf_2[i, j], vlodf_3[i, j]; atol = atol)
     end
 end
 
diff --git a/test/test_equivalent_getters.jl b/test/test_equivalent_getters.jl
index 059a0686..710bc481 100644
--- a/test/test_equivalent_getters.jl
+++ b/test/test_equivalent_getters.jl
@@ -34,6 +34,12 @@
         rating = 150.0,  # rating
         angle_limits = (min = -π / 2, max = π / 2),
     )
+    # Attach the branches so the system-base getters (e.g. the susceptance
+    # weighting in get_impedance_averaged_rating) can resolve the system base.
+    # These lines carry round, illustrative values, so skip data validation.
+    PSY.add_component!(sys, line1; skip_validation = true)
+    PSY.add_component!(sys, line2; skip_validation = true)
+
     # Create BranchesParallel
     bp = PNM.BranchesParallel([line1, line2])
 
@@ -88,7 +94,7 @@ end
     rating3 = PNM.get_equivalent_rating(PNM.ThreeWindingTransformerWinding(trf, 3))
     # Should return winding-specific rating if non-zero, else transformer rating
     expected_rating3 =
-        trf.rating_tertiary == 0.0 ? PSY.get_rating(trf) : trf.rating_tertiary
+        trf.rating_tertiary == 0.0 ? PSY.get_rating(trf, PSY.DU) : trf.rating_tertiary
     @test rating3 == expected_rating3
 
     set_available_secondary!(trf, false)
diff --git a/test/test_lodf_with_radial_branches.jl b/test/test_lodf_with_radial_branches.jl
index 89da0db8..ac3cf8c8 100644
--- a/test/test_lodf_with_radial_branches.jl
+++ b/test/test_lodf_with_radial_branches.jl
@@ -57,7 +57,7 @@
             !PSY.get_available(source) && continue
             bus = PSY.get_bus(source)
             bus_ix = bus_lookup[PSY.get_number(bus)]
-            bus_activepower_injection[bus_ix] += PSY.get_active_power(source)
+            bus_activepower_injection[bus_ix] += PSY.get_active_power(source, PSY.SU)
         end
         bus_activepower_withdrawals = zeros(Float64, n_buses)
         loads = PSY.get_components(x -> !isa(x, PSY.FixedAdmittance), PSY.ElectricLoad, sys)
@@ -65,7 +65,7 @@
             !PSY.get_available(l) && continue
             bus = PSY.get_bus(l)
             bus_ix = bus_lookup[PSY.get_number(bus)]
-            bus_activepower_withdrawals[bus_ix] += PSY.get_active_power(l)
+            bus_activepower_withdrawals[bus_ix] += PSY.get_active_power(l, PSY.SU)
         end
         power_injection =
             deepcopy(bus_activepower_injection - bus_activepower_withdrawals)
diff --git a/test/test_modf_lodf_reductions.jl b/test/test_modf_lodf_reductions.jl
index 8a3caf8d..ed28e230 100644
--- a/test/test_modf_lodf_reductions.jl
+++ b/test/test_modf_lodf_reductions.jl
@@ -11,22 +11,14 @@ pair as a 2-circuit parallel group, enabling N-2 parallel-circuit tests.
 """
 function _clone_transformer_with_name(xfmr::PSY.TapTransformer, new_name::String)
     old_arc = PSY.get_arc(xfmr)
-    new_arc = PSY.Arc(; from = old_arc.from, to = old_arc.to)
-    return PSY.TapTransformer(;
-        name = new_name,
-        available = PSY.get_available(xfmr),
-        active_power_flow = PSY.get_active_power_flow(xfmr),
-        reactive_power_flow = PSY.get_reactive_power_flow(xfmr),
-        arc = new_arc,
-        r = PSY.get_r(xfmr),
-        x = PSY.get_x(xfmr),
-        primary_shunt = PSY.get_primary_shunt(xfmr),
-        tap = PSY.get_tap(xfmr),
-        rating = PSY.get_rating(xfmr),
-        base_power = PSY.get_base_power(xfmr),
-        base_voltage_primary = PSY.get_base_voltage_primary(xfmr),
-        base_voltage_secondary = PSY.get_base_voltage_secondary(xfmr),
-    )
+    # Deep-copy to reproduce every electrical parameter exactly (avoids re-reading
+    # each field in a particular unit basis), then give it a fresh identity, name,
+    # and an arc that shares the original buses rather than the deep-copied ones.
+    clone = deepcopy(xfmr)
+    clone.internal = IS.InfrastructureSystemsInternal()
+    PSY.set_name!(clone, new_name)
+    PSY.set_arc!(clone, PSY.Arc(; from = old_arc.from, to = old_arc.to))
+    return clone
 end
 
 """
@@ -244,7 +236,7 @@ end
     @testset "parallel single-circuit" begin
         arc_tuple, arc_idx = _find_non_islanding_arc(vmodf, nrd.parallel_branch_map)
         parallel = nrd.parallel_branch_map[arc_tuple]
-        b_circuit = PSY.get_series_susceptance(first(parallel.branches))
+        b_circuit = PSY.get_series_susceptance(first(parallel.branches), PSY.SU)
         delta_b = -b_circuit
         verify_modf_lodf_identity(vmodf, vlodf, ptdf, arc_idx, delta_b)
     end
@@ -267,7 +259,7 @@ end
     @testset "parallel single-circuit" begin
         arc_tuple, arc_idx = _find_non_islanding_arc(vmodf, nrd.parallel_branch_map)
         parallel = nrd.parallel_branch_map[arc_tuple]
-        b_circuit = PSY.get_series_susceptance(first(parallel.branches))
+        b_circuit = PSY.get_series_susceptance(first(parallel.branches), PSY.SU)
         delta_b = -b_circuit
         verify_modf_lodf_identity(vmodf, vlodf, ptdf, arc_idx, delta_b)
     end
@@ -291,7 +283,7 @@ end
     @testset "parallel single-circuit" begin
         arc_tuple, arc_idx = _find_non_islanding_arc(vmodf, nrd.parallel_branch_map)
         parallel = nrd.parallel_branch_map[arc_tuple]
-        b_circuit = PSY.get_series_susceptance(first(parallel.branches))
+        b_circuit = PSY.get_series_susceptance(first(parallel.branches), PSY.SU)
         delta_b = -b_circuit
         verify_modf_lodf_identity(vmodf, vlodf, ptdf, arc_idx, delta_b)
     end
@@ -323,7 +315,7 @@ end
     @testset "parallel single-circuit" begin
         arc_tuple, arc_idx = _find_non_islanding_arc(vmodf, nrd.parallel_branch_map)
         parallel = nrd.parallel_branch_map[arc_tuple]
-        b_circuit = PSY.get_series_susceptance(first(parallel.branches))
+        b_circuit = PSY.get_series_susceptance(first(parallel.branches), PSY.SU)
         delta_b = -b_circuit
         verify_modf_lodf_identity(vmodf, vlodf, ptdf, arc_idx, delta_b)
     end
diff --git a/test/test_ptdf.jl b/test/test_ptdf.jl
index fb427491..c9d15373 100644
--- a/test/test_ptdf.jl
+++ b/test/test_ptdf.jl
@@ -198,7 +198,7 @@ end
             r = branch_2.r,
             x = branch_2.x,
             b = branch_2.b,
-            rating = get_rating(branch_2),
+            rating = get_rating(branch_2, PSY.DU),
             angle_limits = get_angle_limits(branch_2),
         ),
     )
diff --git a/test/test_ptdf_with_radial_lines.jl b/test/test_ptdf_with_radial_lines.jl
index 2957d9d0..b4b4dd0e 100644
--- a/test/test_ptdf_with_radial_lines.jl
+++ b/test/test_ptdf_with_radial_lines.jl
@@ -44,7 +44,7 @@
             !PSY.get_available(source) && continue
             bus = PSY.get_bus(source)
             bus_ix = bus_lookup[PSY.get_number(bus)]
-            bus_activepower_injection[bus_ix] += PSY.get_active_power(source)
+            bus_activepower_injection[bus_ix] += PSY.get_active_power(source, PSY.SU)
         end
         bus_activepower_withdrawals = zeros(Float64, n_buses)
         loads = PSY.get_components(x -> !isa(x, PSY.FixedAdmittance), PSY.ElectricLoad, sys)
@@ -52,7 +52,7 @@
             !PSY.get_available(l) && continue
             bus = PSY.get_bus(l)
             bus_ix = bus_lookup[PSY.get_number(bus)]
-            bus_activepower_withdrawals[bus_ix] += PSY.get_active_power(l)
+            bus_activepower_withdrawals[bus_ix] += PSY.get_active_power(l, PSY.SU)
         end
         power_injection =
             deepcopy(bus_activepower_injection - bus_activepower_withdrawals)
@@ -139,7 +139,7 @@ end
                 !PSY.get_available(source) && continue
                 bus = PSY.get_bus(source)
                 bus_ix = bus_lookup[PSY.get_number(bus)]
-                bus_activepower_injection[bus_ix] += PSY.get_active_power(source)
+                bus_activepower_injection[bus_ix] += PSY.get_active_power(source, PSY.SU)
             end
             bus_activepower_withdrawals = zeros(Float64, n_buses)
             loads =
@@ -148,7 +148,7 @@ end
                 !PSY.get_available(l) && continue
                 bus = PSY.get_bus(l)
                 bus_ix = bus_lookup[PSY.get_number(bus)]
-                bus_activepower_withdrawals[bus_ix] += PSY.get_active_power(l)
+                bus_activepower_withdrawals[bus_ix] += PSY.get_active_power(l, PSY.SU)
             end
             power_injection =
                 deepcopy(bus_activepower_injection - bus_activepower_withdrawals)
diff --git a/test/test_ward_reduction.jl b/test/test_ward_reduction.jl
index 2e512ae9..8b9210e2 100644
--- a/test/test_ward_reduction.jl
+++ b/test/test_ward_reduction.jl
@@ -131,9 +131,13 @@ end
     )
     existing_line_susceptance = PSY.get_series_susceptance(
         ptdf_2.network_reduction_data.direct_branch_map[(101, 102)],
+        PSY.SU,
     )
+    # The ward equivalent is a detached synthetic branch storing system-base
+    # impedance; read it back with device base (identity).
     ward_line_susceptance = PSY.get_series_susceptance(
         ptdf_2.network_reduction_data.added_arc_impedance_map[(101, 102)],
+        PSY.DU,
     )
     ward_multiplier =
         existing_line_susceptance / (existing_line_susceptance + ward_line_susceptance)
diff --git a/test/test_ybus_reductions.jl b/test/test_ybus_reductions.jl
index d9508004..9a8cd047 100644
--- a/test/test_ybus_reductions.jl
+++ b/test/test_ybus_reductions.jl
@@ -313,8 +313,8 @@ end
 end
 
 function _set_zero_impedance!(branch)
-    set_r!(branch, 0.0)
-    set_x!(branch, 1e-5)
+    set_r!(branch, 0.0 * PSY.SU)
+    set_x!(branch, 1e-5 * PSY.SU)
 end
 
 @testset "ZeroImpedanceBranchReduction: chained bus merge" begin
@@ -337,8 +337,8 @@ end
 @testset "ZeroImpedanceBranchReduction: transformer arcs are excluded" begin
     sys = PSB.build_system(PSB.PSITestSystems, "c_sys14")
     t = get_component(Transformer2W, sys, "Trans4")  # from=7, to=8
-    set_r!(t, 0.0)
-    set_x!(t, 1e-5)
+    set_r!(t, 0.0 * PSY.SU)
+    set_x!(t, 1e-5 * PSY.SU)
 
     ybus = Ybus(sys)
     nrd = ybus.network_reduction_data
@@ -399,8 +399,8 @@ end
     # makes it qualify and bus 3 merges into bus 2.
     sys = PSB.build_system(PSB.PSITestSystems, "c_sys14")
     line = get_component(Line, sys, "Line3")
-    set_r!(line, 0.0)
-    set_x!(line, 1e-3)  # susceptance ≈ 1 / 1e-3 = 1e3
+    set_r!(line, 0.0 * PSY.SU)
+    set_x!(line, 1e-3 * PSY.SU)  # susceptance ≈ 1 / 1e-3 = 1e3
 
     ybus_default = Ybus(sys)
     @test 3 ∈ PNM.get_bus_axis(ybus_default)
@@ -423,8 +423,8 @@ end
     # drops below the threshold, so the branch is retained instead.
     sys = PSB.build_system(PSB.PSITestSystems, "c_sys14")
     line = get_component(Line, sys, "Line3")  # bus 2 → 3
-    set_r!(line, 0.0)
-    set_x!(line, 0.0)
+    set_r!(line, 0.0 * PSY.SU)
+    set_x!(line, 0.0 * PSY.SU)
 
     ybus_default = Ybus(sys)
     @test 3 ∉ PNM.get_bus_axis(ybus_default)  # merged with the default tiny substitute reactance