diff --git a/.gitignore b/.gitignore
index a2917e1..5aab608 100644
--- a/.gitignore
+++ b/.gitignore
@@ -121,5 +121,8 @@ $RECYCLE.BIN/
 # Windows shortcuts
 *.lnk
 
+# Claude - modify to include claude files 
+.claude/
+
 ## Acknowledgements
 # Many thanks to `https://gitignore.io/`, written and maintained by Joe Blau, which contributed much material to this gitignore file.
diff --git a/Project.toml b/Project.toml
index bb3e04b..05902ef 100644
--- a/Project.toml
+++ b/Project.toml
@@ -16,14 +16,14 @@ PowerSimulations = "e690365d-45e2-57bb-ac84-44ba829e73c4"
 PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
 
 [compat]
-DataStructures = "^0.18"
+DataStructures = "~0.18, ~0.19"
 Dates = "1"
 DocStringExtensions = "~0.8, ~0.9"
-InfrastructureSystems = "2"
+InfrastructureSystems = "3"
 JuMP = "1"
 MPI = "^0.20"
 MathOptInterface = "1"
-PowerNetworkMatrices = "0.12"
-PowerSystems = "4"
-PowerSimulations = "^0.30"
+PowerNetworkMatrices = "^0.20"
+PowerSystems = "5.3"
+PowerSimulations = "0.34"
 julia = "^1.6"
diff --git a/src/algorithms/sequential_algorithm.jl b/src/algorithms/sequential_algorithm.jl
index 9587728..22368a4 100644
--- a/src/algorithms/sequential_algorithm.jl
+++ b/src/algorithms/sequential_algorithm.jl
@@ -20,21 +20,14 @@ function build_main_problem!(
     template::MultiProblemTemplate,
     sys::PSY.System,
 )
-    branch_models_dict = keys(PSI.get_branch_models(template))
     device_models_dict = keys(PSI.get_device_models(template))
-    if :TwoTerminalHVDCLine ∈ branch_models_dict ||
-       :TwoTerminalVSCDCLine ∈ branch_models_dict
-        has_hvdc = true
-    else
-        has_hvdc = false
-    end
     for k in keys(container.subproblems)
-        subsystem_buses = PSY.get_components(PSY.ACBus, sys; subsystem_name=k)
+        subsystem_buses = PSY.get_available_components(PSY.ACBus, sys; subsystem_name=k)
         subsystem_bus_nos = [PSY.get_number(b) for b in subsystem_buses]
         container.subproblem_bus_map[k] = subsystem_bus_nos
         subsystem_static_injectors =
-            PSY.get_components(PSY.StaticInjection, sys; subsystem_name=k)
-        subsystem_buses = PSY.get_components(PSY.ACBus, sys; subsystem_name=k)
+            PSY.get_available_components(PSY.StaticInjection, sys; subsystem_name=k)
+        subsystem_buses = PSY.get_available_components(PSY.ACBus, sys; subsystem_name=k)
         for si in subsystem_static_injectors
             if Symbol(typeof(si)) ∈ device_models_dict
                 if !(PSY.get_bus(si) ∈ subsystem_buses)
