diff --git a/src/core/optimization_container.jl b/src/core/optimization_container.jl index a90fe12..0eeee71 100644 --- a/src/core/optimization_container.jl +++ b/src/core/optimization_container.jl @@ -480,35 +480,13 @@ function compute_conflict!(container::OptimizationContainer) JuMP.unset_silent(jump_model) jump_model.is_model_dirty = false conflict = container.infeasibility_conflict + # Only the `JuMP.compute_conflict!` call itself signals whether the optimizer + # supports IIS/conflict refining. Keep its `try` narrow so that a failure to + # *read back* the conflict of a single constraint container (below) is not + # misreported as "optimizer doesn't support IIS" and does not discard the + # conflict that was successfully computed. try JuMP.compute_conflict!(jump_model) - conflict_status = MOI.get(jump_model, MOI.ConflictStatus()) - if conflict_status != MOI.CONFLICT_FOUND - @error "No conflict could be found for the model. Status: $conflict_status" - if !get_optimizer_solve_log_print(settings) - JuMP.set_silent(jump_model) - end - return conflict_status - end - - for (key, field_container) in get_constraints(container) - conflict_indices = check_conflict_status(jump_model, field_container) - if isempty(conflict_indices) - @info "Conflict Index returned empty for $key" - continue - else - conflict[encode_key(key)] = conflict_indices - end - end - - msg = IOBuffer() - for (k, v) in conflict - PrettyTables.pretty_table(msg, v; header = [k]) - end - - @error "Constraints participating in conflict basis (IIS) \n\n$(String(take!(msg)))" - - return conflict_status catch e jump_model.is_model_dirty = true if isa(e, MethodError) @@ -516,9 +494,46 @@ function compute_conflict!(container::OptimizationContainer) else @error "Can't compute conflict" exception = (e, catch_backtrace()) end + return MOI.NO_CONFLICT_EXISTS + end + + conflict_status = MOI.get(jump_model, MOI.ConflictStatus()) + if conflict_status != MOI.CONFLICT_FOUND + @error "No conflict could be found for the model. Status: $conflict_status" + if !get_optimizer_solve_log_print(settings) + JuMP.set_silent(jump_model) + end + return conflict_status + end + + # Label each constraint container independently: a failure to read the + # conflict status of one container (e.g. an unsupported constraint type) + # must not abort the loop and hide the conflict found for the others. + for (key, field_container) in get_constraints(container) + conflict_indices = try + check_conflict_status(jump_model, field_container) + catch e + @warn "Could not read conflict status for $key; skipping" exception = + (e, catch_backtrace()) + continue + end + if isempty(conflict_indices) + @info "Conflict Index returned empty for $key" + else + conflict[encode_key(key)] = conflict_indices + end end - return MOI.NO_CONFLICT_EXISTS + msg = IOBuffer() + for (k, v) in conflict + PrettyTables.pretty_table(msg, v; header = [k]) + end + @error "Constraints participating in conflict basis (IIS) \n\n$(String(take!(msg)))" + + if !get_optimizer_solve_log_print(settings) + JuMP.set_silent(jump_model) + end + return conflict_status end function write_optimizer_stats!(container::OptimizationContainer) diff --git a/src/objective_function/objective_function_pwl_delta.jl b/src/objective_function/objective_function_pwl_delta.jl index 97fc649..8bce0de 100644 --- a/src/objective_function/objective_function_pwl_delta.jl +++ b/src/objective_function/objective_function_pwl_delta.jl @@ -153,16 +153,18 @@ end # POM overrides these for specific device types and formulations. ################################################## +# Fallbacks accept any formulation (device or service): service formulations +# (e.g. StepwiseCostReserve) carry no min-gen-power offset, so they resolve here. _include_min_gen_power_in_constraint( ::Type, ::Type{<:VariableType}, - ::Type{<:AbstractDeviceFormulation}, + ::Type, ) = false _include_constant_min_gen_power_in_constraint( ::Type, ::Type{<:VariableType}, - ::Type{<:AbstractDeviceFormulation}, + ::Type, ) = false ################################################## @@ -186,16 +188,18 @@ function add_pwl_constraint_delta!( break_points::Vector{<:JuMPOrFloat}, pwl_vars::Vector{JuMP.VariableRef}, period::Int, - ::Type{W}, + ::Type{W}; + meta = CONTAINER_KEY_EMPTY_META, ) where {T <: IS.InfrastructureSystemsComponent, U <: VariableType, - D <: AbstractDeviceFormulation, + D, W <: AbstractPiecewiseLinearBlockOfferConstraint} - variables = get_variable(container, U, T) + variables = get_variable(container, U, T, meta) const_container = lazy_container_addition!( container, W, T, - axes(variables)..., + axes(variables)...; + meta = meta, ) name = get_name(component) diff --git a/src/objective_function/value_curve_cost.jl b/src/objective_function/value_curve_cost.jl index c9a5515..9dab945 100644 --- a/src/objective_function/value_curve_cost.jl +++ b/src/objective_function/value_curve_cost.jl @@ -202,11 +202,12 @@ function _get_raw_pwl_data( ::Type{T}, name::String, cost_data::IS.CostCurve{IS.TimeSeriesPiecewiseIncrementalCurve}, - time::Int, + time::Int; + meta = CONTAINER_KEY_EMPTY_META, ) where {T <: IS.InfrastructureSystemsComponent} SlopeParam = _slope_param(dir) - slope_arr = get_parameter_array(container, SlopeParam, T) - slope_mult = get_parameter_multiplier_array(container, SlopeParam, T) + slope_arr = get_parameter_array(container, SlopeParam, T, meta) + slope_mult = get_parameter_multiplier_array(container, SlopeParam, T, meta) @assert size(slope_arr) == size(slope_mult) seg_axis = axes(slope_arr)[2] slope_cost_component = Vector{Float64}(undef, length(seg_axis)) @@ -215,8 +216,8 @@ function _get_raw_pwl_data( end BreakpointParam = _breakpoint_param(dir) - bp_arr = get_parameter_array(container, BreakpointParam, T) - bp_mult = get_parameter_multiplier_array(container, BreakpointParam, T) + bp_arr = get_parameter_array(container, BreakpointParam, T, meta) + bp_mult = get_parameter_multiplier_array(container, BreakpointParam, T, meta) @assert size(bp_arr) == size(bp_mult) point_axis = axes(bp_arr)[2] breakpoint_cost_component = Vector{Float64}(undef, length(point_axis)) @@ -355,3 +356,79 @@ function _fill_pwl_data_from_arrays!( copyto!(breakpoints, converted_bp) return end + +################################################################################# +# Section 6: Tranche-axis helpers for time-varying PWL parameters (PSY-free) +# +# Time-varying piecewise-linear cost curves are stored in 3-D parameter arrays +# `(component, tranche, time)`. These helpers size the tranche axis to the global +# maximum across all time steps and unwrap each time-series element (an +# `IS.PiecewiseStepData`) into the per-tranche slope/breakpoint vector that fills +# one column of the array. The PSY-touching orchestration (resolving a service's +# time series into raw data, `calc_additional_axes`) lives in the downstream +# package (POM) and calls down into these data-level helpers. +################################################################################# + +# It's nice for debugging to have meaningful labels on the tranche axis. These +# labels are never relied upon numerically. +make_tranche_axis(n_tranches) = "tranche_" .* string.(1:n_tranches) + +""" +Given a parameter array, return any additional axes, i.e. those that aren't the +first (component) or the last (time). +""" +lookup_additional_axes(parameter_array) = axes(parameter_array)[2:(end - 1)] + +# Maximum number of tranches (segments) across a piecewise time series. +get_max_tranches(data::Vector{IS.PiecewiseStepData}) = maximum(length, data) +get_max_tranches(data::TS.TimeArray) = get_max_tranches(TS.values(data)) +get_max_tranches(data::AbstractDict) = maximum(get_max_tranches.(values(data))) + +# Layer of indirection so that parameters whose time series represents multiple +# things (e.g. both slopes and breakpoints come from the same `PiecewiseStepData` +# series) can unwrap each element into the per-tranche vector that fills the +# parameter array. The default is the identity used by scalar time-series +# parameters. +unwrap_for_param(::ParameterType, ts_elem, expected_axs) = ts_elem + +# For piecewise data the number of tranches can vary over time, so the parameter +# container is sized for the maximum number of tranches and shorter curves are +# padded. We pad with "degenerate" tranches at the top end of the curve with +# dx = 0 so their dispatch variables are constrained to 0. The slope of those +# segments shouldn't matter; we use slope = 0 so the term drops trivially from +# the objective. +function unwrap_for_param( + ::AbstractPiecewiseLinearSlopeParameter, + ts_elem::IS.PiecewiseStepData, + expected_axs::Tuple{AbstractVector}, +) + max_len = length(only(expected_axs)) + y_coords = IS.get_y_coords(ts_elem) + length(y_coords) <= max_len || error( + "PiecewiseStepData y-coords ($(length(y_coords))) exceed expected axis length " * + "($max_len) for slope parameter", + ) + fill_value = 0.0 # pad with slope = 0 if necessary (see above) + out = Vector{Float64}(undef, max_len) + copyto!(out, y_coords) + fill!(view(out, (length(y_coords) + 1):max_len), fill_value) + return out +end + +function unwrap_for_param( + ::AbstractPiecewiseLinearBreakpointParameter, + ts_elem::IS.PiecewiseStepData, + expected_axs::Tuple{AbstractVector}, +) + max_len = length(only(expected_axs)) + x_coords = IS.get_x_coords(ts_elem) + length(x_coords) <= max_len || error( + "PiecewiseStepData x-coords ($(length(x_coords))) exceed expected axis length " * + "($max_len) for breakpoint parameter", + ) + fill_value = x_coords[end] # repeat the last breakpoint so dx = 0 (see above) + out = Vector{Float64}(undef, max_len) + copyto!(out, x_coords) + fill!(view(out, (length(x_coords) + 1):max_len), fill_value) + return out +end diff --git a/src/utils/jump_utils.jl b/src/utils/jump_utils.jl index 7380fe8..4918cff 100644 --- a/src/utils/jump_utils.jl +++ b/src/utils/jump_utils.jl @@ -633,10 +633,13 @@ function write_lp_file(jump_model::JuMP.Model, save_path::String) return end -# check_conflict_status functions can't be tested on CI because free solvers don't support IIS +# check_conflict_status functions can't be tested on CI because free solvers don't support IIS. +# The element type is matched covariantly (`<:JuMP.ConstraintRef`) so containers whose eltype is +# a concrete `ConstraintRef{...}` parametrization dispatch here too; an invariant `JuMP.ConstraintRef` +# signature silently misses them and raises a MethodError that aborts the whole conflict report. function check_conflict_status( jump_model::JuMP.Model, - constraint_container::DenseAxisArray{JuMP.ConstraintRef}, + constraint_container::DenseAxisArray{<:JuMP.ConstraintRef}, ) dims = axes(constraint_container) conflict_indices = Vector{Tuple{eltype.(dims)...}}() diff --git a/test/InfrastructureOptimizationModelsTests.jl b/test/InfrastructureOptimizationModelsTests.jl index 750c9a5..7e0e315 100644 --- a/test/InfrastructureOptimizationModelsTests.jl +++ b/test/InfrastructureOptimizationModelsTests.jl @@ -128,6 +128,7 @@ function run_tests() include(joinpath(TEST_DIR, "test_quadratic_curve.jl")) include(joinpath(TEST_DIR, "test_start_up_shut_down.jl")) include(joinpath(TEST_DIR, "test_ts_value_curve_objective.jl")) + include(joinpath(TEST_DIR, "test_tranche_axis_helpers.jl")) # --- common_models/, utils/, initial_conditions/ --- # TODO tests? diff --git a/test/Project.toml b/test/Project.toml index 3b55868..3382525 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -29,8 +29,8 @@ UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" UnoSolver = "1baa60ac-02f7-4b39-a7a8-2f4f58486b05" [sources] +InfrastructureSystems = {url = "https://github.com/Sienna-Platform/InfrastructureSystems.jl", rev = "IS4"} InfrastructureOptimizationModels = {path = ".."} -InfrastructureSystems = {path = "../../InfrastructureSystems.jl"} [compat] HiGHS = "1" diff --git a/test/test_settings.jl b/test/test_settings.jl index 263ed5f..47f8b45 100644 --- a/test/test_settings.jl +++ b/test/test_settings.jl @@ -39,7 +39,8 @@ end @test PSI.get_export_pwl_vars(settings) == false @test PSI.get_allow_fails(settings) == false @test PSI.get_rebuild_model(settings) == false - @test PSI.get_export_optimization_model(settings) == false + @test PSI.get_export_optimization_model(settings) == + PSI.OptimizationModelExportFormat.NONE @test PSI.get_store_variable_names(settings) == false @test PSI.get_check_numerical_bounds(settings) == true @test PSI.get_ext(settings) isa Dict{String, Any} diff --git a/test/test_tranche_axis_helpers.jl b/test/test_tranche_axis_helpers.jl new file mode 100644 index 0000000..7af564e --- /dev/null +++ b/test/test_tranche_axis_helpers.jl @@ -0,0 +1,74 @@ +""" +Unit tests for the PSY-free tranche-axis helpers used to size and fill the 3-D +`(component, tranche, time)` PWL parameter arrays for time-varying ORDC / market-bid +costs. Helpers live in `src/objective_function/value_curve_cost.jl` (Section 6). +""" + +@testset "Tranche-axis helpers (time-varying PWL)" begin + @testset "make_tranche_axis" begin + @test IOM.make_tranche_axis(3) == ["tranche_1", "tranche_2", "tranche_3"] + @test IOM.make_tranche_axis(1) == ["tranche_1"] + @test isempty(IOM.make_tranche_axis(0)) + end + + @testset "lookup_additional_axes" begin + # (component, tranche, time) -> only the middle (tranche) axis is "additional" + arr3 = DenseAxisArray(zeros(1, 2, 4), ["g1"], ["tranche_1", "tranche_2"], 1:4) + @test IOM.lookup_additional_axes(arr3) == (["tranche_1", "tranche_2"],) + + # (component, time) -> no additional axes + arr2 = DenseAxisArray(zeros(1, 4), ["g1"], 1:4) + @test IOM.lookup_additional_axes(arr2) == () + end + + @testset "get_max_tranches" begin + # length(PiecewiseStepData) == number of segments == length(x_coords) - 1 + a = IS.PiecewiseStepData([0.0, 1.0, 2.0], [10.0, 20.0]) # 2 tranches + b = IS.PiecewiseStepData([0.0, 1.0, 2.0, 3.0], [5.0, 6.0, 7.0]) # 3 tranches + + # Vector method: global maximum over time + @test IOM.get_max_tranches([a, b]) == 3 + @test IOM.get_max_tranches([a]) == 2 + + # TimeArray method: unwraps to its values and reuses the Vector method + ts = IOM.TS.TimeArray([DateTime(2024, 1, 1), DateTime(2024, 1, 1, 1)], [a, b]) + @test IOM.get_max_tranches(ts) == 3 + + # Dict method: global maximum across all entries + @test IOM.get_max_tranches(Dict("x" => [a], "y" => [a, b])) == 3 + end + + @testset "unwrap_for_param" begin + # 2 segments / 3 breakpoints + psd = IS.PiecewiseStepData([0.0, 1.0, 2.0], [10.0, 20.0]) + + # default fallback: any non-PWL parameter passes its element through unchanged + @test IOM.unwrap_for_param(TestParameterType(), 3.0, (1:5,)) == 3.0 + + @testset "slope parameter (pads y-coords with 0.0)" begin + slope_param = IOM.DecrementalPiecewiseLinearSlopeParameter() + # exact length: no padding + @test IOM.unwrap_for_param(slope_param, psd, (IOM.make_tranche_axis(2),)) == + [10.0, 20.0] + # short curve: pad slope = 0.0 up to the tranche-axis length + @test IOM.unwrap_for_param(slope_param, psd, (IOM.make_tranche_axis(4),)) == + [10.0, 20.0, 0.0, 0.0] + # too many coords for the axis: error + @test_throws ErrorException IOM.unwrap_for_param( + slope_param, psd, (IOM.make_tranche_axis(1),)) + end + + @testset "breakpoint parameter (repeats last breakpoint)" begin + bp_param = IOM.DecrementalPiecewiseLinearBreakpointParameter() + # exact length (tranches + 1): no padding + @test IOM.unwrap_for_param(bp_param, psd, (IOM.make_tranche_axis(3),)) == + [0.0, 1.0, 2.0] + # short curve: pad by repeating the last breakpoint (so dx = 0) + @test IOM.unwrap_for_param(bp_param, psd, (IOM.make_tranche_axis(5),)) == + [0.0, 1.0, 2.0, 2.0, 2.0] + # too many coords for the axis: error + @test_throws ErrorException IOM.unwrap_for_param( + bp_param, psd, (IOM.make_tranche_axis(2),)) + end + end +end