diff --git a/LICENSE b/LICENSE
index 8b3f7b22d0..f0323d2e2c 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,6 +1,6 @@
 BSD 3-Clause License
 
-Copyright (c) 2018, 2023 Alliance for Sustainable Energy, LLC
+Copyright (c) 2018, 2023 Alliance for Sustainable Energy, LLC and The Regents of the University of California
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
diff --git a/Project.toml b/Project.toml
index 92630cf604..fd1cf39fb9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "PowerSimulations"
 uuid = "e690365d-45e2-57bb-ac84-44ba829e73c4"
 authors = ["Jose Daniel Lara", "Clayton Barrows", "Daniel Thom", "Dheepak Krishnamurthy", "Sourabh Dalvi"]
-version = "0.25.2"
+version = "0.26.0"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
@@ -45,13 +45,13 @@ JuMP = "1"
 LinearAlgebra = "1"
 Logging = "1"
 MathOptInterface = "1"
-PowerModels = "~0.19"
-PowerNetworkMatrices = "^0.9"
+PowerModels = "^0.20"
+PowerNetworkMatrices = "^0.10"
 PowerSystems = "^3"
 PrettyTables = "2"
 ProgressMeter = "^1.5"
 SHA = "0.7"
 Serialization = "1"
-TimeSeries = "~0.23"
+TimeSeries = "~0.23, ~0.24"
 TimerOutputs = "~0.5"
 julia = "^1.6"
diff --git a/src/core/definitions.jl b/src/core/definitions.jl
index b51e41e7db..142967f063 100644
--- a/src/core/definitions.jl
+++ b/src/core/definitions.jl
@@ -53,9 +53,6 @@ const KiB = 1024
 const MiB = KiB * KiB
 const GiB = MiB * KiB
 
-const UNSUPPORTED_POWERMODELS =
-    [PM.SOCBFPowerModel, PM.SOCBFConicPowerModel, PM.IVRPowerModel]
-
 const PSI_NAME_DELIMITER = "__"
 
 const M_VALUE = 1e6
diff --git a/src/core/network_model.jl b/src/core/network_model.jl
index 058d512f77..caa2d8873f 100644
--- a/src/core/network_model.jl
+++ b/src/core/network_model.jl
@@ -20,7 +20,7 @@ Establishes the model for a particular device specified by type.
   - `use_slacks::Bool`: Adds slacks to the network modelings
   - `PTDF::PTDF`: PTDF Array calculated using PowerNetworkMatrices
   - `duals::Vector{DataType}`: Constraint types to calculate the duals
-
+  - `reduce_radial_branches::Bool`: Skips modeling radial branches in the system to reduce problem size
 # Example
 
 ptdf_array = PTDF(system)
@@ -32,21 +32,34 @@ mutable struct NetworkModel{T <: PM.AbstractPowerModel}
     subnetworks::Dict{Int, Set{Int}}
     bus_area_map::Dict{PSY.ACBus, Int}
     duals::Vector{DataType}
+    radial_network_reduction::PNM.RadialNetworkReduction
+    reduce_radial_branches::Bool
 
     function NetworkModel(
         ::Type{T};
         use_slacks = false,
         PTDF_matrix = nothing,
+        reduce_radial_branches = false,
         subnetworks = Dict{Int, Set{Int}}(),
         duals = Vector{DataType}(),
     ) where {T <: PM.AbstractPowerModel}
         _check_pm_formulation(T)
-        new{T}(use_slacks, PTDF_matrix, subnetworks, Dict{PSY.ACBus, Int}(), duals)
+        new{T}(
+            use_slacks,
+            PTDF_matrix,
+            subnetworks,
+            Dict{PSY.ACBus, Int}(),
+            duals,
+            PNM.RadialNetworkReduction(),
+            reduce_radial_branches,
+        )
     end
 end
 
 get_use_slacks(m::NetworkModel) = m.use_slacks
 get_PTDF_matrix(m::NetworkModel) = m.PTDF_matrix
+get_reduce_radial_branches(m::NetworkModel) = m.reduce_radial_branches
+get_radial_network_reduction(m::NetworkModel) = m.radial_network_reduction
 get_duals(m::NetworkModel) = m.duals
 get_network_formulation(::NetworkModel{T}) where {T} = T
 get_reference_buses(m::NetworkModel{T}) where {T <: PM.AbstractPowerModel} =
@@ -62,6 +75,15 @@ function add_dual!(model::NetworkModel, dual)
     return
 end
 
