From 4694aec1767474d8d2254afda01ec12326e06184 Mon Sep 17 00:00:00 2001
From: Luke Kiernan <lkiernan@nrel.gov>
Date: Mon, 30 Mar 2026 15:50:17 -0600
Subject: [PATCH 1/5] Update PSY getter calls for unitful API: add Float64
 target

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
---
 src/BranchesParallel.jl                 |  2 +-
 src/BranchesSeries.jl                   |  8 +++----
 src/ThreeWindingTransformerWinding.jl   | 24 ++++++++++----------
 src/Ybus.jl                             | 30 ++++++++++++-------------
 test/test_A_BA_ABA_with_radial_lines.jl |  4 ++--
 test/test_lodf_with_radial_branches.jl  |  4 ++--
 test/test_ptdf.jl                       |  2 +-
 test/test_ptdf_with_radial_lines.jl     |  8 +++----
 8 files changed, 41 insertions(+), 41 deletions(-)

diff --git a/src/BranchesParallel.jl b/src/BranchesParallel.jl
index e086e6b03..fa6335658 100644
--- a/src/BranchesParallel.jl
+++ b/src/BranchesParallel.jl
@@ -85,7 +85,7 @@ This provides a conservative estimate that accounts for potential overestimation
 """
 function get_equivalent_rating(bp::BranchesParallel)
     # Sum of ratings divided by number of circuits
-    return sum(PSY.get_rating(branch) for branch in bp.branches) / length(bp.branches)
+    return sum(PSY.get_rating(branch, Float64) for branch in bp.branches) / length(bp.branches)
 end
 
 """
diff --git a/src/BranchesSeries.jl b/src/BranchesSeries.jl
index 84101d371..e91a79296 100644
--- a/src/BranchesSeries.jl
+++ b/src/BranchesSeries.jl
@@ -143,7 +143,7 @@ end
 Return the rating for PSY.ACTransmission branches.
 """
 function get_equivalent_rating(bs::PSY.ACTransmission)
-    return PSY.get_rating(bs)
+    return PSY.get_rating(bs, Float64)
 end
 
 """
@@ -163,12 +163,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, Float64))
         @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, Float64)
     end
