diff --git a/src/BA_ABA_matrices.jl b/src/BA_ABA_matrices.jl
index 2351dd24..315d6666 100644
--- a/src/BA_ABA_matrices.jl
+++ b/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,22 +14,38 @@ 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]
@@ -37,7 +53,7 @@ function BA_Matrix(sys::PSY.System)
     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)
diff --git a/src/PowerNetworkMatrices.jl b/src/PowerNetworkMatrices.jl
index 26a27e6a..407e485e 100644
--- a/src/PowerNetworkMatrices.jl
+++ b/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")
diff --git a/src/common.jl b/src/common.jl
index 1ee4fc34..f3d80ef9 100644
--- a/src/common.jl
+++ b/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
diff --git a/src/incedence_matrix.jl b/src/incedence_matrix.jl
index a90451ae..ccd312e6 100644
--- a/src/incedence_matrix.jl
+++ b/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,10 +29,36 @@ 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]
@@ -37,5 +66,15 @@ function IncidenceMatrix(sys::PSY.System)
     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
diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl
index 00897cf1..3f5f7c5c 100644
--- a/src/lodf_calculations.jl
+++ b/src/lodf_calculations.jl
@@ -14,6 +14,9 @@ of how a change in a line's flow affects the flows on other lines in the system.
 - `tol::Base.RefValue{Float64}`:
         tolerance used for sparsifying the matrix (dropping element whose
         absolute value is below this threshold).
+- `radial_branches::RadialBranches`:
+        Structure containing the radial branches and leaf buses that were removed
+        while evaluating the matrix
 """
 struct LODF{Ax, L <: NTuple{2, Dict}, M <: AbstractArray{Float64, 2}} <:
        PowerNetworkMatrix{Float64}
@@ -21,7 +24,7 @@ struct LODF{Ax, L <: NTuple{2, Dict}, M <: AbstractArray{Float64, 2}} <:
     axes::Ax
     lookup::L
     tol::Base.RefValue{Float64}
-    radial_banches::RadialBranches
+    radial_branches::RadialBranches
 end
 
 function _buildlodf(
@@ -257,22 +260,24 @@ Builds the LODF matrix given the data of branches and buses of the system.
         vector of the System AC branches
 - `buses::Vector{PSY.ACBus}`:
         vector of the System buses
+
+# Keyword Arguments
 - `linear_solver`::String
         linear solver to use for matrix evaluation.
         Available options: "KLU", "Dense" (OpenBLAS) or "MKLPardiso".
         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
+- `radial_branches::RadialBranches`:
+        Structure containing the radial branches and leaf buses that were removed
+        while evaluating the ma
 """
 function LODF(
     branches,
     buses::Vector{PSY.ACBus};
     linear_solver::String = "KLU",
     tol::Float64 = eps(),
-    reduce_radial_branches::Bool = false,
+    radial_branches::RadialBranches = RadialBranches(),
 )
 
     # get axis names
@@ -284,12 +289,6 @@ function LODF(
     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)
@@ -307,18 +306,30 @@ end
 
 """
 Builds the LODF matrix from a system.
+Note that `reduce_radial_branches` kwargs is explicitly mentioned because needed inside of the function.
 
 # Arguments
 - `sys::PSY.System`:
         Power Systems system
+
+# Keyword Arguments
+- `reduce_radial_branches::Bool=false`:
+        if True the matrix will be evaluated discarding
+        all the radial branches and leaf buses (optional, default value is false)
 """
 function LODF(
     sys::PSY.System;
+    reduce_radial_branches::Bool = false,
     kwargs...,
 )
-    branches = get_ac_branches(sys)
-    buses = get_buses(sys)
-    return LODF(branches, buses; kwargs...)
+    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)
+    return LODF(branches, buses; radial_branches = rb, kwargs...)
 end
 
 """
@@ -346,25 +357,48 @@ function LODF(
     PTDFm::PTDF;
     linear_solver::String = "KLU",
     tol::Float64 = eps(),
-    reduce_radial_branches::Bool = true,
+    reduce_radial_branches::Bool = false,
 )
     validate_linear_solver(linear_solver)
 
     if PTDFm.tol.x > 1e-15
-        err_msg = string(
-            "The argument `tol` in the PTDF matrix was set to a value dirrent than the default one.\n",
-            "The PTDF matrix used as input of the LODF matrix must have the default `tol` value.\n",
+        warn_msg = string(
+            "The argument `tol` in the PTDF matrix was set to a value different than the default one.\n",
+            "The resulting LODF can include unexpected rounding errors.\n",
         )
-        error(err_msg)
+        @warn(warn_msg)
+        PTDFm_data = Matrix(PTDFm.data)
+    else
+        PTDFm_data = PTDFm.data
     end
+
     if reduce_radial_branches
-        radial_branches = RadialBranches(A)
+        if !isempty(PTDFm.radial_branches)
+            radial_branches = PTDFm.radial_branches
+            @info "Non-empty `radial_branches` field found in PTDF matrix. LODF is evaluated considering radial branches and leaf nodes removed."
+        else
+            error("PTDF has empty `radial_branches` field.")
+        end
+        A_matrix = reduce_A_matrix(
+            A,
+            radial_branches.bus_reduction_map,
+            radial_branches.meshed_branches,
+        )
+        ax_ref = make_ax_ref(sort!(collect(radial_branches.meshed_branches)))
     else
-        radial_branches = RadialBranches()
+        if isempty(PTDFm.radial_branches)
+            radial_branches = RadialBranches()
+            A_matrix = A.data
+            ax_ref = make_ax_ref(A.axes[1])
+        else
+            error(
+                "Field `radial_branches` in PTDF must be empty if `reduce_radial_branches` is not true.",
+            )
+        end
     end
-    ax_ref = make_ax_ref(A.axes[1])
+
     if tol > eps()
-        lodf_t = _buildlodf(A.data, PTDFm.data, linear_solver)
+        lodf_t = _buildlodf(A_matrix, PTDFm_data, linear_solver)
         return LODF(
             sparsify(lodf_t, tol),
             (A.axes[1], A.axes[1]),
@@ -374,7 +408,7 @@ function LODF(
         )
     end
     return LODF(
-        _buildlodf(A.data, PTDFm.data, linear_solver),
+        _buildlodf(A_matrix, PTDFm_data, linear_solver),
         (A.axes[1], A.axes[1]),
         (ax_ref, ax_ref),
         Ref(tol),
@@ -413,13 +447,33 @@ function LODF(
     validate_linear_solver(linear_solver)
     ax_ref = make_ax_ref(A.axes[1])
     if reduce_radial_branches
-        radial_branches = RadialBranches(A)
+        # BA and ABA must contain the same, non-empty RadialBranches stucture
+        if !isempty(BA.radial_branches) && !isempty(ABA.radial_branches) &&
+           isequal(BA.radial_branches, ABA.radial_branches)
+            radial_branches = BA.radial_branches
+            @info "Non-empty `radial_branches` field found in BA and ABA matrix. LODF is evaluated considering radial branches and leaf nodes removed."
+        else
+            error("Mismatch in `radial_branches` field between BA and ABA matrices.")
+        end
+        A_matrix = reduce_A_matrix(
+            A,
+            radial_branches.bus_reduction_map,
+            radial_branches.meshed_branches,
+        )
     else
-        radial_branches = RadialBranches()
+        # BA and ABA must contain the same, empty RadialBranches stucture
+        if isempty(BA.radial_branches) && isempty(ABA.radial_branches)
+            radial_branches = RadialBranches()
+            A_matrix = A.data
+        else
+            error(
+                "Field `radial_branches` in BA and ABA must be empty if `reduce_radial_branches` is not true.",
+            )
+        end
     end
+
     if tol > eps()
-        lodf_t = _buildlodf(A.data, ABA.K, BA.data,
-            A.ref_bus_positions, linear_solver)
+        lodf_t = _buildlodf(A_matrix, ABA.K, BA.data, A.ref_bus_positions, linear_solver)
         return LODF(
             sparsify(lodf_t, tol),
             (A.axes[1], A.axes[1]),
@@ -429,7 +483,7 @@ function LODF(
         )
     end
     return LODF(
-        _buildlodf(A.data, ABA.K, BA.data, A.ref_bus_positions, linear_solver),
+        _buildlodf(A_matrix, ABA.K, BA.data, A.ref_bus_positions, linear_solver),
         (A.axes[1], A.axes[1]),
         (ax_ref, ax_ref),
         Ref(tol),
diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl
index 246e6208..3c50595a 100644
--- a/src/ptdf_calculations.jl
+++ b/src/ptdf_calculations.jl
@@ -21,9 +21,9 @@ The PTDF struct is indexed using the Bus numbers and Branch names.
 - `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
+- `radial_branches::RadialBranches`:
+        Structure containing the radial branches and leaf buses that were removed
+        while evaluating the matrix
 """
 struct PTDF{Ax, L <: NTuple{2, Dict}, M <: AbstractArray{Float64, 2}} <:
        PowerNetworkMatrix{Float64}
