From b7a71b63d8ff406383e81770dab81553f4f8c2f0 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Tue, 14 Apr 2026 12:04:59 -0400 Subject: [PATCH 01/17] saving copilot's code here --- src/bilinear_approximations/bin2.jl | 40 ++++---- src/bilinear_approximations/hybs.jl | 57 ++++++----- src/bilinear_approximations/mccormick.jl | 106 ++++++++++++++++++-- src/bilinear_approximations/nmdt.jl | 77 ++++++-------- src/bilinear_approximations/no_approx.jl | 14 +-- src/quadratic_approximations/common.jl | 17 ++-- src/quadratic_approximations/epigraph.jl | 25 ++--- src/quadratic_approximations/manual_sos2.jl | 21 ++-- src/quadratic_approximations/nmdt.jl | 40 ++++---- src/quadratic_approximations/nmdt_common.jl | 77 +++++++------- src/quadratic_approximations/no_approx.jl | 8 +- src/quadratic_approximations/pwmcc_cuts.jl | 36 +++---- src/quadratic_approximations/sawtooth.jl | 33 +++--- src/quadratic_approximations/solver_sos2.jl | 21 ++-- 14 files changed, 311 insertions(+), 261 deletions(-) diff --git a/src/bilinear_approximations/bin2.jl b/src/bilinear_approximations/bin2.jl index 91634b62..5b77261b 100644 --- a/src/bilinear_approximations/bin2.jl +++ b/src/bilinear_approximations/bin2.jl @@ -34,9 +34,13 @@ Bin2Config(quad_config::QuadraticApproxConfig) = Bin2Config(quad_config, true) # --- Unified bilinear approximation dispatch --- """ - _add_bilinear_approx!(config::Bin2Config, container, C, names, time_steps, x_var, y_var, x_min, x_max, y_min, y_max, meta) + _add_bilinear_approx!(config::Bin2Config, container, C, names, time_steps, x_var, y_var, x_bounds, y_bounds, meta) Standard form: compute x² and y² quadratic approximations, then delegate to precomputed form. + +# Arguments +- `x_bounds::Vector{MinMax}`: per-name lower and upper bounds of x +- `y_bounds::Vector{MinMax}`: per-name lower and upper bounds of y """ function _add_bilinear_approx!( config::Bin2Config, @@ -46,32 +50,34 @@ function _add_bilinear_approx!( time_steps::UnitRange{Int}, x_var, y_var, - x_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} xsq = _add_quadratic_approx!( config.quad_config, container, C, names, time_steps, - x_var, x_min, x_max, meta * "_x", + x_var, x_bounds, meta * "_x", ) ysq = _add_quadratic_approx!( config.quad_config, container, C, names, time_steps, - y_var, y_min, y_max, meta * "_y", + y_var, y_bounds, meta * "_y", ) return _add_bilinear_approx!( config, container, C, names, time_steps, xsq, ysq, x_var, y_var, - x_min, x_max, y_min, y_max, meta, + x_bounds, y_bounds, meta, ) end """ - _add_bilinear_approx!(config::Bin2Config, container, C, names, time_steps, xsq, ysq, x_var, y_var, x_min, x_max, y_min, y_max, meta) + _add_bilinear_approx!(config::Bin2Config, container, C, names, time_steps, xsq, ysq, x_var, y_var, x_bounds, y_bounds, meta) Precomputed form: Bin2 identity z = ½((x+y)² − x² − y²) with optional PWMCC concave cuts. Accepts pre-computed quadratic approximations `xsq` ≈ x² and `ysq` ≈ y². + +# Arguments +- `x_bounds::Vector{MinMax}`: per-name lower and upper bounds of x +- `y_bounds::Vector{MinMax}`: per-name lower and upper bounds of y """ function _add_bilinear_approx!( config::Bin2Config, @@ -83,18 +89,14 @@ function _add_bilinear_approx!( ysq, x_var, y_var, - x_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} # --- Bin2 identity: z = ½((x+y)² − x² − y²) --- - # Bounds for p = x + y - p_min = x_min + y_min - p_max = x_max + y_max - IS.@assert_op p_min <= p_max + # Bounds for p = x + y (per-name) + p_bounds = [MinMax((min = x_bounds[i].min + y_bounds[i].min, max = x_bounds[i].max + y_bounds[i].max)) for i in eachindex(x_bounds)] meta_plus = meta * "_plus" @@ -116,7 +118,7 @@ function _add_bilinear_approx!( # Approximate p² = (x+y)² using the provided quadratic config psq = _add_quadratic_approx!( config.quad_config, container, C, names, time_steps, - p_expr, p_min, p_max, meta_plus, + p_expr, p_bounds, meta_plus, ) result_expr = add_expression_container!( @@ -141,7 +143,7 @@ function _add_bilinear_approx!( _add_reformulated_mccormick!( container, C, names, time_steps, x_var, y_var, psq, xsq, ysq, - x_min, x_max, y_min, y_max, meta, + x_bounds, y_bounds, meta, ) end diff --git a/src/bilinear_approximations/hybs.jl b/src/bilinear_approximations/hybs.jl index be72fbd9..3ab8da61 100644 --- a/src/bilinear_approximations/hybs.jl +++ b/src/bilinear_approximations/hybs.jl @@ -28,9 +28,13 @@ HybSConfig(quad_config::QuadraticApproxConfig, epigraph_depth::Int) = # --- Unified HybS dispatch methods --- """ - _add_bilinear_approx!(config::HybSConfig, container, C, names, time_steps, x_var, y_var, x_min, x_max, y_min, y_max, meta) + _add_bilinear_approx!(config::HybSConfig, container, C, names, time_steps, x_var, y_var, x_bounds, y_bounds, meta) Approximate x·y using HybS (Hybrid Separable) relaxation with config-selected quadratic method. + +# Arguments +- `x_bounds::Vector{MinMax}`: per-name lower and upper bounds of x +- `y_bounds::Vector{MinMax}`: per-name lower and upper bounds of y """ function _add_bilinear_approx!( config::HybSConfig, @@ -40,29 +44,27 @@ function _add_bilinear_approx!( time_steps::UnitRange{Int}, x_var, y_var, - x_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} xsq = _add_quadratic_approx!( config.quad_config, container, C, names, time_steps, - x_var, x_min, x_max, meta * "_x", + x_var, x_bounds, meta * "_x", ) ysq = _add_quadratic_approx!( config.quad_config, container, C, names, time_steps, - y_var, y_min, y_max, meta * "_y", + y_var, y_bounds, meta * "_y", ) return _add_bilinear_approx!( config, container, C, names, time_steps, xsq, ysq, x_var, y_var, - x_min, x_max, y_min, y_max, meta, + x_bounds, y_bounds, meta, ) end """ - _add_bilinear_approx!(config::HybSConfig, container, C, names, time_steps, xsq, ysq, x_var, y_var, x_min, x_max, y_min, y_max, meta) + _add_bilinear_approx!(config::HybSConfig, container, C, names, time_steps, xsq, ysq, x_var, y_var, x_bounds, y_bounds, meta) HybS bilinear approximation with pre-computed quadratic approximations for x² and y². @@ -71,6 +73,10 @@ Combines Bin2 and Bin3 separable identities: - Bin3 upper bound: z ≤ ½(z_x + z_y − z_p2) where z_p2 lower-bounds (x−y)² The cross-terms (x+y)² and (x−y)² always use epigraph Q^{L1} (pure LP). + +# Arguments +- `x_bounds::Vector{MinMax}`: per-name lower and upper bounds of x +- `y_bounds::Vector{MinMax}`: per-name lower and upper bounds of y """ function _add_bilinear_approx!( config::HybSConfig, @@ -82,19 +88,13 @@ function _add_bilinear_approx!( ysq, x_var, y_var, - x_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} - # Bounds for auxiliary variables - p1_min = x_min + y_min - p1_max = x_max + y_max - p2_min = x_min - y_max # p2 = x − y, so min uses −y_max - p2_max = x_max - y_min # and max uses −y_min - IS.@assert_op x_max > x_min - IS.@assert_op y_max > y_min + # Bounds for auxiliary variables (per-name) + p1_bounds = [MinMax((min = x_bounds[i].min + y_bounds[i].min, max = x_bounds[i].max + y_bounds[i].max)) for i in eachindex(x_bounds)] + p2_bounds = [MinMax((min = x_bounds[i].min - y_bounds[i].max, max = x_bounds[i].max - y_bounds[i].min)) for i in eachindex(x_bounds)] jump_model = get_jump_model(container) @@ -139,12 +139,12 @@ function _add_bilinear_approx!( zp1_expr = _add_quadratic_approx!( epi_cfg, container, C, names, time_steps, - p1_expr, p1_min, p1_max, meta_p1, + p1_expr, p1_bounds, meta_p1, ) zp2_expr = _add_quadratic_approx!( epi_cfg, container, C, names, time_steps, - p2_expr, p2_min, p2_max, meta_p2, + p2_expr, p2_bounds, meta_p2, ) # --- Create z variable and two-sided HybS bounds --- @@ -175,11 +175,14 @@ function _add_bilinear_approx!( meta, ) - # Compute valid bounds for z ≈ x·y from variable bounds - z_lo = min(x_min * y_min, x_min * y_max, x_max * y_min, x_max * y_max) - z_hi = max(x_min * y_min, x_min * y_max, x_max * y_min, x_max * y_max) + for (i, name) in enumerate(names), t in time_steps + xb = x_bounds[i] + yb = y_bounds[i] + + # Compute valid bounds for z ≈ x·y from variable bounds + z_lo = min(xb.min * yb.min, xb.min * yb.max, xb.max * yb.min, xb.max * yb.max) + z_hi = max(xb.min * yb.min, xb.min * yb.max, xb.max * yb.min, xb.max * yb.max) - for name in names, t in time_steps z = z_var[name, t] = JuMP.@variable( jump_model, @@ -212,7 +215,7 @@ function _add_bilinear_approx!( _add_mccormick_envelope!( container, C, names, time_steps, x_var, y_var, z_var, - x_min, x_max, y_min, y_max, meta, + x_bounds, y_bounds, meta, ) end diff --git a/src/bilinear_approximations/mccormick.jl b/src/bilinear_approximations/mccormick.jl index 34133ba5..e86fba0e 100644 --- a/src/bilinear_approximations/mccormick.jl +++ b/src/bilinear_approximations/mccormick.jl @@ -173,6 +173,90 @@ function _add_mccormick_envelope!( return end +""" + _add_mccormick_envelope!(container, C, names, time_steps, x_var, y_var, z_var, x_bounds, y_bounds, meta) + +Add McCormick envelope constraints for the bilinear product z ≈ x·y with per-name bounds. + +For each (name, t), adds 4 linear inequalities using bounds looked up by name index. + +# Arguments +- `container::OptimizationContainer`: the optimization container +- `::Type{C}`: component type +- `names::Vector{String}`: component names +- `time_steps::UnitRange{Int}`: time periods +- `x_var`: container of x variables indexed by (name, t) +- `y_var`: container of y variables indexed by (name, t) +- `z_var`: container of z variables indexed by (name, t) +- `x_bounds::Vector{MinMax}`: per-name lower and upper bounds of x +- `y_bounds::Vector{MinMax}`: per-name lower and upper bounds of y +- `meta::String`: identifier for container keys + +# Returns +- Nothing. Constraints are added in-place. +""" +function _add_mccormick_envelope!( + container::OptimizationContainer, + ::Type{C}, + names::Vector{String}, + time_steps::UnitRange{Int}, + x_var, + y_var, + z_var, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, + meta::String; + lower_bounds::Bool = true, +) where {C <: IS.InfrastructureSystemsComponent} + jump_model = get_jump_model(container) + + mc_cons = add_constraints_container!( + container, + McCormickConstraint(), + C, + names, + 1:4, + time_steps; + sparse = true, + meta, + ) + + for (i, name) in enumerate(names), t in time_steps + xb = x_bounds[i] + yb = y_bounds[i] + IS.@assert_op xb.max > xb.min + IS.@assert_op yb.max > yb.min + _add_mccormick_envelope!( + jump_model, mc_cons, (name, t), + x_var[name, t], y_var[name, t], z_var[name, t], + xb.min, xb.max, yb.min, yb.max; + lower_bounds, + ) + end + + return +end + +function _add_mccormick_envelope!( + container::OptimizationContainer, + ::Type{C}, + names::Vector{String}, + time_steps::UnitRange{Int}, + x_var, + z_var, + bounds::Vector{MinMax}, + meta::String; + lower_bounds::Bool = true, +) where {C <: IS.InfrastructureSystemsComponent} + _add_mccormick_envelope!( + container, C, names, time_steps, + x_var, x_var, z_var, + bounds, bounds, + meta; lower_bounds, + ) + return +end + function _add_mccormick_envelope!( jump_model::JuMP.Model, cons, @@ -265,12 +349,16 @@ function _add_reformulated_mccormick_bin2!( end """ - _add_reformulated_mccormick!(container, C, names, time_steps, x_var, y_var, psq, xsq, ysq, x_min, x_max, y_min, y_max, meta) + _add_reformulated_mccormick!(container, C, names, time_steps, x_var, y_var, psq, xsq, ysq, x_bounds, y_bounds, meta) Add 4 reformulated McCormick cuts for Bin2 separable bilinear approximation. Substitutes z = ½(z_p1 − z_x − z_y) into the standard McCormick envelope. `psq`, `xsq`, `ysq` are expression containers for (x+y)², x², y² approximations. + +# Arguments +- `x_bounds::Vector{MinMax}`: per-name lower and upper bounds of x +- `y_bounds::Vector{MinMax}`: per-name lower and upper bounds of y """ function _add_reformulated_mccormick!( container::OptimizationContainer, @@ -282,14 +370,10 @@ function _add_reformulated_mccormick!( psq, xsq, ysq, - x_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} - IS.@assert_op x_max > x_min - IS.@assert_op y_max > y_min jump_model = get_jump_model(container) mc_cons = add_constraints_container!( @@ -303,12 +387,16 @@ function _add_reformulated_mccormick!( meta, ) - for name in names, t in time_steps + for (i, name) in enumerate(names), t in time_steps + xb = x_bounds[i] + yb = y_bounds[i] + IS.@assert_op xb.max > xb.min + IS.@assert_op yb.max > yb.min _add_reformulated_mccormick_bin2!( jump_model, mc_cons, (name, t), x_var[name, t], y_var[name, t], psq[name, t], xsq[name, t], ysq[name, t], - x_min, x_max, y_min, y_max, + xb.min, xb.max, yb.min, yb.max, ) end diff --git a/src/bilinear_approximations/nmdt.jl b/src/bilinear_approximations/nmdt.jl index bed70958..247ca71f 100644 --- a/src/bilinear_approximations/nmdt.jl +++ b/src/bilinear_approximations/nmdt.jl @@ -28,7 +28,7 @@ end # --- DNMDT bilinear approximation --- """ - _add_bilinear_approx!(config::DNMDTBilinearConfig, container, C, names, time_steps, x_disc, y_disc, meta) + _add_bilinear_approx!(config::DNMDTBilinearConfig, container, C, names, time_steps, x_disc, y_disc, x_bounds, y_bounds, meta) Approximate x·y using the DNMDT method from pre-built discretizations. @@ -44,6 +44,8 @@ container. - `time_steps::UnitRange{Int}`: time periods - `x_disc::NMDTDiscretization`: pre-built discretization for x - `y_disc::NMDTDiscretization`: pre-built discretization for y +- `x_bounds::Vector{MinMax}`: per-name lower and upper bounds of x +- `y_bounds::Vector{MinMax}`: per-name lower and upper bounds of y - `meta::String`: identifier encoding the original variable type being approximated """ function _add_bilinear_approx!( @@ -54,10 +56,8 @@ function _add_bilinear_approx!( time_steps::UnitRange{Int}, x_disc::NMDTDiscretization, y_disc::NMDTDiscretization, - x_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} bx_yh_expr = _binary_continuous_product!( @@ -84,13 +84,13 @@ function _add_bilinear_approx!( return _assemble_dnmdt!( container, C, names, time_steps, bx_yh_expr, by_dx_expr, by_xh_expr, bx_dy_expr, - x_disc, y_disc, x_min, x_max, y_min, y_max, + x_disc, y_disc, x_bounds, y_bounds, config.depth, meta; result_type = BilinearProductExpression, ) end """ - _add_bilinear_approx!(config::DNMDTBilinearConfig, container, C, names, time_steps, x_var, y_var, x_min, x_max, y_min, y_max, meta) + _add_bilinear_approx!(config::DNMDTBilinearConfig, container, C, names, time_steps, x_var, y_var, x_bounds, y_bounds, meta) Approximate x·y using the DNMDT method from raw variable inputs. @@ -105,10 +105,8 @@ pre-discretized overload. - `time_steps::UnitRange{Int}`: time periods - `x_var`: container of x variables indexed by (name, t) - `y_var`: container of y variables indexed by (name, t) -- `x_min::Float64`: lower bound of x -- `x_max::Float64`: upper bound of x -- `y_min::Float64`: lower bound of y -- `y_max::Float64`: upper bound of y +- `x_bounds::Vector{MinMax}`: per-name lower and upper bounds of x +- `y_bounds::Vector{MinMax}`: per-name lower and upper bounds of y - `meta::String`: identifier encoding the original variable type being approximated """ function _add_bilinear_approx!( @@ -119,10 +117,8 @@ function _add_bilinear_approx!( time_steps::UnitRange{Int}, x_var, y_var, - x_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} x_disc = _discretize!( @@ -131,8 +127,7 @@ function _add_bilinear_approx!( names, time_steps, x_var, - x_min, - x_max, + x_bounds, config.depth, meta * "_x", ) @@ -142,8 +137,7 @@ function _add_bilinear_approx!( names, time_steps, y_var, - y_min, - y_max, + y_bounds, config.depth, meta * "_y", ) @@ -155,10 +149,8 @@ function _add_bilinear_approx!( time_steps, x_disc, y_disc, - x_min, - x_max, - y_min, - y_max, + x_bounds, + y_bounds, meta, ) end @@ -166,7 +158,7 @@ end # --- NMDT bilinear approximation --- """ - _add_bilinear_approx!(config::NMDTBilinearConfig, container, C, names, time_steps, x_disc, yh_expr, y_min, y_max, meta) + _add_bilinear_approx!(config::NMDTBilinearConfig, container, C, names, time_steps, x_disc, yh_expr, x_bounds, y_bounds, meta) Approximate x·y using the NMDT method from a pre-built x discretization and normalized y. @@ -182,8 +174,8 @@ x·y via `_assemble_product!`. Stores results in a `BilinearProductExpression` c - `time_steps::UnitRange{Int}`: time periods - `x_disc::NMDTDiscretization`: pre-built discretization for x - `yh_expr`: expression container for the normalized variable yh = (y − y_min)/(y_max − y_min) -- `y_min::Float64`: lower bound of y -- `y_max::Float64`: upper bound of y +- `x_bounds::Vector{MinMax}`: per-name lower and upper bounds of x +- `y_bounds::Vector{MinMax}`: per-name lower and upper bounds of y - `meta::String`: identifier encoding the original variable type being approximated """ function _add_bilinear_approx!( @@ -194,10 +186,8 @@ function _add_bilinear_approx!( time_steps::UnitRange{Int}, x_disc::NMDTDiscretization, yh_expr, - x_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, meta::String; ) where {C <: IS.InfrastructureSystemsComponent} bx_y_expr = _binary_continuous_product!( @@ -213,13 +203,13 @@ function _add_bilinear_approx!( return _assemble_product!( container, C, names, time_steps, [bx_y_expr], dz, - x_disc, yh_expr, x_min, x_max, y_min, y_max, + x_disc, yh_expr, x_bounds, y_bounds, meta; result_type = BilinearProductExpression, ) end """ - _add_bilinear_approx!(config::NMDTBilinearConfig, container, C, names, time_steps, x_var, y_var, x_min, x_max, y_min, y_max, meta) + _add_bilinear_approx!(config::NMDTBilinearConfig, container, C, names, time_steps, x_var, y_var, x_bounds, y_bounds, meta) Approximate x·y using the NMDT method from raw variable inputs. @@ -234,10 +224,8 @@ delegates to the pre-discretized overload. - `time_steps::UnitRange{Int}`: time periods - `x_var`: container of x variables indexed by (name, t) - `y_var`: container of y variables indexed by (name, t) -- `x_min::Float64`: lower bound of x -- `x_max::Float64`: upper bound of x -- `y_min::Float64`: lower bound of y -- `y_max::Float64`: upper bound of y +- `x_bounds::Vector{MinMax}`: per-name lower and upper bounds of x +- `y_bounds::Vector{MinMax}`: per-name lower and upper bounds of y - `meta::String`: identifier encoding the original variable type being approximated """ function _add_bilinear_approx!( @@ -248,10 +236,8 @@ function _add_bilinear_approx!( time_steps::UnitRange{Int}, x_var, y_var, - x_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} x_disc = _discretize!( @@ -260,13 +246,12 @@ function _add_bilinear_approx!( names, time_steps, x_var, - x_min, - x_max, + x_bounds, config.depth, meta * "_x", ) yh_expr = - _normed_variable!(container, C, names, time_steps, y_var, y_min, y_max, meta * "_y") + _normed_variable!(container, C, names, time_steps, y_var, y_bounds, meta * "_y") return _add_bilinear_approx!( config, container, @@ -275,10 +260,8 @@ function _add_bilinear_approx!( time_steps, x_disc, yh_expr, - x_min, - x_max, - y_min, - y_max, + x_bounds, + y_bounds, meta, ) end diff --git a/src/bilinear_approximations/no_approx.jl b/src/bilinear_approximations/no_approx.jl index e0ff1f9f..8a7f1007 100644 --- a/src/bilinear_approximations/no_approx.jl +++ b/src/bilinear_approximations/no_approx.jl @@ -5,7 +5,7 @@ struct NoBilinearApproxConfig <: BilinearApproxConfig end """ - _add_bilinear_approx!(::NoBilinearApproxConfig, container, C, names, time_steps, x_var, y_var, x_min, x_max, y_min, y_max, meta) + _add_bilinear_approx!(::NoBilinearApproxConfig, container, C, names, time_steps, x_var, y_var, x_bounds, y_bounds, meta) No-op bilinear approximation: returns exact x·y as a QuadExpr. @@ -17,10 +17,8 @@ No-op bilinear approximation: returns exact x·y as a QuadExpr. - `time_steps::UnitRange{Int}`: time periods - `x_var`: container of x variables indexed by (name, t) - `y_var`: container of y variables indexed by (name, t) -- `x_min::Float64`: lower bound of x domain -- `x_max::Float64`: upper bound of x domain -- `y_min::Float64`: lower bound of y domain -- `y_max::Float64`: upper bound of y domain +- `x_bounds::Vector{MinMax}`: per-name lower and upper bounds of x domain +- `y_bounds::Vector{MinMax}`: per-name lower and upper bounds of y domain - `meta::String`: variable type identifier for the approximation """ function _add_bilinear_approx!( @@ -31,10 +29,8 @@ function _add_bilinear_approx!( time_steps::UnitRange{Int}, x_var, y_var, - x_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} result_expr = add_expression_container!( diff --git a/src/quadratic_approximations/common.jl b/src/quadratic_approximations/common.jl index 8c560966..be37adfa 100644 --- a/src/quadratic_approximations/common.jl +++ b/src/quadratic_approximations/common.jl @@ -10,7 +10,7 @@ struct QuadraticExpression <: ExpressionType end abstract type QuadraticApproxConfig end """ - _normed_variable!(container, C, names, time_steps, x_var, x_min, x_max, meta) + _normed_variable!(container, C, names, time_steps, x_var, bounds, meta) Create an affine expression for the normalized variable xh = (x − x_min) / (x_max − x_min) ∈ [0,1]. @@ -22,8 +22,7 @@ Stores results in a `NormedVariableExpression` expression container. - `names::Vector{String}`: component names - `time_steps::UnitRange{Int}`: time periods - `x_var`: container of variables indexed by (name, t) -- `x_min::Float64`: lower bound of x domain -- `x_max::Float64`: upper bound of x domain +- `bounds::Vector{MinMax}`: per-name lower and upper bounds of x domain - `meta::String`: identifier encoding the original variable type being approximated """ function _normed_variable!( @@ -32,12 +31,9 @@ function _normed_variable!( names::Vector{String}, time_steps::UnitRange{Int}, x_var, - x_min::Float64, - x_max::Float64, + bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} - IS.@assert_op x_max > x_min - lx = x_max - x_min result_expr = add_expression_container!( container, NormedVariableExpression, @@ -47,9 +43,12 @@ function _normed_variable!( meta, ) - for name in names, t in time_steps + for (i, name) in enumerate(names), t in time_steps + b = bounds[i] + IS.@assert_op b.max > b.min + lx = b.max - b.min result = result_expr[name, t] = JuMP.AffExpr(0.0) - add_linear_to_jump_expression!(result, x_var[name, t], 1.0 / lx, -x_min / lx) + add_linear_to_jump_expression!(result, x_var[name, t], 1.0 / lx, -b.min / lx) end return result_expr end diff --git a/src/quadratic_approximations/epigraph.jl b/src/quadratic_approximations/epigraph.jl index 063f5456..2b180b6c 100644 --- a/src/quadratic_approximations/epigraph.jl +++ b/src/quadratic_approximations/epigraph.jl @@ -24,7 +24,7 @@ struct EpigraphQuadConfig <: QuadraticApproxConfig end """ - _add_quadratic_approx!(::EpigraphQuadConfig, container, C, names, time_steps, x_var, x_min, x_max, meta) + _add_quadratic_approx!(::EpigraphQuadConfig, container, C, names, time_steps, x_var, bounds, meta) Create a variable z that lower-bounds x² using tangent-line cuts (Q^{L1} relaxation). @@ -44,8 +44,7 @@ The maximum underestimation gap between the tangent envelope and x² is - `names::Vector{String}`: component names - `time_steps::UnitRange{Int}`: time periods - `x_var`: container of variables indexed by (name, t) -- `x_min::Float64`: lower bound of x domain -- `x_max::Float64`: upper bound of x domain +- `bounds::Vector{MinMax}`: per-name lower and upper bounds of x domain - `meta::String`: variable type identifier for the approximated variable """ function _add_quadratic_approx!( @@ -55,14 +54,11 @@ function _add_quadratic_approx!( names::Vector{String}, time_steps::UnitRange{Int}, x_var, - x_min::Float64, - x_max::Float64, + bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} - IS.@assert_op x_max > x_min IS.@assert_op config.depth >= 1 jump_model = get_jump_model(container) - delta = x_max - x_min g_levels = 0:(config.depth) z_var = add_variable_container!( @@ -126,10 +122,11 @@ function _add_quadratic_approx!( meta, ) - # Upper bound for epigraph variable z ≈ x² - z_ub = max(x_min^2, x_max^2) - - for name in names, t in time_steps + for (i, name) in enumerate(names), t in time_steps + b = bounds[i] + IS.@assert_op b.max > b.min + delta = b.max - b.min + z_ub = max(b.min^2, b.max^2) x = x_var[name, t] # Auxiliary variables g_0,...,g_L ∈ [0, 1] @@ -146,7 +143,7 @@ function _add_quadratic_approx!( # Linking constraint: g_0 = (x - x_min) / Δ link_cons[name, t] = JuMP.@constraint( jump_model, - g0 == (x - x_min) / delta, + g0 == (x - b.min) / delta, ) # T^L constraints for j = 1,...,L @@ -180,13 +177,13 @@ function _add_quadratic_approx!( tangent_cons[(name, j + 1, t)] = JuMP.@constraint( jump_model, z >= - x_min * (2 * delta * g0 + x_min) - fL + delta^2 * (g0 - 2.0^(-2j - 2)) + b.min * (2 * delta * g0 + b.min) - fL + delta^2 * (g0 - 2.0^(-2j - 2)) ) end tangent_cons[name, 1, t] = JuMP.@constraint(jump_model, z >= 0) tangent_cons[name, config.depth + 1, t] = JuMP.@constraint( jump_model, - z >= 2.0 * x_min - 1.0 + 2.0 * delta * g0 + z >= 2.0 * b.min - 1.0 + 2.0 * delta * g0 ) result_expr[name, t] = JuMP.AffExpr(0.0, z => 1.0) diff --git a/src/quadratic_approximations/manual_sos2.jl b/src/quadratic_approximations/manual_sos2.jl index 86b572c1..d061708a 100644 --- a/src/quadratic_approximations/manual_sos2.jl +++ b/src/quadratic_approximations/manual_sos2.jl @@ -24,7 +24,7 @@ end ManualSOS2QuadConfig(depth::Int) = ManualSOS2QuadConfig(depth, 4) """ - _add_quadratic_approx!(config::ManualSOS2QuadConfig, container, C, names, time_steps, x_var, x_min, x_max, meta) + _add_quadratic_approx!(config::ManualSOS2QuadConfig, container, C, names, time_steps, x_var, bounds, meta) Approximate x² using a piecewise linear function with manually-implemented SOS2 constraints. @@ -40,8 +40,7 @@ expression container. - `names::Vector{String}`: component names - `time_steps::UnitRange{Int}`: time periods - `x_var`: container of variables indexed by (name, t) -- `x_min::Float64`: lower bound of x domain -- `x_max::Float64`: upper bound of x domain +- `bounds::Vector{MinMax}`: per-name lower and upper bounds of x domain - `meta::String`: variable type identifier for the approximation (allows multiple approximations per component type) """ function _add_quadratic_approx!( @@ -51,11 +50,9 @@ function _add_quadratic_approx!( names::Vector{String}, time_steps::UnitRange{Int}, x_var, - x_min::Float64, - x_max::Float64, + bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} - lx = x_max - x_min x_bkpts, x_sq_bkpts = _get_breakpoints_for_pwl_function( 0.0, @@ -137,7 +134,9 @@ function _add_quadratic_approx!( meta, ) - for name in names, t in time_steps + for (i, name) in enumerate(names), t in time_steps + b = bounds[i] + lx = b.max - b.min x = x_var[name, t] # Create lambda variables: λ_i ∈ [0, 1] @@ -157,7 +156,7 @@ function _add_quadratic_approx!( for i in eachindex(x_bkpts) add_proportional_to_jump_expression!(link, lambda[i], x_bkpts[i]) end - link_cons[name, t] = JuMP.@constraint(jump_model, (x - x_min) / lx == link) + link_cons[name, t] = JuMP.@constraint(jump_model, (x - b.min) / lx == link) # Σ λ_i = 1 norm = norm_expr[name, t] = JuMP.AffExpr(0.0) @@ -203,15 +202,15 @@ function _add_quadratic_approx!( end x_sq = JuMP.AffExpr(0.0) add_proportional_to_jump_expression!(x_sq, x_hat_sq, lx * lx) - add_proportional_to_jump_expression!(x_sq, x, 2 * x_min) - add_constant_to_jump_expression!(x_sq, -x_min * x_min) + add_proportional_to_jump_expression!(x_sq, x, 2 * b.min) + add_constant_to_jump_expression!(x_sq, -b.min * b.min) result_expr[name, t] = x_sq end if config.pwmcc_segments > 0 _add_pwmcc_concave_cuts!( container, C, names, time_steps, - x_var, result_expr, x_min, x_max, + x_var, result_expr, bounds, config.pwmcc_segments, meta * "_pwmcc", ) end diff --git a/src/quadratic_approximations/nmdt.jl b/src/quadratic_approximations/nmdt.jl index 5a6049ae..343dbaec 100644 --- a/src/quadratic_approximations/nmdt.jl +++ b/src/quadratic_approximations/nmdt.jl @@ -34,7 +34,7 @@ end NMDTQuadConfig(depth::Int) = NMDTQuadConfig(depth, 3 * depth) """ - _add_quadratic_approx!(config::DNMDTQuadConfig, container, C, names, time_steps, x_disc, meta) + _add_quadratic_approx!(config::DNMDTQuadConfig, container, C, names, time_steps, x_disc, bounds, meta) Approximate x² using the Double NMDT (DNMDT) method from a pre-built discretization. @@ -49,6 +49,7 @@ tightens lower bounds with an epigraph relaxation via `_tighten_lower_bounds!`. - `names::Vector{String}`: component names - `time_steps::UnitRange{Int}`: time periods - `x_disc::NMDTDiscretization`: pre-built discretization for x +- `bounds::Vector{MinMax}`: per-name lower and upper bounds of x domain - `meta::String`: identifier encoding the original variable type being approximated """ function _add_quadratic_approx!( @@ -58,8 +59,7 @@ function _add_quadratic_approx!( names::Vector{String}, time_steps::UnitRange{Int}, x_disc::NMDTDiscretization, - x_min::Float64, - x_max::Float64, + bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} tighten = config.epigraph_depth > 0 @@ -77,7 +77,7 @@ function _add_quadratic_approx!( result_expr = _assemble_dnmdt!( container, C, names, time_steps, bx_xh_expr, bx_dx_expr, bx_xh_expr, bx_dx_expr, - x_disc, x_disc, x_min, x_max, x_min, x_max, + x_disc, x_disc, bounds, bounds, config.depth, meta; tighten, result_type = QuadraticExpression, ) @@ -93,7 +93,7 @@ function _add_quadratic_approx!( end """ - _add_quadratic_approx!(config::DNMDTQuadConfig, container, C, names, time_steps, x_var, x_min, x_max, meta) + _add_quadratic_approx!(config::DNMDTQuadConfig, container, C, names, time_steps, x_var, bounds, meta) Approximate x² using the Double NMDT (DNMDT) method from raw variable inputs. @@ -107,8 +107,7 @@ Stores results in a `QuadraticExpression` container. - `names::Vector{String}`: component names - `time_steps::UnitRange{Int}`: time periods - `x_var`: container of variables indexed by (name, t) -- `x_min::Float64`: lower bound of x domain -- `x_max::Float64`: upper bound of x domain +- `bounds::Vector{MinMax}`: per-name lower and upper bounds of x domain - `meta::String`: identifier encoding the original variable type being approximated """ function _add_quadratic_approx!( @@ -118,23 +117,22 @@ function _add_quadratic_approx!( names::Vector{String}, time_steps::UnitRange{Int}, x_var, - x_min::Float64, - x_max::Float64, + bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} x_disc = _discretize!( container, C, names, time_steps, - x_var, x_min, x_max, config.depth, meta, + x_var, bounds, config.depth, meta, ) return _add_quadratic_approx!( config, container, C, names, time_steps, - x_disc, x_min, x_max, meta, + x_disc, bounds, meta, ) end """ - _add_quadratic_approx!(config::NMDTQuadConfig, container, C, names, time_steps, x_disc, meta) + _add_quadratic_approx!(config::NMDTQuadConfig, container, C, names, time_steps, x_disc, bounds, meta) Approximate x² using the NMDT method from a pre-built discretization. @@ -149,6 +147,7 @@ container. Optionally tightens lower bounds with an epigraph relaxation. - `names::Vector{String}`: component names - `time_steps::UnitRange{Int}`: time periods - `x_disc::NMDTDiscretization`: pre-built discretization for x +- `bounds::Vector{MinMax}`: per-name lower and upper bounds of x domain - `meta::String`: identifier encoding the original variable type being approximated """ function _add_quadratic_approx!( @@ -158,8 +157,7 @@ function _add_quadratic_approx!( names::Vector{String}, time_steps::UnitRange{Int}, x_disc::NMDTDiscretization, - x_min::Float64, - x_max::Float64, + bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} tighten = config.epigraph_depth > 0 @@ -177,7 +175,7 @@ function _add_quadratic_approx!( result_expr = _assemble_product!( container, C, names, time_steps, [bx_y_expr], dz, - x_disc, x_disc, x_min, x_max, x_min, x_max, + x_disc, x_disc, bounds, bounds, meta; result_type = QuadraticExpression, ) @@ -192,7 +190,7 @@ function _add_quadratic_approx!( end """ - _add_quadratic_approx!(config::NMDTQuadConfig, container, C, names, time_steps, x_var, x_min, x_max, meta) + _add_quadratic_approx!(config::NMDTQuadConfig, container, C, names, time_steps, x_var, bounds, meta) Approximate x² using the NMDT method from raw variable inputs. @@ -206,8 +204,7 @@ Stores results in a `QuadraticExpression` container. - `names::Vector{String}`: component names - `time_steps::UnitRange{Int}`: time periods - `x_var`: container of variables indexed by (name, t) -- `x_min::Float64`: lower bound of x domain -- `x_max::Float64`: upper bound of x domain +- `bounds::Vector{MinMax}`: per-name lower and upper bounds of x domain - `meta::String`: identifier encoding the original variable type being approximated """ function _add_quadratic_approx!( @@ -217,17 +214,16 @@ function _add_quadratic_approx!( names::Vector{String}, time_steps::UnitRange{Int}, x_var, - x_min::Float64, - x_max::Float64, + bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} x_disc = _discretize!( container, C, names, time_steps, - x_var, x_min, x_max, config.depth, meta, + x_var, bounds, config.depth, meta, ) return _add_quadratic_approx!( config, container, C, names, time_steps, - x_disc, x_min, x_max, meta, + x_disc, bounds, meta, ) end diff --git a/src/quadratic_approximations/nmdt_common.jl b/src/quadratic_approximations/nmdt_common.jl index 3f45695b..3d589ec0 100644 --- a/src/quadratic_approximations/nmdt_common.jl +++ b/src/quadratic_approximations/nmdt_common.jl @@ -36,7 +36,7 @@ struct NMDTDiscretization end """ - _discretize!(container, C, names, time_steps, x_var, x_min, x_max, depth, meta) + _discretize!(container, C, names, time_steps, x_var, bounds, depth, meta) Discretize the normalized variable xh = (x − x_min)/(x_max − x_min) using L binary variables. @@ -50,8 +50,7 @@ returns an `NMDTDiscretization` struct holding all components. - `names::Vector{String}`: component names - `time_steps::UnitRange{Int}`: time periods - `x_var`: container of variables indexed by (name, t) -- `x_min::Float64`: lower bound of x domain -- `x_max::Float64`: upper bound of x domain +- `bounds::Vector{MinMax}`: per-name lower and upper bounds of x domain - `depth::Int`: number of binary discretization levels L - `meta::String`: identifier encoding the original variable type being approximated """ @@ -61,12 +60,10 @@ function _discretize!( names::Vector{String}, time_steps::UnitRange{Int}, x_var, - x_min::Float64, - x_max::Float64, + bounds::Vector{MinMax}, depth::Int, meta::String; ) where {C <: IS.InfrastructureSystemsComponent} - IS.@assert_op x_max > x_min IS.@assert_op depth >= 1 jump_model = get_jump_model(container) @@ -106,7 +103,7 @@ function _discretize!( xh_expr = _normed_variable!( container, C, names, time_steps, - x_var, x_min, x_max, meta, + x_var, bounds, meta, ) for name in names, t in time_steps @@ -259,7 +256,8 @@ function _tighten_lower_bounds!( epi_expr = _add_quadratic_approx!( EpigraphQuadConfig(epigraph_depth), container, C, names, time_steps, - x_disc.norm_expr, 0.0, 1.0, meta * "_epi", + x_disc.norm_expr, fill(MinMax((min = 0.0, max = 1.0)), length(names)), + meta * "_epi", ) epi_cons = add_constraints_container!( container, @@ -341,7 +339,7 @@ function _residual_product!( end """ - _assemble_product!(container, C, names, time_steps, terms, dz_var, x_disc, y_disc, meta; result_type) + _assemble_product!(container, C, names, time_steps, terms, dz_var, xh_norm, yh_norm, x_bounds, y_bounds, meta; result_type) Reconstruct the bilinear product x·y from normalized NMDT components. @@ -363,10 +361,8 @@ Stores results in an expression container of type `result_type`. - `dz_var`: variable container for the residual product δ·y - `xh_norm`: normed expression container for x - `yh_norm`: normed expression container for y -- `x_min::Float64`: original x min -- `x_max::Float64`: original x max -- `y_min::Float64`: original y min -- `y_max::Float64`: original y max +- `x_bounds::Vector{MinMax}`: per-name bounds for x +- `y_bounds::Vector{MinMax}`: per-name bounds for y - `meta::String`: identifier encoding the original variable type being approximated - `result_type`: expression type to store results in (default: `NMDTResultExpression`) """ @@ -379,16 +375,11 @@ function _assemble_product!( dz_var, xh_expr, yh_expr, - x_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, meta::String; result_type = NMDTResultExpression, ) where {C <: IS.InfrastructureSystemsComponent} - lx = x_max - x_min - ly = y_max - y_min - result_expr = add_expression_container!( container, result_type, @@ -398,7 +389,12 @@ function _assemble_product!( meta, ) - for name in names, t in time_steps + for (i, name) in enumerate(names), t in time_steps + xb = x_bounds[i] + yb = y_bounds[i] + lx = xb.max - xb.min + ly = yb.max - yb.min + result = result_expr[name, t] = JuMP.AffExpr(0.0) zh = JuMP.AffExpr(0.0) for term in terms @@ -407,16 +403,16 @@ function _assemble_product!( add_proportional_to_jump_expression!(zh, dz_var[name, t], 1.0) add_proportional_to_jump_expression!(result, zh, lx * ly) - add_proportional_to_jump_expression!(result, xh_expr[name, t], lx * y_min) - add_proportional_to_jump_expression!(result, yh_expr[name, t], ly * x_min) - add_constant_to_jump_expression!(result, x_min * y_min) + add_proportional_to_jump_expression!(result, xh_expr[name, t], lx * yb.min) + add_proportional_to_jump_expression!(result, yh_expr[name, t], ly * xb.min) + add_constant_to_jump_expression!(result, xb.min * yb.min) end return result_expr end """ - _assemble_product!(container, C, names, time_steps, terms, dz_var, x_disc::NMDTDiscretization, y_disc::NMDTDiscretization, x_min, x_max, y_min, y_max, meta; result_type) + _assemble_product!(container, C, names, time_steps, terms, dz_var, x_disc::NMDTDiscretization, y_disc::NMDTDiscretization, x_bounds, y_bounds, meta; result_type) Convenience overload: extracts `norm_expr` from both discretizations and delegates to the core `_assemble_product!`. """ @@ -429,23 +425,21 @@ function _assemble_product!( dz_var, x_disc::NMDTDiscretization, y_disc::NMDTDiscretization, - x_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, meta::String; result_type = NMDTResultExpression, ) where {C <: IS.InfrastructureSystemsComponent} return _assemble_product!( container, C, names, time_steps, terms, dz_var, x_disc.norm_expr, y_disc.norm_expr, - x_min, x_max, y_min, y_max, + x_bounds, y_bounds, meta; result_type, ) end """ - _assemble_product!(container, C, names, time_steps, terms, dz_var, x_disc::NMDTDiscretization, yh_expr, x_min, x_max, y_min, y_max, meta; result_type) + _assemble_product!(container, C, names, time_steps, terms, dz_var, x_disc::NMDTDiscretization, yh_expr, x_bounds, y_bounds, meta; result_type) Convenience overload: extracts `norm_expr` from x_disc and delegates to the core `_assemble_product!`. """ @@ -458,17 +452,15 @@ function _assemble_product!( dz_var, x_disc::NMDTDiscretization, yh_expr, - x_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, meta::String; result_type = NMDTResultExpression, ) where {C <: IS.InfrastructureSystemsComponent} return _assemble_product!( container, C, names, time_steps, terms, dz_var, x_disc.norm_expr, yh_expr, - x_min, x_max, y_min, y_max, + x_bounds, y_bounds, meta; result_type, ) end @@ -497,6 +489,9 @@ container of type `result_type`. - `bx_dy_expr`: expression for β_x·δy binary-continuous products - `x_disc::NMDTDiscretization`: discretization for x - `y_disc::NMDTDiscretization`: discretization for y +- `x_bounds::Vector{MinMax}`: per-name bounds for x +- `y_bounds::Vector{MinMax}`: per-name bounds for y +- `depth::Int`: number of binary discretization levels L - `meta::String`: identifier encoding the original variable type being approximated - `lambda::Float64`: convex combination weight for the two NMDT estimates (default: `DNMDT_LAMBDA` = 0.5) - `result_type`: expression type to store results in (default: `NMDTResultExpression`) @@ -512,10 +507,8 @@ function _assemble_dnmdt!( bx_dy_expr, x_disc::NMDTDiscretization, y_disc::NMDTDiscretization, - x_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, depth::Int, meta::String; lambda::Float64 = DNMDT_LAMBDA, @@ -539,13 +532,13 @@ function _assemble_dnmdt!( z1_expr = _assemble_product!( container, C, names, time_steps, [bx_yh_expr, by_dx_expr], dz, - x_disc, y_disc, x_min, x_max, y_min, y_max, + x_disc, y_disc, x_bounds, y_bounds, meta * "_nmdt1", ) z2_expr = _assemble_product!( container, C, names, time_steps, [by_xh_expr, bx_dy_expr], dz, - y_disc, x_disc, y_min, y_max, x_min, x_max, + y_disc, x_disc, y_bounds, x_bounds, meta * "_nmdt2", ) diff --git a/src/quadratic_approximations/no_approx.jl b/src/quadratic_approximations/no_approx.jl index 93baa057..e2d88498 100644 --- a/src/quadratic_approximations/no_approx.jl +++ b/src/quadratic_approximations/no_approx.jl @@ -5,7 +5,7 @@ struct NoQuadApproxConfig <: QuadraticApproxConfig end """ - _add_quadratic_approx!(::NoQuadApproxConfig, container, C, names, time_steps, x_var, x_min, x_max, meta) + _add_quadratic_approx!(::NoQuadApproxConfig, container, C, names, time_steps, x_var, bounds, meta) No-op quadratic approximation: returns exact x² as a QuadExpr. @@ -16,8 +16,7 @@ No-op quadratic approximation: returns exact x² as a QuadExpr. - `names::Vector{String}`: component names - `time_steps::UnitRange{Int}`: time periods - `x_var`: container of variables indexed by (name, t) -- `x_min::Float64`: lower bound of x domain -- `x_max::Float64`: upper bound of x domain +- `bounds::Vector{MinMax}`: per-name lower and upper bounds of x domain - `meta::String`: variable type identifier for the approximation """ function _add_quadratic_approx!( @@ -27,8 +26,7 @@ function _add_quadratic_approx!( names::Vector{String}, time_steps::UnitRange{Int}, x_var, - x_min::Float64, - x_max::Float64, + bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} result_expr = add_expression_container!( diff --git a/src/quadratic_approximations/pwmcc_cuts.jl b/src/quadratic_approximations/pwmcc_cuts.jl index cded3a5e..3643bc39 100644 --- a/src/quadratic_approximations/pwmcc_cuts.jl +++ b/src/quadratic_approximations/pwmcc_cuts.jl @@ -32,11 +32,11 @@ struct PiecewiseMcCormickTangentLBL <: ConstraintType end struct PiecewiseMcCormickTangentLBR <: ConstraintType end """ - _add_pwmcc_concave_cuts!(container, C, names, time_steps, v_var, q_expr, v_min, v_max, K, meta) + _add_pwmcc_concave_cuts!(container, C, names, time_steps, v_var, q_expr, bounds, K, meta) Add piecewise McCormick cuts on a concave term (-v^2) to tighten its SoS2 LP relaxation. -Partitions [v_min, v_max] into K uniform sub-intervals and adds disaggregated +Partitions each name's [v_min, v_max] into K uniform sub-intervals and adds disaggregated variables, binary interval selectors, and chord/tangent constraints that cut off the interior of the SoS2 relaxation polytope. @@ -47,8 +47,7 @@ the interior of the SoS2 relaxation polytope. - `time_steps::UnitRange{Int}`: time periods - `v_var`: container of the original variable indexed by (name, t) - `q_expr`: expression container for the SoS2 approximation of v^2 (indexed by (name, t)) -- `v_min::Float64`: lower bound of v domain -- `v_max::Float64`: upper bound of v domain +- `bounds::Vector{MinMax}`: per-name lower and upper bounds of v domain - `K::Int`: number of sub-intervals (K=2 is the minimal useful choice) - `meta::String`: unique key prefix, e.g. "pwmcc_x" or "pwmcc_y" """ @@ -59,25 +58,14 @@ function _add_pwmcc_concave_cuts!( time_steps::UnitRange{Int}, v_var, q_expr, - v_min::Float64, - v_max::Float64, + bounds::Vector{MinMax}, K::Int, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} IS.@assert_op K >= 1 - IS.@assert_op v_min < v_max jump_model = get_jump_model(container) - # Pre-compute breakpoints and derived coefficients (depend only on k, not on name/t) - brk = [v_min + k * (v_max - v_min) / K for k in 0:K] - sum_brk = [brk[k] + brk[k + 1] for k in 1:K] - prod_brk = [brk[k] * brk[k + 1] for k in 1:K] - two_brk_l = [2.0 * brk[k] for k in 1:K] - sq_brk_l = [brk[k]^2 for k in 1:K] - two_brk_r = [2.0 * brk[k + 1] for k in 1:K] - sq_brk_r = [brk[k + 1]^2 for k in 1:K] - # Create containers delta_var = add_variable_container!(container, PiecewiseMcCormickBinary, C; meta) vd_var = @@ -145,7 +133,21 @@ function _add_pwmcc_concave_cuts!( delta = Vector{JuMP.VariableRef}(undef, K) vd = Vector{JuMP.VariableRef}(undef, K) - for name in names, t in time_steps + for (idx, name) in enumerate(names), t in time_steps + b = bounds[idx] + IS.@assert_op b.min < b.max + v_min = b.min + v_max = b.max + + # Compute breakpoints and derived coefficients for this name + brk = [v_min + k * (v_max - v_min) / K for k in 0:K] + sum_brk = [brk[k] + brk[k + 1] for k in 1:K] + prod_brk = [brk[k] * brk[k + 1] for k in 1:K] + two_brk_l = [2.0 * brk[k] for k in 1:K] + sq_brk_l = [brk[k]^2 for k in 1:K] + two_brk_r = [2.0 * brk[k + 1] for k in 1:K] + sq_brk_r = [brk[k + 1]^2 for k in 1:K] + v = v_var[name, t] q = q_expr[name, t] diff --git a/src/quadratic_approximations/sawtooth.jl b/src/quadratic_approximations/sawtooth.jl index 1aafed6c..0a0d0d2e 100644 --- a/src/quadratic_approximations/sawtooth.jl +++ b/src/quadratic_approximations/sawtooth.jl @@ -31,7 +31,7 @@ end SawtoothQuadConfig(depth::Int) = SawtoothQuadConfig(depth, 0) """ - _add_quadratic_approx!(config::SawtoothQuadConfig, container, C, names, time_steps, x_var, x_min, x_max, meta) + _add_quadratic_approx!(config::SawtoothQuadConfig, container, C, names, time_steps, x_var, bounds, meta) Approximate x² using the sawtooth MIP formulation. @@ -50,8 +50,7 @@ with maximum overestimation error Δ² · 2^{-2L-2} where Δ = x_max - x_min. - `names::Vector{String}`: component names - `time_steps::UnitRange{Int}`: time periods - `x_var`: container of variables indexed by (name, t) -- `x_min::Float64`: lower bound of x domain -- `x_max::Float64`: upper bound of x domain +- `bounds::Vector{MinMax}`: per-name lower and upper bounds of x domain - `meta::String`: variable type identifier for the approximation (allows multiple approximations per component type) """ function _add_quadratic_approx!( @@ -61,14 +60,11 @@ function _add_quadratic_approx!( names::Vector{String}, time_steps::UnitRange{Int}, x_var, - x_min::Float64, - x_max::Float64, + bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} - IS.@assert_op x_max > x_min IS.@assert_op config.depth >= 1 jump_model = get_jump_model(container) - delta = x_max - x_min # Create containers with known dimensions g_levels = 0:(config.depth) @@ -122,7 +118,7 @@ function _add_quadratic_approx!( lp_expr = _add_quadratic_approx!( EpigraphQuadConfig(config.epigraph_depth), container, C, names, time_steps, - x_var, x_min, x_max, meta * "_lb", + x_var, bounds, meta * "_lb", ) z_var = add_variable_container!( container, @@ -143,14 +139,13 @@ function _add_quadratic_approx!( ) end - # Precompute sawtooth coefficients (invariant across names and time steps) - saw_coeffs = [delta * delta * (2.0^(-2 * j)) for j in alpha_levels] - - # Compute valid bounds for z ≈ x² from variable bounds - z_min = (x_min <= 0.0 <= x_max) ? 0.0 : min(x_min * x_min, x_max * x_max) - z_max = max(x_min * x_min, x_max * x_max) - - for name in names, t in time_steps + for (i, name) in enumerate(names), t in time_steps + b = bounds[i] + IS.@assert_op b.max > b.min + delta = b.max - b.min + saw_coeffs = [delta * delta * (2.0^(-2 * j)) for j in alpha_levels] + z_min = (b.min <= 0.0 <= b.max) ? 0.0 : min(b.min * b.min, b.max * b.max) + z_max = max(b.min * b.min, b.max * b.max) x = x_var[name, t] # Auxiliary variables g_0,...,g_L ∈ [0, 1] @@ -175,7 +170,7 @@ function _add_quadratic_approx!( # Linking constraint: g_0 = (x - x_min) / Δ link_cons[name, t] = JuMP.@constraint( jump_model, - g_var[name, 0, t] == (x - x_min) / delta, + g_var[name, 0, t] == (x - b.min) / delta, ) # S^L constraints for j = 1,...,L @@ -198,11 +193,11 @@ function _add_quadratic_approx!( end # Build x² ≈ x_min² + (2 x_min Δ + Δ²) g_0 - Σ_{j=1}^L Δ² 2^{-2j} g_j - x_sq_approx = JuMP.AffExpr(x_min * x_min) + x_sq_approx = JuMP.AffExpr(b.min * b.min) add_proportional_to_jump_expression!( x_sq_approx, g_var[name, 0, t], - 2.0 * x_min * delta + delta * delta, + 2.0 * b.min * delta + delta * delta, ) for j in alpha_levels add_proportional_to_jump_expression!( diff --git a/src/quadratic_approximations/solver_sos2.jl b/src/quadratic_approximations/solver_sos2.jl index aea0a521..32d34a0a 100644 --- a/src/quadratic_approximations/solver_sos2.jl +++ b/src/quadratic_approximations/solver_sos2.jl @@ -29,7 +29,7 @@ end SolverSOS2QuadConfig(depth::Int) = SolverSOS2QuadConfig(depth, 4) """ - _add_quadratic_approx!(config::SolverSOS2QuadConfig, container, C, names, time_steps, x_var, x_min, x_max, meta) + _add_quadratic_approx!(config::SolverSOS2QuadConfig, container, C, names, time_steps, x_var, bounds, meta) Approximate x² using a piecewise linear function with solver-native SOS2 constraints. @@ -44,8 +44,7 @@ approximating x² in a `QuadraticExpression` expression container. - `names::Vector{String}`: component names - `time_steps::UnitRange{Int}`: time periods - `x_var`: container of variables indexed by (name, t) -- `x_min::Float64`: lower bound of x domain -- `x_max::Float64`: upper bound of x domain +- `bounds::Vector{MinMax}`: per-name lower and upper bounds of x domain - `meta::String`: variable type identifier for the approximation (allows multiple approximations per component type) """ function _add_quadratic_approx!( @@ -55,11 +54,9 @@ function _add_quadratic_approx!( names::Vector{String}, time_steps::UnitRange{Int}, x_var, - x_min::Float64, - x_max::Float64, + bounds::Vector{MinMax}, meta::String, ) where {C <: IS.InfrastructureSystemsComponent} - lx = x_max - x_min x_bkpts, x_sq_bkpts = _get_breakpoints_for_pwl_function( 0.0, @@ -122,7 +119,9 @@ function _add_quadratic_approx!( meta, ) - for name in names, t in time_steps + for (i, name) in enumerate(names), t in time_steps + b = bounds[i] + lx = b.max - b.min x = x_var[name, t] # Create lambda_var variables: λ_i ∈ [0, 1] @@ -142,7 +141,7 @@ function _add_quadratic_approx!( for i in eachindex(x_bkpts) add_proportional_to_jump_expression!(link, lambda[i], x_bkpts[i]) end - link_cons[name, t] = JuMP.@constraint(jump_model, (x - x_min) / lx == link) + link_cons[name, t] = JuMP.@constraint(jump_model, (x - b.min) / lx == link) # Σ λ_i = 1 norm = norm_expr[name, t] = JuMP.AffExpr(0.0) @@ -162,15 +161,15 @@ function _add_quadratic_approx!( end x_sq = JuMP.AffExpr(0.0) add_proportional_to_jump_expression!(x_sq, x_hat_sq, lx * lx) - add_proportional_to_jump_expression!(x_sq, x, 2 * x_min) - add_constant_to_jump_expression!(x_sq, -x_min * x_min) + add_proportional_to_jump_expression!(x_sq, x, 2 * b.min) + add_constant_to_jump_expression!(x_sq, -b.min * b.min) result_expr[name, t] = x_sq end if config.pwmcc_segments > 0 _add_pwmcc_concave_cuts!( container, C, names, time_steps, - x_var, result_expr, x_min, x_max, + x_var, result_expr, bounds, config.pwmcc_segments, meta * "_pwmcc", ) end From 504d8f2e11b3704467e362fb414458bf4354a2c0 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Wed, 15 Apr 2026 18:24:24 -0400 Subject: [PATCH 02/17] converted axis-less containers in some approximations; fixed fallback for unsupported constraint in numerical analysis; added output handling methods for different variable arrays --- src/operation/decision_model_store.jl | 13 +++++++ .../model_numerical_analysis_utils.jl | 5 +-- src/quadratic_approximations/manual_sos2.jl | 25 ++++++++++--- src/quadratic_approximations/pwmcc_cuts.jl | 26 ++++++++++---- src/quadratic_approximations/solver_sos2.jl | 13 +++++-- src/utils/jump_utils.jl | 36 +++++++++++++++++++ 6 files changed, 102 insertions(+), 16 deletions(-) diff --git a/src/operation/decision_model_store.jl b/src/operation/decision_model_store.jl index 89bc5967..e9f65c64 100644 --- a/src/operation/decision_model_store.jl +++ b/src/operation/decision_model_store.jl @@ -112,6 +112,19 @@ function write_output!( return end +function write_output!( + store::DecisionModelStore, + name::Symbol, + key::OptimizationContainerKey, + index::DecisionModelIndexType, + update_timestamp::Dates.DateTime, + array::DenseAxisArray{T, 3, <:Tuple{Vector{String}, UnitRange, UnitRange}}, +) where {T} + container = getfield(store, get_store_container_type(key)) + container[key][index] = array + return +end + function read_outputs( store::DecisionModelStore, key::OptimizationContainerKey; diff --git a/src/operation/model_numerical_analysis_utils.jl b/src/operation/model_numerical_analysis_utils.jl index 0cffa994..a023d49e 100644 --- a/src/operation/model_numerical_analysis_utils.jl +++ b/src/operation/model_numerical_analysis_utils.jl @@ -98,8 +98,9 @@ function update_numerical_bounds(bonuds::NumericalBounds, func::MOI.Interval, id return update_numerical_bounds(bonuds, func.lower, idx) end -# Default fallback for unsupported constraints. -update_numerical_bounds(::NumericalBounds, func, idx) = nothing +# Default fallbacks for unsupported constraints. +update_coefficient_bounds(::ConstraintBounds, func, idx) = nothing +update_rhs_bounds(::ConstraintBounds, func, idx) = nothing function get_constraint_numerical_bounds(model::OperationModel) if !is_built(model) diff --git a/src/quadratic_approximations/manual_sos2.jl b/src/quadratic_approximations/manual_sos2.jl index d061708a..bd63e156 100644 --- a/src/quadratic_approximations/manual_sos2.jl +++ b/src/quadratic_approximations/manual_sos2.jl @@ -65,9 +65,24 @@ function _add_quadratic_approx!( jump_model = get_jump_model(container) # Create all containers upfront - lambda_container = - add_variable_container!(container, QuadraticVariable, C; meta) - z_container = add_variable_container!(container, ManualSOS2BinaryVariable, C; meta) + lambda_container = add_variable_container!( + container, + QuadraticVariable, + C, + names, + 1:n_points, + time_steps; + meta + ) + z_container = add_variable_container!( + container, + ManualSOS2BinaryVariable, + C, + names, + 1:n_bins, + time_steps; + meta + ) link_cons = add_constraints_container!( container, SOS2LinkingConstraint, @@ -143,7 +158,7 @@ function _add_quadratic_approx!( lambda = Vector{JuMP.VariableRef}(undef, n_points) for i in 1:n_points lambda[i] = - lambda_container[(name, i, t)] = JuMP.@variable( + lambda_container[name, i, t] = JuMP.@variable( jump_model, base_name = "QuadraticVariable_$(C)_{$(name), pwl_$(i), $(t)}", lower_bound = 0.0, @@ -169,7 +184,7 @@ function _add_quadratic_approx!( z_vars = Vector{JuMP.VariableRef}(undef, n_bins) for j in 1:n_bins z_vars[j] = - z_container[(name, j, t)] = JuMP.@variable( + z_container[name, j, t] = JuMP.@variable( jump_model, base_name = "ManualSOS2Binary_$(C)_{$(name), $(j), $(t)}", binary = true, diff --git a/src/quadratic_approximations/pwmcc_cuts.jl b/src/quadratic_approximations/pwmcc_cuts.jl index 3643bc39..58aaba58 100644 --- a/src/quadratic_approximations/pwmcc_cuts.jl +++ b/src/quadratic_approximations/pwmcc_cuts.jl @@ -67,10 +67,24 @@ function _add_pwmcc_concave_cuts!( jump_model = get_jump_model(container) # Create containers - delta_var = add_variable_container!(container, PiecewiseMcCormickBinary, C; meta) - vd_var = - add_variable_container!(container, PiecewiseMcCormickDisaggregated, C; meta) - + delta_var = add_variable_container!( + container, + PiecewiseMcCormickBinary, + C, + names, + 1:K, + time_steps; + meta + ) + vd_var = add_variable_container!( + container, + PiecewiseMcCormickDisaggregated, + C, + names, + 1:K, + time_steps; + meta + ) selector_cons = add_constraints_container!( container, PiecewiseMcCormickSelectorSum(), @@ -153,13 +167,13 @@ function _add_pwmcc_concave_cuts!( for k in 1:K delta[k] = - delta_var[(name, k, t)] = JuMP.@variable( + delta_var[name, k, t] = JuMP.@variable( jump_model, base_name = "PwMcCBin_$(C)_{$(name), $(k), $(t)}", binary = true, ) vd[k] = - vd_var[(name, k, t)] = JuMP.@variable( + vd_var[name, k, t] = JuMP.@variable( jump_model, base_name = "PwMcCDis_$(C)_{$(name), $(k), $(t)}", ) diff --git a/src/quadratic_approximations/solver_sos2.jl b/src/quadratic_approximations/solver_sos2.jl index 32d34a0a..1efc7047 100644 --- a/src/quadratic_approximations/solver_sos2.jl +++ b/src/quadratic_approximations/solver_sos2.jl @@ -68,8 +68,15 @@ function _add_quadratic_approx!( jump_model = get_jump_model(container) # Create all containers upfront - lambda_var = - add_variable_container!(container, QuadraticVariable, C; meta) + lambda_var = add_variable_container!( + container, + QuadraticVariable, + C, + names, + 1:n_points, + time_steps; + meta + ) link_cons = add_constraints_container!( container, SOS2LinkingConstraint, @@ -128,7 +135,7 @@ function _add_quadratic_approx!( lambda = Vector{JuMP.VariableRef}(undef, n_points) for i in 1:n_points lambda[i] = - lambda_var[(name, i, t)] = JuMP.@variable( + lambda_var[name, i, t] = JuMP.@variable( jump_model, base_name = "QuadraticVariable_$(C)_{$(name), pwl_$(i), $(t)}", lower_bound = 0.0, diff --git a/src/utils/jump_utils.jl b/src/utils/jump_utils.jl index d8f2d1c5..f2fbc26a 100644 --- a/src/utils/jump_utils.jl +++ b/src/utils/jump_utils.jl @@ -428,6 +428,42 @@ function to_outputs_dataframe( ) end +function to_outputs_dataframe( + array::DenseAxisArray{ + Float64, + 3, + <:Tuple{Vector{String}, UnitRange{Int}, UnitRange{Int}}, + }, + ::Nothing, + ::Val{TableFormat.LONG}, +) + num_rows = length(array.data) + time_col = Vector{Int}(undef, num_rows) + name_col = Vector{String}(undef, num_rows) + name2_col = Vector{Int}(undef, num_rows) + vals = Vector{Float64}(undef, num_rows) + + row_index = 1 + for name in axes(array, 1) + for name2 in axes(array, 2) + for time_index in axes(array, 3) + time_col[row_index] = time_index + name_col[row_index] = name + name2_col[row_index] = name2 + vals[row_index] = array[name, name2, time_index] + row_index += 1 + end + end + end + + return DataFrame( + :time_index => time_col, + :name => name_col, + :name2 => name2_col, + :value => vals, + ) +end + function to_dataframe(array::SparseAxisArray{T, N, K}) where {T, N, K <: NTuple{N, Any}} columns = get_column_names_from_axis_array(array) return DataFrames.DataFrame(_to_matrix(array, columns), columns) From 3b9f1b0ded0ab6da8285cf633f1800f844a0d990 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli <8549957+acostarelli@users.noreply.github.com> Date: Tue, 21 Apr 2026 17:23:05 -0400 Subject: [PATCH 03/17] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/bilinear_approximations/bin2.jl | 7 ++++++- src/bilinear_approximations/hybs.jl | 14 ++++++++++++-- src/quadratic_approximations/manual_sos2.jl | 4 ++-- src/quadratic_approximations/pwmcc_cuts.jl | 4 ++-- src/quadratic_approximations/solver_sos2.jl | 2 +- 5 files changed, 23 insertions(+), 8 deletions(-) diff --git a/src/bilinear_approximations/bin2.jl b/src/bilinear_approximations/bin2.jl index 5b77261b..5c331cd4 100644 --- a/src/bilinear_approximations/bin2.jl +++ b/src/bilinear_approximations/bin2.jl @@ -96,7 +96,12 @@ function _add_bilinear_approx!( # --- Bin2 identity: z = ½((x+y)² − x² − y²) --- # Bounds for p = x + y (per-name) - p_bounds = [MinMax((min = x_bounds[i].min + y_bounds[i].min, max = x_bounds[i].max + y_bounds[i].max)) for i in eachindex(x_bounds)] + p_bounds = [ + MinMax(( + min = x_bounds[i].min + y_bounds[i].min, + max = x_bounds[i].max + y_bounds[i].max, + )) for i in eachindex(x_bounds) + ] meta_plus = meta * "_plus" diff --git a/src/bilinear_approximations/hybs.jl b/src/bilinear_approximations/hybs.jl index 3ab8da61..c705437c 100644 --- a/src/bilinear_approximations/hybs.jl +++ b/src/bilinear_approximations/hybs.jl @@ -93,8 +93,18 @@ function _add_bilinear_approx!( meta::String, ) where {C <: IS.InfrastructureSystemsComponent} # Bounds for auxiliary variables (per-name) - p1_bounds = [MinMax((min = x_bounds[i].min + y_bounds[i].min, max = x_bounds[i].max + y_bounds[i].max)) for i in eachindex(x_bounds)] - p2_bounds = [MinMax((min = x_bounds[i].min - y_bounds[i].max, max = x_bounds[i].max - y_bounds[i].min)) for i in eachindex(x_bounds)] + p1_bounds = [ + MinMax(( + min = x_bounds[i].min + y_bounds[i].min, + max = x_bounds[i].max + y_bounds[i].max, + )) for i in eachindex(x_bounds) + ] + p2_bounds = [ + MinMax(( + min = x_bounds[i].min - y_bounds[i].max, + max = x_bounds[i].max - y_bounds[i].min, + )) for i in eachindex(x_bounds) + ] jump_model = get_jump_model(container) diff --git a/src/quadratic_approximations/manual_sos2.jl b/src/quadratic_approximations/manual_sos2.jl index bd63e156..49f210ad 100644 --- a/src/quadratic_approximations/manual_sos2.jl +++ b/src/quadratic_approximations/manual_sos2.jl @@ -72,7 +72,7 @@ function _add_quadratic_approx!( names, 1:n_points, time_steps; - meta + meta, ) z_container = add_variable_container!( container, @@ -81,7 +81,7 @@ function _add_quadratic_approx!( names, 1:n_bins, time_steps; - meta + meta, ) link_cons = add_constraints_container!( container, diff --git a/src/quadratic_approximations/pwmcc_cuts.jl b/src/quadratic_approximations/pwmcc_cuts.jl index 58aaba58..0de17db9 100644 --- a/src/quadratic_approximations/pwmcc_cuts.jl +++ b/src/quadratic_approximations/pwmcc_cuts.jl @@ -74,7 +74,7 @@ function _add_pwmcc_concave_cuts!( names, 1:K, time_steps; - meta + meta, ) vd_var = add_variable_container!( container, @@ -83,7 +83,7 @@ function _add_pwmcc_concave_cuts!( names, 1:K, time_steps; - meta + meta, ) selector_cons = add_constraints_container!( container, diff --git a/src/quadratic_approximations/solver_sos2.jl b/src/quadratic_approximations/solver_sos2.jl index 1efc7047..02cc3911 100644 --- a/src/quadratic_approximations/solver_sos2.jl +++ b/src/quadratic_approximations/solver_sos2.jl @@ -75,7 +75,7 @@ function _add_quadratic_approx!( names, 1:n_points, time_steps; - meta + meta, ) link_cons = add_constraints_container!( container, From ce70c8bae60748ab08363e97b601f906669b83e7 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Tue, 21 Apr 2026 17:31:40 -0400 Subject: [PATCH 04/17] fix type dispatch --- src/quadratic_approximations/pwmcc_cuts.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/quadratic_approximations/pwmcc_cuts.jl b/src/quadratic_approximations/pwmcc_cuts.jl index 0de17db9..31819aa8 100644 --- a/src/quadratic_approximations/pwmcc_cuts.jl +++ b/src/quadratic_approximations/pwmcc_cuts.jl @@ -95,7 +95,7 @@ function _add_pwmcc_concave_cuts!( ) linking_cons = add_constraints_container!( container, - PiecewiseMcCormickLinking(), + PiecewiseMcCormickLinking, C, names, time_steps; @@ -103,7 +103,7 @@ function _add_pwmcc_concave_cuts!( ) interval_lb_cons = add_constraints_container!( container, - PiecewiseMcCormickIntervalLB(), + PiecewiseMcCormickIntervalLB, C, names, 1:K, @@ -112,7 +112,7 @@ function _add_pwmcc_concave_cuts!( ) interval_ub_cons = add_constraints_container!( container, - PiecewiseMcCormickIntervalUB(), + PiecewiseMcCormickIntervalUB, C, names, 1:K, @@ -121,7 +121,7 @@ function _add_pwmcc_concave_cuts!( ) chord_ub_cons = add_constraints_container!( container, - PiecewiseMcCormickChordUB(), + PiecewiseMcCormickChordUB, C, names, time_steps; @@ -129,7 +129,7 @@ function _add_pwmcc_concave_cuts!( ) tangent_lb_l_cons = add_constraints_container!( container, - PiecewiseMcCormickTangentLBL(), + PiecewiseMcCormickTangentLBL, C, names, time_steps; @@ -137,7 +137,7 @@ function _add_pwmcc_concave_cuts!( ) tangent_lb_r_cons = add_constraints_container!( container, - PiecewiseMcCormickTangentLBR(), + PiecewiseMcCormickTangentLBR, C, names, time_steps; From d1886169ab61592098d355fe3815787fd85296ea Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Tue, 21 Apr 2026 17:38:04 -0400 Subject: [PATCH 05/17] more type dispatch fix --- src/quadratic_approximations/pwmcc_cuts.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/quadratic_approximations/pwmcc_cuts.jl b/src/quadratic_approximations/pwmcc_cuts.jl index 31819aa8..a636e6f8 100644 --- a/src/quadratic_approximations/pwmcc_cuts.jl +++ b/src/quadratic_approximations/pwmcc_cuts.jl @@ -87,7 +87,7 @@ function _add_pwmcc_concave_cuts!( ) selector_cons = add_constraints_container!( container, - PiecewiseMcCormickSelectorSum(), + PiecewiseMcCormickSelectorSum, C, names, time_steps; From d08afe5e49f8817a63b315a44d2b07fa69de09d3 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Wed, 22 Apr 2026 11:54:13 -0400 Subject: [PATCH 06/17] copilot review --- src/bilinear_approximations/hybs.jl | 2 + src/bilinear_approximations/mccormick.jl | 2 +- src/quadratic_approximations/epigraph.jl | 2 +- src/quadratic_approximations/manual_sos2.jl | 1 + src/quadratic_approximations/nmdt_common.jl | 2 + src/quadratic_approximations/sawtooth.jl | 2 +- src/quadratic_approximations/solver_sos2.jl | 1 + src/utils/jump_utils.jl | 42 ++------------------- test/test_jump_utils.jl | 20 ++++++++++ 9 files changed, 32 insertions(+), 42 deletions(-) diff --git a/src/bilinear_approximations/hybs.jl b/src/bilinear_approximations/hybs.jl index c705437c..2a6fe606 100644 --- a/src/bilinear_approximations/hybs.jl +++ b/src/bilinear_approximations/hybs.jl @@ -188,6 +188,8 @@ function _add_bilinear_approx!( for (i, name) in enumerate(names), t in time_steps xb = x_bounds[i] yb = y_bounds[i] + IS.@assert_op xb.max > xb.min "Invalid bounds for $(name): expected max > min, got min=$(xb.min), max=$(xb.max)" + IS.@assert_op yb.max > yb.min "Invalid bounds for $(name): expected max > min, got min=$(yb.min), max=$(yb.max)" # Compute valid bounds for z ≈ x·y from variable bounds z_lo = min(xb.min * yb.min, xb.min * yb.max, xb.max * yb.min, xb.max * yb.max) diff --git a/src/bilinear_approximations/mccormick.jl b/src/bilinear_approximations/mccormick.jl index e86fba0e..8ee2e3d2 100644 --- a/src/bilinear_approximations/mccormick.jl +++ b/src/bilinear_approximations/mccormick.jl @@ -212,7 +212,7 @@ function _add_mccormick_envelope!( mc_cons = add_constraints_container!( container, - McCormickConstraint(), + McCormickConstraint, C, names, 1:4, diff --git a/src/quadratic_approximations/epigraph.jl b/src/quadratic_approximations/epigraph.jl index 2b180b6c..0563debe 100644 --- a/src/quadratic_approximations/epigraph.jl +++ b/src/quadratic_approximations/epigraph.jl @@ -124,7 +124,7 @@ function _add_quadratic_approx!( for (i, name) in enumerate(names), t in time_steps b = bounds[i] - IS.@assert_op b.max > b.min + IS.@assert_op b.max > b.min "Invalid bounds for $(name): expected max > min, got min=$(b.min), max=$(b.max)" delta = b.max - b.min z_ub = max(b.min^2, b.max^2) x = x_var[name, t] diff --git a/src/quadratic_approximations/manual_sos2.jl b/src/quadratic_approximations/manual_sos2.jl index 49f210ad..10309337 100644 --- a/src/quadratic_approximations/manual_sos2.jl +++ b/src/quadratic_approximations/manual_sos2.jl @@ -151,6 +151,7 @@ function _add_quadratic_approx!( for (i, name) in enumerate(names), t in time_steps b = bounds[i] + IS.@assert_op b.max > b.min "Invalid bounds for $(name): expected max > min, got min=$(b.min), max=$(b.max)" lx = b.max - b.min x = x_var[name, t] diff --git a/src/quadratic_approximations/nmdt_common.jl b/src/quadratic_approximations/nmdt_common.jl index 3d589ec0..22113f35 100644 --- a/src/quadratic_approximations/nmdt_common.jl +++ b/src/quadratic_approximations/nmdt_common.jl @@ -392,6 +392,8 @@ function _assemble_product!( for (i, name) in enumerate(names), t in time_steps xb = x_bounds[i] yb = y_bounds[i] + IS.@assert_op xb.max > xb.min "Invalid bounds for $(name): expected max > min, got min=$(xb.min), max=$(xb.max)" + IS.@assert_op yb.max > yb.min "Invalid bounds for $(name): expected max > min, got min=$(yb.min), max=$(yb.max)" lx = xb.max - xb.min ly = yb.max - yb.min diff --git a/src/quadratic_approximations/sawtooth.jl b/src/quadratic_approximations/sawtooth.jl index 0a0d0d2e..4efc3841 100644 --- a/src/quadratic_approximations/sawtooth.jl +++ b/src/quadratic_approximations/sawtooth.jl @@ -141,7 +141,7 @@ function _add_quadratic_approx!( for (i, name) in enumerate(names), t in time_steps b = bounds[i] - IS.@assert_op b.max > b.min + IS.@assert_op b.max > b.min "Invalid bounds for $(name): expected max > min, got min=$(b.min), max=$(b.max)" delta = b.max - b.min saw_coeffs = [delta * delta * (2.0^(-2 * j)) for j in alpha_levels] z_min = (b.min <= 0.0 <= b.max) ? 0.0 : min(b.min * b.min, b.max * b.max) diff --git a/src/quadratic_approximations/solver_sos2.jl b/src/quadratic_approximations/solver_sos2.jl index 02cc3911..86a2b73a 100644 --- a/src/quadratic_approximations/solver_sos2.jl +++ b/src/quadratic_approximations/solver_sos2.jl @@ -128,6 +128,7 @@ function _add_quadratic_approx!( for (i, name) in enumerate(names), t in time_steps b = bounds[i] + IS.@assert_op b.max > b.min "Invalid bounds for $(name): expected max > min, got min=$(b.min), max=$(b.max)" lx = b.max - b.min x = x_var[name, t] diff --git a/src/utils/jump_utils.jl b/src/utils/jump_utils.jl index f2fbc26a..80865062 100644 --- a/src/utils/jump_utils.jl +++ b/src/utils/jump_utils.jl @@ -396,14 +396,14 @@ function to_outputs_dataframe( array::DenseAxisArray{ Float64, 3, - <:Tuple{Vector{String}, Vector{String}, UnitRange{Int}}, + <:Tuple{Vector{String}, T, UnitRange{Int}}, }, ::Nothing, ::Val{TableFormat.LONG}, -) +) where {T <: Union{Vector{String}, UnitRange{Int}}} num_rows = length(array.data) time_col = Vector{Int}(undef, num_rows) - name_col = Vector{String}(undef, num_rows) + name_col = Vector{eltype(T)}(undef, num_rows) name2_col = Vector{String}(undef, num_rows) vals = Vector{Float64}(undef, num_rows) @@ -428,42 +428,6 @@ function to_outputs_dataframe( ) end -function to_outputs_dataframe( - array::DenseAxisArray{ - Float64, - 3, - <:Tuple{Vector{String}, UnitRange{Int}, UnitRange{Int}}, - }, - ::Nothing, - ::Val{TableFormat.LONG}, -) - num_rows = length(array.data) - time_col = Vector{Int}(undef, num_rows) - name_col = Vector{String}(undef, num_rows) - name2_col = Vector{Int}(undef, num_rows) - vals = Vector{Float64}(undef, num_rows) - - row_index = 1 - for name in axes(array, 1) - for name2 in axes(array, 2) - for time_index in axes(array, 3) - time_col[row_index] = time_index - name_col[row_index] = name - name2_col[row_index] = name2 - vals[row_index] = array[name, name2, time_index] - row_index += 1 - end - end - end - - return DataFrame( - :time_index => time_col, - :name => name_col, - :name2 => name2_col, - :value => vals, - ) -end - function to_dataframe(array::SparseAxisArray{T, N, K}) where {T, N, K <: NTuple{N, Any}} columns = get_column_names_from_axis_array(array) return DataFrames.DataFrame(_to_matrix(array, columns), columns) diff --git a/test/test_jump_utils.jl b/test/test_jump_utils.jl index 1085e37e..0d1723d6 100644 --- a/test/test_jump_utils.jl +++ b/test/test_jump_utils.jl @@ -266,4 +266,24 @@ struct MockVariableType <: ISOPT.VariableType end @test :value in propertynames(df) @test nrow(df) == 12 # 2 * 2 * 3 end + + @testset "to_outputs_dataframe - 3D LONG format (String, Int, Int)" begin + data = DenseAxisArray( + zeros(2, 2, 3), + ["gen1", "gen2"], + ["area1", "area2"], + 1:3, + ) + data["gen1", "area1", 1] = 1.0 + data["gen2", "area2", 2] = 2.0 + + timestamps = [DateTime(2024, 1, 1, i) for i in 0:2] + df = IOM.to_outputs_dataframe(data, timestamps, Val(IOM.TableFormat.LONG)) + @test isa(df, DataFrame) + @test :DateTime in propertynames(df) + @test :name in propertynames(df) + @test :name2 in propertynames(df) + @test :value in propertynames(df) + @test nrow(df) == 12 # 2 * 2 * 3 + end end From 442b7dac056a65f7121fc6ddde91160ed306b4ae Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Wed, 22 Apr 2026 12:08:23 -0400 Subject: [PATCH 07/17] fix assert_op class; remove zz tests (again?) --- src/bilinear_approximations/hybs.jl | 4 +- src/quadratic_approximations/epigraph.jl | 2 +- src/quadratic_approximations/manual_sos2.jl | 2 +- src/quadratic_approximations/nmdt_common.jl | 4 +- src/quadratic_approximations/sawtooth.jl | 2 +- src/quadratic_approximations/solver_sos2.jl | 2 +- test/runtests.jl | 8 +- test/test_zzb_approximations.jl | 485 -------------------- test/test_zzi_bilinear.jl | 466 ------------------- 9 files changed, 12 insertions(+), 963 deletions(-) delete mode 100644 test/test_zzb_approximations.jl delete mode 100644 test/test_zzi_bilinear.jl diff --git a/src/bilinear_approximations/hybs.jl b/src/bilinear_approximations/hybs.jl index 2a6fe606..8df74bff 100644 --- a/src/bilinear_approximations/hybs.jl +++ b/src/bilinear_approximations/hybs.jl @@ -188,8 +188,8 @@ function _add_bilinear_approx!( for (i, name) in enumerate(names), t in time_steps xb = x_bounds[i] yb = y_bounds[i] - IS.@assert_op xb.max > xb.min "Invalid bounds for $(name): expected max > min, got min=$(xb.min), max=$(xb.max)" - IS.@assert_op yb.max > yb.min "Invalid bounds for $(name): expected max > min, got min=$(yb.min), max=$(yb.max)" + IS.@assert_op xb.max > xb.min + IS.@assert_op yb.max > yb.min # Compute valid bounds for z ≈ x·y from variable bounds z_lo = min(xb.min * yb.min, xb.min * yb.max, xb.max * yb.min, xb.max * yb.max) diff --git a/src/quadratic_approximations/epigraph.jl b/src/quadratic_approximations/epigraph.jl index 0563debe..2b180b6c 100644 --- a/src/quadratic_approximations/epigraph.jl +++ b/src/quadratic_approximations/epigraph.jl @@ -124,7 +124,7 @@ function _add_quadratic_approx!( for (i, name) in enumerate(names), t in time_steps b = bounds[i] - IS.@assert_op b.max > b.min "Invalid bounds for $(name): expected max > min, got min=$(b.min), max=$(b.max)" + IS.@assert_op b.max > b.min delta = b.max - b.min z_ub = max(b.min^2, b.max^2) x = x_var[name, t] diff --git a/src/quadratic_approximations/manual_sos2.jl b/src/quadratic_approximations/manual_sos2.jl index 10309337..0275424a 100644 --- a/src/quadratic_approximations/manual_sos2.jl +++ b/src/quadratic_approximations/manual_sos2.jl @@ -151,7 +151,7 @@ function _add_quadratic_approx!( for (i, name) in enumerate(names), t in time_steps b = bounds[i] - IS.@assert_op b.max > b.min "Invalid bounds for $(name): expected max > min, got min=$(b.min), max=$(b.max)" + IS.@assert_op b.max > b.min lx = b.max - b.min x = x_var[name, t] diff --git a/src/quadratic_approximations/nmdt_common.jl b/src/quadratic_approximations/nmdt_common.jl index 22113f35..35ae8248 100644 --- a/src/quadratic_approximations/nmdt_common.jl +++ b/src/quadratic_approximations/nmdt_common.jl @@ -392,8 +392,8 @@ function _assemble_product!( for (i, name) in enumerate(names), t in time_steps xb = x_bounds[i] yb = y_bounds[i] - IS.@assert_op xb.max > xb.min "Invalid bounds for $(name): expected max > min, got min=$(xb.min), max=$(xb.max)" - IS.@assert_op yb.max > yb.min "Invalid bounds for $(name): expected max > min, got min=$(yb.min), max=$(yb.max)" + IS.@assert_op xb.max > xb.min + IS.@assert_op yb.max > yb.min lx = xb.max - xb.min ly = yb.max - yb.min diff --git a/src/quadratic_approximations/sawtooth.jl b/src/quadratic_approximations/sawtooth.jl index 4efc3841..0a0d0d2e 100644 --- a/src/quadratic_approximations/sawtooth.jl +++ b/src/quadratic_approximations/sawtooth.jl @@ -141,7 +141,7 @@ function _add_quadratic_approx!( for (i, name) in enumerate(names), t in time_steps b = bounds[i] - IS.@assert_op b.max > b.min "Invalid bounds for $(name): expected max > min, got min=$(b.min), max=$(b.max)" + IS.@assert_op b.max > b.min delta = b.max - b.min saw_coeffs = [delta * delta * (2.0^(-2 * j)) for j in alpha_levels] z_min = (b.min <= 0.0 <= b.max) ? 0.0 : min(b.min * b.min, b.max * b.max) diff --git a/src/quadratic_approximations/solver_sos2.jl b/src/quadratic_approximations/solver_sos2.jl index 86a2b73a..110428fe 100644 --- a/src/quadratic_approximations/solver_sos2.jl +++ b/src/quadratic_approximations/solver_sos2.jl @@ -128,7 +128,7 @@ function _add_quadratic_approx!( for (i, name) in enumerate(names), t in time_steps b = bounds[i] - IS.@assert_op b.max > b.min "Invalid bounds for $(name): expected max > min, got min=$(b.min), max=$(b.max)" + IS.@assert_op b.max > b.min lx = b.max - b.min x = x_var[name, t] diff --git a/test/runtests.jl b/test/runtests.jl index 91774c12..c0e37618 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,10 +9,10 @@ using InfrastructureSystems const IS = InfrastructureSystems # Code Quality Tests -import Aqua -@testset "Code Quality (Aqua.jl)" begin - Aqua.test_all(InfrastructureOptimizationModels; persistent_tasks = false) -end +# import Aqua +# @testset "Code Quality (Aqua.jl)" begin +# Aqua.test_all(InfrastructureOptimizationModels; persistent_tasks = false) +# end # Load the test module include("InfrastructureOptimizationModelsTests.jl") diff --git a/test/test_zzb_approximations.jl b/test/test_zzb_approximations.jl deleted file mode 100644 index 2ebd2ba0..00000000 --- a/test/test_zzb_approximations.jl +++ /dev/null @@ -1,485 +0,0 @@ -const ZZB_META = "ZZBTest" - -@testset "ZZB Quadratic Approximation" begin - - # ==================================================================== - # 1. Pure encoding tests (no JuMP) - # ==================================================================== - @testset "BRGC encoding" begin - @testset "r=1" begin - G = IOM.build_brgc(1) - @test size(G) == (2, 1) - @test G == [0; 1][:, :] - end - - @testset "r=2" begin - G = IOM.build_brgc(2) - @test size(G) == (4, 2) - @test G == [0 0; 0 1; 1 1; 1 0] - end - - @testset "r=3" begin - G = IOM.build_brgc(3) - @test size(G) == (8, 3) - # Verify adjacency property: consecutive rows differ by exactly 1 bit - for i in 1:(size(G, 1) - 1) - diff = sum(abs.(G[i + 1, :] .- G[i, :])) - @test diff == 1 - end - end - - @testset "Adjacency holds for r=1..5" begin - for r in 1:5 - G = IOM.build_brgc(r) - @test size(G) == (2^r, r) - for i in 1:(size(G, 1) - 1) - diff = sum(abs.(G[i + 1, :] .- G[i, :])) - @test diff == 1 - end - end - end - end - - @testset "ZZB coefficient matrices" begin - @testset "r=2 known values" begin - G = IOM.build_brgc(2) - lower_coeffs, upper_coeffs = IOM._build_zzb_coefficients(G, 2) - # G = [0 0; 0 1; 1 1; 1 0], 5 breakpoints, 2 columns - @test size(lower_coeffs) == (5, 2) - @test size(upper_coeffs) == (5, 2) - - # First breakpoint (i=1): adjacent to segment 1 only → G[1,:] = [0,0] - @test lower_coeffs[1, :] == [0, 0] - @test upper_coeffs[1, :] == [0, 0] - - # Last breakpoint (i=5): adjacent to segment 4 only → G[4,:] = [1,0] - @test lower_coeffs[5, :] == [1, 0] - @test upper_coeffs[5, :] == [1, 0] - - # Interior breakpoint i=2: adjacent to segments 1,2 → G[1,:]=[0,0], G[2,:]=[0,1] - @test lower_coeffs[2, :] == [0, 0] # min(0,0), min(0,1) - @test upper_coeffs[2, :] == [0, 1] # max(0,0), max(0,1) - - # Interior breakpoint i=3: adjacent to segments 2,3 → G[2,:]=[0,1], G[3,:]=[1,1] - @test lower_coeffs[3, :] == [0, 1] # min(0,1), min(1,1) - @test upper_coeffs[3, :] == [1, 1] # max(0,1), max(1,1) - - # Interior breakpoint i=4: adjacent to segments 3,4 → G[3,:]=[1,1], G[4,:]=[1,0] - @test lower_coeffs[4, :] == [1, 0] # min(1,1), min(1,0) - @test upper_coeffs[4, :] == [1, 1] # max(1,1), max(1,0) - end - - @testset "Coefficient bounds are valid" begin - for r in 1:4 - G = IOM.build_brgc(r) - lo, hi = IOM._build_zzb_coefficients(G, r) - # All entries must be 0 or 1 - @test all(x -> x == 0 || x == 1, lo) - @test all(x -> x == 0 || x == 1, hi) - # lo ≤ hi elementwise - @test all(lo .<= hi) - end - end - end - - # ==================================================================== - # 2. ZZB quadratic structure test - # ==================================================================== - @testset "Constraint structure" begin - setup = _setup_qa_test(["dev1"], 1:1) - num_segments = 4 # 4 segments → log2(4) = 2 binary variables - - IOM._add_quadratic_approx!( - IOM.ZZBQuadConfig(num_segments), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.var_container, - 0.0, - 4.0, - ZZB_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.QuadraticExpression, - MockThermalGen, - ZZB_META, - ) - - @test expr_container["dev1", 1] isa JuMP.AffExpr - - # log2(num_segments) binary variables - n_bin = count(JuMP.is_binary, JuMP.all_variables(setup.jump_model)) - @test n_bin == Int(log2(num_segments)) - - # Lambda variables exist: num_segments + 1 = 5 lambda vars - @test IOM.has_container_key( - setup.container, - IOM.ZZBLambdaVariable, - MockThermalGen, - ZZB_META, - ) - - # Encoding constraints exist - @test IOM.has_container_key( - setup.container, - IOM.ZZBEncodingConstraint, - MockThermalGen, - ZZB_META, - ) - end - - # ==================================================================== - # 3. ZZB quadratic correctness (solve tests) - # ==================================================================== - @testset "Upper-bounds x^2 on [0,1]" begin - setup = _setup_qa_test(["dev1"], 1:1) - x_var = setup.var_container["dev1", 1] - JuMP.fix(x_var, 0.35; force = true) - - IOM._add_quadratic_approx!( - IOM.ZZBQuadConfig(8), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.var_container, - 0.0, - 1.0, - ZZB_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.QuadraticExpression, - MockThermalGen, - ZZB_META, - ) - z_expr = expr_container["dev1", 1] - - # Minimize: should give value ≤ x² (it's an upper bound, min gives closest) - JuMP.@objective(setup.jump_model, Min, z_expr) - JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) - JuMP.set_silent(setup.jump_model) - JuMP.optimize!(setup.jump_model) - - @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL - obj_val = JuMP.objective_value(setup.jump_model) - # ZZB is a convex combination → upper bound, so min ≥ true - small error - @test obj_val >= 0.35^2 - 1e-6 - # But not too far above - @test obj_val <= 0.35^2 + 0.05 - end - - @testset "Upper-bounds x^2 on non-unit interval [0,2]" begin - setup = _setup_qa_test(["dev1"], 1:1) - x_var = setup.var_container["dev1", 1] - JuMP.fix(x_var, 1.3; force = true) - - IOM._add_quadratic_approx!( - IOM.ZZBQuadConfig(8), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.var_container, - 0.0, - 2.0, - ZZB_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.QuadraticExpression, - MockThermalGen, - ZZB_META, - ) - z_expr = expr_container["dev1", 1] - - JuMP.@objective(setup.jump_model, Min, z_expr) - JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) - JuMP.set_silent(setup.jump_model) - JuMP.optimize!(setup.jump_model) - - @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL - obj_val = JuMP.objective_value(setup.jump_model) - @test obj_val >= 1.3^2 - 1e-6 - @test obj_val <= 1.3^2 + 0.2 - end - - @testset "Tightened ZZB bounds x^2" begin - x0 = 0.35 - true_val = x0^2 - num_segments = 4 # 4 segments → 2 binary variables, epigraph_depth=2 - - # Minimize tightened: should still be ≤ x² - setup_min = _setup_qa_test(["dev1"], 1:1) - JuMP.fix(setup_min.var_container["dev1", 1], x0; force = true) - IOM._add_quadratic_approx!( - IOM.ZZBQuadConfig(num_segments, Int(log2(num_segments))), - setup_min.container, - MockThermalGen, - ["dev1"], - 1:1, - setup_min.var_container, - 0.0, - 1.0, - ZZB_META, - ) - expr_min = IOM.get_expression( - setup_min.container, - IOM.QuadraticExpression, - MockThermalGen, - ZZB_META, - ) - z_min_expr = expr_min["dev1", 1] - - JuMP.@objective(setup_min.jump_model, Min, z_min_expr) - JuMP.set_optimizer(setup_min.jump_model, HiGHS.Optimizer) - JuMP.set_silent(setup_min.jump_model) - JuMP.optimize!(setup_min.jump_model) - @test JuMP.termination_status(setup_min.jump_model) == JuMP.OPTIMAL - min_val = JuMP.objective_value(setup_min.jump_model) - @test min_val <= true_val + 1e-6 - - # Maximize tightened: should be ≥ x² - setup_max = _setup_qa_test(["dev1"], 1:1) - JuMP.fix(setup_max.var_container["dev1", 1], x0; force = true) - IOM._add_quadratic_approx!( - IOM.ZZBQuadConfig(num_segments, Int(log2(num_segments))), - setup_max.container, - MockThermalGen, - ["dev1"], - 1:1, - setup_max.var_container, - 0.0, - 1.0, - ZZB_META, - ) - expr_max = IOM.get_expression( - setup_max.container, - IOM.QuadraticExpression, - MockThermalGen, - ZZB_META, - ) - z_max_expr = expr_max["dev1", 1] - - JuMP.@objective(setup_max.jump_model, Max, z_max_expr) - JuMP.set_optimizer(setup_max.jump_model, HiGHS.Optimizer) - JuMP.set_silent(setup_max.jump_model) - JuMP.optimize!(setup_max.jump_model) - @test JuMP.termination_status(setup_max.jump_model) == JuMP.OPTIMAL - max_val = JuMP.objective_value(setup_max.jump_model) - @test max_val >= true_val - 1e-6 - end - - @testset "Approximation quality improves with segments" begin - errors = Float64[] - true_val = 0.35^2 - for num_segments in [2, 4, 8, 16] - setup = _setup_qa_test(["dev1"], 1:1) - x_var = setup.var_container["dev1", 1] - JuMP.fix(x_var, 0.35; force = true) - - IOM._add_quadratic_approx!( - IOM.ZZBQuadConfig(num_segments), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.var_container, - 0.0, - 1.0, - ZZB_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.QuadraticExpression, - MockThermalGen, - ZZB_META, - ) - z_expr = expr_container["dev1", 1] - - JuMP.@objective(setup.jump_model, Min, z_expr) - JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) - JuMP.set_silent(setup.jump_model) - JuMP.optimize!(setup.jump_model) - - obj_val = JuMP.objective_value(setup.jump_model) - push!(errors, abs(obj_val - true_val)) - end - for i in 2:length(errors) - @test errors[i] <= errors[i - 1] + 1e-10 - end - end - - @testset "Multiple time steps" begin - setup = _setup_qa_test(["dev1"], 1:3) - IOM._add_quadratic_approx!( - IOM.ZZBQuadConfig(4), - setup.container, - MockThermalGen, - ["dev1"], - 1:3, - setup.var_container, - 0.0, - 4.0, - ZZB_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.QuadraticExpression, - MockThermalGen, - ZZB_META, - ) - - for t in 1:3 - @test expr_container["dev1", t] isa JuMP.AffExpr - end - end -end - -# ==================================================================== -# 4. HybS + ZZB bilinear integration -# ==================================================================== -@testset "HybS + ZZB Bilinear Approximation" begin - @testset "Constraint structure" begin - setup = _setup_bilinear_test(["dev1"], 1:1) - num_segments = 4 # 4 segments → 2 binary vars per quadratic - - IOM._add_bilinear_approx!( - IOM.HybSConfig(IOM.ZZBQuadConfig(num_segments), Int(log2(num_segments)), false), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, - 4.0, - 0.0, - 4.0, - ZZB_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - ZZB_META, - ) - - @test expr_container["dev1", 1] isa JuMP.AffExpr - - # Binary count: 2 * log2(num_segments) (one set for x², one for y²), zero from epigraphs - n_bin = count(JuMP.is_binary, JuMP.all_variables(setup.jump_model)) - @test n_bin == 2 * Int(log2(num_segments)) - end - - @testset "Brackets true product at interior points" begin - test_points = [(0.3, 0.7), (0.5, 0.5), (0.1, 0.9), (0.8, 0.2)] - for (x0, y0) in test_points - z_vals = Float64[] - for sense in [JuMP.MIN_SENSE, JuMP.MAX_SENSE] - setup = _setup_bilinear_test(["dev1"], 1:1) - JuMP.fix(setup.x_var_container["dev1", 1], x0; force = true) - JuMP.fix(setup.y_var_container["dev1", 1], y0; force = true) - - IOM._add_bilinear_approx!( - IOM.HybSConfig(IOM.ZZBQuadConfig(4), 2, false), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, - 1.0, - 0.0, - 1.0, - ZZB_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - ZZB_META, - ) - z_expr = expr_container["dev1", 1] - - JuMP.@objective(setup.jump_model, sense, z_expr) - JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) - JuMP.set_silent(setup.jump_model) - JuMP.optimize!(setup.jump_model) - - @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL - push!(z_vals, JuMP.objective_value(setup.jump_model)) - end - true_val = x0 * y0 - @test z_vals[1] <= true_val + 1e-6 - @test z_vals[2] >= true_val - 1e-6 - end - end - - @testset "Fixed-variable correctness" begin - setup = _setup_bilinear_test(["dev1"], 1:1) - x_var = setup.x_var_container["dev1", 1] - y_var = setup.y_var_container["dev1", 1] - JuMP.fix(x_var, 2.0; force = true) - JuMP.fix(y_var, 3.0; force = true) - - IOM._add_bilinear_approx!( - IOM.HybSConfig(IOM.ZZBQuadConfig(8), 3, false), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, - 4.0, - 0.0, - 4.0, - ZZB_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - ZZB_META, - ) - z_expr = expr_container["dev1", 1] - - JuMP.@objective(setup.jump_model, Max, z_expr) - JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) - JuMP.set_silent(setup.jump_model) - JuMP.optimize!(setup.jump_model) - - @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL - @test JuMP.objective_value(setup.jump_model) ≈ 6.0 atol = 0.5 - end - - @testset "HybS+ZZB uses fewer binaries than Bin2+ZZB" begin - for num_segments in [2, 4, 8] - binary_depth = Int(log2(num_segments)) - # HybS - setup_h = _setup_bilinear_test(["dev1"], 1:1) - IOM._add_bilinear_approx!( - IOM.HybSConfig(IOM.ZZBQuadConfig(num_segments), binary_depth, false), - setup_h.container, - MockThermalGen, - ["dev1"], - 1:1, - setup_h.x_var_container, - setup_h.y_var_container, - 0.0, - 1.0, - 0.0, - 1.0, - ZZB_META, - ) - n_bin_hybs = - count(JuMP.is_binary, JuMP.all_variables(setup_h.jump_model)) - - @test n_bin_hybs == 2 * binary_depth - end - end -end diff --git a/test/test_zzi_bilinear.jl b/test/test_zzi_bilinear.jl deleted file mode 100644 index 48c5cc4a..00000000 --- a/test/test_zzi_bilinear.jl +++ /dev/null @@ -1,466 +0,0 @@ -const ZZI_META = "ZZITest" - -@testset "ZZI Bivariate Approximation" begin - - # ==================================================================== - # Section 1: ZZI Encoding (Pure Math, No JuMP) - # ==================================================================== - @testset "ZZI encoding" begin - @testset "d=4 (r=2) known values" begin - C, C_ext = IOM.build_zzi_encoding(4) - @test size(C) == (4, 2) - # build_brgc(2) (prepend convention) = [0 0; 0 1; 1 1; 1 0] - # Cumulative transitions along each column: - # col 1: 0,0,1,1 → C[:,1] = [0, 0, 1, 1] - # col 2: 0,1,1,0 → C[:,2] = [0, 1, 1, 2] - @test C == [0 0; 0 1; 1 1; 1 2] - @test size(C_ext) == (6, 2) - # C_ext: row 1 = C_0 = C[1,:], rows 2..5 = C[1..4,:], row 6 = C[4,:] - @test C_ext[1, :] == [0, 0] # C_0 = C_1 - @test C_ext[6, :] == [1, 2] # C_5 = C_4 - end - - @testset "d=2 (r=1)" begin - C, C_ext = IOM.build_zzi_encoding(2) - @test size(C) == (2, 1) - @test C == reshape([0, 1], 2, 1) - end - - @testset "Monotonicity" begin - for d in [2, 3, 4, 5, 7, 8, 16] - C, _ = IOM.build_zzi_encoding(d) - for k in 1:size(C, 2), i in 2:size(C, 1) - @test C[i, k] >= C[i - 1, k] - end - end - end - end - - @testset "Triangulation chooser" begin - x_bkpts = [0.0, 1.0, 2.0] - y_bkpts = [0.0, 1.0, 2.0] - triang = IOM._choose_triangulation(x_bkpts, y_bkpts) - @test size(triang) == (2, 2) - @test all(t -> t == :U || t == :K, triang) - end - - # ==================================================================== - # Section 2: ZZI Standalone Solve Tests - # ==================================================================== - @testset "ZZI standalone" begin - @testset "Constraint structure" begin - # d1=d2=4 → r1=r2=2 → 4 integers (2 per axis), 6 triangle binaries - setup = _setup_bilinear_test(["dev1"], 1:1) - config = IOM.ZZIBilinearConfig(4, 4, true, false, 0) - IOM._add_bilinear_approx!( - config, - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, - 4.0, - 0.0, - 4.0, - 4, - ZZI_META, - ) - - expr = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - ZZI_META, - ) - @test expr["dev1", 1] isa JuMP.AffExpr - - n_int = count( - v -> JuMP.is_integer(v) && !JuMP.is_binary(v), - JuMP.all_variables(setup.jump_model), - ) - @test n_int == 4 # ceil(log2(4)) + ceil(log2(4)) = 2 + 2 - - n_bin = count(JuMP.is_binary, JuMP.all_variables(setup.jump_model)) - @test n_bin == 6 # triangle selection - - @test IOM.has_container_key( - setup.container, - IOM.ZZILambdaVariable, - MockThermalGen, - ZZI_META, - ) - end - - @testset "Fixed-variable correctness" begin - setup = _setup_bilinear_test(["dev1"], 1:1) - JuMP.fix(setup.x_var_container["dev1", 1], 2.0; force = true) - JuMP.fix(setup.y_var_container["dev1", 1], 3.0; force = true) - - IOM._add_bilinear_approx!( - IOM.ZZIBilinearConfig(4, 4), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, - 4.0, - 0.0, - 4.0, - 4, - ZZI_META, - ) - - z_expr = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - ZZI_META, - )[ - "dev1", - 1, - ] - JuMP.@objective(setup.jump_model, Max, z_expr) - JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) - JuMP.set_silent(setup.jump_model) - JuMP.optimize!(setup.jump_model) - - @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL - @test JuMP.objective_value(setup.jump_model) ≈ 6.0 atol = 1e-4 # exact at grid point - end - - @testset "Vertex optima" begin - setup = _setup_bilinear_test(["dev1"], 1:1) - JuMP.set_lower_bound(setup.x_var_container["dev1", 1], 0.0) - JuMP.set_upper_bound(setup.x_var_container["dev1", 1], 4.0) - JuMP.set_lower_bound(setup.y_var_container["dev1", 1], 0.0) - JuMP.set_upper_bound(setup.y_var_container["dev1", 1], 4.0) - - IOM._add_bilinear_approx!( - IOM.ZZIBilinearConfig(4, 4), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, - 4.0, - 0.0, - 4.0, - 4, - ZZI_META, - ) - - z_expr = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - ZZI_META, - )[ - "dev1", - 1, - ] - JuMP.@objective(setup.jump_model, Min, z_expr) - JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) - JuMP.set_silent(setup.jump_model) - JuMP.optimize!(setup.jump_model) - - @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL - @test JuMP.objective_value(setup.jump_model) ≈ 0.0 atol = 1e-6 - end - - @testset "Approximation quality improves with grid refinement" begin - errors = Float64[] - x0, y0 = 1.3, 2.7 - true_val = x0 * y0 - for d in [2, 4, 8] - setup = _setup_bilinear_test(["dev1"], 1:1) - JuMP.fix(setup.x_var_container["dev1", 1], x0; force = true) - JuMP.fix(setup.y_var_container["dev1", 1], y0; force = true) - - IOM._add_bilinear_approx!( - IOM.ZZIBilinearConfig(d, d), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, - 4.0, - 0.0, - 4.0, - 4, - ZZI_META, - ) - - z_expr = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - ZZI_META, - )[ - "dev1", - 1, - ] - JuMP.@objective(setup.jump_model, Min, z_expr) - JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) - JuMP.set_silent(setup.jump_model) - JuMP.optimize!(setup.jump_model) - - @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL - push!(errors, abs(JuMP.objective_value(setup.jump_model) - true_val)) - end - for i in 2:length(errors) - @test errors[i] <= errors[i - 1] + 1e-10 - end - end - - @testset "Multiple time steps" begin - setup = _setup_bilinear_test(["dev1"], 1:3) - IOM._add_bilinear_approx!( - IOM.ZZIBilinearConfig(4, 4), - setup.container, - MockThermalGen, - ["dev1"], - 1:3, - setup.x_var_container, - setup.y_var_container, - 0.0, - 4.0, - 0.0, - 4.0, - 4, - ZZI_META, - ) - - expr = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - ZZI_META, - ) - for t in 1:3 - @test expr["dev1", t] isa JuMP.AffExpr - end - end - end - - # ==================================================================== - # Section 3: Sawtooth Strengthening - # ==================================================================== - @testset "ZZI with sawtooth strengthening" begin - @testset "Constraint structure" begin - setup = _setup_bilinear_test(["dev1"], 1:1) - config = IOM.ZZIBilinearConfig(4, 4, true, true, 3) - IOM._add_bilinear_approx!( - config, - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, - 4.0, - 0.0, - 4.0, - 4, - ZZI_META, - ) - - @test IOM.has_container_key( - setup.container, - IOM.ZZISawtoothBoundConstraint, - MockThermalGen, - ZZI_META * "_st", - ) - - # 6 triangle binaries + 2*3 sawtooth binaries = 12 binaries - n_bin = count(JuMP.is_binary, JuMP.all_variables(setup.jump_model)) - @test n_bin == 12 - end - - @testset "Fixed-variable correctness at grid point" begin - # Use a grid point to avoid LP infeasibility from ZZI+sawtooth interaction - # (at non-grid points, the sawtooth LP relaxation can conflict with the - # ZZI product equality constraint) - setup = _setup_bilinear_test(["dev1"], 1:1) - JuMP.fix(setup.x_var_container["dev1", 1], 2.0; force = true) - JuMP.fix(setup.y_var_container["dev1", 1], 3.0; force = true) - - IOM._add_bilinear_approx!( - IOM.ZZIBilinearConfig(4, 4, true, true, 3), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, - 4.0, - 0.0, - 4.0, - 4, - ZZI_META, - ) - - z_expr = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - ZZI_META, - )[ - "dev1", - 1, - ] - JuMP.@objective(setup.jump_model, Max, z_expr) - JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) - JuMP.set_silent(setup.jump_model) - JuMP.optimize!(setup.jump_model) - - @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL - @test JuMP.objective_value(setup.jump_model) ≈ 6.0 atol = 1e-4 - end - end - - # ==================================================================== - # Section 4: ZZI Strengthened Wrapper with HybS - # ==================================================================== - @testset "ZZI strengthened HybS" begin - @testset "Constraint structure" begin - setup = _setup_bilinear_test(["dev1"], 1:1) - inner = IOM.HybSConfig(IOM.SawtoothQuadConfig(4), 4) - config = IOM.ZZIStrengthenedConfig(inner, 4, 4) - IOM._add_bilinear_approx!( - config, - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, - 4.0, - 0.0, - 4.0, - 3, - ZZI_META, - ) - - expr = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - ZZI_META, - ) - @test expr["dev1", 1] isa JuMP.AffExpr - - # Both HybS and ZZI containers exist - @test IOM.has_container_key( - setup.container, - IOM.HybSBoundConstraint, - MockThermalGen, - ZZI_META, - ) - meta_zzi = ZZI_META * "_zzir" - @test IOM.has_container_key( - setup.container, - IOM.ZZILambdaVariable, - MockThermalGen, - meta_zzi, - ) - end - - @testset "Fixed-variable correctness" begin - setup = _setup_bilinear_test(["dev1"], 1:1) - JuMP.fix(setup.x_var_container["dev1", 1], 2.0; force = true) - JuMP.fix(setup.y_var_container["dev1", 1], 3.0; force = true) - - inner = IOM.HybSConfig(IOM.SawtoothQuadConfig(4), 4) - config = IOM.ZZIStrengthenedConfig(inner, 4, 4) - IOM._add_bilinear_approx!( - config, - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, - 4.0, - 0.0, - 4.0, - 3, - ZZI_META, - ) - - z_expr = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - ZZI_META, - )[ - "dev1", - 1, - ] - JuMP.@objective(setup.jump_model, Max, z_expr) - JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) - JuMP.set_silent(setup.jump_model) - JuMP.optimize!(setup.jump_model) - - @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL - @test JuMP.objective_value(setup.jump_model) ≈ 6.0 atol = 0.5 - end - end - - # ==================================================================== - # Section 5: ZZI Strengthened Wrapper with Bin2 - # ==================================================================== - @testset "ZZI strengthened Bin2" begin - @testset "Fixed-variable correctness" begin - setup = _setup_bilinear_test(["dev1"], 1:1) - JuMP.fix(setup.x_var_container["dev1", 1], 2.0; force = true) - JuMP.fix(setup.y_var_container["dev1", 1], 3.0; force = true) - - inner = IOM.Bin2Config(IOM.SawtoothQuadConfig(4)) - config = IOM.ZZIStrengthenedConfig(inner, 4, 4) - IOM._add_bilinear_approx!( - config, - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, - 4.0, - 0.0, - 4.0, - 3, - ZZI_META, - ) - - z_expr = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - ZZI_META, - )[ - "dev1", - 1, - ] - JuMP.@objective(setup.jump_model, Min, z_expr) - JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) - JuMP.set_silent(setup.jump_model) - JuMP.optimize!(setup.jump_model) - - @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL - @test JuMP.objective_value(setup.jump_model) ≈ 6.0 atol = 0.5 - end - end -end From ae90871dce367b256fc84f55b2cbb5c513be0dec Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Wed, 22 Apr 2026 13:11:54 -0400 Subject: [PATCH 08/17] tests use correct bounds; pruned pointless tests --- test/test_bilinear_approximations.jl | 289 +++--------------------- test/test_hybs_approximations.jl | 171 ++------------ test/test_nmdt_approximations.jl | 230 +++---------------- test/test_quadratic_approximations.jl | 312 +------------------------- 4 files changed, 82 insertions(+), 920 deletions(-) diff --git a/test/test_bilinear_approximations.jl b/test/test_bilinear_approximations.jl index c517a418..8619eb5d 100644 --- a/test/test_bilinear_approximations.jl +++ b/test/test_bilinear_approximations.jl @@ -2,81 +2,6 @@ const BILINEAR_META = "BilinearTest" @testset "Bilinear Approximations" begin @testset "Bin2 + Solver SOS2" begin - @testset "Constraint structure" begin - setup = _setup_bilinear_test(["dev1"], 1:1) - - IOM._add_bilinear_approx!( - IOM.Bin2Config(IOM.SolverSOS2QuadConfig(4, 0)), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, - BILINEAR_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - BILINEAR_META, - ) - - # add more container tests for the new constraint expressions and p addition stuff - - @test expr_container["dev1", 1] isa JuMP.AffExpr - end - - @testset "Reformulated McCormick constraints present" begin - setup = _setup_bilinear_test(["dev1"], 1:1) - - IOM._add_bilinear_approx!( - IOM.Bin2Config(IOM.SolverSOS2QuadConfig(4, 0), true), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, - BILINEAR_META, - ) - - @test IOM.has_container_key( - setup.container, - IOM.ReformulatedMcCormickConstraint, - MockThermalGen, - BILINEAR_META, - ) - end - - @testset "No reformulated McCormick when disabled" begin - setup = _setup_bilinear_test(["dev1"], 1:1) - - IOM._add_bilinear_approx!( - IOM.Bin2Config(IOM.SolverSOS2QuadConfig(4, 0), false), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, - BILINEAR_META, - ) - - @test !IOM.has_container_key( - setup.container, - IOM.ReformulatedMcCormickConstraint, - MockThermalGen, - BILINEAR_META, - ) - end - @testset "Reformulated McCormick tightens LP relaxation" begin # At interior point x=2.5, y=1.5, true product = 3.75 # Compare max z with and without reformulated McCormick @@ -95,8 +20,8 @@ const BILINEAR_META = "BilinearTest" 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, + [(min = 0.0, max = 4.0)], + [(min = 0.0, max = 4.0)], BILINEAR_META, ) expr_container = IOM.get_expression( @@ -136,8 +61,8 @@ const BILINEAR_META = "BilinearTest" 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, + [(min = 0.0, max = 4.0)], + [(min = 0.0, max = 4.0)], BILINEAR_META, ) expr_container = IOM.get_expression( @@ -169,8 +94,8 @@ const BILINEAR_META = "BilinearTest" 1:1, setup2.x_var_container, setup2.y_var_container, - 0.0, 4.0, - 0.0, 4.0, + [(min = 0.0, max = 4.0)], + [(min = 0.0, max = 4.0)], BILINEAR_META, ) expr_container2 = IOM.get_expression( @@ -207,8 +132,8 @@ const BILINEAR_META = "BilinearTest" 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, + [(min = 0.0, max = 4.0)], + [(min = 0.0, max = 4.0)], BILINEAR_META, ) expr_container = IOM.get_expression( @@ -248,8 +173,8 @@ const BILINEAR_META = "BilinearTest" 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, + [(min = 0.0, max = 4.0)], + [(min = 0.0, max = 4.0)], BILINEAR_META, ) expr_container = IOM.get_expression( @@ -269,32 +194,6 @@ const BILINEAR_META = "BilinearTest" @test JuMP.objective_value(setup.jump_model) ≈ 0.0 atol = 1e-4 end - @testset "Multiple time steps" begin - setup = _setup_bilinear_test(["dev1"], 1:3) - IOM._add_bilinear_approx!( - IOM.Bin2Config(IOM.SolverSOS2QuadConfig(4, 0)), - setup.container, - MockThermalGen, - ["dev1"], - 1:3, - setup.x_var_container, - setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, - BILINEAR_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - BILINEAR_META, - ) - - for t in 1:3 - @test expr_container["dev1", t] isa JuMP.AffExpr - end - end - @testset "Approximation quality improves with segments" begin # Fix x=2.5, y=1.5: true product = 3.75 # Sweep segments, verify gap shrinks @@ -315,8 +214,8 @@ const BILINEAR_META = "BilinearTest" 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, + [(min = 0.0, max = 4.0)], + [(min = 0.0, max = 4.0)], BILINEAR_META, ) expr_container = IOM.get_expression( @@ -342,47 +241,6 @@ const BILINEAR_META = "BilinearTest" end @testset "Bin2 + Manual SOS2" begin - @testset "Constraint structure" begin - setup = _setup_bilinear_test(["dev1"], 1:1) - - IOM._add_bilinear_approx!( - IOM.Bin2Config(IOM.ManualSOS2QuadConfig(4, 0)), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, - BILINEAR_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - BILINEAR_META, - ) - - @test expr_container["dev1", 1] isa JuMP.AffExpr - - # Binary variables should exist for both u² and v² paths - @test IOM.has_container_key( - setup.container, - IOM.ManualSOS2BinaryVariable, - MockThermalGen, - BILINEAR_META * "_plus", - ) - - # NO solver SOS2 constraints - sos2_count = JuMP.num_constraints( - setup.jump_model, - Vector{JuMP.VariableRef}, - MOI.SOS2{Float64}, - ) - @test sos2_count == 0 - end - @testset "Fixed-variable correctness" begin setup = _setup_bilinear_test(["dev1"], 1:1) JuMP.fix(setup.x_var_container["dev1", 1], 2.0; force = true) @@ -396,8 +254,8 @@ const BILINEAR_META = "BilinearTest" 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, + [(min = 0.0, max = 4.0)], + [(min = 0.0, max = 4.0)], BILINEAR_META, ) expr_container = IOM.get_expression( @@ -432,8 +290,8 @@ const BILINEAR_META = "BilinearTest" 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, + [(min = 0.0, max = 4.0)], + [(min = 0.0, max = 4.0)], BILINEAR_META, ) expr_container = IOM.get_expression( @@ -471,8 +329,8 @@ const BILINEAR_META = "BilinearTest" 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, + [(min = 0.0, max = 4.0)], + [(min = 0.0, max = 4.0)], BILINEAR_META, ) expr_container = IOM.get_expression( @@ -492,32 +350,6 @@ const BILINEAR_META = "BilinearTest" @test JuMP.objective_value(setup.jump_model) ≈ 0.0 atol = 1e-4 end - @testset "Multiple time steps" begin - setup = _setup_bilinear_test(["dev1"], 1:3) - IOM._add_bilinear_approx!( - IOM.Bin2Config(IOM.ManualSOS2QuadConfig(4, 0)), - setup.container, - MockThermalGen, - ["dev1"], - 1:3, - setup.x_var_container, - setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, - BILINEAR_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - BILINEAR_META, - ) - - for t in 1:3 - @test expr_container["dev1", t] isa JuMP.AffExpr - end - end - @testset "Approximation quality improves with segments" begin true_product = 2.5 * 1.5 errors = Float64[] @@ -534,8 +366,8 @@ const BILINEAR_META = "BilinearTest" 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, + [(min = 0.0, max = 4.0)], + [(min = 0.0, max = 4.0)], BILINEAR_META, ) expr_container = IOM.get_expression( @@ -561,45 +393,6 @@ const BILINEAR_META = "BilinearTest" end @testset "Bin2 + Sawtooth" begin - @testset "Constraint structure" begin - setup = _setup_bilinear_test(["dev1"], 1:1) - - IOM._add_bilinear_approx!( - IOM.Bin2Config(IOM.SawtoothQuadConfig(2)), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, - BILINEAR_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - BILINEAR_META, - ) - - @test expr_container["dev1", 1] isa JuMP.AffExpr - - # Sawtooth aux/binary variables for both u² and v² paths - @test IOM.has_container_key( - setup.container, - IOM.SawtoothAuxVariable, - MockThermalGen, - BILINEAR_META * "_plus", - ) - @test IOM.has_container_key( - setup.container, - IOM.SawtoothBinaryVariable, - MockThermalGen, - BILINEAR_META * "_plus", - ) - end - @testset "Fixed-variable correctness" begin setup = _setup_bilinear_test(["dev1"], 1:1) JuMP.fix(setup.x_var_container["dev1", 1], 2.0; force = true) @@ -613,8 +406,8 @@ const BILINEAR_META = "BilinearTest" 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, + [(min = 0.0, max = 4.0)], + [(min = 0.0, max = 4.0)], BILINEAR_META, ) expr_container = IOM.get_expression( @@ -649,8 +442,8 @@ const BILINEAR_META = "BilinearTest" 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, + [(min = 0.0, max = 4.0)], + [(min = 0.0, max = 4.0)], BILINEAR_META, ) expr_container = IOM.get_expression( @@ -688,8 +481,8 @@ const BILINEAR_META = "BilinearTest" 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, + [(min = 0.0, max = 4.0)], + [(min = 0.0, max = 4.0)], BILINEAR_META, ) expr_container = IOM.get_expression( @@ -709,32 +502,6 @@ const BILINEAR_META = "BilinearTest" @test JuMP.objective_value(setup.jump_model) ≈ 0.0 atol = 1e-3 end - @testset "Multiple time steps" begin - setup = _setup_bilinear_test(["dev1"], 1:3) - IOM._add_bilinear_approx!( - IOM.Bin2Config(IOM.SawtoothQuadConfig(2)), - setup.container, - MockThermalGen, - ["dev1"], - 1:3, - setup.x_var_container, - setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, - BILINEAR_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - BILINEAR_META, - ) - - for t in 1:3 - @test expr_container["dev1", t] isa JuMP.AffExpr - end - end - @testset "Approximation quality improves with depth" begin true_product = 2.5 * 1.5 errors = Float64[] @@ -751,8 +518,8 @@ const BILINEAR_META = "BilinearTest" 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 4.0, - 0.0, 4.0, + [(min = 0.0, max = 4.0)], + [(min = 0.0, max = 4.0)], BILINEAR_META, ) expr_container = IOM.get_expression( diff --git a/test/test_hybs_approximations.jl b/test/test_hybs_approximations.jl index 2f41340f..da2ab025 100644 --- a/test/test_hybs_approximations.jl +++ b/test/test_hybs_approximations.jl @@ -2,35 +2,6 @@ const HYBS_META = "HybSTest" const HYBS_BILINEAR_META = "BilinearTest" @testset "Epigraph Quadratic Approximation" begin - @testset "Constraint structure" begin - setup = _setup_qa_test(["dev1"], 1:1) - depth = 3 - - IOM._add_quadratic_approx!( - IOM.EpigraphQuadConfig(depth), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.var_container, - 0.0, - 4.0, - HYBS_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.EpigraphExpression, - MockThermalGen, - HYBS_META, - ) - - @test expr_container["dev1", 1] isa JuMP.AffExpr - - # No binary variables (pure LP) - n_bin = count(JuMP.is_binary, JuMP.all_variables(setup.jump_model)) - @test n_bin == 0 - end - @testset "Lower-bounds x^2 on [0,1]" begin setup = _setup_qa_test(["dev1"], 1:1) x_var = setup.var_container["dev1", 1] @@ -43,8 +14,7 @@ const HYBS_BILINEAR_META = "BilinearTest" ["dev1"], 1:1, setup.var_container, - 0.0, - 1.0, + [(min = 0.0, max = 1.0)], HYBS_META, ) expr_container = IOM.get_expression( @@ -77,8 +47,7 @@ const HYBS_BILINEAR_META = "BilinearTest" ["dev1"], 1:1, setup.var_container, - 0.0, - 2.0, + [(min = 0.0, max = 2.0)], HYBS_META, ) expr_container = IOM.get_expression( @@ -99,31 +68,6 @@ const HYBS_BILINEAR_META = "BilinearTest" @test JuMP.objective_value(setup.jump_model) >= 1.69 - 0.1 end - @testset "Multiple time steps" begin - setup = _setup_qa_test(["dev1"], 1:3) - IOM._add_quadratic_approx!( - IOM.EpigraphQuadConfig(2), - setup.container, - MockThermalGen, - ["dev1"], - 1:3, - setup.var_container, - 0.0, - 4.0, - HYBS_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.EpigraphExpression, - MockThermalGen, - HYBS_META, - ) - - for t in 1:3 - @test expr_container["dev1", t] isa JuMP.AffExpr - end - end - @testset "Approximation quality improves with depth" begin errors = Float64[] true_val = 0.35^2 @@ -139,8 +83,7 @@ const HYBS_BILINEAR_META = "BilinearTest" ["dev1"], 1:1, setup.var_container, - 0.0, - 1.0, + [(min = 0.0, max = 1.0)], HYBS_META, ) expr_container = IOM.get_expression( @@ -166,38 +109,6 @@ const HYBS_BILINEAR_META = "BilinearTest" end @testset "HybS Bilinear Approximation" begin - @testset "Constraint structure" begin - setup = _setup_bilinear_test(["dev1"], 1:1) - depth = 2 - - IOM._add_bilinear_approx!( - IOM.HybSConfig(IOM.SawtoothQuadConfig(depth), depth), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.x_var_container, - setup.y_var_container, - 0.0, - 4.0, - 0.0, - 4.0, - HYBS_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - HYBS_META, - ) - - @test expr_container["dev1", 1] isa JuMP.AffExpr - - # Binary count: 2L (L for x², L for y²), zero from epigraphs - n_bin = count(JuMP.is_binary, JuMP.all_variables(setup.jump_model)) - @test n_bin == 2 * depth - end - @testset "Brackets true product at interior points" begin test_points = [(0.3, 0.7), (0.5, 0.5), (0.1, 0.9), (0.8, 0.2)] for (x0, y0) in test_points @@ -215,10 +126,8 @@ end 1:1, setup.x_var_container, setup.y_var_container, - 0.0, - 1.0, - 0.0, - 1.0, + [(min = 0.0, max = 1.0)], + [(min = 0.0, max = 1.0)], HYBS_META, ) expr_container = IOM.get_expression( @@ -258,10 +167,8 @@ end 1:1, setup.x_var_container, setup.y_var_container, - 0.0, - 4.0, - 0.0, - 4.0, + [(min = 0.0, max = 4.0)], + [(min = 0.0, max = 4.0)], HYBS_META, ) expr_container = IOM.get_expression( @@ -298,10 +205,8 @@ end 1:1, setup.x_var_container, setup.y_var_container, - 0.0, - 4.0, - 0.0, - 4.0, + [(min = 0.0, max = 4.0)], + [(min = 0.0, max = 4.0)], HYBS_META, ) expr_container = IOM.get_expression( @@ -338,10 +243,8 @@ end 1:1, setup.x_var_container, setup.y_var_container, - 0.0, - 1.0, - 0.0, - 1.0, + [(min = 0.0, max = 1.0)], + [(min = 0.0, max = 1.0)], HYBS_META, ) expr_container = IOM.get_expression( @@ -381,10 +284,8 @@ end 1:1, setup.x_var_container, setup.y_var_container, - 2.0, - 5.0, - 1.0, - 3.0, + [(min = 2.0, max = 5.0)], + [(min = 1.0, max = 3.0)], HYBS_META, ) expr_container = IOM.get_expression( @@ -407,34 +308,6 @@ end @test z_vals[2] >= true_product - 1e-6 end - @testset "Multiple time steps" begin - setup = _setup_bilinear_test(["dev1"], 1:3) - IOM._add_bilinear_approx!( - IOM.HybSConfig(IOM.SawtoothQuadConfig(2), 2), - setup.container, - MockThermalGen, - ["dev1"], - 1:3, - setup.x_var_container, - setup.y_var_container, - 0.0, - 4.0, - 0.0, - 4.0, - HYBS_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.BilinearProductExpression, - MockThermalGen, - HYBS_META, - ) - - for t in 1:3 - @test expr_container["dev1", t] isa JuMP.AffExpr - end - end - @testset "Vertex optimum" begin setup = _setup_bilinear_test(["dev1"], 1:1) x_var = setup.x_var_container["dev1", 1] @@ -452,10 +325,8 @@ end 1:1, setup.x_var_container, setup.y_var_container, - 0.0, - 4.0, - 0.0, - 4.0, + [(min = 0.0, max = 4.0)], + [(min = 0.0, max = 4.0)], HYBS_META, ) expr_container = IOM.get_expression( @@ -488,10 +359,8 @@ end 1:1, setup_h.x_var_container, setup_h.y_var_container, - 0.0, - 1.0, - 0.0, - 1.0, + [(min = 0.0, max = 1.0)], + [(min = 0.0, max = 1.0)], HYBS_META, ) n_bin_hybs = @@ -507,10 +376,8 @@ end 1:1, setup_b.x_var_container, setup_b.y_var_container, - 0.0, - 1.0, - 0.0, - 1.0, + [(min = 0.0, max = 1.0)], + [(min = 0.0, max = 1.0)], HYBS_BILINEAR_META, ) n_bin_bin2 = diff --git a/test/test_nmdt_approximations.jl b/test/test_nmdt_approximations.jl index cf71c047..a428faf2 100644 --- a/test/test_nmdt_approximations.jl +++ b/test/test_nmdt_approximations.jl @@ -15,7 +15,7 @@ const NMDT_BILINEAR_META = "NMDTBilinearTest" IOM._add_quadratic_approx!( IOM.DNMDTQuadConfig(4, 0), setup.container, MockThermalGen, names, ts, - setup.var_container, 0.0, 1.0, DNMDT_META, + setup.var_container, [(min = 0.0, max = 1.0)], DNMDT_META, ) JuMP.@objective(setup.jump_model, Min, 0) @@ -49,7 +49,7 @@ const NMDT_BILINEAR_META = "NMDTBilinearTest" IOM._add_quadratic_approx!( IOM.DNMDTQuadConfig(3, 0), setup.container, MockThermalGen, ["gen1"], 1:1, - setup.var_container, 0.0, 1.0, DNMDT_META, + setup.var_container, [(min = 0.0, max = 1.0)], DNMDT_META, ) expr = IOM.get_expression( setup.container, IOM.QuadraticExpression, @@ -80,7 +80,7 @@ const NMDT_BILINEAR_META = "NMDTBilinearTest" IOM._add_quadratic_approx!( IOM.DNMDTQuadConfig(2 * L, 0), setup.container, MockThermalGen, ["gen1"], 1:1, - setup.var_container, 0.0, 1.0, DNMDT_META, + setup.var_container, [(min = 0.0, max = 1.0)], DNMDT_META, ) expr = IOM.get_expression( setup.container, IOM.QuadraticExpression, @@ -99,48 +99,6 @@ const NMDT_BILINEAR_META = "NMDTBilinearTest" end end - @testset "Constraint structure" begin - setup = _setup_qa_test(["gen1"], 1:1) - depth = 3 - - IOM._add_quadratic_approx!( - IOM.DNMDTQuadConfig(depth, 0), - setup.container, MockThermalGen, ["gen1"], 1:1, - setup.var_container, 0.0, 1.0, DNMDT_META, - ) - - # L binary variables for univariate - n_bin = count(JuMP.is_binary, JuMP.all_variables(setup.jump_model)) - @test n_bin == depth - - # Container keys exist - @test IOM.has_container_key( - setup.container, IOM.NMDTBinaryVariable, MockThermalGen, DNMDT_META, - ) - @test IOM.has_container_key( - setup.container, IOM.NMDTResidualVariable, MockThermalGen, DNMDT_META, - ) - end - - @testset "Multiple time steps and names" begin - names = ["gen1", "gen2"] - ts = 1:3 - setup = _setup_qa_test(names, ts) - - IOM._add_quadratic_approx!( - IOM.DNMDTQuadConfig(2, 0), - setup.container, MockThermalGen, names, ts, - setup.var_container, 0.0, 1.0, DNMDT_META, - ) - expr = IOM.get_expression( - setup.container, IOM.QuadraticExpression, - MockThermalGen, DNMDT_META, - ) - - for name in names, t in ts - @test expr[name, t] isa JuMP.AffExpr - end - end end @testset "T-D-NMDT Tightening" begin @@ -155,7 +113,7 @@ end IOM._add_quadratic_approx!( (tighten ? IOM.DNMDTQuadConfig(2) : IOM.DNMDTQuadConfig(2, 0)), setup.container, MockThermalGen, ["gen1"], 1:1, - setup.var_container, 0.0, 1.0, DNMDT_META, + setup.var_container, [(min = 0.0, max = 1.0)], DNMDT_META, ) expr = IOM.get_expression( setup.container, IOM.QuadraticExpression, @@ -192,7 +150,7 @@ end IOM._add_quadratic_approx!( IOM.DNMDTQuadConfig(L), setup.container, MockThermalGen, ["gen1"], 1:1, - setup.var_container, 0.0, 1.0, DNMDT_META, + setup.var_container, [(min = 0.0, max = 1.0)], DNMDT_META, ) expr = IOM.get_expression( setup.container, IOM.QuadraticExpression, @@ -227,7 +185,7 @@ end IOM.DNMDTBilinearConfig(2), setup.container, MockThermalGen, ["dev1"], 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 1.0, 0.0, 1.0, DNMDT_META, + [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], DNMDT_META, ) expr = IOM.get_expression( setup.container, IOM.BilinearProductExpression, @@ -261,7 +219,7 @@ end IOM.DNMDTBilinearConfig(2 * L), setup.container, MockThermalGen, ["dev1"], 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 1.0, 0.0, 1.0, DNMDT_META, + [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], DNMDT_META, ) expr = IOM.get_expression( setup.container, IOM.BilinearProductExpression, @@ -296,7 +254,7 @@ end IOM.DNMDTBilinearConfig(8), setup.container, MockThermalGen, ["dev1"], 1:1, setup.x_var_container, setup.y_var_container, - x_min, x_max, y_min, y_max, DNMDT_META, + [(min = x_min, max = x_max)], [(min = y_min, max = y_max)], DNMDT_META, ) expr = IOM.get_expression( setup.container, IOM.BilinearProductExpression, @@ -316,25 +274,6 @@ end end end - @testset "Constraint structure" begin - setup = _setup_bilinear_test(["dev1"], 1:1) - depth = 2 - - IOM._add_bilinear_approx!( - IOM.DNMDTBilinearConfig(depth), - setup.container, MockThermalGen, ["dev1"], 1:1, - setup.x_var_container, setup.y_var_container, - 0.0, 1.0, 0.0, 1.0, DNMDT_META, - ) - # n_bin = count(JuMP.is_binary, JuMP.all_variables(setup.jump_model)) - # @test n_bin == 2 * depth - - # Container keys exist - @test IOM.has_container_key( - setup.container, IOM.BilinearProductExpression, MockThermalGen, - DNMDT_META, - ) - end # @testset "McCormick toggle" begin # setup = _setup_bilinear_test(["dev1"], 1:1) @@ -360,7 +299,7 @@ end IOM.DNMDTBilinearConfig(3), setup.container, MockThermalGen, ["dev1"], 1:1, setup.y_var_container, setup.x_var_container, - 0.0, 4.0, 0.0, 4.0, DNMDT_META, + [(min = 0.0, max = 4.0)], [(min = 0.0, max = 4.0)], DNMDT_META, ) expr = IOM.get_expression( setup.container, IOM.BilinearProductExpression, @@ -376,40 +315,6 @@ end @test JuMP.objective_value(setup.jump_model) ≈ 6.0 atol = 0.5 end - @testset "Multiple time steps" begin - setup = _setup_bilinear_test(["dev1"], 1:3) - - IOM._add_bilinear_approx!( - IOM.DNMDTBilinearConfig(2), - setup.container, MockThermalGen, ["dev1"], 1:3, - setup.x_var_container, setup.y_var_container, - 0.0, 1.0, 0.0, 1.0, DNMDT_META, - ) - expr = IOM.get_expression( - setup.container, IOM.BilinearProductExpression, - MockThermalGen, DNMDT_META, - ) - - for t in 1:3 - @test expr["dev1", t] isa JuMP.AffExpr - end - end - - @testset "Constraint structure: 2L binary variables" begin - setup = _setup_bilinear_test(["dev1"], 1:1) - depth = 3 - - IOM._add_bilinear_approx!( - IOM.DNMDTBilinearConfig(depth), - setup.container, MockThermalGen, ["dev1"], 1:1, - setup.x_var_container, setup.y_var_container, - 0.0, 1.0, 0.0, 1.0, DNMDT_META, - ) - - # D-NMDT discretizes both x and y → 2L binary variables - n_bin = count(JuMP.is_binary, JuMP.all_variables(setup.jump_model)) - @test n_bin == 2 * depth - end @testset "D-NMDT vs HybS comparison" begin true_product = 0.4 * 0.7 @@ -423,7 +328,7 @@ end IOM.DNMDTBilinearConfig(depth), setup_d.container, MockThermalGen, ["dev1"], 1:1, setup_d.x_var_container, setup_d.y_var_container, - 0.0, 1.0, 0.0, 1.0, DNMDT_META, + [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], DNMDT_META, ) expr_d = IOM.get_expression( setup_d.container, IOM.BilinearProductExpression, @@ -445,7 +350,7 @@ end IOM.HybSConfig(IOM.SawtoothQuadConfig(depth), depth), setup_h.container, MockThermalGen, ["dev1"], 1:1, setup_h.x_var_container, setup_h.y_var_container, - 0.0, 1.0, 0.0, 1.0, DNMDT_HYBS_META, + [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], DNMDT_HYBS_META, ) expr_h = IOM.get_expression( setup_h.container, IOM.BilinearProductExpression, @@ -481,7 +386,7 @@ end IOM._add_quadratic_approx!( IOM.NMDTQuadConfig(4, 0), setup.container, MockThermalGen, names, ts, - setup.var_container, 0.0, 1.0, NMDT_META, + setup.var_container, [(min = 0.0, max = 1.0)], NMDT_META, ) JuMP.@objective(setup.jump_model, Min, 0) @@ -515,7 +420,7 @@ end IOM._add_quadratic_approx!( IOM.NMDTQuadConfig(3, 0), setup.container, MockThermalGen, ["gen1"], 1:1, - setup.var_container, 0.0, 1.0, NMDT_META, + setup.var_container, [(min = 0.0, max = 1.0)], NMDT_META, ) expr = IOM.get_expression( setup.container, IOM.QuadraticExpression, @@ -548,7 +453,7 @@ end IOM._add_quadratic_approx!( IOM.NMDTQuadConfig(L, 0), setup.container, MockThermalGen, ["gen1"], 1:1, - setup.var_container, 0.0, 1.0, NMDT_META, + setup.var_container, [(min = 0.0, max = 1.0)], NMDT_META, ) expr = IOM.get_expression( setup.container, IOM.QuadraticExpression, @@ -567,51 +472,6 @@ end end end - @testset "Constraint structure: L binary variables" begin - setup = _setup_qa_test(["gen1"], 1:1) - depth = 3 - - IOM._add_quadratic_approx!( - IOM.NMDTQuadConfig(depth, 0), - setup.container, MockThermalGen, ["gen1"], 1:1, - setup.var_container, 0.0, 1.0, NMDT_META, - ) - - # NMDT uses exactly L binary variables (one per discretization level) - n_bin = count(JuMP.is_binary, JuMP.all_variables(setup.jump_model)) - @test n_bin == depth - - @test IOM.has_container_key( - setup.container, IOM.NMDTBinaryVariable, MockThermalGen, NMDT_META, - ) - @test IOM.has_container_key( - setup.container, IOM.NMDTResidualVariable, MockThermalGen, NMDT_META, - ) - @test IOM.has_container_key( - setup.container, IOM.QuadraticExpression, MockThermalGen, NMDT_META, - ) - end - - @testset "Multiple time steps and names" begin - names = ["gen1", "gen2"] - ts = 1:3 - setup = _setup_qa_test(names, ts) - - IOM._add_quadratic_approx!( - IOM.NMDTQuadConfig(2, 0), - setup.container, MockThermalGen, names, ts, - setup.var_container, 0.0, 1.0, NMDT_META, - ) - expr = IOM.get_expression( - setup.container, IOM.QuadraticExpression, - MockThermalGen, NMDT_META, - ) - - for name in names, t in ts - @test expr[name, t] isa JuMP.AffExpr - end - end - @testset "Tightening improves lower bound" begin for x0 in [0.15, 0.35, 0.65, 0.85] lb_nmdt = NaN @@ -623,7 +483,7 @@ end IOM._add_quadratic_approx!( (tighten ? IOM.NMDTQuadConfig(2) : IOM.NMDTQuadConfig(2, 0)), setup.container, MockThermalGen, ["gen1"], 1:1, - setup.var_container, 0.0, 1.0, NMDT_META, + setup.var_container, [(min = 0.0, max = 1.0)], NMDT_META, ) expr = IOM.get_expression( setup.container, IOM.QuadraticExpression, @@ -669,7 +529,7 @@ end IOM._add_quadratic_approx!( config_fn(L), setup.container, MockThermalGen, ["gen1"], 1:1, - setup.var_container, 0.0, 1.0, NMDT_META, + setup.var_container, [(min = 0.0, max = 1.0)], NMDT_META, ) expr = IOM.get_expression( setup.container, IOM.QuadraticExpression, @@ -704,7 +564,7 @@ end IOM._add_quadratic_approx!( IOM.NMDTQuadConfig(depth, 0), setup_n.container, MockThermalGen, ["gen1"], 1:1, - setup_n.var_container, 0.0, 1.0, NMDT_META, + setup_n.var_container, [(min = 0.0, max = 1.0)], NMDT_META, ) n_bin_nmdt = count(JuMP.is_binary, JuMP.all_variables(setup_n.jump_model)) @@ -712,7 +572,7 @@ end IOM._add_quadratic_approx!( IOM.DNMDTQuadConfig(depth, 0), setup_d.container, MockThermalGen, ["gen1"], 1:1, - setup_d.var_container, 0.0, 1.0, DNMDT_META, + setup_d.var_container, [(min = 0.0, max = 1.0)], DNMDT_META, ) n_bin_dnmdt = count(JuMP.is_binary, JuMP.all_variables(setup_d.jump_model)) @@ -739,7 +599,7 @@ end IOM.NMDTBilinearConfig(3), setup.container, MockThermalGen, ["dev1"], 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 1.0, 0.0, 1.0, NMDT_BILINEAR_META, + [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], NMDT_BILINEAR_META, ) expr = IOM.get_expression( setup.container, IOM.BilinearProductExpression, @@ -774,7 +634,7 @@ end IOM.NMDTBilinearConfig(L), setup.container, MockThermalGen, ["dev1"], 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 1.0, 0.0, 1.0, NMDT_BILINEAR_META, + [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], NMDT_BILINEAR_META, ) expr = IOM.get_expression( setup.container, IOM.BilinearProductExpression, @@ -795,46 +655,6 @@ end end end - @testset "Constraint structure: L binary variables (only x discretized)" begin - setup = _setup_bilinear_test(["dev1"], 1:1) - depth = 3 - - IOM._add_bilinear_approx!( - IOM.NMDTBilinearConfig(depth), - setup.container, MockThermalGen, ["dev1"], 1:1, - setup.x_var_container, setup.y_var_container, - 0.0, 1.0, 0.0, 1.0, NMDT_BILINEAR_META, - ) - - # NMDT only discretizes x → L binary variables (half of D-NMDT's 2L) - n_bin = count(JuMP.is_binary, JuMP.all_variables(setup.jump_model)) - @test n_bin == depth - - @test IOM.has_container_key( - setup.container, IOM.BilinearProductExpression, MockThermalGen, - NMDT_BILINEAR_META, - ) - end - - @testset "Multiple time steps" begin - setup = _setup_bilinear_test(["dev1"], 1:3) - - IOM._add_bilinear_approx!( - IOM.NMDTBilinearConfig(2), - setup.container, MockThermalGen, ["dev1"], 1:3, - setup.x_var_container, setup.y_var_container, - 0.0, 1.0, 0.0, 1.0, NMDT_BILINEAR_META, - ) - expr = IOM.get_expression( - setup.container, IOM.BilinearProductExpression, - MockThermalGen, NMDT_BILINEAR_META, - ) - - for t in 1:3 - @test expr["dev1", t] isa JuMP.AffExpr - end - end - @testset "Fixed-variable correctness" begin setup = _setup_bilinear_test(["dev1"], 1:1) JuMP.fix(setup.x_var_container["dev1", 1], 0.75; force = true) @@ -844,7 +664,7 @@ end IOM.NMDTBilinearConfig(4), setup.container, MockThermalGen, ["dev1"], 1:1, setup.x_var_container, setup.y_var_container, - 0.0, 1.0, 0.0, 1.0, NMDT_BILINEAR_META, + [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], NMDT_BILINEAR_META, ) expr = IOM.get_expression( setup.container, IOM.BilinearProductExpression, @@ -870,7 +690,7 @@ end IOM.NMDTBilinearConfig(L), setup_n.container, MockThermalGen, ["dev1"], 1:1, setup_n.x_var_container, setup_n.y_var_container, - 0.0, 1.0, 0.0, 1.0, NMDT_BILINEAR_META, + [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], NMDT_BILINEAR_META, ) n_bin_nmdt = count(JuMP.is_binary, JuMP.all_variables(setup_n.jump_model)) @@ -879,7 +699,7 @@ end IOM.DNMDTBilinearConfig(L), setup_d.container, MockThermalGen, ["dev1"], 1:1, setup_d.x_var_container, setup_d.y_var_container, - 0.0, 1.0, 0.0, 1.0, DNMDT_META, + [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], DNMDT_META, ) n_bin_dnmdt = count(JuMP.is_binary, JuMP.all_variables(setup_d.jump_model)) @@ -900,7 +720,7 @@ end IOM.NMDTBilinearConfig(L), setup_n.container, MockThermalGen, ["dev1"], 1:1, setup_n.x_var_container, setup_n.y_var_container, - 0.0, 1.0, 0.0, 1.0, NMDT_BILINEAR_META, + [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], NMDT_BILINEAR_META, ) expr_n = IOM.get_expression( setup_n.container, IOM.BilinearProductExpression, @@ -921,7 +741,7 @@ end IOM.DNMDTBilinearConfig(L), setup_d.container, MockThermalGen, ["dev1"], 1:1, setup_d.x_var_container, setup_d.y_var_container, - 0.0, 1.0, 0.0, 1.0, DNMDT_META, + [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], DNMDT_META, ) expr_d = IOM.get_expression( setup_d.container, IOM.BilinearProductExpression, diff --git a/test/test_quadratic_approximations.jl b/test/test_quadratic_approximations.jl index 36c39f75..4bc8f4d0 100644 --- a/test/test_quadratic_approximations.jl +++ b/test/test_quadratic_approximations.jl @@ -3,72 +3,6 @@ const TEST_META = "TestVar" @testset "Quadratic Approximations" begin @testset "Solver SOS2" begin - @testset "Constraint structure" begin - setup = _setup_qa_test(["dev1"], 1:1) - num_segments = 4 - n_points = num_segments + 1 - - IOM._add_quadratic_approx!( - IOM.SolverSOS2QuadConfig(num_segments, 0), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.var_container, - 0.0, - 4.0, - TEST_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.QuadraticExpression, - MockThermalGen, - TEST_META, - ) - x_sq = expr_container["dev1", 1] - - # Returned expression should be AffExpr - @test x_sq isa JuMP.AffExpr - - # Lambda variables should exist - lambda_container = IOM.get_variable( - setup.container, - IOM.QuadraticVariable, - MockThermalGen, - TEST_META, - ) - for i in 1:n_points - @test haskey(lambda_container, ("dev1", i, 1)) - var = lambda_container[("dev1", i, 1)] - @test JuMP.lower_bound(var) == 0.0 - @test JuMP.upper_bound(var) == 1.0 - end - - # Linking constraint should exist - @test IOM.has_container_key( - setup.container, - IOM.SOS2LinkingConstraint, - MockThermalGen, - TEST_META, - ) - - # Normalization constraint should exist - @test IOM.has_container_key( - setup.container, - IOM.SOS2NormConstraint, - MockThermalGen, - TEST_META, - ) - - # SOS2 constraint should exist (solver-native) - sos2_count = JuMP.num_constraints( - setup.jump_model, - Vector{JuMP.VariableRef}, - MOI.SOS2{Float64}, - ) - @test sos2_count == 1 - end - @testset "Solve min x^2 - 4x" begin # Analytic minimum of x^2 - 4x at x=2, value = -4 # With breakpoints at 0,1,2,3,4 the approximation is exact at breakpoints @@ -84,8 +18,7 @@ const TEST_META = "TestVar" ["dev1"], 1:1, setup.var_container, - 0.0, - 4.0, + [(min = 0.0, max = 4.0)], TEST_META, ) expr_container = IOM.get_expression( @@ -121,8 +54,7 @@ const TEST_META = "TestVar" ["dev1"], 1:1, setup.var_container, - 0.0, - 4.0, + [(min = 0.0, max = 4.0)], TEST_META, ) expr_container = IOM.get_expression( @@ -144,43 +76,6 @@ const TEST_META = "TestVar" @test JuMP.value(y) ≈ 1.0 atol = 1e-6 end - @testset "Multiple time steps" begin - setup = _setup_qa_test(["dev1"], 1:3) - IOM._add_quadratic_approx!( - IOM.SolverSOS2QuadConfig(4, 0), - setup.container, - MockThermalGen, - ["dev1"], - 1:3, - setup.var_container, - 0.0, - 4.0, - TEST_META, - ) - - # Verify lambda variables exist for each time step - lambda_container = IOM.get_variable( - setup.container, - IOM.QuadraticVariable, - MockThermalGen, - TEST_META, - ) - for t in 1:3, i in 1:5 - @test haskey(lambda_container, ("dev1", i, t)) - end - - # Expression container should have entries for all (name, t) pairs - expr_container = IOM.get_expression( - setup.container, - IOM.QuadraticExpression, - MockThermalGen, - TEST_META, - ) - for t in 1:3 - @test expr_container["dev1", t] isa JuMP.AffExpr - end - end - @testset "Approximation quality improves with more segments" begin # min (√2)x² - (√3)x on [0, 6], analytic minimum at x=√3/8 analytic_min = sqrt(2) * (sqrt(3 / 8)^2) - sqrt(3) * sqrt(3 / 8) @@ -198,8 +93,7 @@ const TEST_META = "TestVar" ["dev1"], 1:1, setup.var_container, - 0.0, - 6.0, + [(min = 0.0, max = 6.0)], TEST_META, ) expr_container = IOM.get_expression( @@ -226,73 +120,6 @@ const TEST_META = "TestVar" end @testset "Manual SOS2" begin - @testset "Constraint structure" begin - setup = _setup_qa_test(["dev1"], 1:1) - num_segments = 4 - n_points = num_segments + 1 - - IOM._add_quadratic_approx!( - IOM.ManualSOS2QuadConfig(num_segments, 0), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.var_container, - 0.0, - 4.0, - TEST_META, - ) - expr_container = IOM.get_expression( - setup.container, - IOM.QuadraticExpression, - MockThermalGen, - TEST_META, - ) - x_sq = expr_container["dev1", 1] - - # Returned expression should be AffExpr - @test x_sq isa JuMP.AffExpr - - # Lambda variables should exist - lambda_container = IOM.get_variable( - setup.container, - IOM.QuadraticVariable, - MockThermalGen, - TEST_META, - ) - for i in 1:n_points - @test haskey(lambda_container, ("dev1", i, 1)) - end - - # Binary z variables should exist (n_points - 1) - z_container = IOM.get_variable( - setup.container, - IOM.ManualSOS2BinaryVariable, - MockThermalGen, - TEST_META, - ) - for j in 1:(n_points - 1) - @test haskey(z_container, ("dev1", j, 1)) - @test JuMP.is_binary(z_container[("dev1", j, 1)]) - end - - # Segment selection constraint should exist - @test IOM.has_container_key( - setup.container, - IOM.ManualSOS2SegmentSelectionConstraint, - MockThermalGen, - TEST_META, - ) - - # NO solver SOS2 constraints - sos2_count = JuMP.num_constraints( - setup.jump_model, - Vector{JuMP.VariableRef}, - MOI.SOS2{Float64}, - ) - @test sos2_count == 0 - end - @testset "Solve min x^2 - 4x" begin setup = _setup_qa_test(["dev1"], 1:1) x_var = setup.var_container["dev1", 1] @@ -306,8 +133,7 @@ const TEST_META = "TestVar" ["dev1"], 1:1, setup.var_container, - 0.0, - 4.0, + [(min = 0.0, max = 4.0)], TEST_META, ) expr_container = IOM.get_expression( @@ -342,8 +168,7 @@ const TEST_META = "TestVar" ["dev1"], 1:1, setup.var_container, - 0.0, - 4.0, + [(min = 0.0, max = 4.0)], TEST_META, ) expr_container = IOM.get_expression( @@ -366,72 +191,6 @@ const TEST_META = "TestVar" end @testset "Sawtooth" begin - @testset "Constraint structure" begin - setup = _setup_qa_test(["dev1"], 1:1) - depth = 2 - - IOM._add_quadratic_approx!( - IOM.SawtoothQuadConfig(depth), - setup.container, - MockThermalGen, - ["dev1"], - 1:1, - setup.var_container, - 0.0, - 4.0, - TEST_META, - ) - - # Expression container should contain AffExpr for each (name, t) - expr_container = IOM.get_expression( - setup.container, - IOM.QuadraticExpression, - MockThermalGen, - TEST_META, - ) - @test expr_container["dev1", 1] isa JuMP.AffExpr - - # Auxiliary variables g_0, g_1, g_2 should exist - g_container = IOM.get_variable( - setup.container, - IOM.SawtoothAuxVariable, - MockThermalGen, - TEST_META, - ) - for j in 0:depth - var = g_container["dev1", j, 1] - @test JuMP.lower_bound(var) == 0.0 - @test JuMP.upper_bound(var) == 1.0 - end - - # Binary variables α_1, α_2 should exist - alpha_container = IOM.get_variable( - setup.container, - IOM.SawtoothBinaryVariable, - MockThermalGen, - TEST_META, - ) - for j in 1:depth - @test JuMP.is_binary(alpha_container["dev1", j, 1]) - end - - # Linking constraint should exist - @test IOM.has_container_key( - setup.container, - IOM.SawtoothLinkingConstraint, - MockThermalGen, - TEST_META, - ) - - # NO solver SOS2 constraints - sos2_count = JuMP.num_constraints( - setup.jump_model, - Vector{JuMP.VariableRef}, - MOI.SOS2{Float64}, - ) - @test sos2_count == 0 - end - @testset "Solve min x^2 - 4x" begin # depth=2 → breakpoints at 0,1,2,3,4 → exact at x=2 setup = _setup_qa_test(["dev1"], 1:1) @@ -446,8 +205,7 @@ const TEST_META = "TestVar" ["dev1"], 1:1, setup.var_container, - 0.0, - 4.0, + [(min = 0.0, max = 4.0)], TEST_META, ) expr_container = IOM.get_expression( @@ -482,8 +240,7 @@ const TEST_META = "TestVar" ["dev1"], 1:1, setup.var_container, - 0.0, - 4.0, + [(min = 0.0, max = 4.0)], TEST_META, ) expr_container = IOM.get_expression( @@ -504,52 +261,6 @@ const TEST_META = "TestVar" @test JuMP.value(y) ≈ 1.0 atol = 1e-6 end - @testset "Multiple time steps" begin - setup = _setup_qa_test(["dev1"], 1:3) - IOM._add_quadratic_approx!( - IOM.SawtoothQuadConfig(2), - setup.container, - MockThermalGen, - ["dev1"], - 1:3, - setup.var_container, - 0.0, - 4.0, - TEST_META, - ) - - # Verify variables exist for each time step - g_container = IOM.get_variable( - setup.container, - IOM.SawtoothAuxVariable, - MockThermalGen, - TEST_META, - ) - alpha_container = IOM.get_variable( - setup.container, - IOM.SawtoothBinaryVariable, - MockThermalGen, - TEST_META, - ) - for t in 1:3, j in 0:2 - @test JuMP.lower_bound(g_container["dev1", j, t]) == 0.0 - end - for t in 1:3, j in 1:2 - @test JuMP.is_binary(alpha_container["dev1", j, t]) - end - - # Expression container should have entries for all (name, t) pairs - expr_container = IOM.get_expression( - setup.container, - IOM.QuadraticExpression, - MockThermalGen, - TEST_META, - ) - for t in 1:3 - @test expr_container["dev1", t] isa JuMP.AffExpr - end - end - @testset "Approximation quality improves with depth" begin # min (√2)x² - (√3)x on [0, 6], analytic minimum at x=√3/8 analytic_min = sqrt(2) * (sqrt(3 / 8)^2) - sqrt(3) * sqrt(3 / 8) @@ -567,8 +278,7 @@ const TEST_META = "TestVar" ["dev1"], 1:1, setup.var_container, - 0.0, - 6.0, + [(min = 0.0, max = 6.0)], TEST_META, ) expr_container = IOM.get_expression( @@ -610,8 +320,7 @@ const TEST_META = "TestVar" ["dev1"], 1:1, setup.var_container, - 0.0, - 4.0, + [(min = 0.0, max = 4.0)], TEST_META, ) else @@ -622,8 +331,7 @@ const TEST_META = "TestVar" ["dev1"], 1:1, setup.var_container, - 0.0, - 4.0, + [(min = 0.0, max = 4.0)], TEST_META, ) end From 4a75aa4f1defc116257299081acf8bed7cd75ccd Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Wed, 22 Apr 2026 13:13:35 -0400 Subject: [PATCH 09/17] format --- test/test_nmdt_approximations.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/test/test_nmdt_approximations.jl b/test/test_nmdt_approximations.jl index a428faf2..a5c2ecfa 100644 --- a/test/test_nmdt_approximations.jl +++ b/test/test_nmdt_approximations.jl @@ -98,7 +98,6 @@ const NMDT_BILINEAR_META = "NMDTBilinearTest" @test maximum(gaps) <= theoretical_bound + 1e-6 end end - end @testset "T-D-NMDT Tightening" begin @@ -274,7 +273,6 @@ end end end - # @testset "McCormick toggle" begin # setup = _setup_bilinear_test(["dev1"], 1:1) @@ -315,7 +313,6 @@ end @test JuMP.objective_value(setup.jump_model) ≈ 6.0 atol = 0.5 end - @testset "D-NMDT vs HybS comparison" begin true_product = 0.4 * 0.7 for depth in [2, 3] @@ -634,7 +631,8 @@ end IOM.NMDTBilinearConfig(L), setup.container, MockThermalGen, ["dev1"], 1:1, setup.x_var_container, setup.y_var_container, - [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], NMDT_BILINEAR_META, + [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], + NMDT_BILINEAR_META, ) expr = IOM.get_expression( setup.container, IOM.BilinearProductExpression, From cefadf13619007de2af5bff78f612a419bb14b5e Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Wed, 22 Apr 2026 13:43:36 -0600 Subject: [PATCH 10/17] add stubs for 2 edge cases where need PSY bus types --- src/common_models/add_constraint_dual.jl | 3 ++- src/common_models/interfaces.jl | 12 ++++++++++++ src/quadratic_approximations/incremental.jl | 8 ++++++-- 3 files changed, 20 insertions(+), 3 deletions(-) diff --git a/src/common_models/add_constraint_dual.jl b/src/common_models/add_constraint_dual.jl index 643d9e58..fa56fc3e 100644 --- a/src/common_models/add_constraint_dual.jl +++ b/src/common_models/add_constraint_dual.jl @@ -20,7 +20,8 @@ function add_constraint_dual!( model::NetworkModel{T}, ) where {T <: AbstractPowerModel} if !isempty(get_duals(model)) - devices = get_available_components(model, IS.InfrastructureSystemsComponent, sys) + # component is ACBus, but we don't have PSY as a dependency. + devices = get_available_components(model, component_for_network_dual(nothing), sys) for constraint_type in get_duals(model) assign_dual_variable!(container, constraint_type, devices, model) end diff --git a/src/common_models/interfaces.jl b/src/common_models/interfaces.jl index 749a43a6..57879f34 100644 --- a/src/common_models/interfaces.jl +++ b/src/common_models/interfaces.jl @@ -245,3 +245,15 @@ Only called in `emulation_model.jl`: that file's contents and this function shou likely be moved to POM or PSI. """ function update_container_parameter_values! end + +""" +Component type associated with network duals. Returns ACBus when passed in `nothing`. +Only called in a single spot in `add_constraint_dual!` for the network model. +""" +function component_for_network_dual end + +""" +Component type associated with hvdc interpolation constraints. Returns DCBus when passed in +`nothing`. Only called in a single spot in `incremental.jl`. +""" +function component_for_hvdc_interpolation end diff --git a/src/quadratic_approximations/incremental.jl b/src/quadratic_approximations/incremental.jl index 01f31e0a..527256ca 100644 --- a/src/quadratic_approximations/incremental.jl +++ b/src/quadratic_approximations/incremental.jl @@ -165,8 +165,12 @@ function _add_generic_incremental_interpolation_constraint!( # Retrieve all required variables from the optimization container # Retrieve original variable for DCVoltage from the Bus - # FIXME edge case: DCVoltage. W is InterconnectingConverter but key says DCBus. - x_var = get_variable(container, R, W) # Original variable (domain of function) + if R <: DCVoltage + # workaround for the fact that we can't write PSY.DCBus. + x_var = get_variable(container, R, component_for_hvdc_interpolation(nothing)) + else + x_var = get_variable(container, R, W) # Original variable (domain of function) + end y_var = get_variable(container, S, W) # Approximated variable (range of function) δ_var = get_variable(container, T, W) # Interpolation variables (weights for segments) z_var = get_variable(container, U, W) # Binary variables (ordering constraints) From 46c5204ea62fa4d10cc46bb88062dabb64c37ec2 Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Wed, 22 Apr 2026 17:53:19 -0600 Subject: [PATCH 11/17] Route system queries through extension-point stubs (fixes #87) IOM now declares stubs for the system-query interface it needs and calls them unqualified; downstream packages (POM) provide methods that forward to their system type's public API. No more sys.data access in IOM src. Co-Authored-By: Claude Opus 4.7 (1M context) --- src/common_models/interfaces.jl | 34 ++++++++++++ src/core/optimization_problem_outputs.jl | 3 +- src/core/settings.jl | 3 +- src/operation/decision_model.jl | 13 ++--- src/operation/emulation_model.jl | 7 +-- src/utils/component_utils.jl | 16 ++---- test/mocks/mock_system.jl | 71 +++++++----------------- 7 files changed, 70 insertions(+), 77 deletions(-) diff --git a/src/common_models/interfaces.jl b/src/common_models/interfaces.jl index 57879f34..1fa1ad88 100644 --- a/src/common_models/interfaces.jl +++ b/src/common_models/interfaces.jl @@ -145,6 +145,40 @@ For thermals, equivalent to `get_must_run`, but that implementation belongs in P """ skip_proportional_cost(d::IS.InfrastructureSystemsComponent) = false +############################### +#### System query stubs ####### +############################### +# Extension points for querying a system object. POM provides methods for +# PSY.System; tests provide methods for MockSystem. IOM itself never accesses +# sys.data. + +"Extension point: time-series resolutions available on the system." +function get_time_series_resolutions end + +"Extension point: counts summary of time series on the system." +function get_time_series_counts end + +"Extension point: counts by component type of time series on the system." +function get_time_series_counts_by_type end + +"Extension point: forecast interval configured on the system." +function get_forecast_interval end + +"Extension point: forecast horizon configured on the system." +function get_forecast_horizon end + +"Extension point: summary table of forecasts on the system." +function get_forecast_summary_table end + +"Extension point: transform single time series into deterministic forecasts on the system." +function transform_single_time_series! end + +"Extension point: stable UUID for the system (used as a filename identifier)." +function get_system_uuid end + +"Extension point: get components of type `T` in a subsystem of the system." +function get_subsystem_components end + ############################### ###### Start-up Cost ########## ############################### diff --git a/src/core/optimization_problem_outputs.jl b/src/core/optimization_problem_outputs.jl index da4261ec..2db94f1d 100644 --- a/src/core/optimization_problem_outputs.jl +++ b/src/core/optimization_problem_outputs.jl @@ -104,9 +104,8 @@ get_optimizer_stats(res::OptimizationProblemOutputs) = res.optimizer_stats get_parameter_values(res::OptimizationProblemOutputs) = res.parameter_values get_source_data(res::OptimizationProblemOutputs) = res.source_data -# FIXME get_uuid declare as stub make_system_filename(sys::IS.InfrastructureSystemsContainer) = - make_system_filename(IS.get_uuid(sys.data.internal)) + make_system_filename(get_system_uuid(sys)) make_system_filename(sys_uuid::Union{Base.UUID, AbstractString}) = "system-$(sys_uuid).json" """ diff --git a/src/core/settings.jl b/src/core/settings.jl index 37f9be9a..58c38c15 100644 --- a/src/core/settings.jl +++ b/src/core/settings.jl @@ -59,8 +59,7 @@ function Settings( store_variable_names = false, ext = Dict{String, Any}(), ) - # alternatively, declare as stub and implement in POM. - if time_series_cache_size > 0 && IS.stores_time_series_in_memory(sys.data) + if time_series_cache_size > 0 && stores_time_series_in_memory(sys) @info "Overriding time_series_cache_size because time series is stored in memory" time_series_cache_size = 0 end diff --git a/src/operation/decision_model.jl b/src/operation/decision_model.jl index 854a5f8a..5585c422 100644 --- a/src/operation/decision_model.jl +++ b/src/operation/decision_model.jl @@ -1,5 +1,5 @@ function get_deterministic_time_series_type(sys::IS.InfrastructureSystemsContainer) - time_series_types = IS.get_time_series_counts_by_type(sys.data) + time_series_types = get_time_series_counts_by_type(sys) existing_types = Set(d["type"] for d in time_series_types) if ("Deterministic" in existing_types) && ("DeterministicSingleTimeSeries" in existing_types) @@ -242,12 +242,11 @@ function init_model_store_params!(model::DecisionModel) if model_interval != UNSET_INTERVAL interval = model_interval else - interval = IS.get_forecast_interval(system.data) + interval = get_forecast_interval(system) end resolution = get_resolution(model) base_power = get_base_power(system) - # FIXME declare as stub - sys_uuid = IS.get_uuid(system.data.internal) + sys_uuid = get_system_uuid(system) store_params = ModelStoreParams( num_executions, horizon, @@ -264,7 +263,7 @@ end function validate_time_series!(model::DecisionModel{<:DefaultDecisionProblem}) sys = get_system(model) settings = get_settings(model) - available_resolutions = IS.get_time_series_resolutions(sys.data) + available_resolutions = get_time_series_resolutions(sys) if get_resolution(settings) == UNSET_RESOLUTION && length(available_resolutions) != 1 throw( @@ -307,11 +306,11 @@ function validate_time_series!(model::DecisionModel{<:DefaultDecisionProblem}) if get_horizon(settings) == UNSET_HORIZON set_horizon!( settings, - IS.get_forecast_horizon(sys.data; interval = _to_is_interval(model_interval)), + get_forecast_horizon(sys; interval = _to_is_interval(model_interval)), ) end - counts = IS.get_time_series_counts(sys.data) + counts = get_time_series_counts(sys) if counts.forecast_count < 1 error( "The system does not contain forecast data. A DecisionModel can't be built.", diff --git a/src/operation/emulation_model.jl b/src/operation/emulation_model.jl index 33b04284..1898e5b0 100644 --- a/src/operation/emulation_model.jl +++ b/src/operation/emulation_model.jl @@ -211,7 +211,7 @@ end function validate_time_series!(model::EmulationModel{<:DefaultEmulationProblem}) sys = get_system(model) settings = get_settings(model) - available_resolutions = IS.get_time_series_resolutions(sys.data) + available_resolutions = get_time_series_resolutions(sys) if get_resolution(settings) == UNSET_RESOLUTION && length(available_resolutions) != 1 throw( @@ -236,7 +236,7 @@ function validate_time_series!(model::EmulationModel{<:DefaultEmulationProblem}) set_horizon!(settings, get_resolution(settings)) end - counts = IS.get_time_series_counts(sys.data) + counts = get_time_series_counts(sys) if counts.static_time_series_count < 1 error( "The system does not contain Static Time Series data. A EmulationModel can't be built.", @@ -258,8 +258,7 @@ function init_model_store_params!(model::EmulationModel) settings = get_settings(model) horizon = interval = resolution = get_resolution(settings) base_power = get_base_power(system) - # FIXME declare as stub - sys_uuid = IS.get_uuid(system.data.internal) + sys_uuid = get_system_uuid(system) set_store_params!( get_internal(model), ModelStoreParams( diff --git a/src/utils/component_utils.jl b/src/utils/component_utils.jl index 1d240e6b..6dba8400 100644 --- a/src/utils/component_utils.jl +++ b/src/utils/component_utils.jl @@ -66,12 +66,7 @@ function get_available_components( sys::IS.InfrastructureSystemsContainer, ) where {T <: IS.InfrastructureSystemsComponent} subsystem = get_subsystem(model) - # FIXME have to patch thru to sys.data here - return IS.get_components( - T, - sys.data; - subsystem_name = subsystem, - ) + return get_subsystem_components(T, sys; subsystem_name = subsystem) end ################################################## @@ -305,7 +300,7 @@ end is_time_variant(x) = IS.is_time_series_backed(x) function get_forecast_intervals(sys::IS.InfrastructureSystemsContainer) - table = IS.get_forecast_summary_table(sys.data) + table = get_forecast_summary_table(sys) return Set(row.interval for row in eachrow(table) if row.interval !== nothing) end @@ -319,7 +314,7 @@ function auto_transform_time_series!( return end - counts = IS.get_time_series_counts(sys.data) + counts = get_time_series_counts(sys) if counts.static_time_series_count < 1 return end @@ -333,9 +328,8 @@ function auto_transform_time_series!( @info "Auto-transforming SingleTimeSeries to DeterministicSingleTimeSeries" horizon = Dates.canonicalize(model_horizon) interval = Dates.canonicalize(model_interval) - IS.transform_single_time_series!( - sys.data, - IS.DeterministicSingleTimeSeries, + transform_single_time_series!( + sys, model_horizon, model_interval; delete_existing = false, diff --git a/test/mocks/mock_system.jl b/test/mocks/mock_system.jl index 2e375e8e..74b34a88 100644 --- a/test/mocks/mock_system.jl +++ b/test/mocks/mock_system.jl @@ -3,70 +3,39 @@ Minimal mock for PSY.System. Implements only the interface required by OptimizationContainer and models. """ -#= -function stores_time_series_in_memory end -function get_time_series_counts_by_type end -function get_time_series_counts end -function get_forecast_interval end -function get_time_series_resolutions end -function get_forecast_horizon end -function get_uuid end -=# - using InfrastructureSystems const IS = InfrastructureSystems -mutable struct MockSystemData +mutable struct MockSystem <: IS.InfrastructureSystemsContainer base_power::Float64 components::Dict{DataType, Vector{Any}} time_series::Dict{Any, Any} stores_in_memory::Bool - internal::IS.InfrastructureSystemsInternal -end - -mutable struct MockSystem <: IS.InfrastructureSystemsContainer - data::MockSystemData end -MockSystem() = MockSystem(MockSystemData()) -MockSystem(base_power::Float64) = MockSystem(MockSystemData(base_power)) -MockSystem(base_power::Float64, stores_in_memory::Bool) = - MockSystem(MockSystemData(base_power, stores_in_memory)) - # Convenience constructors -MockSystemData() = MockSystemData( - 100.0, - Dict{DataType, Vector{Any}}(), - Dict{Any, Any}(), - false, - IS.InfrastructureSystemsInternal(), -) -MockSystemData(base_power::Float64) = MockSystemData( - base_power, - Dict{DataType, Vector{Any}}(), - Dict{Any, Any}(), - false, - IS.InfrastructureSystemsInternal(), -) -MockSystemData(base_power::Float64, stores_in_memory::Bool) = MockSystemData( - base_power, - Dict{DataType, Vector{Any}}(), - Dict{Any, Any}(), - stores_in_memory, - IS.InfrastructureSystemsInternal(), -) +MockSystem() = MockSystem(100.0, Dict{DataType, Vector{Any}}(), Dict{Any, Any}(), false) +MockSystem(base_power::Float64) = + MockSystem(base_power, Dict{DataType, Vector{Any}}(), Dict{Any, Any}(), false) +MockSystem(base_power::Float64, stores_in_memory::Bool) = + MockSystem( + base_power, + Dict{DataType, Vector{Any}}(), + Dict{Any, Any}(), + stores_in_memory, + ) # Required interface methods - extend InfrastructureOptimizationModels functions for duck-typing -IOM.get_base_power(sys::MockSystem) = sys.data.base_power -IOM.stores_time_series_in_memory(sys::MockSystem) = sys.data.stores_in_memory -IOM.stores_time_series_in_memory(data::MockSystemData) = data.stores_in_memory +InfrastructureOptimizationModels.get_base_power(sys::MockSystem) = sys.base_power +InfrastructureOptimizationModels.stores_time_series_in_memory(sys::MockSystem) = + sys.stores_in_memory function IOM.get_available_components(::NetworkModel, ::Type{T}, sys::MockSystem) where {T} return get_components(T, sys) end function get_components(::Type{T}, sys::MockSystem) where {T} - return get(sys.data.components, T, T[]) + return get(sys.components, T, T[]) end function get_component(::Type{T}, sys::MockSystem, name::String) where {T} @@ -80,10 +49,10 @@ end function add_component!(sys::MockSystem, component) comp_type = typeof(component) - if !haskey(sys.data.components, comp_type) - sys.data.components[comp_type] = [] + if !haskey(sys.components, comp_type) + sys.components[comp_type] = [] end - push!(sys.data.components[comp_type], component) + push!(sys.components[comp_type], component) return end @@ -94,10 +63,10 @@ function get_time_series( args...; kwargs..., ) where {T} - return get(sys.data.time_series, (T, component), nothing) + return get(sys.time_series, (T, component), nothing) end function add_time_series!(sys::MockSystem, component, ts) - sys.data.time_series[(typeof(ts), component)] = ts + sys.time_series[(typeof(ts), component)] = ts return end From c87697b449622304aeee655289c4a6b8dbe30de4 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Thu, 23 Apr 2026 02:46:17 -0400 Subject: [PATCH 12/17] restore aqua tests --- test/runtests.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index c0e37618..91774c12 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,10 +9,10 @@ using InfrastructureSystems const IS = InfrastructureSystems # Code Quality Tests -# import Aqua -# @testset "Code Quality (Aqua.jl)" begin -# Aqua.test_all(InfrastructureOptimizationModels; persistent_tasks = false) -# end +import Aqua +@testset "Code Quality (Aqua.jl)" begin + Aqua.test_all(InfrastructureOptimizationModels; persistent_tasks = false) +end # Load the test module include("InfrastructureOptimizationModelsTests.jl") From cd4f824ebd4185c36f1e9dae7575ea1c87ac78cf Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Thu, 23 Apr 2026 11:32:32 -0400 Subject: [PATCH 13/17] resolved ambiguity --- src/utils/jump_utils.jl | 58 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/src/utils/jump_utils.jl b/src/utils/jump_utils.jl index 80865062..5b7be29d 100644 --- a/src/utils/jump_utils.jl +++ b/src/utils/jump_utils.jl @@ -428,6 +428,64 @@ function to_outputs_dataframe( ) end +function to_outputs_dataframe( + array::DenseAxisArray{ + Float64, + 3, + <:Tuple{Vector{String}, Vector{String}, UnitRange{Int}}, + }, + ::Nothing, + ::Val{TableFormat.LONG}, +) + return to_outputs_dataframe(array, nothing, TableFormat.LONG) +end + +function to_outputs_dataframe( + array::DenseAxisArray{ + Float64, + 3, + <:Tuple{Vector{String}, T, UnitRange{Int}}, + }, + timestamps, + ::Val{TableFormat.LONG}, +) where {T <: Union{Vector{String}, UnitRange{Int}}} + num_timestamps = size(array, 3) + if length(timestamps) != num_timestamps + error( + "The number of timestamps must match the number of rows per component. " * + "timestamps = $(length(timestamps)) " * + "num_timestamps = $num_timestamps", + ) + end + + num_rows = length(array.data) + timestamps_arr = _collect_timestamps(timestamps) + time_col = Vector{Int}(undef, num_rows) + name_col = Vector{eltype(T)}(undef, num_rows) + name2_col = Vector{String}(undef, num_rows) + vals = Vector{Float64}(undef, num_rows) + + row_index = 1 + for name in axes(array, 1) + for name2 in axes(array, 2) + for time_index in axes(array, 3) + time_col[row_index] = timestamps_arr[time_index] + name_col[row_index] = name + name2_col[row_index] = name2 + vals[row_index] = array[name, name2, time_index] + row_index += 1 + end + end + end + + return DataFrame( + :time_index => time_col, + :name => name_col, + :name2 => name2_col, + :value => vals, + ) +end + function to_dataframe(array::SparseAxisArray{T, N, K}) where {T, N, K <: NTuple{N, Any}} columns = get_column_names_from_axis_array(array) return DataFrames.DataFrame(_to_matrix(array, columns), columns) From 6ffcbe7e64deb7ddde5a7ff0d9ed44fdda12db99 Mon Sep 17 00:00:00 2001 From: Luke Kiernan Date: Thu, 23 Apr 2026 15:17:08 -0600 Subject: [PATCH 14/17] helpers for hydro/storage POM fixes --- src/InfrastructureOptimizationModels.jl | 2 +- src/common_models/range_constraint.jl | 25 +++++++++++++++++++--- src/objective_function/proportional.jl | 28 +++++++++++++------------ 3 files changed, 38 insertions(+), 17 deletions(-) diff --git a/src/InfrastructureOptimizationModels.jl b/src/InfrastructureOptimizationModels.jl index 4e11e34b..1fa5f860 100644 --- a/src/InfrastructureOptimizationModels.jl +++ b/src/InfrastructureOptimizationModels.jl @@ -348,7 +348,7 @@ export add_sparse_pwl_interpolation_variables! export JuMPOrFloat # Constraint helpers export add_range_constraints!, add_parameterized_upper_bound_range_constraints -export add_reserve_bound_range_constraints! +export add_reserve_bound_range_constraints!, add_commitment_bound_range_constraints! export add_semicontinuous_range_constraints!, add_semicontinuous_ramp_constraints! # Cost helpers export add_shut_down_cost!, add_start_up_cost! diff --git a/src/common_models/range_constraint.jl b/src/common_models/range_constraint.jl index 2b10439a..8a141acf 100644 --- a/src/common_models/range_constraint.jl +++ b/src/common_models/range_constraint.jl @@ -195,14 +195,18 @@ function add_semicontinuous_range_constraints!( return end -# Generic component version - always uses binary variable +# Generic component version - always uses binary variable. +# `meta_suffix` is appended to the default constraint meta so callers can stack a second +# OnVariable-keyed bound alongside another bound constraint (e.g. a reservation-keyed +# one) on the same `(T, V)` without a meta collision — see `add_commitment_bound_range_constraints!`. function _add_semicontinuous_bound_range_constraints_impl!( container::OptimizationContainer, ::Type{T}, dir::BoundDirection, array, devices::Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, - ::DeviceModel{V, W}, + ::DeviceModel{V, W}; + meta_suffix::String = "", ) where { T <: ConstraintType, V <: IS.InfrastructureSystemsComponent, @@ -212,7 +216,7 @@ function _add_semicontinuous_bound_range_constraints_impl!( names = IS.get_name.(devices) jump_model = get_jump_model(container) con = add_constraints_container!( - container, T, V, names, time_steps; meta = constraint_meta(dir)) + container, T, V, names, time_steps; meta = constraint_meta(dir) * meta_suffix) varbin = get_variable(container, OnVariable, V) for device in devices, t in time_steps @@ -225,6 +229,21 @@ function _add_semicontinuous_bound_range_constraints_impl!( return end +# Exported wrapper: use this from downstream packages to add an OnVariable-keyed bound +# alongside another bound constraint on the same `(T, V)` key — pass `meta_suffix = "_aux"` +# (or similar) to avoid colliding with the default "lb"/"ub" meta. +add_commitment_bound_range_constraints!( + container::OptimizationContainer, + ::Type{T}, + dir::BoundDirection, + array, + devices, + model::DeviceModel; + meta_suffix::String = "", +) where {T <: ConstraintType} = + _add_semicontinuous_bound_range_constraints_impl!( + container, T, dir, array, devices, model; meta_suffix) + # Unified reserve range constraints impl # invert_binary: true for InputActivePower (uses 1-varbin), false for others (uses varbin) function add_reserve_bound_range_constraints!( diff --git a/src/objective_function/proportional.jl b/src/objective_function/proportional.jl index 0860cebf..0abb7cae 100644 --- a/src/objective_function/proportional.jl +++ b/src/objective_function/proportional.jl @@ -5,8 +5,9 @@ # this is only used for ControllableLoads with non-PowerLoadInterruptible formulations. # The rest go through a thin wrapper around the maybe-variant version. """ -Default implementation for proportional cost, where the cost term is not time variant. Anything -time-varying should implement its own method. +Default implementation for proportional cost, where the cost term is not time variant. +See also: `add_proportional_cost_maybe_time_variant!` for a common basis for devices that +might have time-variant proportional costs. """ function add_proportional_cost!( container::OptimizationContainer, @@ -18,7 +19,6 @@ function add_proportional_cost!( U <: VariableType, V <: AbstractDeviceFormulation, } - # NOTE: anything time-varying should implement its own method. multiplier = objective_function_multiplier(U, V) for d in devices op_cost_data = get_operation_cost(d) @@ -26,17 +26,19 @@ function add_proportional_cost!( iszero(cost_term) && continue name = get_name(d) rate = cost_term * multiplier + skip = skip_proportional_cost(d) 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, - ) + if skip + # must-run etc.: bookkeep in ProductionCostExpression but not in objective + add_cost_to_expression!( + container, ProductionCostExpression, rate, T, name, t) + else + variable = get_variable(container, U, T)[name, t] + add_cost_term_invariant!( + container, variable, rate, + ProductionCostExpression, T, name, t, + ) + end end end return From 9ccee6240f9b07ee16bf5a23b144e2a9e4fc4459 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Fri, 24 Apr 2026 13:07:05 -0400 Subject: [PATCH 15/17] check_time_series_consistency --- src/core/optimization_container.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/optimization_container.jl b/src/core/optimization_container.jl index 7dec468b..36f1979b 100644 --- a/src/core/optimization_container.jl +++ b/src/core/optimization_container.jl @@ -309,7 +309,7 @@ function init_optimization_container!( # set_initial_time!(settings, IS.get_forecast_initial_timestamp(sys)) set_initial_time!(settings, temp_get_forecast_initial_timestamp(sys)) elseif get_default_time_series_type(container) <: IS.SingleTimeSeries - ini_time, _ = IS.check_time_series_consistency(sys, IS.SingleTimeSeries) + ini_time, _ = IS.check_time_series_consistency(sys.data, IS.SingleTimeSeries) set_initial_time!(settings, ini_time) end end From 00ff30993461c207cfecd3cc5c0f08fe7f94f4e1 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Mon, 27 Apr 2026 14:01:56 -0400 Subject: [PATCH 16/17] bring back old constraint fallback -- still needed in addition to new ones --- src/operation/model_numerical_analysis_utils.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/operation/model_numerical_analysis_utils.jl b/src/operation/model_numerical_analysis_utils.jl index a023d49e..b27324a9 100644 --- a/src/operation/model_numerical_analysis_utils.jl +++ b/src/operation/model_numerical_analysis_utils.jl @@ -99,6 +99,7 @@ function update_numerical_bounds(bonuds::NumericalBounds, func::MOI.Interval, id end # Default fallbacks for unsupported constraints. +update_numerical_bounds(::NumericalBounds, func, idx) = nothing update_coefficient_bounds(::ConstraintBounds, func, idx) = nothing update_rhs_bounds(::ConstraintBounds, func, idx) = nothing From 10812a6b018077d9a7375aeba0c3e6980b88d19a Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Mon, 27 Apr 2026 14:11:44 -0600 Subject: [PATCH 17/17] simplify to_outputs_dataframe dispatch --- src/utils/jump_utils.jl | 62 ++--------------------------------------- 1 file changed, 2 insertions(+), 60 deletions(-) diff --git a/src/utils/jump_utils.jl b/src/utils/jump_utils.jl index 5b7be29d..a87ab5dc 100644 --- a/src/utils/jump_utils.jl +++ b/src/utils/jump_utils.jl @@ -403,8 +403,8 @@ function to_outputs_dataframe( ) where {T <: Union{Vector{String}, UnitRange{Int}}} num_rows = length(array.data) time_col = Vector{Int}(undef, num_rows) - name_col = Vector{eltype(T)}(undef, num_rows) - name2_col = Vector{String}(undef, num_rows) + name_col = Vector{String}(undef, num_rows) + name2_col = Vector{eltype(T)}(undef, num_rows) vals = Vector{Float64}(undef, num_rows) row_index = 1 @@ -428,64 +428,6 @@ function to_outputs_dataframe( ) end -function to_outputs_dataframe( - array::DenseAxisArray{ - Float64, - 3, - <:Tuple{Vector{String}, Vector{String}, UnitRange{Int}}, - }, - ::Nothing, - ::Val{TableFormat.LONG}, -) - return to_outputs_dataframe(array, nothing, TableFormat.LONG) -end - -function to_outputs_dataframe( - array::DenseAxisArray{ - Float64, - 3, - <:Tuple{Vector{String}, T, UnitRange{Int}}, - }, - timestamps, - ::Val{TableFormat.LONG}, -) where {T <: Union{Vector{String}, UnitRange{Int}}} - num_timestamps = size(array, 3) - if length(timestamps) != num_timestamps - error( - "The number of timestamps must match the number of rows per component. " * - "timestamps = $(length(timestamps)) " * - "num_timestamps = $num_timestamps", - ) - end - - num_rows = length(array.data) - timestamps_arr = _collect_timestamps(timestamps) - time_col = Vector{Int}(undef, num_rows) - name_col = Vector{eltype(T)}(undef, num_rows) - name2_col = Vector{String}(undef, num_rows) - vals = Vector{Float64}(undef, num_rows) - - row_index = 1 - for name in axes(array, 1) - for name2 in axes(array, 2) - for time_index in axes(array, 3) - time_col[row_index] = timestamps_arr[time_index] - name_col[row_index] = name - name2_col[row_index] = name2 - vals[row_index] = array[name, name2, time_index] - row_index += 1 - end - end - end - - return DataFrame( - :time_index => time_col, - :name => name_col, - :name2 => name2_col, - :value => vals, - ) -end - function to_dataframe(array::SparseAxisArray{T, N, K}) where {T, N, K <: NTuple{N, Any}} columns = get_column_names_from_axis_array(array) return DataFrames.DataFrame(_to_matrix(array, columns), columns)