From 0ff1f9414b0d9ea532490aafc36983c55dddc7b7 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Tue, 5 May 2026 15:48:27 -0400 Subject: [PATCH 1/2] claude port --- src/PowerOperationsModels.jl | 22 ++ src/common_models/add_expressions.jl | 121 ++++++++- src/common_models/add_to_expression.jl | 244 +++++++++++------- src/common_models/market_bid_plumbing.jl | 3 +- src/common_models/objective_function.jl | 86 ++++++ .../load_constructor.jl | 12 +- .../renewable_generation.jl | 1 + .../renewablegeneration_constructor.jl | 4 +- .../source_constructor.jl | 4 +- .../thermalgeneration_constructor.jl | 48 ++-- ...evice_renewable_generation_constructors.jl | 33 +++ ..._device_thermal_generation_constructors.jl | 120 +++++---- test/test_model_decision.jl | 6 +- 13 files changed, 506 insertions(+), 198 deletions(-) diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl index 27411c1..cf056e8 100644 --- a/src/PowerOperationsModels.jl +++ b/src/PowerOperationsModels.jl @@ -179,6 +179,19 @@ import InfrastructureOptimizationModels: using InfrastructureOptimizationModels # TODO: use explicit imports. +# Cost expression types are imported explicitly so re-exports below are detected as +# defined bindings in POM (Aqua.test_undefined_exports requirement). +import InfrastructureOptimizationModels: + ProductionCostExpression, + FuelConsumptionExpression, + ConstituentCostExpression, + FuelCostExpression, + StartUpCostExpression, + ShutDownCostExpression, + FixedCostExpression, + VOMCostExpression, + CurtailmentCostExpression + # Note: add_feedforward_arguments!, add_feedforward_constraints!, # get_default_on_variable, get_default_off_variable are defined in POM, not IOM # Note: ABSOLUTE_TOLERANCE is defined in POM's definitions.jl @@ -209,6 +222,8 @@ include("core/initial_conditions.jl") include("common_models/add_expressions.jl") # Device-specific add_to_expression! implementations include("common_models/add_to_expression.jl") +# Device-specific objective function helpers (curtailment cost, compact-form guards) +include("common_models/objective_function.jl") # add_param_container.jl: moved into IOM include("common_models/add_parameters.jl") include("common_models/make_system_expressions.jl") @@ -693,6 +708,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 diff --git a/src/common_models/add_expressions.jl b/src/common_models/add_expressions.jl index 507519c..e4793e9 100644 --- a/src/common_models/add_expressions.jl +++ b/src/common_models/add_expressions.jl @@ -14,6 +14,20 @@ end _get_variable_if_exists(::PSY.MarketBidCost) = nothing _get_variable_if_exists(cost::PSY.OperationalCost) = PSY.get_variable(cost) +# Predicates for fuel-curve detection. Dispatch over the value returned by +# `_get_variable_if_exists` so callers can avoid `isa` checks. +_is_fuel_curve(::Nothing) = false +_is_fuel_curve(::PSY.CostCurve) = false +_is_fuel_curve(::PSY.FuelCurve) = true + +# Predicates for value-curve shape. Used to decide whether a FuelConsumptionExpression +# container needs JuMP.QuadExpr storage (vs the cheaper GAE). +_value_curve_is_quadratic(::PSY.LinearCurve) = false +_value_curve_is_quadratic(::PSY.QuadraticCurve) = true +_value_curve_is_quadratic(::PSY.PiecewisePointCurve) = false +_value_curve_is_quadratic(::PSY.IncrementalCurve) = false +_value_curve_is_quadratic(::PSY.AverageRateCurve) = false + function get_reference_bus( model::NetworkModel{T}, b::PSY.ACBus, @@ -61,14 +75,12 @@ function add_expressions!( names = String[] found_quad_fuel_functions = false for d in devices - op_cost = PSY.get_operation_cost(d) - fuel_curve = _get_variable_if_exists(op_cost) - if fuel_curve isa PSY.FuelCurve - push!(names, PSY.get_name(d)) - if !found_quad_fuel_functions - found_quad_fuel_functions = - PSY.get_value_curve(fuel_curve) isa PSY.QuadraticCurve - end + fuel_curve = _get_variable_if_exists(PSY.get_operation_cost(d)) + _is_fuel_curve(fuel_curve) || continue + push!(names, PSY.get_name(d)) + if !found_quad_fuel_functions + found_quad_fuel_functions = + _value_curve_is_quadratic(PSY.get_value_curve(fuel_curve)) end end @@ -84,6 +96,99 @@ function add_expressions!( return end +################################################################################# +# Cost expression bundles +# add_cost_expressions! sets up the full cost-expression containers for a given +# device type/formulation in one call. ThermalGen and RenewableGen overrides add +# their constituent cost expressions; the default falls back to ProductionCostExpression. +################################################################################# + +""" +Default cost-expression setup: register only `ProductionCostExpression`. +""" +function add_cost_expressions!( + container::OptimizationContainer, + devices::U, + model::DeviceModel{D, W}, +) where { + U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, + W <: AbstractDeviceFormulation, +} where {D <: PSY.Component} + time_steps = get_time_steps(container) + names = PSY.get_name.(devices) + add_expression_container!(container, ProductionCostExpression, D, names, time_steps) + return +end + +""" +Thermal generators get the full constituent decomposition. Constituent expressions +auto-propagate into `ProductionCostExpression` (see IOM `_propagate_to_production_cost!`), +so we register the aggregate as well as the parts. `FuelConsumptionExpression` is added +only when at least one device has a `FuelCurve`, mirroring the existing FuelConsumption +specialization. +""" +function add_cost_expressions!( + container::OptimizationContainer, + devices::U, + model::DeviceModel{D, W}, +) where { + U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, + W <: AbstractThermalFormulation, +} where {D <: PSY.ThermalGen} + time_steps = get_time_steps(container) + n = length(devices) + all_names = Vector{String}(undef, n) + fuel_names = sizehint!(String[], n) + has_quad_fuel = false + for (i, d) in enumerate(devices) + name = PSY.get_name(d) + all_names[i] = name + fuel_curve = _get_variable_if_exists(PSY.get_operation_cost(d)) + _is_fuel_curve(fuel_curve) || continue + push!(fuel_names, name) + if !has_quad_fuel + has_quad_fuel = _value_curve_is_quadratic(PSY.get_value_curve(fuel_curve)) + end + end + if !isempty(fuel_names) + expr_type = has_quad_fuel ? JuMP.QuadExpr : GAE + add_expression_container!( + container, FuelConsumptionExpression, D, fuel_names, time_steps; + expr_type = expr_type, + ) + end + add_expression_container!(container, ProductionCostExpression, D, all_names, time_steps) + add_expression_container!(container, FuelCostExpression, D, all_names, time_steps) + add_expression_container!(container, StartUpCostExpression, D, all_names, time_steps) + add_expression_container!(container, ShutDownCostExpression, D, all_names, time_steps) + add_expression_container!(container, FixedCostExpression, D, all_names, time_steps) + add_expression_container!(container, VOMCostExpression, D, all_names, time_steps) + return +end + +""" +Renewable dispatch formulations track production cost, fixed cost, VOM cost, and +curtailment cost. `CurtailmentCostExpression` is a direct `CostExpressions` subtype +(not a `ConstituentCostExpression`), so it does not propagate into +`ProductionCostExpression` — curtailment is reported standalone. +""" +function add_cost_expressions!( + container::OptimizationContainer, + devices::U, + model::DeviceModel{D, W}, +) where { + U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, + W <: AbstractRenewableDispatchFormulation, +} where {D <: PSY.RenewableGen} + time_steps = get_time_steps(container) + names = PSY.get_name.(devices) + add_expression_container!(container, ProductionCostExpression, D, names, time_steps) + add_expression_container!(container, FixedCostExpression, D, names, time_steps) + add_expression_container!(container, CurtailmentCostExpression, D, names, time_steps) + add_expression_container!(container, VOMCostExpression, D, names, time_steps) + return +end + """ Generic implementation for service models with reserves. """ diff --git a/src/common_models/add_to_expression.jl b/src/common_models/add_to_expression.jl index e8e9c7e..da34ef0 100644 --- a/src/common_models/add_to_expression.jl +++ b/src/common_models/add_to_expression.jl @@ -2609,6 +2609,74 @@ function add_cost_to_expression!( return end +# Per-device fuel consumption term builders, dispatched on the value-curve type so the +# decision of how to translate the curve into JuMP terms is a method-resolution problem +# rather than a runtime branch. + +# LinearCurve fuel: linear in the dispatch variable. +function _add_fuel_consumption_term!( + expression, + variable, + name::String, + var_cost::PSY.FuelCurve, + value_curve::PSY.LinearCurve, + base_power::Float64, + device_base_power::Float64, + dt::Float64, + time_steps, +) + power_units = PSY.get_power_units(var_cost) + proportional_term = PSY.get_proportional_term(value_curve) + prop_term_per_unit = get_proportional_cost_per_system_unit( + proportional_term, power_units, base_power, device_base_power) + for t in time_steps + JuMP.add_to_expression!( + expression[name, t], prop_term_per_unit * dt, variable[name, t]) + end + return +end + +# QuadraticCurve fuel: quadratic in the dispatch variable. +function _add_fuel_consumption_term!( + expression, + variable, + name::String, + var_cost::PSY.FuelCurve, + value_curve::PSY.QuadraticCurve, + base_power::Float64, + device_base_power::Float64, + dt::Float64, + time_steps, +) + power_units = PSY.get_power_units(var_cost) + proportional_term = PSY.get_proportional_term(value_curve) + quadratic_term = PSY.get_quadratic_term(value_curve) + prop_term_per_unit = get_proportional_cost_per_system_unit( + proportional_term, power_units, base_power, device_base_power) + quad_term_per_unit = get_quadratic_cost_per_system_unit( + quadratic_term, power_units, base_power, device_base_power) + for t in time_steps + fuel_expr = ( + variable[name, t] .^ 2 * quad_term_per_unit + + variable[name, t] * prop_term_per_unit + ) * dt + JuMP.add_to_expression!(expression[name, t], fuel_expr) + end + return +end + +# Piecewise/incremental/average-rate value curves are populated through their own +# objective paths; no contribution to FuelConsumptionExpression here. +_add_fuel_consumption_term!( + _, _, ::String, ::PSY.FuelCurve, ::PSY.PiecewisePointCurve, + ::Float64, ::Float64, ::Float64, _) = nothing +_add_fuel_consumption_term!( + _, _, ::String, ::PSY.FuelCurve, ::PSY.IncrementalCurve, + ::Float64, ::Float64, ::Float64, _) = nothing +_add_fuel_consumption_term!( + _, _, ::String, ::PSY.FuelCurve, ::PSY.AverageRateCurve, + ::Float64, ::Float64, ::Float64, _) = nothing + function add_to_expression!( container::OptimizationContainer, ::Type{T}, @@ -2627,62 +2695,85 @@ function add_to_expression!( resolution = get_resolution(container) dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR for d in devices - op_cost = PSY.get_operation_cost(d) - var_cost = _get_variable_if_exists(op_cost) - if !(var_cost isa PSY.FuelCurve) - continue - end + var_cost = _get_variable_if_exists(PSY.get_operation_cost(d)) + _is_fuel_curve(var_cost) || continue expression = get_expression(container, T, V) name = PSY.get_name(d) device_base_power = PSY.get_base_power(d) value_curve = PSY.get_value_curve(var_cost) - if value_curve isa PSY.LinearCurve - power_units = PSY.get_power_units(var_cost) - proportional_term = PSY.get_proportional_term(value_curve) - prop_term_per_unit = get_proportional_cost_per_system_unit( - proportional_term, - power_units, - base_power, - device_base_power, - ) - for t in time_steps - JuMP.add_to_expression!( - expression[name, t], - prop_term_per_unit * dt, - variable[name, t], - ) - end - elseif value_curve isa PSY.QuadraticCurve - power_units = PSY.get_power_units(var_cost) - proportional_term = PSY.get_proportional_term(value_curve) - quadratic_term = PSY.get_quadratic_term(value_curve) - prop_term_per_unit = get_proportional_cost_per_system_unit( - proportional_term, - power_units, - base_power, - device_base_power, - ) - quad_term_per_unit = get_quadratic_cost_per_system_unit( - quadratic_term, - power_units, - base_power, - device_base_power, - ) - for t in time_steps - fuel_expr = - ( - variable[name, t] .^ 2 * quad_term_per_unit + - variable[name, t] * prop_term_per_unit - ) * dt - JuMP.add_to_expression!( - expression[name, t], - fuel_expr, - ) - end + _add_fuel_consumption_term!( + expression, variable, name, var_cost, value_curve, + base_power, device_base_power, dt, time_steps, + ) + end +end + +# Compact formulation: dispatch variable is "above-minimum"; constant P_min term is +# added per-time-step, gated by the SOS status (no_variable / parameter / variable). + +function _add_compact_fuel_consumption_term!( + container::OptimizationContainer, + ::Type{W}, + expression, + variable, + d::V, + var_cost::PSY.FuelCurve, + value_curve::PSY.LinearCurve, + base_power::Float64, + device_base_power::Float64, + dt::Float64, + time_steps, +) where {V <: PSY.ThermalGen, W <: AbstractDeviceFormulation} + name = PSY.get_name(d) + P_min = PSY.get_active_power_limits(d).min + power_units = PSY.get_power_units(var_cost) + proportional_term = PSY.get_proportional_term(value_curve) + prop_term_per_unit = get_proportional_cost_per_system_unit( + proportional_term, power_units, base_power, device_base_power) + for t in time_steps + sos_status = _get_sos_value(container, W, d) + if sos_status == SOSStatusVariable.NO_VARIABLE + JuMP.add_to_expression!(expression[name, t], P_min * prop_term_per_unit * dt) + elseif sos_status == SOSStatusVariable.PARAMETER + param = get_default_on_parameter(d) + bin = get_parameter(container, param, V).parameter_array[name, t] + JuMP.add_to_expression!( + expression[name, t], P_min * prop_term_per_unit * dt, bin) + elseif sos_status == SOSStatusVariable.VARIABLE + var = get_default_on_variable(d) + bin = get_variable(container, var, V)[name, t] + JuMP.add_to_expression!( + expression[name, t], P_min * prop_term_per_unit * dt, bin) + else + @assert false end + JuMP.add_to_expression!( + expression[name, t], prop_term_per_unit * dt, variable[name, t]) end + return end +# Compact formulation does not accept QuadraticCurve fuel — the SOS gating breaks down +# for quadratic terms. +_add_compact_fuel_consumption_term!( + ::OptimizationContainer, ::Type{W}, _, _, ::PSY.ThermalGen, ::PSY.FuelCurve, + ::PSY.QuadraticCurve, ::Float64, ::Float64, ::Float64, _, +) where {W <: AbstractDeviceFormulation} = + error("Quadratic Curves are not accepted with Compact Formulation: $W") + +_add_compact_fuel_consumption_term!( + ::OptimizationContainer, ::Type{<:AbstractDeviceFormulation}, + _, _, ::PSY.ThermalGen, ::PSY.FuelCurve, ::PSY.PiecewisePointCurve, + ::Float64, ::Float64, ::Float64, _) = nothing +_add_compact_fuel_consumption_term!( + ::OptimizationContainer, ::Type{<:AbstractDeviceFormulation}, + _, _, ::PSY.ThermalGen, ::PSY.FuelCurve, ::PSY.IncrementalCurve, + ::Float64, ::Float64, ::Float64, _) = nothing +_add_compact_fuel_consumption_term!( + ::OptimizationContainer, ::Type{<:AbstractDeviceFormulation}, + _, _, ::PSY.ThermalGen, ::PSY.FuelCurve, ::PSY.AverageRateCurve, + ::Float64, ::Float64, ::Float64, _) = nothing + function add_to_expression!( container::OptimizationContainer, ::Type{T}, @@ -2701,60 +2792,15 @@ function add_to_expression!( resolution = get_resolution(container) dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR for d in devices - op_cost = PSY.get_operation_cost(d) - var_cost = _get_variable_if_exists(op_cost) - if !(var_cost isa PSY.FuelCurve) - continue - end + var_cost = _get_variable_if_exists(PSY.get_operation_cost(d)) + _is_fuel_curve(var_cost) || continue expression = get_expression(container, T, V) - name = PSY.get_name(d) device_base_power = PSY.get_base_power(d) value_curve = PSY.get_value_curve(var_cost) - P_min = PSY.get_active_power_limits(d).min - if value_curve isa PSY.LinearCurve - power_units = PSY.get_power_units(var_cost) - proportional_term = PSY.get_proportional_term(value_curve) - prop_term_per_unit = get_proportional_cost_per_system_unit( - proportional_term, - power_units, - base_power, - device_base_power, - ) - for t in time_steps - sos_status = _get_sos_value(container, W, d) - if sos_status == SOSStatusVariable.NO_VARIABLE - JuMP.add_to_expression!( - expression[name, t], - P_min * prop_term_per_unit * dt, - ) - elseif sos_status == SOSStatusVariable.PARAMETER - param = get_default_on_parameter(d) - bin = get_parameter(container, param, V).parameter_array[name, t] - JuMP.add_to_expression!( - expression[name, t], - P_min * prop_term_per_unit * dt, - bin, - ) - elseif sos_status == SOSStatusVariable.VARIABLE - var = get_default_on_variable(d) - bin = get_variable(container, var, V)[name, t] - JuMP.add_to_expression!( - expression[name, t], - P_min * prop_term_per_unit * dt, - bin, - ) - else - @assert false - end - JuMP.add_to_expression!( - expression[name, t], - prop_term_per_unit * dt, - variable[name, t], - ) - end - elseif value_curve isa PSY.QuadraticCurve - error("Quadratic Curves are not accepted with Compact Formulation: $W") - end + _add_compact_fuel_consumption_term!( + container, W, expression, variable, d, var_cost, value_curve, + base_power, device_base_power, dt, time_steps, + ) end end diff --git a/src/common_models/market_bid_plumbing.jl b/src/common_models/market_bid_plumbing.jl index 4d9e562..195ff1f 100644 --- a/src/common_models/market_bid_plumbing.jl +++ b/src/common_models/market_bid_plumbing.jl @@ -550,7 +550,8 @@ function _add_vom_cost_to_objective_helper!( ) where {T <: VariableType, U <: AbstractDeviceFormulation} power_units = IS.get_power_units(cost_data) cost_term = IS.get_proportional_term(IS.get_vom_cost(cost_data)) - IOM.add_proportional_cost_invariant!(container, T, component, cost_term, power_units) + IOM.add_proportional_cost_invariant!( + container, T, component, cost_term, power_units, 1.0, VOMCostExpression) return end diff --git a/src/common_models/objective_function.jl b/src/common_models/objective_function.jl index 14b44b0..1e7620f 100644 --- a/src/common_models/objective_function.jl +++ b/src/common_models/objective_function.jl @@ -15,3 +15,89 @@ function add_variable_cost_to_objective!( ) return end + +################################################################################# +# Curtailment cost for renewable generation +# Renewables can specify a `curtailment_cost` field on their operation cost; this +# captures the dollar value of unprovided available power per timestep: +# cost(t) = price * dt * (offer_max(t) - dispatch(t)) +# Routed to `CurtailmentCostExpression` (a direct `CostExpressions` subtype, not a +# `ConstituentCostExpression`), so it is reported standalone and does NOT propagate +# into `ProductionCostExpression`. +################################################################################# + +# Resolve the operating-point upper bound for a renewable at time t. Falls back to the +# device's static max_active_power when no time-series parameter has been registered. +function _renewable_offer_max( + container::OptimizationContainer, + component::C, + name::String, + t::Int, +) where {C <: PSY.RenewableGen} + has_container_key(container, ActivePowerTimeSeriesParameter, C) || return PSY.get_max_active_power(component) + param_container = get_parameter(container, ActivePowerTimeSeriesParameter, C) + multiplier = get_multiplier_array(param_container) + return get_parameter_column_refs(param_container, name)[t] * multiplier[name, t] +end + +# Cost-curve-shape dispatch. LinearCurve is the supported case; other shapes are +# silently skipped (PSI's PR only implemented the LinearCurve branch). +function _add_curtailment_cost!( + container::OptimizationContainer, + ::Type{T}, + component::C, + cost_function::PSY.CostCurve{PSY.LinearCurve}, + ::Type{U}, +) where {T <: VariableType, C <: PSY.RenewableGen, U <: AbstractDeviceFormulation} + base_power = get_model_base_power(container) + device_base_power = PSY.get_base_power(component) + value_curve = PSY.get_value_curve(cost_function) + power_units = PSY.get_power_units(cost_function) + cost_component = PSY.get_function_data(value_curve) + proportional_term = PSY.get_proportional_term(cost_component) + iszero(proportional_term) && return + + proportional_term_per_unit = get_proportional_cost_per_system_unit( + proportional_term, power_units, base_power, device_base_power, + ) + resolution = get_resolution(container) + dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR + name = PSY.get_name(component) + dispatch_vars = get_variable(container, T, C) + + for t in get_time_steps(container) + offer_max = _renewable_offer_max(container, component, name, t) + dispatch = dispatch_vars[name, t] + curtailment_cost = proportional_term_per_unit * dt * (offer_max - dispatch) + add_to_expression!( + container, CurtailmentCostExpression, curtailment_cost, component, t) + end + return +end + +# Other cost-curve shapes don't have a defined curtailment-cost computation here; +# treat as no-op rather than erroring so renewables with PWL/other costs still build. +_add_curtailment_cost!( + ::OptimizationContainer, ::Type{<:VariableType}, ::PSY.RenewableGen, + ::PSY.OperationalCost, ::Type{<:AbstractDeviceFormulation}) = nothing + +""" +Iterate the device set and route each renewable's `curtailment_cost` into +`CurtailmentCostExpression`. Devices whose operation cost has no `curtailment_cost` +field, or where the field is `nothing`, are skipped. +""" +function add_curtailment_cost!( + container::OptimizationContainer, + ::Type{T}, + devices::IS.FlattenIteratorWrapper{V}, + ::Type{U}, +) where {T <: VariableType, V <: PSY.RenewableGen, U <: AbstractDeviceFormulation} + for d in devices + op_cost_data = PSY.get_operation_cost(d) + hasproperty(op_cost_data, :curtailment_cost) || continue + cost_function = PSY.get_curtailment_cost(op_cost_data) + isnothing(cost_function) && continue + _add_curtailment_cost!(container, T, d, cost_function, U) + end + return +end diff --git a/src/static_injector_models/load_constructor.jl b/src/static_injector_models/load_constructor.jl index 56b5dfc..79076db 100644 --- a/src/static_injector_models/load_constructor.jl +++ b/src/static_injector_models/load_constructor.jl @@ -42,7 +42,7 @@ function construct_device!( add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model) end - add_expressions!(container, ProductionCostExpression, devices, model) + add_cost_expressions!(container, devices, model) add_event_arguments!(container, devices, model, network_model) return end @@ -122,7 +122,7 @@ function construct_device!( process_market_bid_parameters!(container, devices, model, false, true) - add_expressions!(container, ProductionCostExpression, devices, model) + add_cost_expressions!(container, devices, model) add_event_arguments!(container, devices, model, network_model) return end @@ -202,7 +202,7 @@ function construct_device!( process_market_bid_parameters!(container, devices, model, false, true) - add_expressions!(container, ProductionCostExpression, devices, model) + add_cost_expressions!(container, devices, model) add_event_arguments!(container, devices, model, network_model) return end @@ -287,7 +287,7 @@ function construct_device!( add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model) end - add_expressions!(container, ProductionCostExpression, devices, model) + add_cost_expressions!(container, devices, model) add_event_arguments!(container, devices, model, network_model) return end @@ -562,7 +562,7 @@ function construct_device!( network_model, ) - add_expressions!(container, ProductionCostExpression, devices, model) + add_cost_expressions!(container, devices, model) add_event_arguments!(container, devices, model, network_model) return end @@ -672,7 +672,7 @@ function construct_device!( network_model, ) - add_expressions!(container, ProductionCostExpression, devices, model) + add_cost_expressions!(container, devices, model) add_event_arguments!(container, devices, model, network_model) return end diff --git a/src/static_injector_models/renewable_generation.jl b/src/static_injector_models/renewable_generation.jl index c6d2df3..a2e3a6c 100644 --- a/src/static_injector_models/renewable_generation.jl +++ b/src/static_injector_models/renewable_generation.jl @@ -158,5 +158,6 @@ function add_to_objective_function!( ::Type{<:AbstractPowerModel}, ) where {T <: PSY.RenewableGen, U <: AbstractRenewableDispatchFormulation} add_variable_cost!(container, ActivePowerVariable, devices, U) + add_curtailment_cost!(container, ActivePowerVariable, devices, U) return end diff --git a/src/static_injector_models/renewablegeneration_constructor.jl b/src/static_injector_models/renewablegeneration_constructor.jl index 9c42d9a..39ff93c 100644 --- a/src/static_injector_models/renewablegeneration_constructor.jl +++ b/src/static_injector_models/renewablegeneration_constructor.jl @@ -20,7 +20,7 @@ function construct_device!( end process_market_bid_parameters!(container, devices, model) - add_expressions!(container, ProductionCostExpression, devices, model) + add_cost_expressions!(container, devices, model) # Expression add_to_expression!( @@ -139,7 +139,7 @@ function construct_device!( end process_market_bid_parameters!(container, devices, model) - add_expressions!(container, ProductionCostExpression, devices, model) + add_cost_expressions!(container, devices, model) # Expression add_to_expression!( diff --git a/src/static_injector_models/source_constructor.jl b/src/static_injector_models/source_constructor.jl index 7ddd58d..ddc1e04 100644 --- a/src/static_injector_models/source_constructor.jl +++ b/src/static_injector_models/source_constructor.jl @@ -56,7 +56,7 @@ function construct_device!( network_model, ) - add_expressions!(container, ProductionCostExpression, devices, model) + add_cost_expressions!(container, devices, model) add_to_expression!( container, @@ -225,7 +225,7 @@ function construct_device!( model, ) - add_expressions!(container, ProductionCostExpression, devices, model) + add_cost_expressions!(container, devices, model) add_to_expression!( container, diff --git a/src/static_injector_models/thermalgeneration_constructor.jl b/src/static_injector_models/thermalgeneration_constructor.jl index e495f43..5211b79 100644 --- a/src/static_injector_models/thermalgeneration_constructor.jl +++ b/src/static_injector_models/thermalgeneration_constructor.jl @@ -102,8 +102,7 @@ function construct_device!( network_model, ) - add_expressions!(container, FuelConsumptionExpression, devices, device_model) - add_expressions!(container, ProductionCostExpression, devices, device_model) + add_cost_expressions!(container, devices, device_model) add_to_expression!( container, @@ -253,8 +252,7 @@ function construct_device!( network_model, ) - add_expressions!(container, FuelConsumptionExpression, devices, device_model) - add_expressions!(container, ProductionCostExpression, devices, device_model) + add_cost_expressions!(container, devices, device_model) add_to_expression!( container, @@ -389,8 +387,7 @@ function construct_device!( network_model, ) - add_expressions!(container, FuelConsumptionExpression, devices, device_model) - add_expressions!(container, ProductionCostExpression, devices, device_model) + add_cost_expressions!(container, devices, device_model) add_to_expression!( container, @@ -534,8 +531,7 @@ function construct_device!( network_model, ) - add_expressions!(container, FuelConsumptionExpression, devices, device_model) - add_expressions!(container, ProductionCostExpression, devices, device_model) + add_cost_expressions!(container, devices, device_model) add_to_expression!( container, @@ -670,8 +666,7 @@ function construct_device!( network_model, ) - add_expressions!(container, FuelConsumptionExpression, devices, device_model) - add_expressions!(container, ProductionCostExpression, devices, device_model) + add_cost_expressions!(container, devices, device_model) add_to_expression!( container, @@ -795,8 +790,7 @@ function construct_device!( network_model, ) - add_expressions!(container, FuelConsumptionExpression, devices, device_model) - add_expressions!(container, ProductionCostExpression, devices, device_model) + add_cost_expressions!(container, devices, device_model) add_to_expression!( container, @@ -919,8 +913,7 @@ function construct_device!( network_model, ) - add_expressions!(container, FuelConsumptionExpression, devices, device_model) - add_expressions!(container, ProductionCostExpression, devices, device_model) + add_cost_expressions!(container, devices, device_model) add_to_expression!( container, @@ -1028,8 +1021,7 @@ function construct_device!( network_model, ) - add_expressions!(container, FuelConsumptionExpression, devices, device_model) - add_expressions!(container, ProductionCostExpression, devices, device_model) + add_cost_expressions!(container, devices, device_model) add_to_expression!( container, @@ -1168,8 +1160,7 @@ function construct_device!( network_model, ) - add_expressions!(container, FuelConsumptionExpression, devices, device_model) - add_expressions!(container, ProductionCostExpression, devices, device_model) + add_cost_expressions!(container, devices, device_model) add_to_expression!( container, @@ -1348,8 +1339,7 @@ function construct_device!( network_model, ) - add_expressions!(container, FuelConsumptionExpression, devices, device_model) - add_expressions!(container, ProductionCostExpression, devices, device_model) + add_cost_expressions!(container, devices, device_model) add_to_expression!( container, @@ -1526,8 +1516,7 @@ function construct_device!( network_model, ) - add_expressions!(container, FuelConsumptionExpression, devices, device_model) - add_expressions!(container, ProductionCostExpression, devices, device_model) + add_cost_expressions!(container, devices, device_model) add_to_expression!( container, @@ -1681,8 +1670,7 @@ function construct_device!( network_model, ) - add_expressions!(container, FuelConsumptionExpression, devices, device_model) - add_expressions!(container, ProductionCostExpression, devices, device_model) + add_cost_expressions!(container, devices, device_model) add_to_expression!( container, @@ -1834,8 +1822,7 @@ function construct_device!( network_model, ) - add_expressions!(container, FuelConsumptionExpression, devices, device_model) - add_expressions!(container, ProductionCostExpression, devices, device_model) + add_cost_expressions!(container, devices, device_model) add_to_expression!( container, @@ -1986,8 +1973,7 @@ function construct_device!( network_model, ) - add_expressions!(container, FuelConsumptionExpression, devices, device_model) - add_expressions!(container, ProductionCostExpression, devices, device_model) + add_cost_expressions!(container, devices, device_model) add_to_expression!( container, @@ -2113,8 +2099,7 @@ function construct_device!( network_model, ) - add_expressions!(container, FuelConsumptionExpression, devices, device_model) - add_expressions!(container, ProductionCostExpression, devices, device_model) + add_cost_expressions!(container, devices, device_model) add_to_expression!( container, @@ -2261,8 +2246,7 @@ function construct_device!( initial_conditions!(container, devices, ThermalCompactDispatch()) - add_expressions!(container, ProductionCostExpression, devices, device_model) - add_expressions!(container, FuelConsumptionExpression, devices, device_model) + add_cost_expressions!(container, devices, device_model) add_to_expression!( container, diff --git a/test/test_device_renewable_generation_constructors.jl b/test/test_device_renewable_generation_constructors.jl index 7c7add1..9a10e49 100644 --- a/test/test_device_renewable_generation_constructors.jl +++ b/test/test_device_renewable_generation_constructors.jl @@ -75,3 +75,36 @@ end mock_construct_device!(model, device_model; add_event_model = true) moi_tests(model, 0, 0, 0, 0, 0, false) =# end + +@testset "Test Renewable CurtailmentCostExpression nonnegativity" begin + c_sys5_re = PSB.build_system(PSITestSystems, "c_sys5_re") + + template = OperationsProblemTemplate(NetworkModel(CopperPlatePowerModel)) + set_device_model!(template, RenewableDispatch, RenewableFullDispatch) + set_device_model!(template, ThermalStandard, ThermalStandardDispatch) + set_device_model!(template, PowerLoad, StaticPowerLoad) + + model = DecisionModel( + template, + c_sys5_re; + name = "RE_curtailment_cost", + optimizer = HiGHS_optimizer, + optimizer_solve_log_print = true, + ) + + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + outputs = OptimizationProblemOutputs(model) + expr_curt = read_expression( + outputs, + "CurtailmentCostExpression__RenewableDispatch"; + table_format = TableFormat.WIDE, + ) + + tol = 1e-8 + for unit in names(expr_curt)[2:end] + @test all(expr_curt[!, unit] .>= -tol) + end +end diff --git a/test/test_device_thermal_generation_constructors.jl b/test/test_device_thermal_generation_constructors.jl index c784088..bf82664 100644 --- a/test/test_device_thermal_generation_constructors.jl +++ b/test/test_device_thermal_generation_constructors.jl @@ -34,83 +34,111 @@ const TIME1 = DateTime("2024-01-01T00:00:00") "ProductionCostExpression__ThermalStandard"; table_format = TableFormat.WIDE, ) - var_unit_cost = sum(expr[!, "Test Unit"]) + unit = "Test Unit" + var_unit_cost = sum(expr[!, unit]) @test isapprox(var_unit_cost, cost_reference; atol = 1) if thermal_formulation == ThermalBasicUnitCommitment # Tests shut down cost - @test expr[!, "Test Unit"][end] == 0.75 + @test expr[!, unit][end] == 0.75 + + # Decomposition: production == fuel + startup + shutdown + fixed + VOM + expr_fuel = read_expression( + outputs, "FuelCostExpression__ThermalStandard"; + table_format = TableFormat.WIDE, + ) + expr_su = read_expression( + outputs, "StartUpCostExpression__ThermalStandard"; + table_format = TableFormat.WIDE, + ) + expr_sd = read_expression( + outputs, "ShutDownCostExpression__ThermalStandard"; + table_format = TableFormat.WIDE, + ) + expr_fixed = read_expression( + outputs, "FixedCostExpression__ThermalStandard"; + table_format = TableFormat.WIDE, + ) + expr_VOM = read_expression( + outputs, "VOMCostExpression__ThermalStandard"; + table_format = TableFormat.WIDE, + ) + decomp_vec = + expr_fuel[!, unit] .+ expr_su[!, unit] .+ expr_sd[!, unit] .+ + expr_fixed[!, unit] .+ expr_VOM[!, unit] + @test all(isapprox.(decomp_vec, expr[!, unit]; atol = 1e-6)) + @test isapprox(sum(expr[!, unit]), sum(decomp_vec); atol = 1e-6) + + # Nonnegativity (tolerate tiny numerical negatives) + tol = 1e-8 + @test all(expr_fuel[!, unit] .>= -tol) + @test all(expr_su[!, unit] .>= -tol) + @test all(expr_sd[!, unit] .>= -tol) + @test all(expr_fixed[!, unit] .>= -tol) + @test all(expr_VOM[!, unit] .>= -tol) else - @test expr[!, "Test Unit"][end] == 0.0 + @test expr[!, unit][end] == 0.0 end end end - @testset "Test startup cost tracking - compare with and without startup" begin - # Run 1: Normal case (units already ON, no startup) - sys_no_startup = build_system(PSITestSystems, "c_linear_cost_test") - + @testset "Test startup cost tracking" begin template = OperationsProblemTemplate(NetworkModel(CopperPlatePowerModel)) set_device_model!(template, ThermalStandard, ThermalBasicUnitCommitment) set_device_model!(template, PowerLoad, StaticPowerLoad) + unit = "Test Unit" - model_no_startup = DecisionModel( - template, - sys_no_startup; + # Run 1: units initially ON — no startup expected + sys_no_startup = build_system(PSITestSystems, "c_linear_cost_test") + model_no = DecisionModel( + template, sys_no_startup; name = "UC_no_startup", optimizer = HiGHS_optimizer, optimizer_solve_log_print = true, ) - @test build!(model_no_startup; output_dir = test_path) == IOM.ModelBuildStatus.BUILT - @test solve!(model_no_startup) == IOM.RunStatus.SUCCESSFULLY_FINALIZED - - outputs_no_startup = OptimizationProblemOutputs(model_no_startup) - expr_no_startup = read_expression( - outputs_no_startup, - "ProductionCostExpression__ThermalStandard"; + @test build!(model_no; output_dir = test_path) == IOM.ModelBuildStatus.BUILT + @test solve!(model_no) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + out_no = OptimizationProblemOutputs(model_no) + prod_no = read_expression( + out_no, "ProductionCostExpression__ThermalStandard"; table_format = TableFormat.WIDE, ) - cost_no_startup = expr_no_startup[1, "Test Unit"] - - # Run 2: With startup (units initially OFF) - sys_with_startup = build_system(PSITestSystems, "c_linear_cost_test") + cost_no_t1 = prod_no[1, unit] - # Set thermal units to initially be OFF to force a startup - thermal_units = collect(get_components(ThermalStandard, sys_with_startup)) - for unit in thermal_units - set_status!(unit, false) # Set unit to OFF initially - set_time_at_status!(unit, 10.0) # Set time at OFF status > min down time + # Run 2: units initially OFF — startup forced in first timestep + sys_yes = build_system(PSITestSystems, "c_linear_cost_test") + for u in collect(get_components(ThermalStandard, sys_yes)) + set_status!(u, false) + set_time_at_status!(u, 10.0) # > min down time end - - model_with_startup = DecisionModel( - template, - sys_with_startup; + model_yes = DecisionModel( + template, sys_yes; name = "UC_with_startup", optimizer = HiGHS_optimizer, optimizer_solve_log_print = true, ) - @test build!(model_with_startup; output_dir = test_path) == - IOM.ModelBuildStatus.BUILT - @test solve!(model_with_startup) == IOM.RunStatus.SUCCESSFULLY_FINALIZED - - outputs_with_startup = OptimizationProblemOutputs(model_with_startup) - expr_with_startup = read_expression( - outputs_with_startup, - "ProductionCostExpression__ThermalStandard"; + @test build!(model_yes; output_dir = test_path) == IOM.ModelBuildStatus.BUILT + @test solve!(model_yes) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + out_yes = OptimizationProblemOutputs(model_yes) + prod_yes = read_expression( + out_yes, "ProductionCostExpression__ThermalStandard"; table_format = TableFormat.WIDE, ) - cost_with_startup = expr_with_startup[1, "Test Unit"] + cost_yes_t1 = prod_yes[1, unit] - # Verify startup actually occurred start_vars = read_variable( - outputs_with_startup, - "StartVariable__ThermalStandard"; + out_yes, "StartVariable__ThermalStandard"; table_format = TableFormat.WIDE, ) - @test start_vars[1, "Test Unit"] > 0.5 # Startup in first timestep + @test start_vars[1, unit] > 0.5 - # The key test: cost with startup should be greater than without startup - # because it includes the startup cost in addition to generation costs - @test cost_with_startup > cost_no_startup + expr_su = read_expression( + out_yes, "StartUpCostExpression__ThermalStandard"; + table_format = TableFormat.WIDE, + ) + startup_cost_t1 = expr_su[1, unit] + @test startup_cost_t1 > 0.0 + @test cost_yes_t1 > cost_no_t1 + @test isapprox(cost_yes_t1 - cost_no_t1, startup_cost_t1; atol = 1e-6) end end diff --git a/test/test_model_decision.jl b/test/test_model_decision.jl index 0e4ccd2..aa67785 100644 --- a/test/test_model_decision.jl +++ b/test/test_model_decision.jl @@ -68,7 +68,7 @@ end @test length(read_variables(res; table_format = TableFormat.WIDE)) == 4 @test length(read_parameters(res; table_format = TableFormat.WIDE)) == 1 @test length(read_duals(res; table_format = TableFormat.WIDE)) == 0 - @test length(read_expressions(res; table_format = TableFormat.WIDE)) == 2 + @test length(read_expressions(res; table_format = TableFormat.WIDE)) == 7 @test read_variables( res, ["StartVariable__ThermalStandard"]; @@ -416,7 +416,9 @@ end # Manually Multiply by the base power var1_a has natural units and export writes directly from the solver @test var1_a.value == var4.value .* 100.0 - @test length(readdir(IOM.export_realized_outputs(outputs1))) === 7 + exported = readdir(IOM.export_realized_outputs(outputs1)) + @test length(exported) >= 12 + @test any(contains.(exported, "ProductionCostExpression")) end @testset "Test Numerical Stability of Constraints" begin From 70a49677404b8a3ba9335790a1e9d720eaa48224 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Tue, 5 May 2026 17:03:50 -0400 Subject: [PATCH 2/2] formatting --- src/common_models/add_to_expression.jl | 9 +++++---- src/common_models/objective_function.jl | 3 ++- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/common_models/add_to_expression.jl b/src/common_models/add_to_expression.jl index da34ef0..07c14c3 100644 --- a/src/common_models/add_to_expression.jl +++ b/src/common_models/add_to_expression.jl @@ -2656,10 +2656,11 @@ function _add_fuel_consumption_term!( quad_term_per_unit = get_quadratic_cost_per_system_unit( quadratic_term, power_units, base_power, device_base_power) for t in time_steps - fuel_expr = ( - variable[name, t] .^ 2 * quad_term_per_unit + - variable[name, t] * prop_term_per_unit - ) * dt + fuel_expr = + ( + variable[name, t] .^ 2 * quad_term_per_unit + + variable[name, t] * prop_term_per_unit + ) * dt JuMP.add_to_expression!(expression[name, t], fuel_expr) end return diff --git a/src/common_models/objective_function.jl b/src/common_models/objective_function.jl index 1e7620f..94d9ad1 100644 --- a/src/common_models/objective_function.jl +++ b/src/common_models/objective_function.jl @@ -34,7 +34,8 @@ function _renewable_offer_max( name::String, t::Int, ) where {C <: PSY.RenewableGen} - has_container_key(container, ActivePowerTimeSeriesParameter, C) || return PSY.get_max_active_power(component) + has_container_key(container, ActivePowerTimeSeriesParameter, C) || + return PSY.get_max_active_power(component) param_container = get_parameter(container, ActivePowerTimeSeriesParameter, C) multiplier = get_multiplier_array(param_container) return get_parameter_column_refs(param_container, name)[t] * multiplier[name, t]