+function check_radial_branch_reduction_compatibility(
+    ::Type{T},
+) where {T <: PM.AbstractPowerModel}
+    if T ∈ INCOMPATIBLE_WITH_RADIAL_BRANCHES_POWERMODELS
+        error("Network Model $T is not compatible with radial branch reduction")
+    end
+    return
+end
+
 function instantiate_network_model(
     model::NetworkModel{T},
     sys::PSY.System,
@@ -69,6 +91,10 @@ function instantiate_network_model(
     if isempty(model.subnetworks)
         model.subnetworks = PNM.find_subnetworks(sys)
     end
+    if model.reduce_radial_branches
+        check_radial_branch_reduction_compatibility(T)
+        model.radial_network_reduction = PNM.RadialNetworkReduction(sys)
+    end
     return
 end
 
@@ -90,7 +116,12 @@ end
 function instantiate_network_model(model::NetworkModel{StandardPTDFModel}, sys::PSY.System)
     if get_PTDF_matrix(model) === nothing
         @info "PTDF Matrix not provided. Calculating using PowerNetworkMatrices.PTDF"
-        model.PTDF_matrix = PNM.PTDF(sys)
+        model.PTDF_matrix =
+            PNM.PTDF(sys; reduce_radial_branches = model.reduce_radial_branches)
+    end
+    if model.reduce_radial_branches
+        @assert !isempty(model.PTDF_matrix.radial_network_reduction)
+        model.radial_network_reduction = model.PTDF_matrix.radial_network_reduction
     end
     get_PTDF_matrix(model).subnetworks
     model.subnetworks = deepcopy(get_PTDF_matrix(model).subnetworks)
@@ -107,19 +138,27 @@ function _assign_subnetworks_to_buses(
 ) where {T <: Union{CopperPlatePowerModel, StandardPTDFModel}}
     subnetworks = model.subnetworks
     temp_bus_map = Dict{Int, Int}()
-    for bus in get_available_components(PSY.ACBus, sys)
+    radial_network_reduction = PSI.get_radial_network_reduction(model)
+    for bus in PSI.get_available_components(PSY.ACBus, sys)
         bus_no = PSY.get_number(bus)
+        mapped_bus_no = PNM.get_mapped_bus_number(radial_network_reduction, bus)
         if haskey(temp_bus_map, bus_no)
             model.bus_area_map[bus] = temp_bus_map[bus_no]
+            continue
         else
+            bus_mapped = false
             for (subnet, bus_set) in subnetworks
-                if bus_no ∈ bus_set
+                if mapped_bus_no ∈ bus_set
                     temp_bus_map[bus_no] = subnet
                     model.bus_area_map[bus] = subnet
+                    bus_mapped = true
                     break
                 end
             end
         end
+        if !bus_mapped
+            error("Bus $(PSY.summary(bus)) not mapped to any reference bus")
+        end
     end
     return
 end
diff --git a/src/core/optimization_container.jl b/src/core/optimization_container.jl
index be85ae1a9c..9686a1d06d 100644
--- a/src/core/optimization_container.jl
+++ b/src/core/optimization_container.jl
@@ -431,9 +431,14 @@ function _make_system_expressions!(
     subnetworks::Dict{Int, Set{Int}},
     dc_bus_numbers::Vector{Int},
     ::Type{<:PM.AbstractPowerModel},
+    bus_reduction_map::Dict{Int64, Set{Int64}},
 )
     time_steps = get_time_steps(container)
-    ac_bus_numbers = collect(Iterators.flatten(values(subnetworks)))
+    if isempty(bus_reduction_map)
+        ac_bus_numbers = collect(Iterators.flatten(values(subnetworks)))
+    else
+        ac_bus_numbers = collect(keys(bus_reduction_map))
+    end
     container.expressions = Dict(
         ExpressionKey(ActivePowerBalance, PSY.ACBus) =>
             _make_container_array(ac_bus_numbers, time_steps),
@@ -450,9 +455,14 @@ function _make_system_expressions!(
     subnetworks::Dict{Int, Set{Int}},
     dc_bus_numbers::Vector{Int},
     ::Type{<:PM.AbstractActivePowerModel},
+    bus_reduction_map::Dict{Int64, Set{Int64}},
 )
     time_steps = get_time_steps(container)
-    ac_bus_numbers = collect(Iterators.flatten(values(subnetworks)))
+    if isempty(bus_reduction_map)
+        ac_bus_numbers = collect(Iterators.flatten(values(subnetworks)))
+    else
+        ac_bus_numbers = collect(keys(bus_reduction_map))
+    end
     container.expressions = Dict(
         ExpressionKey(ActivePowerBalance, PSY.ACBus) =>
             _make_container_array(ac_bus_numbers, time_steps),
@@ -467,6 +477,7 @@ function _make_system_expressions!(
     subnetworks::Dict{Int, Set{Int}},
     ::Vector{Int},
     ::Type{CopperPlatePowerModel},
+    bus_reduction_map::Dict{Int64, Set{Int64}},
 )
     time_steps = get_time_steps(container)
     subnetworks_ref_buses = collect(keys(subnetworks))
@@ -482,9 +493,14 @@ function _make_system_expressions!(
     subnetworks::Dict{Int, Set{Int}},
     dc_bus_numbers::Vector{Int},
     ::Type{T},
+    bus_reduction_map::Dict{Int64, Set{Int64}},
 ) where {T <: Union{PTDFPowerModel, StandardPTDFModel}}
     time_steps = get_time_steps(container)
-    ac_bus_numbers = sort!(collect(Iterators.flatten(values(subnetworks))))
+    if isempty(bus_reduction_map)
+        ac_bus_numbers = collect(Iterators.flatten(values(subnetworks)))
+    else
+        ac_bus_numbers = collect(keys(bus_reduction_map))
+    end
     subnetworks = collect(keys(subnetworks))
     container.expressions = Dict(
         ExpressionKey(ActivePowerBalance, PSY.System) =>
@@ -492,7 +508,9 @@ function _make_system_expressions!(
         ExpressionKey(ActivePowerBalance, PSY.DCBus) =>
             _make_container_array(dc_bus_numbers, time_steps),
         ExpressionKey(ActivePowerBalance, PSY.ACBus) =>
-            _make_container_array(ac_bus_numbers, time_steps),
+        # Bus numbers are sorted to guarantee consistency in the order between the
+        # containers
+            _make_container_array(sort!(ac_bus_numbers), time_steps),
     )
     return
 end
@@ -502,9 +520,10 @@ function initialize_system_expressions!(
     ::Type{T},
     subnetworks::Dict{Int, Set{Int}},
     system::PSY.System,
+    bus_reduction_map::Dict{Int64, Set{Int64}},
 ) where {T <: PM.AbstractPowerModel}
     dc_bus_numbers = [PSY.get_number(b) for b in PSY.get_components(PSY.DCBus, system)]
-    _make_system_expressions!(container, subnetworks, dc_bus_numbers, T)
+    _make_system_expressions!(container, subnetworks, dc_bus_numbers, T, bus_reduction_map)
     return
 end
 
@@ -520,7 +539,7 @@ function build_impl!(
         transmission,
         transmission_model.subnetworks,
         sys,
-    )
+        transmission_model.radial_network_reduction.bus_reduction_map)
 
     # Order is required
     for device_model in values(template.devices)
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 1176f83700..e14215cec1 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_network_reduction = get_radial_network_reduction(network_model)
+    radial_branches_names = PNM.get_radial_branches(radial_network_reduction)
+
     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_network_reduction = get_radial_network_reduction(network_model)
+    radial_branches_names = PNM.get_radial_branches(radial_network_reduction)
+
     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
@@ -168,14 +183,55 @@ function add_constraints!(
     container::OptimizationContainer,
     cons_type::Type{RateLimitConstraint},
     devices::IS.FlattenIteratorWrapper{T},
-    model::DeviceModel{T, U},
-    ::NetworkModel{X},
+    ::DeviceModel{T, U},
+    network_model::NetworkModel{V},
 ) where {
     T <: PSY.ACBranch,
     U <: AbstractBranchFormulation,
-    X <: PM.AbstractActivePowerModel,
+    V <: PM.AbstractActivePowerModel,
 }
-    add_range_constraints!(container, cons_type, FlowActivePowerVariable, devices, model, X)
+    time_steps = get_time_steps(container)
+    radial_network_reduction = get_radial_network_reduction(network_model)
+    if isempty(radial_network_reduction)
+        device_names = [PSY.get_name(d) for d in devices]
+    else
+        device_names = PNM.get_meshed_branches(radial_network_reduction)
+    end
+
+    con_lb =
+        add_constraints_container!(
+            container,
+            cons_type(),
+            T,
+            device_names,
+            time_steps;
+            meta = "lb",
+        )
+    con_ub =
+        add_constraints_container!(
+            container,
+            cons_type(),
+            T,
+            device_names,
+            time_steps;
+            meta = "ub",
+        )
+
+    array = get_variable(container, FlowActivePowerVariable(), T)
+
+    for device in devices
+        ci_name = PSY.get_name(device)
+        if ci_name ∈ PNM.get_radial_branches(radial_network_reduction)
+            continue
+        end
+        limits = get_min_max_limits(device, RateLimitConstraint, U) # depends on constraint type and formulation type
+        for t in time_steps
+            con_ub[ci_name, t] =
+                JuMP.@constraint(container.JuMPmodel, array[ci_name, t] <= limits.max)
+            con_lb[ci_name, t] =
+                JuMP.@constraint(container.JuMPmodel, array[ci_name, t] >= limits.min)
+        end
+    end
     return
 end
 
@@ -215,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]
 
@@ -231,7 +287,13 @@ function add_constraints!(
     )
     constraint = get_constraint(container, cons_type(), B)
 
+    radial_network_reduction = get_radial_network_reduction(network_model)
+    radial_branches_names = PNM.get_radial_branches(radial_network_reduction)
+
     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),
@@ -250,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]
 
@@ -266,7 +328,13 @@ function add_constraints!(
     )
     constraint = get_constraint(container, cons_type(), B)
 
+    radial_network_reduction = get_radial_network_reduction(network_model)
+    radial_branches_names = PNM.get_radial_branches(radial_network_reduction)
+
     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/devices_models/devices/common/add_to_expression.jl b/src/devices_models/devices/common/add_to_expression.jl
index 54aa14a590..c9c7e1cb68 100644
--- a/src/devices_models/devices/common/add_to_expression.jl
+++ b/src/devices_models/devices/common/add_to_expression.jl
@@ -79,7 +79,7 @@ function add_to_expression!(
     ::Type{U},
     devices::IS.FlattenIteratorWrapper{V},
     model::DeviceModel{V, W},
-    ::NetworkModel{X},
+    network_model::NetworkModel{X},
 ) where {
     T <: SystemBalanceExpressions,
     U <: TimeSeriesParameter,
@@ -89,11 +89,12 @@ function add_to_expression!(
 }
     param_container = get_parameter(container, U(), V)
     multiplier = get_multiplier_array(param_container)
+    radial_network_reduction = get_radial_network_reduction(network_model)
     for d in devices, t in get_time_steps(container)
-        bus_number = PSY.get_number(PSY.get_bus(d))
+        bus_no = PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_bus(d))
         name = PSY.get_name(d)
         _add_to_jump_expression!(
-            get_expression(container, T(), PSY.ACBus)[bus_number, t],
+            get_expression(container, T(), PSY.ACBus)[bus_no, t],
             get_parameter_column_refs(param_container, name)[t],
             multiplier[name, t],
         )
@@ -107,7 +108,7 @@ function add_to_expression!(
     ::Type{U},
     devices::IS.FlattenIteratorWrapper{V},
     ::DeviceModel{V, W},
-    ::NetworkModel{X},
+    network_model::NetworkModel{X},
 ) where {
     T <: ActivePowerBalance,
     U <: OnStatusParameter,
@@ -116,13 +117,13 @@ function add_to_expression!(
     X <: PM.AbstractPowerModel,
 }
     parameter = get_parameter_array(container, U(), V)
-
+    radial_network_reduction = get_radial_network_reduction(network_model)
     for d in devices, t in get_time_steps(container)
-        bus_number = PSY.get_number(PSY.get_bus(d))
+        bus_no = PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_bus(d))
         name = PSY.get_name(d)
         mult = get_expression_multiplier(U(), T(), d, W())
         _add_to_jump_expression!(
-            get_expression(container, T(), PSY.ACBus)[bus_number, t],
+            get_expression(container, T(), PSY.ACBus)[bus_no, t],
             parameter[name, t],
             mult,
         )
@@ -139,7 +140,7 @@ function add_to_expression!(
     ::Type{U},
     devices::IS.FlattenIteratorWrapper{V},
     ::DeviceModel{V, W},
-    ::NetworkModel{X},
+    network_model::NetworkModel{X},
 ) where {
     T <: SystemBalanceExpressions,
     U <: VariableType,
@@ -149,11 +150,12 @@ function add_to_expression!(
 }
     variable = get_variable(container, U(), V)
     expression = get_expression(container, T(), PSY.ACBus)
+    radial_network_reduction = get_radial_network_reduction(network_model)
     for d in devices, t in get_time_steps(container)
         name = PSY.get_name(d)
-        bus_number = PSY.get_number(PSY.get_bus(d))
+        bus_no = PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_bus(d))
         _add_to_jump_expression!(
-            expression[bus_number, t],
+            expression[bus_no, t],
             variable[name, t],
             get_variable_multiplier(U(), V, W()),
         )
@@ -215,8 +217,9 @@ function add_to_expression!(
     var = get_variable(container, U(), V)
     nodal_expr = get_expression(container, T(), PSY.ACBus)
     sys_expr = get_expression(container, T(), PSY.System)
+    radial_network_reduction = get_radial_network_reduction(network_model)
     for d in devices
-        bus_no_to = PSY.get_number(PSY.get_arc(d).to)
+        bus_no_to = PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).to)
         ref_bus_from = get_reference_bus(network_model, PSY.get_arc(d).from)
         ref_bus_to = get_reference_bus(network_model, PSY.get_arc(d).to)
         for t in get_time_steps(container)