@@ -46,12 +39,9 @@ function build_main_problem!(
                 end
             end
         end
-        subsystem_hvdcs = PSY.get_components(
-            Union{PSY.TwoTerminalHVDCLine, PSY.TwoTerminalVSCDCLine},
-            sys;
-            subsystem_name=k,
-        )
-        if has_hvdc
+        subsystem_hvdcs =
+            PSY.get_available_components(PSY.TwoTerminalHVDC, sys; subsystem_name=k)
+        if _has_hvdc_model(template)
             for hvdc in subsystem_hvdcs
                 from_bus_no = PSY.get_number(PSY.get_from(PSY.get_arc(hvdc)))
                 to_bus_no = PSY.get_number(PSY.get_to(PSY.get_arc(hvdc)))
@@ -67,6 +57,20 @@ function build_main_problem!(
     end
 end
 
+function _has_hvdc_model(template::MultiProblemTemplate)
+    branch_models_dict = keys(PSI.get_branch_models(template))
+    for hvdc_type in CONCRETE_HVDC_TYPES
+        if hvdc_type in branch_models_dict
+            return true
+        end
+    end
+    return false
+end
+
+# Note: With the addition of 3D results processing, we can eliminate this design of writing
+# all subsystem results to a signle container. This way we avoid overwriting the results from 
+# different subsystem problems and can still access all relevant results using the API. 
+
 # The drawback of this approach is that it will loop over the results twice
 # once to write into the main container and a second time when writing into the
 # store. The upside of this approach is that doesn't require overloading write_model_XXX_results!
@@ -79,23 +83,30 @@ function write_results_to_main_container(container::MultiOptimizationContainer)
             subproblem_data_field = getproperty(subproblem, field)
             main_container_data_field = getproperty(container, field)
             for (key, src) in subproblem_data_field
+                @debug "writing $key to main container"
                 if src isa JuMP.Containers.SparseAxisArray
                     @debug "Skip SparseAxisArray" field key
                     continue
                 end
                 num_dims = ndims(src)
-                num_dims > 2 && error("ndims = $(num_dims) is not supported yet")
-                data = nothing
-                data = PSI.jump_value.(src)
                 dst = main_container_data_field[key]
                 if num_dims == 1
-                    dst[1:length(axes(src)[1])] = data
+                    columns = _get_main_container_columns(container, k, key, src)
+                    for col in columns
+                        if isassigned(src, col)
+                            dst[col] = PSI.jump_value.(src[col])
+                        end
+                    end
                 elseif num_dims == 2
                     columns = _get_main_container_columns(container, k, key, src)
-                    len = length(axes(src)[2])
-                    dst[columns, 1:len] = PSI.jump_value.(src[columns, :])
+                    for col in columns
+                        for t in axes(src)[2]
+                            if isassigned(src, col, t)
+                                dst[col, t] = PSI.jump_value.(src[col, t])
+                            end
+                        end
+                    end
                 elseif num_dims == 3
-                    # TODO: untested
                     axis1 = axes(src)[1]
                     axis2 = axes(src)[2]
                     len = length(axes(src)[3])
diff --git a/src/core/definitions.jl b/src/core/definitions.jl
index a7d701b..1018903 100644
--- a/src/core/definitions.jl
+++ b/src/core/definitions.jl
@@ -1,3 +1,6 @@
 const CONTAINER_FIELDS = [:variables, :aux_variables, :constraints, :expressions, :duals]
 const ALL_CONTAINER_FIELDS =
     [:variables, :aux_variables, :constraints, :expressions, :duals, :parameters]
+
+const CONCRETE_HVDC_TYPES =
+    [:TwoTerminalGenericHVDCLine, :TwoTerminalLCCLine, :TwoTerminalGenericHVDCLine]
diff --git a/src/models/branch_models.jl b/src/models/branch_models.jl
index 210798d..34d8a0d 100644
--- a/src/models/branch_models.jl
+++ b/src/models/branch_models.jl
@@ -6,13 +6,6 @@ function PSI.construct_device!(
     network_model::PSI.NetworkModel{<:PSI.AbstractPTDFModel},
 )
     devices = PSI.get_available_components(model, sys)
-    PSI.add_variables!(
-        container,
-        PSI.FlowActivePowerVariable,
-        network_model,
-        devices,
-        PSI.StaticBranchUnbounded(),
-    )
     PSI.add_parameters!(container, StateEstimationFlows, devices, model)
     PSI.add_feedforward_arguments!(container, model, devices)
     return
@@ -30,7 +23,7 @@ function PSI.add_parameters!(
         container,
         StateEstimationFlows(),
         D,
-        ISOPT.VariableKey{PSI.FlowActivePowerVariable, D}(""),
+        ISOPT.ExpressionKey{PSI.PTDFBranchFlow, D}(""),
         branch_names,
         time_steps,
     )
@@ -54,7 +47,7 @@ function PSI._make_flow_expressions!(
 )
     # NOTE - can alternatively avoid passing the system by ensuring the parameter container for state estimation injections
     # is ordered by bus number as is the case for the Active Power Balance expressions
-    all_buses = PSY.get_components(
+    all_buses = PSY.get_available_components(
         x -> PSY.get_bustype(x) != PSY.ACBusTypes.ISOLATED,
         PSY.ACBus,
         sys;
@@ -118,20 +111,14 @@ function PSI._make_flow_expressions!(
     container::PSI.OptimizationContainer,
     branches::Vector{String},
     time_steps::UnitRange{Int},
-    ptdf::PSI.ValidPTDFS,
-    nodal_balance_expressions::PSI.JuMPAffineExpressionDArray,
+    ptdf::Union{PNM.PTDF, PNM.VirtualPTDF},
+    nodal_balance_expressions::PSI.JuMPAffineExpressionDArrayIntInt,
     state_estimation_injections,
     state_estimation_flows,
     branch_Type::DataType,
     sys::PSY.System,
 )
-    branch_flow_expr = PSI.add_expression_container!(
-        container,
-        PSI.PTDFBranchFlow(),
-        branch_Type,
-        branches,
-        time_steps,
-    )
+    branch_flow_expr = PSI.get_expression(container, PSI.PTDFBranchFlow(), branch_Type)
 
     jump_model = PSI.get_jump_model(container)
 
@@ -183,8 +170,8 @@ function PSI.add_constraints!(
     ptdf = PSI.get_PTDF_matrix(network_model)
     # This is a workaround to not call the same list comprehension to find
     # The subset of branches of type B in the PTDF
-    flow_variables = PSI.get_variable(container, PSI.FlowActivePowerVariable(), B)
-    branches = flow_variables.axes[1]
+    flows = PSI.get_expression(container, PSI.PTDFBranchFlow(), B)
+    branches = flows.axes[1]
     time_steps = PSI.get_time_steps(container)
     branch_flow = PSI.add_constraints_container!(
         container,
@@ -199,7 +186,6 @@ function PSI.add_constraints!(
         PSI.get_parameter_array(container, StateEstimationInjections(), PSY.ACBus)
 
     state_estimation_flows = PSI.get_parameter(container, StateEstimationFlows(), B)
-    flow_variables = PSI.get_variable(container, PSI.FlowActivePowerVariable(), B)
     branch_flow_expr = PSI._make_flow_expressions!(
         container,
         branches,
@@ -216,7 +202,7 @@ function PSI.add_constraints!(
         for t in time_steps
             branch_flow[name, t] = JuMP.@constraint(
                 jump_model,
-                branch_flow_expr[name, t] - flow_variables[name, t] == 0.0
+                branch_flow_expr[name, t] - flows[name, t] == 0.0
             )
         end
     end
@@ -231,6 +217,7 @@ function PSI.construct_device!(
     network_model::PSI.NetworkModel{<:PSI.AbstractPTDFModel},
 )
     devices = PSI.get_available_components(model, sys)
+    PSI.add_expressions!(container, PSI.PTDFBranchFlow, devices, model, network_model)
     # NOTE - changes required in handling of feedforwards? 
     PSI.add_feedforward_constraints!(container, model, devices)
     PSI.add_constraints!(
diff --git a/src/models/network_models.jl b/src/models/network_models.jl
index e61b4bf..287fa01 100644
--- a/src/models/network_models.jl
+++ b/src/models/network_models.jl
@@ -42,7 +42,6 @@ function PSI.construct_network!(
         model,
     )
     PSI.add_constraints!(container, PSI.CopperPlateBalanceConstraint, sys, model)
-    PSI.add_constraints!(container, PSI.NodalBalanceActiveConstraint, sys, model)
     PSI.add_constraint_dual!(container, sys, model)
     return
 end
@@ -79,7 +78,7 @@ function PSI.add_to_expression!(
     parameter_array =
         PSI.get_parameter_array(container, StateEstimationInjections(), PSY.ACBus)
     subsys = PSI.get_subsystem(network_model)
-    all_buses = PSY.get_components(
+    all_buses = PSY.get_available_components(
         x -> PSY.get_bustype(x) != PSY.ACBusTypes.ISOLATED,
         PSY.ACBus,
         sys;
@@ -88,12 +87,12 @@ function PSI.add_to_expression!(
     # These are the buses not in the same subsystem as the one being built
 
     expression = PSI.get_expression(container, PSI.ActivePowerBalance(), PSY.ACBus)
-    radial_network_reduction = PSI.get_radial_network_reduction(network_model)
+    nrd = PSI.get_network_reduction(network_model)
     for b in all_buses, t in PSI.get_time_steps(container)
         if PSY.has_component(sys, subsys, b)
             continue
         end
-        bus_no = PNM.get_mapped_bus_number(radial_network_reduction, b)
+        bus_no = PNM.get_mapped_bus_number(nrd, b)
         PSI._add_to_jump_expression!(
             expression[bus_no, t],
             parameter_array[string(bus_no), t],
@@ -130,13 +129,7 @@ function PSI.add_parameters!(
     time_steps = PSI.get_time_steps(container)
     subsys = PSI.get_subsystem(network_model)
 
-    all_buses = PSY.get_components(
-        x -> PSY.get_bustype(x) != PSY.ACBusTypes.ISOLATED,
-        PSY.ACBus,
-        sys;
-    )
-
-    bus_numbers = [string(PSY.get_number(b)) for b in all_buses]
+    bus_numbers = string.(PNM.get_bus_axis(PSI.get_PTDF_matrix(network_model)))
     @assert !isempty(bus_numbers)
 
     parameter_container = PSI.add_param_container!(
@@ -162,6 +155,7 @@ function PSI.initialize_system_expressions!(
     container::PSI.OptimizationContainer,
     network_model::PSI.NetworkModel{SplitAreaPTDFPowerModel},
     subnetworks::Dict{Int, Set{Int}},
+    ::PSI.BranchModelContainer,
     system::PSY.System,
     bus_reduction_map::Dict{Int64, Set{Int64}},
 )
@@ -198,12 +192,12 @@ function PSI.add_to_expression!(
     variable = PSI.get_variable(container, U(), V)
     area_expr = PSI.get_expression(container, T(), PSY.Area)
     nodal_expr = PSI.get_expression(container, T(), PSY.ACBus)
-    radial_network_reduction = PSI.get_radial_network_reduction(network_model)
+    nrd = PSI.get_network_reduction(network_model)
     for d in devices
         name = PSY.get_name(d)
         device_bus = PSY.get_bus(d)
         area_name = PSY.get_name(PSY.get_area(device_bus))
-        bus_no = PNM.get_mapped_bus_number(radial_network_reduction, device_bus)
+        bus_no = PNM.get_mapped_bus_number(nrd, device_bus)
         for t in PSI.get_time_steps(container)
             PSI._add_to_jump_expression!(
                 area_expr[area_name, t],
@@ -281,19 +275,17 @@ function _update_parameter_values!(
 ) where {T <: Union{JuMP.VariableRef, Float64}}
     state = PSI.get_system_states(simulation_state)
     state_values = PSI.get_dataset_values(state, PSI.get_attribute_key(attributes))
-
     if !isfinite(first(state_values))
         @warn "first value not present, updating state estimation injections from decision state"
         state = PSI.get_decision_states(simulation_state)
         state_values = PSI.get_dataset_values(state, PSI.get_attribute_key(attributes))
+    elseif size(parameter_array)[2] > size(state_values)[2]
+        @warn "Cannot update; state estimation injection parameter has more timesteps than the system state, updating state estimation injections from decision state"
+        state = PSI.get_decision_states(simulation_state)
+        state_values = PSI.get_dataset_values(state, PSI.get_attribute_key(attributes))
+    else
+        @info "Updating state estimation injection parameters from system state"
     end
-
-    if size(parameter_array)[2] > size(state_values)[2]
-        error(
-            "Cannot update: state estimation injection parameter has more timesteps than the state used for updating.",
-        )
-    end
-
     component_names, time = axes(parameter_array)
     for t in time
         for name in component_names
diff --git a/src/multi_optimization_container.jl b/src/multi_optimization_container.jl
index 6362ec3..343502b 100644
--- a/src/multi_optimization_container.jl
+++ b/src/multi_optimization_container.jl
@@ -6,7 +6,6 @@ Base.@kwdef mutable struct MultiOptimizationContainer{T <: DecompositionAlgorith
     time_steps::UnitRange{Int}
     resolution::Dates.TimePeriod
     settings::PSI.Settings
-    settings_copy::PSI.Settings
     variables::Dict{ISOPT.VariableKey, AbstractArray}
     aux_variables::Dict{ISOPT.AuxVarKey, AbstractArray}
     duals::Dict{ISOPT.ConstraintKey, AbstractArray}
@@ -59,7 +58,6 @@ function MultiOptimizationContainer(
         time_steps=1:1,
         resolution=IS.time_period_conversion(resolution),
         settings=settings,
-        settings_copy=PSI.copy_for_serialization(settings),
         variables=Dict{ISOPT.VariableKey, AbstractArray}(),
         aux_variables=Dict{ISOPT.AuxVarKey, AbstractArray}(),
         duals=Dict{ISOPT.ConstraintKey, AbstractArray}(),
@@ -187,3 +185,18 @@ function PSI.serialize_optimization_model(
     container::MultiOptimizationContainer,
     save_path::String,
 ) end
+
+function PSI.get_column_names(
+    ::MultiOptimizationContainer,
+    field::Symbol,
+    subcontainer,
+    key::PSI.OptimizationContainerKey,
+)
+    return if field == :parameters
+        # Parameters are stored in ParameterContainer.
+        PSI.get_column_names(key, subcontainer)
+    else
+        # The others are in DenseAxisArrays.
+        PSI.get_column_names_from_axis_array(key, subcontainer)
+    end
+end
diff --git a/src/multiproblem_template.jl b/src/multiproblem_template.jl
index 4113930..afa37c2 100644
--- a/src/multiproblem_template.jl
+++ b/src/multiproblem_template.jl
@@ -215,11 +215,79 @@ function PSI.set_service_model!(
     return
 end
 
+# Set service model for only one subsystem:
+function PSI.set_service_model!(
+    template::MultiProblemTemplate,
+    model::PSI.ServiceModel{<:PSY.Service, <:PSI.AbstractServiceFormulation},
+    subsystem_id::String,
+)
+    PSI.set_service_model!(template.base_template, model)
+    for (id, sub_template) in get_sub_templates(template)
+        if id == subsystem_id
+            PSI.set_subsystem!(model, id)
+            PSI.set_service_model!(sub_template, deepcopy(model))
+        end
+    end
+    return
+end
+
 function finalize_template!(template::MultiProblemTemplate, sys::PSY.System)
     PSI.finalize_template!(template.base_template, sys)
     for (ix, sub_template) in get_sub_templates(template)
-        @debug "Finalizing template for sub probem $ix"
-        PSI.finalize_template!(sub_template, sys)
+        finalize_template!(sub_template, sys, ix)
+    end
+    return
+end
+
+function finalize_template!(
+    template::PSI.ProblemTemplate,
+    sys::PSY.System,
+    subsystem::String,
+)
+    _add_modeled_ac_branches!(template, sys, subsystem)
+    _check_for_empty_device_models!(template, sys, subsystem)
+    PSI._populate_aggregated_service_model!(template, sys)
+    PSI._populate_contributing_devices!(template, sys)
+    PSI._add_services_to_device_model!(template)
+    return
+end
+
+function _check_for_empty_device_models!(
+    template::PSI.ProblemTemplate,
+    sys::PSY.System,
+    subsystem::String,
+)
+    for device_model in values(PSI.get_device_models(template))
+        component_type = PSI.get_component_type(device_model)
+        components =
+            PSY.get_available_components(component_type, sys; subsystem_name=subsystem)
+        if isempty(components)
+            pop!(PSI.get_device_models(template), Symbol(component_type))
+            @warn "Device model for $component_type is included in the main template but there are no available components of this type in subsystem $subsystem and will be removed from the template"
+        end
+    end
+    return
+end
+
+function _add_modeled_ac_branches!(
+    template::PSI.ProblemTemplate,
+    sys::PSY.System,
+    subsystem::String,
+)
+    network_model = PSI.get_network_model(template)
+    branch_models = PSI.get_branch_models(template)
+    for v in values(branch_models)
+        component_type = PSI.get_component_type(v)
+        if isempty(
+            PSY.get_available_components(component_type, sys; subsystem_name=subsystem),
+        )
+            pop!(branch_models, Symbol(component_type))
+            @warn "$component_type is modeled but not in subsystem $subsystem so removed from the template"
+        else
+            if (component_type <: PSY.ACTransmission)
+                push!(network_model.modeled_ac_branch_types, component_type)
+            end
+        end
     end
     return
 end
diff --git a/src/print.jl b/src/print.jl
index e1de52b..f9938be 100644
--- a/src/print.jl
+++ b/src/print.jl
@@ -1,7 +1,15 @@
-function Base.show(io::IO, ::MIME"text/plain", input::MultiProblemTemplate)
-    println(io, "Print somenthing clever here. Template")
+function Base.show(io::IO, x::MIME"text/plain", input::MultiProblemTemplate)
+    for (k, template) in input.sub_templates
+        println(io, "SUBTEMPLATE $k:")
+        Base.show(io::IO, x, template)
+        println(io, "\n")
+    end
 end
 
-function Base.show(io::IO, ::MIME"text/plain", input::PSI.DecisionModel{MultiRegionProblem})
-    println(io, "Print somenthing clever here. Problem")
+function Base.show(
+    io::IO,
+    x::MIME"text/plain",
+    input::PSI.DecisionModel{MultiRegionProblem},
+)
+    Base.show(io::IO, x, input.template)
 end
diff --git a/src/problems/multi_region_problem.jl b/src/problems/multi_region_problem.jl
index c651cbb..45986de 100644
--- a/src/problems/multi_region_problem.jl
+++ b/src/problems/multi_region_problem.jl
@@ -8,6 +8,7 @@ function PSI.DecisionModel{MultiRegionProblem}(
 )
     name = Symbol(get(kwargs, :name, nameof(MultiRegionProblem)))
     settings = PSI.Settings(sys; [k for k in kwargs if first(k) ∉ [:name]]...)
+    PSI.auto_transform_time_series!(sys, settings)
     internal = ISOPT.ModelInternal(
         MultiOptimizationContainer(
             SequentialAlgorithm,
@@ -18,7 +19,6 @@ function PSI.DecisionModel{MultiRegionProblem}(
         ),
     )
     template_ = deepcopy(template)
-
     finalize_template!(template_, sys)
 
     model = PSI.DecisionModel{MultiRegionProblem}(
@@ -104,6 +104,14 @@ function _get_axes!(
     return
 end
 
+function _make_joint_axes!(
+    dim1::Set{String},
+    dim2::Set{T},
+    dim3::Set{UnitRange{Int}},
+) where {T <: Union{Int, String}}
+    return (collect(dim1), collect(dim2), first(dim3))
+end
+
 function _make_joint_axes!(
     dim1::Set{T},
     dim2::Set{UnitRange{Int}},
@@ -116,7 +124,6 @@ function _make_joint_axes!(dim1::Set{UnitRange{Int}})
 end
 
 function _make_joint_axes!(dim1::Set{String})
-    @error dim1
     return (collect(dim1),)
 end
 
@@ -129,7 +136,6 @@ function _map_containers(model::PSI.DecisionModel{MultiRegionProblem})
     for subproblem_container in values(container.subproblems)
         _get_axes!(common_axes, subproblem_container)
     end
-
     for (field, vals) in common_axes
         field_data = getproperty(container, field)
         for (key, axes_data) in vals
@@ -276,19 +282,22 @@ function handle_initial_conditions!(model::PSI.DecisionModel{MultiRegionProblem}
 
 function instantiate_network_model(model::PSI.DecisionModel{MultiRegionProblem})
     template = PSI.get_template(model)
+    sys = PSI.get_system(model)
     for (id, sub_template) in get_sub_templates(template)
         network_model = PSI.get_network_model(sub_template)
         PSI.set_subsystem!(network_model, id)
-        PSI.instantiate_network_model(network_model, PSI.get_system(model))
+        branch_models = PSI.get_branch_models(sub_template)
+        number_of_steps = PSI.get_time_steps(PSI.get_optimization_container(model))[end]
+        for model in values(branch_models)
+            top_level_filter = get(model.attributes, "filter_function", x -> true)
+            model.attributes["filter_function"] =
+                x -> top_level_filter(x) && PSY.has_component(sys, id, x)
+        end
+        PSI.instantiate_network_model!(network_model, branch_models, number_of_steps, sys)
     end
     return
 end
 
-function PSI.serialize_problem(
-    model::PSI.DecisionModel{MultiRegionProblem};
-    optimizer::Nothing,
-) end
-
 function PSI.build_model!(model::PSI.DecisionModel{MultiRegionProblem})
     build_impl!(
         PSI.get_optimization_container(model),
diff --git a/test/Project.toml b/test/Project.toml
index 2b9808f..6ca88c2 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,9 +1,12 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 HydroPowerSimulations = "fc1677e0-6ad7-4515-bf3a-bd6bf20a0b1b"
 InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
+PowerFlows = "94fada2c-fd9a-4e89-8d82-81405f5cb4f6"
 PowerNetworkMatrices = "bed98974-b02a-5e2f-9fe0-a103f5c450dd"
 PowerSimulations = "e690365d-45e2-57bb-ac84-44ba829e73c4"
 PowerSimulationsDecomposition = "bed98974-b02a-5e2f-9ee0-a103f5c450dd"
@@ -11,4 +14,5 @@ PowerSystemCaseBuilder = "f00506e0-b84f-492a-93c2-c0a9afc4364e"
 PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 StorageSystemsSimulations = "e2f1a126-19d0-4674-9252-42b2384f8e3c"
-Xpress = "9e70acf3-d6c9-5be6-b5bd-4e2c73e3e054"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TimeSeries = "9e3dc215-6440-5c97-bce1-76c03772f85e"
diff --git a/test/first_test.jl b/test/first_test.jl
deleted file mode 100644
index b15931f..0000000
--- a/test/first_test.jl
+++ /dev/null
@@ -1,148 +0,0 @@
-# build function able to exclude elements belonging to a certain sub region
-# THIS IS A HACK!!!
-
-# TODO: add field `ext` to all the topology types (Area, LoadZone, Arc)
-
-# SIIP Packages
-using Revise
-
-using HydroPowerSimulations
-using PowerSimulations
-using PowerSystemCaseBuilder
-using InfrastructureSystems
-using StorageSystemsSimulations
-using PowerSystems
-using Xpress
-using PowerSimulationsDecomposition
-using Logging
-# using Plots
-
-# using Xpress
-# using JuMP
-# using Logging
-# using Dates
-# using DataFrames
-
-const PSY = PowerSystems
-const PSI = PowerSimulations
-const PSB = PowerSystemCaseBuilder
-
-# consider the use of custom system used for GDO case
-name_ = "AC_inter"
-sys_twin_rts_DA =
-    PSY.System("GDO systems/saved_main_RTS_GMLC_DA_final_sys_" * name_ * ".json")   # day ahead
-
-# ! check reserves
-
-# modify the system -> add features in the "ext" field
-for d in PSY.get_components(PSY.Component, sys_twin_rts_DA)
-    if typeof(d) <: PSY.Bus || :available in fieldnames(typeof(d))
-        if occursin("twin", PSY.get_name(d))
-            PSY.set_ext!(d, Dict("subregion" => Set(["2"])))
-        else
-            PSY.set_ext!(d, Dict("subregion" => Set(["1"])))
-        end
-    end
-end
-
-# interconnection is shared between the two regions
-if name_ == "AC_inter"
-    br = get_component(PSY.MonitoredLine, sys_twin_rts_DA, "AC_interconnection")
-elseif name_ == "HVDC_inter"
-    br = get_component(PSY.TwoTerminalHVDCLine, sys_twin_rts_DA, "HVDC_interconnection")
-end
-
-PSY.set_ext!(br, Dict("subregion" => Set(["1", "2"])))
-arc_ = get_arc(br)
-for b in [get_from(arc_), get_to(arc_)]
-    PSY.set_ext!(b, Dict("subregion" => Set(["1", "2"])))
-end
-
-if name_ == "HVDC_inter"
-    HVDC_inter = true
-else
-    HVDC_inter = false
-end
-
-# DC power flow reference solution
-
-# define battery model
-storage_model = DeviceModel(
-    EnergyReservoirStorage,
-    StorageDispatchWithReserves;
-    attributes=Dict(
-        "reservation" => false,
-        "cycling_limits" => false,
-        "energy_target" => false,
-        "complete_coverage" => false,
-        "regularization" => true,
-    ),
-)
-
-# UC model
-template_uc = MultiProblemTemplate(
-    # NetworkModel(StandardPTDFModel; PTDF_matrix = PTDF(sys_twin_rts)),
-    NetworkModel(DCPPowerModel; use_slacks=true),
-    ["1", "2"],
-)
-set_device_model!(template_uc, ThermalStandard, ThermalStandardUnitCommitment)
-set_device_model!(template_uc, RenewableDispatch, RenewableFullDispatch)
-set_device_model!(template_uc, RenewableFix, FixedOutput)
-set_device_model!(template_uc, PowerLoad, StaticPowerLoad)
-set_device_model!(template_uc, Line, StaticBranch)
-set_device_model!(template_uc, Transformer2W, StaticBranchUnbounded)
-set_device_model!(template_uc, TapTransformer, StaticBranchUnbounded)
-set_device_model!(template_uc, HydroDispatch, FixedOutput)
-set_device_model!(template_uc, HydroEnergyReservoir, FixedOutput)
-set_device_model!(template_uc, storage_model)
-set_service_model!(
-    template_uc,
-    ServiceModel(VariableReserve{ReserveUp}, RangeReserve; use_slacks=true),
-)
-set_service_model!(
-    template_uc,
-    ServiceModel(VariableReserve{ReserveDown}, RangeReserve; use_slacks=true),
-)
-
-# add the HVDC line in case is present
-if HVDC_inter == "true"
-    set_device_model!(template_uc, TwoTerminalHVDCLine, HVDCTwoTerminalDispatch)
-else
-    set_device_model!(template_uc, MonitoredLine, StaticBranch)
-end
-
-model = DecisionModel(
-    MultiRegionProblem,
-    template_uc,
-    sys_twin_rts_DA;
-    name="UC",
-    optimizer=optimizer_with_attributes(
-        Xpress.Optimizer,
-        "MIPRELSTOP" => 0.01,       # Set the relative mip gap tolerance
-        "MAXMEMORYSOFT" => 600000,   # Set the maximum amount of memory the solver can use (in MB)
-    ),
-    system_to_file=false,
-    initialize_model=true,
-    optimizer_solve_log_print=true,
-    direct_mode_optimizer=true,
-    rebuild_model=false,
-    store_variable_names=true,
-    calculate_conflict=true,
-)
-
-# for b in get_components(ACBus, model.sys)
-#     @show get_bustype(b)
-#     @show get_bustype(b) == PSY.ACBusTypes.ISOLATED
-#     # if get_bustype(b) == PSY.ACBusTypes.ISOLATED
-#     #     @show get_name(b)
-#     # end
-# end
-
-# b = get_component(ACBus, model.sys, "Caesar")
-# get_bustype(b) == ACBusTypes.ISOLATED
-# get_name(b)
-
-# PowerSimulationsDecomposition.instantiate_network_model(model)
-
-build!(model; console_level=Logging.Info, output_dir=mktempdir())
-solve!(model; console_level=Logging.Info)
diff --git a/test/runtests.jl b/test/runtests.jl
index 290a1a3..20003a2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,11 +9,13 @@ using HydroPowerSimulations
 import PowerSystemCaseBuilder: PSITestSystems
 using PowerNetworkMatrices
 using StorageSystemsSimulations
+using TimeSeries
 using Dates
 using HiGHS
 using JuMP
 const IS = InfrastructureSystems
 const PSI = PowerSimulations
+const PSY = PowerSystems
 
 # Test Packages
 using Test
@@ -30,13 +32,16 @@ Aqua.test_undefined_exports(PowerSimulationsDecomposition)
 Aqua.test_stale_deps(PowerSimulationsDecomposition)
 Aqua.test_deps_compat(PowerSimulationsDecomposition)
 
-LOG_FILE = "power-systems.log"
-LOG_LEVELS = Dict(
+const LOG_FILE = "power-systems.log"
+const LOG_LEVELS = Dict(
     "Debug" => Logging.Debug,
     "Info" => Logging.Info,
     "Warn" => Logging.Warn,
     "Error" => Logging.Error,
 )
+const DISABLED_TEST_FILES = [  # Can generate with ls -1 test | grep "test_.*.jl"
+    "test_reserve_deliverability.jl",
+]
 
 """
 Copied @includetests from https://github.com/ssfrr/TestSetExtensions.jl.
@@ -72,7 +77,11 @@ macro includetests(testarg...)
             tests = map(f -> string(f, ".jl"), tests)
         end
         println()
+        if !isempty(DISABLED_TEST_FILES)
+            @warn("Some tests are disabled $DISABLED_TEST_FILES")
+        end
         for test in tests
+            (test ∈ DISABLED_TEST_FILES) && continue
             print(splitext(test)[1], ": ")
             include(test)
             println()
diff --git a/test/test_component_assignment.jl b/test/test_component_assignment.jl
new file mode 100644
index 0000000..8c3b635
--- /dev/null
+++ b/test/test_component_assignment.jl
@@ -0,0 +1,64 @@
+@testset "Test branches are assigned to subsystems correctly" begin
+    sys = build_system(PSISystems, "two_area_pjm_DA"; add_reserves=true)
+
+    area_subsystem_map = Dict("Area1" => "a", "Area2" => "b")
+    make_subsystems!(sys, area_subsystem_map)
+    # From bus determines subsystem of ACBranch
+    a_buses = get_components(ACBus, sys; subsystem_name="a")
+    b_buses = get_components(ACBus, sys; subsystem_name="b")
+    for b in get_components(ACTransmission, sys)
+        from_bus = get_from(get_arc(b))
+        if from_bus ∈ a_buses
+            add_component_to_subsystem!(sys, "a", b)
+        elseif from_bus ∈ b_buses
+            add_component_to_subsystem!(sys, "b", b)
+        end
+    end
+    template = MultiProblemTemplate(NetworkModel(SplitAreaPTDFPowerModel), ["a", "b"])
+    set_device_model!(template, ThermalStandard, ThermalBasicDispatch)
+    set_device_model!(template, PowerLoad, StaticPowerLoad)
+    set_device_model!(template, DeviceModel(MonitoredLine, StaticBranch))
+    set_device_model!(template, DeviceModel(Line, StaticBranch))
+
+    problem = DecisionModel(
+        MultiRegionProblem,
+        template,
+        sys;
+        horizon=Hour(24),
+        interval=Hour(1),
+        resolution=Hour(1),
+        name="UC_Subsystem",
+        optimizer=HiGHS_optimizer_small_gap,
+    )
+
+    build_out = build!(problem; output_dir=mktempdir())
+
+    @test size(
+        problem.internal.container.subproblems["a"].constraints[PSI.ConstraintKey{
+            FlowRateConstraint,
+            Line,
+        }(
+            "lb",
+        )],
+    )[1] == 6
+    @test size(
+        problem.internal.container.subproblems["b"].constraints[PSI.ConstraintKey{
+            FlowRateConstraint,
+            Line,
+        }(
+            "lb",
+        )],
+    )[1] == 6
+    @test size(
+        problem.internal.container.subproblems["a"].constraints[PSI.ConstraintKey{
+            FlowRateConstraint,
+            MonitoredLine,
+        }(
+            "lb",
+        )],
+    )[1] == 1
+    @test !haskey(
+        problem.internal.container.subproblems["b"].constraints,
+        PSI.ConstraintKey{FlowRateConstraint, MonitoredLine}("lb"),
+    )
+end
diff --git a/test/test_core.jl b/test/test_core.jl
deleted file mode 100644
index 9aec8f9..0000000
--- a/test/test_core.jl
+++ /dev/null
@@ -1,3 +0,0 @@
-@test 1 == 1
-
-#sys = build_system(PSISystems, "sys10_pjm_ac_dc")
diff --git a/test/test_horizontal_passing.jl b/test/test_horizontal_passing.jl
index dcf79f0..230b076 100644
--- a/test/test_horizontal_passing.jl
+++ b/test/test_horizontal_passing.jl
@@ -1,9 +1,8 @@
 
 @testset "Test horizontal passing without emulator" begin
     sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
-    sys2 = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
     results, _ = run_rts_multi_stage_decomposition_simulation(
-        [sys, sys2];
+        sys;
         NT=5,
         mode="horizontal",
         monitored_line_formulations=[StaticBranchUnbounded, StaticBranchUnbounded],
@@ -15,31 +14,55 @@
     #Test "horizontal passing": ActivePowerBalance__ACBus(t-1) =  StateEstimationInjections__ACBus(t).
     #For the real time problem, the state estimation comes from the previous time interval. 
     for b in [string(get_number(x)) for x in get_components(ACBus, sys)]
-        apb = read_realized_variable(results_rt, "ActivePowerBalance__ACBus")[!, b]
-        sei = read_realized_variable(results_rt, "StateEstimationInjections__ACBus")[!, b]
+        apb = read_realized_variable(
+            results_rt,
+            "ActivePowerBalance__ACBus";
+            table_format=TableFormat.WIDE,
+        )[
+            !,
+            b,
+        ]
+        sei = read_realized_variable(
+            results_rt,
+            "StateEstimationInjections__ACBus";
+            table_format=TableFormat.WIDE,
+        )[
+            !,
+            b,
+        ]
         @test isapprox(sei[2:end], apb[1:(end - 1)])
     end
     # Test values to ensure implementation changes aren't causing unexpected changes in results
-    @test read_realized_variable(results_rt, "ActivePowerBalance__ACBus")[1, "116"] ==
-          -0.3456209797192982
-    @test read_realized_variable(results_rt, "ActivePowerBalance__ACBus")[1, "119"] ==
-          -0.6255739732919298
+    @test read_realized_variable(
+        results_rt,
+        "ActivePowerBalance__ACBus";
+        table_format=TableFormat.WIDE,
+    )[
+        1,
+        "116",
+    ] == -0.3456209797192982
+    @test read_realized_variable(
+        results_rt,
+        "ActivePowerBalance__ACBus";
+        table_format=TableFormat.WIDE,
+    )[
+        1,
+        "119",
+    ] == -0.6255739732919298
 end
 
 @testset "Horizontal passing; compare branch models without emulator" begin
     sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
-    sys2 = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
     results_original, _ = run_rts_multi_stage_decomposition_simulation(
-        [sys, sys2];
+        sys;
         NT=5,
         mode="horizontal",
         monitored_line_formulations=[StaticBranchUnbounded, StaticBranchUnbounded],
         use_emulator=false,
     )
     sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
-    sys2 = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
     results_se_line, _ = run_rts_multi_stage_decomposition_simulation(
-        [sys, sys2];
+        sys;
         NT=5,
         mode="horizontal",
         monitored_line_formulations=[
@@ -52,11 +75,13 @@ end
     results_sub_se_line = get_decision_problem_results(results_se_line, "UC_Subsystem")
     flow_sub_original = read_realized_variable(
         results_sub_original,
-        "FlowActivePowerVariable__MonitoredLine",
+        "PTDFBranchFlow__MonitoredLine";
+        table_format=TableFormat.WIDE,
     )
     flow_sub_se_line = read_realized_variable(
         results_sub_se_line,
-        "FlowActivePowerVariable__MonitoredLine",
+        "PTDFBranchFlow__MonitoredLine";
+        table_format=TableFormat.WIDE,
     )
 
     # WITHOUT the emulator, we expect some difference in the flows outside of the first timestep:
@@ -66,10 +91,8 @@ end
 
 @testset "Horizontal passing; compare branch models with emulator" begin
     sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
-    sys2 = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
-    sys3 = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
     results_original, _ = run_rts_multi_stage_decomposition_simulation(
-        [sys, sys2, sys3];
+        sys;
         NT=5,
         mode="horizontal",
         monitored_line_formulations=[
@@ -80,10 +103,8 @@ end
         use_emulator=true,
     )
     sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
-    sys2 = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
-    sys3 = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
     results_se_line, _ = run_rts_multi_stage_decomposition_simulation(
-        [sys, sys2, sys3];
+        sys;
         NT=5,
         mode="horizontal",
         monitored_line_formulations=[
@@ -97,11 +118,13 @@ end
     results_sub_se_line = get_decision_problem_results(results_se_line, "UC_Subsystem")
     flow_sub_original = read_realized_variable(
         results_sub_original,
-        "FlowActivePowerVariable__MonitoredLine",
+        "PTDFBranchFlow__MonitoredLine";
+        table_format=TableFormat.WIDE,
     )
     flow_sub_se_line = read_realized_variable(
         results_sub_se_line,
-        "FlowActivePowerVariable__MonitoredLine",
+        "PTDFBranchFlow__MonitoredLine";
+        table_format=TableFormat.WIDE,
     )
 
     # WITH the emulator, we expect the formulations to be equivalent:
diff --git a/test/test_hvdc.jl b/test/test_hvdc.jl
index 0ca2d9b..c2255f5 100644
--- a/test/test_hvdc.jl
+++ b/test/test_hvdc.jl
@@ -2,7 +2,7 @@
     sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
     area_subsystem_map = Dict("1" => "a", "2" => "b", "3" => "b")
     make_subsystems!(sys, area_subsystem_map)
-    hvdc = get_component(TwoTerminalHVDCLine, sys, "DC1")
+    hvdc = get_component(TwoTerminalGenericHVDCLine, sys, "DC1")
     PowerSystems.add_component_to_subsystem!(sys, "a", hvdc)
     template_uc2 = MultiProblemTemplate(NetworkModel(SplitAreaPTDFPowerModel), ["a", "b"])
     problem = DecisionModel(
@@ -15,34 +15,35 @@
     @test build!(problem; output_dir=mktempdir()) == PowerSimulations.ModelBuildStatus.BUILT
 end
 
-sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
-area_subsystem_map = Dict("1" => "a", "2" => "b", "3" => "b")
-make_subsystems!(sys, area_subsystem_map)
-hvdc = get_component(TwoTerminalHVDCLine, sys, "DC1")
-PowerSystems.add_component_to_subsystem!(sys, "a", hvdc)
-template_uc2 = MultiProblemTemplate(NetworkModel(SplitAreaPTDFPowerModel), ["a", "b"])
-set_device_model!(template_uc2, TwoTerminalHVDCLine, HVDCTwoTerminalLossless)
-problem = DecisionModel(
-    MultiRegionProblem,
-    template_uc2,
-    sys;
-    name="UC_Subsystem",
-    optimizer=optimizer_with_attributes(HiGHS.Optimizer),
-)
-build_out = build!(problem; console_level=Logging.AboveMaxLevel, output_dir=mktempdir())
-@test build_out == PowerSimulations.ModelBuildStatus.FAILED
-
+@testset "HVDC spans subsystems and is modeled (build fails)" begin
+    sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
+    area_subsystem_map = Dict("1" => "a", "2" => "b", "3" => "b")
+    make_subsystems!(sys, area_subsystem_map)
+    hvdc = get_component(TwoTerminalGenericHVDCLine, sys, "DC1")
+    PowerSystems.add_component_to_subsystem!(sys, "a", hvdc)
+    template_uc2 = MultiProblemTemplate(NetworkModel(SplitAreaPTDFPowerModel), ["a", "b"])
+    set_device_model!(template_uc2, TwoTerminalGenericHVDCLine, HVDCTwoTerminalLossless)
+    problem = DecisionModel(
+        MultiRegionProblem,
+        template_uc2,
+        sys;
+        name="UC_Subsystem",
+        optimizer=optimizer_with_attributes(HiGHS.Optimizer),
+    )
+    build_out = build!(problem; console_level=Logging.AboveMaxLevel, output_dir=mktempdir())
+    @test build_out == PowerSimulations.ModelBuildStatus.FAILED
+end
 
 @testset "HVDC spans subsystems and is modeled but both terminal buses belong to same subsytem (build suceeds)" begin
     sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
     area_subsystem_map = Dict("1" => "a", "2" => "b", "3" => "b")
     make_subsystems!(sys, area_subsystem_map)
-    hvdc = get_component(TwoTerminalHVDCLine, sys, "DC1")
+    hvdc = get_component(TwoTerminalGenericHVDCLine, sys, "DC1")
     PowerSystems.add_component_to_subsystem!(sys, "a", hvdc)
     PowerSystems.remove_component_from_subsystem!(sys, "b", hvdc.arc.to)
     PowerSystems.add_component_to_subsystem!(sys, "a", hvdc.arc.to)
     template_uc2 = MultiProblemTemplate(NetworkModel(SplitAreaPTDFPowerModel), ["a", "b"])
-    set_device_model!(template_uc2, TwoTerminalHVDCLine, HVDCTwoTerminalLossless)
+    set_device_model!(template_uc2, TwoTerminalGenericHVDCLine, HVDCTwoTerminalLossless)
     problem = DecisionModel(
         MultiRegionProblem,
         template_uc2,
@@ -58,12 +59,12 @@ end
     sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
     area_subsystem_map = Dict("1" => "a", "2" => "b", "3" => "b")
     make_subsystems!(sys, area_subsystem_map)
-    hvdc = get_component(TwoTerminalHVDCLine, sys, "DC1")
+    hvdc = get_component(TwoTerminalGenericHVDCLine, sys, "DC1")
     PowerSystems.add_component_to_subsystem!(sys, "a", hvdc)
     PowerSystems.remove_component_from_subsystem!(sys, "b", hvdc.arc.to)
     PowerSystems.add_component_to_subsystem!(sys, "a", hvdc.arc.to)
     template_uc2 = MultiProblemTemplate(NetworkModel(SplitAreaPTDFPowerModel), ["a", "b"])
-    set_device_model!(template_uc2, TwoTerminalHVDCLine, HVDCTwoTerminalLossless)
+    set_device_model!(template_uc2, TwoTerminalGenericHVDCLine, HVDCTwoTerminalLossless)
     problem = DecisionModel(
         MultiRegionProblem,
         template_uc2,
@@ -79,12 +80,12 @@ end
     sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
     area_subsystem_map = Dict("1" => "a", "2" => "b", "3" => "b")
     make_subsystems!(sys, area_subsystem_map)
-    hvdc = get_component(TwoTerminalHVDCLine, sys, "DC1")
+    hvdc = get_component(TwoTerminalGenericHVDCLine, sys, "DC1")
     PowerSystems.add_component_to_subsystem!(sys, "a", hvdc)
     PowerSystems.remove_component_from_subsystem!(sys, "b", hvdc.arc.to)
     PowerSystems.add_component_to_subsystem!(sys, "a", hvdc.arc.to)
     template_uc2 = MultiProblemTemplate(NetworkModel(SplitAreaPTDFPowerModel), ["a", "b"])
-    set_device_model!(template_uc2, TwoTerminalHVDCLine, HVDCTwoTerminalLossless)
+    set_device_model!(template_uc2, TwoTerminalGenericHVDCLine, HVDCTwoTerminalLossless)
     problem = DecisionModel(
         MultiRegionProblem,
         template_uc2,
@@ -114,6 +115,6 @@ end
         name="UC_Subsystem",
         optimizer=optimizer_with_attributes(HiGHS.Optimizer),
     )
-    build_out = build!(problem; output_dir=mktempdir())
+    build_out = build!(problem; console_level=Logging.AboveMaxLevel, output_dir=mktempdir())
     @test build_out == PowerSimulations.ModelBuildStatus.FAILED
 end
diff --git a/test/test_hydro.jl b/test/test_hydro.jl
new file mode 100644
index 0000000..fd5644c
--- /dev/null
+++ b/test/test_hydro.jl
@@ -0,0 +1,14 @@
+@testset "Test decomposition with hydro budget (1D results)" begin
+    sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
+    add_hydro_budget_time_series_to_rts!(sys)
+    results, _ = run_rts_multi_stage_decomposition_simulation(
+        sys;
+        NT=5,
+        mode="horizontal",
+        monitored_line_formulations=[StaticBranchUnbounded, StaticBranchUnbounded],
+        use_emulator=false,
+        add_reserves=false,
+        add_hydro=true,
+    )
+    @test isa(results, SimulationResults)
+end
diff --git a/test/test_reserve_deliverability.jl b/test/test_reserve_deliverability.jl
new file mode 100644
index 0000000..5d2a63d
--- /dev/null
+++ b/test/test_reserve_deliverability.jl
@@ -0,0 +1,347 @@
+#Utility function for testing for adding outage data to generators and reserves:
+function add_outages_to_systems!(
+    systems::Vector{System},
+    outage_specifications::Vector{
+        @NamedTuple{
+            outage_generators::Vector{String},
+            responding_reserves::Dict{DataType, String},
+        }
+    },
+)
+    for outage_specification in outage_specifications
+        outage_generators = outage_specification.outage_generators
+        responding_reserves = outage_specification.responding_reserves
+        transition_data = GeometricDistributionForcedOutage(;
+            mean_time_to_recovery=10,  # Units of hours - This value does not have any influence for G-1 formulation
+            outage_transition_probability=1.0,
+        )
+        for sys in systems
+            for gen_name in outage_generators
+                gen = get_component(PSY.Generator, sys, gen_name)
+                add_supplemental_attribute!(sys, gen, transition_data)
+            end
+            for (reserve, reserve_name) in responding_reserves
+                reserve_product = get_component(reserve, sys, reserve_name)
+                add_supplemental_attribute!(sys, reserve_product, transition_data)
+            end
+        end
+    end
+end
+
+function _get_reserve_name_for_results(reserve_type::DataType, reserve_name::String)
+    a = string(Base.typename(reserve_type).wrapper)
+    b = reserve_type.parameters[1]
+    c = reserve_name
+    return "__" * string(a) * "__" * string(b) * "__" * string(c)
+end
+
+@testset "10 bus; RampReserveWithDeliverabilityConstraints + AreaPTDFPowerModel (reference): separate reserves" begin
+    sys = build_system(PSISystems, "two_area_pjm_DA"; add_reserves=true)
+    # Make feasible by ignoring ramp limits: 
+    for g in get_components(ThermalStandard, sys)
+        set_ramp_limits!(g, nothing)
+    end
+
+    # Make Sundance must run with a minimum active power: 
+    set_must_run!(get_component(ThermalStandard, sys, "Sundance_1"), true)
+    set_active_power_limits!(
+        get_component(ThermalStandard, sys, "Sundance_1"),
+        (min=1.0, max=2.0),
+    )
+    set_must_run!(get_component(ThermalStandard, sys, "Sundance_2"), true)
+    set_active_power_limits!(
+        get_component(ThermalStandard, sys, "Sundance_2"),
+        (min=1.0, max=2.0),
+    )
+
+    outages_specifications = [
+        (
+            outage_generators=["Sundance_1"],
+            responding_reserves=Dict(PSY.VariableReserve{ReserveUp} => "Reserve1_1"),
+        )
+        (
+            outage_generators=["Sundance_2"],
+            responding_reserves=Dict(PSY.VariableReserve{ReserveUp} => "Reserve1_2"),
+        )
+    ]
+    add_outages_to_systems!([sys], outages_specifications)
+
+    template = ProblemTemplate(NetworkModel(AreaPTDFPowerModel))
+    set_device_model!(template, ThermalStandard, ThermalBasicDispatch)
+    set_device_model!(template, PowerLoad, StaticPowerLoad)
+    set_device_model!(template, DeviceModel(MonitoredLine, StaticBranch))  # Only assign Monitored Line to subsystem A template
+    set_device_model!(template, Line, StaticBranch)
+    set_service_model!(
+        template,
+        ServiceModel(
+            VariableReserve{ReserveUp},
+            RampReserveWithDeliverabilityConstraints,
+            "Reserve1_1",
+        ),
+    )
+    set_service_model!(
+        template,
+        ServiceModel(
+            VariableReserve{ReserveUp},
+            RampReserveWithDeliverabilityConstraints,
+            "Reserve1_2",
+        ),
+    )
+    problem = DecisionModel(
+        template,
+        sys;
+        horizon=Hour(24),
+        interval=Hour(1),
+        resolution=Hour(1),
+        name="UC_Subsystem",
+        optimizer=HiGHS_optimizer_small_gap,
+    )
+
+    build_out = build!(problem; output_dir=mktempdir())
+    solve!(problem)
+    res = OptimizationProblemResults(problem)
+
+    post_contingency_deployment_1 = read_variable(
+        res,
+        "PostContingencyActivePowerReserveDeploymentVariable__VariableReserve__ReserveUp__Reserve1_1",
+    )
+    post_contingency_deployment_2 = read_variable(
+        res,
+        "PostContingencyActivePowerReserveDeploymentVariable__VariableReserve__ReserveUp__Reserve1_2",
+    )
+    # Sum of the reserve deployments is equivalent to the outaged generator: 
+    for post_contingency_deployment in
+        [post_contingency_deployment_1, post_contingency_deployment_2]
+        for i in unique(post_contingency_deployment[!, :DateTime])
+            @test isapprox(
+                sum(
+                    filter(row -> row.DateTime == i, post_contingency_deployment)[
+                        :,
+                        "value",
+                    ],
+                ),
+                100.0,
+            )
+        end
+    end
+
+    flows = read_expression(res, "PTDFBranchFlow__Line", table_format=TableFormat.WIDE)
+    post_contingency_flows_1 = read_expression(
+        res,
+        "PostContingencyBranchFlow__VariableReserve__ReserveUp__Reserve1_1",
+    )
+    post_contingency_flows_2 = read_expression(
+        res,
+        "PostContingencyBranchFlow__VariableReserve__ReserveUp__Reserve1_2",
+    )
+    for line_name in get_name.(get_components(Line, sys))
+        flow = flows[:, line_name]
+        post_contingency_flow_1 =
+            filter(row -> row.name2 == line_name, post_contingency_flows_1)[:, :value]
+        post_contingency_flow_2 =
+            filter(row -> row.name2 == line_name, post_contingency_flows_2)[:, :value]
+        # Post contingency flow in the area of the outage is different, post contingency flow from the other area is the same:
+        if occursin("_1", line_name)
+            @test all(isapprox.(flow, post_contingency_flow_2))
+            @test !all(isapprox.(flow, post_contingency_flow_1))
+        elseif occursin("_2", line_name)
+            @test all(isapprox.(flow, post_contingency_flow_1))
+            @test !all(isapprox.(flow, post_contingency_flow_2))
+        end
+    end
+end
+
+@testset "10 bus; RampReserveWithDeliverabilityConstraints + SplitAreaPTDFPowerModel: separate reserves" begin
+    sys = build_system(PSISystems, "two_area_pjm_DA"; add_reserves=true)
+    # Make feasible by ignoring ramp limits: 
+    for g in get_components(ThermalStandard, sys)
+        set_ramp_limits!(g, nothing)
+    end
+
+    # Make Sundance must run with a minimum active power: 
+    set_must_run!(get_component(ThermalStandard, sys, "Sundance_1"), true)
+    set_active_power_limits!(
+        get_component(ThermalStandard, sys, "Sundance_1"),
+        (min=1.0, max=2.0),
+    )
+    set_must_run!(get_component(ThermalStandard, sys, "Sundance_2"), true)
+    set_active_power_limits!(
+        get_component(ThermalStandard, sys, "Sundance_2"),
+        (min=1.0, max=2.0),
+    )
+
+    outages_specifications = [
+        (
+            outage_generators=["Sundance_1"],
+            responding_reserves=Dict(PSY.VariableReserve{ReserveUp} => "Reserve1_1"),
+        )
+        (
+            outage_generators=["Sundance_2"],
+            responding_reserves=Dict(PSY.VariableReserve{ReserveUp} => "Reserve1_2"),
+        )
+    ]
+    add_outages_to_systems!([sys], outages_specifications)
+
+    area_subsystem_map = Dict("Area1" => "a", "Area2" => "b")
+    make_subsystems!(sys, area_subsystem_map)
+    r1_1 = get_component(VariableReserve, sys, "Reserve1_1")
+    add_component_to_subsystem!(sys, "a", r1_1)
+    r1_2 = get_component(VariableReserve, sys, "Reserve1_2")
+    add_component_to_subsystem!(sys, "b", r1_2)
+
+    # From bus determines subsystem of ACBranch
+    a_buses = get_components(ACBus, sys; subsystem_name="a")
+    b_buses = get_components(ACBus, sys; subsystem_name="b")
+    for b in get_components(ACTransmission, sys)
+        from_bus = get_from(get_arc(b))
+        if from_bus ∈ a_buses
+            add_component_to_subsystem!(sys, "a", b)
+        elseif from_bus ∈ b_buses
+            add_component_to_subsystem!(sys, "b", b)
+        end
+    end
+    template = MultiProblemTemplate(NetworkModel(SplitAreaPTDFPowerModel), ["a", "b"])
+    set_device_model!(template, ThermalStandard, ThermalBasicDispatch)
+    set_device_model!(template, PowerLoad, StaticPowerLoad)
+    set_device_model!(template, DeviceModel(MonitoredLine, StaticBranch))
+    set_device_model!(template, Line, StaticBranch)
+    set_service_model!(
+        template,
+        ServiceModel(
+            VariableReserve{ReserveUp},
+            RampReserveWithDeliverabilityConstraints,
+            "Reserve1_1",
+        ),
+        "a",
+    )
+    set_service_model!(
+        template,
+        ServiceModel(
+            VariableReserve{ReserveUp},
+            RampReserveWithDeliverabilityConstraints,
+            "Reserve1_2",
+        ),
+        "b",
+    )
+    problem = DecisionModel(
+        MultiRegionProblem,
+        template,
+        sys;
+        horizon=Hour(24),
+        interval=Hour(1),
+        resolution=Hour(1),
+        name="UC_Subsystem",
+        optimizer=HiGHS_optimizer_small_gap,
+    )
+
+    build_out = build!(problem; output_dir=mktempdir())
+    solve!(problem)
+    res = OptimizationProblemResults(problem)
+    post_contingency_deployment_1 = get_variable_values(res)[PSI.VariableKey{
+        PostContingencyActivePowerReserveDeploymentVariable,
+        VariableReserve{ReserveUp},
+    }(
+        "Reserve1_1",
+    )]
+    post_contingency_deployment_2 = get_variable_values(res)[PSI.VariableKey{
+        PostContingencyActivePowerReserveDeploymentVariable,
+        VariableReserve{ReserveUp},
+    }(
+        "Reserve1_2",
+    )]
+    # Sum of the reserve deployments is equivalent to the outaged generator: 
+    for post_contingency_deployment in
+        [post_contingency_deployment_1, post_contingency_deployment_2]
+        for i in unique(post_contingency_deployment[!, :time_index])
+            @test isapprox(
+                sum(
+                    filter(row -> row.time_index == i, post_contingency_deployment)[
+                        :,
+                        "value",
+                    ],
+                ),
+                1.0,
+            )
+        end
+    end
+
+    flows = get_expression_values(res)[PSI.ExpressionKey{PTDFBranchFlow, Line}("")]
+    post_contingency_flows_1 = get_expression_values(res)[PSI.ExpressionKey{
+        PostContingencyBranchFlow,
+        VariableReserve{ReserveUp},
+    }(
+        "Reserve1_1",
+    )]
+    post_contingency_flows_2 = get_expression_values(res)[PSI.ExpressionKey{
+        PostContingencyBranchFlow,
+        VariableReserve{ReserveUp},
+    }(
+        "Reserve1_2",
+    )]
+    for line_name in get_name.(get_components(Line, sys))
+        flow = filter(row -> row.name == line_name, flows)[:, :value]
+        post_contingency_flow_1 =
+            filter(row -> row.name2 == line_name, post_contingency_flows_1)[:, :value]
+        post_contingency_flow_2 =
+            filter(row -> row.name2 == line_name, post_contingency_flows_2)[:, :value]
+        # Post contingency flow in the area of the outage is different
+        if occursin("_1", line_name)
+            @test !all(isapprox.(flow, post_contingency_flow_1; atol=1e-6))
+        elseif occursin("_2", line_name)
+            @test !all(isapprox.(flow, post_contingency_flow_2; atol=1e-6))
+        end
+    end
+end
+
+@testset "RTS multi-stage sim w/ RampReserveWithDeliverabilityConstraints and SplitAreaPTDFPowerModel" begin
+    sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
+    outages_specifications = [(
+        outage_generators=["313_CC_1"],
+        responding_reserves=Dict(PSY.VariableReserve{ReserveUp} => "Spin_Up_R3"),
+    )]
+    add_outages_to_systems!([sys], outages_specifications)
+    results, sim = run_rts_multi_stage_decomposition_simulation(
+        sys;
+        NT=5,
+        mode="vertical",
+        monitored_line_formulations=[StaticBranchUnbounded, StaticBranchUnbounded],
+        use_emulator=false,
+        add_reserves=true,
+        in_memory=true, # fails with in_memory = false
+    )
+    results_uc0 = get_decision_problem_results(results, "UC0")
+    results_uc_subsystem = get_decision_problem_results(results, "UC_Subsystem")
+    active_power = read_realized_variable(
+        results_uc_subsystem,
+        "ActivePowerVariable__ThermalStandard";
+        table_format=TableFormat.WIDE,
+    )
+    # Test sum of the responding reserve deployments is equivalent to the total outaged generation: 
+    for outages_specification in outages_specifications
+        outage_generators = outages_specification.outage_generators
+        responding_reserves = outages_specification.responding_reserves
+        for i in unique(active_power[!, :DateTime])
+            total_reserve_deployment = 0.0
+            total_lost_generation = 0.0
+            for (reserve_type, reserve_name) in responding_reserves
+                # NOTE: workaround because read_realized_variable broken for 3D results: https://github.com/NREL-Sienna/PowerSimulations.jl/issues/1390
+                post_contingency_deployment = read_variable(
+                    results_uc_subsystem,
+                    "PostContingencyActivePowerReserveDeploymentVariable" *
+                    _get_reserve_name_for_results(reserve_type, reserve_name),
+                )[DateTime("2020-01-01T00:00:00")]
+                total_reserve_deployment += sum(
+                    filter(row -> row.DateTime == i, post_contingency_deployment)[
+                        :,
+                        "value",
+                    ],
+                )
+            end
+            for generator in outage_generators
+                total_lost_generation +=
+                    filter(row -> row.DateTime == i, active_power)[1, generator]
+            end
+            @test isapprox(total_lost_generation, total_reserve_deployment)
+        end
+    end
+end
diff --git a/test/test_reserves.jl b/test/test_reserves.jl
index d62269c..8ac5bf8 100644
--- a/test/test_reserves.jl
+++ b/test/test_reserves.jl
@@ -23,7 +23,6 @@ end
     make_subsystems!(sys, area_subsystem_map)
     template = MultiProblemTemplate(NetworkModel(SplitAreaPTDFPowerModel), ["a", "b"])
     set_device_model!(template, ThermalStandard, ThermalBasicUnitCommitment)
-    set_service_model!(template, VariableReserve{ReserveUp}, RangeReserve)
     problem = DecisionModel(
         MultiRegionProblem,
         template,
@@ -34,8 +33,8 @@ end
     build_out = build!(problem; output_dir=mktempdir())
     @test build_out == PowerSimulations.ModelBuildStatus.BUILT
     jump_problem_dict = get_jump_models(problem)
-    moi_tests(jump_problem_dict["a"], 10434, 0, 1728, 864, 2640, true)
-    moi_tests(jump_problem_dict["b"], 17364, 0, 3456, 1728, 5280, true)
+    moi_tests(jump_problem_dict["a"], 10470, 0, 1728, 864, 2640, true)
+    moi_tests(jump_problem_dict["b"], 17436, 0, 3456, 1728, 5280, true)
 end
 
 @testset "MOI test - reserves in A" begin
@@ -49,6 +48,8 @@ end
             end
         end
     end
+    remove_component!(sys, get_component(VariableReserve, sys, "Reg_Up"))
+    remove_component!(sys, get_component(VariableReserve, sys, "Reg_Down"))
     area_subsystem_map = Dict("1" => "a", "2" => "b", "3" => "b")
     make_subsystems!(sys, area_subsystem_map)
     r1 = get_component(VariableReserve, sys, "Spin_Up_R1")
@@ -66,8 +67,8 @@ end
     build_out = build!(problem; output_dir=mktempdir())
     @test build_out == PowerSimulations.ModelBuildStatus.BUILT
     jump_problem_dict = get_jump_models(problem)
-    moi_tests(jump_problem_dict["a"], 11250, 0, 1728, 912, 2640, true)
-    moi_tests(jump_problem_dict["b"], 17364, 0, 3456, 1728, 5280, true)
+    moi_tests(jump_problem_dict["a"], 11286, 0, 1728, 912, 2640, true)
+    moi_tests(jump_problem_dict["b"], 17436, 0, 3456, 1728, 5280, true)
 end
 
 @testset "MOI test - reserves in B" begin
@@ -81,6 +82,8 @@ end
             end
         end
     end
+    remove_component!(sys, get_component(VariableReserve, sys, "Reg_Up"))
+    remove_component!(sys, get_component(VariableReserve, sys, "Reg_Down"))
     area_subsystem_map = Dict("1" => "a", "2" => "b", "3" => "b")
     make_subsystems!(sys, area_subsystem_map)
     r2 = get_component(VariableReserve, sys, "Spin_Up_R2")
@@ -98,6 +101,6 @@ end
     build_out = build!(problem; output_dir=mktempdir())
     @test build_out == PowerSimulations.ModelBuildStatus.BUILT
     jump_problem_dict = get_jump_models(problem)
-    moi_tests(jump_problem_dict["a"], 10434, 0, 1728, 864, 2640, true)
-    moi_tests(jump_problem_dict["b"], 18276, 0, 3456, 1776, 5280, true)
+    moi_tests(jump_problem_dict["a"], 10470, 0, 1728, 864, 2640, true)
+    moi_tests(jump_problem_dict["b"], 18348, 0, 3456, 1776, 5280, true)
 end
diff --git a/test/test_utils/model_checks.jl b/test/test_utils/model_checks.jl
index f8cef48..3a5da8f 100644
--- a/test/test_utils/model_checks.jl
+++ b/test/test_utils/model_checks.jl
@@ -1,9 +1,17 @@
 const GAEVF = JuMP.GenericAffExpr{Float64, VariableRef}
 const GQEVF = JuMP.GenericQuadExpr{Float64, VariableRef}
 
+HiGHS_optimizer_small_gap = JuMP.optimizer_with_attributes(
+    HiGHS.Optimizer,
+    "time_limit" => 100.0,
+    "random_seed" => 12345,
+    "mip_rel_gap" => 0.001,
+    "log_to_console" => false,
+)
+
 function get_jump_models(model::PSI.DecisionModel{MultiRegionProblem})
     jump_model_dict = Dict()
-    subproblem_keys = keys(problem.internal.container.subproblems)
+    subproblem_keys = keys(model.internal.container.subproblems)
     for k in subproblem_keys
         jump_model_dict[k] = PSI.get_jump_model(
             IS.Optimization.get_container(PSI.get_internal(model)).subproblems[k],
diff --git a/test/test_utils/simulation_runs.jl b/test/test_utils/simulation_runs.jl
index 1059b48..67eb4f7 100644
--- a/test/test_utils/simulation_runs.jl
+++ b/test/test_utils/simulation_runs.jl
@@ -1,42 +1,39 @@
 function run_rts_multi_stage_decomposition_simulation(
-    systems;
+    sys;
     NT=5,
     mode="vertical",
     monitored_line_formulations=[StaticBranchUnbounded, StaticBranchUnbounded],
     use_emulator=false,
+    add_reserves=false,
+    add_hydro=false,
+    in_memory=false,
 )
     modeled_lines = ["CA-1", "CB-1", "AB1", "A28"]
     convert_to_monitored_line = [(name="A28", flow_limit=20.0)]
-    sys = systems[1]
-    sys2 = systems[2]
-    transform_single_time_series!(sys, Hour(NT), Hour(NT))
-    if mode == "vertical"
-        transform_single_time_series!(sys2, Hour(NT), Hour(NT))
-    elseif mode == "horizontal"
-        transform_single_time_series!(sys2, Hour(1), Hour(1))
-    end
-    if use_emulator
-        @assert length(systems) == 3
-        sys3 = systems[3]
-        transform_single_time_series!(sys3, Hour(1), Hour(1))
-    end
     area_subsystem_map = Dict("1" => "a", "2" => "b", "3" => "a")
 
-    make_subsystems!(sys2, area_subsystem_map)
+    make_subsystems!(sys, area_subsystem_map)
 
     for b in modeled_lines
-        l = get_component(ACBranch, sys2, b)
-        PowerSystems.add_component_to_subsystem!(sys2, "a", l)
-        #PowerSystems.add_component_to_subsystem!(sys2, "b", l)
+        l = get_component(ACBranch, sys, b)
+        to_bus_area = get_name(get_area(get_to(get_arc(l))))
+        PowerSystems.add_component_to_subsystem!(sys, area_subsystem_map[to_bus_area], l)
+        @info "Assigning line $(get_name(l)) to subsystem $(area_subsystem_map[to_bus_area])"
     end
-    for sys in systems
-        add_interchanges!(sys)
-        add_monitored_lines!(sys, convert_to_monitored_line)
+    #Add all area interchanges and include in all subsystems
+    add_interchanges!(sys)
+    for b in get_components(AreaInterchange, sys)
+        for subsystem in unique(values(area_subsystem_map))
+            add_component_to_subsystem!(sys, subsystem, b)
+        end
     end
+    add_monitored_lines!(sys, convert_to_monitored_line)
 
     template_uc = ProblemTemplate(NetworkModel(AreaPTDFPowerModel; use_slacks=true))
     set_device_model!(template_uc, ThermalStandard, ThermalBasicUnitCommitment)
     set_device_model!(template_uc, PowerLoad, StaticPowerLoad)
+    add_hydro &&
+        set_device_model!(template_uc, HydroDispatch, HydroDispatchRunOfRiverBudget)
     set_device_model!(template_uc, AreaInterchange, StaticBranch)
     set_device_model!(
         template_uc,
@@ -56,6 +53,33 @@ function run_rts_multi_stage_decomposition_simulation(
             attributes=Dict("filter_function" => x -> get_name(x) in modeled_lines),
         ),
     )
+    if add_reserves
+        # add reserve formulations:
+        set_service_model!(
+            template_uc,
+            ServiceModel(
+                VariableReserve{ReserveUp},
+                RampReserveWithDeliverabilityConstraints,
+                "Spin_Up_R1",
+            ),
+        )
+        set_service_model!(
+            template_uc,
+            ServiceModel(
+                VariableReserve{ReserveUp},
+                RampReserveWithDeliverabilityConstraints,
+                "Spin_Up_R2",
+            ),
+        )
+        set_service_model!(
+            template_uc,
+            ServiceModel(
+                VariableReserve{ReserveUp},
+                RampReserveWithDeliverabilityConstraints,
+                "Spin_Up_R3",
+            ),
+        )
+    end
     # Set up Model 2 (MultiProblem)
     template_uc2 = MultiProblemTemplate(
         NetworkModel(SplitAreaPTDFPowerModel; use_slacks=true),
@@ -63,6 +87,8 @@ function run_rts_multi_stage_decomposition_simulation(
     )
     set_device_model!(template_uc2, ThermalStandard, ThermalBasicUnitCommitment)
     set_device_model!(template_uc2, PowerLoad, StaticPowerLoad)
+    add_hydro &&
+        set_device_model!(template_uc2, HydroDispatch, HydroDispatchRunOfRiverBudget)
     set_device_model!(template_uc2, AreaInterchange, StaticBranch)
     set_device_model!(
         template_uc2,
@@ -82,10 +108,50 @@ function run_rts_multi_stage_decomposition_simulation(
             attributes=Dict("filter_function" => x -> get_name(x) in modeled_lines),
         ),
     )
+    if add_reserves
+        # add modeled reserve components to subsystems: 
+        r1 = get_component(VariableReserve{ReserveUp}, sys, "Spin_Up_R1")
+        add_component_to_subsystem!(sys, "a", r1)
+        r2 = get_component(VariableReserve{ReserveUp}, sys, "Spin_Up_R2")
+        add_component_to_subsystem!(sys, "b", r2)
+        r3 = get_component(VariableReserve{ReserveUp}, sys, "Spin_Up_R3")
+        add_component_to_subsystem!(sys, "a", r3)
+
+        # add reserve formulations (also to specific subsystems):
+        set_service_model!(
+            template_uc2,
+            ServiceModel(
+                VariableReserve{ReserveUp},
+                RampReserveWithDeliverabilityConstraints,
+                "Spin_Up_R1",
+            ),
+            "a",
+        )
+        set_service_model!(
+            template_uc2,
+            ServiceModel(
+                VariableReserve{ReserveUp},
+                RampReserveWithDeliverabilityConstraints,
+                "Spin_Up_R2",
+            ),
+            "b",
+        )
+        set_service_model!(
+            template_uc2,
+            ServiceModel(
+                VariableReserve{ReserveUp},
+                RampReserveWithDeliverabilityConstraints,
+                "Spin_Up_R3",
+            ),
+            "a",
+        )
+    end
     if use_emulator
         template_em = ProblemTemplate(NetworkModel(AreaPTDFPowerModel; use_slacks=true))
         set_device_model!(template_em, ThermalStandard, ThermalBasicUnitCommitment)
         set_device_model!(template_em, PowerLoad, StaticPowerLoad)
+        add_hydro &&
+            set_device_model!(template_em, HydroDispatch, HydroDispatchRunOfRiverBudget)
         set_device_model!(template_em, AreaInterchange, StaticBranch)
         set_device_model!(
             template_em,
@@ -106,6 +172,13 @@ function run_rts_multi_stage_decomposition_simulation(
             ),
         )
     end
+    if mode == "vertical"
+        d2_horizon = Hour(NT)
+        d2_interval = Hour(NT)
+    elseif mode == "horizontal"
+        d2_horizon = Hour(1)
+        d2_interval = Hour(1)
+    end
     models = SimulationModels(;
         decision_models=[
             DecisionModel(
@@ -113,7 +186,9 @@ function run_rts_multi_stage_decomposition_simulation(
                 sys;
                 name="UC0",
                 optimizer=optimizer_with_attributes(HiGHS.Optimizer),
-                system_to_file=false,
+                horizon=Hour(NT),
+                interval=Hour(NT),
+                resolution=Hour(1),
                 optimizer_solve_log_print=false,
                 direct_mode_optimizer=true,
                 store_variable_names=true,
@@ -122,10 +197,12 @@ function run_rts_multi_stage_decomposition_simulation(
             DecisionModel(
                 MultiRegionProblem,
                 template_uc2,
-                sys2;
+                sys;
                 name="UC_Subsystem",
                 optimizer=optimizer_with_attributes(HiGHS.Optimizer),
-                system_to_file=false,
+                horizon=d2_horizon,
+                interval=d2_interval,
+                resolution=Hour(1),
                 initialize_model=true,
                 optimizer_solve_log_print=false,
                 direct_mode_optimizer=true,
@@ -137,9 +214,10 @@ function run_rts_multi_stage_decomposition_simulation(
         emulation_model=use_emulator ?
                         EmulationModel(
             template_em,
-            sys3;
+            sys;
             name="EM",
             optimizer=optimizer_with_attributes(HiGHS.Optimizer),
+            resolution=Hour(1),
             store_variable_names=true,
         ) : nothing,
     )
@@ -168,8 +246,8 @@ function run_rts_multi_stage_decomposition_simulation(
         simulation_folder=mktempdir(),
     )
 
-    build_out = build!(sim; console_level=Logging.Info, serialize=false)
-    execute_status = execute!(sim; enable_progress_bar=true)
+    build_out = build!(sim; console_level=Logging.Info)
+    execute_status = execute!(sim; in_memory=in_memory, enable_progress_bar=true)
 
     return SimulationResults(sim), sim
 end
diff --git a/test/test_utils/system_modifications.jl b/test/test_utils/system_modifications.jl
index b2a69d9..1b8d27a 100644
--- a/test/test_utils/system_modifications.jl
+++ b/test/test_utils/system_modifications.jl
@@ -1,3 +1,11 @@
+function add_hydro_budget_time_series_to_rts!(sys)
+    tstamp = range(DateTime("2020-01-01T00:00:00"); step=Dates.Hour(1), length=8784)
+    data = ones(length(tstamp))
+    ts = SingleTimeSeries("hydro_budget", TimeArray(tstamp, data))
+    add_time_series!(sys, first(get_components(HydroDispatch, sys)), ts)
+    transform_single_time_series!(sys, Day(2), Day(1))
+    return
+end
 
 function make_subsystems!(sys, area_subsystem_map)
     subsystems = unique([v for (_, v) in area_subsystem_map])
diff --git a/test/test_vertical_passing.jl b/test/test_vertical_passing.jl
index b61637b..69fcc7e 100644
--- a/test/test_vertical_passing.jl
+++ b/test/test_vertical_passing.jl
@@ -1,8 +1,7 @@
 @testset "Test vertical passing without emulator" begin
     sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
-    sys2 = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
     results, sim = run_rts_multi_stage_decomposition_simulation(
-        [sys, sys2];
+        sys;
         NT=5,
         mode="vertical",
         monitored_line_formulations=[StaticBranchUnbounded, StaticBranchUnbounded],
@@ -11,57 +10,80 @@
     results_uc0 = get_decision_problem_results(results, "UC0")
     results_ucsub = get_decision_problem_results(results, "UC_Subsystem")
     results_em = get_emulation_problem_results(results)
-    read_realized_variable(results_uc0, "ActivePowerBalance__ACBus")
+    read_realized_variable(
+        results_uc0,
+        "ActivePowerBalance__ACBus";
+        table_format=TableFormat.WIDE,
+    )
     # Tests "vertical passing": ActivePowerBalance__ACBus from the full system problem are passed as StateEstimationInjections__ACBus
     # for the MultiProblem (for the same timesteps)
     for b in [string(get_number(x)) for x in get_components(ACBus, sys)]
-        apb = read_realized_variable(results_uc0, "ActivePowerBalance__ACBus")[!, b]
-        sei =
-            read_realized_variable(results_ucsub, "StateEstimationInjections__ACBus")[!, b]
+        apb = read_realized_variable(
+            results_uc0,
+            "ActivePowerBalance__ACBus";
+            table_format=TableFormat.WIDE,
+        )[
+            !,
+            b,
+        ]
+        sei = read_realized_variable(
+            results_ucsub,
+            "StateEstimationInjections__ACBus";
+            table_format=TableFormat.WIDE,
+        )[
+            !,
+            b,
+        ]
         @test isapprox(apb, sei)
     end
     # Test values to ensure implementation changes aren't causing unexpected changes in results
-    @test read_realized_variable(results_uc0, "ActivePowerBalance__ACBus")[1, "116"] ==
-          -0.3456209797192982
-    @test read_realized_variable(results_uc0, "ActivePowerBalance__ACBus")[1, "119"] ==
-          -0.6255739732919298
+    @test read_realized_variable(
+        results_uc0,
+        "ActivePowerBalance__ACBus";
+        table_format=TableFormat.WIDE,
+    )[
+        1,
+        "116",
+    ] == -0.3456209797192982
+    @test read_realized_variable(
+        results_uc0,
+        "ActivePowerBalance__ACBus";
+        table_format=TableFormat.WIDE,
+    )[
+        1,
+        "119",
+    ] == -0.6255739732919298
+
+    state_estimation_injection = read_realized_variable(
+        results_ucsub,
+        "StateEstimationInjections__ACBus";
+        table_format=TableFormat.WIDE,
+    )
+    active_power_balance = read_realized_variable(
+        results_uc0,
+        "ActivePowerBalance__ACBus";
+        table_format=TableFormat.WIDE,
+    )
 
-    # NOTE - open issue for results processing (this is only testing a single value at t=0): https://github.com/NREL-Sienna/PowerSimulations.jl/issues/1307
-    # We can manually go in and grab the results from the container to test all 5 values: 
-    param = sim.models.decision_models[2].internal.container.parameters
-    state_estimation_injection =
-        param[InfrastructureSystems.Optimization.ParameterKey{
-            PowerSimulationsDecomposition.StateEstimationInjections,
-            ACBus,
-        }(
-            "",
-        )].parameter_array
-    expr = sim.models.decision_models[1].internal.container.expressions
-    active_power_balance =
-        expr[InfrastructureSystems.Optimization.ExpressionKey{ActivePowerBalance, ACBus}(
-            "",
-        )]
     for b_number in [get_number(x) for x in get_components(ACBus, sys)]
-        apb = value.(active_power_balance[b_number, :]).data
-        sei = state_estimation_injection[string(b_number), :].data
+        apb = value.(active_power_balance[:, string(b_number)])
+        sei = state_estimation_injection[:, string(b_number)]
         @test isapprox(apb, sei)
     end
 end
 
 @testset "Vertical passing; compare branch models without emulator" begin
     sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
-    sys2 = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
     results_original, sim_original = run_rts_multi_stage_decomposition_simulation(
-        [sys, sys2];
+        sys;
         NT=5,
         mode="vertical",
         monitored_line_formulations=[StaticBranchUnbounded, StaticBranchUnbounded],
         use_emulator=false,
     )
     sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
-    sys2 = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
     results_se_line, sim_se_line = run_rts_multi_stage_decomposition_simulation(
-        [sys, sys2];
+        sys;
         NT=5,
         mode="vertical",
         monitored_line_formulations=[
@@ -74,34 +96,14 @@ end
     results_sub_se_line = get_decision_problem_results(results_se_line, "UC_Subsystem")
     flow_sub_original = read_realized_variable(
         results_sub_original,
-        "FlowActivePowerVariable__MonitoredLine",
+        "PTDFBranchFlow__MonitoredLine";
+        table_format=TableFormat.WIDE,
     )
     flow_sub_se_line = read_realized_variable(
         results_sub_se_line,
-        "FlowActivePowerVariable__MonitoredLine",
+        "PTDFBranchFlow__MonitoredLine";
+        table_format=TableFormat.WIDE,
     )
 
-    @test isapprox(flow_sub_original[1, "A28"], flow_sub_se_line[1, "A28"])
-
-    # NOTE - open issue for results processing (this is only testing a single value at t=0): https://github.com/NREL-Sienna/PowerSimulations.jl/issues/1307
-    # We can manually go in and grab the results from the container to test all 5 values: 
-    vars_original = sim_original.models.decision_models[2].internal.container.variables
-    flows_original = Vector(
-        vars_original[PowerSimulations.VariableKey{FlowActivePowerVariable, MonitoredLine}(
-            "",
-        )][
-            "A28",
-            :,
-        ],
-    )
-    vars_se_line = sim_se_line.models.decision_models[2].internal.container.variables
-    flows_se_line = Vector(
-        vars_se_line[PowerSimulations.VariableKey{FlowActivePowerVariable, MonitoredLine}(
-            "",
-        )][
-            "A28",
-            :,
-        ],
-    )
-    @test isapprox(flows_original, flows_se_line)
-
+    @test isapprox(flow_sub_original[:, "A28"], flow_sub_se_line[:, "A28"])
+end