Merge pull request #67 from NREL-Sienna/afc/radial_branches_matrices

afc/radial_branches_matrices

src/BA_ABA_matrices.jl

@@ -3,7 +3,7 @@ Structure containing the BA matrix and other relevant data.
 
 # Arguments
 - `data::SparseArrays.SparseMatrixCSC{Float64, Int}`:
-        the transposed BA matrix coming from the product between the Incidence 
+        the transposed BA matrix coming from the product between the Incidence
         Matrix A and the Matrix of Susceptance B
 - `axes<:NTuple{2, Dict}`:
         Tuple containing two vectors, the first one contains the names of each
@@ -14,30 +14,46 @@ Structure containing the BA matrix and other relevant data.
         Tuple containing 2 Dictionaries mapping the number of rows and columns
         with the names of buses and branches
 - `ref_bus_positions::Set{Int}`:
-        Vector containing the indexes of the columns of the BA matrix corresponding
+        Set containing the indexes of the columns of the BA matrix corresponding
         to the refence buses
+- `radial_branches::RadialBranches`:
+        Structure containing the radial branches and leaf buses that were removed
+        while evaluating the matrix
 """
 struct BA_Matrix{Ax, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{Float64}
     data::SparseArrays.SparseMatrixCSC{Float64, Int}
     axes::Ax
     lookup::L
     ref_bus_positions::Set{Int}
+    radial_branches::RadialBranches
 end
 
 """
 Build the BA matrix from a given System
+
+# Arguments
+- `sys::PSY.System`:
+        PSY system for which the matrix is constructed
+- `reduce_radial_branches::Bool`:
+        if True the matrix is build considering radial branches removed from
+        the system
 """
-function BA_Matrix(sys::PSY.System)
-    branches = get_ac_branches(sys)
-    buses = get_buses(sys)
+function BA_Matrix(sys::PSY.System; reduce_radial_branches::Bool = false)
+    if reduce_radial_branches
+        rb = RadialBranches(IncidenceMatrix(sys))
+    else
+        rb = RadialBranches()
+    end
+    branches = get_ac_branches(sys, rb.radial_branches)
+    buses = get_buses(sys, rb.bus_reduction_map)
     ref_bus_positions = find_slack_positions(buses)
     bus_lookup = make_ax_ref(buses)
     line_ax = [PSY.get_name(branch) for branch in branches]
     bus_ax = [PSY.get_number(bus) for bus in setdiff(buses, ref_bus_positions)]
     axes = (bus_ax, line_ax)
     lookup = (make_ax_ref(bus_ax), make_ax_ref(line_ax))
     data = calculate_BA_matrix(branches, bus_lookup)
-    return BA_Matrix(data, axes, lookup, ref_bus_positions)
+    return BA_Matrix(data, axes, lookup, ref_bus_positions, rb)
 end
 
 """
@@ -48,7 +64,7 @@ Structure containing the ABA matrix and other relevant data.
         the ABA matrix coming from the product between the Incidence Matrix A and
         the Matrix BA.
 - `axes<:NTuple{2, Dict}`:
-        Tuple containing two identical vectors, both containing the number of 
+        Tuple containing two identical vectors, both containing the number of
         each bus of the network (each one related to a row/column of the Matrix
         in "data"), excluding the slack buses.
 - `lookup<:NTuple{2, Dict}`:
@@ -59,6 +75,9 @@ Structure containing the ABA matrix and other relevant data.
         to the refence buses
 - `K<:Union{Nothing, KLU.KLUFactorization{Float64, Int}}`:
         either nothing or a container for KLU factorization matrices (LU factorization)
+- `radial_branches::RadialBranches`:
+        Structure containing the radial branches and leaf buses that were removed
+        while evaluating the matrix
 """
 struct ABA_Matrix{
     Ax,
@@ -70,6 +89,7 @@ struct ABA_Matrix{
     lookup::L
     ref_bus_positions::Set{Int}
     K::F
+    radial_branches::RadialBranches
 end
 
 """