-    return PSY.get_rating_b(branch)
+    return PSY.get_rating_b(branch, Float64)
 end
 
 """
diff --git a/src/ThreeWindingTransformerWinding.jl b/src/ThreeWindingTransformerWinding.jl
index f63723535..e3662bdbb 100644
--- a/src/ThreeWindingTransformerWinding.jl
+++ b/src/ThreeWindingTransformerWinding.jl
@@ -48,11 +48,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, Float64)
     elseif winding_num == 2
-        return PSY.get_r_secondary(tfw)
+        return PSY.get_r_secondary(tfw, Float64)
     elseif winding_num == 3
-        return PSY.get_r_tertiary(tfw)
+        return PSY.get_r_tertiary(tfw, Float64)
     else
         throw(ArgumentError("Invalid winding number: $winding_num"))
     end
@@ -69,11 +69,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, Float64)
     elseif winding_num == 2
-        return PSY.get_x_secondary(tfw)
+        return PSY.get_x_secondary(tfw, Float64)
     elseif winding_num == 3
-        return PSY.get_x_tertiary(tfw)
+        return PSY.get_x_tertiary(tfw, Float64)
     else
         throw(ArgumentError("Invalid winding number: $winding_num"))
     end
@@ -92,7 +92,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, Float64), 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 +112,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, Float64)
     elseif winding_num == 2
-        PSY.get_rating_secondary(tfw)
+        PSY.get_rating_secondary(tfw, Float64)
     elseif winding_num == 3
-        PSY.get_rating_tertiary(tfw)
+        PSY.get_rating_tertiary(tfw, Float64)
     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, Float64))
         return 0.0
     else
-        return PSY.get_rating(tfw)
+        return PSY.get_rating(tfw, Float64)
     end
 end
 
diff --git a/src/Ybus.jl b/src/Ybus.jl
index e68472f02..b7510633c 100644
--- a/src/Ybus.jl
+++ b/src/Ybus.jl
@@ -379,16 +379,16 @@ 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, Float64)[node] + 1im * PSY.get_b(br, Float64)[node]
 _get_shunt(::PSY.DiscreteControlledACBranch, ::Symbol) = zero(YBUS_ELTYPE)
 
 """Ybus entries for a `Line` or a `DiscreteControlledACBranch`."""
 function ybus_branch_entries(br::PSY.ACTransmission)
-    Y_l = (1 / (PSY.get_r(br) + PSY.get_x(br) * 1im))
+    Y_l = (1 / (PSY.get_r(br, Float64) + PSY.get_x(br, Float64) * 1im))
     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, Float64)), x = $(PSY.get_x(br, Float64))",
         )
     end
     Y12 = -Y_l
@@ -422,10 +422,10 @@ _get_tap(br::PSY.TwoWindingTransformer) = PSY.get_tap(br)
 
 """Ybus entries for a `Transformer2W`, `TapTransformer`, or `PhaseShiftingTransformer`."""
 function ybus_branch_entries(br::PSY.TwoWindingTransformer)
-    Y_t = 1 / (PSY.get_r(br) + PSY.get_x(br) * 1im)
+    Y_t = 1 / (PSY.get_r(br, Float64) + PSY.get_x(br, Float64) * 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, Float64)
     Y11 = Y_t / abs2(tap)
     if !isfinite(Y11) || !isfinite(Y_t) || !isfinite(y_shunt * c_tap)
         error(
@@ -443,39 +443,39 @@ function ybus_branch_entries(tp::ThreeWindingTransformerWinding)
     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, Float64) + PSY.get_x_primary(br, Float64) * 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, Float64) + im * PSY.get_b(br, Float64)
         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, Float64)), x_p = $(PSY.get_x_primary(br, Float64)), 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, Float64) + PSY.get_x_secondary(br, Float64) * 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, Float64)), x_s = $(PSY.get_x_secondary(br, Float64)), 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, Float64) + PSY.get_x_tertiary(br, Float64) * 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, Float64)), x_t = $(PSY.get_x_tertiary(br, Float64)), tertiary_turns_ratio = $(PSY.get_tertiary_turns_ratio(br))",
             )
         end
     end
@@ -643,7 +643,7 @@ 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, Float64) - im * PSY.get_impedance_reactive_power(fa, Float64)
     if !isfinite(Y)
         error(
             "Data in $(PSY.get_name(fa)) is incorrect. Y = $(Y)",
@@ -821,8 +821,8 @@ function Ybus(
     breaker_switches = Vector{PSY.DiscreteControlledACBranch}()
     for br in PSY.get_components(PSY.DiscreteControlledACBranch, sys)
         !PSY.get_available(br) && continue
-        r = PSY.get_r(br)
-        x = PSY.get_x(br)
+        r = PSY.get_r(br, Float64)
+        x = PSY.get_x(br, Float64)
         status = PSY.get_branch_status(br)
         if status == PSY.DiscreteControlledBranchStatus.CLOSED
             if r == 0.0 && x < ZERO_IMPEDANCE_LINE_REACTANCE_THRESHOLD
diff --git a/test/test_A_BA_ABA_with_radial_lines.jl b/test/test_A_BA_ABA_with_radial_lines.jl
index 3ec326464..19c3900f0 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, Float64)
         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, Float64)
         end
         power_injection =
             deepcopy(bus_activepower_injection - bus_activepower_withdrawals)
diff --git a/test/test_lodf_with_radial_branches.jl b/test/test_lodf_with_radial_branches.jl
index 89da0db80..5b319cff9 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, Float64)
         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, Float64)
         end
         power_injection =
             deepcopy(bus_activepower_injection - bus_activepower_withdrawals)
diff --git a/test/test_ptdf.jl b/test/test_ptdf.jl
index 749cbf7a0..9109de23f 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, Float64),
             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 2957d9d03..34efe3c45 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, Float64)
         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, Float64)
         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, Float64)
             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, Float64)
             end
             power_injection =
                 deepcopy(bus_activepower_injection - bus_activepower_withdrawals)

From abde495c3026e550496ce417ac216a74149841b3 Mon Sep 17 00:00:00 2001
From: Luke Kiernan <lkiernan@nrel.gov>
Date: Tue, 26 May 2026 20:57:17 -0600
Subject: [PATCH 2/5] Migrate to unitful PowerSystems getter API; pin
 units-branch sources

Replace the placeholder `Float64` getter target with explicit unit systems:
`SU` (system base) for impedances, admittances, and shunts used in matrix
construction, and `DU` (device base) for equivalent ratings, which must work
on detached branches. Thread `units::IS.AbstractUnitSystem` through the
`get_series_susceptance` accessors (PSY now requires it).

Pin InfrastructureSystems/PowerSystems (and test deps PowerSystemCaseBuilder/
PowerFlowFileParser) to their units-aware branches via `[sources]` so CI can
resolve the as-yet-unreleased upstream changes.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 Project.toml                            |  4 ++++
 src/BA_ABA_matrices.jl                  |  6 ++---
 src/BranchesParallel.jl                 | 11 +++++----
 src/BranchesSeries.jl                   | 13 +++++-----
 src/ThreeWindingTransformerWinding.jl   | 31 +++++++++++++-----------
 src/Ybus.jl                             | 32 +++++++++++++------------
 test/Project.toml                       |  8 +++++++
 test/test_A_BA_ABA_with_radial_lines.jl |  4 ++--
 test/test_equivalent_getters.jl         |  2 +-
 test/test_lodf_with_radial_branches.jl  |  4 ++--
 test/test_ptdf.jl                       |  2 +-
 test/test_ptdf_with_radial_lines.jl     |  8 +++----
 test/test_ward_reduction.jl             |  1 +
 13 files changed, 73 insertions(+), 53 deletions(-)

diff --git a/Project.toml b/Project.toml
index c3065fa2d..9ba1d1353 100644
--- a/Project.toml
+++ b/Project.toml
@@ -22,6 +22,10 @@ Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 AppleAccelerateExt = "AppleAccelerate"
 MKLPardisoExt = "Pardiso"
 
+[sources]
+InfrastructureSystems = {url = "https://github.com/NREL-Sienna/InfrastructureSystems.jl", rev = "lk/units-domain-agnostic-is4"}
+PowerSystems = {url = "https://github.com/NREL-Sienna/PowerSystems.jl", rev = "lk/units-fold-psu"}
+
 [compat]
 AppleAccelerate = "^0.4"
 DataStructures = "^0.19"
diff --git a/src/BA_ABA_matrices.jl b/src/BA_ABA_matrices.jl
index ffe7dffa3..9141f841f 100644
--- a/src/BA_ABA_matrices.jl
+++ b/src/BA_ABA_matrices.jl
@@ -109,7 +109,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
@@ -135,8 +135,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 fa6335658..012edb7e8 100644
--- a/src/BranchesParallel.jl
+++ b/src/BranchesParallel.jl
@@ -51,15 +51,15 @@ 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::BranchesParallel)
-    return sum(get_series_susceptance(branch) for branch in segment.branches)
+function get_series_susceptance(segment::BranchesParallel, units::IS.AbstractUnitSystem)
+    return sum(get_series_susceptance(branch, units) for branch in segment.branches)
 end
 
 function get_equivalent_physical_branch_parameters(bp::BranchesParallel)
@@ -85,7 +85,8 @@ This provides a conservative estimate that accounts for potential overestimation
 """
 function get_equivalent_rating(bp::BranchesParallel)
     # Sum of ratings divided by number of circuits
