diff --git a/src/InfrastructureOptimizationModels.jl b/src/InfrastructureOptimizationModels.jl
index 1fa5f86..b7b500e 100644
--- a/src/InfrastructureOptimizationModels.jl
+++ b/src/InfrastructureOptimizationModels.jl
@@ -59,9 +59,7 @@ import InfrastructureSystems.Optimization:
     AbstractStorageFormulation,
     AbstractLoadFormulation,
     AbstractHVDCNetworkModel,
-    AbstractPowerModel,
-    AbstractPowerFlowEvaluationModel,
-    AbstractPowerFlowEvaluationData
+    AbstractPowerModel
 
 import InfrastructureSystems:
     @scoped_enum,
@@ -170,7 +168,7 @@ export InitialCondition
 export NetworkModel
 export get_PTDF_matrix, get_LODF_matrix, get_reduce_radial_branches
 export get_duals, get_reference_buses, get_subnetworks, get_bus_area_map
-export get_power_flow_evaluation, has_subnetworks, get_subsystem
+export get_evaluations, has_subnetworks, get_subsystem
 export set_subsystem!, add_dual!
 export requires_all_branch_models, supports_branch_filtering, ignores_branch_filtering
 export validate_network_model
@@ -223,6 +221,7 @@ export add_constant_to_jump_expression!
 export add_proportional_to_jump_expression!
 export add_linear_to_jump_expression!
 # Cost term helpers (generic objective function building blocks)
+export add_cost_term_to_expression!
 export add_cost_term_invariant!
 export add_cost_term_variant!
 export add_pwl_variables_delta!
@@ -376,7 +375,7 @@ export get_incompatible_devices
 
 # Bulk export: symbols POM needs that weren't previously exported
 # Core types
-export OptimizationContainer, OperationModel, AbstractPowerFlowEvaluationModel
+export OptimizationContainer, OperationModel
 export ArgumentConstructStage, ModelConstructStage
 export EmulationModelStore, DeviceModelForBranches
 export SOSStatusVariable
@@ -443,6 +442,13 @@ export EmergencyUp
 export EmergencyDown
 export RawACE
 export ProductionCostExpression
+export ConstituentCostExpression
+export FuelCostExpression
+export StartUpCostExpression
+export ShutDownCostExpression
+export FixedCostExpression
+export VOMCostExpression
+export CurtailmentCostExpression
 export FuelConsumptionExpression
 export ActivePowerRangeExpressionLB
 export ActivePowerRangeExpressionUB
@@ -468,8 +474,25 @@ export SparseVariableType, InterpolationVariableType, BinaryInterpolationVariabl
 export UpperBound, LowerBound, BoundDirection, get_bound_direction
 export EventParameter
 
-# Abstract types for extensions (from InfrastructureSystems.Optimization)
-export AbstractPowerFlowEvaluationData
+# External evaluation abstraction (replaces the PowerFlows-specific shims)
+export AbstractEvaluator
+export AbstractEvaluationData
+export EvaluationContainer
+export initialize_evaluation_data
+export evaluate!
+export reset!
+export is_solved
+export get_evaluation_data
+export get_inner_data
+export get_evaluators
+export get_evaluator
+export add_evaluator!
+export add_evaluation_data!
+export reset_evaluations!
+export is_from_evaluator
+export lookup_value
+export get_entry_type
+export get_component_names
 
 # Status Enums (from InfrastructureSystems)
 export ModelBuildStatus
@@ -549,6 +572,7 @@ include("core/operation_model_abstract_types.jl")
 include("core/network_reductions.jl")
 include("core/service_model.jl")
 include("core/device_model.jl")
+include("core/external_evaluation.jl")
 include("core/network_model.jl")
 include("core/initial_conditions.jl")
 include("core/settings.jl")
diff --git a/src/common_models/interfaces.jl b/src/common_models/interfaces.jl
index 1fa1ad8..538793a 100644
--- a/src/common_models/interfaces.jl
+++ b/src/common_models/interfaces.jl
@@ -203,12 +203,6 @@ function start_up_cost(
     )
 end
 
-"""
-Extension point: Solve the power flow model.
-Default: error. Concrete implementations require PowerFlows integration.
-"""
-function solve_powerflow! end
-
 """
 Extension point: Calculate auxiliary variable values.
 Concrete implementations in PowerOperationsModels for specific aux variable types.
@@ -216,10 +210,11 @@ Concrete implementations in PowerOperationsModels for specific aux variable type
 function calculate_aux_variable_value! end
 
 """
-Extension point: Check if an auxiliary variable type comes from power flow evaluation.
-Default: false. Override in POM for PowerFlowAuxVariableType subtypes.
+Extension point: Check if an auxiliary variable type is produced by an external
+evaluator (see `AbstractEvaluator`). Default: false. Override in POM for
+evaluator-bound aux variable subtypes.
 """
-is_from_power_flow(::Type{<:AuxVariableType}) = false
+is_from_evaluator(::Type{<:AuxVariableType}) = false
 
 """
 Extension point: Get minimum and maximum limits for a given component, constraint type, and device formulation.
diff --git a/src/core/external_evaluation.jl b/src/core/external_evaluation.jl
new file mode 100644
index 0000000..30acc09
--- /dev/null
+++ b/src/core/external_evaluation.jl
@@ -0,0 +1,141 @@
+"""
+Abstract type for evaluator configurations attached to a `NetworkModel`.
+
+An `AbstractEvaluator` is stateless configuration that, when used to build an
+`OptimizationContainer`, produces an `AbstractEvaluationData` (its runtime
+counterpart) via [`initialize_evaluation_data`](@ref).
+
+Concrete subtypes (e.g. `PowerFlowEvaluationModel`) are defined in downstream
+packages such as PowerOperationsModels and its PowerFlows extension.
+"""
+abstract type AbstractEvaluator end
+
+"""
+Abstract type for runtime evaluator state stored inside the
+`OptimizationContainer`'s [`EvaluationContainer`](@ref).
+
+Concrete subtypes wrap whatever domain-specific data (power-flow solver state,
+exporter handles, etc.) the evaluator needs across iterations of
+`calculate_aux_variables!`.
+"""
+abstract type AbstractEvaluationData end
+
+"""
+Holds the evaluators registered on a `NetworkModel` and their associated
+runtime data, keyed by concrete evaluator type.
+
+`evaluators` is populated by the user when constructing a `NetworkModel`;
+`evaluation_data` is populated during `OptimizationContainer` construction by
+calling [`initialize_evaluation_data`](@ref) for each registered evaluator.
+
+The split between the two dictionaries may evolve; for now we keep them as
+parallel keyed dictionaries so that `evaluation_data[T]` is the runtime state
+produced from `evaluators[T]`.
+"""
+mutable struct EvaluationContainer
+    evaluators::Dict{DataType, Any}
+    evaluation_data::Dict{DataType, AbstractEvaluationData}
+end
+
+function EvaluationContainer()
+    return EvaluationContainer(
+        Dict{DataType, Any}(),
+        Dict{DataType, AbstractEvaluationData}(),
+    )
+end
+
+get_evaluators(ec::EvaluationContainer) = ec.evaluators
+get_evaluation_data(ec::EvaluationContainer) = ec.evaluation_data
+
+get_evaluator(ec::EvaluationContainer, T::DataType) = ec.evaluators[T]
+get_evaluation_data(ec::EvaluationContainer, T::DataType) = ec.evaluation_data[T]
+
+add_evaluator!(ec::EvaluationContainer, T::DataType, ev) =
+    (ec.evaluators[T] = ev)
+add_evaluation_data!(
+    ec::EvaluationContainer,
+    T::DataType,
+    d::AbstractEvaluationData,
+) = (ec.evaluation_data[T] = d)
+
+Base.isempty(ec::EvaluationContainer) = isempty(ec.evaluators)
+Base.length(ec::EvaluationContainer) = length(ec.evaluators)
+Base.haskey(ec::EvaluationContainer, T::DataType) = haskey(ec.evaluators, T)
+
+#=
+=========================================================================
+AbstractEvaluator interface (config-side)
+=========================================================================
+=#
+
+"""
+Build the runtime state for `ev` and return an `AbstractEvaluationData`.
+Called once during `OptimizationContainer` construction. Concrete methods are
+registered in downstream packages.
+"""
+function initialize_evaluation_data(
+    ev::AbstractEvaluator,
+    container,
+    system,
+)
+    error(
+        "initialize_evaluation_data not implemented for evaluator type $(typeof(ev)). " *
+        "Define `initialize_evaluation_data(::$(typeof(ev)), container, system) ::AbstractEvaluationData` " *
+        "in the package that owns the concrete type.",
+    )
+end
+
+#=
+=========================================================================
+AbstractEvaluationData interface (runtime-side)
+=========================================================================
+=#
+
+"""
+Run the evaluation: read state from `container`/`system`, write aux-variable
+inputs back into `container`, and mark `data` as solved. Concrete methods are
+registered in downstream packages.
+"""
+function evaluate!(data::AbstractEvaluationData, container, system)
+    error(
+        "evaluate! not implemented for evaluation data type $(typeof(data)). " *
+        "Define `evaluate!(::$(typeof(data)), container, system)` " *
+        "in the package that owns the concrete type.",
+    )
+end
+
+"""
+Mark `data` as not-yet-solved for the current iteration. Called by
+`calculate_aux_variables!` at the top of each pass.
+"""
+function reset!(data::AbstractEvaluationData)
+    error(
+        "reset! not implemented for evaluation data type $(typeof(data)). " *
+        "Define `reset!(::$(typeof(data)))` " *
+        "in the package that owns the concrete type.",
+    )
+end
+
+"""
+Return `true` if `data` has been solved in the current iteration.
+"""
+function is_solved(data::AbstractEvaluationData)
+    error(
+        "is_solved not implemented for evaluation data type $(typeof(data)). " *
+        "Define `is_solved(::$(typeof(data))) ::Bool` " *
+        "in the package that owns the concrete type.",
+    )
+end
+
+"""
+Return the underlying domain-specific data wrapped by `data` (e.g., a
+`PowerFlowData` for a power-flow evaluator). Used by aux-variable calculators
+that need direct read access to the evaluator's output.
+"""
+function get_inner_data(data::AbstractEvaluationData)
+    error(
+        "get_inner_data not implemented for evaluation data type $(typeof(data)). " *
+        "Define `get_inner_data(::$(typeof(data)))` " *
+        "in the package that owns the concrete type.",
+    )
+end
diff --git a/src/core/network_model.jl b/src/core/network_model.jl
index 928fa85..f619460 100644
--- a/src/core/network_model.jl
+++ b/src/core/network_model.jl
@@ -12,11 +12,6 @@ function _check_pm_formulation(::Type{T}) where {T <: AbstractPowerModel}
     end
 end
 