@@ -250,8 +253,10 @@ function add_to_expression!(
     var = get_variable(container, U(), V)
     nodal_expr = get_expression(container, T(), PSY.ACBus)
     sys_expr = get_expression(container, T(), PSY.System)
+    radial_network_reduction = get_radial_network_reduction(network_model)
     for d in devices
-        bus_no_from = PSY.get_number(PSY.get_arc(d).from)
+        bus_no_from =
+            PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).from)
         ref_bus_to = get_reference_bus(network_model, PSY.get_arc(d).to)
         ref_bus_from = get_reference_bus(network_model, PSY.get_arc(d).from)
         for t in get_time_steps(container)
@@ -352,12 +357,14 @@ function add_to_expression!(
 }
     variable = get_variable(container, U(), V)
     expression = get_expression(container, T(), PSY.ACBus)
+    radial_network_reduction = get_radial_network_reduction(network_model)
     for d in devices
         name = PSY.get_name(d)
-        bus_number = PSY.get_number(PSY.get_arc(d).from)
+        bus_no_ = PSY.get_number(PSY.get_arc(d).from)
+        bus_no = PNM.get_mapped_bus_number(radial_network_reduction, bus_no_)
         for t in get_time_steps(container)
             _add_to_jump_expression!(
-                expression[bus_number, t],
+                expression[bus_no, t],
                 variable[name, t],
                 get_variable_multiplier(U(), V, W()),
             )