@@ -33,16 +33,7 @@ struct PTDF{Ax, L <: NTuple{2, Dict}, M <: AbstractArray{Float64, 2}} <:
     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)
+    radial_branches::RadialBranches
 end
 
 """
@@ -86,23 +77,24 @@ function _buildptdf(
 end
 
 function _buildptdf_from_matrices(
-    A::IncidenceMatrix,
+    A::SparseArrays.SparseMatrixCSC{Int8, Int},
     BA::SparseArrays.SparseMatrixCSC{T, Int} where {T <: Union{Float32, Float64}},
+    ref_bus_positions::Set{Int},
     dist_slack::Vector{Float64},
     linear_solver::String)
     if linear_solver == "KLU"
-        PTDFm = _calculate_PTDF_matrix_KLU(A.data, BA, A.ref_bus_positions, dist_slack)
+        PTDFm = _calculate_PTDF_matrix_KLU(A, BA, ref_bus_positions, dist_slack)
     elseif linear_solver == "Dense"
         # Convert SparseMatrices to Dense
         PTDFm = _calculate_PTDF_matrix_DENSE(
-            Matrix(A.data),
+            Matrix(A),
             Matrix(BA),
             A.ref_bus_positions,
             dist_slack,
         )
     elseif linear_solver == "MKLPardiso"
         PTDFm =
-            _calculate_PTDF_matrix_MKLPardiso(A.data, BA, A.ref_bus_positions, dist_slack)
+            _calculate_PTDF_matrix_MKLPardiso(A, BA, ref_bus_positions, dist_slack)
     end
 
     return PTDFm
@@ -376,6 +368,8 @@ Builds the PTDF matrix from a group of branches and buses. The return is a PTDF
         vector of the System AC branches
 - `buses::Vector{PSY.ACBus}`:
         vector of the System buses
+
+# Keyword Arguments
 - `dist_slack::Vector{Float64}`:
         vector of weights to be used as distributed slack bus.
         The distributed slack vector has to be the same length as the number of buses
@@ -383,9 +377,9 @@ Builds the PTDF matrix from a group of branches and buses. The return is a PTDF
         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
+- `radial_branches::RadialBranches`:
+        Structure containing the radial branches and leaf buses that were removed
+        while evaluating the ma
 """
 function PTDF(
     branches,
@@ -393,8 +387,10 @@ function PTDF(
     dist_slack::Vector{Float64} = Float64[],
     linear_solver::String = "KLU",
     tol::Float64 = eps(),
-    reduce_radial_branches::Bool = false)
+    radial_branches::RadialBranches = RadialBranches(),
+)
     validate_linear_solver(linear_solver)
+
     #Get axis names
     line_ax = [PSY.get_name(branch) for branch in branches]
     bus_ax = [PSY.get_number(bus) for bus in buses]
