Jd/radial branches

Project.toml

@@ -18,8 +18,10 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 DocStringExtensions = "~0.8, ~0.9"
 HDF5 = "0.17"
 InfrastructureSystems = "^1.20"
+LinearAlgebra = "1"
 MKL = "0.6"
 KLU = "~0.4"
 Pardiso = "~0.5.4"
 PowerSystems = "3"
+SparseArrays = "1"
 julia = "^1.6"

src/PowerNetworkMatrices.jl

@@ -8,10 +8,13 @@ export find_subnetworks
 export from_hdf5
 export get_ptdf_data
 export get_lodf_data
+export get_bus_reduction_map
+export get_radial_branches
 export IncidenceMatrix
 export is_factorized
 export LODF
 export PTDF
+export RadialBranches
 export to_hdf5
 export validate_connectivity
 export VirtualLODF
@@ -49,6 +52,7 @@ include("BA_ABA_matrices.jl")
 include("incedence_matrix.jl")
 include("adjacency_matrix.jl")
 include("common.jl")
+include("radial_braches.jl")
 include("definitions.jl")
 include("ptdf_calculations.jl")
 include("ybus_calculations.jl")

src/lodf_calculations.jl

@@ -21,6 +21,7 @@
     axes::Ax
     lookup::L
     tol::Base.RefValue{Float64}
+    radial_banches::RadialBranches
 end
 
 function _buildlodf(
@@ -262,31 +263,44 @@
         Default solver: "KLU".
 - `tol::Float64`:
         Tolerance to eliminate entries in the LODF matrix (default eps())
+- `reduce_radial_branches::Bool`:
+        True to reduce the network by simplifying the radial branches and mapping the
+        eliminate buses
 """
 function LODF(
     branches,
     buses::Vector{PSY.ACBus};
     linear_solver::String = "KLU",
     tol::Float64 = eps(),
+    reduce_radial_branches::Bool = false,
 )
 
     # get axis names
     line_ax = [branch.name for branch in branches]
     axes = (line_ax, line_ax)
-    look_up = (make_ax_ref(line_ax), make_ax_ref(line_ax))
+    line_map = make_ax_ref(line_ax)
+    look_up = (line_map, line_map)
     bus_ax = [PSY.get_number(bus) for bus in buses]
     bus_lookup = make_ax_ref(bus_ax)
     # get network matrices
     ptdf_t, a = calculate_PTDF_matrix_KLU(branches, buses, bus_lookup, Float64[])
+    if reduce_radial_branches
+        data, ref_bus_positions = calculate_A_matrix(branches, buses)
+        radial_branches = RadialBranches(data, bus_lookup, line_map, ref_bus_positions)
+    else
+        radial_branches = RadialBranches()
+    end
+
     if tol > eps()
         lodf_t = _buildlodf(a, ptdf_t, linear_solver)
-        return LODF(sparsify(lodf_t, tol), axes, look_up, Ref(tol))
+        return LODF(sparsify(lodf_t, tol), axes, look_up, Ref(tol), radial_branches)
     else
         return LODF(
             _buildlodf(a, ptdf_t, linear_solver),
             axes,
             look_up,
             Ref(tol),
+            radial_branches,
         )
     end
 end
@@ -323,12 +337,16 @@
         Linear solver to be used. Options are "Dense" and "KLU".
 - `tol::Float64`:
         Tolerance to eliminate entries in the LODF matrix (default eps()).
+- `reduce_radial_branches::Bool`:
+        True to reduce the network by simplifying the radial branches and mapping the
+        eliminate buses
 """
 function LODF(
     A::IncidenceMatrix,
     PTDFm::PTDF;
     linear_solver::String = "KLU",
     tol::Float64 = eps(),
+    reduce_radial_branches::Bool = true,
 )
     validate_linear_solver(linear_solver)
 
@@ -339,7 +357,11 @@
         )
         error(err_msg)
     end
-
+    if reduce_radial_branches
+        radial_branches = RadialBranches(A)
+    else
+        radial_branches = RadialBranches()
+    end
     ax_ref = make_ax_ref(A.axes[1])
     if tol > eps()
         lodf_t = _buildlodf(A.data, PTDFm.data, linear_solver)
@@ -348,13 +370,15 @@
             (A.axes[1], A.axes[1]),
             (ax_ref, ax_ref),
             Ref(tol),
+            radial_branches,
         )
     end
     return LODF(
         _buildlodf(A.data, PTDFm.data, linear_solver),
         (A.axes[1], A.axes[1]),
         (ax_ref, ax_ref),
         Ref(tol),
+        radial_branches,
     )
 end
 
@@ -374,16 +398,25 @@
         Linear solver to be used. Options are "Dense" and "KLU".
 - `tol::Float64`:
         Tolerance to eliminate entries in the LODF matrix (default eps()).
