diff --git a/src/core/network_model.jl b/src/core/network_model.jl
index 32a6f2be32..8c5a21abae 100644
--- a/src/core/network_model.jl
+++ b/src/core/network_model.jl
@@ -33,6 +33,7 @@ mutable struct NetworkModel{T <: PM.AbstractPowerModel}
     bus_area_map::Dict{PSY.ACBus, Int}
     duals::Vector{DataType}
     radial_branches::PNM.RadialBranches
+    bus_aggregation_map::Dict{Int, Int}
     reduce_radial_branches::Bool
 
     function NetworkModel(
@@ -51,6 +52,7 @@ mutable struct NetworkModel{T <: PM.AbstractPowerModel}
             Dict{PSY.ACBus, Int}(),
             duals,
             PNM.RadialBranches(),
+            Dict{Int, Int}(),
             reduce_radial_branches,
         )
     end
diff --git a/src/core/optimization_container.jl b/src/core/optimization_container.jl
index 89ba392b50..82e5fb88b2 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) =>
@@ -502,9 +518,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
 
@@ -519,7 +536,8 @@ function build_impl!(
         container,
         transmission,
         transmission_model.subnetworks,
-        sys,)
+        sys,
+        transmission_model.radial_branches.bus_reduction_map)
 
     # Order is required
     for device_model in values(template.devices)