@@ -404,7 +400,7 @@ function PTDF(
     subnetworks = find_subnetworks(M, bus_ax)
     if length(subnetworks) > 1
         @info "Network is not connected, using subnetworks"
-        subnetworks = assing_reference_buses(subnetworks, ref_bus_positions)
+        subnetworks = assign_reference_buses(subnetworks, ref_bus_positions)
     end
     look_up = (bus_ax_ref, make_ax_ref(line_ax))
     S, _ = _buildptdf(
@@ -414,12 +410,6 @@ function PTDF(
         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),
@@ -436,18 +426,35 @@ end
 
 """
 Builds the PTDF matrix from a system. The return is a PTDF array indexed with the bus numbers.
+Note that `dist_slack` and `reduce_radial_branches` kwargs are explicitly mentioned because needed inside of the function.
 
 # Arguments
 - `sys::PSY.System`:
-        Power Systems system
+        PSY system for which the matrix is constructed
+
+# Keyword Arguments
+- `dist_slack::Vector{Float64}=Float64[]`:
+        vector of weights to be used as distributed slack bus.
+        The distributed slack vector has to be the same length as the number of buse
+- `reduce_radial_branches::Bool=false`:
+        if True the matrix will be evaluated discarding
+        all the radial branches and leaf buses (optional, default value is false)
 """
 function PTDF(
     sys::PSY.System;
+    dist_slack::Vector{Float64} = Float64[],
+    reduce_radial_branches::Bool = false,
     kwargs...,
 )
-    branches = get_ac_branches(sys)
-    buses = get_buses(sys)
-    return PTDF(branches, buses; kwargs...)
+    if reduce_radial_branches
+        A = IncidenceMatrix(sys)
+        dist_slack, rb = redistribute_dist_slack(dist_slack, A)
+    else
+        rb = RadialBranches()
+    end
+    branches = get_ac_branches(sys, rb.radial_branches)
+    buses = get_buses(sys, rb.bus_reduction_map)
+    return PTDF(branches, buses; dist_slack = dist_slack, radial_branches = rb, kwargs...)
 end
 
 """
@@ -458,6 +465,8 @@ Builds the PTDF matrix from a system. The return is a PTDF array indexed with th
         Incidence Matrix (full structure)
 - `BA::BA_Matrix`:
         BA matrix (full structure)
+
+# Keyword Arguments
 - `dist_slack::Vector{Float64}`:
         Vector of weights to be used as distributed slack bus.
         The distributed slack vector has to be the same length as the number of buses.
@@ -475,18 +484,43 @@ function PTDF(
     dist_slack::Vector{Float64} = Float64[],
     linear_solver = "KLU",
     tol::Float64 = eps(),
-    reduce_radial_branches::Bool = false)
+    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)
+        if !isempty(BA.radial_branches)
+            radial_branches = BA.radial_branches
+            @info "Non-empty `radial_branches` field found in BA matrix. PTDF is evaluated considering radial branches and leaf nodes removed."
+        else
+            error("BA has empty `radial_branches` field.")
+        end
+        A_matrix = reduce_A_matrix(
+            A,
+            radial_branches.bus_reduction_map,
+            radial_branches.meshed_branches,
+        )
+        axes = BA.axes
+        lookup = BA.lookup
     else
-        radial_branches = RadialBranches()
+        if isempty(BA.radial_branches)
+            radial_branches = RadialBranches()
+            A_matrix = A.data
+            axes = (A.axes[2], A.axes[1])
+            lookup = (A.lookup[2], A.lookup[1])
+        else
+            error(
+                "Field `radial_branches` in BA must be empty if `reduce_radial_branches` is not true.",
+            )
+        end
     end
+    S = _buildptdf_from_matrices(
+        A_matrix,
+        BA.data,
+        BA.ref_bus_positions,
+        dist_slack,
+        linear_solver,
+    )
     if tol > eps()
         return PTDF(
             sparsify(S, tol),
@@ -510,6 +544,10 @@ function PTDF(
     end
 end
 
+##############################################################################
+########################### Auxiliary functions ##############################
+##############################################################################
+
 # PTDF stores the transposed matrix. Overload indexing and how data is exported.
 function Base.getindex(A::PTDF, line, bus)
     i, j = to_index(A, bus, line)
@@ -539,3 +577,26 @@ end
 function get_tol(ptdf::PTDF)
     return ptdf.tol
 end
+
+function redistribute_dist_slack(
+    dist_slack::Vector{Float64},
+    A::IncidenceMatrix,
+)
+    dist_slack1 = deepcopy(dist_slack)
+    rb = RadialBranches(A)
+    # if original length of dist_slack is correct
+    if length(dist_slack) == size(A.data, 2)
+        for i in keys(rb.bus_reduction_map)
+            for j in rb.bus_reduction_map[i]
+                dist_slack1[A.lookup[2][i]] += dist_slack1[A.lookup[2][j]]
+                dist_slack1[A.lookup[2][j]] = -9999
+            end
+        end
+        # redefine dist_slack
+        return dist_slack1[dist_slack1 .!= -9999], rb
+        # otherwise throw an error
+    elseif !isempty(dist_slack) && length(dist_slack) != size(A.data, 2)
+        error("Distributed bus specification doesn't match the number of the buses.")
+    end
+    return dist_slack, rb
+end
diff --git a/src/radial_braches.jl b/src/radial_braches.jl
index 75b60d67..91ef19f4 100644
--- a/src/radial_braches.jl
+++ b/src/radial_braches.jl
@@ -1,17 +1,65 @@
 struct RadialBranches
     bus_reduction_map::Dict{Int, Set{Int}}
+    reverse_bus_search_map::Dict{Int, Int}
     radial_branches::Set{String}
+    meshed_branches::Set{String}
 end
 
 get_bus_reduction_map(rb::RadialBranches) = rb.bus_reduction_map
+get_reverse_bus_search_map(rb::RadialBranches) = rb.reverse_bus_search_map
 get_radial_branches(rb::RadialBranches) = rb.radial_branches
-Base.isempty(rb::RadialBranches) =
-    isempty(rb.bus_reduction_map) && isempty(rb.bus_reduction_map)
+get_meshed_branches(rb::RadialBranches) = rb.meshed_branches
+
+function Base.isempty(rb::RadialBranches)
+    if !isempty(rb.bus_reduction_map)
+        return false
+    end
+
+    if !isempty(rb.bus_reduction_map)
+        return false
+    end
+
+    if !isempty(rb.reverse_bus_search_map)
+        return false
+    end
+    return true
+end
 
 function RadialBranches(;
     bus_reduction_map::Dict{Int, Set{Int}} = Dict{Int, Set{Int}}(),
-    radial_branches::Set{String} = Set{String}())
-    return RadialBranches(bus_reduction_map, radial_branches)
+    reverse_bus_search_map::Dict{Int, Int} = Dict{Int, Int}(),
+    radial_branches::Set{String} = Set{String}(),
+    meshed_branches::Set{String} = Set{String}())
+    return RadialBranches(
+        bus_reduction_map,
+        reverse_bus_search_map,
+        radial_branches,
+        meshed_branches,
+    )
+end
+
+"""
+Structure to collect information the branches in the system that are radial and can be
+ignored in the models.
+
+# Arguments
+- `A::IncidenceMatrix`: IncidenceMatrix
+
+"""
+function RadialBranches(A::IncidenceMatrix)
+    return calculate_radial_branches(A.data, A.lookup[1], A.lookup[2], A.ref_bus_positions)
+end
+
+"""
+Structure to collect information the branches in the system that are radial and can be
+ignored in the models.
+
+# Arguments
+- `sys::PSY.System`: System Data
+"""
+
+function RadialBranches(sys::PSY.System)
+    return RadialBranches(IncidenceMatrix(sys))
 end
 
 function _find_upstream_bus(
@@ -52,7 +100,7 @@ function _new_parent(
         )
         new_parent_bus_number = reverse_bus_map[new_parent_val]
         new_set = push!(pop!(bus_reduction_map_index, parent_bus_number), parent_bus_number)
-        bus_reduction_map_index[new_parent_bus_number] = new_set
+        union!(bus_reduction_map_index[new_parent_bus_number], new_set)
         _new_parent(
             A,
             new_parent_val,
@@ -74,6 +122,7 @@ function _reverse_search(
     reverse_bus_map::Dict{Int, Int},
 )
     j_bus_number = reverse_bus_map[j]
+    pop!(bus_reduction_map_index, j_bus_number)
     reducion_set = Set{Int}(j_bus_number)
     parent = _find_upstream_bus(
         A,
@@ -84,7 +133,7 @@ function _reverse_search(
         reverse_bus_map,
     )
     parent_bus_number = reverse_bus_map[parent]
-    bus_reduction_map_index[parent_bus_number] = reducion_set
+    union!(bus_reduction_map_index[parent_bus_number], reducion_set)
     _new_parent(
         A,
         parent,
@@ -96,11 +145,30 @@ function _reverse_search(
     return
 end
 
-function RadialBranches(A::IncidenceMatrix)
-    return RadialBranches(A.data, A.lookup[1], A.lookup[2], A.ref_bus_positions)
+function _make_reverse_bus_search_map(bus_reduction_map::Dict{Int, Set{Int}}, n_buses::Int)
+    map = Dict{Int, Int}()
+    sizehint!(map, n_buses)
+    for (parent, children) in bus_reduction_map
+        for bus in children
+            map[bus] = parent
+        end
+    end
+    return map
 end
 
-function RadialBranches(
+"""
+used to calculate the branches in the system that are radial and can be
+ignored in the models by exploring the structure of the incidence matrix
+
+# Arguments
+- `A::SparseArrays.SparseMatrixCSC{Int8, Int64}`: Data from the IncidenceMatrix
+- `line_map::Dict{String, Int}`: Map of Line Name to Matrix Index
+- `bus_map::Dict{Int, Int}`: Map of Bus Name to Matrix Index
+- `ref_bus_positions::Set{Int}`:
+        Set containing the indexes of the columns of the BA matrix corresponding
+        to the refence buses
+"""
+function calculate_radial_branches(
     A::SparseArrays.SparseMatrixCSC{Int8, Int64},
     line_map::Dict{String, Int},
     bus_map::Dict{Int, Int},
@@ -108,15 +176,15 @@ function RadialBranches(
 )
     lk = ReentrantLock()
     buscount = length(bus_map)
-    bus_reduction_map_index = Dict{Int, Set{Int}}()
     radial_branches = Set{String}()
     reverse_line_map = Dict(reverse(kv) for kv in line_map)
     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
-            if j ∈ ref_bus_positions
-                continue
-            end
             lock(lk) do
                 _reverse_search(
                     A,
@@ -129,6 +197,48 @@ function RadialBranches(
             end
         end
     end
+    reverse_bus_search_map = _make_reverse_bus_search_map(bus_reduction_map_index, buscount)
+    meshed_branches = Set{String}()
+    for k in keys(line_map)
+        if k in radial_branches
+            continue
+        end
+        push!(meshed_branches, k)
+    end
+    return RadialBranches(
+        bus_reduction_map_index,
+        reverse_bus_search_map,
+        radial_branches,
+        meshed_branches,
+    )
+end
+
+"""
+Interface to obtain the parent bus number of a reduced bus when radial branches are eliminated
+
+# Arguments
+- `rb::RadialBranches`: RadialBranches object
+- `bus_number::Int`: Bus number of the reduced bus
+"""
+function get_mapped_bus_number(rb::RadialBranches, bus_number::Int)
+    if isempty(rb)
+        return bus_number
+    end
+    return get(rb.reverse_bus_search_map, bus_number, bus_number)
+end
 
-    return RadialBranches(bus_reduction_map_index, radial_branches)
+##############################################################################
+########################### Auxiliary functions ##############################
+##############################################################################
+
+function isequal(
+    rb1::RadialBranches,
+    rb2::RadialBranches,
+)
+    for field in fieldnames(typeof(rb1))
+        if getfield(rb1, field) != getfield(rb2, field)
+            return false
+        end
+    end
+    return true
 end
diff --git a/src/system_utils.jl b/src/system_utils.jl
index af322dc8..1a3df072 100644
--- a/src/system_utils.jl
+++ b/src/system_utils.jl
@@ -21,5 +21,5 @@ by the reference bus of the subnetworks if they exist
 """
 function find_subnetworks(sys::PSY.System)
     sbn, ref_bus_positions = _find_subnetworks(sys)
-    return assing_reference_buses(sbn, ref_bus_positions)
+    return assign_reference_buses(sbn, ref_bus_positions)
 end
diff --git a/src/virtual_lodf_calculations.jl b/src/virtual_lodf_calculations.jl
index 1659d65e..be9c1662 100644
--- a/src/virtual_lodf_calculations.jl
+++ b/src/virtual_lodf_calculations.jl
@@ -38,6 +38,9 @@ The VirtualLODF struct is indexed using branch names.
         Dictionary containing the subsets of buses defining the different subnetwork of the system.
 - `tol::Base.RefValue{Float64}`:
         Tolerance related to scarification and values to drop.
+- `radial_branches::RadialBranches`:
+        Structure containing the radial branches and leaf buses that were removed
+        while evaluating the matrix
 """
 struct VirtualLODF{Ax, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{Float64}
     K::KLU.KLUFactorization{Float64, Int}
@@ -52,7 +55,7 @@ struct VirtualLODF{Ax, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{Float64}
     cache::RowCache
     subnetworks::Dict{Int, Set{Int}}
     tol::Base.RefValue{Float64}
-    radial_banches::RadialBranches
+    radial_branches::RadialBranches
 end
 
 function Base.show(io::IO, ::MIME{Symbol("text/plain")}, array::VirtualLODF)
@@ -91,22 +94,64 @@ function _get_PTDF_A_diag(
     return diag_
 end
 
+"""
+Builds the Virtual LODF matrix from a system. The return is a VirtualLODF
+struct with an empty cache.
+
+# Arguments
+- `sys::PSY.System`:
+        PSY system for which the matrix is constructed
+
+# Keyword Arguments
+- `reduce_radial_branches::Bool=false`:
+        if True the matrix will be evaluated discarding
+        all the radial branches and leaf buses (optional, default value is false)
+- `kwargs...`:
+        other keyword arguments used by VirtualPTDF
+"""
 function VirtualLODF(
-    sys::PSY.System; kwargs...,
+    sys::PSY.System;
+    reduce_radial_branches::Bool = false,
+    kwargs...,
 )
-    branches = get_ac_branches(sys)
-    buses = get_buses(sys)
-
-    return VirtualLODF(branches, buses; kwargs...)
+    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)
+    return VirtualLODF(branches, buses; radial_branches = rb, kwargs...)
 end
 
+"""
+Builds the LODF matrix from a group of branches and buses. The return is a
+VirtualLODF struct with an empty cache.
+
+# Arguments
+- `branches`:
+        Vector of the system's AC branches.
+- `buses::Vector{PSY.ACBus}`:
+        Vector of the system's buses.
+
+# Keyword Arguments
+- `tol::Float64 = eps()`:
+        Tolerance related to sparsification and values to drop.
+- `max_cache_size::Int`:
+        max cache size in MiB (inizialized as MAX_CACHE_SIZE_MiB).
+- `persistent_lines::Vector{String}`:
+        line to be evaluated as soon as the VirtualPTDF is created (initialized as empty vector of strings).
+- `radial_branches::RadialBranches`:
+        Structure containing the radial branches and leaf buses that were removed
+        while evaluating the matrix
+"""
 function VirtualLODF(
     branches,
     buses::Vector{PSY.ACBus};
     tol::Float64 = eps(),
     max_cache_size::Int = MAX_CACHE_SIZE_MiB,
     persistent_lines::Vector{String} = String[],
-    reduce_radial_branches::Bool = false,
+    radial_branches::RadialBranches = RadialBranches(),
 )
 
     #Get axis names and lookups
@@ -126,7 +171,7 @@ function VirtualLODF(
     # check subnetworks
     if length(subnetworks) > 1
         @info "Network is not connected, using subnetworks"
-        subnetworks = assing_reference_buses(subnetworks, ref_bus_positions)
+        subnetworks = assign_reference_buses(subnetworks, ref_bus_positions)
     end
     # get diagonal of PTDF
     temp_data = zeros(length(bus_ax))
@@ -151,13 +196,6 @@ function VirtualLODF(
                 length(bus_ax) * sizeof(Float64),
             )
     end
-    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
-
     return VirtualLODF(
         K,
         BA,
diff --git a/src/virtual_ptdf_calculations.jl b/src/virtual_ptdf_calculations.jl
index ddf27f35..b3137070 100644
--- a/src/virtual_ptdf_calculations.jl
+++ b/src/virtual_ptdf_calculations.jl
@@ -42,6 +42,9 @@ matrix.
         Dictionary containing the subsets of buses defining the different subnetwork of the system.
 - `tol::Base.RefValue{Float64}`:
         Tolerance related to scarification and values to drop.
+- `radial_branches::RadialBranches`:
+        Structure containing the radial branches and leaf buses that were removed
+        while evaluating the matrix
 """
 struct VirtualPTDF{Ax, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{Float64}
     K::KLU.KLUFactorization{Float64, Int}
@@ -55,7 +58,7 @@ struct VirtualPTDF{Ax, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{Float64}
     cache::RowCache
     subnetworks::Dict{Int, Set{Int}}
     tol::Base.RefValue{Float64}
-    radial_banches::RadialBranches
+    radial_branches::RadialBranches
 end
 
 function Base.show(io::IO, ::MIME{Symbol("text/plain")}, array::VirtualPTDF)
@@ -68,13 +71,15 @@ end
 
 """
 Builds the PTDF matrix from a group of branches and buses. The return is a
-PTDF array indexed with the branch numbers.
+VirtualPTDF struct with an empty cache.
 
 # Arguments
 - `branches`:
         Vector of the system's AC branches.
 - `buses::Vector{PSY.ACBus}`:
         Vector of the system's buses.
+
+# Keyword Arguments
 - `dist_slack::Vector{Float64} = Float64[]`:
         Vector of weights to be used as distributed slack bus.
         The distributed slack vector has to be the same length as the number of buses.
@@ -84,9 +89,9 @@ PTDF array indexed with the branch numbers.
         max cache size in MiB (inizialized as MAX_CACHE_SIZE_MiB).
 - `persistent_lines::Vector{String}`:
         line to be evaluated as soon as the VirtualPTDF is created (initialized as empty vector of strings).
-- `reduce_radial_branches::Bool`:
-        True to reduce the network by simplifying the radial branches and mapping the
-        eliminate buses
+- `radial_branches::RadialBranches`:
+        Structure containing the radial branches and leaf buses that were removed
+        while evaluating the matrix
 """
 function VirtualPTDF(
     branches,
@@ -95,7 +100,8 @@ function VirtualPTDF(
     tol::Float64 = eps(),
     max_cache_size::Int = MAX_CACHE_SIZE_MiB,
     persistent_lines::Vector{String} = String[],
-    reduce_radial_branches::Bool = false)
+    radial_branches::RadialBranches = RadialBranches(),
+)
     if length(dist_slack) != 0
         @info "Distributed bus"
     end
@@ -114,7 +120,7 @@ function VirtualPTDF(
     subnetworks = find_subnetworks(M, bus_ax)
     if length(subnetworks) > 1
         @info "Network is not connected, using subnetworks"
-        subnetworks = assing_reference_buses(subnetworks, ref_bus_positions)
+        subnetworks = assign_reference_buses(subnetworks, ref_bus_positions)
     end
     temp_data = zeros(length(bus_ax))
 
@@ -131,13 +137,6 @@ function VirtualPTDF(
             )
     end
 
-    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
-
     return VirtualPTDF(
         klu(ABA),
         BA,
@@ -157,13 +156,42 @@ end
 """
 Builds the Virtual PTDF matrix from a system. The return is a VirtualPTDF
 struct with an empty cache.
+
+# Arguments
+- `sys::PSY.System`:
+        PSY system for which the matrix is constructed
+
+# Keyword Arguments
+- `dist_slack::Vector{Float64}=Float64[]`:
+        vector of weights to be used as distributed slack bus.
+        The distributed slack vector has to be the same length as the number of buse
+- `reduce_radial_branches::Bool=false`:
+        if True the matrix will be evaluated discarding
+        all the radial branches and leaf buses (optional, default value is false)
+- `kwargs...`:
+        other keyword arguments used by VirtualPTDF
 """
 function VirtualPTDF(
-    sys::PSY.System; kwargs...)
-    branches = get_ac_branches(sys)
-    buses = get_buses(sys)
-
-    return VirtualPTDF(branches, buses; kwargs...)
+    sys::PSY.System;
+    dist_slack::Vector{Float64} = Float64[],
+    reduce_radial_branches::Bool = false,
+    kwargs...,
+)
+    if reduce_radial_branches
+        A = IncidenceMatrix(sys)
+        dist_slack, rb = redistribute_dist_slack(dist_slack, A)
+    else
+        rb = RadialBranches()
+    end
+    branches = get_ac_branches(sys, rb.radial_branches)
+    buses = get_buses(sys, rb.bus_reduction_map)
+    return VirtualPTDF(
+        branches,
+        buses;
+        dist_slack = dist_slack,
+        radial_branches = rb,
+        kwargs...,
+    )
 end
 
 # Overload Base functions
diff --git a/test/test_A_BA_ABA_with_radial_lines.jl b/test/test_A_BA_ABA_with_radial_lines.jl
new file mode 100644
index 00000000..c82a6c22
--- /dev/null
+++ b/test/test_A_BA_ABA_with_radial_lines.jl
@@ -0,0 +1,117 @@
+@testset "Test BA matrix with radial lines" begin
+    for name in ["c_sys14", "test_RTS_GMLC_sys"]
+        # load the system
+        sys = PSB.build_system(PSB.PSITestSystems, name)
+        # get the incidence matrix
+        BA = BA_Matrix(sys)
+        # ... and with radial lines
+        BA_rad = BA_Matrix(sys; reduce_radial_branches = true)
+        # get inidices for the leaf nodes
+        rb = BA_rad.radial_branches
+        bus_numbers = []
+        for i in keys(rb.bus_reduction_map)
+            append!(bus_numbers, collect(rb.bus_reduction_map[i]))
+        end
+        bus_idx = setdiff(1:size(BA.data, 1), [BA.lookup[1][i] for i in bus_numbers])
+        # ... and radial branches
+        br_idx = setdiff(1:size(BA.data, 2), [BA.lookup[2][i] for i in rb.radial_branches])
+        # now extract A matrix anc compare
+        @test all(isapprox.(BA.data[bus_idx, br_idx], BA_rad.data))
+    end
+end
+
+@testset "Test ABA matrix with radial lines" begin
+    # to test the if the ABA matrix contains the same information, the power
+    # flows must be evaluated
+    for name in ["c_sys14", "test_RTS_GMLC_sys"]
+        # load the system
+        sys = PSB.build_system(PSB.PSITestSystems, name)
+
+        # get the RadialBranches struct
+        rb = RadialBranches(IncidenceMatrix(sys))
+
+        # get the original and reduced IncidenceMatrix, BA and ABA
+        A = IncidenceMatrix(sys)
+        A_rad = IncidenceMatrix(sys)
+        BA = BA_Matrix(sys)
+        BA_rad = BA_Matrix(sys; reduce_radial_branches = true)
+        ABA = ABA_Matrix(sys; factorize = true)
+        ABA_rad = ABA_Matrix(sys; factorize = true, reduce_radial_branches = true)
+
+        # check if the same angles and flows are coputed with the matrices of the reduced systems
+        # get the indices for the reduced system
+        bus_numbers = []
+        for i in keys(rb.bus_reduction_map)
+            append!(bus_numbers, collect(rb.bus_reduction_map[i]))
+        end
+        bus_idx = setdiff(
+            1:size(A.data, 2),
+            append!([A.lookup[2][i] for i in bus_numbers], A.ref_bus_positions),
+        )
+        br_idx = setdiff(1:size(A.data, 1), [A.lookup[1][i] for i in rb.radial_branches])
+
+        # now get the injuctions from the system
+        n_buses = length(axes(BA, 1))
+        bus_lookup = BA.lookup[1]
+        branch_lookup = BA.lookup[2]
+        bus_angles = zeros(Float64, n_buses)
+        branch_flow_values = zeros(Float64, length(axes(BA, 2)))
+        temp_bus_map = Dict{Int, String}(
+            PSY.get_number(b) => PSY.get_name(b) for
+            b in PSY.get_components(PSY.ACBus, sys)
+        )
+        for (bus_no, ix) in bus_lookup
+            bus_name = temp_bus_map[bus_no]
+            bus = PSY.get_component(PSY.Bus, sys, bus_name)
+            if PSY.get_bustype(bus) == PSY.ACBusTypes.REF
+                bus_angles[ix] = 0.0
+            else
+                bus_angles[ix] = PSY.get_angle(bus)
+            end
+        end
+        bus_activepower_injection = zeros(Float64, n_buses)
+        sources =
+            PSY.get_components(d -> !isa(d, PSY.ElectricLoad), PSY.StaticInjection, sys)
+        for source in sources
+            !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)
+        end
+        bus_activepower_withdrawals = zeros(Float64, n_buses)
+        loads = PSY.get_components(x -> !isa(x, PSY.FixedAdmittance), PSY.ElectricLoad, sys)
+        for l in loads
+            !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)
+        end
+        power_injection =
+            deepcopy(bus_activepower_injection - bus_activepower_withdrawals)
+        valid_ix = setdiff(1:length(power_injection), BA.ref_bus_positions)
+        ref_bus_angles = deepcopy(bus_angles)
+        ref_flow_values = deepcopy(branch_flow_values)
+
+        # now get the reference power flows and angle evaluated with BA and ABA
+        ref_bus_angles[valid_ix] = ABA.K \ power_injection[valid_ix]
+        ref_flow_values = transpose(BA.data) * ref_bus_angles
+
+        # evalaute according to the matrices with no radial branches
+        reduced_bus_angles = zeros((length(bus_idx) + length(A.ref_bus_positions),))
+        reduce_flow_values = zeros((length(br_idx),))
+        # change power injection for affrefated leaf buses
+        power_injection2 = deepcopy(power_injection)
+        for i in keys(rb.bus_reduction_map)
+            for j in rb.bus_reduction_map[i]
+                power_injection2[BA.lookup[1][i]] += power_injection[BA.lookup[1][j]]
+            end
+        end
+        valid_ix2 = setdiff(1:size(BA_rad.data, 1), BA_rad.ref_bus_positions)
+        reduced_bus_angles[valid_ix2] = ABA_rad.K \ power_injection2[bus_idx]
+        reduced_flow_values = transpose(BA_rad.data) * reduced_bus_angles
+
+        # now check if flows and angles are the same
+        @test isapprox(ref_bus_angles[bus_idx], reduced_bus_angles[valid_ix2])
+        @test isapprox(ref_flow_values[br_idx], reduced_flow_values)
+    end
+end
diff --git a/test/test_lodf.jl b/test/test_lodf.jl
index 0dec2892..3fb8ca9f 100644
--- a/test/test_lodf.jl
+++ b/test/test_lodf.jl
@@ -57,8 +57,10 @@
 
     # test if error is thrown in case `tol` is defined in PTDF
     P5 = PTDF(sys5; tol = 1e-3)
-    test_value = false
-    @test_throws ErrorException L5NS_from_ptdf = LODF(A, P5)
+    @test_logs (
+        :warn,
+        "The argument `tol` in the PTDF matrix was set to a value different than the default one.\nThe resulting LODF can include unexpected rounding errors.\n",
+    ) match_mode = :any LODF(A, P5)
 
     # get 14 bus system
     buses_14 = nodes14()
diff --git a/test/test_lodf_with_radial_branches.jl b/test/test_lodf_with_radial_branches.jl
new file mode 100644
index 00000000..948d7ecb
--- /dev/null
+++ b/test/test_lodf_with_radial_branches.jl
@@ -0,0 +1,138 @@
+@testset "Test LODF with radial branches" begin
+    # check if the flows obtained with original LODF and reduced one are the same
+    for name in ["c_sys14", "test_RTS_GMLC_sys"]
+        # load the system
+        sys = PSB.build_system(PSB.PSITestSystems, name)
+
+        # get A, BA and ABA matrix with radial branches
+        A_rad = IncidenceMatrix(sys)
+        BA_rad = BA_Matrix(sys, ; reduce_radial_branches = true)
+        ABA_rad = ABA_Matrix(sys; reduce_radial_branches = true, factorize = true)
+        ptdf = PTDF(sys)
+        ptdf_rad = PTDF(sys; reduce_radial_branches = true)
+
+        # get the original and reduced PTDF matrices (consider different methods)
+        lodf = LODF(sys)
+        lodf_rad = LODF(sys; reduce_radial_branches = true)
+        lodf_rad_A_BA_ABA = LODF(A_rad, ABA_rad, BA_rad; reduce_radial_branches = true)
+        lodf_rad_A_PTDF = LODF(A_rad, ptdf_rad; reduce_radial_branches = true)
+
+        rb = RadialBranches(IncidenceMatrix(sys))
+
+        # at first check if all the matrices are the same
+        @test isapprox(lodf_rad.data, lodf_rad_A_BA_ABA.data, atol = 1e-10)
+        @test isapprox(lodf_rad_A_PTDF.data, lodf_rad_A_BA_ABA.data, atol = 1e-10)
+
+        # now check if the flows are the same
+
+        # first evaluate the flows on the original and reduced system
+        bus_numbers = []
+        for i in keys(rb.bus_reduction_map)
+            append!(bus_numbers, collect(rb.bus_reduction_map[i]))
+        end
+        bus_idx = setdiff(
+            1:size(ptdf.data, 1),
+            append!([ptdf.lookup[1][i] for i in bus_numbers]),
+        )
+        br_idx =
+            setdiff(1:size(ptdf.data, 2), [ptdf.lookup[2][i] for i in rb.radial_branches])
+        # now get the injections from the system
+        n_buses = length(axes(ptdf, 1))
+        bus_lookup = ptdf.lookup[1]
+        branch_flow_values = zeros(Float64, length(axes(ptdf, 2)))
+        bus_activepower_injection = zeros(Float64, n_buses)
+        sources =
+            PSY.get_components(d -> !isa(d, PSY.ElectricLoad), PSY.StaticInjection, sys)
+        for source in sources
+            !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)
+        end
+        bus_activepower_withdrawals = zeros(Float64, n_buses)
+        loads = PSY.get_components(x -> !isa(x, PSY.FixedAdmittance), PSY.ElectricLoad, sys)
+        for l in loads
+            !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)
+        end
+        power_injection =
+            deepcopy(bus_activepower_injection - bus_activepower_withdrawals)
+        # get the flows with the PTDF matrix
+        ref_flow_values = transpose(ptdf.data) * power_injection
+        # evalaute according to the matrix with no radial branches
+        reduce_flow_values = zeros((length(br_idx),))
+        # change power injection for affrefated leaf buses
+        power_injection2 = deepcopy(power_injection)
+        for i in keys(rb.bus_reduction_map)
+            for j in rb.bus_reduction_map[i]
+                power_injection2[ptdf.lookup[1][i]] += power_injection[ptdf.lookup[1][j]]
+            end
+        end
+        reduced_flow_values = transpose(ptdf_rad.data) * power_injection2[bus_idx]
+
+        # now evaluate the flows with the reduced LODF matrix
+        ref_lodf_flows = transpose(lodf.data) * ref_flow_values
+        ref_lodf_flows_reduced1 = transpose(lodf_rad.data) * reduced_flow_values
+        ref_lodf_flows_reduced2 = transpose(lodf_rad.data) * ref_flow_values[br_idx]
+
+        @test isapprox(ref_lodf_flows_reduced1, ref_lodf_flows_reduced2, atol = 1e-10)
+        @test isapprox(ref_lodf_flows[br_idx], ref_lodf_flows_reduced1, atol = 1e-10)
+    end
+end
+
+@testset "Test LODF errors" begin
+
+    # load the system
+    sys = PSB.build_system(PSB.PSITestSystems, "c_sys14")
+
+    # get A, BA and ABA matrix with radial branches
+    A = IncidenceMatrix(sys)
+    BA_rad = BA_Matrix(sys, ; reduce_radial_branches = true)
+    ABA = ABA_Matrix(sys; factorize = true)
+    ptdf = PTDF(sys)
+    ptdf_rad = PTDF(sys; reduce_radial_branches = true)
+
+    # test LODF from A, ABA and BA
+    test_value = false
+    try
+        lodf_rad_A_BA_ABA = LODF(A, ABA, BA_rad; reduce_radial_branches = true)
+    catch err
+        if err isa Exception
+            test_value = true
+        end
+    end
+    @test test_value
+
+    test_value = false
+    try
+        lodf_rad_A_BA_ABA = LODF(A, ABA, BA_rad)
+    catch err
+        if err isa Exception
+            test_value = true
+        end
+    end
+    @test test_value
+
+    # test LODF from A, PTDF
+    test_value = false
+    try
+        lodf_rad_A_PTDF = LODF(A_rad, ptdf; reduce_radial_branches = true)
+    catch err
+        if err isa Exception
+            test_value = true
+        end
+    end
+    @test test_value
+
+    test_value = false
+    try
+        lodf_rad_A_PTDF = LODF(A_rad, ptdf_rad)
+    catch err
+        if err isa Exception
+            test_value = true
+        end
+    end
+    @test test_value
+end
diff --git a/test/test_ptdf_with_radial_lines.jl b/test/test_ptdf_with_radial_lines.jl
new file mode 100644
index 00000000..ddb154ed
--- /dev/null
+++ b/test/test_ptdf_with_radial_lines.jl
@@ -0,0 +1,183 @@
+@testset "Test PTDF matrix with radial lines" begin
+    for name in ["c_sys14", "test_RTS_GMLC_sys"]
+        # load the system
+        sys = PSB.build_system(PSB.PSITestSystems, name)
+
+        # get the RadialBranches struct
+        rb = RadialBranches(IncidenceMatrix(sys))
+
+        # get the A and BA matrices without radial lines
+        A_rad = IncidenceMatrix(sys)
+        BA_rad = BA_Matrix(sys; reduce_radial_branches = true)
+
+        # get the original and reduced PTDF matrices (consider different methods)
+        ptdf = PTDF(sys)
+        ptdf_rad = PTDF(sys; reduce_radial_branches = true)
+        ptdf_rad_A_BA = PTDF(A_rad, BA_rad; reduce_radial_branches = true)
+
+        # check if the same angles and flows are coputed with the matrices of the reduced systems
+        # get the indices for the reduced system
+        bus_numbers = []
+        for i in keys(rb.bus_reduction_map)
+            append!(bus_numbers, collect(rb.bus_reduction_map[i]))
+        end
+        bus_idx = setdiff(
+            1:size(ptdf.data, 1),
+            append!([ptdf.lookup[1][i] for i in bus_numbers]),
+        )
+        br_idx =
+            setdiff(1:size(ptdf.data, 2), [ptdf.lookup[2][i] for i in rb.radial_branches])
+
+        # now get the injections from the system
+        n_buses = length(axes(ptdf, 1))
+        bus_lookup = ptdf.lookup[1]
+        branch_flow_values = zeros(Float64, length(axes(ptdf, 2)))
+        bus_activepower_injection = zeros(Float64, n_buses)
+        sources =
+            PSY.get_components(d -> !isa(d, PSY.ElectricLoad), PSY.StaticInjection, sys)
+        for source in sources
+            !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)
+        end
+        bus_activepower_withdrawals = zeros(Float64, n_buses)
+        loads = PSY.get_components(x -> !isa(x, PSY.FixedAdmittance), PSY.ElectricLoad, sys)
+        for l in loads
+            !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)
+        end
+        power_injection =
+            deepcopy(bus_activepower_injection - bus_activepower_withdrawals)
+
+        # get the flows with the PTDF matrix
+        ref_flow_values = transpose(ptdf.data) * power_injection
+
+        # evalaute according to the matrix with no radial branches
+        reduce_flow_values = zeros((length(br_idx),))
+        # change power injection for affrefated leaf buses
+        power_injection2 = deepcopy(power_injection)
+        for i in keys(rb.bus_reduction_map)
+            for j in rb.bus_reduction_map[i]
+                power_injection2[ptdf.lookup[1][i]] += power_injection[ptdf.lookup[1][j]]
+            end
+        end
+        reduced_flow_values = transpose(ptdf_rad.data) * power_injection2[bus_idx]
+
+        # now check if flows are the same
+        @test isapprox(ref_flow_values[br_idx], reduced_flow_values)
+        # for the PTDF from A and BA matrices just need to check the elements
+        @test isapprox(ptdf_rad.data, ptdf_rad_A_BA.data, atol = 1e-5)
+    end
+end
+
+@testset "Test PTDF with radial lines and distributed slack" begin
+    for name in ["c_sys14", "test_RTS_GMLC_sys"]
+        # load the system
+        sys = PSB.build_system(PSB.PSITestSystems, name)
+        # get the radial branches
+        A = IncidenceMatrix(sys)
+        rb = RadialBranches(A)
+        # get number of buses
+        buscount = length(PNM.get_buses(sys))
+        # now compute distributed slack vector
+        dist_slack = 1 / buscount * ones(buscount)
+        slack_array = dist_slack / sum(dist_slack)
+        # adjust to have the same vector with and without leaf node reduction
+        bus_numbers = reduce(
+            vcat,
+            [collect(rb.bus_reduction_map[i]) for i in keys(rb.bus_reduction_map)],
+        )
+        bus_idx = setdiff(
+            1:size(A.data, 2),
+            append!([A.lookup[2][i] for i in bus_numbers]),
+        )
+        br_idx = setdiff(1:size(A.data, 1), [A.lookup[1][i] for i in rb.radial_branches])
+        for i in keys(rb.bus_reduction_map)
+            for j in rb.bus_reduction_map[i]
+                slack_array[A.lookup[2][i]] += slack_array[A.lookup[2][j]]
+                slack_array[A.lookup[2][j]] = -9999
+            end
+        end
+        # redefine slack_array
+        slack_array[slack_array .== -9999] .= 0
+        slack_array1 = slack_array[slack_array .!= -9999]
+        # now get the PTDF matrices
+        ptdf = PTDF(sys; dist_slack = slack_array)
+        ptdf_rad = PTDF(sys; reduce_radial_branches = true, dist_slack = slack_array1)
+
+        # now get the injections from the system
+        n_buses = length(axes(ptdf, 1))
+        bus_lookup = ptdf.lookup[1]
+        branch_flow_values = zeros(Float64, length(axes(ptdf, 2)))
+        bus_activepower_injection = zeros(Float64, n_buses)
+        sources =
+            PSY.get_components(d -> !isa(d, PSY.ElectricLoad), PSY.StaticInjection, sys)
+        for source in sources
+            !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)
+        end
+        bus_activepower_withdrawals = zeros(Float64, n_buses)
+        loads = PSY.get_components(x -> !isa(x, PSY.FixedAdmittance), PSY.ElectricLoad, sys)
+        for l in loads
+            !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)
+        end
+        power_injection =
+            deepcopy(bus_activepower_injection - bus_activepower_withdrawals)
+
+        # get the flows with the PTDF matrix
+        ref_flow_values = transpose(ptdf.data) * power_injection
+
+        # evalaute according to the matrix with no radial branches
+        reduce_flow_values = zeros((length(br_idx),))
+        # change power injection for affrefated leaf buses
+        power_injection2 = deepcopy(power_injection)
+        for i in keys(rb.bus_reduction_map)
+            for j in rb.bus_reduction_map[i]
+                power_injection2[ptdf.lookup[1][i]] += power_injection[ptdf.lookup[1][j]]
+            end
+        end
+        reduced_flow_values = transpose(ptdf_rad.data) * power_injection2[bus_idx]
+
+        # now check if flows are the same
+        @test isapprox(ref_flow_values[br_idx], reduced_flow_values)
+    end
+end
+
+@testset "Test PTDF errors" begin
+
+    # load the system
+    sys = PSB.build_system(PSB.PSITestSystems, "c_sys14")
+
+    # get the A and BA matrices without radial lines
+    A = IncidenceMatrix(sys)
+    BA = BA_Matrix(sys)
+    BA_rad = BA_Matrix(sys; reduce_radial_branches = true)
+
+    test_value = false
+    try
+        ptdf = PTDF(A, BA; reduce_radial_branches = true)
+    catch err
+        if err isa Exception
+            test_value = true
+        end
+    end
+    @test test_value
+
+    test_value = false
+    try
+        ptdf = PTDF(A, BA_rad)
+    catch err
+        if err isa Exception
+            test_value = true
+        end
+    end
+    @test test_value
+end
diff --git a/test/test_radial_branches.jl b/test/test_radial_branches.jl
index abdc96c1..76519152 100644
--- a/test/test_radial_branches.jl
+++ b/test/test_radial_branches.jl
@@ -37,3 +37,16 @@ end
         @test k ∉ v
     end
 end