@@ -82,9 +102,19 @@ Builds the ABA matrix from a System
 # Keyword arguments
 - `factorize`: if true populates ABA_Matrix.K with KLU factorization matrices
 """
-function ABA_Matrix(sys::PSY.System; factorize = false)
-    branches = get_ac_branches(sys)
-    buses = get_buses(sys)
+function ABA_Matrix(
+    sys::PSY.System;
+    factorize = false,
+    reduce_radial_branches::Bool = false,
+)
+    if reduce_radial_branches
+        rb = RadialBranches(IncidenceMatrix(sys))
+    else
+        rb = RadialBranches()
+    end
+
+    branches = get_ac_branches(sys, rb.radial_branches)
+    buses = get_buses(sys, rb.bus_reduction_map)
     bus_lookup = make_ax_ref(buses)
 
     A, ref_bus_positions = calculate_A_matrix(branches, buses)
@@ -107,6 +137,7 @@ function ABA_Matrix(sys::PSY.System; factorize = false)
         lookup,
         ref_bus_positions,
         K,
+        rb,
     )
 end
 
@@ -124,6 +155,7 @@ function factorize(ABA::ABA_Matrix{Ax, L, Nothing}) where {Ax, L <: NTuple{2, Di
         deepcopy(ABA.lookup),
         deepcopy(ABA.ref_bus_positions),
         klu(ABA.data),
+        deepcopy(ABA.radial_branches),
     )
     return ABA_lu
 end
@@ -149,7 +181,7 @@ function Base.getindex(
     return A.data[bus_number, line_number]
 end
 
-# get_index functions: ABA_Matrix stores a square matrix whose number of rows 
+# get_index functions: ABA_Matrix stores a square matrix whose number of rows
 # and column is equal to the number of the system's buses minus the slack ones,
 # NOTE: bus_1, bus_2 are bus numbers not row and column indices!
 function Base.getindex(A::ABA_Matrix, bus_1, bus_2)

src/PowerNetworkMatrices.jl

@@ -48,11 +48,11 @@ import Pardiso
 
 # network calculations
 include("PowerNetworkMatrix.jl")
-include("BA_ABA_matrices.jl")
 include("incedence_matrix.jl")
 include("adjacency_matrix.jl")
-include("common.jl")
 include("radial_braches.jl")
+include("common.jl")
+include("BA_ABA_matrices.jl")
 include("definitions.jl")
 include("ptdf_calculations.jl")
 include("ybus_calculations.jl")

src/common.jl

@@ -13,7 +13,10 @@ end
 """
 Gets the AC branches from a given Systems.
 """
-function get_ac_branches(sys::PSY.System)
+function get_ac_branches(
+    sys::PSY.System,
+    radial_branches::Set{String} = Set{String}(),
+)::Vector{PSY.ACBranch}
     collection = Vector{PSY.ACBranch}()
     for br in PSY.get_components(PSY.get_available, PSY.ACBranch, sys)
         arc = PSY.get_arc(br)
@@ -31,7 +34,9 @@ function get_ac_branches(sys::PSY.System)
                 ),
             )
         end
-        _add_to_collection!(collection, br)
+        if PSY.get_name(br) ∉ radial_branches
+            _add_to_collection!(collection, br)
+        end
     end
     return sort!(collection;
         by = x -> (PSY.get_number(PSY.get_arc(x).from), PSY.get_number(PSY.get_arc(x).to)),
@@ -41,18 +46,35 @@ end
 """
 Gets the non-isolated buses from a given System
 """
-function get_buses(sys::PSY.System)::Vector{PSY.ACBus}
-    return sort!(
-        collect(
-            PSY.get_components(
-                x -> PSY.get_bustype(x) != ACBusTypes.ISOLATED,
-                PSY.ACBus,
-                sys,
-            ),
-        );
-        by = x -> PSY.get_number(x),
-    )
+function get_buses(
+    sys::PSY.System,
+    bus_reduction_map::Dict{Int64, Set{Int64}} = Dict{Int64, Set{Int64}}(),
+)::Vector{PSY.ACBus}
+    leaf_buses = Set{PSY.Int64}()
+    if !isempty(bus_reduction_map)
+        for vals in values(bus_reduction_map)
+            union!(leaf_buses, vals)
+        end
+    end
+
+    count_i = 1
+    all_buses = PSY.get_components(PSY.ACBus, sys)
+    buses = Vector{PSY.ACBus}(undef, length(all_buses))
+    for b in all_buses
+        if PSY.get_bustype(b) == ACBusTypes.ISOLATED
+            continue
+        end
+
+        if PSY.get_number(b) ∈ leaf_buses
+            continue
+        end
+        buses[count_i] = b
+        count_i += 1
+    end
+
+    return sort!(deleteat!(buses, count_i:length(buses)); by = x -> PSY.get_number(x))
 end
+
 """
 Gets the indices  of the reference (slack) buses.
 NOTE:
@@ -157,11 +179,17 @@ function calculate_adjacency(
     a = SparseArrays.spzeros(Int8, buscount, buscount)
 
     for b in branches
-        (fr_b, to_b) = get_bus_indices(b, bus_lookup)
+        fr_b, to_b = get_bus_indices(b, bus_lookup)
         a[fr_b, to_b] = 1
         a[to_b, fr_b] = -1
-        a[fr_b, fr_b] = 1
-        a[to_b, to_b] = 1
+    end
+
+    # If a line is disconnected needs to check for the buses correctly
+    for i in 1:buscount
+        if PSY.get_bustype(buses[i]) == ACBusTypes.ISOLATED
+            continue
+        end
+        a[i, i] = 1
     end
 
     # Return both for type stability
@@ -191,6 +219,10 @@ function calculate_BA_matrix(
         (fr_b, to_b) = get_bus_indices(b, bus_lookup)
         b_val = PSY.get_series_susceptance(b)
 
+        if !isfinite(b_val)
+            error("Invalid value for branch $(PSY.summary(b)), $b_val")
+        end
+
         push!(BA_I, fr_b)
         push!(BA_J, ix)
         push!(BA_V, b_val)
@@ -275,7 +307,7 @@ end
 """
 !!! MISSING DOCUMENTATION !!!
 """
-function assing_reference_buses(
+function assign_reference_buses(
     subnetworks::Dict{Int, Set{Int}},
     ref_bus_positions::Set{Int},
 )
@@ -320,7 +352,7 @@ function find_subnetworks(M::SparseArrays.SparseMatrixCSC, bus_numbers::Vector{I
     subnetworks = Dict{Int, Set{Int}}()
     for (ix, bus_number) in enumerate(bus_numbers)
         neighbors = SparseArrays.nzrange(M, ix)
-        if length(neighbors) < 1
+        if length(neighbors) <= 1
             @warn "Bus $bus_number is islanded"
             subnetworks[bus_number] = Set{Int}(bus_number)
             continue

src/incedence_matrix.jl

@@ -11,7 +11,10 @@ Incidence matrix: shows connection between buses, defining lines
         Tuple containing two dictionaries, the first mapping the branches
         and buses with their enumerated indexes.
 - `ref_bus_positions::Set{Int}`:
-        vector containing the indices of the reference slack buses.
+        Vector containing the indices of the reference slack buses.
+- `radial_branches::RadialBranches`:
+        Structure containing the radial branches and leaf buses that were removed
+        while evaluating the matrix
 """
 struct IncidenceMatrix{Ax, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{Int8}
     data::SparseArrays.SparseMatrixCSC{Int8, Int}
@@ -26,16 +29,52 @@ get_lookup(A::IncidenceMatrix) = A.lookup
 get_slack_position(A::IncidenceMatrix) = A.ref_bus_positions
 
 """
-Builds the Incidence matrix by evaluating the actual matrix and other relevant
+Builds the Incidence matrix of the system by evaluating the actual matrix and other relevant
 values.
+
+# Arguments
+- `sys::PSY.System`:
+        the PowerSystem system to consider
+- `reduce_radial_branches::Bool`:
+        if True the matrix will be evaluated discarding
+        all the radial branches and leaf buses (optional, default value is false)
 """
-function IncidenceMatrix(sys::PSY.System)
+function IncidenceMatrix(
+    sys::PSY.System,
+)
+    data, axes, lookup, ref_bus_positions = evaluate_A_matrix_values(sys)
+    return IncidenceMatrix(data, axes, lookup, ref_bus_positions)
+end
+
+"""
+Builds the Incidence matrix of the system by evaluating the actual matrix and other relevant
+values.
+
+# Arguments
+- `sys::PSY.System`: the PowerSystem system to consider
+- `radial_branches::RadialBranches`:
+        Structure containing the radial branches and leaf buses that were removed
+        while evaluating the matrix
+"""
+function evaluate_A_matrix_values(
+    sys::PSY.System,
+)
     branches = get_ac_branches(sys)
     buses = get_buses(sys)
     line_ax = [PSY.get_name(branch) for branch in branches]
     bus_ax = [PSY.get_number(bus) for bus in buses]
     data, ref_bus_positions = calculate_A_matrix(branches, buses)
     axes = (line_ax, bus_ax)
     lookup = (make_ax_ref(line_ax), make_ax_ref(bus_ax))
-    return IncidenceMatrix(data, axes, lookup, ref_bus_positions)
+    return data, axes, lookup, ref_bus_positions
+end
+
+function reduce_A_matrix(
+    A::IncidenceMatrix,
+    bus_reduction_map::Dict{Int, Set{Int}},
+    meshed_branches::Set{String},
+)
+    branch_ixs = sort!([A.lookup[1][k] for k in meshed_branches])
+    bus_ixs = sort!([A.lookup[2][k] for k in keys(bus_reduction_map)])
+    return A.data[branch_ixs, bus_ixs]
 end