+- `reduce_radial_branches::Bool`:
+        True to reduce the network by simplifying the radial branches and mapping the
+        eliminate buses
 """
 function LODF(
     A::IncidenceMatrix,
     ABA::ABA_Matrix,
     BA::BA_Matrix;
     linear_solver::String = "KLU",
     tol::Float64 = eps(),
+    reduce_radial_branches::Bool = false,
 )
     validate_linear_solver(linear_solver)
     ax_ref = make_ax_ref(A.axes[1])
+    if reduce_radial_branches
+        radial_branches = RadialBranches(A)
+    else
+        radial_branches = RadialBranches()
+    end
     if tol > eps()
         lodf_t = _buildlodf(A.data, ABA.K, BA.data,
             A.ref_bus_positions, linear_solver)
@@ -392,13 +425,15 @@
             (A.axes[1], A.axes[1]),
             (ax_ref, ax_ref),
             Ref(tol),
+            radial_branches,
         )
     end
     return LODF(
         _buildlodf(A.data, ABA.K, BA.data, A.ref_bus_positions, linear_solver),
         (A.axes[1], A.axes[1]),
         (ax_ref, ax_ref),
         Ref(tol),
+        radial_branches,
     )
 end
 

src/ptdf_calculations.jl

@@ -21,6 +21,9 @@
 - `tol::Base.RefValue{Float64}`:
         tolerance used for sparsifying the matrix (dropping element whose
         absolute value is below this threshold).
+- `reduce_radial_branches::Bool`:
+        True to reduce the network by simplifying the radial branches and mapping the
+        eliminate buses
 """
 struct PTDF{Ax, L <: NTuple{2, Dict}, M <: AbstractArray{Float64, 2}} <:
        PowerNetworkMatrix{Float64}
@@ -30,6 +33,16 @@
     subnetworks::Dict{Int, Set{Int}}
     ref_bus_positions::Set{Int}
     tol::Base.RefValue{Float64}
+    radial_banches::RadialBranches
+end
+
+function PTDF(data::AbstractArray{Float64, 2},
+    axes,
+    lookup::NTuple{2, Dict},
+    subnetworks::Dict{Int, Set{Int}},
+    ref_bus_positions::Set{Int},
+    tol::Base.RefValue{Float64})
+    return PTDF(data, axes, lookup, subnetworks, ref_bus_positions, tol, RadialBranches)
 end
 
 """
@@ -370,13 +383,17 @@
         Linear solver to be used. Options are "Dense", "KLU" and "MKLPardiso
 - `tol::Float64`:
         Tolerance to eliminate entries in the PTDF matrix (default eps())
+- `reduce_radial_branches::Bool`:
+        True to reduce the network by simplifying the radial branches and mapping the
+        eliminate buses
 """
 function PTDF(
     branches,
     buses::Vector{PSY.ACBus};
     dist_slack::Vector{Float64} = Float64[],
     linear_solver::String = "KLU",
-    tol::Float64 = eps())
+    tol::Float64 = eps(),
+    reduce_radial_branches::Bool = false)
     validate_linear_solver(linear_solver)
     #Get axis names
     line_ax = [PSY.get_name(branch) for branch in branches]
@@ -397,6 +414,12 @@
         dist_slack,
         linear_solver,
     )
+    if reduce_radial_branches
+        data, _ = calculate_A_matrix(branches, buses)
+        radial_branches = RadialBranches(data, bus_lookup, line_map, ref_bus_positions)
+    else
+        radial_branches = RadialBranches()
+    end
     if tol > eps()
         return PTDF(
             sparsify(S, tol),
@@ -405,9 +428,10 @@
             subnetworks,
             ref_bus_positions,
             Ref(tol),
+            radial_branches,
         )
     end
-    return PTDF(S, axes, look_up, subnetworks, ref_bus_positions, Ref(tol))
+    return PTDF(S, axes, look_up, subnetworks, ref_bus_positions, Ref(tol), radial_branches)
 end
 
 """
@@ -441,18 +465,28 @@
         Linear solver to be used. Options are "Dense", "KLU" and "MKLPardiso.
 - `tol::Float64`:
         Tolerance to eliminate entries in the PTDF matrix (default eps()).
+- `reduce_radial_branches::Bool`:
+        True to reduce the network by simplifying the radial branches and mapping the
+        eliminate buses
 """
 function PTDF(
     A::IncidenceMatrix,
     BA::BA_Matrix;
     dist_slack::Vector{Float64} = Float64[],
     linear_solver = "KLU",
-    tol::Float64 = eps())
+    tol::Float64 = eps(),
+    reduce_radial_branches::Bool = false)
     validate_linear_solver(linear_solver)
     S = _buildptdf_from_matrices(A, BA.data, dist_slack, linear_solver)
     axes = (A.axes[2], A.axes[1])
     lookup = (A.lookup[2], A.lookup[1])
     @warn "PTDF creates via other matrices doesn't compute the subnetworks"
+    if reduce_radial_branches
+        data, ref_bus_positions = calculate_A_matrix(branches, buses)
+        radial_branches = RadialBranches(data, bus_lookup, line_map, ref_bus_positions)
+    else
+        radial_branches = RadialBranches()
+    end
     if tol > eps()
         return PTDF(
             sparsify(S, tol),
@@ -461,9 +495,18 @@
             Dict{Int, Set{Int}}(),
             BA.ref_bus_positions,
             Ref(tol),
+            radial_branches,
         )
     else
-        return PTDF(S, axes, lookup, Dict{Int, Set{Int}}(), BA.ref_bus_positions, Ref(tol))
+        return PTDF(
+            S,
+            axes,
+            lookup,
+            Dict{Int, Set{Int}}(),
+            BA.ref_bus_positions,
+            Ref(tol),
+            radial_branches,
+        )
     end
 end