-# Note: PowerFlowEvaluationModel support has been moved to PowerOperationsModels
-# These stub functions maintain API compatibility
-_maybe_flatten_pfem(pfem::Vector) = pfem
-_maybe_flatten_pfem(pfem) = [pfem]
-
 """
 Establishes the NetworkModel for a given AC network formulation type.
 
@@ -39,8 +34,9 @@ Establishes the NetworkModel for a given AC network formulation type.
     subnetworks are inferred from PTDF/VirtualPTDF or discovered from the system.
 - `duals::Vector{DataType}` = Vector{DataType}()
     Constraint types for which duals should be recorded.
-- `power_flow_evaluation::Union{AbstractPowerFlowEvaluationModel, Vector{<:AbstractPowerFlowEvaluationModel}}`
-    Power-flow evaluation model(s). A single model is flattened to a vector internally.
+- `evaluations::EvaluationContainer`
+    External evaluators (e.g. power-flow) keyed by concrete evaluator type.
+    Default is an empty container — no evaluator runs.
 
 # Notes
 - `modeled_branch_types` and `reduced_branch_tracker` are internal fields managed by the model.
@@ -51,8 +47,10 @@ Establishes the NetworkModel for a given AC network formulation type.
 
 # Examples (concrete types like PTDFPowerModel, CopperPlatePowerModel are defined in PowerSimulations)
 # ptdf = PNM.VirtualPTDF(system)
+# ec = EvaluationContainer()
+# add_evaluator!(ec, PFS.PowerFlowEvaluationModel, PFS.PowerFlowEvaluationModel())
 # nw = NetworkModel(PTDFPowerModel; PTDF_matrix = ptdf, reduce_radial_branches = true,
-#                   power_flow_evaluation = PFS.PowerFlowEvaluationModel())
+#                   evaluations = ec)
 #
 # nw2 = NetworkModel(CopperPlatePowerModel; subnetworks = Dict(1 => Set([1,2,3])))
 """
@@ -66,7 +64,7 @@ mutable struct NetworkModel{T <: AbstractPowerModel}
     network_reduction::PNM.NetworkReductionData
     reduce_radial_branches::Bool
     reduce_degree_two_branches::Bool
-    power_flow_evaluation::Vector{<:AbstractPowerFlowEvaluationModel}
+    evaluations::EvaluationContainer
     subsystem::Union{Nothing, String}
     hvdc_network_model::Union{Nothing, AbstractHVDCNetworkModel}
     modeled_branch_types::Vector{DataType}
@@ -81,10 +79,7 @@ mutable struct NetworkModel{T <: AbstractPowerModel}
         reduce_degree_two_branches = false,
         subnetworks = Dict{Int, Set{Int}}(),
         duals = Vector{DataType}(),
-        power_flow_evaluation::Union{
-            AbstractPowerFlowEvaluationModel,
-            Vector{<:AbstractPowerFlowEvaluationModel},
-        } = AbstractPowerFlowEvaluationModel[],
+        evaluations = EvaluationContainer(),
         hvdc_network_model = nothing,
     ) where {T <: AbstractPowerModel}
         _check_pm_formulation(T)
@@ -98,7 +93,7 @@ mutable struct NetworkModel{T <: AbstractPowerModel}
             PNM.NetworkReductionData(),
             reduce_radial_branches,
             reduce_degree_two_branches,
-            _maybe_flatten_pfem(power_flow_evaluation),
+            evaluations,
             nothing,
             hvdc_network_model,
             Vector{DataType}(),
@@ -119,7 +114,7 @@ get_reference_buses(m::NetworkModel{T}) where {T <: AbstractPowerModel} =
     collect(keys(m.subnetworks))
 get_subnetworks(m::NetworkModel) = m.subnetworks
 get_bus_area_map(m::NetworkModel) = m.bus_area_map
-get_power_flow_evaluation(m::NetworkModel) = m.power_flow_evaluation
+get_evaluations(m::NetworkModel) = m.evaluations
 has_subnetworks(m::NetworkModel) = !isempty(m.bus_area_map)
 get_subsystem(m::NetworkModel) = m.subsystem
 get_hvdc_network_model(m::NetworkModel) = m.hvdc_network_model
diff --git a/src/core/optimization_container.jl b/src/core/optimization_container.jl
index 36f1979..598fbfd 100644
--- a/src/core/optimization_container.jl
+++ b/src/core/optimization_container.jl
@@ -79,11 +79,11 @@ mutable struct OptimizationContainer <: AbstractOptimizationContainer
     model_base_power::Float64
     optimizer_stats::OptimizerStats
     built_for_recurrent_solves::Bool
-    pf_aux_var_keys::Vector{AuxVarKey}
-    non_pf_aux_var_keys::Vector{AuxVarKey}
+    evaluator_aux_var_keys::Vector{AuxVarKey}
+    standalone_aux_var_keys::Vector{AuxVarKey}
     metadata::OptimizationContainerMetadata
     default_time_series_type::Type