-    return sum(PSY.get_rating(branch, Float64) for branch in bp.branches) / length(bp.branches)
+    return sum(PSY.get_rating(branch, PSY.DU) for branch in bp.branches) /
+           length(bp.branches)
 end
 
 """
diff --git a/src/BranchesSeries.jl b/src/BranchesSeries.jl
index e91a79296..f67a00039 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
@@ -143,7 +144,7 @@ end
 Return the rating for PSY.ACTransmission branches.
 """
 function get_equivalent_rating(bs::PSY.ACTransmission)
-    return PSY.get_rating(bs, Float64)
+    return PSY.get_rating(bs, PSY.DU)
 end
 
 """
@@ -163,12 +164,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, Float64))
+    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, Float64)
+        return PSY.get_rating(branch, PSY.DU)
     end
-    return PSY.get_rating_b(branch, Float64)
+    return PSY.get_rating_b(branch, PSY.DU)
 end
 
 """
diff --git a/src/ThreeWindingTransformerWinding.jl b/src/ThreeWindingTransformerWinding.jl
index e3662bdbb..146333e38 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, Float64)
+        return PSY.get_r_primary(tfw, PSY.SU)
     elseif winding_num == 2
-        return PSY.get_r_secondary(tfw, Float64)
+        return PSY.get_r_secondary(tfw, PSY.SU)
     elseif winding_num == 3