@@ -375,7 +382,7 @@ function add_to_expression!(
     ::Type{U},
     devices::IS.FlattenIteratorWrapper{V},
     ::DeviceModel{V, W},
-    ::NetworkModel{X},
+    network_model::NetworkModel{X},
 ) where {
     T <: ActivePowerBalance,
     U <: FlowActivePowerToFromVariable,
@@ -385,12 +392,14 @@ function add_to_expression!(
 }
     variable = get_variable(container, U(), V)
     expression = get_expression(container, T(), PSY.ACBus)
+    radial_network_reduction = get_radial_network_reduction(network_model)
     for d in devices
         name = PSY.get_name(d)
-        bus_number = PSY.get_number(PSY.get_arc(d).to)
+        bus_no_ = PSY.get_number(PSY.get_arc(d).to)
+        bus_no = PNM.get_mapped_bus_number(radial_network_reduction, bus_no_)
         for t in get_time_steps(container)
             _add_to_jump_expression!(
-                expression[bus_number, t],
+                expression[bus_no, t],
                 variable[name, t],
                 get_variable_multiplier(U(), V, W()),
             )
@@ -405,7 +414,7 @@ function add_to_expression!(
     ::Type{U},
     devices::IS.FlattenIteratorWrapper{V},
     ::DeviceModel{V, W},
-    ::NetworkModel{X},
+    network_model::NetworkModel{X},
 ) where {
     T <: SystemBalanceExpressions,
     U <: OnVariable,
@@ -415,11 +424,13 @@ function add_to_expression!(
 }
     variable = get_variable(container, U(), V)
     expression = get_expression(container, T(), PSY.ACBus)
+    radial_network_reduction = get_radial_network_reduction(network_model)
     for d in devices, t in get_time_steps(container)
         name = PSY.get_name(d)
-        bus_number = PSY.get_number(PSY.get_bus(d))
+        bus_no_ = PSY.get_number(PSY.get_bus(d))
+        bus_no = PNM.get_mapped_bus_number(radial_network_reduction, bus_no_)
         _add_to_jump_expression!(
-            expression[bus_number, t],
+            expression[bus_no, t],
             variable[name, t],
             get_variable_multiplier(U(), d, W()),
         )
@@ -574,10 +585,12 @@ function add_to_expression!(
     multiplier = get_multiplier_array(param_container)
     sys_expr = get_expression(container, T(), PSY.System)
     nodal_expr = get_expression(container, T(), PSY.ACBus)
+    radial_network_reduction = get_radial_network_reduction(network_model)
     for d in devices
         name = PSY.get_name(d)
         device_bus = PSY.get_bus(d)
-        bus_no = PSY.get_number(device_bus)
+        bus_no_ = PSY.get_number(device_bus)
+        bus_no = PNM.get_mapped_bus_number(radial_network_reduction, bus_no_)
         ref_bus = get_reference_bus(network_model, device_bus)
         param = get_parameter_column_refs(param_container, name)
         for t in get_time_steps(container)
@@ -605,9 +618,11 @@ function add_to_expression!(
     parameter = get_parameter_array(container, U(), V)
     sys_expr = get_expression(container, T(), PSY.System)
     nodal_expr = get_expression(container, T(), PSY.ACBus)
+    radial_network_reduction = get_radial_network_reduction(network_model)
     for d in devices, t in get_time_steps(container)
         name = PSY.get_name(d)
-        bus_no = PSY.get_number(PSY.get_bus(d))
+        bus_no_ = PSY.get_number(PSY.get_bus(d))
+        bus_no = PNM.get_mapped_bus_number(radial_network_reduction, bus_no_)
         mult = get_expression_multiplier(U(), T(), d, W())
         device_bus = PSY.get_bus(d)
         ref_bus = get_reference_bus(network_model, device_bus)
@@ -637,10 +652,11 @@ function add_to_expression!(
     variable = get_variable(container, U(), V)
     sys_expr = get_expression(container, T(), PSY.System)
     nodal_expr = get_expression(container, T(), PSY.ACBus)
+    radial_network_reduction = get_radial_network_reduction(network_model)
     for d in devices
         name = PSY.get_name(d)
         device_bus = PSY.get_bus(d)
-        bus_no = PSY.get_number(device_bus)
+        bus_no = PNM.get_mapped_bus_number(radial_network_reduction, device_bus)
         ref_bus = get_reference_bus(network_model, device_bus)
         for t in get_time_steps(container)
             _add_to_jump_expression!(
@@ -675,11 +691,11 @@ function add_to_expression!(
     variable = get_variable(container, U(), V)
     sys_expr = get_expression(container, T(), PSY.System)
     nodal_expr = get_expression(container, T(), PSY.ACBus)
+    radial_network_reduction = get_radial_network_reduction(network_model)
     for d in devices
         name = PSY.get_name(d)
-        device_bus = PSY.get_bus(d)
-        bus_no = PSY.get_number(device_bus)
-        ref_bus = get_reference_bus(network_model, device_bus)
+        bus_no = PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_bus(d))
+        ref_bus = get_reference_bus(network_model, PSY.get_bus(d))
         for t in get_time_steps(container)
             _add_to_jump_expression!(
                 sys_expr[ref_bus, t],
@@ -715,16 +731,20 @@ function add_to_expression!(
 }
     var = get_variable(container, U(), V)
     expression = get_expression(container, T(), PSY.ACBus)
+    radial_network_reduction = get_radial_network_reduction(network_model)
     for d in devices
+        bus_no_from =
+            PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).from)
+        bus_no_to = PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).to)
         for t in get_time_steps(container)
             flow_variable = var[PSY.get_name(d), t]
             _add_to_jump_expression!(
-                expression[PSY.get_number(PSY.get_arc(d).from), t],
+                expression[bus_no_from, t],
                 flow_variable,
                 -1.0,
             )
             _add_to_jump_expression!(
-                expression[PSY.get_number(PSY.get_arc(d).to), t],
+                expression[bus_no_to, t],
                 flow_variable,
                 1.0,
             )
@@ -752,9 +772,11 @@ function add_to_expression!(
     var = get_variable(container, U(), V)
     nodal_expr = get_expression(container, T(), PSY.ACBus)
     sys_expr = get_expression(container, T(), PSY.System)
+    radial_network_reduction = get_radial_network_reduction(network_model)
     for d in devices
-        bus_no_from = PSY.get_number(PSY.get_arc(d).from)
-        bus_no_to = PSY.get_number(PSY.get_arc(d).to)
+        bus_no_from =
+            PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).from)
+        bus_no_to = PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).to)
         ref_bus_from = get_reference_bus(network_model, PSY.get_arc(d).from)
         ref_bus_to = get_reference_bus(network_model, PSY.get_arc(d).to)
         for t in get_time_steps(container)
@@ -826,16 +848,20 @@ function add_to_expression!(
 ) where {T <: ActivePowerBalance, U <: PhaseShifterAngle, V <: PhaseAngleControl}
     var = get_variable(container, U(), PSY.PhaseShiftingTransformer)
     expression = get_expression(container, T(), PSY.ACBus)
+    radial_network_reduction = get_radial_network_reduction(network_model)
     for d in devices
+        bus_no_from =
+            PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).from)
+        bus_no_to = PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).to)
         for t in get_time_steps(container)
             flow_variable = var[PSY.get_name(d), t]
             _add_to_jump_expression!(
-                expression[PSY.get_number(PSY.get_arc(d).from), t],
+                expression[bus_no_from, t],
                 flow_variable,
                 -get_variable_multiplier(U(), d, V()),
             )
             _add_to_jump_expression!(
-                expression[PSY.get_number(PSY.get_arc(d).to), t],
+                expression[bus_no_to, t],
                 flow_variable,
                 get_variable_multiplier(U(), d, V()),
             )
