Corner cases catches

src/common.jl

@@ -307,25 +307,22 @@
 """
 !!! MISSING DOCUMENTATION !!!
 """
-function assign_reference_buses(
+function assign_reference_buses!(
     subnetworks::Dict{Int, Set{Int}},
-    ref_bus_positions::Set{Int},
+    ref_buses::Vector{Int},
 )
-    if isempty(ref_bus_positions) || length(ref_bus_positions) != length(subnetworks)
+    if isempty(ref_buses) || length(ref_buses) != length(subnetworks)
         @warn "The reference bus positions are not consistent with the subnetworks. References buses will be assigned arbitrarily"
         return deepcopy(subnetworks)
     end
     bus_groups = Dict{Int, Set{Int}}()
     for (bus_key, subnetwork_buses) in subnetworks
-        ref_bus = intersect(ref_bus_positions, subnetwork_buses)
+        ref_bus = intersect(ref_buses, subnetwork_buses)
         if length(ref_bus) == 1
             bus_groups[first(ref_bus)] = pop!(subnetworks, bus_key)
-            continue
         elseif length(ref_bus) == 0
-            @warn "No reference bus in the subnetwork associated with bus $bus_key. References buses will be assigned arbitrarily"
-            return subnetworks
+            @warn "No reference bus in the subnetwork associated with bus $bus_key. Reference bus assigned arbitrarily"
         elseif length(ref_bus) > 1
-            # TODO: still to implement
             error(
                 "More than one reference bus in the subnetwork associated with bus $bus_key",
             )
@@ -336,6 +333,15 @@
     return bus_groups
 end
 
+function assign_reference_buses!(
+    subnetworks::Dict{Int, Set{Int}},
+    ref_bus_positions::Set{Int},
+    bus_lookup::Dict{Int, Int},
+)
+    ref_buses = [k for (k, v) in bus_lookup if v in ref_bus_positions]
+    assign_reference_buses!(subnetworks, ref_buses)
+end
+
 """
 Finds the subnetworks present in the considered System. This is evaluated by taking
 a the ABA or Adjacency Matrix.

src/network_radial_reduction.jl

@@ -99,6 +99,13 @@
             reverse_bus_map,
         )
         new_parent_bus_number = reverse_bus_map[new_parent_val]
+        # This check is meant to capture cases of a full radial network which can happen in
+        # system with small islands that represent larger interconnected areas.
+        if length(SparseArrays.nzrange(A, new_parent_val)) < 2
+            @warn "Bus $parent_bus_number Parent $new_parent_bus_number is a leaf node. Indicating there is an island."
+            push!(bus_reduction_map_index[parent_bus_number], new_parent_bus_number)
+            return
+        end
         new_set = push!(pop!(bus_reduction_map_index, parent_bus_number), parent_bus_number)
         union!(bus_reduction_map_index[new_parent_bus_number], new_set)
         _new_parent(
@@ -120,7 +127,11 @@
     reverse_line_map::Dict{Int64, String},
     radial_branches::Set{String},
     reverse_bus_map::Dict{Int, Int},
+    ref_bus_positions::Set{Int},
 )
+    if j ∈ ref_bus_positions
+        return
+    end
     j_bus_number = reverse_bus_map[j]
     pop!(bus_reduction_map_index, j_bus_number)
     reducion_set = Set{Int}(j_bus_number)
@@ -134,6 +145,9 @@
     )
     parent_bus_number = reverse_bus_map[parent]
     union!(bus_reduction_map_index[parent_bus_number], reducion_set)
+    if parent ∈ ref_bus_positions
+        return
+    end
     _new_parent(
         A,
         parent,
@@ -181,9 +195,6 @@
     reverse_bus_map = Dict(reverse(kv) for kv in bus_map)
     bus_reduction_map_index = Dict{Int, Set{Int}}(k => Set{Int}() for k in keys(bus_map))
     Threads.@threads for j in 1:buscount
-        if j ∈ ref_bus_positions
-            continue
-        end
         if length(SparseArrays.nzrange(A, j)) == 1
             lock(lk) do
                 _reverse_search(
@@ -193,6 +204,7 @@
                     reverse_line_map,
                     radial_branches,
                     reverse_bus_map,
+                    ref_bus_positions,
                 )
             end
         end

src/ptdf_calculations.jl

@@ -400,7 +400,7 @@ function PTDF(
     subnetworks = find_subnetworks(M, bus_ax)
     if length(subnetworks) > 1
         @info "Network is not connected, using subnetworks"
-        subnetworks = assign_reference_buses(subnetworks, ref_bus_positions)
+        subnetworks = assign_reference_buses!(subnetworks, ref_bus_positions, bus_ax_ref)
     end
     look_up = (bus_ax_ref, make_ax_ref(line_ax))
     S, _ = _buildptdf(

src/system_utils.jl

@@ -11,15 +11,16 @@ function _find_subnetworks(sys::PSY.System)
     branches = get_ac_branches(sys)
     @info "Validating connectivity with depth first search (network traversal)"
     M, bus_lookup = calculate_adjacency(branches, buses)
-    slack_positions = find_slack_positions(buses, bus_lookup)
-    return find_subnetworks(M, PSY.get_number.(buses)), slack_positions
+    ref_positions = find_slack_positions(buses, bus_lookup)
+    ref_buses = [k for (k, v) in bus_lookup if v in ref_positions]
+    return find_subnetworks(M, PSY.get_number.(buses)), ref_buses
 end
 
 """
 Finds the subnetworks in a system using Depth First Search (DFS). Returns a dictionary keyed
 by the reference bus of the subnetworks if they exist
 """
 function find_subnetworks(sys::PSY.System)
-    sbn, ref_bus_positions = _find_subnetworks(sys)
-    return assign_reference_buses(sbn, ref_bus_positions)
+    sbn, ref_buses = _find_subnetworks(sys)
+    return assign_reference_buses!(sbn, ref_buses)
 end

src/virtual_lodf_calculations.jl

@@ -171,7 +171,7 @@ function VirtualLODF(
     # check subnetworks
     if length(subnetworks) > 1
         @info "Network is not connected, using subnetworks"
-        subnetworks = assign_reference_buses(subnetworks, ref_bus_positions)
+        subnetworks = assign_reference_buses!(subnetworks, ref_bus_positions, bus_ax_ref)
     end
     # get diagonal of PTDF
     temp_data = zeros(length(bus_ax))

src/virtual_ptdf_calculations.jl

@@ -120,7 +120,7 @@ function VirtualPTDF(
     subnetworks = find_subnetworks(M, bus_ax)
     if length(subnetworks) > 1
         @info "Network is not connected, using subnetworks"
-        subnetworks = assign_reference_buses(subnetworks, ref_bus_positions)
+        subnetworks = assign_reference_buses!(subnetworks, ref_bus_positions, bus_ax_ref)
     end
     temp_data = zeros(length(bus_ax))