diff --git a/src/devices_models/device_constructors/branch_constructor.jl b/src/devices_models/device_constructors/branch_constructor.jl
index ec1b7d6d9d..616500cafb 100644
--- a/src/devices_models/device_constructors/branch_constructor.jl
+++ b/src/devices_models/device_constructors/branch_constructor.jl
@@ -106,7 +106,7 @@ function construct_device!(
             container,
             devices,
             device_model,
-            get_network_formulation(network_model),
+            network_model,
         )
     end
     return
@@ -273,7 +273,7 @@ function construct_device!(
         container,
         devices,
         device_model,
-        get_network_formulation(network_model),
+        network_model,
     )
     add_constraint_dual!(container, sys, device_model)
     return
@@ -331,7 +331,7 @@ function construct_device!(
 ) where {T <: PSY.ACBranch}
     devices =
         get_available_components(T, sys, get_attribute(model, "filter_function"))
-    branch_rate_bounds!(container, devices, model, get_network_formulation(network_model))
+    branch_rate_bounds!(container, devices, model, network_model)
 
     add_constraints!(container, RateLimitConstraintFromTo, devices, model, network_model)
     add_constraints!(container, RateLimitConstraintToFrom, devices, model, network_model)
@@ -356,7 +356,7 @@ function construct_device!(
 ) where {T <: PSY.ACBranch}
     devices =
         get_available_components(T, sys, get_attribute(model, "filter_function"))
-    branch_rate_bounds!(container, devices, model, get_network_formulation(network_model))
+    branch_rate_bounds!(container, devices, model, network_model)
     add_constraint_dual!(container, sys, model)
     return
 end
diff --git a/src/devices_models/devices/AC_branches.jl b/src/devices_models/devices/AC_branches.jl
index bf24ae3c3c..de74ab915c 100644
--- a/src/devices_models/devices/AC_branches.jl
+++ b/src/devices_models/devices/AC_branches.jl
@@ -104,11 +104,18 @@ function branch_rate_bounds!(
     container::OptimizationContainer,
     devices::IS.FlattenIteratorWrapper{B},
     ::DeviceModel{B, <:AbstractBranchFormulation},
-    ::Type{<:PM.AbstractDCPModel},
+    network_model::NetworkModel{<:PM.AbstractDCPModel},
 ) where {B <: PSY.ACBranch}
     var = get_variable(container, FlowActivePowerVariable(), B)
+
+    radial_branches = get_radial_branches(network_model)
+    radial_branches_names = PNM.get_radial_branches(radial_branches)
+
     for d in devices
         name = PSY.get_name(d)
+        if name ∈ radial_branches_names
+            continue
+        end
         for t in get_time_steps(container)
             JuMP.set_upper_bound(var[name, t], PSY.get_rate(d))
             JuMP.set_lower_bound(var[name, t], -1.0 * PSY.get_rate(d))