@@ -1079,7 +1105,8 @@ function add_to_expression!(
     variable = get_variable(container, U(), PSY.ACBus)
     expression = get_expression(container, T(), PSY.ACBus)
     @assert_op length(axes(variable, 1)) == length(axes(expression, 1))
-    for t in get_time_steps(container), n in axes(variable, 1)
+    # We uses axis here to avoid double addition of the slacks to the aggregated buses
+    for t in get_time_steps(container), n in axes(expression, 1)
         _add_to_jump_expression!(
             expression[n, t],
             variable[n, t],
@@ -1102,8 +1129,8 @@ function add_to_expression!(
 }
     variable = get_variable(container, U(), PSY.ACBus, "P")
     expression = get_expression(container, T(), PSY.ACBus)
-    bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys))
-    for t in get_time_steps(container), n in bus_numbers
+    # We uses axis here to avoid double addition of the slacks to the aggregated buses
+    for t in get_time_steps(container), n in axes(expression, 1)
         _add_to_jump_expression!(
             expression[n, t],
             variable[n, t],
@@ -1126,8 +1153,8 @@ function add_to_expression!(
 }
     variable = get_variable(container, U(), PSY.ACBus, "Q")
     expression = get_expression(container, T(), PSY.ACBus)
-    bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys))
-    for t in get_time_steps(container), n in bus_numbers
+    # We uses axis here to avoid double addition of the slacks to the aggregated buses
+    for t in get_time_steps(container), n in axes(expression, 1)
         _add_to_jump_expression!(
             expression[n, t],
             variable[n, t],
diff --git a/src/devices_models/devices/thermal_generation.jl b/src/devices_models/devices/thermal_generation.jl
index b58d01f965..95e3733aac 100644
--- a/src/devices_models/devices/thermal_generation.jl
+++ b/src/devices_models/devices/thermal_generation.jl
@@ -341,11 +341,21 @@ end
 Startup and shutdown active power limits for Compact Unit Commitment
 """
 function get_startup_shutdown_limits(
-    device,
+    device::PSY.ThermalMultiStart,
     ::Type{ActivePowerVariableLimitsConstraint},
-    ::Type{<:ThermalMultiStartUnitCommitment},
+    ::Type{ThermalMultiStartUnitCommitment},
 )
-    return PSY.get_power_trajectory(device)
+    startup_shutdown = PSY.get_power_trajectory(device)
+    if isnothing(startup_shutdown)
+        @warn(
+            "Generator $(summary(device)) has a Nothing startup_shutdown property. Using active power limits."
+        )
+        return (
+            startup = PSY.get_active_power_limits(device).max,
+            shutdown = PSY.get_active_power_limits(device).max,
+        )
+    end
+    return startup_shutdown
 end
 
 """
@@ -355,7 +365,7 @@ function get_min_max_limits(
     device,
     ::Type{ActivePowerVariableLimitsConstraint},
     ::Type{<:AbstractCompactUnitCommitment},
-) #  -> Union{Nothing, NamedTuple{(:startup, :shutdown), Tuple{Float64, Float64}}}
+) #  -> Union{Nothing, NamedTuple{(:min, :max), Tuple{Float64, Float64}}}
     return (
         min = 0,
         max = PSY.get_active_power_limits(device).max -
@@ -562,6 +572,7 @@ function add_constraints!(
         name = PSY.get_name(device)
         limits = get_min_max_limits(device, T, W) # depends on constraint type and formulation type
         startup_shutdown_limits = get_startup_shutdown_limits(device, T, W)
+        @assert !isnothing(startup_shutdown_limits) "$(name)"
         if JuMP.has_lower_bound(varp[name, t])
             JuMP.set_lower_bound(varp[name, t], 0.0)
         end
diff --git a/src/initial_conditions/initialization.jl b/src/initial_conditions/initialization.jl
index 80ff0f1fd7..e5222f6bfc 100644
--- a/src/initial_conditions/initialization.jl
+++ b/src/initial_conditions/initialization.jl
@@ -5,7 +5,12 @@ function get_initial_conditions_template(model::OperationModel)
         get_network_formulation(model.template);
         use_slacks = get_use_slacks(get_network_model(model.template)),
         PTDF_matrix = get_PTDF_matrix(get_network_model(model.template)),
+        reduce_radial_branches = get_reduce_radial_branches(
+            get_network_model(model.template),
+        ),
     )
+    network_model.radial_network_reduction =
+        get_radial_network_reduction(get_network_model(model.template))
     network_model.subnetworks = get_subnetworks(get_network_model(model.template))
     bus_area_map = get_bus_area_map(get_network_model(model.template))
     if !isempty(bus_area_map)
diff --git a/src/network_models/network_constructor.jl b/src/network_models/network_constructor.jl
index e2d4a5d850..7f53454439 100644
--- a/src/network_models/network_constructor.jl
+++ b/src/network_models/network_constructor.jl
@@ -202,15 +202,14 @@ function construct_network!(
     end
 
     if get_use_slacks(model)
-        add_variables!(container, SystemBalanceSlackUp, sys, T)
-        add_variables!(container, SystemBalanceSlackDown, sys, T)
+        add_variables!(container, SystemBalanceSlackUp, sys, model)
+        add_variables!(container, SystemBalanceSlackDown, sys, model)
         add_to_expression!(
             container,
             ActivePowerBalance,
             SystemBalanceSlackUp,
             sys,
             model,
-            T,
         )
         add_to_expression!(
             container,
@@ -218,7 +217,6 @@ function construct_network!(
             SystemBalanceSlackDown,
             sys,
             model,
-            T,
         )
         add_to_expression!(
             container,
@@ -226,7 +224,6 @@ function construct_network!(
             SystemBalanceSlackUp,
             sys,
             model,
-            T,
         )
         add_to_expression!(
             container,
@@ -234,7 +231,6 @@ function construct_network!(
             SystemBalanceSlackDown,
             sys,
             model,
-            T,
         )
         objective_function!(container, PSY.ACBus, model)
     end
diff --git a/src/network_models/network_slack_variables.jl b/src/network_models/network_slack_variables.jl
index 619aeab047..82e1638e49 100644
--- a/src/network_models/network_slack_variables.jl
+++ b/src/network_models/network_slack_variables.jl
@@ -37,7 +37,13 @@ function add_variables!(
     U <: PM.AbstractActivePowerModel,
 }
     time_steps = get_time_steps(container)
-    bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys))
+    radial_network_reduction = get_radial_network_reduction(network_model)
+    if isempty(radial_network_reduction)
+        bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys))
+    else
+        bus_numbers = collect(keys(PNM.get_bus_reduction_map(radial_network_reduction)))
+    end
+
     variable = add_variable_container!(container, T(), PSY.ACBus, bus_numbers, time_steps)
     for t in time_steps, n in bus_numbers
         variable[n, t] = JuMP.@variable(
@@ -59,7 +65,12 @@ function add_variables!(
     U <: PM.AbstractPowerModel,
 }
     time_steps = get_time_steps(container)
-    bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys))
+    radial_network_reduction = get_radial_network_reduction(network_model)
+    if isempty(radial_network_reduction)
+        bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys))
+    else
+        bus_numbers = collect(keys(PNM.get_bus_reduction_map(radial_network_reduction)))
+    end
     variable_active =
         add_variable_container!(container, T(), PSY.ACBus, "P", bus_numbers, time_steps)
     variable_reactive =