-        return PSY.get_r_tertiary(tfw, Float64)
+        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, Float64)
+        return PSY.get_x_primary(tfw, PSY.SU)
     elseif winding_num == 2
-        return PSY.get_x_secondary(tfw, Float64)
+        return PSY.get_x_secondary(tfw, PSY.SU)
     elseif winding_num == 3
-        return PSY.get_x_tertiary(tfw, Float64)
+        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, Float64), 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, Float64)
+        PSY.get_rating_primary(tfw, PSY.DU)
     elseif winding_num == 2
-        PSY.get_rating_secondary(tfw, Float64)
+        PSY.get_rating_secondary(tfw, PSY.DU)
     elseif winding_num == 3
-        PSY.get_rating_tertiary(tfw, Float64)
+        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, Float64))
+    elseif isnothing(PSY.get_rating(tfw, PSY.DU))
         return 0.0
     else
-        return PSY.get_rating(tfw, Float64)
+        return PSY.get_rating(tfw, PSY.DU)
     end
 end
 
diff --git a/src/Ybus.jl b/src/Ybus.jl
index b7510633c..9ec35a939 100644
--- a/src/Ybus.jl
+++ b/src/Ybus.jl
@@ -379,16 +379,16 @@ function add_branch_entries_to_indexing_maps!(
 end
 
 _get_shunt(br::PSY.ACTransmission, node::Symbol) =
-    PSY.get_g(br, Float64)[node] + 1im * PSY.get_b(br, Float64)[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 a `DiscreteControlledACBranch`."""
 function ybus_branch_entries(br::PSY.ACTransmission)
-    Y_l = (1 / (PSY.get_r(br, Float64) + PSY.get_x(br, Float64) * 1im))
+    Y_l = (1 / (PSY.get_r(br, PSY.SU) + PSY.get_x(br, PSY.SU) * 1im))
     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, Float64)), x = $(PSY.get_x(br, Float64))",
+            "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
@@ -422,10 +422,10 @@ _get_tap(br::PSY.TwoWindingTransformer) = PSY.get_tap(br)
 
 """Ybus entries for a `Transformer2W`, `TapTransformer`, or `PhaseShiftingTransformer`."""
 function ybus_branch_entries(br::PSY.TwoWindingTransformer)
-    Y_t = 1 / (PSY.get_r(br, Float64) + PSY.get_x(br, Float64) * 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, Float64)
+    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(
@@ -443,39 +443,39 @@ function ybus_branch_entries(tp::ThreeWindingTransformerWinding)
     br = get_transformer(tp)
     winding_number = get_winding_number(tp)
     if winding_number == 1
-        Y_t = 1 / (PSY.get_r_primary(br, Float64) + PSY.get_x_primary(br, Float64) * 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, Float64) + im * PSY.get_b(br, Float64)
+        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, Float64)), x_p = $(PSY.get_x_primary(br, Float64)), 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, Float64) + PSY.get_x_secondary(br, Float64) * 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, Float64)), x_s = $(PSY.get_x_secondary(br, Float64)), 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, Float64) + PSY.get_x_tertiary(br, Float64) * 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, Float64)), x_t = $(PSY.get_x_tertiary(br, Float64)), 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
@@ -643,7 +643,9 @@ function _ybus!(
     nr::NetworkReductionData,
 )
     bus_no = get_bus_index(fa, num_bus, nr)
-    Y = PSY.get_impedance_active_power(fa, Float64) - im * PSY.get_impedance_reactive_power(fa, Float64)
+    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)",
@@ -821,8 +823,8 @@ function Ybus(
     breaker_switches = Vector{PSY.DiscreteControlledACBranch}()
     for br in PSY.get_components(PSY.DiscreteControlledACBranch, sys)
         !PSY.get_available(br) && continue
-        r = PSY.get_r(br, Float64)
-        x = PSY.get_x(br, Float64)
+        r = PSY.get_r(br, PSY.SU)
+        x = PSY.get_x(br, PSY.SU)
         status = PSY.get_branch_status(br)
         if status == PSY.DiscreteControlledBranchStatus.CLOSED
             if r == 0.0 && x < ZERO_IMPEDANCE_LINE_REACTANCE_THRESHOLD
diff --git a/test/Project.toml b/test/Project.toml
index 8196a9053..e71baf919 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 = "lk/units-domain-agnostic-is4"}
+PowerFlowFileParser = {url = "https://github.com/luke-kiernan/PowerFlowFileParser.jl", rev = "lk/units-fixes"}
+PowerSystemCaseBuilder = {url = "https://github.com/NREL-Sienna/PowerSystemCaseBuilder.jl", rev = "lk/units-table-driven"}
+PowerSystems = {url = "https://github.com/NREL-Sienna/PowerSystems.jl", rev = "lk/units-fold-psu"}
+
 [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 19c3900f0..39d4c575b 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, Float64)
+            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, Float64)
+            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_equivalent_getters.jl b/test/test_equivalent_getters.jl
index a2f615772..3ed6998b7 100644
--- a/test/test_equivalent_getters.jl
+++ b/test/test_equivalent_getters.jl
@@ -79,7 +79,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 5b319cff9..ac3cf8c8a 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, Float64)
+            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, Float64)
+            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_ptdf.jl b/test/test_ptdf.jl
index 9109de23f..ae8c68340 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, Float64),
+            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 34efe3c45..b4b4dd0ed 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, Float64)
+            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, Float64)
+            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, Float64)
+                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, Float64)
+                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 75506d371..6f948b715 100644
--- a/test/test_ward_reduction.jl
+++ b/test/test_ward_reduction.jl
@@ -113,6 +113,7 @@ end
     )
     existing_line_susceptance = PSY.get_series_susceptance(
         ptdf_2.network_reduction_data.direct_branch_map[(101, 102)],
+        PSY.SU,
     )
     ward_line_susceptance =
         1 / imag(1 / (ptdf_2.network_reduction_data.added_branch_map[(101, 102)]))

From f7c1bfa0a85c9f44175d2ad7742d5a2c344fe9d0 Mon Sep 17 00:00:00 2001
From: Luke Kiernan <lkiernan@nrel.gov>
Date: Fri, 29 May 2026 12:49:03 -0600
Subject: [PATCH 3/5] Re-pin sources to IS4/psy6 release branches; fix units
 call sites

Re-point [sources] from the units feature branches to the next-release
branches now that they are merged upstream:
  - InfrastructureSystems: lk/units-domain-agnostic-is4 -> IS4
  - PowerSystems:          lk/units-fold-psu            -> psy6
  - PowerSystemCaseBuilder lk/units-table-driven        -> psy6
  - PowerFlowFileParser:   fork lk/units-fixes          -> NREL-Sienna main
    (the fix is merged to main but predates the only release, v0.1.0)

Fix two breakages surfaced by the bump to psy6:
  - virtual_lodf_calculations.jl: pass PSY.SU to get_series_susceptance
    in _extract_branch_susceptances_by_arc (missed call sites; the
    accessor now requires an explicit unit system).
  - test_arc_types_and_reductions.jl: attach the MixedBranchesParallel
    test's Line/TapTransformer to a System so the system-base getters in
    ybus_branch_entries resolve (psy6 now throws on detached components).
  - test_connectivity.jl: relax the VirtualPTDF/VirtualLODF cross-method
    equality checks from == to isapprox; the on-demand paths now differ
    by ULPs on the IS4/psy6 stack.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 Project.toml                          | 4 ++--
 src/virtual_lodf_calculations.jl      | 4 ++--
 test/Project.toml                     | 8 ++++----
 test/test_arc_types_and_reductions.jl | 8 ++++++++
 test/test_connectivity.jl             | 4 ++--
 5 files changed, 18 insertions(+), 10 deletions(-)

diff --git a/Project.toml b/Project.toml
index d8842e8a3..4cb4f39c4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -21,8 +21,8 @@ Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 MKLPardisoExt = "Pardiso"
 
 [sources]
-InfrastructureSystems = {url = "https://github.com/NREL-Sienna/InfrastructureSystems.jl", rev = "lk/units-domain-agnostic-is4"}
-PowerSystems = {url = "https://github.com/NREL-Sienna/PowerSystems.jl", rev = "lk/units-fold-psu"}
+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"
diff --git a/src/virtual_lodf_calculations.jl b/src/virtual_lodf_calculations.jl
index 3417393a4..10749e5a8 100644
--- a/src/virtual_lodf_calculations.jl
+++ b/src/virtual_lodf_calculations.jl
@@ -240,12 +240,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 e71baf919..ecdeb90cc 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -21,10 +21,10 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TimeSeries = "9e3dc215-6440-5c97-bce1-76c03772f85e"
 
 [sources]
-InfrastructureSystems = {url = "https://github.com/NREL-Sienna/InfrastructureSystems.jl", rev = "lk/units-domain-agnostic-is4"}
-PowerFlowFileParser = {url = "https://github.com/luke-kiernan/PowerFlowFileParser.jl", rev = "lk/units-fixes"}
-PowerSystemCaseBuilder = {url = "https://github.com/NREL-Sienna/PowerSystemCaseBuilder.jl", rev = "lk/units-table-driven"}
-PowerSystems = {url = "https://github.com/NREL-Sienna/PowerSystems.jl", rev = "lk/units-fold-psu"}
+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_arc_types_and_reductions.jl b/test/test_arc_types_and_reductions.jl
index c734d29ab..059d71f14 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 0fe8f01cf..ed4e21cd1 100644
--- a/test/test_connectivity.jl
+++ b/test/test_connectivity.jl
@@ -134,10 +134,10 @@ end
         @test lodf_1[i, j] == lodf_2[i, j] == lodf_3[i, j]
     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 vptdf_1[i, j] ≈ vptdf_2[i, j] ≈ vptdf_3[i, j]
     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 vlodf_1[i, j] ≈ vlodf_2[i, j] ≈ vlodf_3[i, j]
     end
 end
 

From f4dc6fab996a6e52eb0837911d54f9ddc5d01f9c Mon Sep 17 00:00:00 2001
From: Luke Kiernan <lkiernan@nrel.gov>
Date: Fri, 29 May 2026 14:30:27 -0600
Subject: [PATCH 4/5] Complete units migration in connectivity/contingency/ward
 code for psy6

Finish threading the unit-aware PSY getter API through the subsystems the
original migration missed, surfaced by the bump to the psy6/IS4 release
branches:

- connectivity_checks.jl: reimplement find_connected_components with PNM's
  own union-find instead of PowerModels.calc_connected_components, which
  PowerSystems no longer re-exports (PowerModels dependency dropped).
- common.jl, network_modification.jl: pass PSY.SU to get_series_susceptance
  and get_impedance_active/reactive_power in the contingency/MODF delta
  computations (all operate on system-base, system-attached components).
- Ybus.jl, BranchesSeries.jl: the synthetic ward-equivalent
  GenericArcImpedance is detached by design and stores system-base impedance;
  read it back with device base (DU, identity) since system base would need a
  system base power a detached component cannot resolve.
- BranchesParallel.jl: document why get_impedance_averaged_rating uses system
  base (consistent weighting across the group; requires attached branches).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 src/BranchesParallel.jl     |  5 +++++
 src/BranchesSeries.jl       |  5 +++--
 src/Ybus.jl                 |  7 +++++--
 src/common.jl               |  6 +++---
 src/connectivity_checks.jl  | 40 +++++++++++++++++++++----------------
 src/network_modification.jl | 18 ++++++++---------
 6 files changed, 48 insertions(+), 33 deletions(-)

diff --git a/src/BranchesParallel.jl b/src/BranchesParallel.jl
index 2e1b203f9..909f0e509 100644
--- a/src/BranchesParallel.jl
+++ b/src/BranchesParallel.jl
@@ -132,6 +132,11 @@ 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)
+    # 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(
diff --git a/src/BranchesSeries.jl b/src/BranchesSeries.jl
index 38239fc87..1fa87c783 100644
--- a/src/BranchesSeries.jl
+++ b/src/BranchesSeries.jl
@@ -155,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
 
 """
@@ -190,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/Ybus.jl b/src/Ybus.jl
index cbce5b47e..b23ae42e2 100644
--- a/src/Ybus.jl
+++ b/src/Ybus.jl
@@ -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
diff --git a/src/common.jl b/src/common.jl
index 660d0779e..16c48eac2 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 3ffe5e3da..7a8ce2542 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 08b94b076..b11e33095 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

From e3b7ed9749d3883f8591771fabb39137f7fff8b1 Mon Sep 17 00:00:00 2001
From: Luke Kiernan <lkiernan@nrel.gov>
Date: Fri, 29 May 2026 14:30:52 -0600
Subject: [PATCH 5/5] Adapt tests to the psy6 unit-aware getter/setter API

- test_connectivity.jl: compare the dense PTDF/LODF cross-method checks with
  isapprox + an absolute tolerance; the dense and on-demand paths drift by a
  few ULPs run-to-run on the IS4/psy6 stack and many entries are physically
  zero (relative tolerance is meaningless there).
- test_equivalent_getters.jl: attach the parallel branches to the system
  (skip_validation) so the system-base susceptance weighting can resolve.
- test_modf_lodf_reductions.jl: clone transformers via deepcopy + fresh
  identity instead of re-reading each field, and pass PSY.SU to
  get_series_susceptance on attached branches.
- test_ward_reduction.jl: read the detached ward-equivalent's susceptance with
  device base (DU).
- test_ybus_reductions.jl: setters now require unitful values (val * PSY.SU).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 test/test_connectivity.jl         | 17 ++++++++++++----
 test/test_equivalent_getters.jl   |  6 ++++++
 test/test_modf_lodf_reductions.jl | 32 ++++++++++++-------------------
 test/test_ward_reduction.jl       |  3 +++
 test/test_ybus_reductions.jl      | 16 ++++++++--------
 5 files changed, 42 insertions(+), 32 deletions(-)

diff --git a/test/test_connectivity.jl b/test/test_connectivity.jl
index ed4e21cd1..92354fddc 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 95a77dda2..710bc4811 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])
 
diff --git a/test/test_modf_lodf_reductions.jl b/test/test_modf_lodf_reductions.jl
index 8a3caf8d7..ed28e230e 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_ward_reduction.jl b/test/test_ward_reduction.jl
index ac004daeb..8b9210e21 100644
--- a/test/test_ward_reduction.jl
+++ b/test/test_ward_reduction.jl
@@ -133,8 +133,11 @@ end
         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 d95080048..9a8cd047c 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