-    power_flow_evaluation_data::Vector{<:AbstractPowerFlowEvaluationData}
+    evaluations::EvaluationContainer
 end
 
 function OptimizationContainer(
@@ -127,7 +127,7 @@ function OptimizationContainer(
         AuxVarKey[],
         OptimizationContainerMetadata(),
         T,
-        AbstractPowerFlowEvaluationData[],
+        EvaluationContainer(),
     )
 end
 
@@ -165,8 +165,7 @@ get_jump_model(container::OptimizationContainer) = container.JuMPmodel
 get_metadata(container::OptimizationContainer) = container.metadata
 get_optimizer_stats(container::OptimizationContainer) = container.optimizer_stats
 get_parameters(container::OptimizationContainer) = container.parameters
-get_power_flow_evaluation_data(container::OptimizationContainer) =
-    container.power_flow_evaluation_data
+get_evaluations(container::OptimizationContainer) = container.evaluations
 get_resolution(container::OptimizationContainer) = get_resolution(container.settings)
 get_settings(container::OptimizationContainer) = container.settings
 get_time_steps(container::OptimizationContainer) = container.time_steps
@@ -180,9 +179,9 @@ is_synchronized(container::OptimizationContainer) =
 set_time_steps!(container::OptimizationContainer, time_steps::UnitRange{Int64}) =
     container.time_steps = time_steps
 
-function reset_power_flow_is_solved!(container::OptimizationContainer)
-    for pf_e_data in get_power_flow_evaluation_data(container)
-        pf_e_data.is_solved = false
+function reset_evaluations!(container::OptimizationContainer)
+    for data in values(get_evaluation_data(get_evaluations(container)))
+        reset!(data)
     end
 end
 
@@ -352,8 +351,8 @@ function reset_optimization_model!(container::OptimizationContainer)
     for field in [:variables, :aux_variables, :constraints, :expressions, :duals]
         empty!(getfield(container, field))
     end
-    empty!(container.pf_aux_var_keys)
-    empty!(container.non_pf_aux_var_keys)
+    empty!(container.evaluator_aux_var_keys)
+    empty!(container.standalone_aux_var_keys)
     container.initial_conditions_data = InitialConditionsData()
     container.objective_function = ObjectiveFunction()
     container.primal_values_cache = PrimalValuesCache()
@@ -1230,9 +1229,9 @@ end
 
 function _cache_aux_variable_key_partitions!(container::OptimizationContainer)
     aux_var_keys = keys(get_aux_variables(container))
-    pf_keys = filter(is_from_power_flow ∘ get_entry_type, aux_var_keys)
-    container.pf_aux_var_keys = collect(pf_keys)
-    container.non_pf_aux_var_keys = collect(setdiff(aux_var_keys, pf_keys))
+    evaluator_keys = filter(is_from_evaluator ∘ get_entry_type, aux_var_keys)
+    container.evaluator_aux_var_keys = collect(evaluator_keys)
+    container.standalone_aux_var_keys = collect(setdiff(aux_var_keys, evaluator_keys))
     return
 end
 
@@ -1240,28 +1239,30 @@ function calculate_aux_variables!(
     container::OptimizationContainer,
     system::IS.InfrastructureSystemsContainer,
 )
-    if isempty(container.pf_aux_var_keys) && isempty(container.non_pf_aux_var_keys)
+    if isempty(container.evaluator_aux_var_keys) &&
+       isempty(container.standalone_aux_var_keys)
         _cache_aux_variable_key_partitions!(container)
     end
-    pf_aux_var_keys = container.pf_aux_var_keys
-    non_pf_aux_var_keys = container.non_pf_aux_var_keys
-    # We should only have power flow aux vars if we have power flow evaluators
-    @assert isempty(pf_aux_var_keys) || !isempty(get_power_flow_evaluation_data(container))
-
-    TimerOutputs.@timeit RUN_SIMULATION_TIMER "Power Flow Evaluation" begin
-        reset_power_flow_is_solved!(container)
-        # Power flow-related aux vars get calculated once per power flow
-        for (i, pf_e_data) in enumerate(get_power_flow_evaluation_data(container))
-            @debug "Processing power flow $i"
-            solve_powerflow!(pf_e_data, container)
-            for key in pf_aux_var_keys
+    evaluator_aux_var_keys = container.evaluator_aux_var_keys
+    standalone_aux_var_keys = container.standalone_aux_var_keys
+    evaluation_data = get_evaluation_data(get_evaluations(container))
+    # We should only have evaluator-bound aux vars if we have evaluators registered
+    @assert isempty(evaluator_aux_var_keys) || !isempty(evaluation_data)
+
+    TimerOutputs.@timeit RUN_SIMULATION_TIMER "External Evaluation" begin
+        reset_evaluations!(container)
+        # Evaluator-bound aux vars get calculated once per evaluator
+        for (T, data) in evaluation_data
+            @debug "Processing evaluator $(T)"
+            evaluate!(data, container, system)
+            for key in evaluator_aux_var_keys
                 calculate_aux_variable_value!(container, key, system)
             end
         end
     end
 
     # Other aux vars get calculated once at the end
-    for key in non_pf_aux_var_keys
+    for key in standalone_aux_var_keys
         calculate_aux_variable_value!(container, key, system)
     end
     return RunStatus.SUCCESSFULLY_FINALIZED
diff --git a/src/core/parameter_container.jl b/src/core/parameter_container.jl
index d362c46..aa71385 100644
--- a/src/core/parameter_container.jl
+++ b/src/core/parameter_container.jl
@@ -108,7 +108,15 @@ function ParameterContainer(parameter_array, multiplier_array)
 end
 
 get_parameter_array(c::ParameterContainer) = c.parameter_array
+# Underlying dense storage of the parameter array. `parent` on a JuMP `DenseAxisArray`
+# returns the array itself, so reach for `.data` directly to bypass the axis-keyed lookup.
+# Restricted to dense storage: `SparseAxisArray` has no `.data` field.
+get_parameter_array_data(c::ParameterContainer{<:DenseAxisArray}) =
+    get_parameter_array(c).data
 get_multiplier_array(c::ParameterContainer) = c.multiplier_array
+# Same shortcut for the multiplier array — used by the integer-indexed fast path.
+get_multiplier_array_data(c::ParameterContainer{<:Any, <:DenseAxisArray}) =
+    get_multiplier_array(c).data
 get_attributes(c::ParameterContainer) = c.attributes
 Base.length(c::ParameterContainer) = length(c.parameter_array)
 Base.size(c::ParameterContainer) = size(c.parameter_array)
@@ -281,3 +289,178 @@ function set_parameter!(
     assign_maybe_broadcast!(get_parameter_array(container), parameter, ixs)
     return
 end
+
+# Fast-path setters that skip DenseAxisArray's string-keyed axis lookup. Callers pass
+# `get_parameter_array_data(container)` once, then write into the underlying Array
+# by integer indices. The (i, t) layout matches the canonical (component, time) axis
+# order produced by `add_param_container!`.
+#
+# 2D scalar path: covers Float64 and Tuple{Vararg{Float64}} eltypes (the latter is
+# used by piecewise-cost MarketBid parameters whose storage is a Matrix of tuples).
+@inline function _set_parameter_at!(
+    parent_param::Array{T, 2},
+    ::JuMP.Model,
+    value::T,
+    i::Int,
+    t::Int,
+) where {T <: ValidDataParamEltypes}
+    parent_param[i, t] = value
+    return
+end
+
+# 2D recurrent-rebuild paths: param storage is `Array{JuMP.VariableRef, 2}`. Either we
+# need a fresh JuMP parameter (Float64 input) or we reuse one created by an earlier
+# parallel branch type (VariableRef input).
+@inline function _set_parameter_at!(
+    parent_param::Array{JuMP.VariableRef, 2},
+    jump_model::JuMP.Model,
+    value::Float64,
+    i::Int,
+    t::Int,
+)
+    parent_param[i, t] = add_jump_parameter(jump_model, value)
+    return
+end
+
+@inline function _set_parameter_at!(
+    parent_param::Array{JuMP.VariableRef, 2},
+    ::JuMP.Model,
+    parameter::JuMP.VariableRef,
+    i::Int,
+    t::Int,
+)
+    parent_param[i, t] = parameter
+    return
+end
+
+# 3D fast paths (parameter container with a middle additional axis, e.g. piecewise
+# tranches). The supplied `value` is a length-(size(parent_param, 2)) vector that fills
+# the middle axis at position (i, :, t). Eltype constrained to `ValidDataParamEltypes`
+# so tuples-of-floats are also accepted (piecewise breakpoint storage).
+@inline function _set_parameter_at!(
+    parent_param::Array{T, 3},
+    ::JuMP.Model,
+    value::AbstractVector{<:T},
+    i::Int,
+    t::Int,
+) where {T <: ValidDataParamEltypes}
+    @views parent_param[i, :, t] .= value
+    return
+end
+
+@inline function _set_parameter_at!(
+    parent_param::Array{JuMP.VariableRef, 3},
+    jump_model::JuMP.Model,
+    value::AbstractVector{Float64},
+    i::Int,
+    t::Int,
+)
+    for k in 1:size(parent_param, 2)
+        parent_param[i, k, t] = add_jump_parameter(jump_model, value[k])
+    end
+    return
+end
+
+# Fast-path setters for the multiplier array, mirroring `_set_parameter_at!`.
+# Multipliers are always Float64-valued (or tuples-of-floats for piecewise
+# parameters), so a single typed family covers every call site. Callers should
+# hoist `parent_mult = get_multiplier_array_data(parameter_container)` once
+# above the device loop and pass the integer device row index.
+
+# 2D row fill: assigns `value` across the whole time slice for component `i`
+# (the canonical pattern at parameter-creation time, where the multiplier is
+# constant per device).
+@inline function _set_multiplier_at!(
+    parent_mult::Array{T, 2},
+    value::T,
+    i::Int,
+) where {T <: ValidDataParamEltypes}
+    @views parent_mult[i, :] .= value
+    return
+end
+
+# 2D scalar write at a single (component, time) cell.
+@inline function _set_multiplier_at!(
+    parent_mult::Array{T, 2},
+    value::T,
+    i::Int,
+    t::Int,
+) where {T <: ValidDataParamEltypes}
+    parent_mult[i, t] = value
+    return
+end
+
+# 3D row fill: assigns `value` across all (tranche, time) for component `i`.
+@inline function _set_multiplier_at!(
+    parent_mult::Array{T, 3},
+    value::T,
+    i::Int,
+) where {T <: ValidDataParamEltypes}
+    @views parent_mult[i, :, :] .= value
+    return
+end
+
+# 3D point write at a single (component, tranche, time) cell.
+@inline function _set_multiplier_at!(
+    parent_mult::Array{T, 3},
+    value::T,
+    i::Int,
+    j::Int,
+    t::Int,
+) where {T <: ValidDataParamEltypes}
+    parent_mult[i, j, t] = value
+    return
+end
+
+# Fast-path setters for the simulation-step parameter-VALUE update path.
+# Used by `_update_parameter_values!` and `_fix_parameter_value!` overloads
+# where the parameter container already exists and we are pushing new values
+# into it from upstream model results or time series. Caller hoists
+# `parent_param = get_parameter_array_data(parameter_container)` (and an
+# optional `parent_mult` / `parent_var`) above the device loop, then writes
+# by integer (i, t) — bypassing DenseAxisArray's string-keyed axis lookup.
+#
+# For Float64-typed storage we write directly; for `JuMP.VariableRef` storage
+# we update the JuMP parameter's bound via `JuMP.fix(...; force=true)`.
+@inline function _set_param_value_at!(
+    parent_param::Array{T, 2},
+    value::T,
+    i::Int,
+    t::Int,
+) where {T <: ValidDataParamEltypes}
+    @inbounds parent_param[i, t] = value
+    return
+end
+
+@inline function _set_param_value_at!(
+    parent_param::Array{JuMP.VariableRef, 2},
+    value::Float64,
+    i::Int,
+    t::Int,
+)
+    @inbounds JuMP.fix(parent_param[i, t], value; force = true)
+    return
+end
+
+# 3D paths for piecewise-tranche updates.
+@inline function _set_param_value_at!(
+    parent_param::Array{T, 3},
+    value::AbstractVector{<:T},
+    i::Int,
+    t::Int,
+) where {T <: ValidDataParamEltypes}
+    @inbounds @views parent_param[i, :, t] .= value
+    return
+end
+
+@inline function _set_param_value_at!(
+    parent_param::Array{JuMP.VariableRef, 3},
+    value::AbstractVector{Float64},
+    i::Int,
+    t::Int,
+)
+    @inbounds for k in 1:size(parent_param, 2)
+        JuMP.fix(parent_param[i, k, t], value[k]; force = true)
+    end
+    return
+end
diff --git a/src/core/standard_variables_expressions.jl b/src/core/standard_variables_expressions.jl
index 23aa953..58cf876 100644
--- a/src/core/standard_variables_expressions.jl
+++ b/src/core/standard_variables_expressions.jl
@@ -50,6 +50,13 @@ abstract type PostContingencyExpressions <: ExpressionType end
 
 # Concrete expression types used in IOM code
 struct ProductionCostExpression <: CostExpressions end
+abstract type ConstituentCostExpression <: CostExpressions end
+struct FuelCostExpression <: ConstituentCostExpression end
+struct StartUpCostExpression <: ConstituentCostExpression end
+struct ShutDownCostExpression <: ConstituentCostExpression end
+struct FixedCostExpression <: ConstituentCostExpression end
+struct VOMCostExpression <: ConstituentCostExpression end
+struct CurtailmentCostExpression <: CostExpressions end
 struct FuelConsumptionExpression <: ExpressionType end
 struct ActivePowerRangeExpressionLB <: RangeConstraintLBExpressions end
 struct ActivePowerRangeExpressionUB <: RangeConstraintUBExpressions end
@@ -75,8 +82,11 @@ should_write_resulting_value(::Type{ActivePowerBalance}) = true
 should_write_resulting_value(::Type{ReactivePowerBalance}) = true
 should_write_resulting_value(::Type{DCCurrentBalance}) = true
 
-# ProductionCostExpression-specific container method (moved here from optimization_container.jl
-# because it requires ProductionCostExpression to be defined first)
+# CostExpressions container method (moved here from optimization_container.jl
+# because it requires the cost expression types to be defined first). Covers
+# ProductionCostExpression, ConstituentCostExpression subtypes, and
+# CurtailmentCostExpression — all use a JuMP.QuadExpr container so quadratic
+# fuel terms can be stored in the same expression.
 function add_expression_container!(
     container::OptimizationContainer,
     ::Type{T},
@@ -85,7 +95,7 @@ function add_expression_container!(
     sparse = false,
     meta = CONTAINER_KEY_EMPTY_META,
 ) where {
-    T <: ProductionCostExpression,
+    T <: CostExpressions,
     U <: Union{IS.InfrastructureSystemsComponent, IS.InfrastructureSystemsContainer},
 }
     expr_container =
diff --git a/src/objective_function/common.jl b/src/objective_function/common.jl
index ae3c7fe..ff110db 100644
--- a/src/objective_function/common.jl
+++ b/src/objective_function/common.jl
@@ -20,6 +20,8 @@ function add_cost_to_expression!(
             cost_expression,
         )
     end
+    _propagate_to_production_cost!(
+        container, S, T, component_name, time_period, cost_expression)
     return
 end
 
@@ -69,7 +71,8 @@ function _add_vom_cost_to_objective!(
     variable_cost_data = variable_cost(op_cost, T, C, U)
     power_units = IS.get_power_units(variable_cost_data)
     cost_term = IS.get_proportional_term(IS.get_vom_cost(variable_cost_data))
-    add_proportional_cost_invariant!(container, T, component, cost_term, power_units)
+    add_proportional_cost_invariant!(
+        container, T, component, cost_term, power_units, 1.0, VOMCostExpression)
     return
 end
 
@@ -114,7 +117,7 @@ function _add_time_varying_fuel_variable_cost!(
             container,
             expression[name, t],
             FuelCostParameter,
-            ProductionCostExpression,
+            FuelCostExpression,
             V,
             name,
             t,
diff --git a/src/objective_function/cost_term_helpers.jl b/src/objective_function/cost_term_helpers.jl
index 2f3a85f..dd27d2d 100644
--- a/src/objective_function/cost_term_helpers.jl
+++ b/src/objective_function/cost_term_helpers.jl
@@ -7,26 +7,61 @@
 # - No concrete expression or parameter types - caller passes them
 # - Single timestep only - looping stays in PSI/POM
 
+#######################################
+###### Constituent Propagation ########
+#######################################
+
+# Default no-op: any expression type that is not a ConstituentCostExpression
+# does not propagate into the aggregate ProductionCostExpression.
+_propagate_to_production_cost!(
+    ::OptimizationContainer,
+    ::Type{<:ExpressionType},
+    ::Type{<:IS.InfrastructureSystemsComponent},
+    ::String,
+    ::Int,
+    ::Any,
+) = nothing
+
+# ConstituentCostExpression overload: also write the same cost into
+# ProductionCostExpression so the aggregate stays consistent with its parts.
+function _propagate_to_production_cost!(
+    container::OptimizationContainer,
+    ::Type{<:ConstituentCostExpression},
+    ::Type{C},
+    name::String,
+    t::Int,
+    cost,
+) where {C <: IS.InfrastructureSystemsComponent}
+    has_container_key(container, ProductionCostExpression, C) || return
+    prod_expr = get_expression(container, ProductionCostExpression, C)
+    JuMP.add_to_expression!(prod_expr[name, t], cost)
+    return
+end
+
 #######################################
 ######## Linear Cost Helpers ##########
 #######################################
 
 """
-Add cost term to expression and invariant objective.
+Add cost term to a target expression (no objective hook).
 
-Computes `cost = quantity * rate`, adds to target expression (if present),
-and adds to the time-invariant part of the objective.
+Computes `cost = quantity * rate`, adds it to expression `E` for component `C` at
+time `t` if that expression exists, and propagates into `ProductionCostExpression`
+when `E <: ConstituentCostExpression`. Use this when the caller wants to record
+the cost in an expression container without adding to the JuMP objective (e.g.,
+fuel consumption, where the term is a downstream quantity rather than a cost
+that should be minimized).
 
 # Arguments
 - `container`: the optimization container
 - `quantity`: the value being costed (e.g., variable value, expression value)
 - `rate`: scalar cost rate (e.g., \$/MWh, \$/MMBTU)
-- `E`: target expression type (caller provides, e.g., ProductionCostExpression)
+- `E`: target expression type (caller provides)
 - `C`: component type
 - `name`: component name
 - `t`: time period
 """
-function add_cost_term_invariant!(
+function add_cost_term_to_expression!(
     container::OptimizationContainer,
     quantity::JuMPOrFloat,
     rate::Float64,
@@ -40,6 +75,36 @@ function add_cost_term_invariant!(
         expr = get_expression(container, E, C)
         JuMP.add_to_expression!(expr[name, t], cost)
     end
+    _propagate_to_production_cost!(container, E, C, name, t, cost)
+    return cost
+end
+
+"""
+Add cost term to expression and invariant objective.
+
+Computes `cost = quantity * rate`, adds to target expression (if present),
+propagates to `ProductionCostExpression` (when applicable), and adds to the
+time-invariant part of the objective.
+
+# Arguments
+- `container`: the optimization container
+- `quantity`: the value being costed (e.g., variable value, expression value)
+- `rate`: scalar cost rate (e.g., \$/MWh, \$/MMBTU)
+- `E`: target expression type (caller provides, e.g., ProductionCostExpression)
+- `C`: component type
+- `name`: component name
+- `t`: time period
+"""
+function add_cost_term_invariant!(
+    container::OptimizationContainer,
+    quantity::JuMPOrFloat,
+    rate::Float64,
+    ::Type{E},
+    ::Type{C},
+    name::String,
+    t::Int,
+) where {E <: ExpressionType, C <: IS.InfrastructureSystemsComponent}
+    cost = add_cost_term_to_expression!(container, quantity, rate, E, C, name, t)
     add_to_objective_invariant_expression!(container, cost)
     return cost
 end
@@ -76,6 +141,7 @@ function add_cost_term_variant!(
         expr = get_expression(container, E, C)
         JuMP.add_to_expression!(expr[name, t], cost)
     end
+    _propagate_to_production_cost!(container, E, C, name, t, cost)
     add_to_objective_variant_expression!(container, cost)
     return cost
 end
@@ -113,6 +179,7 @@ function add_cost_term_variant!(
         expr = get_expression(container, E, C)
         JuMP.add_to_expression!(expr[name, t], cost)
     end
+    _propagate_to_production_cost!(container, E, C, name, t, cost)
     add_to_objective_variant_expression!(container, cost)
     return cost
 end
@@ -122,7 +189,7 @@ end
 Add a proportional (linear) cost to the invariant objective across all time steps.
 
 Normalizes `cost_term` from `power_units` to system per-unit, multiplies by `dt` and
-`multiplier`, then adds `variable * rate` to `ProductionCostExpression` and the invariant
+`multiplier`, then adds `variable * rate` to the target expression `E` and the invariant
 objective for each time step.
 
 # Arguments
@@ -132,6 +199,9 @@ objective for each time step.
 - `cost_term`: raw proportional cost (e.g., \$/MWh before normalization)
 - `power_units`: unit system of `cost_term`
 - `multiplier`: additional scalar (e.g., `objective_function_multiplier`, fuel cost)
+- `E`: target cost expression type (e.g., `FuelCostExpression`, `VOMCostExpression`).
+  Constituent types auto-propagate into `ProductionCostExpression` via
+  `_propagate_to_production_cost!`.
 """
 function add_proportional_cost_invariant!(
     container::OptimizationContainer,
@@ -139,8 +209,9 @@ function add_proportional_cost_invariant!(
     component::C,
     cost_term::Float64,
     power_units::IS.UnitSystem,
-    multiplier::Float64 = 1.0,
-) where {T <: VariableType, C <: IS.InfrastructureSystemsComponent}
+    multiplier::Float64,
+    ::Type{E},
+) where {T <: VariableType, C <: IS.InfrastructureSystemsComponent, E <: CostExpressions}
     iszero(cost_term) && return
     base_power = get_model_base_power(container)
     device_base_power = get_base_power(component)
@@ -151,8 +222,7 @@ function add_proportional_cost_invariant!(
     rate = cost_per_unit * multiplier * dt
     for t in get_time_steps(container)
         variable = get_variable(container, T, C)[name, t]
-        add_cost_term_invariant!(
-            container, variable, rate, ProductionCostExpression, C, name, t)
+        add_cost_term_invariant!(container, variable, rate, E, C, name, t)
     end
     return
 end
diff --git a/src/objective_function/linear_curve.jl b/src/objective_function/linear_curve.jl
index 0faf1ad..c4881bc 100644
--- a/src/objective_function/linear_curve.jl
+++ b/src/objective_function/linear_curve.jl
@@ -20,7 +20,8 @@ function add_variable_cost_to_objective!(
     proportional_term = get_proportional_term(get_function_data(value_curve))
     multiplier = objective_function_multiplier(T, U)
     add_proportional_cost_invariant!(
-        container, T, component, proportional_term, power_units, multiplier)
+        container, T, component, proportional_term, power_units, multiplier,
+        FuelCostExpression)
     return
 end
 
@@ -59,7 +60,8 @@ function add_variable_cost_to_objective!(
     # Multiplier is not necessary here. There is no negative cost for fuel curves.
     if fuel_cost isa Float64
         add_proportional_cost_invariant!(
-            container, T, component, proportional_term, power_units, fuel_cost)
+            container, T, component, proportional_term, power_units, fuel_cost,
+            FuelCostExpression)
     else
         # Time-varying fuel cost: normalize, then delegate to variant path
         base_power = get_model_base_power(container)
diff --git a/src/objective_function/objective_function_pwl_lambda.jl b/src/objective_function/objective_function_pwl_lambda.jl
index 8a55a59..d852814 100644
--- a/src/objective_function/objective_function_pwl_lambda.jl
+++ b/src/objective_function/objective_function_pwl_lambda.jl
@@ -126,7 +126,9 @@ function _determine_bin_lhs(
     sos_status::SOSStatusVariable,
     ::Type{T},
     name::String,
-    period::Int) where {T <: IS.InfrastructureSystemsComponent}
+    period::Int;
+    on_var_type::Type{<:VariableType} = OnVariable,
+) where {T <: IS.InfrastructureSystemsComponent}
     if sos_status == SOSStatusVariable.NO_VARIABLE
         @debug "Using Piecewise Linear cost function but no variable/parameter ref for ON status is passed. Default status will be set to online (1.0)" _group =
             LOG_GROUP_COST_FUNCTIONS
@@ -136,9 +138,9 @@ function _determine_bin_lhs(
             LOG_GROUP_COST_FUNCTIONS
         return get_parameter(container, OnStatusParameter, T).parameter_array[name, period]
     elseif sos_status == SOSStatusVariable.VARIABLE
-        @debug "Using Piecewise Linear cost function with variable OnVariable $T" _group =
+        @debug "Using Piecewise Linear cost function with variable $on_var_type for $T" _group =
             LOG_GROUP_COST_FUNCTIONS
-        return get_variable(container, OnVariable, T)[name, period]
+        return get_variable(container, on_var_type, T)[name, period]
     else
         @assert false
     end
diff --git a/src/objective_function/piecewise_linear.jl b/src/objective_function/piecewise_linear.jl
index cc1a06a..46fd82f 100644
--- a/src/objective_function/piecewise_linear.jl
+++ b/src/objective_function/piecewise_linear.jl
@@ -178,7 +178,7 @@ function add_variable_cost_to_objective!(
     for t in get_time_steps(container)
         add_cost_to_expression!(
             container,
-            ProductionCostExpression,
+            FuelCostExpression,
             pwl_cost_expressions[t],
             C,
             component_name,
@@ -227,7 +227,7 @@ function add_variable_cost_to_objective!(
         pwl_cost_expression = pwl_fuel_consumption_expressions[t] * fuel_cost_value
         add_cost_to_expression!(
             container,
-            ProductionCostExpression,
+            FuelCostExpression,
             pwl_cost_expression,
             V,
             component_name,
diff --git a/src/objective_function/proportional.jl b/src/objective_function/proportional.jl
index 0abb7ca..7cf2e75 100644
--- a/src/objective_function/proportional.jl
+++ b/src/objective_function/proportional.jl
@@ -29,14 +29,14 @@ function add_proportional_cost!(
         skip = skip_proportional_cost(d)
         for t in get_time_steps(container)
             if skip
-                # must-run etc.: bookkeep in ProductionCostExpression but not in objective
+                # must-run etc.: bookkeep in FixedCostExpression but not in objective
                 add_cost_to_expression!(
-                    container, ProductionCostExpression, rate, T, name, t)
+                    container, FixedCostExpression, rate, T, name, t)
             else
                 variable = get_variable(container, U, T)[name, t]
                 add_cost_term_invariant!(
                     container, variable, rate,
-                    ProductionCostExpression, T, name, t,
+                    FixedCostExpression, T, name, t,
                 )
             end
         end
@@ -75,7 +75,7 @@ function add_proportional_cost_maybe_time_variant!(
                 # Only add to expression, not objective
                 add_cost_to_expression!(
                     container,
-                    ProductionCostExpression,
+                    FixedCostExpression,
                     rate,
                     T,
                     name,
@@ -85,10 +85,10 @@ function add_proportional_cost_maybe_time_variant!(
                 variable = get_variable(container, U, T)[name, t]
                 if add_as_time_variant
                     add_cost_term_variant!(
-                        container, variable, rate, ProductionCostExpression, T, name, t)
+                        container, variable, rate, FixedCostExpression, T, name, t)
                 else
                     add_cost_term_invariant!(
-                        container, variable, rate, ProductionCostExpression, T, name, t)
+                        container, variable, rate, FixedCostExpression, T, name, t)
                 end
             end
         end
diff --git a/src/objective_function/quadratic_curve.jl b/src/objective_function/quadratic_curve.jl
index 65623a4..82f16cb 100644
--- a/src/objective_function/quadratic_curve.jl
+++ b/src/objective_function/quadratic_curve.jl
@@ -12,23 +12,19 @@ function _add_quadraticcurve_variable_term_to_model!(
     name = get_name(component)
     var = get_variable(container, T, V)[name, time_period]
 
-    cost = if quadratic_term_per_unit >= eps()
+    if abs(quadratic_term_per_unit) >= eps()
         @debug "$name Quadratic Variable Cost" _group = LOG_GROUP_COST_FUNCTIONS name
         q_cost =
             (var .^ 2 * quadratic_term_per_unit + var * proportional_term_per_unit) * dt
         add_to_objective_invariant_expression!(container, q_cost)
-        q_cost
+        # add_cost_to_expression! handles ConstituentCostExpression -> ProductionCostExpression
+        # propagation via _propagate_to_production_cost!.
+        add_cost_to_expression!(
+            container, FuelCostExpression, q_cost, V, name, time_period)
     else
         add_cost_term_invariant!(
             container, var, proportional_term_per_unit * dt,
-            ProductionCostExpression, V, name, time_period)
-    end
-
-    # For quadratic case, still need to add to expression (linear case handled by helper)
-    if quadratic_term_per_unit >= eps() &&
-       has_container_key(container, ProductionCostExpression, V)
-        expr = get_expression(container, ProductionCostExpression, V)
-        JuMP.add_to_expression!(expr[name, time_period], cost)
+            FuelCostExpression, V, name, time_period)
     end
     return
 end
diff --git a/src/objective_function/start_up_shut_down.jl b/src/objective_function/start_up_shut_down.jl
index cb98860..e210f00 100644
--- a/src/objective_function/start_up_shut_down.jl
+++ b/src/objective_function/start_up_shut_down.jl
@@ -55,7 +55,7 @@ function _add_shut_down_cost_per_device!(
             rate = cost_term * multiplier
             variable = get_variable(container, U, T)[name, t]
             add_cost_term_variant!(
-                container, variable, rate, ProductionCostExpression, T, name, t)
+                container, variable, rate, ShutDownCostExpression, T, name, t)
         end
     else
         cost_term = _shutdown_cost_value(get_shut_down(op_cost))
@@ -64,7 +64,7 @@ function _add_shut_down_cost_per_device!(
         for t in get_time_steps(container)
             variable = get_variable(container, U, T)[name, t]
             add_cost_term_invariant!(
-                container, variable, rate, ProductionCostExpression, T, name, t)
+                container, variable, rate, ShutDownCostExpression, T, name, t)
         end
     end
 end
@@ -117,7 +117,7 @@ function _add_start_up_cost_to_objective!(
             rate = cost_term * multiplier
             variable = get_variable(container, T, C)[name, t]
             add_cost_term_variant!(
-                container, variable, rate, ProductionCostExpression, C, name, t)
+                container, variable, rate, StartUpCostExpression, C, name, t)
         end
     else
         raw_startup_cost = get_start_up(op_cost)
@@ -127,7 +127,7 @@ function _add_start_up_cost_to_objective!(
             rate = cost_term * multiplier
             variable = get_variable(container, T, C)[name, t]
             add_cost_term_invariant!(
-                container, variable, rate, ProductionCostExpression, C, name, t)
+                container, variable, rate, StartUpCostExpression, C, name, t)
         end
     end
     return
diff --git a/src/utils/indexing.jl b/src/utils/indexing.jl
index 542a7c8..ce8138c 100644
--- a/src/utils/indexing.jl
+++ b/src/utils/indexing.jl
@@ -1,50 +1,62 @@
-# If `ixs` does not index all dimensions of `dest`, add a `:` for the rest (like Python's
-# `...`) to prepare for broadcast-assigning.
-function expand_ixs(ixs::NTuple{1, T}, dest::AbstractArray) where {T}
-    if length(ixs) <= ndims(dest)
-        return (ixs[1], fill(:, ndims(dest) - 1)...)
-    else
-        throw(ArgumentError("`ixs` must not index more dimensions than `dest` has"))
-    end
+# Pad `ixs` with `:` for any unindexed middle dimensions of `dest` (Python `...`-style).
+# Fast path for AbstractArrays where `N` is known at compile time → allocation-free `Val(N - K)`.
+@inline function expand_ixs(
+    ixs::NTuple{K, Any},
+    dest::AbstractArray{<:Any, N},
+) where {K, N}
+    K <= N || throw(ArgumentError("`ixs` must not index more dimensions than `dest` has"))
+    K == N && return ixs
+    # Single-element ixs is the leading axis; multi-element preserves first..last with `:` filling the middle.
+    K == 1 && return (only(ixs), ntuple(_ -> Colon(), Val(N - 1))...)
+    return (Base.front(ixs)..., ntuple(_ -> Colon(), Val(N - K))..., last(ixs))
 end
 
-function expand_ixs(ixs::Tuple{T, U}, dest::AbstractArray) where {T, U}
-    if length(ixs) <= ndims(dest)
-        return (ixs[1], fill(:, ndims(dest) - 2)..., ixs[end])
-    else
-        throw(ArgumentError("`ixs` must not index more dimensions than `dest` has"))
-    end
+# Fallback for non-AbstractArray containers (e.g. `HDF5.Dataset`) — `ndims` resolved at runtime.
+@inline function expand_ixs(ixs::NTuple{K, Any}, dest) where {K}
+    N = ndims(dest)
+    K <= N || throw(ArgumentError("`ixs` must not index more dimensions than `dest` has"))
+    K == N && return ixs
+    K == 1 && return (only(ixs), ntuple(_ -> Colon(), N - 1)...)
+    return (Base.front(ixs)..., ntuple(_ -> Colon(), N - K)..., last(ixs))
 end
 
-function expand_ixs(ixs::Tuple, dest::AbstractArray)
-    if length(ixs) <= ndims(dest)
-        return (ixs[1:(end - 1)]..., fill(:, ndims(dest) - length(ixs))..., ixs[end])
-    else
-        throw(ArgumentError("`ixs` must not index more dimensions than `dest` has"))
-    end
+# Concrete fast path: scalar `src` with a fully-specified index tuple goes through `setindex!`.
+@inline function assign_maybe_broadcast!(
+    dest::DenseAxisArray{T, N},
+    src::T,
+    ixs::NTuple{N, Any},
+) where {T, N}
+    dest[ixs...] = src
+    return
 end
 
-# If `src` is an array, we want to set a slice of `dest` equal to `src`. Broadcast
-# assignment from integer-indexed `src` to DenseAxisArray `dest` slice of same shape doesn't
-# work when `dest`'s axes being broadcast across aren't 1:n, but standard assignment does
-# the trick in that case and (PERF) seems to not appreciably affect simulation performance
-function assign_maybe_broadcast!(dest::AbstractArray, src::AbstractArray, ixs::Tuple)
-    expanded_axs = expand_ixs(ixs, dest)
-    dest[expanded_axs...] = src
+# Array `src`: assign a slice of `dest` from `src` (standard assignment handles non-1:n axes).
+# `dest` is left untyped so non-`AbstractArray` containers (e.g. `HDF5.Dataset`) are also accepted.
+@inline function assign_maybe_broadcast!(dest, src::AbstractArray, ixs::Tuple)
+    dest[expand_ixs(ixs, dest)...] = src
     return
 end
 
-# If `src` is a tuple or scalar, we want to set all values across a slice of `dest` equal to `src`
-function assign_maybe_broadcast!(dest::AbstractArray, src, ixs::Tuple)
-    expanded_axs = expand_ixs(ixs, dest)
-    dest[expanded_axs...] .= Ref(src)
+# Scalar/tuple `src`: broadcast the value across the indexed slice of `dest`.
+@inline function assign_maybe_broadcast!(dest, src, ixs::Tuple)
+    expanded = expand_ixs(ixs, dest)
+    @views dest[expanded...] .= Ref(src)
     return
 end
 
 # Similar to assign_maybe_broadcast! but for fixing JuMP VariableRefs
-fix_expand(dest::AbstractArray, src, ixs::Tuple) =
+@inline function fix_maybe_broadcast!(
+    dest::DenseAxisArray{JuMP.VariableRef, N},
+    src::Float64,
+    ixs::NTuple{N, Any},
+) where {N}
+    JuMP.fix(dest[ixs...], src; force = true)
+    return
+end
+
+fix_expand(dest, src, ixs::Tuple) =
     fix_parameter_value.(dest[expand_ixs(ixs, dest)...], src)
-fix_maybe_broadcast!(dest::AbstractArray, src::AbstractArray, ixs::Tuple) =
+fix_maybe_broadcast!(dest, src::AbstractArray, ixs::Tuple) =
     fix_expand(dest, src, ixs)
-fix_maybe_broadcast!(dest::AbstractArray, src, ixs::Tuple) =
+fix_maybe_broadcast!(dest, src, ixs::Tuple) =
     fix_expand(dest, Ref(src), ixs)
diff --git a/test/InfrastructureOptimizationModelsTests.jl b/test/InfrastructureOptimizationModelsTests.jl
index f218b68..97968b5 100644
--- a/test/InfrastructureOptimizationModelsTests.jl
+++ b/test/InfrastructureOptimizationModelsTests.jl
@@ -106,6 +106,7 @@ function run_tests()
                     include(joinpath(TEST_DIR, "test_optimization_container_metadata.jl"))
                     # optimization_container_types.jl: no need for tests
                     include(joinpath(TEST_DIR, "test_optimization_container.jl"))
+                    include(joinpath(TEST_DIR, "test_external_evaluation.jl"))
                     # optimization_problem_outputs_export.jl: low-complexity
                     include(joinpath(TEST_DIR, "test_optimization_outputs.jl"))
                     include(joinpath(TEST_DIR, "test_optimizer_stats.jl"))
diff --git a/test/test_basic_model_structs.jl b/test/test_basic_model_structs.jl
index 519e6f0..ecbae94 100644
--- a/test/test_basic_model_structs.jl
+++ b/test/test_basic_model_structs.jl
@@ -6,17 +6,25 @@ end
 
 @testset "NetworkModel Tests" begin
     @test_throws ArgumentError NetworkModel(AbstractPowerModel)
+    ec_multi = EvaluationContainer()
+    add_evaluator!(ec_multi, DCPowerFlow, DCPowerFlow())
+    add_evaluator!(
+        ec_multi,
+        PSSEExportPowerFlow,
+        PSSEExportPowerFlow(:v33, "exports"),
+    )
     @test NetworkModel(
         PTDFPowerModel;
         use_slacks = true,
-        power_flow_evaluation = [DCPowerFlow(), PSSEExportPowerFlow(:v33, "exports")],
+        evaluations = ec_multi,
     ) isa NetworkModel
-    @test NetworkModel(
-        PTDFPowerModel;
-        use_slacks = true,
-        power_flow_evaluation = ACPowerFlow(;
-            exporter =
-            PSSEExportPowerFlow(
+
+    ec_ac = EvaluationContainer()
+    add_evaluator!(
+        ec_ac,
+        ACPowerFlow,
+        ACPowerFlow(;
+            exporter = PSSEExportPowerFlow(
                 :v33,
                 "exports";
                 name = "my_export_name",
@@ -24,6 +32,11 @@ end
                 overwrite = true,
             ),
         ),
+    )
+    @test NetworkModel(
+        PTDFPowerModel;
+        use_slacks = true,
+        evaluations = ec_ac,
     ) isa NetworkModel
 end
 
diff --git a/test/test_cost_term_helpers.jl b/test/test_cost_term_helpers.jl
index ccfb828..f4a286f 100644
--- a/test/test_cost_term_helpers.jl
+++ b/test/test_cost_term_helpers.jl
@@ -547,4 +547,241 @@ Test types defined in test_utils/test_types.jl.
             @test JuMP.coefficient(cost_expr, pwl_vars[2]) ≈ 15.0
         end
     end
+
+    # Helpers for the ProductionCostExpression propagation testset
+    function _setup_prop_test_container(time_steps, names = ["gen1"])
+        container = make_test_container(time_steps)
+        var_container = IOM.add_variable_container!(
+            container, TestCostVariable, MockThermalGen, names, time_steps,
+        )
+        jump_model = IOM.get_jump_model(container)
+        for name in names, t in time_steps
+            var_container[name, t] =
+                JuMP.@variable(jump_model, base_name = "v_$(name)_$(t)")
+        end
+        return container, var_container
+    end
+
+    _expr_coef(container, ::Type{E}, var, name, t) where {E} = JuMP.coefficient(
+        IOM.get_expression(container, E, MockThermalGen)[name, t], var,
+    )
+    _inv_coef(container, var) = JuMP.coefficient(
+        IOM.get_invariant_terms(IOM.get_objective_expression(container)), var,
+    )
+    _var_coef(container, var) = JuMP.coefficient(
+        IOM.get_variant_terms(IOM.get_objective_expression(container)), var,
+    )
+
+    @testset "ProductionCostExpression propagation" begin
+        time_steps = 1:1
+
+        @testset "add_cost_term_to_expression! propagates constituent → ProductionCost" begin
+            container, vars = _setup_prop_test_container(time_steps)
+            for E in (IOM.FuelCostExpression, IOM.ProductionCostExpression)
+                add_test_expression!(container, E, MockThermalGen, ["gen1"], time_steps)
+            end
+            r = 11.0
+            IOM.add_cost_term_to_expression!(
+                container, vars["gen1", 1], r,
+                IOM.FuelCostExpression, MockThermalGen, "gen1", 1,
+            )
+            @test _expr_coef(
+                container, IOM.FuelCostExpression, vars["gen1", 1], "gen1", 1) ≈ r
+            @test _expr_coef(
+                container, IOM.ProductionCostExpression, vars["gen1", 1], "gen1", 1) ≈ r
+            # No objective hook on this entry point.
+            @test _inv_coef(container, vars["gen1", 1]) ≈ 0.0
+            @test _var_coef(container, vars["gen1", 1]) ≈ 0.0
+        end
+
+        @testset "add_cost_term_invariant! → constituent, ProductionCost, and invariant obj" begin
+            container, vars = _setup_prop_test_container(time_steps)
+            for E in (IOM.FuelCostExpression, IOM.ProductionCostExpression)
+                add_test_expression!(container, E, MockThermalGen, ["gen1"], time_steps)
+            end
+            r = 13.0
+            IOM.add_cost_term_invariant!(
+                container, vars["gen1", 1], r,
+                IOM.FuelCostExpression, MockThermalGen, "gen1", 1,
+            )
+            @test _expr_coef(
+                container, IOM.FuelCostExpression, vars["gen1", 1], "gen1", 1) ≈ r
+            @test _expr_coef(
+                container, IOM.ProductionCostExpression, vars["gen1", 1], "gen1", 1) ≈ r
+            @test _inv_coef(container, vars["gen1", 1]) ≈ r
+            @test _var_coef(container, vars["gen1", 1]) ≈ 0.0
+        end
+
+        @testset "add_cost_term_invariant! direct ProductionCost write does not recurse" begin
+            container, vars = _setup_prop_test_container(time_steps)
+            add_test_expression!(
+                container, IOM.ProductionCostExpression, MockThermalGen,
+                ["gen1"], time_steps,
+            )
+            r = 17.0
+            IOM.add_cost_term_invariant!(
+                container, vars["gen1", 1], r,
+                IOM.ProductionCostExpression, MockThermalGen, "gen1", 1,
+            )
+            # ProductionCostExpression ⊀ ConstituentCostExpression, so the propagation
+            # hook is a no-op — exactly one write to ProductionCostExpression.
+            @test _expr_coef(
+                container, IOM.ProductionCostExpression, vars["gen1", 1], "gen1", 1) ≈ r
+            @test _inv_coef(container, vars["gen1", 1]) ≈ r
+        end
+
+        @testset "add_cost_term_invariant! to non-constituent reaches obj, not ProductionCost" begin
+            container, vars = _setup_prop_test_container(time_steps)
+            add_test_expression!(
+                container, TestCostExpression, MockThermalGen, ["gen1"], time_steps,
+            )
+            add_test_expression!(
+                container, IOM.ProductionCostExpression, MockThermalGen,
+                ["gen1"], time_steps,
+            )
+            r = 19.0
+            IOM.add_cost_term_invariant!(
+                container, vars["gen1", 1], r,
+                TestCostExpression, MockThermalGen, "gen1", 1,
+            )
+            @test _expr_coef(
+                container, TestCostExpression, vars["gen1", 1], "gen1", 1) ≈ r
+            # TestCostExpression is not a ConstituentCostExpression — ProductionCost untouched.
+            @test _expr_coef(
+                container, IOM.ProductionCostExpression,
+                vars["gen1", 1], "gen1", 1) ≈ 0.0
+            @test _inv_coef(container, vars["gen1", 1]) ≈ r
+        end
+
+        @testset "add_cost_term_variant! (param-rate) propagates and hits variant obj" begin
+            container, vars = _setup_prop_test_container(time_steps)
+            for E in (IOM.FuelCostExpression, IOM.ProductionCostExpression)
+                add_test_expression!(container, E, MockThermalGen, ["gen1"], time_steps)
+            end
+            param_rate = 23.0
+            add_test_parameter!(
+                container, TestCostParameter, MockThermalGen,
+                ["gen1"], time_steps, fill(param_rate, 1, length(time_steps)),
+            )
+            IOM.add_cost_term_variant!(
+                container, vars["gen1", 1], TestCostParameter,
+                IOM.FuelCostExpression, MockThermalGen, "gen1", 1,
+            )
+            @test _expr_coef(
+                container, IOM.FuelCostExpression, vars["gen1", 1], "gen1", 1) ≈ param_rate
+            @test _expr_coef(
+                container, IOM.ProductionCostExpression,
+                vars["gen1", 1], "gen1", 1) ≈ param_rate
+            @test _var_coef(container, vars["gen1", 1]) ≈ param_rate
+            @test _inv_coef(container, vars["gen1", 1]) ≈ 0.0
+        end
+
+        @testset "add_cost_term_variant! (explicit-rate) propagates and hits variant obj" begin
+            container, vars = _setup_prop_test_container(time_steps)
+            for E in (IOM.FuelCostExpression, IOM.ProductionCostExpression)
+                add_test_expression!(container, E, MockThermalGen, ["gen1"], time_steps)
+            end
+            r = 29.0
+            IOM.add_cost_term_variant!(
+                container, vars["gen1", 1], r,
+                IOM.FuelCostExpression, MockThermalGen, "gen1", 1,
+            )
+            @test _expr_coef(
+                container, IOM.FuelCostExpression, vars["gen1", 1], "gen1", 1) ≈ r
+            @test _expr_coef(
+                container, IOM.ProductionCostExpression, vars["gen1", 1], "gen1", 1) ≈ r
+            @test _var_coef(container, vars["gen1", 1]) ≈ r
+            @test _inv_coef(container, vars["gen1", 1]) ≈ 0.0
+        end
+
+        @testset "add_cost_to_expression! propagates constituent → ProductionCost" begin
+            container, vars = _setup_prop_test_container(time_steps)
+            for E in (IOM.FuelCostExpression, IOM.ProductionCostExpression)
+                add_test_expression!(container, E, MockThermalGen, ["gen1"], time_steps)
+            end
+            r = 31.0
+            IOM.add_cost_to_expression!(
+                container, IOM.FuelCostExpression, r * vars["gen1", 1],
+                MockThermalGen, "gen1", 1,
+            )
+            @test _expr_coef(
+                container, IOM.FuelCostExpression, vars["gen1", 1], "gen1", 1) ≈ r
+            @test _expr_coef(
+                container, IOM.ProductionCostExpression, vars["gen1", 1], "gen1", 1) ≈ r
+            # No objective hook on this entry point.
+            @test _inv_coef(container, vars["gen1", 1]) ≈ 0.0
+            @test _var_coef(container, vars["gen1", 1]) ≈ 0.0
+        end
+
+        @testset "add_proportional_cost_invariant! propagates per-time-step" begin
+            ts = 1:3
+            container, vars = _setup_prop_test_container(ts)
+            for E in (IOM.FuelCostExpression, IOM.ProductionCostExpression)
+                add_test_expression!(container, E, MockThermalGen, ["gen1"], ts)
+            end
+            device = make_mock_thermal("gen1"; base_power = 100.0)
+            cost_term = 7.0
+            multiplier = 1.0
+            # SYSTEM_BASE → no normalization; dt = 1 hour → rate = cost_term * multiplier.
+            IOM.add_proportional_cost_invariant!(
+                container, TestCostVariable, device, cost_term,
+                IS.UnitSystem.SYSTEM_BASE, multiplier, IOM.FuelCostExpression,
+            )
+            expected = cost_term * multiplier
+            for t in ts
+                @test _expr_coef(
+                    container, IOM.FuelCostExpression,
+                    vars["gen1", t], "gen1", t) ≈ expected
+                @test _expr_coef(
+                    container, IOM.ProductionCostExpression,
+                    vars["gen1", t], "gen1", t) ≈ expected
+                @test _inv_coef(container, vars["gen1", t]) ≈ expected
+            end
+        end
+
+        @testset "multiple constituents sum into ProductionCost without double-count" begin
+            container, vars = _setup_prop_test_container(time_steps)
+            for E in (IOM.FuelCostExpression, IOM.VOMCostExpression,
+                IOM.StartUpCostExpression, IOM.ProductionCostExpression)
+                add_test_expression!(container, E, MockThermalGen, ["gen1"], time_steps)
+            end
+            r_fuel, r_vom, r_su = 3.0, 5.0, 7.0
+            for (E, r) in ((IOM.FuelCostExpression, r_fuel),
+                (IOM.VOMCostExpression, r_vom),
+                (IOM.StartUpCostExpression, r_su))
+                IOM.add_cost_term_invariant!(
+                    container, vars["gen1", 1], r, E, MockThermalGen, "gen1", 1,
+                )
+            end
+            @test _expr_coef(
+                container, IOM.FuelCostExpression, vars["gen1", 1], "gen1", 1) ≈ r_fuel
+            @test _expr_coef(
+                container, IOM.VOMCostExpression, vars["gen1", 1], "gen1", 1) ≈ r_vom
+            @test _expr_coef(
+                container, IOM.StartUpCostExpression, vars["gen1", 1], "gen1", 1) ≈ r_su
+            # Each constituent contributes once; no double-count, no missing terms.
+            @test _expr_coef(
+                container, IOM.ProductionCostExpression,
+                vars["gen1", 1], "gen1", 1) ≈ r_fuel + r_vom + r_su
+            @test _inv_coef(container, vars["gen1", 1]) ≈ r_fuel + r_vom + r_su
+        end
+
+        @testset "ProductionCost not registered: constituent write is a no-op for prod" begin
+            container, vars = _setup_prop_test_container(time_steps)
+            # Only the constituent container is registered.
+            add_test_expression!(
+                container, IOM.FuelCostExpression, MockThermalGen, ["gen1"], time_steps,
+            )
+            r = 41.0
+            IOM.add_cost_term_invariant!(
+                container, vars["gen1", 1], r,
+                IOM.FuelCostExpression, MockThermalGen, "gen1", 1,
+            )
+            @test _expr_coef(
+                container, IOM.FuelCostExpression, vars["gen1", 1], "gen1", 1) ≈ r
+            @test !IOM.has_container_key(
+                container, IOM.ProductionCostExpression, MockThermalGen)
+            @test _inv_coef(container, vars["gen1", 1]) ≈ r
+        end
+    end
 end
diff --git a/test/test_external_evaluation.jl b/test/test_external_evaluation.jl
new file mode 100644
index 0000000..55a1536
--- /dev/null
+++ b/test/test_external_evaluation.jl
@@ -0,0 +1,132 @@
+import InfrastructureOptimizationModels:
+    AbstractEvaluator,
+    AbstractEvaluationData,
+    EvaluationContainer,
+    get_evaluators,
+    get_evaluation_data,
+    get_evaluator,
+    add_evaluator!,
+    add_evaluation_data!,
+    initialize_evaluation_data,
+    evaluate!,
+    reset!,
+    is_solved,
+    get_inner_data,
+    is_from_evaluator
+
+# Concrete-but-minimal subtypes for interface-stub coverage. No methods registered.
+struct DummyEvaluator <: AbstractEvaluator end
+struct DummyEvaluationData <: AbstractEvaluationData end
+
+# Full mock implementation used to exercise the dispatch path end-to-end.
+struct MockEvaluator <: AbstractEvaluator
+    tag::Symbol
+end
+
+mutable struct MockEvaluationData <: AbstractEvaluationData
+    evaluator::MockEvaluator
+    solved::Bool
+    payload::Vector{Float64}
+end
+
+IOM.initialize_evaluation_data(ev::MockEvaluator, _container, _system) =
+    MockEvaluationData(ev, false, Float64[])
+
+function IOM.evaluate!(data::MockEvaluationData, _container, _system)
+    push!(data.payload, length(data.payload) + 1.0)
+    data.solved = true
+    return
+end
+
+IOM.reset!(data::MockEvaluationData) = (data.solved = false; return)
+IOM.is_solved(data::MockEvaluationData) = data.solved
+IOM.get_inner_data(data::MockEvaluationData) = data.payload
+
+@testset "EvaluationContainer CRUD" begin
+    ec = EvaluationContainer()
+    @test ec isa EvaluationContainer
+    @test isempty(ec)
+    @test length(ec) == 0
+    @test !haskey(ec, MockEvaluator)
+    @test isempty(get_evaluators(ec))
+    @test isempty(get_evaluation_data(ec))
+
+    ev = MockEvaluator(:primary)
+    add_evaluator!(ec, MockEvaluator, ev)
+    @test !isempty(ec)
+    @test length(ec) == 1
+    @test haskey(ec, MockEvaluator)
+    @test get_evaluator(ec, MockEvaluator) === ev
+
+    data = MockEvaluationData(ev, false, Float64[])
+    add_evaluation_data!(ec, MockEvaluator, data)
+    @test get_evaluation_data(ec, MockEvaluator) === data
+    @test length(get_evaluation_data(ec)) == 1
+end
+
+@testset "AbstractEvaluator interface stubs throw with type names" begin
+    container = nothing
+    system = nothing
+    dummy_ev = DummyEvaluator()
+    dummy_data = DummyEvaluationData()
+
+    err1 = @test_throws ErrorException initialize_evaluation_data(
+        dummy_ev,
+        container,
+        system,
+    )
+    @test occursin("DummyEvaluator", err1.value.msg)
+
+    err2 = @test_throws ErrorException evaluate!(dummy_data, container, system)
+    @test occursin("DummyEvaluationData", err2.value.msg)
+
+    err3 = @test_throws ErrorException reset!(dummy_data)
+    @test occursin("DummyEvaluationData", err3.value.msg)
+
+    err4 = @test_throws ErrorException is_solved(dummy_data)
+    @test occursin("DummyEvaluationData", err4.value.msg)
+
+    err5 = @test_throws ErrorException get_inner_data(dummy_data)
+    @test occursin("DummyEvaluationData", err5.value.msg)
+end
+
+@testset "MockEvaluator end-to-end dispatch" begin
+    ec = EvaluationContainer()
+    ev = MockEvaluator(:e1)
+    add_evaluator!(ec, MockEvaluator, ev)
+    data = initialize_evaluation_data(ev, nothing, nothing)
+    add_evaluation_data!(ec, MockEvaluator, data)
+
+    @test !is_solved(data)
+    evaluate!(data, nothing, nothing)
+    @test is_solved(data)
+    @test get_inner_data(data) == [1.0]
+
+    # Iterate via the dict view, as `calculate_aux_variables!` does
+    for (T, d) in get_evaluation_data(ec)
+        reset!(d)
+    end
+    @test !is_solved(data)
+
+    evaluate!(data, nothing, nothing)
+    @test is_solved(data)
+    @test get_inner_data(data) == [1.0, 2.0]
+end
+
+@testset "is_from_evaluator default is false" begin
+    @test is_from_evaluator(MockAuxVariable) === false
+end
+
+@testset "NetworkModel evaluations field defaults to empty" begin
+    nw = IOM.NetworkModel(TestPowerModel)
+    @test IOM.get_evaluations(nw) isa EvaluationContainer
+    @test isempty(IOM.get_evaluations(nw))
+end
+
+@testset "NetworkModel evaluations accepts a populated container" begin
+    ec = EvaluationContainer()
+    add_evaluator!(ec, MockEvaluator, MockEvaluator(:populated))
+    nw = IOM.NetworkModel(TestPowerModel; evaluations = ec)
+    @test IOM.get_evaluations(nw) === ec
+    @test haskey(IOM.get_evaluations(nw), MockEvaluator)
+end
diff --git a/test/test_linear_curve.jl b/test/test_linear_curve.jl
index b83661e..3cb6545 100644
--- a/test/test_linear_curve.jl
+++ b/test/test_linear_curve.jl
@@ -73,6 +73,8 @@ end
             device,
             25.0,
             IS.UnitSystem.NATURAL_UNITS,
+            1.0,
+            IOM.ProductionCostExpression,
         )
 
         expected_coef = 25.0 * 100.0 * 1.0
@@ -104,6 +106,7 @@ end
             20.0,
             IS.UnitSystem.SYSTEM_BASE,
             2.0,
+            IOM.ProductionCostExpression,
         )
 
         expected_coef = 20.0 * 2.0 * 0.25
@@ -128,6 +131,8 @@ end
             device,
             0.0,
             IS.UnitSystem.NATURAL_UNITS,
+            1.0,
+            IOM.ProductionCostExpression,
         )
 
         # Should have zero coefficients since cost_term is zero