diff --git a/src/network_models/pm_translator.jl b/src/network_models/pm_translator.jl
index 627a5b7cf0..ee86e0544b 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_network_reduction = get_radial_network_reduction(network_model)
+    radial_branches_names = PNM.get_radial_branches(radial_network_reduction)
     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 it is 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..ea0b6bb991 100644
--- a/src/network_models/powermodels_interface.jl
+++ b/src/network_models/powermodels_interface.jl
@@ -6,6 +6,14 @@
 #################################################################################
 # Model Definitions
 
+const UNSUPPORTED_POWERMODELS =
+    [PM.SOCBFPowerModel, PM.SOCBFConicPowerModel, PM.IVRPowerModel]
+
+const INCOMPATIBLE_WITH_RADIAL_BRANCHES_POWERMODELS = [
+    PM.SparseSDPWRMPowerModel,
+    PTDFPowerModel,
+]
+
 function instantiate_nip_expr_model(data::Dict{String, Any}, model_constructor; kwargs...)
     return PM.instantiate_model(data, model_constructor, instantiate_nip_expr; kwargs...)
 end
@@ -13,7 +21,6 @@ end
 # replicates PM.build_mn_opf
 function instantiate_nip_expr(pm::PM.AbstractPowerModel)
     for n in eachindex(PM.nws(pm))
-        @assert !PM.ismulticonductor(pm; nw = n)
         PM.variable_bus_voltage(pm; nw = n)
         PM.variable_branch_power(pm; nw = n, bounded = false)
         PM.variable_dcline_power(pm; nw = n, bounded = false)
@@ -60,7 +67,6 @@ end
 # replicates PM.build_opf_ptdf
 function instantiate_nip_ptdf_expr(pm::PM.AbstractPowerModel)
     for n in eachindex(PM.nws(pm))
-        @assert !PM.ismulticonductor(pm; nw = n)
 
         #PM.variable_gen_power(pm) #connect P__* with these
 
@@ -104,7 +110,6 @@ end
 # replicates PM.build_mn_opf_bf_strg
 function instantiate_bfp_expr(pm::PM.AbstractPowerModel)
     for n in eachindex(PM.nws(pm))
-        @assert !PM.ismulticonductor(pm; nw = n)
         PM.variable_bus_voltage(pm; nw = n)
         PM.variable_branch_power(pm; nw = n, bounded = false)
         PM.variable_dcline_power(pm; nw = n, bounded = false)
