diff --git a/src/PowerSimulations.jl b/src/PowerSimulations.jl
index e3a3eb9fb9..f4d882d98d 100644
--- a/src/PowerSimulations.jl
+++ b/src/PowerSimulations.jl
@@ -513,7 +513,6 @@ include("feedforward/feedforward_arguments.jl")
 include("feedforward/feedforward_constraints.jl")
 
 include("parameters/add_parameters.jl")
-include("parameters/update_parameters.jl")
 
 include("simulation/optimization_output_cache.jl")
 include("simulation/optimization_output_caches.jl")
@@ -534,6 +533,10 @@ include("simulation/simulation_internal.jl")
 include("simulation/simulation.jl")
 include("simulation/simulation_results_export.jl")
 include("simulation/simulation_results.jl")
+include("operation/operation_model_simulation_interface.jl")
+include("parameters/update_container_parameter_values.jl")
+include("parameters/update_cost_parameters.jl")
+include("parameters/update_parameters.jl")
 
 include("devices_models/devices/common/objective_function/common.jl")
 include("devices_models/devices/common/objective_function/linear_curve.jl")
diff --git a/src/core/parameters.jl b/src/core/parameters.jl
index 7e2e60afe2..2946dd68fe 100644
--- a/src/core/parameters.jl
+++ b/src/core/parameters.jl
@@ -6,6 +6,7 @@ struct TimeSeriesAttributes{T <: PSY.TimeSeriesData} <: ParameterAttributes
     name::String
     multiplier_id::Base.RefValue{Int}
     component_name_to_ts_uuid::Dict{String, String}