@@ -121,15 +128,23 @@ function branch_rate_bounds!(
     container::OptimizationContainer,
     devices::IS.FlattenIteratorWrapper{B},
     ::DeviceModel{B, <:AbstractBranchFormulation},
-    ::Type{<:PM.AbstractPowerModel},
+    network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {B <: PSY.ACBranch}
     vars = [
         get_variable(container, FlowActivePowerFromToVariable(), B),
         get_variable(container, FlowActivePowerToFromVariable(), B),
     ]
+
+    time_steps = get_time_steps(container)
+    radial_branches = get_radial_branches(network_model)
+    radial_branches_names = PNM.get_radial_branches(radial_branches)
+
     for d in devices
         name = PSY.get_name(d)
-        for t in get_time_steps(container), var in vars
+        if name ∈ radial_branches_names
+            continue
+        end
+        for t in time_steps, var in vars
             JuMP.set_upper_bound(var[name, t], PSY.get_rate(d))
             JuMP.set_lower_bound(var[name, t], -1.0 * PSY.get_rate(d))
         end
@@ -256,7 +271,7 @@ function add_constraints!(
     cons_type::Type{RateLimitConstraintFromTo},
     devices::IS.FlattenIteratorWrapper{B},
     ::DeviceModel{B, <:AbstractBranchFormulation},
-    ::NetworkModel{T},
+    network_model::NetworkModel{T},
 ) where {B <: PSY.ACBranch, T <: PM.AbstractPowerModel}
     rating_data = [(PSY.get_name(h), PSY.get_rate(h)) for h in devices]
 
@@ -272,7 +287,13 @@ function add_constraints!(
     )
     constraint = get_constraint(container, cons_type(), B)
 
+    radial_branches = get_radial_branches(network_model)
+    radial_branches_names = PNM.get_radial_branches(radial_branches)
+
     for r in rating_data
+        if r[1] ∈ radial_branches_names
+            continue
+        end
         for t in time_steps
             constraint[r[1], t] = JuMP.@constraint(
                 get_jump_model(container),
@@ -291,7 +312,7 @@ function add_constraints!(
     cons_type::Type{RateLimitConstraintToFrom},
     devices::IS.FlattenIteratorWrapper{B},
     ::DeviceModel{B, <:AbstractBranchFormulation},
-    ::NetworkModel{T},
+    network_model::NetworkModel{T},
 ) where {B <: PSY.ACBranch, T <: PM.AbstractPowerModel}
     rating_data = [(PSY.get_name(h), PSY.get_rate(h)) for h in devices]
 
@@ -307,7 +328,13 @@ function add_constraints!(
     )
     constraint = get_constraint(container, cons_type(), B)
 
+    radial_branches = get_radial_branches(network_model)
+    radial_branches_names = PNM.get_radial_branches(radial_branches)
+
     for r in rating_data
+        if r[1] ∈ radial_branches_names
+            continue
+        end
         for t in time_steps
             constraint[r[1], t] = JuMP.@constraint(
                 get_jump_model(container),
diff --git a/src/network_models/pm_translator.jl b/src/network_models/pm_translator.jl
index 627a5b7cf0..ec012ff097 100644
--- a/src/network_models/pm_translator.jl
+++ b/src/network_models/pm_translator.jl
@@ -381,7 +381,7 @@ end
 
 function get_branches_to_pm(
     sys::PSY.System,
-    ::Type{S},
+    network_model::NetworkModel{S},
     ::Type{T},
     branch_template::BranchModelContainer,
     start_idx = 0,
@@ -389,6 +389,8 @@ function get_branches_to_pm(
     PM_branches = Dict{String, Any}()
     PMmap_br = Dict{PM_MAP_TUPLE, T}()
 
+    radial_branches = get_radial_branches(network_model)
+    radial_branches_names = PNM.get_radial_branches(radial_branches)
     for (d, device_model) in branch_template
         comp_type = get_component_type(device_model)
         if comp_type <: TwoTerminalHVDCTypes
@@ -398,6 +400,10 @@ function get_branches_to_pm(
         start_idx += length(PM_branches)
         filter_func = get_attribute(device_model, "filter_function")
         for (i, branch) in enumerate(get_available_components(comp_type, sys, filter_func))
+            if PSY.get_name(branch) ∈ radial_branches_names
+                @debug "Skipping branch $(PSY.get_name(branch)) since its radial"
+                continue
+            end
             ix = i + start_idx
             PM_branches["$(ix)"] =
                 get_branch_to_pm(ix, branch, get_formulation(device_model), S)
@@ -413,7 +419,7 @@ end
 
 function get_branches_to_pm(
     sys::PSY.System,
-    ::Type{S},
+    network_model::NetworkModel{S},
     ::Type{T},
     branch_template::BranchModelContainer,
     start_idx = 0,
@@ -442,7 +448,7 @@ end
 
 function get_branches_to_pm(
     ::PSY.System,
-    ::Type{PTDFPowerModel},
+    network_model::NetworkModel{PTDFPowerModel},
     ::Type{T},
     branch_template::BranchModelContainer,
     start_idx = 0,
@@ -485,13 +491,13 @@ end
 function pass_to_pm(sys::PSY.System, template::ProblemTemplate, time_periods::Int)
     ac_lines, PMmap_ac = get_branches_to_pm(
         sys,
-        get_network_formulation(template),
+        get_network_model(template),
         PSY.ACBranch,
         template.branches,
     )
     two_terminal_dc_lines, PMmap_dc = get_branches_to_pm(
         sys,
-        get_network_formulation(template),
+        get_network_model(template),
         TwoTerminalHVDCTypes,
         template.branches,
         length(ac_lines),
diff --git a/src/network_models/powermodels_interface.jl b/src/network_models/powermodels_interface.jl
index 932f2c158f..942a735985 100644
--- a/src/network_models/powermodels_interface.jl
+++ b/src/network_models/powermodels_interface.jl
@@ -291,17 +291,25 @@ function powermodels_network!(
 ) where {S <: PM.AbstractPowerModel}
     time_steps = get_time_steps(container)
     pm_data, PM_map = pass_to_pm(sys, template, time_steps[end])
-    buses = get_available_components(PSY.ACBus, sys)
 
-    for t in time_steps, bus in buses
-        pm_data["nw"]["$(t)"]["bus"]["$(bus.number)"]["inj_p"] =
+    network_model = get_network_model(template)
+    radial_branches = get_radial_branches(network_model)
+    if isempty(radial_branches)
+        ac_bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys))
+    else
+        bus_reduction_map = radial_branches.bus_reduction_map
+        ac_bus_numbers = collect(keys(bus_reduction_map))
+    end
+
+    for t in time_steps, bus_no in ac_bus_numbers
+        pm_data["nw"]["$(t)"]["bus"]["$bus_no"]["inj_p"] =
             container.expressions[ExpressionKey(ActivePowerBalance, PSY.ACBus)][
-                bus.number,
+                bus_no,
                 t,
             ]
-        pm_data["nw"]["$(t)"]["bus"]["$(bus.number)"]["inj_q"] =
+        pm_data["nw"]["$(t)"]["bus"]["$bus_no"]["inj_q"] =
             container.expressions[ExpressionKey(ReactivePowerBalance, PSY.ACBus)][
-                bus.number,
+                bus_no,
                 t,
             ]
     end
@@ -324,10 +332,19 @@ function powermodels_network!(
     pm_data, PM_map = pass_to_pm(sys, template, time_steps[end])
     buses = get_available_components(PSY.ACBus, sys)
 
-    for t in time_steps, bus in buses
-        pm_data["nw"]["$(t)"]["bus"]["$(PSY.get_number(bus))"]["inj_p"] =
+    network_model = get_network_model(template)
+    radial_branches = get_radial_branches(network_model)
+    if isempty(radial_branches)
+        ac_bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys))
+    else
+        bus_reduction_map = radial_branches.bus_reduction_map
+        ac_bus_numbers = collect(keys(bus_reduction_map))
+    end
+
+    for t in time_steps, bus_no in ac_bus_numbers
+        pm_data["nw"]["$(t)"]["bus"]["$bus_no"]["inj_p"] =
             container.expressions[ExpressionKey(ActivePowerBalance, PSY.ACBus)][
-                PSY.get_number(bus),
+                bus_no,
                 t,
             ]
         # pm_data["nw"]["$(t)"]["bus"]["$(bus.number)"]["inj_q"] = 0.0