@@ -291,17 +296,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_network_reduction = get_radial_network_reduction(network_model)
+    if isempty(radial_network_reduction)
+        ac_bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys))
+    else
+        bus_reduction_map = PNM.get_bus_reduction_map(radial_network_reduction)
+        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 +337,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_network_reduction = get_radial_network_reduction(network_model)
+    if isempty(radial_network_reduction)
+        ac_bus_numbers = PSY.get_number.(get_available_components(PSY.ACBus, sys))
+    else
+        bus_reduction_map = PNM.get_bus_reduction_map(radial_network_reduction)
+        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
diff --git a/src/operation/decision_model.jl b/src/operation/decision_model.jl
index 8287346b53..a0cde57e3d 100644
--- a/src/operation/decision_model.jl
+++ b/src/operation/decision_model.jl
@@ -9,6 +9,15 @@ Generic PowerSimulations Operation Problem Type for unspecified models
 """
 struct GenericOpProblem <: DefaultDecisionProblem end
 
+mutable struct DecisionModel{M <: DecisionProblem} <: OperationModel
+    name::Symbol
+    template::ProblemTemplate
+    sys::PSY.System
+    internal::Union{Nothing, ModelInternal}
+    store::DecisionModelStore
+    ext::Dict{String, Any}
+end
+
 """
     DecisionModel{M}(
         template::ProblemTemplate,
@@ -51,40 +60,31 @@ template = ProblemTemplate(CopperPlatePowerModel, devices, branches, services)
 OpModel = DecisionModel(MockOperationProblem, template, system)
 ```
 """
-mutable struct DecisionModel{M <: DecisionProblem} <: OperationModel
-    name::Symbol
-    template::ProblemTemplate
-    sys::PSY.System
-    internal::Union{Nothing, ModelInternal}
-    store::DecisionModelStore
-    ext::Dict{String, Any}
-
-    function DecisionModel{M}(
-        template::ProblemTemplate,
-        sys::PSY.System,
-        settings::Settings,
-        jump_model::Union{Nothing, JuMP.Model} = nothing;
-        name = nothing,
-    ) where {M <: DecisionProblem}
-        if name === nothing
-            name = nameof(M)
-        elseif name isa String
-            name = Symbol(name)
-        end
-        internal = ModelInternal(
-            OptimizationContainer(sys, settings, jump_model, PSY.Deterministic),
-        )
-        template_ = deepcopy(template)
-        finalize_template!(template_, sys)
-        return new{M}(
-            name,
-            template_,
-            sys,
-            internal,
-            DecisionModelStore(),
-            Dict{String, Any}(),
-        )
+function DecisionModel{M}(
+    template::ProblemTemplate,
+    sys::PSY.System,
+    settings::Settings,
+    jump_model::Union{Nothing, JuMP.Model} = nothing;
+    name = nothing,
+) where {M <: DecisionProblem}
+    if name === nothing
+        name = nameof(M)
+    elseif name isa String
+        name = Symbol(name)
     end
+    internal = ModelInternal(
+        OptimizationContainer(sys, settings, jump_model, PSY.Deterministic),
+    )
+    template_ = deepcopy(template)
+    finalize_template!(template_, sys)
+    return DecisionModel{M}(
+        name,
+        template_,
+        sys,
+        internal,
+        DecisionModelStore(),
+        Dict{String, Any}(),
+    )
 end
 
 function DecisionModel{M}(
diff --git a/src/operation/decision_model_store.jl b/src/operation/decision_model_store.jl
index b8c2591164..686a5b9586 100644
--- a/src/operation/decision_model_store.jl
+++ b/src/operation/decision_model_store.jl
@@ -27,7 +27,7 @@ end
 
 function initialize_storage!(
     store::DecisionModelStore,
-    container::OptimizationContainer,
+    container::AbstractModelContainer,
     params::ModelStoreParams,
 )
     num_of_executions = get_num_executions(params)
diff --git a/src/operation/model_internal.jl b/src/operation/model_internal.jl
index f00b690259..4753e6f648 100644
--- a/src/operation/model_internal.jl
+++ b/src/operation/model_internal.jl
@@ -12,9 +12,9 @@ mutable struct SimulationInfo
     sequence_uuid::Base.UUID
 end
 
-mutable struct ModelInternal
-    container::OptimizationContainer
-    ic_model_container::Union{Nothing, OptimizationContainer}
+mutable struct ModelInternal{T <: AbstractModelContainer}
+    container::T
+    ic_model_container::Union{Nothing, T}
     status::BuildStatus
     run_status::RunStatus
     base_conversion::Bool
@@ -31,11 +31,11 @@ mutable struct ModelInternal
 end
 
 function ModelInternal(
-    container::OptimizationContainer;
+    container::T;
     ext = Dict{String, Any}(),
     recorders = [],
-)
-    return ModelInternal(
+) where {T <: AbstractModelContainer}
+    return ModelInternal{T}(
         container,
         nothing,
         BuildStatus.EMPTY,
diff --git a/src/operation/operation_model_interface.jl b/src/operation/operation_model_interface.jl
index 688f6013f7..faa323ae6b 100644
--- a/src/operation/operation_model_interface.jl
+++ b/src/operation/operation_model_interface.jl
@@ -10,7 +10,7 @@ get_execution_count(model::OperationModel) = get_internal(model).execution_count
 get_executions(model::OperationModel) = get_internal(model).executions
 get_initial_time(model::OperationModel) = get_initial_time(get_settings(model))
 get_internal(model::OperationModel) = model.internal
-get_jump_model(model::OperationModel) = get_internal(model).container.JuMPmodel
+get_jump_model(model::OperationModel) = get_jump_model(get_internal(model).container)
 get_name(model::OperationModel) = model.name
 get_store(model::OperationModel) = model.store
 is_synchronized(model::OperationModel) = is_synchronized(get_optimization_container(model))
@@ -415,7 +415,7 @@ end
 
 function serialize_optimization_model(model::OperationModel, save_path::String)
     serialize_jump_optimization_model(
-        get_optimization_container(model),
+        get_jump_model(get_optimization_container(model)),
         save_path,
     )
     return
diff --git a/test/test_device_hvdc.jl b/test/test_device_hvdc.jl
index 42ce383876..2db47195ba 100644
--- a/test/test_device_hvdc.jl
+++ b/test/test_device_hvdc.jl
@@ -15,7 +15,7 @@
     set_device_model!(template_uc, DeviceModel(TModelHVDCLine, LossLessLine))
     model = DecisionModel(template_uc, sys_5; name = "UC", optimizer = HiGHS_optimizer)
     @test build!(model; output_dir = mktempdir()) == PSI.BuildStatus.BUILT
-    moi_tests(model, 1656, 0, 1536, 816, 888, true)
+    moi_tests(model, 1656, 288, 1248, 528, 888, true)
     @test solve!(model) == RunStatus.SUCCESSFUL
 
     template_uc = ProblemTemplate(NetworkModel(
diff --git a/test/test_network_constructors.jl b/test/test_network_constructors.jl
index 4de655700b..519b19d28f 100644
--- a/test/test_network_constructors.jl
+++ b/test/test_network_constructors.jl
@@ -1,29 +1,30 @@
 # Note to devs. Use GLPK or Cbc for models with linear constraints and linear cost functions
 # Use OSQP for models with quadratic cost function and linear constraints and ipopt otherwise
+const NETWORKS_FOR_TESTING = [
+    (PM.ACPPowerModel, fast_ipopt_optimizer),
+    (PM.ACRPowerModel, fast_ipopt_optimizer),
+    (PM.ACTPowerModel, fast_ipopt_optimizer),
+    #(PM.IVRPowerModel, fast_ipopt_optimizer), #instantiate_ivp_expr_model not implemented
+    (PM.DCPPowerModel, fast_ipopt_optimizer),
+    (PM.DCMPPowerModel, fast_ipopt_optimizer),
+    (PM.NFAPowerModel, fast_ipopt_optimizer),
+    (PM.DCPLLPowerModel, fast_ipopt_optimizer),
+    (PM.LPACCPowerModel, fast_ipopt_optimizer),
+    (PM.SOCWRPowerModel, fast_ipopt_optimizer),
+    (PM.SOCWRConicPowerModel, scs_solver),
+    (PM.QCRMPowerModel, fast_ipopt_optimizer),
+    (PM.QCLSPowerModel, fast_ipopt_optimizer),
+    #(PM.SOCBFPowerModel, fast_ipopt_optimizer), # not implemented
+    (PM.BFAPowerModel, fast_ipopt_optimizer),
+    #(PM.SOCBFConicPowerModel, fast_ipopt_optimizer), # not implemented
+    (PM.SDPWRMPowerModel, scs_solver),
+    (PM.SparseSDPWRMPowerModel, scs_solver),
+    (PTDFPowerModel, fast_ipopt_optimizer),
+]
+
 @testset "All PowerModels models construction" begin
-    networks = [
-        (PM.ACPPowerModel, fast_ipopt_optimizer),
-        (PM.ACRPowerModel, fast_ipopt_optimizer),
-        (PM.ACTPowerModel, fast_ipopt_optimizer),
-        #(PM.IVRPowerModel, fast_ipopt_optimizer), #instantiate_ivp_expr_model not implemented
-        (PM.DCPPowerModel, fast_ipopt_optimizer),
-        (PM.DCMPPowerModel, fast_ipopt_optimizer),
-        (PM.NFAPowerModel, fast_ipopt_optimizer),
-        (PM.DCPLLPowerModel, fast_ipopt_optimizer),
-        (PM.LPACCPowerModel, fast_ipopt_optimizer),
-        (PM.SOCWRPowerModel, fast_ipopt_optimizer),
-        (PM.SOCWRConicPowerModel, scs_solver),
-        (PM.QCRMPowerModel, fast_ipopt_optimizer),
-        (PM.QCLSPowerModel, fast_ipopt_optimizer),
-        #(PM.SOCBFPowerModel, fast_ipopt_optimizer), # not implemented
-        (PM.BFAPowerModel, fast_ipopt_optimizer),
-        #(PM.SOCBFConicPowerModel, fast_ipopt_optimizer), # not implemented
-        (PM.SDPWRMPowerModel, scs_solver),
-        (PM.SparseSDPWRMPowerModel, scs_solver),
-        (PTDFPowerModel, fast_ipopt_optimizer),
-    ]
     c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
-    for (network, solver) in networks
+    for (network, solver) in NETWORKS_FOR_TESTING
         template = get_thermal_dispatch_template_network(
             NetworkModel(network; PTDF_matrix = PTDF(c_sys5)),
         )
@@ -265,9 +266,9 @@ end
         PSI.ConstraintKey(PSI.NodalBalanceActiveConstraint, PSY.ACBus),
     ]
     test_results = IdDict{System, Vector{Int}}(
-        c_sys5 => [384, 0, 408, 408, 288],
-        c_sys14 => [936, 0, 1080, 1080, 840],
-        c_sys14_dc => [984, 0, 1080, 984, 840],
+        c_sys5 => [384, 144, 264, 264, 288],
+        c_sys14 => [936, 480, 600, 600, 840],
+        c_sys14_dc => [984, 432, 648, 552, 840],
     )
     test_obj_values = IdDict{System, Float64}(
         c_sys5 => 342000.0,
@@ -313,9 +314,9 @@ end
         PSI.ConstraintKey(PSI.NodalBalanceReactiveConstraint, PSY.ACBus),
     ]
     test_results = IdDict{System, Vector{Int}}(
-        c_sys5 => [1056, 0, 384, 384, 264],
-        c_sys14 => [2832, 0, 720, 720, 696],
-        c_sys14_dc => [2832, 0, 768, 672, 744],
+        c_sys5 => [1056, 144, 240, 240, 264],
+        c_sys14 => [2832, 480, 240, 240, 696],
+        c_sys14_dc => [2832, 432, 336, 240, 744],
     )
     test_obj_values = IdDict{System, Float64}(
         c_sys5 => 340000.0,
@@ -409,9 +410,9 @@ end
         c_sys14_dc => [2832, 0, 336, 240, 744],
     )
     ACT_test_results = Dict{System, Vector{Int}}(
-        c_sys5 => [1344, 0, 384, 384, 840],
-        c_sys14 => [3792, 0, 720, 720, 2616],
-        c_sys14_dc => [3696, 0, 768, 672, 2472],
+        c_sys5 => [1344, 144, 240, 240, 840],
+        c_sys14 => [3792, 480, 240, 240, 2616],
+        c_sys14_dc => [3696, 432, 336, 240, 2472],
     )
     test_results = Dict(zip(networks, [ACR_test_results, ACT_test_results]))
     for network in networks, sys in systems
@@ -448,14 +449,14 @@ end
         c_sys14_dc => 142000.0,
     )
     DCPLL_test_results = Dict{System, Vector{Int}}(
-        c_sys5 => [528, 0, 408, 408, 288],
-        c_sys14 => [1416, 0, 1080, 1080, 840],
-        c_sys14_dc => [1416, 0, 1080, 984, 840],
+        c_sys5 => [528, 144, 264, 264, 288],
+        c_sys14 => [1416, 480, 600, 600, 840],
+        c_sys14_dc => [1416, 432, 648, 552, 840],
     )
     LPACC_test_results = Dict{System, Vector{Int}}(
-        c_sys5 => [1200, 0, 384, 384, 840],
-        c_sys14 => [3312, 0, 720, 720, 2616],
-        c_sys14_dc => [3264, 0, 768, 672, 2472],
+        c_sys5 => [1200, 144, 240, 240, 840],
+        c_sys14 => [3312, 480, 240, 240, 2616],
+        c_sys14_dc => [3264, 432, 336, 240, 2472],
     )
     test_results = Dict(zip(networks, [DCPLL_test_results, LPACC_test_results]))
     for network in networks, (ix, sys) in enumerate(systems)
@@ -696,3 +697,86 @@ end
     )
     @test all(isapprox.(sum(zone_2_gen .+ zone_2_load; dims = 2), 0.0; atol = 1e-3))
 end
+
+# These models are easier to test due to their lossless nature
+@testset "StandardPTDF/DCPPowerModel Radial Branches Test" begin
+    new_sys = PSB.build_system(PSITestSystems, "c_sys5_radial")
+
+    for net_model in [DCPPowerModel, StandardPTDFModel]
+        template_uc = template_unit_commitment(;
+            network = NetworkModel(net_model;
+                reduce_radial_branches = true,
+                use_slacks = false,
+            ),
+        )
+        thermal_model = ThermalStandardUnitCommitment
+        set_device_model!(template_uc, ThermalStandard, thermal_model)
+
+        ##### Solve Reduced Model ####
+        solver = GLPK_optimizer
+        uc_model_red = DecisionModel(
+            template_uc,
+            new_sys;
+            optimizer = solver,
+            name = "UC_RED",
+            store_variable_names = true,
+        )
+
+        @test build!(uc_model_red; output_dir = mktempdir(; cleanup = true)) ==
+              PSI.BuildStatus.BUILT
+        solve!(uc_model_red)
+
+        res_red = ProblemResults(uc_model_red)
+
+        flow_lines = read_variable(res_red, "FlowActivePowerVariable__Line")
+        line_names = DataFrames.names(flow_lines)[2:end]
+
+        ##### Solve Original Model ####
+        template_uc_orig = template_unit_commitment(;
+            network = NetworkModel(net_model;
+                reduce_radial_branches = false,
+                use_slacks = false,
+            ),
+        )
+        set_device_model!(template_uc_orig, ThermalStandard, thermal_model)
+
+        uc_model_orig = DecisionModel(
+            template_uc_orig,
+            new_sys;
+            optimizer = solver,
+            name = "UC_ORIG",
+            store_variable_names = true,
+        )
+
+        @test build!(uc_model_orig; output_dir = mktempdir(; cleanup = true)) ==
+              PSI.BuildStatus.BUILT
+        solve!(uc_model_orig)
+
+        res_orig = ProblemResults(uc_model_orig)
+
+        flow_lines_orig = read_variable(res_orig, "FlowActivePowerVariable__Line")
+
+        for line in line_names
+            @test isapprox(flow_lines[!, line], flow_lines_orig[!, line])
+        end
+    end
+end
+
+@testset "All PowerModels models construction with reduced radial branches" begin
+    new_sys = PSB.build_system(PSITestSystems, "c_sys5_radial")
+    for (network, solver) in NETWORKS_FOR_TESTING
+        if network ∈ PSI.INCOMPATIBLE_WITH_RADIAL_BRANCHES_POWERMODELS
+            continue
+        end
+        template = get_thermal_dispatch_template_network(
+            NetworkModel(network;
+                PTDF_matrix = PTDF(new_sys),
+                reduce_radial_branches = true,
+                use_slacks = true),
+        )
+        ps_model = DecisionModel(template, new_sys; optimizer = solver)
+        @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) ==
+              PSI.BuildStatus.BUILT
+        @test ps_model.internal.container.pm !== nothing
+    end
+end
diff --git a/test/test_utils/mock_operation_models.jl b/test/test_utils/mock_operation_models.jl
index a41128efb5..759684109c 100644
--- a/test/test_utils/mock_operation_models.jl
+++ b/test/test_utils/mock_operation_models.jl
@@ -117,6 +117,7 @@ function mock_construct_device!(
         PSI.get_network_formulation(template),
         PSI.get_network_model(template).subnetworks,
         PSI.get_system(problem),
+        Dict{Int64, Set{Int64}}(),
     )
     if PSI.validate_available_devices(model, PSI.get_system(problem))
         PSI.construct_device!(