+
+@testset "Check reference bus in Radial Branches" begin
+    for name in ["matpower_ACTIVSg2000_sys", "matpower_ACTIVSg10k_sys"]
+        sys = build_system(MatpowerTestSystems, name; add_forecasts = false)
+        a_mat = IncidenceMatrix(sys)
+        rb = RadialBranches(IncidenceMatrix(sys))
+        leaf_buses = Int64[]
+        for i in keys(rb.bus_reduction_map)
+            append!(leaf_buses, collect(rb.bus_reduction_map[i]))
+        end
+        @test all(a_mat.ref_bus_positions .∉ leaf_buses)
+    end
+end
diff --git a/test/test_virtual_lodf_with_radial_branches.jl b/test/test_virtual_lodf_with_radial_branches.jl
new file mode 100644
index 00000000..69698c61
--- /dev/null
+++ b/test/test_virtual_lodf_with_radial_branches.jl
@@ -0,0 +1,18 @@
+@testset "Test VirtualLODF with radial lines" begin
+    # get the system
+    sys = PSB.build_system(PSB.PSITestSystems, "test_RTS_GMLC_sys")
+    # get the PTDF matrix for reference
+    lodf_rad = LODF(sys; reduce_radial_branches = true)
+    vlodf_rad = VirtualLODF(sys; reduce_radial_branches = true)
+
+    for i in axes(lodf_rad, 2)
+        virtual = vlodf_rad[i, :]
+        for j in axes(lodf_rad, 1)
+            # check values using PTDFs axes
+            @test isapprox(lodf_rad[i, j], vlodf_rad[i, j]; atol = 1e-10)
+        end
+    end
+    # Check the cache is populated
+    @test length(vlodf_rad.cache) == length(vlodf_rad.axes[1])
+    @test length(vlodf_rad.cache[1]) == length(vlodf_rad.axes[2])
+end
diff --git a/test/test_virtual_ptdf_with_radial_branches.jl b/test/test_virtual_ptdf_with_radial_branches.jl
new file mode 100644
index 00000000..5444a89e
--- /dev/null
+++ b/test/test_virtual_ptdf_with_radial_branches.jl
@@ -0,0 +1,36 @@
+@testset "Test VirtualPTDF with radial lines" begin
+    # get the system
+    sys = PSB.build_system(PSB.PSITestSystems, "test_RTS_GMLC_sys")
+    # get the PTDF matrix for reference
+    ptdf_rad = PTDF(sys; reduce_radial_branches = true)
+    vptdf_rad = VirtualPTDF(sys; reduce_radial_branches = true)
+
+    for i in axes(ptdf_rad, 2)
+        virtual = vptdf_rad[i, :]
+        for j in axes(ptdf_rad, 1)
+            # check values using PTDFs axes
+            @test isapprox(ptdf_rad[i, j], vptdf_rad[i, j]; atol = 1e-10)
+        end
+    end
+    # Check the cache is populated
+    @test length(vptdf_rad.cache) == length(vptdf_rad.axes[1])
+    @test length(vptdf_rad.cache[1]) == length(vptdf_rad.axes[2])
+end
+
+@testset "Test VirtualPTDF with radial lines and distributed slack" begin
+    # check if VirtualPTDF have same values as PTDF row-wise
+    # sys = PSB.build_system(PSB.PSITestSystems, "test_RTS_GMLC_sys")
+    sys = PSB.build_system(PSB.PSITestSystems, "c_sys14")   # ! remove
+    buscount = length(PNM.get_buses(sys))
+    dist_slack = 1 / buscount * ones(buscount)
+    slack_array = dist_slack / sum(dist_slack)
+    ptdf_rad = PTDF(sys; reduce_radial_branches = true, dist_slack = slack_array)
+    vptdf_rad = VirtualPTDF(sys; reduce_radial_branches = true, dist_slack = slack_array)
+    for i in axes(ptdf_rad, 2)
+        virtual = vptdf_rad[i, :]
+        for j in axes(ptdf_rad, 1)
+            # check values using PTDFs axes
+            @test isapprox(ptdf_rad[i, j], vptdf_rad[i, j]; atol = 1e-10)
+        end
+    end
+end