+    subsystem::Base.RefValue{String}
 end
 
 function TimeSeriesAttributes(
@@ -18,6 +19,7 @@ function TimeSeriesAttributes(
         name,
         Base.RefValue{Int}(multiplier_id),
         component_name_to_ts_uuid,
+        Base.RefValue{String}(""),
     )
 end
 
@@ -29,6 +31,13 @@ function set_time_series_multiplier_id!(attr::TimeSeriesAttributes, val::Int)
     return
 end
 
+get_subsystem(attr::TimeSeriesAttributes) = attr.subsystem[]
+function set_subsystem!(attr::TimeSeriesAttributes, val::String)
+    attr.subsystem[] = val
+    return
+end
+set_subsystem!(::TimeSeriesAttributes, ::Nothing) = nothing
+
 function add_component_name!(attr::TimeSeriesAttributes, name::String, uuid::String)
     if haskey(attr.component_name_to_ts_uuid, name)
         throw(ArgumentError("$name is already stored"))
diff --git a/src/core/settings.jl b/src/core/settings.jl
index d9dd095cbc..063ad1cce0 100644
--- a/src/core/settings.jl
+++ b/src/core/settings.jl
@@ -16,6 +16,7 @@ struct Settings
     export_pwl_vars::Bool
     allow_fails::Bool
     rebuild_model::Bool
+    export_optimization_model::Bool
     store_variable_names::Bool
     check_numerical_bounds::Bool
     ext::Dict{String, Any}
@@ -41,6 +42,7 @@ function Settings(
     allow_fails::Bool = false,
     check_numerical_bounds = true,
     rebuild_model = false,
+    export_optimization_model = false,
     store_variable_names = false,
     ext = Dict{String, Any}(),
 )
@@ -77,6 +79,7 @@ function Settings(
         export_pwl_vars,
         allow_fails,
         rebuild_model,
+        export_optimization_model,
         store_variable_names,
         check_numerical_bounds,
         ext,
@@ -151,6 +154,7 @@ get_detailed_optimizer_stats(settings::Settings) = settings.detailed_optimizer_s
 get_direct_mode_optimizer(settings::Settings) = settings.direct_mode_optimizer
 get_store_variable_names(settings::Settings) = settings.store_variable_names
 get_rebuild_model(settings::Settings) = settings.rebuild_model
+get_export_optimization_model(settings::Settings) = settings.export_optimization_model
 use_time_series_cache(settings::Settings) = settings.time_series_cache_size > 0
 
 function set_horizon!(settings::Settings, horizon::Dates.TimePeriod)
diff --git a/src/devices_models/device_constructors/branch_constructor.jl b/src/devices_models/device_constructors/branch_constructor.jl
index b1926fefbb..a4a517ce7d 100644
--- a/src/devices_models/device_constructors/branch_constructor.jl
+++ b/src/devices_models/device_constructors/branch_constructor.jl
@@ -10,9 +10,8 @@ function construct_device!(
         NetworkModel{AreaBalancePowerModel},
     },
 ) where {T <: PSY.ACBranch}
+    devices = get_available_components(model, sys)
     if has_subnetworks(network_model)
-        devices = get_available_components(model, sys)
-
         if get_use_slacks(model)
             add_variables!(
                 container,
@@ -46,6 +45,7 @@ function construct_device!(
             network_model,
         )
     end
+    add_feedforward_arguments!(container, model, devices)
     return
 end
 
@@ -59,6 +59,7 @@ function construct_device!(
         NetworkModel{AreaBalancePowerModel},
     },
 ) where {T <: PSY.ACBranch}
+    devices = get_available_components(model, sys)
     if has_subnetworks(network_model)
         devices =
             get_available_components(model, sys)
@@ -71,6 +72,7 @@ function construct_device!(
         )
         add_constraint_dual!(container, sys, model)
     end
+    add_feedforward_constraints!(container, model, devices)
     return
 end
 
@@ -87,9 +89,8 @@ function construct_device!(
     if get_use_slacks(model)
         throw(ArgumentError("StaticBranchBounds is not compatible with the use of slacks"))
     end
+    devices = get_available_components(model, sys)
     if has_subnetworks(network_model)
-        devices =
-            get_available_components(model, sys)
         add_variables!(
             container,
             FlowActivePowerVariable,
@@ -106,6 +107,7 @@ function construct_device!(
             network_model,
         )
     end
+    add_feedforward_arguments!(container, model, devices)
     return
 end
 
@@ -119,9 +121,8 @@ function construct_device!(
         NetworkModel{AreaBalancePowerModel},
     },
 ) where {T <: PSY.ACBranch}
+    devices = get_available_components(model, sys)
     if has_subnetworks(network_model)
-        devices =
-            get_available_components(model, sys)
         branch_rate_bounds!(
             container,
             devices,
@@ -129,6 +130,7 @@ function construct_device!(
             network_model,
         )
     end
+    add_feedforward_constraints!(container, model, devices)
     return
 end
 
@@ -142,9 +144,8 @@ function construct_device!(
         NetworkModel{AreaBalancePowerModel},
     },
 ) where {T <: PSY.ACBranch}
+    devices = get_available_components(model, sys)
     if has_subnetworks(network_model)
-        devices =
-            get_available_components(model, sys)
         add_variables!(
             container,
             FlowActivePowerVariable,
@@ -161,19 +162,22 @@ function construct_device!(
             network_model,
         )
     end
+    add_feedforward_arguments!(container, model, devices)
     return
 end
 
 function construct_device!(
     ::OptimizationContainer,
-    ::PSY.System,
+    sys::PSY.System,
     ::ModelConstructStage,
-    ::DeviceModel{<:PSY.ACBranch, StaticBranchUnbounded},
+    model::DeviceModel{<:PSY.ACBranch, StaticBranchUnbounded},
     network_model::Union{
         NetworkModel{CopperPlatePowerModel},
         NetworkModel{AreaBalancePowerModel},
     },
 )
+    devices = get_available_components(model, sys)
+    add_feedforward_constraints!(container, model, devices)
     return
 end
 
@@ -217,8 +221,7 @@ function construct_device!(
         )
     end
     add_feedforward_arguments!(container, device_model, devices)
-
-    add_feedforward_arguments!(container, device_model, devices)
+    return
 end
 
 # For DC Power only. Implements constraints
@@ -272,6 +275,7 @@ function construct_device!(
         devices,
         StaticBranch(),
     )
+    add_feedforward_arguments!(container, model, devices)
     return
 end
 
@@ -285,6 +289,7 @@ function construct_device!(
     devices = get_available_components(model, sys)
     add_constraints!(container, NetworkFlowConstraint, devices, model, network_model)
     add_constraints!(container, RateLimitConstraint, devices, model, network_model)
+    add_feedforward_constraints!(container, model, devices)
     objective_function!(container, devices, model, PTDFPowerModel)
     add_constraint_dual!(container, sys, model)
     return
@@ -310,6 +315,7 @@ function construct_device!(
         devices,
         StaticBranchBounds(),
     )
+    add_feedforward_arguments!(container, model, devices)
     return
 end
 
@@ -320,8 +326,7 @@ function construct_device!(
     model::DeviceModel{T, StaticBranchBounds},
     network_model::NetworkModel{<:AbstractPTDFModel},
 ) where {T <: PSY.ACBranch}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
     add_constraints!(container, NetworkFlowConstraint, devices, model, network_model)
     branch_rate_bounds!(
         container,
@@ -329,6 +334,7 @@ function construct_device!(
         model,
         network_model,
     )
+    add_feedforward_constraints!(container, model, devices)
     add_constraint_dual!(container, sys, model)
     return
 end
@@ -340,8 +346,7 @@ function construct_device!(
     model::DeviceModel{T, StaticBranchUnbounded},
     network_model::NetworkModel{<:AbstractPTDFModel},
 ) where {T <: PSY.ACBranch}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
     add_variables!(
         container,
         FlowActivePowerVariable,
@@ -349,6 +354,7 @@ function construct_device!(
         devices,
         StaticBranchUnbounded(),
     )
+    add_feedforward_arguments!(container, model, devices)
     return
 end
 
@@ -359,9 +365,8 @@ function construct_device!(
     model::DeviceModel{T, StaticBranchUnbounded},
     network_model::NetworkModel{<:AbstractPTDFModel},
 ) where {T <: PSY.ACBranch}
-    devices =
-        get_available_components(model, sys)
-
+    devices = get_available_components(model, sys)
+    add_feedforward_constraints!(container, model, devices)
     add_constraints!(container, NetworkFlowConstraint, devices, model, network_model)
     add_constraint_dual!(container, sys, model)
     return
@@ -386,7 +391,7 @@ function construct_device!(
             StaticBranch(),
         )
     end
-
+    add_feedforward_arguments!(container, device_model, devices)
     return
 end
 
@@ -397,10 +402,9 @@ function construct_device!(
     model::DeviceModel{T, StaticBranch},
     network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {T <: PSY.ACBranch}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
     branch_rate_bounds!(container, devices, model, network_model)
-
+    add_feedforward_constraints!(container, model, devices)
     add_constraints!(container, RateLimitConstraintFromTo, devices, model, network_model)
     add_constraints!(container, RateLimitConstraintToFrom, devices, model, network_model)
     add_constraint_dual!(container, sys, model)
@@ -448,6 +452,7 @@ function construct_device!(
         network_model,
     )
     add_constraint_dual!(container, sys, device_model)
+
     return
 end
 
@@ -544,8 +549,7 @@ function construct_device!(
     model::DeviceModel{T, HVDCTwoTerminalUnbounded},
     network_model::NetworkModel{<:AbstractPTDFModel},
 ) where {T <: TwoTerminalHVDCTypes}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
     add_variables!(container, FlowActivePowerVariable, devices, HVDCTwoTerminalUnbounded())
     add_to_expression!(
         container,
@@ -587,8 +591,7 @@ function construct_device!(
     model::DeviceModel{T, HVDCTwoTerminalLossless},
     network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {T <: TwoTerminalHVDCTypes}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
     add_constraints!(container, FlowRateConstraint, devices, model, network_model)
     add_constraint_dual!(container, sys, model)
     return
@@ -602,8 +605,7 @@ function construct_device!(
     model::DeviceModel{T, HVDCTwoTerminalLossless},
     network_model::NetworkModel{<:AbstractPTDFModel},
 ) where {T <: TwoTerminalHVDCTypes}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
     add_variables!(container, FlowActivePowerVariable, devices, HVDCTwoTerminalLossless())
     add_to_expression!(
         container,
@@ -626,8 +628,7 @@ function construct_device!(
 ) where {
     T <: TwoTerminalHVDCTypes,
 }
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
     add_constraints!(container, FlowRateConstraint, devices, model, network_model)
     add_constraint_dual!(container, sys, model)
     return
@@ -640,8 +641,7 @@ function construct_device!(
     model::DeviceModel{T, HVDCTwoTerminalDispatch},
     network_model::NetworkModel{<:AbstractPTDFModel},
 ) where {T <: TwoTerminalHVDCTypes}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
     add_variables!(
         container,
         FlowActivePowerToFromVariable,
@@ -690,8 +690,7 @@ function construct_device!(
     model::DeviceModel{T, HVDCTwoTerminalDispatch},
     network_model::NetworkModel{<:AbstractPTDFModel},
 ) where {T <: TwoTerminalHVDCTypes}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
     add_constraints!(container, FlowRateConstraintFromTo, devices, model, network_model)
     add_constraints!(container, FlowRateConstraintToFrom, devices, model, network_model)
     add_constraints!(container, HVDCPowerBalance, devices, model, network_model)
@@ -705,8 +704,7 @@ function construct_device!(
     model::DeviceModel{T, HVDCTwoTerminalDispatch},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {T <: TwoTerminalHVDCTypes}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
     add_variables!(
         container,
         FlowActivePowerToFromVariable,
@@ -747,8 +745,7 @@ function construct_device!(
     model::DeviceModel{T, HVDCTwoTerminalDispatch},
     network_model::NetworkModel{CopperPlatePowerModel},
 ) where {T <: TwoTerminalHVDCTypes}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
     @warn "CopperPlatePowerModel models with HVDC ignores inter-area losses"
     add_constraints!(container, FlowRateConstraintFromTo, devices, model, network_model)
     add_constraints!(container, FlowRateConstraintToFrom, devices, model, network_model)
@@ -892,8 +889,7 @@ function construct_device!(
     model::DeviceModel{T, HVDCTwoTerminalDispatch},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {T <: TwoTerminalHVDCTypes}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
     add_constraints!(container, FlowRateConstraintFromTo, devices, model, network_model)
     add_constraints!(container, FlowRateConstraintToFrom, devices, model, network_model)
     add_constraints!(container, HVDCPowerBalance, devices, model, network_model)
@@ -926,6 +922,7 @@ function construct_device!(
         model,
         network_model,
     )
+    add_feedforward_arguments!(container, model, devices)
     return
 end
 
@@ -950,6 +947,7 @@ function construct_device!(
         model,
         network_model,
     )
+    add_feedforward_arguments!(container, model, devices)
     return
 end
 
@@ -968,6 +966,7 @@ function construct_device!(
     add_constraints!(container, PhaseAngleControlLimit, devices, model, network_model)
     add_constraints!(container, NetworkFlowConstraint, devices, model, network_model)
     add_constraint_dual!(container, sys, model)
+    add_feedforward_constraints!(container, model, devices)
     return
 end
 
@@ -986,18 +985,32 @@ function construct_device!(
     add_constraints!(container, PhaseAngleControlLimit, devices, model, network_model)
     add_constraints!(container, NetworkFlowConstraint, devices, model, network_model)
     add_constraint_dual!(container, sys, model)
+    add_feedforward_constraints!(container, model, devices)
     return
 end
 
 ################################# AreaInterchange Models ################################
 function construct_device!(
-    ::OptimizationContainer,
-    ::PSY.System,
+    container::OptimizationContainer,
+    sys::PSY.System,
     ::ArgumentConstructStage,
     model::DeviceModel{PSY.AreaInterchange, U},
+    network_model::NetworkModel{CopperPlatePowerModel},
+) where {U <: Union{StaticBranchUnbounded, StaticBranch}}
+    devices = get_available_components(model, sys)
+    add_feedforward_arguments!(container, model, devices)
+    return
+end
+
+function construct_device!(
+    container::OptimizationContainer,
+    sys::PSY.System,
+    ::ModelConstructStage,
+    model::DeviceModel{PSY.AreaInterchange, StaticBranchUnbounded},
     network_model::NetworkModel{T},
-) where {T <: PM.AbstractPowerModel, U <: Union{StaticBranchUnbounded, StaticBranch}}
-    error("AreaInterchange is not yet implemented for $T")
+) where {T <: PM.AbstractActivePowerModel}
+    devices = get_available_components(model, sys)
+    add_feedforward_constraints!(container, model, devices)
     return
 end
 
@@ -1009,8 +1022,10 @@ function construct_device!(
     network_model::NetworkModel{U},
 ) where {
     T <: Union{StaticBranchUnbounded, StaticBranch},
-    U <: Union{AreaBalancePowerModel, AreaPTDFPowerModel},
+    U <: PM.AbstractActivePowerModel,
 }
+    devices = get_available_components(model, sys)
+    has_ts = PSY.has_time_series.(devices)
     if get_use_slacks(model)
         add_variables!(
             container,
@@ -1027,8 +1042,6 @@ function construct_device!(
             T(),
         )
     end
-    devices = get_available_components(model, sys)
-    has_ts = PSY.has_time_series.(devices)
     if any(has_ts) && !all(has_ts)
         error(
             "Not all AreaInterchange devices have time series. Check data to complete (or remove) time series.",
@@ -1062,6 +1075,7 @@ function construct_device!(
         add_parameters!(container, FromToFlowLimitParameter, devices, model)
         add_parameters!(container, ToFromFlowLimitParameter, devices, model)
     end
+    add_feedforward_arguments!(container, model, devices)
     return
 end
 
@@ -1071,15 +1085,16 @@ function construct_device!(
     ::ModelConstructStage,
     model::DeviceModel{PSY.AreaInterchange, StaticBranch},
     network_model::NetworkModel{T},
-) where {T <: AreaBalancePowerModel}
+) where {T <: PM.AbstractActivePowerModel}
     devices = get_available_components(model, sys)
     add_constraints!(container, FlowLimitConstraint, devices, model, network_model)
+    add_feedforward_constraints!(container, model, devices)
     return
 end
 
 function _get_branch_map(
     container::OptimizationContainer,
-    network_model::NetworkModel{AreaPTDFPowerModel},
+    network_model::NetworkModel,
     sys::PSY.System,
 )
     @assert !isempty(network_model.modeled_branch_types)
@@ -1122,7 +1137,7 @@ function construct_device!(
     ::ModelConstructStage,
     model::DeviceModel{PSY.AreaInterchange, StaticBranch},
     network_model::NetworkModel{T},
-) where {T <: AreaPTDFPowerModel}
+) where {T <: PSI.AbstractPTDFModel}
     devices = get_available_components(model, sys)
     add_constraints!(container, FlowLimitConstraint, devices, model, network_model)
     # Not ideal to do this here, but it is a not terrible workaround
@@ -1138,6 +1153,7 @@ function construct_device!(
         network_model,
         inter_area_branch_map,
     )
+    add_feedforward_constraints!(container, model, devices)
     return
 end
 
@@ -1148,6 +1164,7 @@ function construct_device!(
     model::DeviceModel{PSY.AreaInterchange, StaticBranchUnbounded},
     network_model::NetworkModel{AreaBalancePowerModel},
 )
+    add_feedforward_constraints!(container, model, devices)
     return
 end
 
@@ -1158,6 +1175,7 @@ function construct_device!(
     model::DeviceModel{PSY.AreaInterchange, StaticBranchUnbounded},
     network_model::NetworkModel{AreaPTDFPowerModel},
 )
+    devices = get_available_components(model, sys)
     inter_area_branch_map = _get_branch_map(container, network_model, sys)
     # Not ideal to do this here, but it is a not terrible workaround
     # The area interchanges are like a services/device mix.
@@ -1170,5 +1188,6 @@ function construct_device!(
         network_model,
         inter_area_branch_map,
     )
+    add_feedforward_constraints!(container, model, devices)
     return
 end
diff --git a/src/devices_models/device_constructors/renewablegeneration_constructor.jl b/src/devices_models/device_constructors/renewablegeneration_constructor.jl
index 7c05d30ceb..10939d52ea 100644
--- a/src/devices_models/device_constructors/renewablegeneration_constructor.jl
+++ b/src/devices_models/device_constructors/renewablegeneration_constructor.jl
@@ -8,8 +8,7 @@ function construct_device!(
     R <: PSY.RenewableGen,
     D <: AbstractRenewableDispatchFormulation,
 }
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_variables!(container, ActivePowerVariable, devices, D())
     add_variables!(container, ReactivePowerVariable, devices, D())
@@ -62,8 +61,7 @@ function construct_device!(
     model::DeviceModel{R, <:AbstractRenewableDispatchFormulation},
     network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {R <: PSY.RenewableGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     if has_service_model(model)
         add_constraints!(
@@ -119,8 +117,7 @@ function construct_device!(
     R <: PSY.RenewableGen,
     D <: AbstractRenewableDispatchFormulation,
 }
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_variables!(container, ActivePowerVariable, devices, D())
 
@@ -164,8 +161,7 @@ function construct_device!(
     model::DeviceModel{R, <:AbstractRenewableDispatchFormulation},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {R <: PSY.RenewableGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     if has_service_model(model)
         add_constraints!(
@@ -210,8 +206,7 @@ function construct_device!(
     model::DeviceModel{R, FixedOutput},
     network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {R <: PSY.RenewableGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)
     add_parameters!(container, ReactivePowerTimeSeriesParameter, devices, model)
@@ -242,8 +237,7 @@ function construct_device!(
     model::DeviceModel{R, FixedOutput},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {R <: PSY.RenewableGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)
     add_to_expression!(
diff --git a/src/devices_models/device_constructors/thermalgeneration_constructor.jl b/src/devices_models/device_constructors/thermalgeneration_constructor.jl
index f24d154b08..5225586c43 100644
--- a/src/devices_models/device_constructors/thermalgeneration_constructor.jl
+++ b/src/devices_models/device_constructors/thermalgeneration_constructor.jl
@@ -5,8 +5,7 @@ function construct_device!(
     device_model::DeviceModel{T, FixedOutput},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
     add_parameters!(container, ActivePowerTimeSeriesParameter, devices, device_model)
     add_to_expression!(
         container,
@@ -44,8 +43,7 @@ function construct_device!(
     T <: PSY.ThermalGen,
     D <: AbstractStandardUnitCommitment,
 }
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_variables!(container, ActivePowerVariable, devices, D())
     add_variables!(container, ReactivePowerVariable, devices, D())
@@ -107,8 +105,7 @@ function construct_device!(
     model::DeviceModel{T, <:AbstractStandardUnitCommitment},
     network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_constraints!(
         container,
@@ -155,8 +152,7 @@ function construct_device!(
     model::DeviceModel{T, D},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {T <: PSY.ThermalGen, D <: AbstractStandardUnitCommitment}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_variables!(container, ActivePowerVariable, devices, D())
     add_variables!(container, OnVariable, devices, D())
@@ -210,8 +206,7 @@ function construct_device!(
     model::DeviceModel{T, <:AbstractStandardUnitCommitment},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
     add_constraints!(
         container,
         ActivePowerVariableLimitsConstraint,
@@ -251,8 +246,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalBasicUnitCommitment},
     network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_variables!(container, ActivePowerVariable, devices, ThermalBasicUnitCommitment())
     add_variables!(container, ReactivePowerVariable, devices, ThermalBasicUnitCommitment())
@@ -312,8 +306,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalBasicUnitCommitment},
     network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_constraints!(
         container,
@@ -359,8 +352,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalBasicUnitCommitment},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_variables!(container, ActivePowerVariable, devices, ThermalBasicUnitCommitment())
     add_variables!(container, OnVariable, devices, ThermalBasicUnitCommitment())
@@ -411,8 +403,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalBasicUnitCommitment},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_constraints!(
         container,
@@ -450,8 +441,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalStandardDispatch},
     network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_variables!(container, ActivePowerVariable, devices, ThermalStandardDispatch())
     add_variables!(container, ReactivePowerVariable, devices, ThermalStandardDispatch())
@@ -508,8 +498,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalStandardDispatch},
     network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_constraints!(
         container,
@@ -555,8 +544,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalStandardDispatch},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_variables!(container, ActivePowerVariable, devices, ThermalStandardDispatch())
 
@@ -604,8 +592,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalStandardDispatch},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_constraints!(
         container,
@@ -643,8 +630,7 @@ function construct_device!(
     T <: PSY.ThermalGen,
     D <: AbstractThermalDispatchFormulation,
 }
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_variables!(container, ActivePowerVariable, devices, D())
     add_variables!(container, ReactivePowerVariable, devices, D())
@@ -696,8 +682,7 @@ function construct_device!(
     model::DeviceModel{T, <:AbstractThermalDispatchFormulation},
     network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_constraints!(
         container,
@@ -742,8 +727,7 @@ function construct_device!(
     T <: PSY.ThermalGen,
     D <: AbstractThermalDispatchFormulation,
 }
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_variables!(container, ActivePowerVariable, devices, D())
 
@@ -786,8 +770,7 @@ function construct_device!(
     model::DeviceModel{T, <:AbstractThermalDispatchFormulation},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_constraints!(
         container,
@@ -1085,8 +1068,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalCompactUnitCommitment},
     network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_variables!(
         container,
@@ -1157,8 +1139,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalCompactUnitCommitment},
     network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_constraints!(
         container,
@@ -1203,8 +1184,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalCompactUnitCommitment},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_variables!(
         container,
@@ -1268,8 +1248,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalCompactUnitCommitment},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_constraints!(
         container,
@@ -1306,8 +1285,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalBasicCompactUnitCommitment},
     network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_variables!(
         container,
@@ -1376,8 +1354,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalBasicCompactUnitCommitment},
     network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_constraints!(
         container,
@@ -1420,8 +1397,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalBasicCompactUnitCommitment},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_variables!(
         container,
@@ -1483,8 +1459,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalBasicCompactUnitCommitment},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_constraints!(
         container,
@@ -1519,8 +1494,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalCompactDispatch},
     network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_variables!(container, PowerAboveMinimumVariable, devices, ThermalCompactDispatch())
     add_variables!(container, ReactivePowerVariable, devices, ThermalCompactDispatch())
@@ -1560,23 +1534,6 @@ function construct_device!(
         model,
         network_model,
     )
-    add_to_expression!(
-        container,
-        ReactivePowerBalance,
-        ReactivePowerVariable,
-        devices,
-        model,
-        network_model,
-    )
-    add_to_expression!(
-        container,
-        ActivePowerBalance,
-        OnStatusParameter,
-        devices,
-        model,
-        network_model,
-    )
-
     add_to_expression!(
         container,
         ActivePowerRangeExpressionLB,
@@ -1603,8 +1560,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalCompactDispatch},
     network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_constraints!(
         container,
@@ -1647,8 +1603,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalCompactDispatch},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_variables!(container, PowerAboveMinimumVariable, devices, ThermalCompactDispatch())
 
@@ -1706,8 +1661,7 @@ function construct_device!(
     model::DeviceModel{T, ThermalCompactDispatch},
     network_model::NetworkModel{<:PM.AbstractActivePowerModel},
 ) where {T <: PSY.ThermalGen}
-    devices =
-        get_available_components(model, sys)
+    devices = get_available_components(model, sys)
 
     add_constraints!(
         container,
diff --git a/src/devices_models/devices/AC_branches.jl b/src/devices_models/devices/AC_branches.jl
index 8046f17242..48767cd794 100644
--- a/src/devices_models/devices/AC_branches.jl
+++ b/src/devices_models/devices/AC_branches.jl
@@ -21,9 +21,14 @@
 get_variable_binary(::FlowActivePowerVariable, ::Type{<:PSY.ACBranch}, ::AbstractBranchFormulation,) = false
 get_variable_binary(::PhaseShifterAngle, ::Type{PSY.PhaseShiftingTransformer}, ::AbstractBranchFormulation,) = false
 
-get_parameter_multiplier(::FixValueParameter, ::PSY.ACBranch, ::StaticBranch) = 1.0
+get_parameter_multiplier(::FixValueParameter, ::PSY.ACBranch, ::AbstractBranchFormulation) = 1.0
+get_parameter_multiplier(::LowerBoundValueParameter, ::PSY.ACBranch, ::AbstractBranchFormulation) = 1.0
+get_parameter_multiplier(::UpperBoundValueParameter, ::PSY.ACBranch, ::AbstractBranchFormulation) = 1.0
+
 get_variable_multiplier(::PhaseShifterAngle, d::PSY.PhaseShiftingTransformer, ::PhaseAngleControl) = 1.0/PSY.get_x(d)
 
+
+
 get_initial_conditions_device_model(::OperationModel, ::DeviceModel{T, U}) where {T <: PSY.ACBranch, U <: AbstractBranchFormulation} = DeviceModel(T, U)
 
 #### Properties of slack variables
diff --git a/src/devices_models/devices/area_interchange.jl b/src/devices_models/devices/area_interchange.jl
index 4233f9c7c0..84c953e54b 100644
--- a/src/devices_models/devices/area_interchange.jl
+++ b/src/devices_models/devices/area_interchange.jl
@@ -1,6 +1,16 @@
 #! format: off
 get_multiplier_value(::FromToFlowLimitParameter, d::PSY.AreaInterchange, ::AbstractBranchFormulation) = -1.0 * PSY.get_from_to_flow_limit(d)
 get_multiplier_value(::ToFromFlowLimitParameter, d::PSY.AreaInterchange, ::AbstractBranchFormulation) = PSY.get_to_from_flow_limit(d)
+
+get_parameter_multiplier(::FixValueParameter, ::PSY.AreaInterchange, ::AbstractBranchFormulation) = 1.0
+get_parameter_multiplier(::LowerBoundValueParameter, ::PSY.AreaInterchange, ::AbstractBranchFormulation) = 1.0
+get_parameter_multiplier(::UpperBoundValueParameter, ::PSY.AreaInterchange, ::AbstractBranchFormulation) = 1.0
+
+get_initial_conditions_device_model(
+    ::OperationModel,
+    model::DeviceModel{PSY.AreaInterchange, T},
+) where {T <: AbstractBranchFormulation} = DeviceModel(PSY.AreaInterchange, T)
+
 #! format: on
 
 function get_default_time_series_names(
@@ -26,7 +36,7 @@ function add_variables!(
     model::NetworkModel{T},
     devices::IS.FlattenIteratorWrapper{PSY.AreaInterchange},
     formulation::AbstractBranchFormulation,
-) where {T <: Union{AreaBalancePowerModel, AreaPTDFPowerModel}}
+) where {T <: PM.AbstractPowerModel}
     time_steps = get_time_steps(container)
 
     variable = add_variable_container!(
@@ -47,6 +57,18 @@ function add_variables!(
     return
 end
 
+function add_variables!(
+    container::OptimizationContainer,
+    ::Type{FlowActivePowerVariable},
+    model::NetworkModel{CopperPlatePowerModel},
+    devices::IS.FlattenIteratorWrapper{PSY.AreaInterchange},
+    formulation::AbstractBranchFormulation,
+)
+    @warn(
+        "CopperPlatePowerModel ignores AreaInterchanges. Instead use AreaBalancePowerModel."
+    )
+end
+
 """
 Add flow constraints for area interchanges
 """
@@ -56,7 +78,7 @@ function add_constraints!(
     devices::IS.FlattenIteratorWrapper{PSY.AreaInterchange},
     model::DeviceModel{PSY.AreaInterchange, StaticBranch},
     ::NetworkModel{T},
-) where {T <: Union{AreaBalancePowerModel, AreaPTDFPowerModel}}
+) where {T <: PM.AbstractActivePowerModel}
     time_steps = get_time_steps(container)
     device_names = [PSY.get_name(d) for d in devices]
 
@@ -139,13 +161,14 @@ function add_constraints!(
     ::Type{LineFlowBoundConstraint},
     devices::IS.FlattenIteratorWrapper{PSY.AreaInterchange},
     model::DeviceModel{PSY.AreaInterchange, <:AbstractBranchFormulation},
-    network_model::NetworkModel{AreaPTDFPowerModel},
+    network_model::NetworkModel{T},
     inter_area_branch_map::Dict{
         Tuple{PSY.Area, PSY.Area},
         Dict{DataType, Vector{<:PSY.ACBranch}},
     },
-)
+) where {T <: AbstractPTDFModel}
     @assert !isempty(inter_area_branch_map)
+
     time_steps = get_time_steps(container)
     device_names = [PSY.get_name(d) for d in devices]
 
@@ -173,13 +196,16 @@ function add_constraints!(
         inter_change_name = PSY.get_name(area_interchange)
         area_from = PSY.get_from_area(area_interchange)
         area_to = PSY.get_to_area(area_interchange)
+        direction_branch_map = Dict{Float64, Dict{DataType, Vector{<:PSY.ACBranch}}}()
         if haskey(inter_area_branch_map, (area_from, area_to))
-            inter_area_branches = inter_area_branch_map[(area_from, area_to)]
-            mult = 1.0
-        elseif haskey(inter_area_branch_map, (area_to, area_from))
-            inter_area_branches = inter_area_branch_map[(area_to, area_from)]
-            mult = -1.0
-        else
+            # 1 is the multiplier
+            direction_branch_map[1.0] = inter_area_branch_map[(area_from, area_to)]
+        end
+        if haskey(inter_area_branch_map, (area_to, area_from))
+            # -1 is the multiplier because the direction is reversed
+            direction_branch_map[-1.0] = inter_area_branch_map[(area_to, area_from)]
+        end
+        if isempty(direction_branch_map)
             @warn(
                 "There are no branches modeled in Area InterChange $(summary(area_interchange)) \
           LineFlowBoundConstraint not created"
@@ -189,11 +215,13 @@ function add_constraints!(
 
         for t in time_steps
             sum_of_flows = JuMP.AffExpr()
-            for (type, branches) in inter_area_branches
-                flow_vars = get_variable(container, FlowActivePowerVariable(), type)
-                for b in branches
-                    b_name = PSY.get_name(b)
-                    _add_to_jump_expression!(sum_of_flows, flow_vars[b_name, t], mult)
+            for (mult, inter_area_branches) in direction_branch_map
+                for (type, branches) in inter_area_branches
+                    flow_vars = get_variable(container, FlowActivePowerVariable(), type)
+                    for b in branches
+                        b_name = PSY.get_name(b)
+                        _add_to_jump_expression!(sum_of_flows, flow_vars[b_name, t], mult)
+                    end
                 end
             end
             con_ub[inter_change_name, t] =
diff --git a/src/devices_models/devices/common/add_to_expression.jl b/src/devices_models/devices/common/add_to_expression.jl
index 51c33b245b..93978aa22c 100644
--- a/src/devices_models/devices/common/add_to_expression.jl
+++ b/src/devices_models/devices/common/add_to_expression.jl
@@ -1011,7 +1011,7 @@ function add_to_expression!(
     network_model::NetworkModel{U},
 ) where {
     T <: ActivePowerBalance,
-    U <: Union{AreaBalancePowerModel, AreaPTDFPowerModel},
+    U <: PM.AbstractActivePowerModel,
     W <: AbstractBranchFormulation,
 }
     flow_variable = get_variable(container, FlowActivePowerVariable(), PSY.AreaInterchange)
diff --git a/src/devices_models/devices/common/objective_function/piecewise_linear.jl b/src/devices_models/devices/common/objective_function/piecewise_linear.jl
index 37ae04b62d..bc184fb972 100644
--- a/src/devices_models/devices/common/objective_function/piecewise_linear.jl
+++ b/src/devices_models/devices/common/objective_function/piecewise_linear.jl
@@ -408,13 +408,25 @@ function _add_pwl_term!(
     end
 
     # Compact PWL data does not exists anymore
-
-    if slopes[1] != 0.0
-        @debug "PWL has no 0.0 intercept for generator $(component_name)"
-        # adds a first intercept a x = 0.0 and y below the intercept of the first tuple to make convex equivalent
-        intercept_point = (x = 0.0, y = first(data).y - COST_EPSILON)
-        data = PSY.PiecewiseLinearData(vcat(intercept_point, get_points(data)))
-        @assert PSY.is_convex(slopes)
+    x_coords = PSY.get_x_coords(data)
+    if x_coords[1] != 0.0
+        y_coords = PSY.get_y_coords(data)
+        x_first = round(x_coords[1]; digits = 3)
+        y_first = round(y_coords[1]; digits = 3)
+        slope_first = round(slopes[1]; digits = 3)
+        guess_y_zero = y_coords[1] - slopes[1] * x_coords[1]
+        @warn(
+            "PWL has no 0.0 intercept for generator $(name). First point is given at (x = $(x_first), y = $(y_first)). Adding a first intercept at (x = 0.0, y = $(round(guess_y_zero, digits = 3)) to have equal initial slope $(slope_first)"
+        )
+        if guess_y_zero < 0.0
+            error(
+                "Added zero intercept has negative cost for generator $(name). Consider using other formulation or improve data.",
+            )
+        end
+        # adds a first intercept a x = 0.0 and y above the intercept of the first tuple to make convex equivalent (avoid floating point issues of almost equal slopes)
+        intercept_point = (x = 0.0, y = guess_y_zero + COST_EPSILON)
+        data = PSY.PiecewiseLinearData(vcat(intercept_point, PSY.get_points(data)))
+        @assert PSY.is_convex(data)
     end
 
     time_steps = get_time_steps(container)
@@ -422,7 +434,7 @@ function _add_pwl_term!(
     break_points = PSY.get_x_coords(data)
     sos_val = _get_sos_value(container, V, component)
     for t in time_steps
-        _add_pwl_variables!(container, T, component_name, t, data)
+        _add_pwl_variables!(container, T, name, t, data)
         _add_pwl_constraint!(container, component, U(), break_points, sos_val, t)
         pwl_cost =
             _get_pwl_cost_expression(container, component, t, cost_function, U(), V())
diff --git a/src/devices_models/devices/electric_loads.jl b/src/devices_models/devices/electric_loads.jl
index 2452b18f7d..6725b2bd3e 100644
--- a/src/devices_models/devices/electric_loads.jl
+++ b/src/devices_models/devices/electric_loads.jl
@@ -7,14 +7,14 @@ get_variable_multiplier(_, ::Type{<:PSY.ElectricLoad}, ::AbstractLoadFormulation
 
 get_variable_binary(::ActivePowerVariable, ::Type{<:PSY.ElectricLoad}, ::AbstractLoadFormulation) = false
 get_variable_lower_bound(::ActivePowerVariable, d::PSY.ElectricLoad, ::AbstractLoadFormulation) = 0.0
-get_variable_upper_bound(::ActivePowerVariable, d::PSY.ElectricLoad, ::AbstractLoadFormulation) = PSY.get_active_power(d)
+get_variable_upper_bound(::ActivePowerVariable, d::PSY.ElectricLoad, ::AbstractLoadFormulation) = PSY.get_max_active_power(d)
 
 ########################### ReactivePowerVariable, ElectricLoad ####################################
 
 get_variable_binary(::ReactivePowerVariable, ::Type{<:PSY.ElectricLoad}, ::AbstractLoadFormulation) = false
 
 get_variable_lower_bound(::ReactivePowerVariable, d::PSY.ElectricLoad, ::AbstractLoadFormulation) = 0.0
-get_variable_upper_bound(::ReactivePowerVariable, d::PSY.ElectricLoad, ::AbstractLoadFormulation) = PSY.get_reactive_power(d)
+get_variable_upper_bound(::ReactivePowerVariable, d::PSY.ElectricLoad, ::AbstractLoadFormulation) = PSY.get_max_reactive_power(d)
 
 ########################### ReactivePowerVariable, ElectricLoad ####################################
 
diff --git a/src/feedforward/feedforward_arguments.jl b/src/feedforward/feedforward_arguments.jl
index 20fdd57767..25e40081b9 100644
--- a/src/feedforward/feedforward_arguments.jl
+++ b/src/feedforward/feedforward_arguments.jl
@@ -109,9 +109,10 @@ function _add_feedforward_slack_variables!(
         var_type = get_entry_type(var)
         variable = add_variable_container!(
             container,
-            T,
+            T(),
             U,
-            [PSY.get_name(d) for d in devices];
+            [PSY.get_name(d) for d in devices],
+            time_steps;
             meta = "$(var_type)",
         )
 
@@ -135,11 +136,12 @@ function _add_feedforward_arguments!(
     parameter_type = get_default_parameter_type(ff, T)
     add_parameters!(container, parameter_type, ff, model, devices)
     if get_slacks(ff)
-        add_feedforward_slack_variables!(
+        _add_feedforward_slack_variables!(
             container,
-            UpperBoundFeedForwardSlack,
-            devices,
+            UpperBoundFeedForwardSlack(),
+            ff,
             model,
+            devices,
         )
     end
     return
@@ -154,11 +156,12 @@ function _add_feedforward_arguments!(
     parameter_type = get_default_parameter_type(ff, SR)
     add_parameters!(container, parameter_type, ff, model, contributing_devices)
     if get_slacks(ff)
-        add_feedforward_slack_variables!(
+        _add_feedforward_slack_variables!(
             container,
-            UpperBoundFeedForwardSlack,
-            contributing_devices,
+            UpperBoundFeedForwardSlack(),
+            ff,
             model,
+            contributing_devices,
         )
     end
     return
@@ -175,7 +178,7 @@ function _add_feedforward_arguments!(
     if get_slacks(ff)
         _add_feedforward_slack_variables!(
             container,
-            LowerBoundFeedForwardSlack,
+            LowerBoundFeedForwardSlack(),
             ff,
             model,
             devices,
diff --git a/src/feedforward/feedforward_constraints.jl b/src/feedforward/feedforward_constraints.jl
index 2acac1dfb7..a85333b85f 100644
--- a/src/feedforward/feedforward_constraints.jl
+++ b/src/feedforward/feedforward_constraints.jl
@@ -456,7 +456,6 @@ function add_feedforward_constraints!(
     devices::Union{Vector{T}, IS.FlattenIteratorWrapper{T}},
     ff::FixValueFeedforward,
 ) where {T <: PSY.Component}
-    time_steps = get_time_steps(container)
     parameter_type = get_default_parameter_type(ff, T)
     source_key = get_optimization_container_key(ff)
     var_type = get_entry_type(source_key)
@@ -466,7 +465,6 @@ function add_feedforward_constraints!(
         variable = get_variable(container, var)
         set_name, set_time = JuMP.axes(variable)
         IS.@assert_op set_name == [PSY.get_name(d) for d in devices]
-        #IS.@assert_op set_time == time_steps
 
         for t in set_time, name in set_name
             JuMP.fix(variable[name, t], param[name, t] * multiplier[name, t]; force = true)
diff --git a/src/initial_conditions/initialization.jl b/src/initial_conditions/initialization.jl
index e375b25c12..cff6e81ed0 100644
--- a/src/initial_conditions/initialization.jl
+++ b/src/initial_conditions/initialization.jl
@@ -13,10 +13,12 @@ function get_initial_conditions_template(model::OperationModel)
         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)
         network_model.bus_area_map = get_bus_area_map(get_network_model(model.template))
     end
-
+    network_model.modeled_branch_types =
+        get_network_model(model.template).modeled_branch_types
     ic_template = ProblemTemplate(network_model)
     for device_model in values(model.template.devices)
         base_model = get_initial_conditions_device_model(model, device_model)
diff --git a/src/operation/decision_model.jl b/src/operation/decision_model.jl
index b0718de645..785742ac72 100644
--- a/src/operation/decision_model.jl
+++ b/src/operation/decision_model.jl
@@ -114,6 +114,7 @@ function DecisionModel{M}(
     direct_mode_optimizer = false,
     store_variable_names = false,
     rebuild_model = false,
+    export_optimization_model = false,
     check_numerical_bounds = true,
     initial_time = UNSET_INI_TIME,
     time_series_cache_size::Int = IS.TIME_SERIES_CACHE_SIZE_BYTES,
@@ -139,6 +140,7 @@ function DecisionModel{M}(
         check_numerical_bounds = check_numerical_bounds,
         store_variable_names = store_variable_names,
         rebuild_model = rebuild_model,
+        export_optimization_model = export_optimization_model,
     )
     return DecisionModel{M}(template, sys, settings, jump_model; name = name)
 end
@@ -444,7 +446,7 @@ keyword arguments to that function.
   - `console_level = Logging.Error`:
   - `file_level = Logging.Info`:
   - `disable_timer_outputs = false` : Enable/Disable timing outputs
-  - `serialize::Bool = true`: If true, serialize the model to a file to allow re-execution later.
+  - `export_optimization_problem::Bool = true`: If true, serialize the model to a file to allow re-execution later.
 
 # Examples
 
@@ -459,7 +461,7 @@ function solve!(
     console_level = Logging.Error,
     file_level = Logging.Info,
     disable_timer_outputs = false,
-    serialize = true,
+    export_optimization_problem = true,
     kwargs...,
 )
     build_if_not_already_built!(
@@ -500,7 +502,7 @@ function solve!(
                         current_time,
                     )
                 end
-                if serialize
+                if export_optimization_problem
                     TimerOutputs.@timeit RUN_OPERATION_MODEL_TIMER "Serialize" begin
                         serialize_problem(model; optimizer = optimizer)
                         serialize_optimization_model(model)
@@ -558,19 +560,3 @@ function solve!(
     end
     return get_run_status(model)
 end
-
-function update_parameters!(
-    model::DecisionModel,
-    decision_states::DatasetContainer{InMemoryDataset},
-)
-    cost_function_unsynch(get_optimization_container(model))
-    for key in keys(get_parameters(model))
-        update_parameter_values!(model, key, decision_states)
-    end
-    if !is_synchronized(model)
-        update_objective_function!(get_optimization_container(model))
-        obj_func = get_objective_expression(get_optimization_container(model))
-        set_synchronized_status!(obj_func, true)
-    end
-    return
-end
diff --git a/src/operation/emulation_model.jl b/src/operation/emulation_model.jl
index c635cf35b6..cd42d60ec9 100644
--- a/src/operation/emulation_model.jl
+++ b/src/operation/emulation_model.jl
@@ -454,6 +454,30 @@ function update_model!(
     return
 end
 
+"""
+Update parameter function an OperationModel
+"""
+function update_parameter_values!(
+    model::EmulationModel,
+    key::ParameterKey{T, U},
+    input::DatasetContainer{InMemoryDataset},
+) where {T <: ParameterType, U <: PSY.Component}
+    # Enable again for detailed debugging
+    # TimerOutputs.@timeit RUN_SIMULATION_TIMER "$T $U Parameter Update" begin
+    optimization_container = get_optimization_container(model)
+    update_container_parameter_values!(optimization_container, model, key, input)
+    parameter_attributes = get_parameter_attributes(optimization_container, key)
+    IS.@record :execution ParameterUpdateEvent(
+        T,
+        U,
+        parameter_attributes,
+        get_current_timestamp(model),
+        get_name(model),
+    )
+    #end
+    return
+end
+
 function update_model!(model::EmulationModel)
     update_model!(model, get_store(model), InterProblemChronology())
     return
@@ -507,7 +531,7 @@ keyword arguments to that function.
   - `export_problem_results::Bool`: If true, export OptimizationProblemResults DataFrames to CSV files.
   - `output_dir::String`: Required if the model is not already built, otherwise ignored
   - `enable_progress_bar::Bool`: Enables/Disable progress bar printing
-  - `serialize::Bool`: If true, serialize the model to a file to allow re-execution later.
+  - `export_optimization_model::Bool`: If true, serialize the model to a file to allow re-execution later.
 
 # Examples
 
@@ -522,7 +546,7 @@ function run!(
     console_level = Logging.Error,
     file_level = Logging.Info,
     disable_timer_outputs = false,
-    serialize = true,
+    export_optimization_model = true,
     kwargs...,
 )
     build_if_not_already_built!(
@@ -555,7 +579,7 @@ function run!(
                     run_impl!(model; kwargs...)
                     set_run_status!(model, RunStatus.SUCCESSFULLY_FINALIZED)
                 end
-                if serialize
+                if export_optimization_model
                     TimerOutputs.@timeit RUN_OPERATION_MODEL_TIMER "Serialize" begin
                         optimizer = get(kwargs, :optimizer, nothing)
                         serialize_problem(model; optimizer = optimizer)
diff --git a/src/operation/operation_model_interface.jl b/src/operation/operation_model_interface.jl
index 1eb0e8ecc3..d2157ee967 100644
--- a/src/operation/operation_model_interface.jl
+++ b/src/operation/operation_model_interface.jl
@@ -104,13 +104,22 @@ end
 
 function solve_impl!(model::OperationModel)
     container = get_optimization_container(model)
+    model_name = get_name(model)
+    ts = get_current_timestamp(model)
+    output_dir = get_output_dir(model)
+
+    if get_export_optimization_model(get_settings(model))
+        model_output_dir = joinpath(output_dir, "optimization_model_exports")
+        mkpath(model_output_dir)
+        tss = replace("$(ts)", ":" => "_")
+        model_export_path = joinpath(model_output_dir, "exported_$(model_name)_$(tss).json")
+        serialize_optimization_model(container, model_export_path)
+    end
+
     status = solve_impl!(container, get_system(model))
     set_run_status!(model, status)
     if status != RunStatus.SUCCESSFULLY_FINALIZED
         settings = get_settings(model)
-        model_name = get_name(model)
-        ts = get_current_timestamp(model)
-        output_dir = get_output_dir(model)
         infeasible_opt_path = joinpath(output_dir, "infeasible_$(model_name).json")
         @error("Serializing Infeasible Problem at $(infeasible_opt_path)")
         serialize_optimization_model(container, infeasible_opt_path)
@@ -432,16 +441,6 @@ function serialize_optimization_model(model::OperationModel)
     return
 end
 
-function update_model!(model::OperationModel, source, ini_cond_chronology)
-    TimerOutputs.@timeit RUN_SIMULATION_TIMER "Parameter Updates" begin
-        update_parameters!(model, get_decision_states(source))
-    end
-    TimerOutputs.@timeit RUN_SIMULATION_TIMER "Ini Cond Updates" begin
-        update_initial_conditions!(model, source, ini_cond_chronology)
-    end
-    return
-end
-
 function instantiate_network_model(model::OperationModel)
     template = get_template(model)
     network_model = get_network_model(template)
diff --git a/src/operation/operation_model_simulation_interface.jl b/src/operation/operation_model_simulation_interface.jl
new file mode 100644
index 0000000000..3ec50b2c92
--- /dev/null
+++ b/src/operation/operation_model_simulation_interface.jl
@@ -0,0 +1,31 @@
+function update_model!(model::OperationModel, source::SimulationState, ini_cond_chronology)
+    TimerOutputs.@timeit RUN_SIMULATION_TIMER "Parameter Updates" begin
+        update_parameters!(model, source)
+    end
+    TimerOutputs.@timeit RUN_SIMULATION_TIMER "Ini Cond Updates" begin
+        update_initial_conditions!(model, source, ini_cond_chronology)
+    end
+    return
+end
+
+function update_parameters!(model::EmulationModel, state::SimulationState)
+    data = get_decision_states(state)
+    update_parameters!(model, data)
+    return
+end
+
+function update_parameters!(
+    model::DecisionModel,
+    simulation_state::SimulationState,
+)
+    cost_function_unsynch(get_optimization_container(model))
+    for key in keys(get_parameters(model))
+        update_parameter_values!(model, key, simulation_state)
+    end
+    if !is_synchronized(model)
+        update_objective_function!(get_optimization_container(model))
+        obj_func = get_objective_expression(get_optimization_container(model))
+        set_synchronized_status!(obj_func, true)
+    end
+    return
+end
diff --git a/src/operation/time_series_interface.jl b/src/operation/time_series_interface.jl
index 040f7b2e1d..990e6dff7f 100644
--- a/src/operation/time_series_interface.jl
+++ b/src/operation/time_series_interface.jl
@@ -2,9 +2,9 @@ function get_time_series_values!(
     time_series_type::Type{T},
     model::DecisionModel,
     component,
-    name,
+    name::String,
     multiplier_id::Int,
-    initial_time,
+    initial_time::Dates.DateTime,
     horizon::Int;
     ignore_scaling_factors = true,
 ) where {T <: PSY.Forecast}
diff --git a/src/parameters/add_parameters.jl b/src/parameters/add_parameters.jl
index b581893ed0..6261ded49c 100644
--- a/src/parameters/add_parameters.jl
+++ b/src/parameters/add_parameters.jl
@@ -216,6 +216,7 @@ function _add_time_series_parameters!(
         time_steps,
     )
     set_time_series_multiplier_id!(get_attributes(param_container), time_series_mult_id)
+    set_subsystem!(get_attributes(param_container), get_subsystem(model))
     jump_model = get_jump_model(container)
 
     for (ts_uuid, ts_values) in initial_values
@@ -282,6 +283,7 @@ function _add_parameters!(
     )
 
     set_time_series_multiplier_id!(get_attributes(parameter_container), time_series_mult_id)
+    set_subsystem!(get_attributes(parameter_container), get_subsystem(model))
     jump_model = get_jump_model(container)
     ts_vector = get_time_series(container, service, T(), name)
     multiplier = get_multiplier_value(T(), service, V())
diff --git a/src/parameters/update_container_parameter_values.jl b/src/parameters/update_container_parameter_values.jl
new file mode 100644
index 0000000000..5784c0964f
--- /dev/null
+++ b/src/parameters/update_container_parameter_values.jl
@@ -0,0 +1,455 @@
+function _update_parameter_values!(
+    ::AbstractArray{T},
+    ::NoAttributes,
+    args...,
+) where {T <: Union{Float64, JuMP.VariableRef}} end
+
+######################## Methods to update Parameters from Time Series #####################
+function _set_param_value!(param::JuMPVariableMatrix, value::Float64, name::String, t::Int)
+    fix_parameter_value(param[name, t], value)
+    return
+end
+
+function _set_param_value!(
+    param::DenseAxisArray{Vector{NTuple{2, Float64}}},
+    value::Vector{NTuple{2, Float64}},
+    name::String,
+    t::Int,
+)
+    param[name, t] = value
+    return
+end
+
+function _set_param_value!(param::JuMPFloatArray, value::Float64, name::String, t::Int)
+    param[name, t] = value
+    return
+end
+
+function _update_parameter_values!(
+    parameter_array::AbstractArray{T},
+    attributes::TimeSeriesAttributes{U},
+    ::Type{V},
+    model::DecisionModel,
+    ::DatasetContainer{InMemoryDataset},
+) where {
+    T <: Union{JuMP.VariableRef, Float64},
+    U <: PSY.AbstractDeterministic,
+    V <: PSY.Component,
+}
+    initial_forecast_time = get_current_time(model) # Function not well defined for DecisionModels
+    horizon = get_time_steps(get_optimization_container(model))[end]
+    ts_name = get_time_series_name(attributes)
+    multiplier_id = get_time_series_multiplier_id(attributes)
+    subsystem = get_subsystem(attributes)
+    template = get_template(model)
+    if isempty(subsystem)
+        device_model = get_model(template, V)
+    else
+        device_model = get_model(template, V, subsystem)
+    end
+    components = get_available_components(device_model, get_system(model))
+    ts_uuids = Set{String}()
+    for component in components
+        ts_uuid = string(IS.get_time_series_uuid(U, component, ts_name))
+        if !(ts_uuid in ts_uuids)
+            ts_vector = get_time_series_values!(
+                U,
+                model,
+                component,
+                ts_name,
+                multiplier_id,
+                initial_forecast_time,
+                horizon,
+            )
+            for (t, value) in enumerate(ts_vector)
+                if !isfinite(value)
+                    error("The value for the time series $(ts_name) is not finite. \
+                          Check that the data in the time series is valid.")
+                end
+                _set_param_value!(parameter_array, value, ts_uuid, t)
+            end
+            push!(ts_uuids, ts_uuid)
+        end
+    end
+end
+
+function _update_parameter_values!(
+    parameter_array::AbstractArray{T},
+    attributes::TimeSeriesAttributes{U},
+    service::V,
+    model::DecisionModel,
+    ::DatasetContainer{InMemoryDataset},
+) where {
+    T <: Union{JuMP.VariableRef, Float64},
+    U <: PSY.AbstractDeterministic,
+    V <: PSY.Service,
+}
+    initial_forecast_time = get_current_time(model) # Function not well defined for DecisionModels
+    horizon = get_time_steps(get_optimization_container(model))[end]
+    ts_name = get_time_series_name(attributes)
+    ts_uuid = string(IS.get_time_series_uuid(U, service, ts_name))
+    ts_vector = get_time_series_values!(
+        U,
+        model,
+        service,
+        get_time_series_name(attributes),
+        get_time_series_multiplier_id(attributes),
+        initial_forecast_time,
+        horizon,
+    )
+    for (t, value) in enumerate(ts_vector)
+        if !isfinite(value)
+            error("The value for the time series $(ts_name) is not finite. \
+                  Check that the data in the time series is valid.")
+        end
+        _set_param_value!(parameter_array, value, ts_uuid, t)
+    end
+end
+
+function _update_parameter_values!(
+    parameter_array::AbstractArray{T},
+    attributes::TimeSeriesAttributes{U},
+    ::Type{V},
+    model::EmulationModel,
+    ::DatasetContainer{InMemoryDataset},
+) where {T <: Union{JuMP.VariableRef, Float64}, U <: PSY.SingleTimeSeries, V <: PSY.Device}
+    initial_forecast_time = get_current_time(model)
+    template = get_template(model)
+    device_model = get_model(template, V)
+    components = get_available_components(device_model, get_system(model))
+    ts_name = get_time_series_name(attributes)
+    ts_uuids = Set{String}()
+    for component in components
+        ts_uuid = string(IS.get_time_series_uuid(U, component, ts_name))
+        if !(ts_uuid in ts_uuids)
+            # Note: This interface reads one single value per component at a time.
+            value = get_time_series_values!(
+                U,
+                model,
+                component,
+                get_time_series_name(attributes),
+                get_time_series_multiplier_id(attributes),
+                initial_forecast_time,
+            )[1]
+            if !isfinite(value)
+                error("The value for the time series $(ts_name) is not finite. \
+                      Check that the data in the time series is valid.")
+            end
+            _set_param_value!(parameter_array, value, ts_uuid, 1)
+            push!(ts_uuids, ts_uuid)
+        end
+    end
+    return
+end
+
+function _update_parameter_values!(
+    parameter_array::AbstractArray{T},
+    attributes::VariableValueAttributes,
+    ::Type{<:PSY.Device},
+    model::DecisionModel,
+    state::DatasetContainer{InMemoryDataset},
+) where {T <: Union{JuMP.VariableRef, Float64}}
+    current_time = get_current_time(model)
+    state_values = get_dataset_values(state, get_attribute_key(attributes))
+    component_names, time = axes(parameter_array)
+    model_resolution = get_resolution(model)
+    state_data = get_dataset(state, get_attribute_key(attributes))
+    state_timestamps = state_data.timestamps
+    max_state_index = get_num_rows(state_data)
+    if model_resolution < state_data.resolution
+        t_step = 1
+    else
+        t_step = model_resolution ÷ state_data.resolution
+    end
+    state_data_index = find_timestamp_index(state_timestamps, current_time)
+    sim_timestamps = range(current_time; step = model_resolution, length = time[end])
+    for t in time
+        timestamp_ix = min(max_state_index, state_data_index + t_step)
+        @debug "parameter horizon is over the step" max_state_index > state_data_index + 1
+        if state_timestamps[timestamp_ix] <= sim_timestamps[t]
+            state_data_index = timestamp_ix
+        end
+        for name in component_names
+            # Pass indices in this way since JuMP DenseAxisArray don't support view()
+            state_value = state_values[name, state_data_index]
+            if !isfinite(state_value)
+                error(
+                    "The value for the system state used in $(encode_key_as_string(get_attribute_key(attributes))) is not a finite value $(state_value) \
+                     This is commonly caused by referencing a state value at a time when such decision hasn't been made. \
+                     Consider reviewing your models' horizon and interval definitions",
+                )
+            end
+            _set_param_value!(parameter_array, state_value, name, t)
+        end
+    end
+    return
+end
+
+function _update_parameter_values!(
+    parameter_array::AbstractArray{T},
+    attributes::VariableValueAttributes,
+    ::PSY.Reserve,
+    model::DecisionModel,
+    state::DatasetContainer{InMemoryDataset},
+) where {T <: Union{JuMP.VariableRef, Float64}}
+    current_time = get_current_time(model)
+    state_values = get_dataset_values(state, get_attribute_key(attributes))
+    component_names, time = axes(parameter_array)
+    model_resolution = get_resolution(model)
+    state_data = get_dataset(state, get_attribute_key(attributes))
+    state_timestamps = state_data.timestamps
+    max_state_index = get_num_rows(state_data)
+    if model_resolution < state_data.resolution
+        t_step = 1
+    else
+        t_step = model_resolution ÷ state_data.resolution
+    end
+    state_data_index = find_timestamp_index(state_timestamps, current_time)
+    sim_timestamps = range(current_time; step = model_resolution, length = time[end])
+    for t in time
+        timestamp_ix = min(max_state_index, state_data_index + t_step)
+        @debug "parameter horizon is over the step" max_state_index > state_data_index + 1
+        if state_timestamps[timestamp_ix] <= sim_timestamps[t]
+            state_data_index = timestamp_ix
+        end
+        for name in component_names
+            # Pass indices in this way since JuMP DenseAxisArray don't support view()
+            state_value = state_values[name, state_data_index]
+            if !isfinite(state_value)
+                error(
+                    "The value for the system state used in $(encode_key_as_string(get_attribute_key(attributes))) is not a finite value $(state_value) \
+                     This is commonly caused by referencing a state value at a time when such decision hasn't been made. \
+                     Consider reviewing your models' horizon and interval definitions",
+                )
+            end
+            _set_param_value!(parameter_array, state_value, name, t)
+        end
+    end
+    return
+end
+
+function _update_parameter_values!(
+    parameter_array::AbstractArray{T},
+    attributes::VariableValueAttributes{VariableKey{OnVariable, U}},
+    ::Type{U},
+    model::DecisionModel,
+    state::DatasetContainer{InMemoryDataset},
+) where {T <: Union{JuMP.VariableRef, Float64}, U <: PSY.Device}
+    current_time = get_current_time(model)
+    state_values = get_dataset_values(state, get_attribute_key(attributes))
+    component_names, time = axes(parameter_array)
+    model_resolution = get_resolution(model)
+    state_data = get_dataset(state, get_attribute_key(attributes))
+    state_timestamps = state_data.timestamps
+    max_state_index = get_num_rows(state_data)
+    if model_resolution < state_data.resolution
+        t_step = 1
+    else
+        t_step = model_resolution ÷ state_data.resolution
+    end
+    state_data_index = find_timestamp_index(state_timestamps, current_time)
+
+    sim_timestamps = range(current_time; step = model_resolution, length = time[end])
+    for t in time
+        timestamp_ix = min(max_state_index, state_data_index + t_step)
+        @debug "parameter horizon is over the step" max_state_index > state_data_index + 1
+        if state_timestamps[timestamp_ix] <= sim_timestamps[t]
+            state_data_index = timestamp_ix
+        end
+        for name in component_names
+            # Pass indices in this way since JuMP DenseAxisArray don't support view()
+            value = round(state_values[name, state_data_index])
+            if !isfinite(value)
+                error(
+                    "The value for the system state used in $(encode_key_as_string(get_attribute_key(attributes))) is not a finite value $(value) \
+                     This is commonly caused by referencing a state value at a time when such decision hasn't been made. \
+                     Consider reviewing your models' horizon and interval definitions",
+                )
+            end
+            if 0.0 > value || value > 1.0
+                error(
+                    "The value for the system state used in $(encode_key_as_string(get_attribute_key(attributes))): $(value) is out of the [0, 1] range",
+                )
+            end
+            _set_param_value!(parameter_array, value, name, t)
+        end
+    end
+    return
+end
+
+function _update_parameter_values!(
+    parameter_array::AbstractArray{T},
+    attributes::VariableValueAttributes,
+    ::Type{<:PSY.Component},
+    model::EmulationModel,
+    state::DatasetContainer{InMemoryDataset},
+) where {T <: Union{JuMP.VariableRef, Float64}}
+    current_time = get_current_time(model)
+    state_values = get_dataset_values(state, get_attribute_key(attributes))
+    component_names, _ = axes(parameter_array)
+    state_data = get_dataset(state, get_attribute_key(attributes))
+    state_timestamps = state_data.timestamps
+    state_data_index = find_timestamp_index(state_timestamps, current_time)
+    for name in component_names
+        # Pass indices in this way since JuMP DenseAxisArray don't support view()
+        _set_param_value!(parameter_array, state_values[name, state_data_index], name, 1)
+    end
+    return
+end
+
+function _update_parameter_values!(
+    parameter_array::AbstractArray{T},
+    attributes::VariableValueAttributes{VariableKey{OnVariable, U}},
+    ::Type{<:PSY.Component},
+    model::EmulationModel,
+    state::DatasetContainer{InMemoryDataset},
+) where {T <: Union{JuMP.VariableRef, Float64}, U <: PSY.Component}
+    current_time = get_current_time(model)
+    state_values = get_dataset_values(state, get_attribute_key(attributes))
+    component_names, _ = axes(parameter_array)
+    state_data = get_dataset(state, get_attribute_key(attributes))
+    state_timestamps = state_data.timestamps
+    state_data_index = find_timestamp_index(state_timestamps, current_time)
+    for name in component_names
+        # Pass indices in this way since JuMP DenseAxisArray don't support view()
+        value = round(state_values[name, state_data_index])
+        if !isfinite(value)
+            error(
+                "The value for the system state used in $(encode_key_as_string(get_attribute_key(attributes))) is not a finite value $(value) \
+                 This is commonly caused by referencing a state value at a time when such decision hasn't been made. \
+                 Consider reviewing your models' horizon and interval definitions",
+            )
+        end
+        if 0.0 > value || value > 1.0
+            error(
+                "The value for the system state used in $(encode_key_as_string(get_attribute_key(attributes))): $(value) is out of the [0, 1] range",
+            )
+        end
+        _set_param_value!(parameter_array, value, name, 1)
+    end
+    return
+end
+
+function _update_parameter_values!(
+    ::AbstractArray{T},
+    ::VariableValueAttributes,
+    ::Type{<:PSY.Component},
+    ::EmulationModel,
+    ::EmulationModelStore,
+) where {T <: Union{JuMP.VariableRef, Float64}}
+    error("The emulation model has parameters that can't be updated from its results")
+    return
+end
+
+"""
+Update parameter function an OperationModel
+"""
+function update_container_parameter_values!(
+    optimization_container::OptimizationContainer,
+    model::OperationModel,
+    key::ParameterKey{T, U},
+    input::DatasetContainer{InMemoryDataset},
+) where {T <: ParameterType, U <: PSY.Component}
+    # Enable again for detailed debugging
+    # TimerOutputs.@timeit RUN_SIMULATION_TIMER "$T $U Parameter Update" begin
+    # Note: Do not instantite a new key here because it might not match the param keys in the container
+    # if the keys have strings in the meta fields
+    parameter_array = get_parameter_array(optimization_container, key)
+    parameter_attributes = get_parameter_attributes(optimization_container, key)
+    _update_parameter_values!(parameter_array, parameter_attributes, U, model, input)
+    return
+end
+
+function update_container_parameter_values!(
+    optimization_container::OptimizationContainer,
+    model::OperationModel,
+    key::ParameterKey{T, U},
+    input::DatasetContainer{InMemoryDataset},
+) where {T <: ObjectiveFunctionParameter, U <: PSY.Component}
+    # Note: Do not instantite a new key here because it might not match the param keys in the container
+    # if the keys have strings in the meta fields
+    parameter_array = get_parameter_array(optimization_container, key)
+    # Multiplier is only needed for the objective function since `_update_parameter_values!` also updates the objective function
+    parameter_multiplier = get_parameter_multiplier_array(optimization_container, key)
+    parameter_attributes = get_parameter_attributes(optimization_container, key)
+    _update_parameter_values!(
+        parameter_array,
+        parameter_multiplier,
+        parameter_attributes,
+        U,
+        model,
+        input,
+    )
+    return
+end
+
+function update_container_parameter_values!(
+    optimization_container::OptimizationContainer,
+    model::OperationModel,
+    key::ParameterKey{T, U},
+    input::DatasetContainer{InMemoryDataset},
+) where {T <: ObjectiveFunctionParameter, U <: PSY.Service}
+    # Note: Do not instantite a new key here because it might not match the param keys in the container
+    # if the keys have strings in the meta fields
+    parameter_array = get_parameter_array(optimization_container, key)
+    # Multiplier is only needed for the objective function since `_update_parameter_values!` also updates the objective function
+    parameter_multiplier = get_parameter_multiplier_array(optimization_container, key)
+    parameter_attributes = get_parameter_attributes(optimization_container, key)
+    _update_parameter_values!(
+        parameter_array,
+        parameter_multiplier,
+        parameter_attributes,
+        U,
+        model,
+        input,
+    )
+    return
+end
+
+function update_container_parameter_values!(
+    optimization_container::OptimizationContainer,
+    model::OperationModel,
+    key::ParameterKey{FixValueParameter, U},
+    input::DatasetContainer{InMemoryDataset},
+) where {U <: PSY.Component}
+    # Note: Do not instantite a new key here because it might not match the param keys in the container
+    # if the keys have strings in the meta fields
+    parameter_array = get_parameter_array(optimization_container, key)
+    parameter_attributes = get_parameter_attributes(optimization_container, key)
+    _update_parameter_values!(parameter_array, parameter_attributes, U, model, input)
+    _fix_parameter_value!(optimization_container, parameter_array, parameter_attributes)
+    return
+end
+
+function update_container_parameter_values!(
+    optimization_container::OptimizationContainer,
+    model::OperationModel,
+    key::ParameterKey{FixValueParameter, U},
+    input::DatasetContainer{InMemoryDataset},
+) where {U <: PSY.Service}
+    # Note: Do not instantite a new key here because it might not match the param keys in the container
+    # if the keys have strings in the meta fields
+    parameter_array = get_parameter_array(optimization_container, key)
+    parameter_attributes = get_parameter_attributes(optimization_container, key)
+    service = PSY.get_component(U, get_system(model), key.meta)
+    @assert service !== nothing
+    _update_parameter_values!(parameter_array, parameter_attributes, U, model, input)
+    _fix_parameter_value!(optimization_container, parameter_array, parameter_attributes)
+    return
+end
+
+function update_container_parameter_values!(
+    optimization_container::OptimizationContainer,
+    model::OperationModel,
+    key::ParameterKey{T, U},
+    input::DatasetContainer{InMemoryDataset},
+) where {T <: ParameterType, U <: PSY.Service}
+    # Note: Do not instantite a new key here because it might not match the param keys in the container
+    # if the keys have strings in the meta fields
+    parameter_array = get_parameter_array(optimization_container, key)
+    parameter_attributes = get_parameter_attributes(optimization_container, key)
+    service = PSY.get_component(U, get_system(model), key.meta)
+    @assert service !== nothing
+    _update_parameter_values!(parameter_array, parameter_attributes, service, model, input)
+    return
+end
diff --git a/src/parameters/update_cost_parameters.jl b/src/parameters/update_cost_parameters.jl
new file mode 100644
index 0000000000..600eca675c
--- /dev/null
+++ b/src/parameters/update_cost_parameters.jl
@@ -0,0 +1,125 @@
+function _update_parameter_values!(
+    parameter_array::DenseAxisArray,
+    parameter_multiplier::JuMPFloatArray,
+    attributes::CostFunctionAttributes,
+    ::Type{V},
+    model::DecisionModel,
+    ::DatasetContainer{InMemoryDataset},
+) where {V <: PSY.Component}
+    initial_forecast_time = get_current_time(model) # Function not well defined for DecisionModels
+    time_steps = get_time_steps(get_optimization_container(model))
+    horizon = time_steps[end]
+    container = get_optimization_container(model)
+    @assert !is_synchronized(container)
+    template = get_template(model)
+    device_model = get_model(template, V)
+    components = get_available_components(device_model, get_system(model))
+
+    for component in components
+        if _has_variable_cost_parameter(component)
+            name = PSY.get_name(component)
+            ts_vector = PSY.get_variable_cost(
+                component,
+                PSY.get_operation_cost(component);
+                start_time = initial_forecast_time,
+                len = horizon,
+            )
+            variable_cost_forecast_values = TimeSeries.values(ts_vector)
+            for (t, value) in enumerate(variable_cost_forecast_values)
+                if attributes.uses_compact_power
+                    # TODO implement this
+                    value, _ = _convert_variable_cost(value)
+                end
+                # TODO removed an apparently unused block of code here?
+                _set_param_value!(parameter_array, value, name, t)
+                update_variable_cost!(
+                    container,
+                    parameter_array,
+                    parameter_multiplier,
+                    attributes,
+                    component,
+                    t,
+                )
+            end
+        end
+    end
+    return
+end
+
+_has_variable_cost_parameter(component::PSY.Component) =
+    _has_variable_cost_parameter(PSY.get_operation_cost(component))
+_has_variable_cost_parameter(::PSY.MarketBidCost) = true
+_has_variable_cost_parameter(::T) where {T <: PSY.OperationalCost} = false
+
+function _update_pwl_cost_expression(
+    container::OptimizationContainer,
+    ::Type{T},
+    component_name::String,
+    time_period::Int,
+    cost_data::PSY.PiecewiseLinearData,
+) where {T <: PSY.Component}
+    pwl_var_container = get_variable(container, PieceWiseLinearCostVariable(), T)
+    resolution = get_resolution(container)
+    dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR
+    gen_cost = JuMP.AffExpr(0.0)
+    slopes = PSY.get_slopes(cost_data)
+    upb = get_breakpoint_upper_bounds(cost_data)
+    for i in 1:length(cost_data)
+        JuMP.add_to_expression!(
+            gen_cost,
+            slopes[i] * upb[i] * dt * pwl_var_container[(component_name, i, time_period)],
+        )
+    end
+    return gen_cost
+end
+
+function update_variable_cost!(
+    container::OptimizationContainer,
+    parameter_array::JuMPFloatArray,
+    parameter_multiplier::JuMPFloatArray,
+    attributes::CostFunctionAttributes{Float64},
+    component::T,
+    time_period::Int,
+) where {T <: PSY.Component}
+    resolution = get_resolution(container)
+    dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR
+    base_power = get_base_power(container)
+    component_name = PSY.get_name(component)
+    cost_data = parameter_array[component_name, time_period]  # TODO is this a new-style cost?
+    if iszero(cost_data)
+        return
+    end
+    mult_ = parameter_multiplier[component_name, time_period]
+    variable = get_variable(container, get_variable_type(attributes)(), T)
+    gen_cost = variable[component_name, time_period] * mult_ * cost_data * base_power * dt
+    add_to_objective_variant_expression!(container, gen_cost)
+    set_expression!(container, ProductionCostExpression, gen_cost, component, time_period)
+    return
+end
+
+function update_variable_cost!(
+    container::OptimizationContainer,
+    parameter_array::DenseAxisArray{Vector{NTuple{2, Float64}}},
+    parameter_multiplier::JuMPFloatArray,
+    ::CostFunctionAttributes{Vector{NTuple{2, Float64}}},
+    component::T,
+    time_period::Int,
+) where {T <: PSY.Component}
+    component_name = PSY.get_name(component)
+    cost_data = parameter_array[component_name, time_period]
+    if all(iszero.(last.(cost_data)))
+        return
+    end
+    mult_ = parameter_multiplier[component_name, time_period]
+    gen_cost =
+        _update_pwl_cost_expression(
+            container,
+            T,
+            component_name,
+            time_period,
+            PSY.PiecewiseLinearData(cost_data),
+        )
+    add_to_objective_variant_expression!(container, mult_ * gen_cost)
+    set_expression!(container, ProductionCostExpression, gen_cost, component, time_period)
+    return
+end
diff --git a/src/parameters/update_parameters.jl b/src/parameters/update_parameters.jl
index 39e4b8f996..53e3f4d35b 100644
--- a/src/parameters/update_parameters.jl
+++ b/src/parameters/update_parameters.jl
@@ -1,455 +1,15 @@
-function _update_parameter_values!(
-    ::AbstractArray{T},
-    ::NoAttributes,
-    args...,
-) where {T <: Union{Float64, JuMP.VariableRef}} end
-
-######################## Methods to update Parameters from Time Series #####################
-function _set_param_value!(param::JuMPVariableMatrix, value::Float64, name::String, t::Int)
-    fix_parameter_value(param[name, t], value)
-    return
-end
-
-function _set_param_value!(
-    param::DenseAxisArray{Vector{NTuple{2, Float64}}},
-    value::Vector{NTuple{2, Float64}},
-    name::String,
-    t::Int,
-)
-    param[name, t] = value
-    return
-end
-
-function _set_param_value!(param::JuMPFloatArray, value::Float64, name::String, t::Int)
-    param[name, t] = value
-    return
-end
-
-function _update_parameter_values!(
-    parameter_array::AbstractArray{T},
-    attributes::TimeSeriesAttributes{U},
-    ::Type{V},
-    model::DecisionModel,
-    ::DatasetContainer{InMemoryDataset},
-) where {
-    T <: Union{JuMP.VariableRef, Float64},
-    U <: PSY.AbstractDeterministic,
-    V <: PSY.Component,
-}
-    initial_forecast_time = get_current_time(model) # Function not well defined for DecisionModels
-    horizon = get_time_steps(get_optimization_container(model))[end]
-    ts_name = get_time_series_name(attributes)
-    multiplier_id = get_time_series_multiplier_id(attributes)
-    template = get_template(model)
-    device_model = get_model(template, V)
-    components = get_available_components(device_model, get_system(model))
-    ts_uuids = Set{String}()
-    for component in components
-        ts_uuid = string(IS.get_time_series_uuid(U, component, ts_name))
-        if !(ts_uuid in ts_uuids)
-            ts_vector = get_time_series_values!(
-                U,
-                model,
-                component,
-                ts_name,
-                multiplier_id,
-                initial_forecast_time,
-                horizon,
-            )
-            for (t, value) in enumerate(ts_vector)
-                if !isfinite(value)
-                    error("The value for the time series $(ts_name) is not finite. \
-                          Check that the data in the time series is valid.")
-                end
-                _set_param_value!(parameter_array, value, ts_uuid, t)
-            end
-            push!(ts_uuids, ts_uuid)
-        end
-    end
-end
-
-function _update_parameter_values!(
-    parameter_array::AbstractArray{T},
-    attributes::TimeSeriesAttributes{U},
-    service::V,
-    model::DecisionModel,
-    ::DatasetContainer{InMemoryDataset},
-) where {
-    T <: Union{JuMP.VariableRef, Float64},
-    U <: PSY.AbstractDeterministic,
-    V <: PSY.Service,
-}
-    initial_forecast_time = get_current_time(model) # Function not well defined for DecisionModels
-    horizon = get_time_steps(get_optimization_container(model))[end]
-    ts_name = get_time_series_name(attributes)
-    ts_uuid = string(IS.get_time_series_uuid(U, service, ts_name))
-    ts_vector = get_time_series_values!(
-        U,
-        model,
-        service,
-        get_time_series_name(attributes),
-        get_time_series_multiplier_id(attributes),
-        initial_forecast_time,
-        horizon,
-    )
-    for (t, value) in enumerate(ts_vector)
-        if !isfinite(value)
-            error("The value for the time series $(ts_name) is not finite. \
-                  Check that the data in the time series is valid.")
-        end
-        _set_param_value!(parameter_array, value, ts_uuid, t)
-    end
-end
-
-function _update_parameter_values!(
-    parameter_array::AbstractArray{T},
-    attributes::TimeSeriesAttributes{U},
-    ::Type{V},
-    model::EmulationModel,
-    ::DatasetContainer{InMemoryDataset},
-) where {T <: Union{JuMP.VariableRef, Float64}, U <: PSY.SingleTimeSeries, V <: PSY.Device}
-    initial_forecast_time = get_current_time(model)
-    template = get_template(model)
-    device_model = get_model(template, V)
-    components = get_available_components(device_model, get_system(model))
-    ts_name = get_time_series_name(attributes)
-    ts_uuids = Set{String}()
-    for component in components
-        ts_uuid = string(IS.get_time_series_uuid(U, component, ts_name))
-        if !(ts_uuid in ts_uuids)
-            # Note: This interface reads one single value per component at a time.
-            value = get_time_series_values!(
-                U,
-                model,
-                component,
-                get_time_series_name(attributes),
-                get_time_series_multiplier_id(attributes),
-                initial_forecast_time,
-            )[1]
-            if !isfinite(value)
-                error("The value for the time series $(ts_name) is not finite. \
-                      Check that the data in the time series is valid.")
-            end
-            _set_param_value!(parameter_array, value, ts_uuid, 1)
-            push!(ts_uuids, ts_uuid)
-        end
-    end
-    return
-end
-
-function _update_parameter_values!(
-    parameter_array::AbstractArray{T},
-    attributes::VariableValueAttributes,
-    ::Type{<:PSY.Device},
-    model::DecisionModel,
-    state::DatasetContainer{InMemoryDataset},
-) where {T <: Union{JuMP.VariableRef, Float64}}
-    current_time = get_current_time(model)
-    state_values = get_dataset_values(state, get_attribute_key(attributes))
-    component_names, time = axes(parameter_array)
-    model_resolution = get_resolution(model)
-    state_data = get_dataset(state, get_attribute_key(attributes))
-    state_timestamps = state_data.timestamps
-    max_state_index = get_num_rows(state_data)
-    if model_resolution < state_data.resolution
-        t_step = 1
-    else
-        t_step = model_resolution ÷ state_data.resolution
-    end
-    state_data_index = find_timestamp_index(state_timestamps, current_time)
-    sim_timestamps = range(current_time; step = model_resolution, length = time[end])
-    for t in time
-        timestamp_ix = min(max_state_index, state_data_index + t_step)
-        @debug "parameter horizon is over the step" max_state_index > state_data_index + 1
-        if state_timestamps[timestamp_ix] <= sim_timestamps[t]
-            state_data_index = timestamp_ix
-        end
-        for name in component_names
-            # Pass indices in this way since JuMP DenseAxisArray don't support view()
-            state_value = state_values[name, state_data_index]
-            if !isfinite(state_value)
-                error(
-                    "The value for the system state used in $(encode_key_as_string(get_attribute_key(attributes))) is not a finite value $(state_value) \
-                     This is commonly caused by referencing a state value at a time when such decision hasn't been made. \
-                     Consider reviewing your models' horizon and interval definitions",
-                )
-            end
-            _set_param_value!(parameter_array, state_value, name, t)
-        end
-    end
-    return
-end
-
-function _update_parameter_values!(
-    parameter_array::AbstractArray{T},
-    attributes::VariableValueAttributes,
-    ::PSY.Reserve,
-    model::DecisionModel,
-    state::DatasetContainer{InMemoryDataset},
-) where {T <: Union{JuMP.VariableRef, Float64}}
-    current_time = get_current_time(model)
-    state_values = get_dataset_values(state, get_attribute_key(attributes))
-    component_names, time = axes(parameter_array)
-    model_resolution = get_resolution(model)
-    state_data = get_dataset(state, get_attribute_key(attributes))
-    state_timestamps = state_data.timestamps
-    max_state_index = get_num_rows(state_data)
-    if model_resolution < state_data.resolution
-        t_step = 1
-    else
-        t_step = model_resolution ÷ state_data.resolution
-    end
-    state_data_index = find_timestamp_index(state_timestamps, current_time)
-    sim_timestamps = range(current_time; step = model_resolution, length = time[end])
-    for t in time
-        timestamp_ix = min(max_state_index, state_data_index + t_step)
-        @debug "parameter horizon is over the step" max_state_index > state_data_index + 1
-        if state_timestamps[timestamp_ix] <= sim_timestamps[t]
-            state_data_index = timestamp_ix
-        end
-        for name in component_names
-            # Pass indices in this way since JuMP DenseAxisArray don't support view()
-            state_value = state_values[name, state_data_index]
-            if !isfinite(state_value)
-                error(
-                    "The value for the system state used in $(encode_key_as_string(get_attribute_key(attributes))) is not a finite value $(state_value) \
-                     This is commonly caused by referencing a state value at a time when such decision hasn't been made. \
-                     Consider reviewing your models' horizon and interval definitions",
-                )
-            end
-            _set_param_value!(parameter_array, state_value, name, t)
-        end
-    end
-    return
-end
-
-function _update_parameter_values!(
-    parameter_array::AbstractArray{T},
-    attributes::VariableValueAttributes{VariableKey{OnVariable, U}},
-    ::Type{U},
-    model::DecisionModel,
-    state::DatasetContainer{InMemoryDataset},
-) where {T <: Union{JuMP.VariableRef, Float64}, U <: PSY.Device}
-    current_time = get_current_time(model)
-    state_values = get_dataset_values(state, get_attribute_key(attributes))
-    component_names, time = axes(parameter_array)
-    model_resolution = get_resolution(model)
-    state_data = get_dataset(state, get_attribute_key(attributes))
-    state_timestamps = state_data.timestamps
-    max_state_index = get_num_rows(state_data)
-    if model_resolution < state_data.resolution
-        t_step = 1
-    else
-        t_step = model_resolution ÷ state_data.resolution
-    end
-    state_data_index = find_timestamp_index(state_timestamps, current_time)
-
-    sim_timestamps = range(current_time; step = model_resolution, length = time[end])
-    for t in time
-        timestamp_ix = min(max_state_index, state_data_index + t_step)
-        @debug "parameter horizon is over the step" max_state_index > state_data_index + 1
-        if state_timestamps[timestamp_ix] <= sim_timestamps[t]
-            state_data_index = timestamp_ix
-        end
-        for name in component_names
-            # Pass indices in this way since JuMP DenseAxisArray don't support view()
-            value = round(state_values[name, state_data_index])
-            @assert 0.0 <= value <= 1.0
-            if !isfinite(value)
-                error(
-                    "The value for the system state used in $(encode_key_as_string(get_attribute_key(attributes))) is not a finite value $(value) \
-                     This is commonly caused by referencing a state value at a time when such decision hasn't been made. \
-                     Consider reviewing your models' horizon and interval definitions",
-                )
-            end
-            _set_param_value!(parameter_array, value, name, t)
-        end
-    end
-    return
-end
-
-function _update_parameter_values!(
-    parameter_array::AbstractArray{T},
-    attributes::VariableValueAttributes,
-    ::Type{<:PSY.Component},
-    model::EmulationModel,
-    state::DatasetContainer{InMemoryDataset},
-) where {T <: Union{JuMP.VariableRef, Float64}}
-    current_time = get_current_time(model)
-    state_values = get_dataset_values(state, get_attribute_key(attributes))
-    component_names, _ = axes(parameter_array)
-    state_data = get_dataset(state, get_attribute_key(attributes))
-    state_timestamps = state_data.timestamps
-    state_data_index = find_timestamp_index(state_timestamps, current_time)
-    for name in component_names
-        # Pass indices in this way since JuMP DenseAxisArray don't support view()
-        _set_param_value!(parameter_array, state_values[name, state_data_index], name, 1)
-    end
-    return
-end
-
-function _update_parameter_values!(
-    parameter_array::AbstractArray{T},
-    attributes::VariableValueAttributes{VariableKey{OnVariable, U}},
-    ::Type{<:PSY.Component},
-    model::EmulationModel,
-    state::DatasetContainer{InMemoryDataset},
-) where {T <: Union{JuMP.VariableRef, Float64}, U <: PSY.Component}
-    current_time = get_current_time(model)
-    state_values = get_dataset_values(state, get_attribute_key(attributes))
-    component_names, _ = axes(parameter_array)
-    state_data = get_dataset(state, get_attribute_key(attributes))
-    state_timestamps = state_data.timestamps
-    state_data_index = find_timestamp_index(state_timestamps, current_time)
-    for name in component_names
-        # Pass indices in this way since JuMP DenseAxisArray don't support view()
-        value = round(state_values[name, state_data_index])
-        @assert 0.0 <= value <= 1.0
-        if !isfinite(value)
-            error(
-                "The value for the system state used in $(encode_key_as_string(get_attribute_key(attributes))) is not a finite value $(value) \
-                 This is commonly caused by referencing a state value at a time when such decision hasn't been made. \
-                 Consider reviewing your models' horizon and interval definitions",
-            )
-        end
-        _set_param_value!(parameter_array, value, name, 1)
-    end
-    return
-end
-
-function _update_parameter_values!(
-    ::AbstractArray{T},
-    ::VariableValueAttributes,
-    ::Type{<:PSY.Component},
-    ::EmulationModel,
-    ::EmulationModelStore,
-) where {T <: Union{JuMP.VariableRef, Float64}}
-    error("The emulation model has parameters that can't be updated from its results")
-    return
-end
-
-"""
-Update parameter function an OperationModel
-"""
-function update_container_parameter_values!(
-    optimization_container::OptimizationContainer,
-    model::OperationModel,
-    key::ParameterKey{T, U},
-    input::DatasetContainer{InMemoryDataset},
-) where {T <: ParameterType, U <: PSY.Component}
-    # Enable again for detailed debugging
-    # TimerOutputs.@timeit RUN_SIMULATION_TIMER "$T $U Parameter Update" begin
-    # Note: Do not instantite a new key here because it might not match the param keys in the container
-    # if the keys have strings in the meta fields
-    parameter_array = get_parameter_array(optimization_container, key)
-    parameter_attributes = get_parameter_attributes(optimization_container, key)
-    _update_parameter_values!(parameter_array, parameter_attributes, U, model, input)
-    return
-end
-
-function update_container_parameter_values!(
-    optimization_container::OptimizationContainer,
-    model::OperationModel,
-    key::ParameterKey{T, U},
-    input::DatasetContainer{InMemoryDataset},
-) where {T <: ObjectiveFunctionParameter, U <: PSY.Component}
-    # Note: Do not instantite a new key here because it might not match the param keys in the container
-    # if the keys have strings in the meta fields
-    parameter_array = get_parameter_array(optimization_container, key)
-    # Multiplier is only needed for the objective function since `_update_parameter_values!` also updates the objective function
-    parameter_multiplier = get_parameter_multiplier_array(optimization_container, key)
-    parameter_attributes = get_parameter_attributes(optimization_container, key)
-    _update_parameter_values!(
-        parameter_array,
-        parameter_multiplier,
-        parameter_attributes,
-        U,
-        model,
-        input,
-    )
-    return
-end
-
-function update_container_parameter_values!(
-    optimization_container::OptimizationContainer,
-    model::OperationModel,
-    key::ParameterKey{T, U},
-    input::DatasetContainer{InMemoryDataset},
-) where {T <: ObjectiveFunctionParameter, U <: PSY.Service}
-    # Note: Do not instantite a new key here because it might not match the param keys in the container
-    # if the keys have strings in the meta fields
-    parameter_array = get_parameter_array(optimization_container, key)
-    # Multiplier is only needed for the objective function since `_update_parameter_values!` also updates the objective function
-    parameter_multiplier = get_parameter_multiplier_array(optimization_container, key)
-    parameter_attributes = get_parameter_attributes(optimization_container, key)
-    _update_parameter_values!(
-        parameter_array,
-        parameter_multiplier,
-        parameter_attributes,
-        U,
-        model,
-        input,
-    )
-    return
-end
-
-function update_container_parameter_values!(
-    optimization_container::OptimizationContainer,
-    model::OperationModel,
-    key::ParameterKey{FixValueParameter, U},
-    input::DatasetContainer{InMemoryDataset},
-) where {U <: PSY.Component}
-    # Note: Do not instantite a new key here because it might not match the param keys in the container
-    # if the keys have strings in the meta fields
-    parameter_array = get_parameter_array(optimization_container, key)
-    parameter_attributes = get_parameter_attributes(optimization_container, key)
-    _update_parameter_values!(parameter_array, parameter_attributes, T, model, input)
-    _fix_parameter_value!(optimization_container, parameter_array, parameter_attributes)
-    return
-end
-
-function update_container_parameter_values!(
-    optimization_container::OptimizationContainer,
-    model::OperationModel,
-    key::ParameterKey{FixValueParameter, U},
-    input::DatasetContainer{InMemoryDataset},
-) where {U <: PSY.Service}
-    # Note: Do not instantite a new key here because it might not match the param keys in the container
-    # if the keys have strings in the meta fields
-    parameter_array = get_parameter_array(optimization_container, key)
-    parameter_attributes = get_parameter_attributes(optimization_container, key)
-    _update_parameter_values!(parameter_array, parameter_attributes, T, model, input)
-    _fix_parameter_value!(optimization_container, parameter_array, parameter_attributes)
-    return
-end
-
-function update_container_parameter_values!(
-    optimization_container::OptimizationContainer,
-    model::OperationModel,
-    key::ParameterKey{T, U},
-    input::DatasetContainer{InMemoryDataset},
-) where {T <: ParameterType, U <: PSY.Service}
-    # Note: Do not instantite a new key here because it might not match the param keys in the container
-    # if the keys have strings in the meta fields
-    parameter_array = get_parameter_array(optimization_container, key)
-    parameter_attributes = get_parameter_attributes(optimization_container, key)
-    service = PSY.get_component(U, get_system(model), key.meta)
-    @assert service !== nothing
-    _update_parameter_values!(parameter_array, parameter_attributes, service, model, input)
-    return
-end
-
 """
 Update parameter function an OperationModel
 """
 function update_parameter_values!(
     model::OperationModel,
     key::ParameterKey{T, U},
-    input::DatasetContainer{InMemoryDataset},
+    simulation_state::SimulationState,
 ) where {T <: ParameterType, U <: PSY.Component}
     # Enable again for detailed debugging
     # TimerOutputs.@timeit RUN_SIMULATION_TIMER "$T $U Parameter Update" begin
     optimization_container = get_optimization_container(model)
+    input = get_decision_states(simulation_state)
     update_container_parameter_values!(optimization_container, model, key, input)
     parameter_attributes = get_parameter_attributes(optimization_container, key)
     IS.@record :execution ParameterUpdateEvent(
@@ -483,7 +43,7 @@ end
 function update_parameter_values!(
     model::OperationModel,
     key::ParameterKey{FixValueParameter, T},
-    input::DatasetContainer{InMemoryDataset},
+    simulation_state::SimulationState,
 ) where {T <: PSY.Service}
     # Enable again for detailed debugging
     # TimerOutputs.@timeit RUN_SIMULATION_TIMER "$T $U Parameter Update" begin
@@ -494,7 +54,14 @@ function update_parameter_values!(
     parameter_attributes = get_parameter_attributes(optimization_container, key)
     service = PSY.get_component(T, get_system(model), key.meta)
     @assert service !== nothing
-    _update_parameter_values!(parameter_array, parameter_attributes, service, model, input)
+    input = get_decision_states(simulation_state)
+    _update_parameter_values!(
+        parameter_array,
+        parameter_attributes,
+        service,
+        model,
+        input,
+    )
     _fix_parameter_value!(optimization_container, parameter_array, parameter_attributes)
     IS.@record :execution ParameterUpdateEvent(
         FixValueParameter,
@@ -506,129 +73,3 @@ function update_parameter_values!(
     #end
     return
 end
-
-function _update_parameter_values!(
-    parameter_array::DenseAxisArray,
-    parameter_multiplier::JuMPFloatArray,
-    attributes::CostFunctionAttributes,
-    ::Type{V},
-    model::DecisionModel,
-    ::DatasetContainer{InMemoryDataset},
-) where {V <: PSY.Component}
-    initial_forecast_time = get_current_time(model) # Function not well defined for DecisionModels
-    time_steps = get_time_steps(get_optimization_container(model))
-    horizon = time_steps[end]
-    container = get_optimization_container(model)
-    @assert !is_synchronized(container)
-    template = get_template(model)
-    device_model = get_model(template, V)
-    components = get_available_components(device_model, get_system(model))
-
-    for component in components
-        if _has_variable_cost_parameter(component)
-            name = PSY.get_name(component)
-            ts_vector = PSY.get_variable_cost(
-                component,
-                PSY.get_operation_cost(component);
-                start_time = initial_forecast_time,
-                len = horizon,
-            )
-            variable_cost_forecast_values = TimeSeries.values(ts_vector)
-            for (t, value) in enumerate(variable_cost_forecast_values)
-                if attributes.uses_compact_power
-                    # TODO implement this
-                    value, _ = _convert_variable_cost(value)
-                end
-                # TODO removed an apparently unused block of code here?
-                _set_param_value!(parameter_array, value, name, t)
-                update_variable_cost!(
-                    container,
-                    parameter_array,
-                    parameter_multiplier,
-                    attributes,
-                    component,
-                    t,
-                )
-            end
-        end
-    end
-    return
-end
-
-_has_variable_cost_parameter(component::PSY.Component) =
-    _has_variable_cost_parameter(PSY.get_operation_cost(component))
-_has_variable_cost_parameter(::PSY.MarketBidCost) = true
-_has_variable_cost_parameter(::T) where {T <: PSY.OperationalCost} = false
-
-function _update_pwl_cost_expression(
-    container::OptimizationContainer,
-    ::Type{T},
-    component_name::String,
-    time_period::Int,
-    cost_data::PSY.PiecewiseLinearData,
-) where {T <: PSY.Component}
-    pwl_var_container = get_variable(container, PieceWiseLinearCostVariable(), T)
-    resolution = get_resolution(container)
-    dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR
-    gen_cost = JuMP.AffExpr(0.0)
-    slopes = PSY.get_slopes(cost_data)
-    upb = get_breakpoint_upper_bounds(cost_data)
-    for i in 1:length(cost_data)
-        JuMP.add_to_expression!(
-            gen_cost,
-            slopes[i] * upb[i] * dt * pwl_var_container[(component_name, i, time_period)],
-        )
-    end
-    return gen_cost
-end
-
-function update_variable_cost!(
-    container::OptimizationContainer,
-    parameter_array::JuMPFloatArray,
-    parameter_multiplier::JuMPFloatArray,
-    attributes::CostFunctionAttributes{Float64},
-    component::T,
-    time_period::Int,
-) where {T <: PSY.Component}
-    resolution = get_resolution(container)
-    dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR
-    base_power = get_base_power(container)
-    component_name = PSY.get_name(component)
-    cost_data = parameter_array[component_name, time_period]  # TODO is this a new-style cost?
-    if iszero(cost_data)
-        return
-    end
-    mult_ = parameter_multiplier[component_name, time_period]
-    variable = get_variable(container, get_variable_type(attributes)(), T)
-    gen_cost = variable[component_name, time_period] * mult_ * cost_data * base_power * dt
-    add_to_objective_variant_expression!(container, gen_cost)
-    set_expression!(container, ProductionCostExpression, gen_cost, component, time_period)
-    return
-end
-
-function update_variable_cost!(
-    container::OptimizationContainer,
-    parameter_array::DenseAxisArray{Vector{NTuple{2, Float64}}},
-    parameter_multiplier::JuMPFloatArray,
-    ::CostFunctionAttributes{Vector{NTuple{2, Float64}}},
-    component::T,
-    time_period::Int,
-) where {T <: PSY.Component}
-    component_name = PSY.get_name(component)
-    cost_data = parameter_array[component_name, time_period]
-    if all(iszero.(last.(cost_data)))
-        return
-    end
-    mult_ = parameter_multiplier[component_name, time_period]
-    gen_cost =
-        _update_pwl_cost_expression(
-            container,
-            T,
-            component_name,
-            time_period,
-            PSY.PiecewiseLinearData(cost_data),
-        )
-    add_to_objective_variant_expression!(container, mult_ * gen_cost)
-    set_expression!(container, ProductionCostExpression, gen_cost, component, time_period)
-    return
-end
diff --git a/src/services_models/services_constructor.jl b/src/services_models/services_constructor.jl
index 09d9c9d2ec..6deb80d300 100644
--- a/src/services_models/services_constructor.jl
+++ b/src/services_models/services_constructor.jl
@@ -104,6 +104,7 @@ function construct_service!(
 ) where {SR <: PSY.Reserve}
     name = get_service_name(model)
     service = PSY.get_component(SR, sys, name)
+    !PSY.get_available(service) && return
     add_parameters!(container, RequirementTimeSeriesParameter, service, model)
     contributing_devices = get_contributing_devices(model)
 
@@ -130,6 +131,7 @@ function construct_service!(
 ) where {SR <: PSY.Reserve}
     name = get_service_name(model)
     service = PSY.get_component(SR, sys, name)
+    !PSY.get_available(service) && return
     contributing_devices = get_contributing_devices(model)
 
     add_constraints!(container, RequirementConstraint, service, contributing_devices, model)
@@ -160,6 +162,7 @@ function construct_service!(
 ) where {SR <: PSY.ConstantReserve}
     name = get_service_name(model)
     service = PSY.get_component(SR, sys, name)
+    !PSY.get_available(service) && return
     contributing_devices = get_contributing_devices(model)
 
     add_variables!(
@@ -185,6 +188,7 @@ function construct_service!(
 ) where {SR <: PSY.ConstantReserve}
     name = get_service_name(model)
     service = PSY.get_component(SR, sys, name)
+    !PSY.get_available(service) && return
     contributing_devices = get_contributing_devices(model)
 
     add_constraints!(container, RequirementConstraint, service, contributing_devices, model)
@@ -214,6 +218,7 @@ function construct_service!(
 ) where {SR <: PSY.Reserve}
     name = get_service_name(model)
     service = PSY.get_component(SR, sys, name)
+    !PSY.get_available(service) && return
     contributing_devices = get_contributing_devices(model)
     add_variable!(container, ServiceRequirementVariable(), service, StepwiseCostReserve())
     add_variables!(
@@ -239,6 +244,7 @@ function construct_service!(
 ) where {SR <: PSY.Reserve}
     name = get_service_name(model)
     service = PSY.get_component(SR, sys, name)
+    !PSY.get_available(service) && return
     contributing_devices = get_contributing_devices(model)
 
     add_constraints!(container, RequirementConstraint, service, contributing_devices, model)
@@ -251,6 +257,7 @@ function construct_service!(
     return
 end
 
+#=
 function construct_service!(
     container::OptimizationContainer,
     sys::PSY.System,
@@ -343,6 +350,7 @@ function construct_service!(
     objective_function!(container, services, model)
     return
 end
+=#
 
 """
     Constructs a service for ConstantReserveGroup.
@@ -358,6 +366,7 @@ function construct_service!(
 ) where {SR <: PSY.ConstantReserveGroup}
     name = get_service_name(model)
     service = PSY.get_component(SR, sys, name)
+    !PSY.get_available(service) && return
     contributing_services = PSY.get_contributing_services(service)
     # check if variables exist
     check_activeservice_variables(container, contributing_services)
@@ -376,6 +385,7 @@ function construct_service!(
 ) where {SR <: PSY.ConstantReserveGroup}
     name = get_service_name(model)
     service = PSY.get_component(SR, sys, name)
+    !PSY.get_available(service) && return
     contributing_services = PSY.get_contributing_services(service)
 
     add_constraints!(
@@ -401,6 +411,7 @@ function construct_service!(
 ) where {SR <: PSY.Reserve}
     name = get_service_name(model)
     service = PSY.get_component(SR, sys, name)
+    !PSY.get_available(service) && return
     contributing_devices = get_contributing_devices(model)
     add_parameters!(container, RequirementTimeSeriesParameter, service, model)
 
@@ -427,6 +438,7 @@ function construct_service!(
 ) where {SR <: PSY.Reserve}
     name = get_service_name(model)
     service = PSY.get_component(SR, sys, name)
+    !PSY.get_available(service) && return
     contributing_devices = get_contributing_devices(model)
 
     add_constraints!(container, RequirementConstraint, service, contributing_devices, model)
@@ -458,6 +470,7 @@ function construct_service!(
 ) where {SR <: PSY.ReserveNonSpinning}
     name = get_service_name(model)
     service = PSY.get_component(SR, sys, name)
+    !PSY.get_available(service) && return
     contributing_devices = get_contributing_devices(model)
     add_parameters!(container, RequirementTimeSeriesParameter, service, model)
 
@@ -483,6 +496,7 @@ function construct_service!(
 ) where {SR <: PSY.ReserveNonSpinning}
     name = get_service_name(model)
     service = PSY.get_component(SR, sys, name)
+    !PSY.get_available(service) && return
     contributing_devices = get_contributing_devices(model)
 
     add_constraints!(container, RequirementConstraint, service, contributing_devices, model)
@@ -549,6 +563,7 @@ function construct_service!(
 ) where {T <: PSY.TransmissionInterface}
     name = get_service_name(model)
     service = PSY.get_component(T, sys, name)
+    !PSY.get_available(service) && return
 
     add_to_expression!(
         container,
@@ -640,6 +655,7 @@ function construct_service!(
 ) where {T <: PSY.TransmissionInterface}
     name = get_service_name(model)
     service = PSY.get_component(T, sys, name)
+    !PSY.get_available(service) && return
 
     add_to_expression!(
         container,
diff --git a/src/simulation/initial_condition_update_simulation.jl b/src/simulation/initial_condition_update_simulation.jl
index 7d444cf609..4f4a46559e 100644
--- a/src/simulation/initial_condition_update_simulation.jl
+++ b/src/simulation/initial_condition_update_simulation.jl
@@ -89,7 +89,11 @@ function update_initial_conditions!(
                       $(comp_type)-$(comp_name) is out of bounds [$(min), $(max)].")
             end
         else
-            @assert isapprox(var_val, 0.0, atol = ABSOLUTE_TOLERANCE) "status and power don't match"
+            if !isapprox(var_val, 0.0; atol = ABSOLUTE_TOLERANCE)
+                error("Status and Power variables don't match for $comp_name. \\
+                ActivePowerVariable: $(var_val)\\
+                Status value: $(status_val) for OnVariable")
+            end
             set_ic_quantity!(ic, 0.0)
         end
     end
diff --git a/test/test_model_emulation.jl b/test/test_model_emulation.jl
index f824c7ec4b..9688073310 100644
--- a/test/test_model_emulation.jl
+++ b/test/test_model_emulation.jl
@@ -223,7 +223,7 @@ end
             model;
             executions = 10,
             output_dir = mktempdir(; cleanup = true),
-            serialize = serialize,
+            export_optimization_model = serialize,
         ) == PSI.RunStatus.SUCCESSFULLY_FINALIZED
     end
 end
diff --git a/test/test_network_constructors.jl b/test/test_network_constructors.jl
index 62be7f13f6..28e40011c8 100644
--- a/test/test_network_constructors.jl
+++ b/test/test_network_constructors.jl
@@ -964,7 +964,7 @@ end
         PSI.get_constraint(opt_container, CopperPlateBalanceConstraint(), PSY.Area)
     @test size(copper_plate_constraints) == (2, 24)
 
-    psi_checksolve_test(ps_model, [MOI.OPTIMAL], 662467, 1)
+    psi_checksolve_test(ps_model, [MOI.OPTIMAL], 671937, 1)
 
     results = OptimizationProblemResults(ps_model)
     interarea_flow = read_variable(results, "FlowActivePowerVariable__AreaInterchange")
diff --git a/test/test_simulation_build.jl b/test/test_simulation_build.jl
index 0b3f2e4d7c..db4b9e6be2 100644
--- a/test/test_simulation_build.jl
+++ b/test/test_simulation_build.jl
@@ -212,3 +212,188 @@ end
     )
     @test !isempty(c)
 end
+
+@testset "Test Upper/Lower Bound Feedforwards" begin
+    template_uc = get_template_basic_uc_simulation()
+    set_network_model!(template_uc, NetworkModel(PTDFPowerModel; use_slacks = true))
+    set_device_model!(template_uc, DeviceModel(Line, StaticBranchUnbounded))
+    template_ed =
+        get_template_nomin_ed_simulation(NetworkModel(PTDFPowerModel; use_slacks = true))
+    set_device_model!(template_ed, DeviceModel(Line, StaticBranchUnbounded))
+    c_sys5_hy_uc = PSB.build_system(PSITestSystems, "c_sys5_hy_uc")
+    c_sys5_hy_ed = PSB.build_system(PSITestSystems, "c_sys5_hy_ed")
+    models = SimulationModels(;
+        decision_models = [
+            DecisionModel(
+                template_uc,
+                c_sys5_hy_uc;
+                name = "UC",
+                optimizer = HiGHS_optimizer,
+                initialize_model = false,
+            ),
+            DecisionModel(
+                template_ed,
+                c_sys5_hy_ed;
+                name = "ED",
+                optimizer = ipopt_optimizer,
+                initialize_model = false,
+            ),
+        ],
+    )
+
+    sequence = SimulationSequence(;
+        models = models,
+        feedforwards = Dict(
+            "ED" => [
+                SemiContinuousFeedforward(;
+                    component_type = ThermalStandard,
+                    source = OnVariable,
+                    affected_values = [ActivePowerVariable],
+                ),
+                LowerBoundFeedforward(;
+                    component_type = Line,
+                    source = FlowActivePowerVariable,
+                    affected_values = [FlowActivePowerVariable],
+                    add_slacks = true,
+                ),
+                UpperBoundFeedforward(;
+                    component_type = Line,
+                    source = FlowActivePowerVariable,
+                    affected_values = [FlowActivePowerVariable],
+                    add_slacks = true,
+                ),
+            ],
+        ),
+        ini_cond_chronology = InterProblemChronology(),
+    )
+
+    sim = Simulation(;
+        name = "reactive_feedforward",
+        steps = 2,
+        models = models,
+        sequence = sequence,
+        simulation_folder = mktempdir(; cleanup = true),
+    )
+    build_out = build!(sim)
+    @test build_out == PSI.SimulationBuildStatus.BUILT
+    ed_power_model = PSI.get_simulation_model(PSI.get_models(sim), :ED)
+    c = PSI.get_constraint(
+        PSI.get_optimization_container(ed_power_model),
+        FeedforwardSemiContinuousConstraint(),
+        ThermalStandard,
+        "ActivePowerVariable_ub",
+    )
+    @test !isempty(c)
+    c = PSI.get_constraint(
+        PSI.get_optimization_container(ed_power_model),
+        FeedforwardSemiContinuousConstraint(),
+        ThermalStandard,
+        "ActivePowerVariable_lb",
+    )
+    @test !isempty(c)
+    c = PSI.get_constraint(
+        PSI.get_optimization_container(ed_power_model),
+        FeedforwardLowerBoundConstraint(),
+        Line,
+        "FlowActivePowerVariablelb",
+    )
+    @test !isempty(c)
+    c = PSI.get_constraint(
+        PSI.get_optimization_container(ed_power_model),
+        FeedforwardUpperBoundConstraint(),
+        Line,
+        "FlowActivePowerVariableub",
+    )
+    @test !isempty(c)
+    c = PSI.get_variable(
+        PSI.get_optimization_container(ed_power_model),
+        UpperBoundFeedForwardSlack(),
+        Line,
+        "FlowActivePowerVariable",
+    )
+    @test !isempty(c)
+    c = PSI.get_variable(
+        PSI.get_optimization_container(ed_power_model),
+        LowerBoundFeedForwardSlack(),
+        Line,
+        "FlowActivePowerVariable",
+    )
+    @test !isempty(c)
+end
+
+@testset "Test FixValue Feedforwards" begin
+    template_uc = get_template_basic_uc_simulation()
+    set_network_model!(template_uc, NetworkModel(PTDFPowerModel; use_slacks = true))
+    set_device_model!(template_uc, DeviceModel(Line, StaticBranchUnbounded))
+    set_service_model!(template_uc, ServiceModel(VariableReserve{ReserveUp}, RangeReserve))
+    template_ed =
+        get_template_nomin_ed_simulation(NetworkModel(PTDFPowerModel; use_slacks = true))
+    set_device_model!(template_ed, DeviceModel(Line, StaticBranchUnbounded))
+    set_service_model!(template_ed, ServiceModel(VariableReserve{ReserveUp}, RangeReserve))
+    c_sys5_hy_uc = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true)
+    c_sys5_hy_ed = PSB.build_system(PSITestSystems, "c_sys5_ed"; add_reserves = true)
+    models = SimulationModels(;
+        decision_models = [
+            DecisionModel(
+                template_uc,
+                c_sys5_hy_uc;
+                name = "UC",
+                optimizer = HiGHS_optimizer,
+                initialize_model = false,
+            ),
+            DecisionModel(
+                template_ed,
+                c_sys5_hy_ed;
+                name = "ED",
+                optimizer = ipopt_optimizer,
+                initialize_model = false,
+            ),
+        ],
+    )
+
+    sequence = SimulationSequence(;
+        models = models,
+        feedforwards = Dict(
+            "ED" => [
+                SemiContinuousFeedforward(;
+                    component_type = ThermalStandard,
+                    source = OnVariable,
+                    affected_values = [ActivePowerVariable],
+                ),
+                FixValueFeedforward(;
+                    component_type = VariableReserve{ReserveUp},
+                    source = ActivePowerReserveVariable,
+                    affected_values = [ActivePowerReserveVariable],
+                ),
+            ],
+        ),
+        ini_cond_chronology = InterProblemChronology(),
+    )
+
+    sim = Simulation(;
+        name = "reserve_feedforward",
+        steps = 2,
+        models = models,
+        sequence = sequence,
+        simulation_folder = mktempdir(; cleanup = true),
+    )
+    build_out = build!(sim)
+    @test build_out == PSI.SimulationBuildStatus.BUILT
+    ed_power_model = PSI.get_simulation_model(PSI.get_models(sim), :ED)
+    c = PSI.get_parameter(
+        PSI.get_optimization_container(ed_power_model),
+        FixValueParameter(),
+        VariableReserve{ReserveUp},
+        "Reserve1",
+    )
+    @test !isempty(c.multiplier_array)
+    @test !isempty(c.parameter_array)
+    c = PSI.get_parameter(
+        PSI.get_optimization_container(ed_power_model),
+        FixValueParameter(),
+        VariableReserve{ReserveUp},
+        "Reserve11",
+    )
+    @test !isempty(c.multiplier_array)
+    @test !isempty(c.parameter_array)
+end
diff --git a/test/test_simulation_execute.jl b/test/test_simulation_execute.jl
index 6bc4581769..22b8cbb251 100644
--- a/test/test_simulation_execute.jl
+++ b/test/test_simulation_execute.jl
@@ -1,4 +1,5 @@
-function test_single_stage_sequential(in_memory, rebuild)
+function test_single_stage_sequential(in_memory, rebuild, export_model)
+    tmp_dir = mktempdir(; cleanup = true)
     template_ed = get_template_nomin_ed_simulation()
     c_sys = PSB.build_system(PSITestSystems, "c_sys5_uc")
     models = SimulationModels([
@@ -8,6 +9,7 @@ function test_single_stage_sequential(in_memory, rebuild)
             name = "ED",
             optimizer = ipopt_optimizer,
             rebuild_model = rebuild,
+            export_optimization_model = export_model,
         ),
     ])
     test_sequence =
@@ -20,20 +22,29 @@ function test_single_stage_sequential(in_memory, rebuild)
         steps = 2,
         models = models,
         sequence = test_sequence,
-        simulation_folder = mktempdir(; cleanup = true),
+        simulation_folder = tmp_dir,
     )
     build_out = build!(sim_single)
     @test build_out == PSI.SimulationBuildStatus.BUILT
     execute_out = execute!(sim_single; in_memory = in_memory)
     @test execute_out == PSI.RunStatus.SUCCESSFULLY_FINALIZED
+    return tmp_dir
 end
 
 @testset "Single stage sequential tests" begin
     for in_memory in (true, false), rebuild in (true, false)
-        test_single_stage_sequential(in_memory, rebuild)
+        test_single_stage_sequential(in_memory, rebuild, false)
     end
 end
 
+@testset "Test model export at each solve" begin
+    folder = test_single_stage_sequential(true, false, true)
+    test_path =
+        joinpath(folder, "consecutive", "problems", "ED", "optimization_model_exports")
+    @test ispath(test_path)
+    @test length(readdir(test_path)) == 2
+end
+
 function test_2_stage_decision_models_with_feedforwards(in_memory)
     template_uc = get_template_basic_uc_simulation()
     template_ed = get_template_nomin_ed_simulation()