diff --git a/src/InfrastructureOptimizationModels.jl b/src/InfrastructureOptimizationModels.jl index b7b500e2..f0b44ffb 100644 --- a/src/InfrastructureOptimizationModels.jl +++ b/src/InfrastructureOptimizationModels.jl @@ -624,25 +624,27 @@ include("objective_function/objective_function_pwl_delta.jl") # delta/increment include("objective_function/piecewise_linear.jl") # CostCurve/FuelCurve → lambda PWL include("objective_function/value_curve_cost.jl") # ValueCurve → delta PWL -# Quadratic approximations (PWL via SOS2) -include("quadratic_approximations/common.jl") -include("quadratic_approximations/no_approx.jl") -include("quadratic_approximations/pwl_utils.jl") -include("quadratic_approximations/incremental.jl") -include("quadratic_approximations/solver_sos2.jl") -include("quadratic_approximations/manual_sos2.jl") -include("quadratic_approximations/sawtooth.jl") -include("quadratic_approximations/epigraph.jl") -include("quadratic_approximations/nmdt_common.jl") -include("quadratic_approximations/nmdt.jl") -include("quadratic_approximations/pwmcc_cuts.jl") - -# Bilinear approximations (x·y via Bin2/HybS decomposition) -include("bilinear_approximations/mccormick.jl") -include("bilinear_approximations/bin2.jl") -include("bilinear_approximations/no_approx.jl") -include("bilinear_approximations/hybs.jl") -include("bilinear_approximations/nmdt.jl") +# Quadratic and bilinear approximations. +# Each method ships a scalar `build_*` (pure JuMP) and an `add_*_approx!` +# IOM adapter (allocate, loop, write) in the same file. +include("approximations/common.jl") +include("approximations/pwl_utils.jl") +include("approximations/mccormick.jl") +include("approximations/pwmcc_cuts.jl") +include("approximations/epigraph.jl") # must precede sawtooth and NMDT (tightening) +include("approximations/nmdt_discretization.jl") # must precede NMDT quad/bilinear +# Quadratic methods (each file is self-contained: config + scalar build + IOM adapter) +include("approximations/no_approx_quadratic.jl") +include("approximations/solver_sos2.jl") +include("approximations/manual_sos2.jl") +include("approximations/sawtooth.jl") +include("approximations/nmdt_quadratic.jl") +include("approximations/incremental.jl") +# Bilinear methods (compose with quadratic — must follow) +include("approximations/no_approx_bilinear.jl") +include("approximations/nmdt_bilinear.jl") +include("approximations/bin2.jl") +include("approximations/hybs.jl") # add_param_container! wrappers — must come after piecewise_linear.jl # (which defines VariableValueParameter and FixValueParameter) diff --git a/src/approximations/bin2.jl b/src/approximations/bin2.jl new file mode 100644 index 00000000..51c35642 --- /dev/null +++ b/src/approximations/bin2.jl @@ -0,0 +1,226 @@ +# Bin2 separable approximation of bilinear products z = x·y. +# Uses the identity x·y = ½·((x+y)² − x² − y²). +# Composes a quadratic approximation (chosen via `quad_config`) for x², y², +# and (x+y)². Optionally adds reformulated McCormick cuts to tighten the LP +# relaxation in terms of the three quadratic approximations. + +""" +Config for Bin2 bilinear approximation using z = ½·((x+y)² − x² − y²). + +# Fields +- `quad_config::QuadraticApproxConfig`: quadratic method used for x², y², and (x+y)². +- `add_mccormick::Bool`: whether to add reformulated McCormick cuts (default true). +""" +struct Bin2Config{QC <: QuadraticApproxConfig} <: BilinearApproxConfig + quad_config::QC + add_mccormick::Bool +end +function Bin2Config(quad_config::QuadraticApproxConfig) + return Bin2Config(quad_config, true) +end + +""" + build_bilinear_approx(config::Bin2Config, model, x, y, x_min, x_max, y_min, y_max) + +Scalar form: build x², y², (x+y)² via the chosen quadratic method, combine +via z = ½·(psq − xsq − ysq). If `config.add_mccormick`, also build the +four reformulated McCormick cuts. + +Returns `(; approximation, xsq, ysq, psq, sum_expression, mccormick_constraints)` +where `mccormick_constraints` is `nothing` or a NamedTuple `(c1, c2, c3, c4)`. +""" +function build_bilinear_approx( + config::Bin2Config, + model::JuMP.Model, + x::JuMP.AbstractJuMPScalar, + y::JuMP.AbstractJuMPScalar, + x_min::Float64, + x_max::Float64, + y_min::Float64, + y_max::Float64, +) + xsq = build_quadratic_approx(config.quad_config, model, x, x_min, x_max) + ysq = build_quadratic_approx(config.quad_config, model, y, y_min, y_max) + p_expr = JuMP.@expression(model, x + y) + psq = build_quadratic_approx( + config.quad_config, model, p_expr, x_min + y_min, x_max + y_max, + ) + approximation = JuMP.@expression( + model, + 0.5 * (psq.approximation - xsq.approximation - ysq.approximation), + ) + mc = if config.add_mccormick + build_reformulated_mccormick( + model, x, y, + psq.approximation, xsq.approximation, ysq.approximation, + x_min, x_max, y_min, y_max, + ) + else + nothing + end + return (; + approximation, + xsq, + ysq, + psq, + sum_expression = p_expr, + mccormick_constraints = mc, + ) +end + +""" + add_bilinear_approx!(config::Bin2Config, container, ::Type{C}, x_var, y_var, x_bounds, y_bounds, meta) + +Build x² and y² via `add_quadratic_approx!(config.quad_config, ...)`, +build the (x+y) expression container and its psq via the same quad +adapter, then assemble z = ½·(psq − xsq − ysq) and (optionally) the +reformulated McCormick cuts. +""" +function add_bilinear_approx!( + config::Bin2Config, + container::OptimizationContainer, + ::Type{C}, + x_var, + y_var, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + name_axis = axes(x_var, 1) + time_axis = axes(x_var, 2) + @assert length(name_axis) == length(x_bounds) + @assert length(name_axis) == length(y_bounds) + model = get_jump_model(container) + + p_target = add_expression_container!( + container, VariableSumExpression, C, name_axis, time_axis; + meta = meta * "_plus", + ) + for (i, name) in enumerate(name_axis) + for t in time_axis + p_target[name, t] = x_var[name, t] + y_var[name, t] + end + end + p_bounds = [ + (min = x_bounds[i].min + y_bounds[i].min, + max = x_bounds[i].max + y_bounds[i].max) + for i in eachindex(x_bounds) + ] + + xsq = add_quadratic_approx!( + config.quad_config, + container, + C, + x_var, + x_bounds, + meta * "_x", + ) + ysq = add_quadratic_approx!( + config.quad_config, + container, + C, + y_var, + y_bounds, + meta * "_y", + ) + psq = add_quadratic_approx!( + config.quad_config, container, C, p_target, p_bounds, meta * "_plus", + ) + + return _bin2_assemble_and_mccormick!( + container, C, name_axis, time_axis, model, + x_var, y_var, xsq, ysq, psq, x_bounds, y_bounds, meta; + add_mccormick = config.add_mccormick, + ) +end + +""" + add_bilinear_approx!(config::Bin2Config, container, ::Type{C}, xsq, ysq, x_var, y_var, x_bounds, y_bounds, meta) + +Precomputed-form: accepts already-built `xsq` ≈ x² and `ysq` ≈ y² 2D +expression containers (rather than rebuilding them). Builds only the +(x+y)² approximation on top, the Bin2 assembly, and the optional cuts. +""" +function add_bilinear_approx!( + config::Bin2Config, + container::OptimizationContainer, + ::Type{C}, + xsq, + ysq, + x_var, + y_var, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + name_axis = axes(x_var, 1) + time_axis = axes(x_var, 2) + @assert length(name_axis) == length(x_bounds) + @assert length(name_axis) == length(y_bounds) + model = get_jump_model(container) + + p_target = add_expression_container!( + container, VariableSumExpression, C, name_axis, time_axis; + meta = meta * "_plus", + ) + for (i, name) in enumerate(name_axis) + for t in time_axis + p_target[name, t] = x_var[name, t] + y_var[name, t] + end + end + p_bounds = [ + (min = x_bounds[i].min + y_bounds[i].min, + max = x_bounds[i].max + y_bounds[i].max) + for i in eachindex(x_bounds) + ] + psq = add_quadratic_approx!( + config.quad_config, container, C, p_target, p_bounds, meta * "_plus", + ) + + return _bin2_assemble_and_mccormick!( + container, C, name_axis, time_axis, model, + x_var, y_var, xsq, ysq, psq, x_bounds, y_bounds, meta; + add_mccormick = config.add_mccormick, + ) +end + +# Allocate the BilinearProductExpression result + optional ReformulatedMcCormick +# container, then loop (name, t) to assemble z and the McCormick cuts. +function _bin2_assemble_and_mccormick!( + container, ::Type{C}, name_axis, time_axis, model, + x_var, y_var, xsq, ysq, psq, x_bounds, y_bounds, meta; + add_mccormick::Bool, +) where {C <: IS.InfrastructureSystemsComponent} + approx_target = add_expression_container!( + container, BilinearProductExpression, C, name_axis, time_axis; meta, + ) + mc_target = if add_mccormick + add_constraints_container!( + container, ReformulatedMcCormickConstraint, C, + name_axis, 1:4, time_axis; meta, + ) + else + nothing + end + + for (i, name) in enumerate(name_axis) + xmn, xmx = x_bounds[i].min, x_bounds[i].max + ymn, ymx = y_bounds[i].min, y_bounds[i].max + for t in time_axis + approx_target[name, t] = + 0.5 * (psq[name, t] - xsq[name, t] - ysq[name, t]) + if add_mccormick + r = build_reformulated_mccormick( + model, + x_var[name, t], y_var[name, t], + psq[name, t], xsq[name, t], ysq[name, t], + xmn, xmx, ymn, ymx, + ) + for (k, ref) in enumerate(r) + mc_target[name, k, t] = ref + end + end + end + end + return approx_target +end diff --git a/src/approximations/common.jl b/src/approximations/common.jl new file mode 100644 index 00000000..6026a588 --- /dev/null +++ b/src/approximations/common.jl @@ -0,0 +1,34 @@ +# Shared infrastructure for quadratic and bilinear approximations. +# +# Each method ships a scalar `build_` function (pure-JuMP, no IOM +# dependencies) and an `add__approx!` IOM adapter that allocates +# containers with known `(name, t, ...)` axes, loops over (name, t), calls +# the scalar build per cell, and slots refs into the container. + +# --- Abstract config supertypes --- + +"Abstract supertype for quadratic approximation method configurations." +abstract type QuadraticApproxConfig end + +"Abstract supertype for bilinear approximation method configurations." +abstract type BilinearApproxConfig end + +# --- Shared expression / variable key types --- + +"Expression container for the normalized variable xh = (x − x_min) / (x_max − x_min) ∈ [0,1]." +struct NormedVariableExpression <: ExpressionType end + +"Expression container for quadratic (x²) approximation results." +struct QuadraticExpression <: ExpressionType end + +"Expression container for bilinear product (x·y) approximation results." +struct BilinearProductExpression <: ExpressionType end + +"Variable container for bilinear product (x·y) approximation results." +struct BilinearProductVariable <: VariableType end + +"Expression container for sums of two variables, x + y." +struct VariableSumExpression <: ExpressionType end + +"Expression container for differences of two variables, x − y." +struct VariableDifferenceExpression <: ExpressionType end diff --git a/src/approximations/epigraph.jl b/src/approximations/epigraph.jl new file mode 100644 index 00000000..57e79a7f --- /dev/null +++ b/src/approximations/epigraph.jl @@ -0,0 +1,222 @@ +# Epigraph (Q^{L1}) LP-only lower bound for x² using tangent-line cuts. +# Pure LP — zero binary variables. Creates a variable z that lower-bounds +# x² (approximately) bounded from below by supporting hyperplanes of the +# parabola. +# Reference: Beach, Burlacu, Hager, Hildebrand (2024), Q^{L1} relaxation. + +"Expression container for epigraph quadratic approximation results." +struct EpigraphExpression <: ExpressionType end + +"Auxiliary continuous variables (g₀, …, g_L) for tooth-based PWL approximations." +struct SawtoothAuxVariable <: VariableType end + +"LP relaxation constraints (g_j ≤ 2g_{j-1}, g_j ≤ 2(1−g_{j-1}))." +struct SawtoothLPConstraint <: ConstraintType end + +"Links g₀ to the normalized x value." +struct SawtoothLinkingConstraint <: ConstraintType end + +"Variable representing a lower-bounded approximation of x² in epigraph relaxation." +struct EpigraphVariable <: VariableType end + +"Tangent-line lower-bound constraints in epigraph relaxation." +struct EpigraphTangentConstraint <: ConstraintType end + +"Tangent-line lower-bound expression fL used in the epigraph formulation." +struct EpigraphTangentExpression <: ExpressionType end + +""" + scale_back_g_basis(x_min, delta, g_var, levels) + +Build the affine "scale back to actual dimensions" expression at one cell: + x² ≈ x_min² + (2·x_min·δ + δ²)·g₀ − Σ_{j ∈ levels} δ²·2^{−2j}·g_j +where `g_var` is a 1D `DenseAxisArray` over the levels axis for a single +(name, t). Shared by sawtooth (PWL approximation) and epigraph (tangent cuts). +""" +@inline function scale_back_g_basis(x_min, delta, g_var, levels) + return x_min^2 + + (2.0 * x_min * delta + delta^2) * g_var[0] - + sum(delta^2 * 2.0^(-2j) * g_var[j] for j in levels) +end + +""" +Config for epigraph (Q^{L1}) LP-only lower-bound quadratic approximation. + +# Fields +- `depth::Int`: number of tangent-line breakpoints (2^depth + 1 tangent lines); + pure LP, zero binary variables. +""" +struct EpigraphQuadConfig <: QuadraticApproxConfig + depth::Int +end + +""" + build_quadratic_approx(config::EpigraphQuadConfig, model, x, x_min, x_max) + +Scalar form: build the epigraph relaxation for a single JuMP scalar `x` +with bounds `[x_min, x_max]`. Creates the per-cell g basis (j ∈ 0:depth), +LP relaxation constraints, z variable, tangent expression `fL`, and the +`depth + 2` tangent cuts. + +Returns a NamedTuple with `(approximation, z_var, g_var, link_constraint, +lp_constraints, tangent_expression, tangent_constraints)`. +""" +function build_quadratic_approx( + config::EpigraphQuadConfig, + model::JuMP.Model, + x::JuMP.AbstractJuMPScalar, + x_min::Float64, + x_max::Float64, +) + @assert config.depth >= 1 + @assert x_max > x_min + + depth = config.depth + delta = x_max - x_min + z_ub = max(x_min^2, x_max^2) + + g_var = JuMP.@variable( + model, [j = 0:depth], + lower_bound = 0.0, upper_bound = 1.0, + base_name = "SawtoothAux", + ) + + link_con = JuMP.@constraint(model, g_var[0] == (x - x_min) / delta) + + lp_a = JuMP.@constraint(model, [j = 1:depth], g_var[j] <= 2.0 * g_var[j - 1]) + lp_b = JuMP.@constraint( + model, [j = 1:depth], g_var[j] <= 2.0 * (1.0 - g_var[j - 1]), + ) + lp_cons = JuMP.Containers.DenseAxisArray{eltype(lp_a)}( + undef, 1:depth, 1:2, + ) + @views lp_cons.data[:, 1] .= lp_a + @views lp_cons.data[:, 2] .= lp_b + + z_var = JuMP.@variable( + model, lower_bound = 0.0, upper_bound = z_ub, base_name = "EpigraphVar", + ) + + fL_expr = JuMP.@expression( + model, + sum(delta^2 * 2.0^(-2j) * g_var[j] for j in 1:depth), + ) + + tangent_zero = JuMP.@constraint(model, z_var >= 0.0) + tangent_anchor = JuMP.@constraint( + model, z_var >= 2.0 * x_min - 1.0 + 2.0 * delta * g_var[0], + ) + tangent_levels = JuMP.@constraint( + model, [j = 1:depth], + z_var >= + scale_back_g_basis(x_min, delta, g_var, 1:j) - + delta^2 * 2.0^(-2j - 2), + ) + + tangent_cons = JuMP.Containers.DenseAxisArray{typeof(tangent_zero)}( + undef, 1:(depth + 2), + ) + tangent_cons[1] = tangent_zero + tangent_cons[2] = tangent_anchor + @views tangent_cons.data[3:end] .= tangent_levels + + approximation = JuMP.@expression(model, 1.0 * z_var) + + return (; + approximation, + z_var, + g_var, + link_constraint = link_con, + lp_constraints = lp_cons, + tangent_expression = fL_expr, + tangent_constraints = tangent_cons, + ) +end + +# --- IOM allocation + per-cell write helpers (shared with sawtooth tightening +# and NMDT tightening) --- + +function _alloc_epigraph_targets!( + container::OptimizationContainer, ::Type{C}, name_axis, time_axis, depth::Int, meta, +) where {C <: IS.InfrastructureSystemsComponent} + return ( + z = add_variable_container!( + container, EpigraphVariable, C, name_axis, time_axis; meta, + ), + g = add_variable_container!( + container, SawtoothAuxVariable, C, name_axis, 0:depth, time_axis; meta, + ), + link = add_constraints_container!( + container, SawtoothLinkingConstraint, C, name_axis, time_axis; meta, + ), + fL = add_expression_container!( + container, EpigraphTangentExpression, C, name_axis, time_axis; meta, + ), + approx = add_expression_container!( + container, EpigraphExpression, C, name_axis, time_axis; meta, + ), + lp = add_constraints_container!( + container, SawtoothLPConstraint, C, name_axis, 1:depth, 1:2, time_axis; + meta, + ), + tangent = add_constraints_container!( + container, EpigraphTangentConstraint, C, name_axis, 1:(depth + 2), time_axis; + meta, + ), + ) +end + +function _write_epigraph_cell!(targets, name, t, r, depth::Int) + targets.z[name, t] = r.z_var + for j in 0:depth + targets.g[name, j, t] = r.g_var[j] + end + targets.link[name, t] = r.link_constraint + targets.fL[name, t] = r.tangent_expression + targets.approx[name, t] = r.approximation + for j in 1:depth, k in 1:2 + targets.lp[name, j, k, t] = r.lp_constraints[j, k] + end + for j in 1:(depth + 2) + targets.tangent[name, j, t] = r.tangent_constraints[j] + end + return +end + +""" + add_quadratic_approx!(config::EpigraphQuadConfig, container, ::Type{C}, x_var, x_bounds, meta) + +Allocate all output containers (z, g, link/lp/tangent constraints, fL and +approximation expressions) with axes drawn from `x_var`'s `(name, t)` plus +the internal `(depth)` axes, then loop `(name, t)` calling the scalar +build per cell. +""" +function add_quadratic_approx!( + config::EpigraphQuadConfig, + container::OptimizationContainer, + ::Type{C}, + x_var, + x_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + name_axis = axes(x_var, 1) + time_axis = axes(x_var, 2) + depth = config.depth + @assert depth >= 1 + @assert length(name_axis) == length(x_bounds) + for b in x_bounds + @assert b.max > b.min + end + + model = get_jump_model(container) + targets = _alloc_epigraph_targets!(container, C, name_axis, time_axis, depth, meta) + + for (i, name) in enumerate(name_axis) + xmn, xmx = x_bounds[i].min, x_bounds[i].max + for t in time_axis + r = build_quadratic_approx(config, model, x_var[name, t], xmn, xmx) + _write_epigraph_cell!(targets, name, t, r, depth) + end + end + return targets.approx +end diff --git a/src/approximations/hybs.jl b/src/approximations/hybs.jl new file mode 100644 index 00000000..8b02c63e --- /dev/null +++ b/src/approximations/hybs.jl @@ -0,0 +1,270 @@ +# HybS (Hybrid Separable) MIP relaxation for bilinear products z = x·y. +# Combines a Bin2 lower bound and a Bin3 upper bound with shared quadratic +# approximations of x², y² and pure-LP epigraph relaxations of (x+y)², (x−y)². +# Uses 2L binaries instead of Bin2's 3L. +# Reference: Beach, Burlacu, Bärmann, Hager, Hildebrand (2024), Definition 10. + +"Two-sided HybS bound constraints: Bin2 lower + Bin3 upper." +struct HybSBoundConstraint <: ConstraintType end + +""" +Config for HybS bilinear approximation. + +# Fields +- `quad_config::QuadraticApproxConfig`: quadratic method used for x² and y². +- `epigraph_depth::Int`: depth of the epigraph Q^{L1} approximation of (x±y)². +- `add_mccormick::Bool`: whether to add a standard McCormick envelope on z (default false). +""" +struct HybSConfig{QC <: QuadraticApproxConfig} <: BilinearApproxConfig + quad_config::QC + epigraph_depth::Int + add_mccormick::Bool +end +function HybSConfig(quad_config::QuadraticApproxConfig, epigraph_depth::Int) + return HybSConfig(quad_config, epigraph_depth, false) +end + +""" + build_bilinear_approx(config::HybSConfig, model, x, y, x_min, x_max, y_min, y_max) + +Scalar form: build x² and y² via the chosen quadratic method, build +(x+y)² and (x−y)² via the epigraph Q^{L1} relaxation, and constrain a +fresh product variable z with two-sided bounds derived from the Bin2 lower +/ Bin3 upper identities. + +Returns `(; approximation, xsq, ysq, sum_expression, diff_expression, +sum_epigraph, diff_epigraph, z_var, bound_constraints, mccormick_constraints)`. +""" +function build_bilinear_approx( + config::HybSConfig, + model::JuMP.Model, + x::JuMP.AbstractJuMPScalar, + y::JuMP.AbstractJuMPScalar, + x_min::Float64, + x_max::Float64, + y_min::Float64, + y_max::Float64, +) + xsq = build_quadratic_approx(config.quad_config, model, x, x_min, x_max) + ysq = build_quadratic_approx(config.quad_config, model, y, y_min, y_max) + return _build_hybs_scalar( + config, model, x, y, xsq, ysq, x_min, x_max, y_min, y_max, + ) +end + +# Shared math between the standard and precomputed-form scalar entrypoints. +# `xsq`/`ysq` are NamedTuples (or any object with an `.approximation` field). +function _build_hybs_scalar( + config::HybSConfig, + model::JuMP.Model, + x::JuMP.AbstractJuMPScalar, + y::JuMP.AbstractJuMPScalar, + xsq, + ysq, + x_min::Float64, + x_max::Float64, + y_min::Float64, + y_max::Float64, +) + p1_expr = JuMP.@expression(model, x + y) + p2_expr = JuMP.@expression(model, x - y) + p1_min, p1_max = x_min + y_min, x_max + y_max + p2_min, p2_max = x_min - y_max, x_max - y_min + + epi_cfg = EpigraphQuadConfig(config.epigraph_depth) + zp1 = build_quadratic_approx(epi_cfg, model, p1_expr, p1_min, p1_max) + zp2 = build_quadratic_approx(epi_cfg, model, p2_expr, p2_min, p2_max) + + 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) + z_var = JuMP.@variable( + model, lower_bound = z_lo, upper_bound = z_hi, base_name = "HybSProduct", + ) + + bound_1 = JuMP.@constraint( + model, + z_var >= 0.5 * (zp1.approximation - xsq.approximation - ysq.approximation), + ) + bound_2 = JuMP.@constraint( + model, + z_var <= 0.5 * (xsq.approximation + ysq.approximation - zp2.approximation), + ) + bound_cons = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, 1:2) + bound_cons[1] = bound_1 + bound_cons[2] = bound_2 + + approximation = JuMP.@expression(model, 1.0 * z_var) + + mc = if config.add_mccormick + build_mccormick_envelope( + model, x, y, z_var, x_min, x_max, y_min, y_max, + ) + else + nothing + end + + return (; + approximation, + xsq, + ysq, + sum_expression = p1_expr, + diff_expression = p2_expr, + sum_epigraph = zp1, + diff_epigraph = zp2, + z_var, + bound_constraints = bound_cons, + mccormick_constraints = mc, + ) +end + +""" + add_bilinear_approx!(config::HybSConfig, container, ::Type{C}, x_var, y_var, x_bounds, y_bounds, meta) + +Build x² and y² via `add_quadratic_approx!(config.quad_config, ...)`, +build the (x+y) and (x−y) expression containers, build their epigraphs via +`add_quadratic_approx!(EpigraphQuadConfig(...), ...)`, then allocate the +HybS product variable + bound constraints and assemble per cell. +""" +function add_bilinear_approx!( + config::HybSConfig, + container::OptimizationContainer, + ::Type{C}, + x_var, + y_var, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + xsq = add_quadratic_approx!( + config.quad_config, + container, + C, + x_var, + x_bounds, + meta * "_x", + ) + ysq = add_quadratic_approx!( + config.quad_config, + container, + C, + y_var, + y_bounds, + meta * "_y", + ) + return _add_hybs_adapter!( + container, C, config, x_var, y_var, xsq, ysq, x_bounds, y_bounds, meta, + ) +end + +""" + add_bilinear_approx!(config::HybSConfig, container, ::Type{C}, xsq, ysq, x_var, y_var, x_bounds, y_bounds, meta) + +Precomputed-form: accepts already-built `xsq` ≈ x² and `ysq` ≈ y² 2D +expression containers; builds the HybS pieces on top. +""" +function add_bilinear_approx!( + config::HybSConfig, + container::OptimizationContainer, + ::Type{C}, + xsq, + ysq, + x_var, + y_var, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + return _add_hybs_adapter!( + container, C, config, x_var, y_var, xsq, ysq, x_bounds, y_bounds, meta, + ) +end + +# Allocate the HybS-specific containers (sum/diff exprs, two epigraphs, z var, +# bounds, approximation, optional McCormick) and loop (name, t) to assemble. +function _add_hybs_adapter!( + container::OptimizationContainer, ::Type{C}, config::HybSConfig, + x_var, y_var, xsq, ysq, x_bounds, y_bounds, meta, +) where {C <: IS.InfrastructureSystemsComponent} + name_axis = axes(x_var, 1) + time_axis = axes(x_var, 2) + @assert length(name_axis) == length(x_bounds) + @assert length(name_axis) == length(y_bounds) + model = get_jump_model(container) + + p1_target = add_expression_container!( + container, VariableSumExpression, C, name_axis, time_axis; + meta = meta * "_plus", + ) + p2_target = add_expression_container!( + container, VariableDifferenceExpression, C, name_axis, time_axis; + meta = meta * "_diff", + ) + for (i, name) in enumerate(name_axis) + for t in time_axis + p1_target[name, t] = x_var[name, t] + y_var[name, t] + p2_target[name, t] = x_var[name, t] - y_var[name, t] + end + end + p1_bounds = [ + (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 = [ + (min = x_bounds[i].min - y_bounds[i].max, + max = x_bounds[i].max - y_bounds[i].min) + for i in eachindex(x_bounds) + ] + + epi_cfg = EpigraphQuadConfig(config.epigraph_depth) + zp1 = add_quadratic_approx!(epi_cfg, container, C, p1_target, p1_bounds, meta * "_plus") + zp2 = add_quadratic_approx!(epi_cfg, container, C, p2_target, p2_bounds, meta * "_diff") + + z_target = add_variable_container!( + container, BilinearProductVariable, C, name_axis, time_axis; meta, + ) + bound_target = add_constraints_container!( + container, HybSBoundConstraint, C, name_axis, 1:2, time_axis; meta, + ) + approx_target = add_expression_container!( + container, BilinearProductExpression, C, name_axis, time_axis; meta, + ) + mc_target = if config.add_mccormick + add_constraints_container!( + container, McCormickConstraint, C, name_axis, 1:4, time_axis; meta, + ) + else + nothing + end + + for (i, name) in enumerate(name_axis) + xmn, xmx = x_bounds[i].min, x_bounds[i].max + ymn, ymx = y_bounds[i].min, y_bounds[i].max + z_lo = min(xmn * ymn, xmn * ymx, xmx * ymn, xmx * ymx) + z_hi = max(xmn * ymn, xmn * ymx, xmx * ymn, xmx * ymx) + for t in time_axis + z_v = JuMP.@variable( + model, lower_bound = z_lo, upper_bound = z_hi, + base_name = "HybSProduct", + ) + z_target[name, t] = z_v + bound_target[name, 1, t] = JuMP.@constraint( + model, z_v >= 0.5 * (zp1[name, t] - xsq[name, t] - ysq[name, t]), + ) + bound_target[name, 2, t] = JuMP.@constraint( + model, z_v <= 0.5 * (xsq[name, t] + ysq[name, t] - zp2[name, t]), + ) + approx_target[name, t] = JuMP.@expression(model, 1.0 * z_v) + if config.add_mccormick + r = build_mccormick_envelope( + model, x_var[name, t], y_var[name, t], z_v, + xmn, xmx, ymn, ymx, + ) + for (k, ref) in enumerate(r) + mc_target[name, k, t] = ref + end + end + end + end + return approx_target +end diff --git a/src/quadratic_approximations/incremental.jl b/src/approximations/incremental.jl similarity index 100% rename from src/quadratic_approximations/incremental.jl rename to src/approximations/incremental.jl diff --git a/src/approximations/manual_sos2.jl b/src/approximations/manual_sos2.jl new file mode 100644 index 00000000..3ec589a5 --- /dev/null +++ b/src/approximations/manual_sos2.jl @@ -0,0 +1,212 @@ +# SOS2 piecewise linear approximation of x² with manually-implemented adjacency +# via binary segment-selector variables. Useful when the solver does not +# natively support MOI.SOS2. + +"Binary segment-selection variables (z) for manual SOS2 quadratic approximation." +struct ManualSOS2BinaryVariable <: SparseVariableType end + +"Ensures exactly one segment is active (∑ z_j = 1) in manual SOS2 quadratic approximation." +struct ManualSOS2SegmentSelectionConstraint <: ConstraintType end + +"Expression for the segment selection sum Σ z_j in manual SOS2 quadratic approximation." +struct ManualSOS2SegmentSelectionExpression <: ExpressionType end + +"Links active segment to lambda variables." +struct ManualSOS2AdjacencyConstraint <: ConstraintType end + +""" +Config for manual binary-variable SOS2 quadratic approximation. + +# Fields +- `depth::Int`: number of PWL segments (breakpoints = depth + 1). +- `pwmcc_segments::Int`: number of piecewise McCormick cut partitions; + 0 to disable (default 4). +""" +struct ManualSOS2QuadConfig <: QuadraticApproxConfig + depth::Int + pwmcc_segments::Int +end +function ManualSOS2QuadConfig(depth::Int) + return ManualSOS2QuadConfig(depth, 4) +end + +""" + build_quadratic_approx(config::ManualSOS2QuadConfig, model, x, x_min, x_max) + +Scalar form: PWL approximation of x² with manually-enforced SOS2 adjacency +via binary segment-selectors z_j and constraints λ_i ≤ z_{i-1} + z_i (with +boundary cases at i=1 and i=n_points). If `config.pwmcc_segments > 0`, +also adds piecewise McCormick concave cuts per cell. + +Returns a NamedTuple with `(approximation, lambda, z_var, link_constraint, +norm_constraint, segment_sum_constraint, adjacency_constraints, +link_expression, norm_expression, segment_sum_expression, pwmcc)`. +""" +function build_quadratic_approx( + config::ManualSOS2QuadConfig, + model::JuMP.Model, + x::JuMP.AbstractJuMPScalar, + x_min::Float64, + x_max::Float64, +) + @assert x_max > x_min + n_points = config.depth + 1 + n_bins = n_points - 1 + x_bkpts, x_sq_bkpts = _get_breakpoints_for_pwl_function( + 0.0, 1.0, _square; num_segments = config.depth, + ) + lx = x_max - x_min + + lambda = JuMP.@variable( + model, [i = 1:n_points], + lower_bound = 0.0, upper_bound = 1.0, + base_name = "QuadraticVariable", + ) + z_var = JuMP.@variable( + model, [j = 1:n_bins], binary = true, base_name = "ManualSOS2Binary", + ) + + link_expr = JuMP.@expression( + model, sum(x_bkpts[i] * lambda[i] for i in 1:n_points), + ) + link_con = JuMP.@constraint(model, (x - x_min) / lx == link_expr) + norm_expr = JuMP.@expression(model, sum(lambda[i] for i in 1:n_points)) + norm_con = JuMP.@constraint(model, norm_expr == 1.0) + seg_expr = JuMP.@expression(model, sum(z_var[j] for j in 1:n_bins)) + seg_con = JuMP.@constraint(model, seg_expr == 1) + + adj_first = JuMP.@constraint(model, lambda[1] <= z_var[1]) + adj_interior = JuMP.@constraint( + model, [i = 2:(n_points - 1)], lambda[i] <= z_var[i - 1] + z_var[i], + ) + adj_last = JuMP.@constraint(model, lambda[n_points] <= z_var[n_bins]) + + adj_cons = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, 1:n_points) + adj_cons[1] = adj_first + if n_points >= 3 + for i in 2:(n_points - 1) + adj_cons[i] = adj_interior[i] + end + end + adj_cons[n_points] = adj_last + + approximation = JuMP.@expression( + model, + lx * lx * sum(x_sq_bkpts[i] * lambda[i] for i in 1:n_points) + + 2.0 * x_min * x - x_min * x_min, + ) + + pwmcc = if config.pwmcc_segments > 0 + build_pwmcc_concave_cuts( + model, x, approximation, x_min, x_max, config.pwmcc_segments, + ) + else + nothing + end + + return (; + approximation, + lambda, + z_var, + link_constraint = link_con, + norm_constraint = norm_con, + segment_sum_constraint = seg_con, + adjacency_constraints = adj_cons, + link_expression = link_expr, + norm_expression = norm_expr, + segment_sum_expression = seg_expr, + pwmcc, + ) +end + +""" + add_quadratic_approx!(config::ManualSOS2QuadConfig, container, ::Type{C}, x_var, x_bounds, meta) + +Allocate manual-SOS2 containers (λ, z, link/norm/seg/adjacency, expressions, +approximation) plus, when `config.pwmcc_segments > 0`, the PWMCC containers +under `meta * "_pwmcc"`. Loop `(name, t)`. +""" +function add_quadratic_approx!( + config::ManualSOS2QuadConfig, + container::OptimizationContainer, + ::Type{C}, + x_var, + x_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + name_axis = axes(x_var, 1) + time_axis = axes(x_var, 2) + n_points = config.depth + 1 + n_bins = n_points - 1 + @assert length(name_axis) == length(x_bounds) + for b in x_bounds + @assert b.max > b.min + end + + model = get_jump_model(container) + + lambda_target = add_variable_container!( + container, QuadraticVariable, C, name_axis, 1:n_points, time_axis; meta, + ) + z_target = add_variable_container!( + container, ManualSOS2BinaryVariable, C, name_axis, 1:n_bins, time_axis; meta, + ) + link_cons_target = add_constraints_container!( + container, SOS2LinkingConstraint, C, name_axis, time_axis; meta, + ) + norm_cons_target = add_constraints_container!( + container, SOS2NormConstraint, C, name_axis, time_axis; meta, + ) + seg_cons_target = add_constraints_container!( + container, ManualSOS2SegmentSelectionConstraint, C, name_axis, time_axis; meta, + ) + link_expr_target = add_expression_container!( + container, SOS2LinkingExpression, C, name_axis, time_axis; meta, + ) + norm_expr_target = add_expression_container!( + container, SOS2NormExpression, C, name_axis, time_axis; meta, + ) + seg_expr_target = add_expression_container!( + container, ManualSOS2SegmentSelectionExpression, C, name_axis, time_axis; meta, + ) + approx_target = add_expression_container!( + container, QuadraticExpression, C, name_axis, time_axis; meta, + ) + adj_target = add_constraints_container!( + container, ManualSOS2AdjacencyConstraint, C, name_axis, 1:n_points, time_axis; + meta, + ) + + use_pwmcc = config.pwmcc_segments > 0 + K = config.pwmcc_segments + pwmcc_targets = if use_pwmcc + _alloc_pwmcc_targets!(container, C, name_axis, time_axis, K, meta * "_pwmcc") + else + nothing + end + + for (i, name) in enumerate(name_axis) + xmn, xmx = x_bounds[i].min, x_bounds[i].max + for t in time_axis + r = build_quadratic_approx(config, model, x_var[name, t], xmn, xmx) + for ip in 1:n_points + lambda_target[name, ip, t] = r.lambda[ip] + adj_target[name, ip, t] = r.adjacency_constraints[ip] + end + for j in 1:n_bins + z_target[name, j, t] = r.z_var[j] + end + link_cons_target[name, t] = r.link_constraint + norm_cons_target[name, t] = r.norm_constraint + seg_cons_target[name, t] = r.segment_sum_constraint + link_expr_target[name, t] = r.link_expression + norm_expr_target[name, t] = r.norm_expression + seg_expr_target[name, t] = r.segment_sum_expression + approx_target[name, t] = r.approximation + if use_pwmcc + _write_pwmcc_cell!(pwmcc_targets, name, t, r.pwmcc, K) + end + end + end + return approx_target +end diff --git a/src/approximations/mccormick.jl b/src/approximations/mccormick.jl new file mode 100644 index 00000000..db0e940a --- /dev/null +++ b/src/approximations/mccormick.jl @@ -0,0 +1,256 @@ +# McCormick envelope for bilinear products z = x·y. +# +# Three scalar build_* functions form the pure-JuMP math layer (no IOM deps): +# build_mccormick_envelope — 4 constraints (upper_1, upper_2, lower_1, lower_2) +# build_mccormick_upper — 2 constraints (upper_1, upper_2 only) +# build_reformulated_mccormick — 4 constraints for the Bin2 separable identity +# +# Each has a matching IOM adapter (add_*_approx!) that allocates a container +# with a 1:N dimension matching the math's return shape, loops over (name, t), +# and writes the scalar refs into the slots. + +# --- Container key types --- + +"McCormick envelope upper-bound constraints: z ≤ ... (axis 1:2 in container)." +struct McCormickUpperConstraint <: ConstraintType end + +"Combined McCormick envelope constraints for the full 4-constraint envelope (axis 1:4)." +struct McCormickConstraint <: ConstraintType end + +"Reformulated McCormick constraints on Bin2 separable variables (axis 1:4)." +struct ReformulatedMcCormickConstraint <: ConstraintType end + +# --- Scalar build_* (pure JuMP) --- + +""" + build_mccormick_envelope(model, x, y, z, x_min, x_max, y_min, y_max) + +Build the four McCormick inequalities for z ≈ x·y at a single cell. +Returns a flat NamedTuple `(upper_1, upper_2, lower_1, lower_2)` of scalar +constraint refs. +""" +function build_mccormick_envelope( + model::JuMP.Model, + x::JuMP.AbstractJuMPScalar, + y::JuMP.AbstractJuMPScalar, + z::JuMP.AbstractJuMPScalar, + x_min::Float64, + x_max::Float64, + y_min::Float64, + y_max::Float64, +) + upper_1 = JuMP.@constraint(model, z <= x_max * y + x * y_min - x_max * y_min) + upper_2 = JuMP.@constraint(model, z <= x_min * y + x * y_max - x_min * y_max) + lower_1 = JuMP.@constraint(model, z >= x_min * y + x * y_min - x_min * y_min) + lower_2 = JuMP.@constraint(model, z >= x_max * y + x * y_max - x_max * y_max) + return (; upper_1, upper_2, lower_1, lower_2) +end + +""" + build_mccormick_upper(model, x, y, z, x_min, x_max, y_min, y_max) + +Build only the upper-envelope McCormick inequalities (z ≤ ...) at a single +cell. Used when a tighter lower bound on z is supplied elsewhere. + +Returns a flat NamedTuple `(upper_1, upper_2)` of scalar constraint refs. +""" +function build_mccormick_upper( + model::JuMP.Model, + x::JuMP.AbstractJuMPScalar, + y::JuMP.AbstractJuMPScalar, + z::JuMP.AbstractJuMPScalar, + x_min::Float64, + x_max::Float64, + y_min::Float64, + y_max::Float64, +) + upper_1 = JuMP.@constraint(model, z <= x_max * y + x * y_min - x_max * y_min) + upper_2 = JuMP.@constraint(model, z <= x_min * y + x * y_max - x_min * y_max) + return (; upper_1, upper_2) +end + +""" + build_reformulated_mccormick(model, x, y, zp1, zx, zy, x_min, x_max, y_min, y_max) + +Build the four reformulated-McCormick inequalities for the Bin2 separable +identity (zp1 ≈ (x+y)², zx ≈ x², zy ≈ y²) at a single cell. + +Returns a flat NamedTuple `(c1, c2, c3, c4)` of scalar constraint refs. +""" +function build_reformulated_mccormick( + model::JuMP.Model, + x::JuMP.AbstractJuMPScalar, + y::JuMP.AbstractJuMPScalar, + zp1::JuMP.AbstractJuMPScalar, + zx::JuMP.AbstractJuMPScalar, + zy::JuMP.AbstractJuMPScalar, + x_min::Float64, + x_max::Float64, + y_min::Float64, + y_max::Float64, +) + c1 = JuMP.@constraint( + model, + zp1 - zx - zy >= 2.0 * (x_min * y + x * y_min - x_min * y_min), + ) + c2 = JuMP.@constraint( + model, + zp1 - zx - zy >= 2.0 * (x_max * y + x * y_max - x_max * y_max), + ) + c3 = JuMP.@constraint( + model, + zp1 - zx - zy <= 2.0 * (x_max * y + x * y_min - x_max * y_min), + ) + c4 = JuMP.@constraint( + model, + zp1 - zx - zy <= 2.0 * (x_min * y + x * y_max - x_min * y_max), + ) + return (; c1, c2, c3, c4) +end + +# --- IOM adapters --- + +""" + add_mccormick_approx!(container, ::Type{C}, x_var, y_var, z_var, x_bounds, y_bounds, meta) + +Allocate a `McCormickConstraint` container with axes `(name, 1:4, time)`, +loop `(name, t)`, call `build_mccormick_envelope` per cell, and slot the +four refs. +""" +function add_mccormick_approx!( + container::OptimizationContainer, + ::Type{C}, + x_var, + y_var, + z_var, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + name_axis = axes(x_var, 1) + time_axis = axes(x_var, 2) + @assert length(name_axis) == length(x_bounds) + @assert length(name_axis) == length(y_bounds) + for i in eachindex(x_bounds) + @assert x_bounds[i].max > x_bounds[i].min + @assert y_bounds[i].max > y_bounds[i].min + end + + model = get_jump_model(container) + target = add_constraints_container!( + container, McCormickConstraint, C, name_axis, 1:4, time_axis; meta, + ) + + for (i, name) in enumerate(name_axis) + xmn, xmx = x_bounds[i].min, x_bounds[i].max + ymn, ymx = y_bounds[i].min, y_bounds[i].max + for t in time_axis + r = build_mccormick_envelope( + model, + x_var[name, t], y_var[name, t], z_var[name, t], + xmn, xmx, ymn, ymx, + ) + for (k, ref) in enumerate(r) + target[name, k, t] = ref + end + end + end + return target +end + +""" + add_mccormick_upper_approx!(container, ::Type{C}, x_var, y_var, z_var, x_bounds, y_bounds, meta) + +Allocate `McCormickUpperConstraint` `(name, 1:2, time)`; loop, call +`build_mccormick_upper` per cell, slot two refs. +""" +function add_mccormick_upper_approx!( + container::OptimizationContainer, + ::Type{C}, + x_var, + y_var, + z_var, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + name_axis = axes(x_var, 1) + time_axis = axes(x_var, 2) + @assert length(name_axis) == length(x_bounds) + @assert length(name_axis) == length(y_bounds) + for i in eachindex(x_bounds) + @assert x_bounds[i].max > x_bounds[i].min + @assert y_bounds[i].max > y_bounds[i].min + end + + model = get_jump_model(container) + target = add_constraints_container!( + container, McCormickUpperConstraint, C, name_axis, 1:2, time_axis; meta, + ) + + for (i, name) in enumerate(name_axis) + xmn, xmx = x_bounds[i].min, x_bounds[i].max + ymn, ymx = y_bounds[i].min, y_bounds[i].max + for t in time_axis + r = build_mccormick_upper( + model, + x_var[name, t], y_var[name, t], z_var[name, t], + xmn, xmx, ymn, ymx, + ) + for (k, ref) in enumerate(r) + target[name, k, t] = ref + end + end + end + return target +end + +""" + add_reformulated_mccormick_approx!(container, ::Type{C}, x_var, y_var, zp1, zx, zy, x_bounds, y_bounds, meta) + +Allocate `ReformulatedMcCormickConstraint` `(name, 1:4, time)`; loop, call +`build_reformulated_mccormick` per cell, slot four refs. +""" +function add_reformulated_mccormick_approx!( + container::OptimizationContainer, + ::Type{C}, + x_var, + y_var, + zp1_expr, + zx_expr, + zy_expr, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + name_axis = axes(x_var, 1) + time_axis = axes(x_var, 2) + @assert length(name_axis) == length(x_bounds) + @assert length(name_axis) == length(y_bounds) + for i in eachindex(x_bounds) + @assert x_bounds[i].max > x_bounds[i].min + @assert y_bounds[i].max > y_bounds[i].min + end + + model = get_jump_model(container) + target = add_constraints_container!( + container, ReformulatedMcCormickConstraint, C, name_axis, 1:4, time_axis; meta, + ) + + for (i, name) in enumerate(name_axis) + xmn, xmx = x_bounds[i].min, x_bounds[i].max + ymn, ymx = y_bounds[i].min, y_bounds[i].max + for t in time_axis + r = build_reformulated_mccormick( + model, + x_var[name, t], y_var[name, t], + zp1_expr[name, t], zx_expr[name, t], zy_expr[name, t], + xmn, xmx, ymn, ymx, + ) + for (k, ref) in enumerate(r) + target[name, k, t] = ref + end + end + end + return target +end diff --git a/src/approximations/nmdt_bilinear.jl b/src/approximations/nmdt_bilinear.jl new file mode 100644 index 00000000..b11dacb9 --- /dev/null +++ b/src/approximations/nmdt_bilinear.jl @@ -0,0 +1,292 @@ +# NMDT bilinear approximations of z = x·y. +# NMDTBilinearConfig — discretize x only; y is just normalized. +# DNMDTBilinearConfig — discretize both x and y, combine two estimates. + +""" +Config for single-NMDT bilinear approximation (discretizes x only). + +# Fields +- `depth::Int`: number of binary discretization levels L for x. +""" +struct NMDTBilinearConfig <: BilinearApproxConfig + depth::Int +end + +""" +Config for double-NMDT bilinear approximation (discretizes both x and y). + +# Fields +- `depth::Int`: number of binary discretization levels L for both x and y. +""" +struct DNMDTBilinearConfig <: BilinearApproxConfig + depth::Int +end + +# --- Scalar build (pure JuMP) --- + +""" + build_bilinear_approx(config::NMDTBilinearConfig, model, x, y, x_min, x_max, y_min, y_max) + +Scalar form: approximate x·y via NMDT for one cell. Discretize x, normalize +y to yh ∈ [0,1], build the binary-continuous product β·yh and residual +δ·yh, reassemble x·y. + +Returns `(; approximation, x_discretization, yh_expression, bx_yh_product, +residual_product)`. +""" +function build_bilinear_approx( + config::NMDTBilinearConfig, + model::JuMP.Model, + x::JuMP.AbstractJuMPScalar, + y::JuMP.AbstractJuMPScalar, + x_min::Float64, + x_max::Float64, + y_min::Float64, + y_max::Float64, +) + x_disc = build_discretization(model, x, x_min, x_max, config.depth) + yh_expr = (y - y_min) / (y_max - y_min) + bx_yh = build_binary_continuous_product( + model, x_disc.beta_var, yh_expr, 0.0, 1.0, config.depth, + ) + dz = build_residual_product( + model, x_disc.delta_var, yh_expr, 1.0, config.depth, + ) + approximation = build_assembled_product( + model, [bx_yh.result_expression], dz.z_var, + x_disc.norm_expr, yh_expr, + x_min, x_max, y_min, y_max, + ) + return (; + approximation, + x_discretization = x_disc, + yh_expression = yh_expr, + bx_yh_product = bx_yh, + residual_product = dz, + ) +end + +""" + build_bilinear_approx(config::DNMDTBilinearConfig, model, x, y, x_min, x_max, y_min, y_max) + +Scalar form: DNMDT bilinear approximation at one cell. Discretize both x +and y, form all four cross binary-continuous products, convex-combine. +""" +function build_bilinear_approx( + config::DNMDTBilinearConfig, + model::JuMP.Model, + x::JuMP.AbstractJuMPScalar, + y::JuMP.AbstractJuMPScalar, + x_min::Float64, + x_max::Float64, + y_min::Float64, + y_max::Float64, +) + x_disc = build_discretization(model, x, x_min, x_max, config.depth) + y_disc = build_discretization(model, y, y_min, y_max, config.depth) + bx_yh = build_binary_continuous_product( + model, x_disc.beta_var, y_disc.norm_expr, 0.0, 1.0, config.depth, + ) + by_dx = build_binary_continuous_product( + model, y_disc.beta_var, x_disc.delta_var, + 0.0, 2.0^(-config.depth), config.depth, + ) + by_xh = build_binary_continuous_product( + model, y_disc.beta_var, x_disc.norm_expr, 0.0, 1.0, config.depth, + ) + bx_dy = build_binary_continuous_product( + model, x_disc.beta_var, y_disc.delta_var, + 0.0, 2.0^(-config.depth), config.depth, + ) + dz = build_residual_product( + model, x_disc.delta_var, y_disc.delta_var, + 2.0^(-config.depth), config.depth, + ) + approximation = build_assembled_dnmdt( + model, + bx_yh.result_expression, + by_dx.result_expression, + by_xh.result_expression, + bx_dy.result_expression, + dz.z_var, + x_disc.norm_expr, y_disc.norm_expr, + x_min, x_max, y_min, y_max, + ) + return (; + approximation, + x_discretization = x_disc, + y_discretization = y_disc, + bx_yh_product = bx_yh, + by_dx_product = by_dx, + by_xh_product = by_xh, + bx_dy_product = bx_dy, + residual_product = dz, + ) +end + +# --- IOM adapter --- + +""" + add_bilinear_approx!(config::NMDTBilinearConfig, container, ::Type{C}, x_var, y_var, x_bounds, y_bounds, meta) + +Allocate x discretization + yh expression + binary-continuous-product + +residual-product + BilinearProductExpression containers. Loop `(name, t)`. +""" +function add_bilinear_approx!( + config::NMDTBilinearConfig, + container::OptimizationContainer, + ::Type{C}, + x_var, + y_var, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + name_axis = axes(x_var, 1) + time_axis = axes(x_var, 2) + depth = config.depth + @assert length(name_axis) == length(x_bounds) + @assert length(name_axis) == length(y_bounds) + for i in eachindex(x_bounds) + @assert x_bounds[i].max > x_bounds[i].min + @assert y_bounds[i].max > y_bounds[i].min + end + + model = get_jump_model(container) + + x_disc_targets = _alloc_discretization_targets!( + container, C, name_axis, time_axis, depth, meta * "_x", + ) + yh_target = add_expression_container!( + container, NormedVariableExpression, C, name_axis, time_axis; + meta = meta * "_y", + ) + bx_yh_targets = _alloc_binary_continuous_product_targets!( + container, C, name_axis, time_axis, depth, meta; tighten = false, + ) + res_targets = _alloc_residual_product_targets!( + container, C, name_axis, time_axis, meta; tighten = false, + ) + approx_target = add_expression_container!( + container, BilinearProductExpression, C, name_axis, time_axis; meta, + ) + + for (i, name) in enumerate(name_axis) + xmn, xmx = x_bounds[i].min, x_bounds[i].max + ymn, ymx = y_bounds[i].min, y_bounds[i].max + for t in time_axis + r = build_bilinear_approx( + config, model, x_var[name, t], y_var[name, t], xmn, xmx, ymn, ymx, + ) + _write_discretization_cell!(x_disc_targets, name, t, r.x_discretization, depth) + yh_target[name, t] = r.yh_expression + _write_binary_continuous_product_cell!( + bx_yh_targets, + name, + t, + r.bx_yh_product, + depth, + ) + _write_residual_product_cell!(res_targets, name, t, r.residual_product) + approx_target[name, t] = r.approximation + end + end + return approx_target +end + +""" + add_bilinear_approx!(config::DNMDTBilinearConfig, container, ::Type{C}, x_var, y_var, x_bounds, y_bounds, meta) + +Allocate two discretizations + four binary-continuous-product + one +residual-product + BilinearProductExpression containers. Loop `(name, t)`. +""" +function add_bilinear_approx!( + config::DNMDTBilinearConfig, + container::OptimizationContainer, + ::Type{C}, + x_var, + y_var, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + name_axis = axes(x_var, 1) + time_axis = axes(x_var, 2) + depth = config.depth + @assert length(name_axis) == length(x_bounds) + @assert length(name_axis) == length(y_bounds) + for i in eachindex(x_bounds) + @assert x_bounds[i].max > x_bounds[i].min + @assert y_bounds[i].max > y_bounds[i].min + end + + model = get_jump_model(container) + + x_disc_targets = _alloc_discretization_targets!( + container, C, name_axis, time_axis, depth, meta * "_x", + ) + y_disc_targets = _alloc_discretization_targets!( + container, C, name_axis, time_axis, depth, meta * "_y", + ) + bx_yh_targets = _alloc_binary_continuous_product_targets!( + container, C, name_axis, time_axis, depth, meta * "_bx_yh"; tighten = false, + ) + by_dx_targets = _alloc_binary_continuous_product_targets!( + container, C, name_axis, time_axis, depth, meta * "_by_dx"; tighten = false, + ) + by_xh_targets = _alloc_binary_continuous_product_targets!( + container, C, name_axis, time_axis, depth, meta * "_by_xh"; tighten = false, + ) + bx_dy_targets = _alloc_binary_continuous_product_targets!( + container, C, name_axis, time_axis, depth, meta * "_bx_dy"; tighten = false, + ) + res_targets = _alloc_residual_product_targets!( + container, C, name_axis, time_axis, meta; tighten = false, + ) + approx_target = add_expression_container!( + container, BilinearProductExpression, C, name_axis, time_axis; meta, + ) + + for (i, name) in enumerate(name_axis) + xmn, xmx = x_bounds[i].min, x_bounds[i].max + ymn, ymx = y_bounds[i].min, y_bounds[i].max + for t in time_axis + r = build_bilinear_approx( + config, model, x_var[name, t], y_var[name, t], xmn, xmx, ymn, ymx, + ) + _write_discretization_cell!(x_disc_targets, name, t, r.x_discretization, depth) + _write_discretization_cell!(y_disc_targets, name, t, r.y_discretization, depth) + _write_binary_continuous_product_cell!( + bx_yh_targets, + name, + t, + r.bx_yh_product, + depth, + ) + _write_binary_continuous_product_cell!( + by_dx_targets, + name, + t, + r.by_dx_product, + depth, + ) + _write_binary_continuous_product_cell!( + by_xh_targets, + name, + t, + r.by_xh_product, + depth, + ) + _write_binary_continuous_product_cell!( + bx_dy_targets, + name, + t, + r.bx_dy_product, + depth, + ) + _write_residual_product_cell!(res_targets, name, t, r.residual_product) + approx_target[name, t] = r.approximation + end + end + return approx_target +end diff --git a/src/approximations/nmdt_discretization.jl b/src/approximations/nmdt_discretization.jl new file mode 100644 index 00000000..977ca2ec --- /dev/null +++ b/src/approximations/nmdt_discretization.jl @@ -0,0 +1,340 @@ +# Shared NMDT machinery used by both nmdt_quadratic and nmdt_bilinear. +# +# NMDT (Normalized Multiparametric Disaggregation Technique) discretizes a +# normalized variable xh ∈ [0,1] as xh = Σᵢ 2^{−i}·β_i + δ, with β_i ∈ {0,1} +# binary digits and δ ∈ [0, 2^{−L}] a residual. The discretization is then +# combined with another normalized variable (for bilinear products) or with +# itself (for quadratic) via McCormick-linearized binary-continuous products. + +# --- Container key types --- + +"Binary discretization variables β_i ∈ {0,1} in the NMDT decomposition of xh." +struct NMDTBinaryVariable <: VariableType end +"Residual variable δ ∈ [0, 2^{−L}] capturing the NMDT discretization error." +struct NMDTResidualVariable <: VariableType end +"McCormick linearization variables u_i ≈ β_i·y in NMDT binary-continuous products." +struct NMDTBinaryContinuousProductVariable <: VariableType end +"Variable z ≈ δ · y linearizing the residual-continuous product in NMDT." +struct NMDTResidualProductVariable <: VariableType end + +"Expression container for the NMDT binary discretization: Σ 2^{−i}·β_i + δ ≈ xh." +struct NMDTDiscretizationExpression <: ExpressionType end +"Expression container for the NMDT binary-continuous product: Σ 2^{−i}·u_i ≈ β·y." +struct NMDTBinaryContinuousProductExpression <: ExpressionType end + +"Constraint enforcing xh = Σ 2^{−i}·β_i + δ in the NMDT discretization." +struct NMDTEDiscretizationConstraint <: ConstraintType end +"Epigraph lower-bound tightening constraint on the NMDT quadratic result." +struct NMDTTightenConstraint <: ConstraintType end + +# --- Scalar build helpers (pure JuMP) --- + +""" + build_discretization(model, x, x_min, x_max, depth) + +Scalar form: build the NMDT binary discretization of the normalized +variable xh = (x − x_min)/(x_max − x_min) for a single cell. Creates +`depth` binary variables β_i and one residual δ ∈ [0, 2^{−depth}], +enforcing xh = Σ 2^{−i}·β_i + δ. + +Returns `(; norm_expr, beta_var, delta_var, disc_constraint, disc_expression)`. +""" +function build_discretization( + model::JuMP.Model, + x::JuMP.AbstractJuMPScalar, + x_min::Float64, + x_max::Float64, + depth::Int, +) + @assert depth >= 1 + @assert x_max > x_min + norm_expr = (x - x_min) / (x_max - x_min) + beta_var = JuMP.@variable( + model, [i = 1:depth], binary = true, base_name = "NMDTBinary", + ) + delta_var = JuMP.@variable( + model, + lower_bound = 0.0, upper_bound = 2.0^(-depth), + base_name = "NMDTResidual", + ) + disc_expr = JuMP.@expression( + model, sum(2.0^(-i) * beta_var[i] for i in 1:depth) + delta_var, + ) + disc_con = JuMP.@constraint(model, norm_expr == disc_expr) + return (; + norm_expr, + beta_var, + delta_var, + disc_constraint = disc_con, + disc_expression = disc_expr, + ) +end + +""" + build_binary_continuous_product(model, beta_var, cont_var, cont_min, cont_max, depth; tighten=false) + +Scalar form: build Σᵢ 2^{−i}·u_i ≈ β·y for a single cell. `beta_var` is a +1D `DenseAxisArray` over `1:depth` (binary), `cont_var` is a JuMP scalar. + +Returns `(; u_var, mccormick_lower, mccormick_upper, result_expression)` where +`mccormick_lower` is `nothing` when `tighten = true`. +""" +function build_binary_continuous_product( + model::JuMP.Model, + beta_var::AbstractVector, + cont_var::JuMP.AbstractJuMPScalar, + cont_min::Float64, + cont_max::Float64, + depth::Int; + tighten::Bool = false, +) + u_var = JuMP.@variable( + model, [i = 1:depth], + lower_bound = cont_min, upper_bound = cont_max, + base_name = "NMDTBinContProd", + ) + upper_1 = JuMP.@constraint( + model, [i = 1:depth], u_var[i] <= cont_max * beta_var[i], + ) + upper_2 = JuMP.@constraint( + model, [i = 1:depth], + u_var[i] <= cont_min * beta_var[i] + cont_var - cont_min, + ) + mccormick_lower = if tighten + nothing + else + lower_1 = JuMP.@constraint( + model, [i = 1:depth], u_var[i] >= cont_min * beta_var[i], + ) + lower_2 = JuMP.@constraint( + model, [i = 1:depth], + u_var[i] >= cont_max * beta_var[i] + cont_var - cont_max, + ) + (c1 = lower_1, c2 = lower_2) + end + result_expr = JuMP.@expression( + model, sum(2.0^(-i) * u_var[i] for i in 1:depth), + ) + return (; + u_var, + mccormick_lower, + mccormick_upper = (c1 = upper_1, c2 = upper_2), + result_expression = result_expr, + ) +end + +""" + build_residual_product(model, delta_var, cont_var, cont_max, depth; tighten=false) + +Scalar form: build z ≈ δ·y for a single cell, where δ ∈ [0, 2^{−depth}] +and y ∈ [0, cont_max], using McCormick envelopes. When `tighten`, only the +upper envelopes are added. + +Returns `(; z_var, mccormick_constraints)` where `mccormick_constraints` is +the NamedTuple returned by `build_mccormick_upper` (`tighten = true`) or +`build_mccormick_envelope` (`tighten = false`). +""" +function build_residual_product( + model::JuMP.Model, + delta_var::JuMP.AbstractJuMPScalar, + cont_var::JuMP.AbstractJuMPScalar, + cont_max::Float64, + depth::Int; + tighten::Bool = false, +) + delta_max = 2.0^(-depth) + z_var = JuMP.@variable( + model, + lower_bound = 0.0, upper_bound = delta_max * cont_max, + base_name = "NMDTResidualProduct", + ) + mc = if tighten + build_mccormick_upper( + model, delta_var, cont_var, z_var, 0.0, delta_max, 0.0, cont_max, + ) + else + build_mccormick_envelope( + model, delta_var, cont_var, z_var, 0.0, delta_max, 0.0, cont_max, + ) + end + return (; z_var, mccormick_constraints = mc) +end + +""" + build_assembled_product(model, terms, dz, xh_expr, yh_expr, x_min, x_max, y_min, y_max) + +Scalar form: affine reassembly of x·y from normalized NMDT pieces: + x·y = lx·ly·zh + lx·y_min·xh + ly·x_min·yh + x_min·y_min +where `zh = sum(terms) + dz`. `terms` is a list of scalar AffExpr values. +""" +function build_assembled_product( + model::JuMP.Model, + terms::AbstractVector, + dz::JuMP.AbstractJuMPScalar, + xh_expr, + yh_expr, + x_min::Float64, + x_max::Float64, + y_min::Float64, + y_max::Float64, +) + lx = x_max - x_min + ly = y_max - y_min + return JuMP.@expression( + model, + lx * ly * (sum(terms) + dz) + + lx * y_min * xh_expr + + ly * x_min * yh_expr + + x_min * y_min, + ) +end + +""" + build_assembled_dnmdt(model, bx_yh, by_dx, by_xh, bx_dy, dz, xh_expr, yh_expr, x_min, x_max, y_min, y_max; lambda) + +Scalar form: convex combination of two NMDT estimates of x·y at one cell. +""" +function build_assembled_dnmdt( + model::JuMP.Model, + bx_yh::JuMP.AbstractJuMPScalar, + by_dx::JuMP.AbstractJuMPScalar, + by_xh::JuMP.AbstractJuMPScalar, + bx_dy::JuMP.AbstractJuMPScalar, + dz::JuMP.AbstractJuMPScalar, + xh_expr, + yh_expr, + x_min::Float64, + x_max::Float64, + y_min::Float64, + y_max::Float64; + lambda::Float64 = DNMDT_LAMBDA, +) + z1 = build_assembled_product( + model, [bx_yh, by_dx], dz, xh_expr, yh_expr, + x_min, x_max, y_min, y_max, + ) + z2 = build_assembled_product( + model, [by_xh, bx_dy], dz, yh_expr, xh_expr, + y_min, y_max, x_min, x_max, + ) + return JuMP.@expression(model, lambda * z1 + (1.0 - lambda) * z2) +end + +# --- IOM allocation + per-cell write helpers (used by NMDT quad/bilinear adapters) --- + +function _alloc_discretization_targets!( + container::OptimizationContainer, ::Type{C}, name_axis, time_axis, depth::Int, meta, +) where {C <: IS.InfrastructureSystemsComponent} + return ( + norm = add_expression_container!( + container, NormedVariableExpression, C, name_axis, time_axis; meta, + ), + beta = add_variable_container!( + container, NMDTBinaryVariable, C, name_axis, 1:depth, time_axis; meta, + ), + delta = add_variable_container!( + container, NMDTResidualVariable, C, name_axis, time_axis; meta, + ), + disc_expr = add_expression_container!( + container, NMDTDiscretizationExpression, C, name_axis, time_axis; meta, + ), + disc_cons = add_constraints_container!( + container, NMDTEDiscretizationConstraint, C, name_axis, time_axis; meta, + ), + ) +end + +function _write_discretization_cell!(targets, name, t, r, depth::Int) + targets.norm[name, t] = r.norm_expr + for i in 1:depth + targets.beta[name, i, t] = r.beta_var[i] + end + targets.delta[name, t] = r.delta_var + targets.disc_expr[name, t] = r.disc_expression + targets.disc_cons[name, t] = r.disc_constraint + return +end + +function _alloc_binary_continuous_product_targets!( + container::OptimizationContainer, ::Type{C}, name_axis, time_axis, depth::Int, meta; + tighten::Bool, +) where {C <: IS.InfrastructureSystemsComponent} + u_target = add_variable_container!( + container, NMDTBinaryContinuousProductVariable, C, + name_axis, 1:depth, time_axis; meta, + ) + mc_meta = meta * "_bc" + mc_upper_target = add_constraints_container!( + container, McCormickUpperConstraint, C, + name_axis, 1:depth, 1:2, time_axis; meta = mc_meta, + ) + mc_lower_target = if tighten + nothing + else + add_constraints_container!( + container, McCormickUpperConstraint, C, + name_axis, 1:depth, 1:2, time_axis; meta = mc_meta * "_lb", + ) + end + result_expr_target = add_expression_container!( + container, NMDTBinaryContinuousProductExpression, C, + name_axis, time_axis; meta, + ) + return ( + u = u_target, + mc_upper = mc_upper_target, + mc_lower = mc_lower_target, + result_expr = result_expr_target, + ) +end + +function _write_binary_continuous_product_cell!(targets, name, t, r, depth::Int) + for i in 1:depth + targets.u[name, i, t] = r.u_var[i] + targets.mc_upper[name, i, 1, t] = r.mccormick_upper.c1[i] + targets.mc_upper[name, i, 2, t] = r.mccormick_upper.c2[i] + if targets.mc_lower !== nothing + targets.mc_lower[name, i, 1, t] = r.mccormick_lower.c1[i] + targets.mc_lower[name, i, 2, t] = r.mccormick_lower.c2[i] + end + end + targets.result_expr[name, t] = r.result_expression + return +end + +function _alloc_residual_product_targets!( + container::OptimizationContainer, ::Type{C}, name_axis, time_axis, meta; + tighten::Bool, +) where {C <: IS.InfrastructureSystemsComponent} + z_target = add_variable_container!( + container, NMDTResidualProductVariable, C, name_axis, time_axis; meta, + ) + res_meta = meta * "_res" + mc_target = if tighten + add_constraints_container!( + container, McCormickUpperConstraint, C, + name_axis, 1:2, time_axis; meta = res_meta, + ) + else + add_constraints_container!( + container, McCormickConstraint, C, + name_axis, 1:4, time_axis; meta = res_meta, + ) + end + return (z = z_target, mc = mc_target, tighten = tighten) +end + +function _write_residual_product_cell!(targets, name, t, r) + targets.z[name, t] = r.z_var + mc = r.mccormick_constraints + if targets.tighten + targets.mc[name, 1, t] = mc.upper_1 + targets.mc[name, 2, t] = mc.upper_2 + else + targets.mc[name, 1, t] = mc.upper_1 + targets.mc[name, 2, t] = mc.upper_2 + targets.mc[name, 3, t] = mc.lower_1 + targets.mc[name, 4, t] = mc.lower_2 + end + return +end diff --git a/src/approximations/nmdt_quadratic.jl b/src/approximations/nmdt_quadratic.jl new file mode 100644 index 00000000..43f1ae6d --- /dev/null +++ b/src/approximations/nmdt_quadratic.jl @@ -0,0 +1,315 @@ +# NMDT (Normalized Multiparametric Disaggregation Technique) approximations +# for x². Two variants: +# NMDTQuadConfig — single discretization on x. +# DNMDTQuadConfig — double NMDT: convex combination of (x discretized) and +# (x discretized) again, with shared β·δ residual. +# +# Both normalize x to xh ∈ [0,1], discretize via L binary digits + residual, +# linearize the binary-continuous products with McCormick envelopes, and +# reassemble x² from the normalized components. Optional epigraph Q^{L1} +# lower-bound tightening on xh². +# Reference: Teles, Castro, Matos (2013), Multiparametric disaggregation +# technique for global optimization of polynomial programming problems. + +""" +Config for single-NMDT quadratic approximation. + +# Fields +- `depth::Int`: number of binary discretization levels L. +- `epigraph_depth::Int`: LP tightening depth via epigraph Q^{L1} lower bound; + 0 to disable (default 3·depth). +""" +struct NMDTQuadConfig <: QuadraticApproxConfig + depth::Int + epigraph_depth::Int +end +function NMDTQuadConfig(depth::Int) + return NMDTQuadConfig(depth, 3 * depth) +end + +""" +Config for double-NMDT quadratic approximation. + +# Fields +- `depth::Int`: number of binary discretization levels L. +- `epigraph_depth::Int`: LP tightening depth via epigraph Q^{L1} lower bound; + 0 to disable (default 3·depth). +""" +struct DNMDTQuadConfig <: QuadraticApproxConfig + depth::Int + epigraph_depth::Int +end +function DNMDTQuadConfig(depth::Int) + return DNMDTQuadConfig(depth, 3 * depth) +end + +# --- Scalar build (pure JuMP) --- + +""" + build_quadratic_approx(config::NMDTQuadConfig, model, x, x_min, x_max) + +Scalar form: approximate x² via single NMDT for one cell. Discretize xh, +build the binary-continuous product β·xh, the residual product δ·xh, and +reassemble x². When `epigraph_depth > 0`, also build an epigraph lower +bound on xh² and tighten with `x²_approx ≥ epi`. + +Returns `(; approximation, discretization, bx_xh_product, residual_product, +tightening)`. +""" +function build_quadratic_approx( + config::NMDTQuadConfig, + model::JuMP.Model, + x::JuMP.AbstractJuMPScalar, + x_min::Float64, + x_max::Float64, +) + tighten = config.epigraph_depth > 0 + disc = build_discretization(model, x, x_min, x_max, config.depth) + bx_xh = build_binary_continuous_product( + model, disc.beta_var, disc.norm_expr, 0.0, 1.0, config.depth; tighten, + ) + dz = build_residual_product( + model, disc.delta_var, disc.norm_expr, 1.0, config.depth; tighten, + ) + approximation = build_assembled_product( + model, [bx_xh.result_expression], dz.z_var, + disc.norm_expr, disc.norm_expr, x_min, x_max, x_min, x_max, + ) + tightening = if tighten + epi = build_quadratic_approx( + EpigraphQuadConfig(config.epigraph_depth), model, disc.norm_expr, 0.0, + 1.0, + ) + tcon = JuMP.@constraint(model, approximation >= epi.approximation) + (; epigraph = epi, constraint = tcon) + else + nothing + end + return (; + approximation, + discretization = disc, + bx_xh_product = bx_xh, + residual_product = dz, + tightening, + ) +end + +""" + build_quadratic_approx(config::DNMDTQuadConfig, model, x, x_min, x_max) + +Scalar form: approximate x² via double NMDT for one cell — convex +combination of two NMDT estimates with shared discretization on x. + +Returns the same NamedTuple shape as `NMDTQuadConfig` plus an additional +`bx_dx_product` field for the second binary-continuous product step. +""" +function build_quadratic_approx( + config::DNMDTQuadConfig, + model::JuMP.Model, + x::JuMP.AbstractJuMPScalar, + x_min::Float64, + x_max::Float64, +) + tighten = config.epigraph_depth > 0 + disc = build_discretization(model, x, x_min, x_max, config.depth) + bx_xh = build_binary_continuous_product( + model, disc.beta_var, disc.norm_expr, 0.0, 1.0, config.depth; tighten, + ) + bx_dx = build_binary_continuous_product( + model, disc.beta_var, disc.delta_var, + 0.0, 2.0^(-config.depth), config.depth; tighten, + ) + dz = build_residual_product( + model, disc.delta_var, disc.delta_var, + 2.0^(-config.depth), config.depth; tighten, + ) + approximation = build_assembled_dnmdt( + model, + bx_xh.result_expression, + bx_dx.result_expression, + bx_xh.result_expression, + bx_dx.result_expression, + dz.z_var, + disc.norm_expr, disc.norm_expr, + x_min, x_max, x_min, x_max, + ) + tightening = if tighten + epi = build_quadratic_approx( + EpigraphQuadConfig(config.epigraph_depth), model, disc.norm_expr, 0.0, + 1.0, + ) + tcon = JuMP.@constraint(model, approximation >= epi.approximation) + (; epigraph = epi, constraint = tcon) + else + nothing + end + return (; + approximation, + discretization = disc, + bx_xh_product = bx_xh, + bx_dx_product = bx_dx, + residual_product = dz, + tightening, + ) +end + +# --- IOM allocation + per-cell write helpers --- + +function _alloc_nmdt_tightening_targets!( + container::OptimizationContainer, ::Type{C}, name_axis, time_axis, epi_depth::Int, meta, +) where {C <: IS.InfrastructureSystemsComponent} + epi_targets = _alloc_epigraph_targets!( + container, C, name_axis, time_axis, epi_depth, meta * "_epi", + ) + tighten_cons = add_constraints_container!( + container, NMDTTightenConstraint, C, name_axis, time_axis; meta, + ) + return (epi = epi_targets, tighten = tighten_cons) +end + +function _write_nmdt_tightening_cell!(targets, name, t, tightening, epi_depth::Int) + _write_epigraph_cell!(targets.epi, name, t, tightening.epigraph, epi_depth) + targets.tighten[name, t] = tightening.constraint + return +end + +""" + add_quadratic_approx!(config::NMDTQuadConfig, container, ::Type{C}, x_var, x_bounds, meta) + +Allocate discretization + binary-continuous-product + residual-product +containers, plus when tightening is enabled the epigraph + tightening +constraint containers. Loop `(name, t)`. +""" +function add_quadratic_approx!( + config::NMDTQuadConfig, + container::OptimizationContainer, + ::Type{C}, + x_var, + x_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + name_axis = axes(x_var, 1) + time_axis = axes(x_var, 2) + depth = config.depth + tighten = config.epigraph_depth > 0 + @assert length(name_axis) == length(x_bounds) + for b in x_bounds + @assert b.max > b.min + end + + model = get_jump_model(container) + + disc_targets = _alloc_discretization_targets!( + container, C, name_axis, time_axis, depth, meta, + ) + bx_targets = _alloc_binary_continuous_product_targets!( + container, C, name_axis, time_axis, depth, meta; tighten, + ) + res_targets = _alloc_residual_product_targets!( + container, C, name_axis, time_axis, meta; tighten, + ) + approx_target = add_expression_container!( + container, QuadraticExpression, C, name_axis, time_axis; meta, + ) + tighten_targets = if tighten + _alloc_nmdt_tightening_targets!( + container, C, name_axis, time_axis, config.epigraph_depth, meta, + ) + else + nothing + end + + for (i, name) in enumerate(name_axis) + xmn, xmx = x_bounds[i].min, x_bounds[i].max + for t in time_axis + r = build_quadratic_approx(config, model, x_var[name, t], xmn, xmx) + _write_discretization_cell!(disc_targets, name, t, r.discretization, depth) + _write_binary_continuous_product_cell!( + bx_targets, + name, + t, + r.bx_xh_product, + depth, + ) + _write_residual_product_cell!(res_targets, name, t, r.residual_product) + approx_target[name, t] = r.approximation + if tighten + _write_nmdt_tightening_cell!( + tighten_targets, name, t, r.tightening, config.epigraph_depth, + ) + end + end + end + return approx_target +end + +""" + add_quadratic_approx!(config::DNMDTQuadConfig, container, ::Type{C}, x_var, x_bounds, meta) + +Allocate two binary-continuous-product container sets +(`meta * "_bx_xh"` and `meta * "_bx_dx"`) plus the rest of the NMDT pieces. +""" +function add_quadratic_approx!( + config::DNMDTQuadConfig, + container::OptimizationContainer, + ::Type{C}, + x_var, + x_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + name_axis = axes(x_var, 1) + time_axis = axes(x_var, 2) + depth = config.depth + tighten = config.epigraph_depth > 0 + @assert length(name_axis) == length(x_bounds) + for b in x_bounds + @assert b.max > b.min + end + + model = get_jump_model(container) + + disc_targets = _alloc_discretization_targets!( + container, C, name_axis, time_axis, depth, meta, + ) + bx_xh_targets = _alloc_binary_continuous_product_targets!( + container, C, name_axis, time_axis, depth, meta * "_bx_xh"; tighten, + ) + bx_dx_targets = _alloc_binary_continuous_product_targets!( + container, C, name_axis, time_axis, depth, meta * "_bx_dx"; tighten, + ) + res_targets = _alloc_residual_product_targets!( + container, C, name_axis, time_axis, meta; tighten, + ) + approx_target = add_expression_container!( + container, QuadraticExpression, C, name_axis, time_axis; meta, + ) + tighten_targets = if tighten + _alloc_nmdt_tightening_targets!( + container, C, name_axis, time_axis, config.epigraph_depth, meta, + ) + else + nothing + end + + for (i, name) in enumerate(name_axis) + xmn, xmx = x_bounds[i].min, x_bounds[i].max + for t in time_axis + r = build_quadratic_approx(config, model, x_var[name, t], xmn, xmx) + _write_discretization_cell!(disc_targets, name, t, r.discretization, depth) + _write_binary_continuous_product_cell!( + bx_xh_targets, name, t, r.bx_xh_product, depth, + ) + _write_binary_continuous_product_cell!( + bx_dx_targets, name, t, r.bx_dx_product, depth, + ) + _write_residual_product_cell!(res_targets, name, t, r.residual_product) + approx_target[name, t] = r.approximation + if tighten + _write_nmdt_tightening_cell!( + tighten_targets, name, t, r.tightening, config.epigraph_depth, + ) + end + end + end + return approx_target +end diff --git a/src/approximations/no_approx_bilinear.jl b/src/approximations/no_approx_bilinear.jl new file mode 100644 index 00000000..48b5999e --- /dev/null +++ b/src/approximations/no_approx_bilinear.jl @@ -0,0 +1,89 @@ +# No-op bilinear approximation: returns the exact x·y as a JuMP.QuadExpr. +# For NLP-capable solvers or testing. + +"No-op config for bilinear approximation: returns exact x·y as a QuadExpr." +struct NoBilinearApproxConfig <: BilinearApproxConfig end + +""" + build_bilinear_approx(::NoBilinearApproxConfig, model, x, y, x_min, x_max, y_min, y_max) + +Scalar form: return `(; approximation = x*y)` for a single JuMP scalar pair. +Bounds accepted for signature parity, unused. +""" +function build_bilinear_approx( + ::NoBilinearApproxConfig, + model::JuMP.Model, + x::JuMP.AbstractJuMPScalar, + y::JuMP.AbstractJuMPScalar, + x_min::Float64, + x_max::Float64, + y_min::Float64, + y_max::Float64, +) + return (; approximation = x * y) +end + +""" + add_bilinear_approx!(::NoBilinearApproxConfig, container, ::Type{C}, x_var, y_var, x_bounds, y_bounds, meta) + +Allocate a `BilinearProductExpression` container with axes `(name, t)`, +loop, and write the exact `x*y` per cell. +""" +function add_bilinear_approx!( + ::NoBilinearApproxConfig, + container::OptimizationContainer, + ::Type{C}, + x_var, + y_var, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + name_axis = axes(x_var, 1) + time_axis = axes(x_var, 2) + @assert length(name_axis) == length(x_bounds) + @assert length(name_axis) == length(y_bounds) + model = get_jump_model(container) + target = add_expression_container!( + container, BilinearProductExpression, C, name_axis, time_axis; + meta, expr_type = JuMP.QuadExpr, + ) + for (i, name) in enumerate(name_axis) + xmn, xmx = x_bounds[i].min, x_bounds[i].max + ymn, ymx = y_bounds[i].min, y_bounds[i].max + for t in time_axis + r = build_bilinear_approx( + NoBilinearApproxConfig(), model, + x_var[name, t], y_var[name, t], + xmn, xmx, ymn, ymx, + ) + target[name, t] = r.approximation + end + end + return target +end + +""" + add_bilinear_approx!(::NoBilinearApproxConfig, container, ::Type{C}, xsq, ysq, x_var, y_var, x_bounds, y_bounds, meta) + +Precomputed-form entrypoint: signature-compatible with the precomputed-form +of `Bin2Config` / `HybSConfig`, so a caller can swap configs without +changing the call site. `xsq` and `ysq` are accepted but ignored — the +no-op approximation just returns the exact `x*y` product. +""" +function add_bilinear_approx!( + ::NoBilinearApproxConfig, + container::OptimizationContainer, + ::Type{C}, + xsq, + ysq, + x_var, + y_var, + x_bounds::Vector{MinMax}, + y_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + return add_bilinear_approx!( + NoBilinearApproxConfig(), container, C, x_var, y_var, x_bounds, y_bounds, meta, + ) +end diff --git a/src/approximations/no_approx_quadratic.jl b/src/approximations/no_approx_quadratic.jl new file mode 100644 index 00000000..bffbda37 --- /dev/null +++ b/src/approximations/no_approx_quadratic.jl @@ -0,0 +1,59 @@ +# No-op quadratic approximation: returns the exact x² as a JuMP.QuadExpr. +# For NLP-capable solvers or testing. + +"No-op config for quadratic approximation: returns exact x² as a QuadExpr." +struct NoQuadApproxConfig <: QuadraticApproxConfig end + +""" + build_quadratic_approx(::NoQuadApproxConfig, model, x, x_min, x_max) + +Scalar form: return `(; approximation = x*x)` for a single JuMP scalar. +Bounds accepted for signature parity with other quadratic methods, unused. +""" +function build_quadratic_approx( + ::NoQuadApproxConfig, + model::JuMP.Model, + x::JuMP.AbstractJuMPScalar, + x_min::Float64, + x_max::Float64, +) + return (; approximation = x * x) +end + +""" + add_quadratic_approx!(::NoQuadApproxConfig, container, ::Type{C}, x_var, x_bounds, meta) + +Allocate a `QuadraticExpression` container with axes `(name, t)`, loop, +call the scalar build per cell, write the exact `x*x` per cell. +""" +function add_quadratic_approx!( + ::NoQuadApproxConfig, + container::OptimizationContainer, + ::Type{C}, + x_var, + x_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + name_axis = axes(x_var, 1) + time_axis = axes(x_var, 2) + @assert length(name_axis) == length(x_bounds) + model = get_jump_model(container) + target = add_expression_container!( + container, QuadraticExpression, C, name_axis, time_axis; + meta, expr_type = JuMP.QuadExpr, + ) + for (i, name) in enumerate(name_axis) + xmn, xmx = x_bounds[i].min, x_bounds[i].max + for t in time_axis + r = build_quadratic_approx( + NoQuadApproxConfig(), + model, + x_var[name, t], + xmn, + xmx, + ) + target[name, t] = r.approximation + end + end + return target +end diff --git a/src/approximations/pwl_utils.jl b/src/approximations/pwl_utils.jl new file mode 100644 index 00000000..03b8ea6f --- /dev/null +++ b/src/approximations/pwl_utils.jl @@ -0,0 +1,30 @@ +# Pure helpers for piecewise linear approximation breakpoint generation. + +""" + _get_breakpoints_for_pwl_function(min_val, max_val, f; num_segments = DEFAULT_INTERPOLATION_LENGTH) + +Generate `num_segments + 1` equally-spaced breakpoints over `[min_val, max_val]` +and evaluate `f` at each. Returns `(x_bkpts, y_bkpts)`. +""" +function _get_breakpoints_for_pwl_function( + min_val::Float64, + max_val::Float64, + f; + num_segments = DEFAULT_INTERPOLATION_LENGTH, +) + num_bkpts = num_segments + 1 + step = (max_val - min_val) / num_segments + x_bkpts = Vector{Float64}(undef, num_bkpts) + y_bkpts = Vector{Float64}(undef, num_bkpts) + x_bkpts[1] = min_val + y_bkpts[1] = f(min_val) + for i in 1:num_segments + x = min_val + step * i + x_bkpts[i + 1] = x + y_bkpts[i + 1] = f(x) + end + return x_bkpts, y_bkpts +end + +"Returns x² (used as the default function for PWL breakpoint generation)." +_square(x::Float64) = x * x diff --git a/src/approximations/pwmcc_cuts.jl b/src/approximations/pwmcc_cuts.jl new file mode 100644 index 00000000..b54627c4 --- /dev/null +++ b/src/approximations/pwmcc_cuts.jl @@ -0,0 +1,159 @@ +# Piecewise McCormick (PWMCC) cuts for concave terms in PWL approximations of x². +# Adds K local chord upper bounds on the v² PWL approximation by partitioning +# v's domain into K sub-intervals. The LP gap shrinks from Delta²/4 to +# Delta²/(4·K²). These cuts supplement (do not replace) the underlying PWL +# (SOS2 or manual-SOS2) constraints. + +# --- Container key types --- + +"Binary interval selector for piecewise McCormick cuts." +struct PiecewiseMcCormickBinary <: SparseVariableType end + +"Disaggregated variable for piecewise McCormick cuts." +struct PiecewiseMcCormickDisaggregated <: SparseVariableType end + +"Selector sum constraint: sum_k delta_k = 1." +struct PiecewiseMcCormickSelectorSum <: ConstraintType end + +"Disaggregation linking constraint: v = sum_k v^d_k." +struct PiecewiseMcCormickLinking <: ConstraintType end + +"Interval activation lower bound: t_{k-1} * delta_k <= v^d_k." +struct PiecewiseMcCormickIntervalLB <: ConstraintType end + +"Interval activation upper bound: v^d_k <= t_k * delta_k." +struct PiecewiseMcCormickIntervalUB <: ConstraintType end + +"Piecewise McCormick chord upper-bound constraint on v² approximation." +struct PiecewiseMcCormickChordUB <: ConstraintType end + +"Piecewise McCormick tangent lower-bound constraint (left endpoint)." +struct PiecewiseMcCormickTangentLBL <: ConstraintType end + +"Piecewise McCormick tangent lower-bound constraint (right endpoint)." +struct PiecewiseMcCormickTangentLBR <: ConstraintType end + +""" + build_pwmcc_concave_cuts(model, v, q_expr, v_min, v_max, K) + +Scalar form: build the K-segment PWMCC cuts at a single cell. `v` is the +JuMP scalar variable and `q_expr` is the scalar expression for the +existing PWL v² approximation at this cell. + +Returns a NamedTuple with `(delta_var, vd_var, selector_constraint, +linking_constraint, interval_lb_constraints, interval_ub_constraints, +chord_ub_constraint, tangent_lb_l_constraint, tangent_lb_r_constraint)`. +""" +function build_pwmcc_concave_cuts( + model::JuMP.Model, + v::JuMP.AbstractJuMPScalar, + q_expr, + v_min::Float64, + v_max::Float64, + K::Int, +) + @assert K >= 1 + @assert v_min < v_max + + brk = [v_min + k * (v_max - v_min) / K for k in 0:K] + + delta_var = JuMP.@variable( + model, [k = 1:K], binary = true, base_name = "PwMcCBin", + ) + vd_var = JuMP.@variable(model, [k = 1:K], base_name = "PwMcCDis") + + selector_con = JuMP.@constraint(model, sum(delta_var[k] for k in 1:K) == 1.0) + linking_con = JuMP.@constraint(model, sum(vd_var[k] for k in 1:K) == v) + + interval_lb = JuMP.@constraint( + model, [k = 1:K], brk[k] * delta_var[k] <= vd_var[k], + ) + interval_ub = JuMP.@constraint( + model, [k = 1:K], vd_var[k] <= brk[k + 1] * delta_var[k], + ) + chord_ub = JuMP.@constraint( + model, + q_expr <= sum( + (brk[k] + brk[k + 1]) * vd_var[k] - + brk[k] * brk[k + 1] * delta_var[k] for k in 1:K + ), + ) + tangent_lb_l = JuMP.@constraint( + model, + q_expr >= sum( + 2.0 * brk[k] * vd_var[k] - brk[k]^2 * delta_var[k] for k in 1:K + ), + ) + tangent_lb_r = JuMP.@constraint( + model, + q_expr >= sum( + 2.0 * brk[k + 1] * vd_var[k] - brk[k + 1]^2 * delta_var[k] for k in 1:K + ), + ) + + return (; + delta_var, + vd_var, + selector_constraint = selector_con, + linking_constraint = linking_con, + interval_lb_constraints = interval_lb, + interval_ub_constraints = interval_ub, + chord_ub_constraint = chord_ub, + tangent_lb_l_constraint = tangent_lb_l, + tangent_lb_r_constraint = tangent_lb_r, + ) +end + +# --- Allocation + per-cell write helpers (used by SOS2 adapters) --- + +function _alloc_pwmcc_targets!( + container::OptimizationContainer, ::Type{C}, name_axis, time_axis, K::Int, meta, +) where {C <: IS.InfrastructureSystemsComponent} + return ( + delta = add_variable_container!( + container, PiecewiseMcCormickBinary, C, name_axis, 1:K, time_axis; meta, + ), + vd = add_variable_container!( + container, PiecewiseMcCormickDisaggregated, C, + name_axis, 1:K, time_axis; meta, + ), + selector = add_constraints_container!( + container, PiecewiseMcCormickSelectorSum, C, name_axis, time_axis; meta, + ), + linking = add_constraints_container!( + container, PiecewiseMcCormickLinking, C, name_axis, time_axis; meta, + ), + interval_lb = add_constraints_container!( + container, PiecewiseMcCormickIntervalLB, C, + name_axis, 1:K, time_axis; meta, + ), + interval_ub = add_constraints_container!( + container, PiecewiseMcCormickIntervalUB, C, + name_axis, 1:K, time_axis; meta, + ), + chord = add_constraints_container!( + container, PiecewiseMcCormickChordUB, C, name_axis, time_axis; meta, + ), + tangent_l = add_constraints_container!( + container, PiecewiseMcCormickTangentLBL, C, name_axis, time_axis; meta, + ), + tangent_r = add_constraints_container!( + container, PiecewiseMcCormickTangentLBR, C, name_axis, time_axis; meta, + ), + ) +end + +function _write_pwmcc_cell!(targets, name, t, r, K::Int) + for k in 1:K + targets.delta[name, k, t] = r.delta_var[k] + targets.vd[name, k, t] = r.vd_var[k] + targets.interval_lb[name, k, t] = r.interval_lb_constraints[k] + targets.interval_ub[name, k, t] = r.interval_ub_constraints[k] + end + targets.selector[name, t] = r.selector_constraint + targets.linking[name, t] = r.linking_constraint + targets.chord[name, t] = r.chord_ub_constraint + targets.tangent_l[name, t] = r.tangent_lb_l_constraint + targets.tangent_r[name, t] = r.tangent_lb_r_constraint + return +end diff --git a/src/approximations/sawtooth.jl b/src/approximations/sawtooth.jl new file mode 100644 index 00000000..772da1f5 --- /dev/null +++ b/src/approximations/sawtooth.jl @@ -0,0 +1,216 @@ +# Sawtooth MIP approximation of x² for use in constraints. +# Uses recursive tooth function compositions with O(log(1/ε)) binary variables. +# Reference: Beach, Burlacu, Hager, Hildebrand (2024). + +"Binary variables (α₁, …, α_L) for sawtooth quadratic approximation." +struct SawtoothBinaryVariable <: VariableType end + +"Variable result in tightened version." +struct SawtoothTightenedVariable <: VariableType end + +"Constrains g_j based on g_{j-1}." +struct SawtoothMIPConstraint <: ConstraintType end + +"Bounds tightened sawtooth variable." +struct SawtoothTightenedConstraint <: ConstraintType end + +""" +Config for sawtooth MIP quadratic approximation. + +# Fields +- `depth::Int`: recursion depth L; uses L binary variables for 2^L + 1 breakpoints. +- `epigraph_depth::Int`: LP tightening depth via epigraph Q^{L1} lower bound; + 0 to disable (default 0). +""" +struct SawtoothQuadConfig <: QuadraticApproxConfig + depth::Int + epigraph_depth::Int +end +function SawtoothQuadConfig(depth::Int) + return SawtoothQuadConfig(depth, 0) +end + +""" + build_quadratic_approx(config::SawtoothQuadConfig, model, x, x_min, x_max) + +Scalar form: PWL approximation of x² for a single JuMP scalar `x` with +bounds `[x_min, x_max]`, using `depth` binary variables. If +`config.epigraph_depth > 0`, also builds an epigraph Q^{L1} lower bound and +tightens with z ≤ sawtooth_upper and z ≥ epigraph. + +Returns a NamedTuple with `(approximation, g_var, alpha_var, link_constraint, +mip_constraints, tightening)` where `tightening` is `nothing` or +`(; z_var, constraints :: 1D over 1:2, epigraph)`. +""" +function build_quadratic_approx( + config::SawtoothQuadConfig, + model::JuMP.Model, + x::JuMP.AbstractJuMPScalar, + x_min::Float64, + x_max::Float64, +) + @assert config.depth >= 1 + @assert x_max > x_min + + depth = config.depth + delta = x_max - x_min + + g_var = JuMP.@variable( + model, [j = 0:depth], + lower_bound = 0.0, upper_bound = 1.0, + base_name = "SawtoothAux", + ) + alpha_var = JuMP.@variable( + model, [j = 1:depth], + binary = true, + base_name = "SawtoothBin", + ) + + link_con = JuMP.@constraint(model, g_var[0] == (x - x_min) / delta) + + mip_a = JuMP.@constraint( + model, [j = 1:depth], g_var[j] <= 2.0 * g_var[j - 1], + ) + mip_b = JuMP.@constraint( + model, [j = 1:depth], g_var[j] <= 2.0 * (1.0 - g_var[j - 1]), + ) + mip_c = JuMP.@constraint( + model, [j = 1:depth], g_var[j] >= 2.0 * (g_var[j - 1] - alpha_var[j]), + ) + mip_d = JuMP.@constraint( + model, [j = 1:depth], g_var[j] >= 2.0 * (alpha_var[j] - g_var[j - 1]), + ) + mip_cons = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, 1:depth, 1:4) + @views mip_cons.data[:, 1] .= mip_a + @views mip_cons.data[:, 2] .= mip_b + @views mip_cons.data[:, 3] .= mip_c + @views mip_cons.data[:, 4] .= mip_d + + x_sq_approx = JuMP.@expression( + model, + scale_back_g_basis(x_min, delta, g_var, 1:depth), + ) + + if config.epigraph_depth > 0 + epi = build_quadratic_approx( + EpigraphQuadConfig(config.epigraph_depth), model, x, x_min, x_max, + ) + z_min = (x_min <= 0.0 <= x_max) ? 0.0 : min(x_min^2, x_max^2) + z_max = max(x_min^2, x_max^2) + z_var = JuMP.@variable( + model, lower_bound = z_min, upper_bound = z_max, + base_name = "TightenedSawtooth", + ) + tight_a = JuMP.@constraint(model, z_var <= x_sq_approx) + tight_b = JuMP.@constraint(model, z_var >= epi.approximation) + tight_cons = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, 1:2) + tight_cons[1] = tight_a + tight_cons[2] = tight_b + approximation = JuMP.@expression(model, 1.0 * z_var) + tightening = (; z_var, constraints = tight_cons, epigraph = epi) + return (; + approximation, + g_var, + alpha_var, + link_constraint = link_con, + mip_constraints = mip_cons, + tightening, + ) + end + + return (; + approximation = x_sq_approx, + g_var, + alpha_var, + link_constraint = link_con, + mip_constraints = mip_cons, + tightening = nothing, + ) +end + +""" + add_quadratic_approx!(config::SawtoothQuadConfig, container, ::Type{C}, x_var, x_bounds, meta) + +Allocate sawtooth containers (g, α, link, mip, approximation) plus, when +`config.epigraph_depth > 0`, the tightened-z + 2-constraint containers AND +the full set of epigraph containers under `meta * "_lb"`. Loop `(name, t)`. +""" +function add_quadratic_approx!( + config::SawtoothQuadConfig, + container::OptimizationContainer, + ::Type{C}, + x_var, + x_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + name_axis = axes(x_var, 1) + time_axis = axes(x_var, 2) + depth = config.depth + @assert depth >= 1 + @assert length(name_axis) == length(x_bounds) + for b in x_bounds + @assert b.max > b.min + end + + model = get_jump_model(container) + + g_target = add_variable_container!( + container, SawtoothAuxVariable, C, name_axis, 0:depth, time_axis; meta, + ) + alpha_target = add_variable_container!( + container, SawtoothBinaryVariable, C, name_axis, 1:depth, time_axis; meta, + ) + link_target = add_constraints_container!( + container, SawtoothLinkingConstraint, C, name_axis, time_axis; meta, + ) + mip_target = add_constraints_container!( + container, SawtoothMIPConstraint, C, name_axis, 1:depth, 1:4, time_axis; meta, + ) + approx_target = add_expression_container!( + container, QuadraticExpression, C, name_axis, time_axis; meta, + ) + + tighten = config.epigraph_depth > 0 + local st_z_target, st_tight_target, epi_targets + local epi_depth::Int + if tighten + st_z_target = add_variable_container!( + container, SawtoothTightenedVariable, C, name_axis, time_axis; meta, + ) + st_tight_target = add_constraints_container!( + container, SawtoothTightenedConstraint, C, name_axis, 1:2, time_axis; meta, + ) + epi_depth = config.epigraph_depth + epi_targets = _alloc_epigraph_targets!( + container, C, name_axis, time_axis, epi_depth, meta * "_lb", + ) + end + + for (i, name) in enumerate(name_axis) + xmn, xmx = x_bounds[i].min, x_bounds[i].max + for t in time_axis + r = build_quadratic_approx(config, model, x_var[name, t], xmn, xmx) + for j in 0:depth + g_target[name, j, t] = r.g_var[j] + end + for j in 1:depth + alpha_target[name, j, t] = r.alpha_var[j] + end + link_target[name, t] = r.link_constraint + for j in 1:depth, k in 1:4 + mip_target[name, j, k, t] = r.mip_constraints[j, k] + end + approx_target[name, t] = r.approximation + + if tighten + tt = r.tightening + st_z_target[name, t] = tt.z_var + for k in 1:2 + st_tight_target[name, k, t] = tt.constraints[k] + end + _write_epigraph_cell!(epi_targets, name, t, tt.epigraph, epi_depth) + end + end + end + return approx_target +end diff --git a/src/approximations/solver_sos2.jl b/src/approximations/solver_sos2.jl new file mode 100644 index 00000000..89e5886c --- /dev/null +++ b/src/approximations/solver_sos2.jl @@ -0,0 +1,176 @@ +# Solver-native SOS2 piecewise linear approximation of x² for use in constraints. +# Uses solver-native MOI.SOS2 constraints for adjacency enforcement among λ weights. + +"Lambda (λ) convex-combination weight variables for SOS2 quadratic approximation." +struct QuadraticVariable <: SparseVariableType end + +"Links x to the weighted sum of breakpoints in SOS2 quadratic approximation." +struct SOS2LinkingConstraint <: ConstraintType end + +"Expression for the weighted sum of breakpoints Σ λ_i · x_i linking x to λ variables." +struct SOS2LinkingExpression <: ExpressionType end + +"Ensures the sum of λ weights equals 1 in SOS2 quadratic approximation." +struct SOS2NormConstraint <: ConstraintType end + +"Expression for the normalization sum Σ λ_i in SOS2 quadratic approximation." +struct SOS2NormExpression <: ExpressionType end + +"Solver-native MOI.SOS2 adjacency constraint on lambda variables." +struct SolverSOS2Constraint <: ConstraintType end + +""" +Config for solver-native SOS2 quadratic approximation (MOI.SOS2 adjacency). + +# Fields +- `depth::Int`: number of PWL segments (breakpoints = depth + 1). +- `pwmcc_segments::Int`: number of piecewise McCormick cut partitions; + 0 to disable (default 4). +""" +struct SolverSOS2QuadConfig <: QuadraticApproxConfig + depth::Int + pwmcc_segments::Int +end +function SolverSOS2QuadConfig(depth::Int) + return SolverSOS2QuadConfig(depth, 4) +end + +""" + build_quadratic_approx(config::SolverSOS2QuadConfig, model, x, x_min, x_max) + +Scalar form: PWL approximation of x² for a single JuMP scalar `x` using +solver-native MOI.SOS2 adjacency. If `config.pwmcc_segments > 0`, also +adds piecewise McCormick concave cuts per cell. + +Returns a NamedTuple with `(approximation, lambda, link_constraint, +norm_constraint, sos_constraint, link_expression, norm_expression, pwmcc)`. +""" +function build_quadratic_approx( + config::SolverSOS2QuadConfig, + model::JuMP.Model, + x::JuMP.AbstractJuMPScalar, + x_min::Float64, + x_max::Float64, +) + @assert x_max > x_min + n_points = config.depth + 1 + x_bkpts, x_sq_bkpts = _get_breakpoints_for_pwl_function( + 0.0, 1.0, _square; num_segments = config.depth, + ) + lx = x_max - x_min + + lambda = JuMP.@variable( + model, [i = 1:n_points], + lower_bound = 0.0, upper_bound = 1.0, + base_name = "QuadraticVariable", + ) + link_expr = JuMP.@expression( + model, sum(x_bkpts[i] * lambda[i] for i in 1:n_points), + ) + link_con = JuMP.@constraint(model, (x - x_min) / lx == link_expr) + norm_expr = JuMP.@expression(model, sum(lambda[i] for i in 1:n_points)) + norm_con = JuMP.@constraint(model, norm_expr == 1.0) + sos_con = JuMP.@constraint( + model, [lambda[i] for i in 1:n_points] in MOI.SOS2(collect(1:n_points)), + ) + approximation = JuMP.@expression( + model, + lx * lx * sum(x_sq_bkpts[i] * lambda[i] for i in 1:n_points) + + 2.0 * x_min * x - x_min * x_min, + ) + + pwmcc = if config.pwmcc_segments > 0 + build_pwmcc_concave_cuts( + model, x, approximation, x_min, x_max, config.pwmcc_segments, + ) + else + nothing + end + + return (; + approximation, + lambda, + link_constraint = link_con, + norm_constraint = norm_con, + sos_constraint = sos_con, + link_expression = link_expr, + norm_expression = norm_expr, + pwmcc, + ) +end + +""" + add_quadratic_approx!(config::SolverSOS2QuadConfig, container, ::Type{C}, x_var, x_bounds, meta) + +Allocate SOS2 containers (λ, link/norm/sos constraints, expressions, +approximation) plus, when `config.pwmcc_segments > 0`, the full set of +PWMCC containers under `meta * "_pwmcc"`. Loop `(name, t)`. +""" +function add_quadratic_approx!( + config::SolverSOS2QuadConfig, + container::OptimizationContainer, + ::Type{C}, + x_var, + x_bounds::Vector{MinMax}, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + name_axis = axes(x_var, 1) + time_axis = axes(x_var, 2) + n_points = config.depth + 1 + @assert length(name_axis) == length(x_bounds) + for b in x_bounds + @assert b.max > b.min + end + + model = get_jump_model(container) + + lambda_target = add_variable_container!( + container, QuadraticVariable, C, name_axis, 1:n_points, time_axis; meta, + ) + link_cons_target = add_constraints_container!( + container, SOS2LinkingConstraint, C, name_axis, time_axis; meta, + ) + norm_cons_target = add_constraints_container!( + container, SOS2NormConstraint, C, name_axis, time_axis; meta, + ) + sos_cons_target = add_constraints_container!( + container, SolverSOS2Constraint, C, name_axis, time_axis; meta, + ) + link_expr_target = add_expression_container!( + container, SOS2LinkingExpression, C, name_axis, time_axis; meta, + ) + norm_expr_target = add_expression_container!( + container, SOS2NormExpression, C, name_axis, time_axis; meta, + ) + approx_target = add_expression_container!( + container, QuadraticExpression, C, name_axis, time_axis; meta, + ) + + use_pwmcc = config.pwmcc_segments > 0 + K = config.pwmcc_segments + pwmcc_targets = if use_pwmcc + _alloc_pwmcc_targets!(container, C, name_axis, time_axis, K, meta * "_pwmcc") + else + nothing + end + + for (i, name) in enumerate(name_axis) + xmn, xmx = x_bounds[i].min, x_bounds[i].max + for t in time_axis + r = build_quadratic_approx(config, model, x_var[name, t], xmn, xmx) + for ip in 1:n_points + lambda_target[name, ip, t] = r.lambda[ip] + end + link_cons_target[name, t] = r.link_constraint + norm_cons_target[name, t] = r.norm_constraint + sos_cons_target[name, t] = r.sos_constraint + link_expr_target[name, t] = r.link_expression + norm_expr_target[name, t] = r.norm_expression + approx_target[name, t] = r.approximation + if use_pwmcc + _write_pwmcc_cell!(pwmcc_targets, name, t, r.pwmcc, K) + end + end + end + return approx_target +end diff --git a/src/bilinear_approximations/bin2.jl b/src/bilinear_approximations/bin2.jl deleted file mode 100644 index 5c331cd4..00000000 --- a/src/bilinear_approximations/bin2.jl +++ /dev/null @@ -1,156 +0,0 @@ -# Bin2 separable approximation of bilinear products z = x·y. -# Uses the identity: x·y = (1/2)*((x+y)² − x² - y²). -# Calls existing quadratic approximation functions for p²=(x+y)² - -"Expression container for bilinear product (x·y) approximation results." -struct BilinearProductExpression <: ExpressionType end -"Variable container for bilinear product (x ̇y) approximation results." -struct BilinearProductVariable <: VariableType end -"Expression container for adding variables." -struct VariableSumExpression <: ExpressionType end -"Expression container for subtracting variables." -struct VariableDifferenceExpression <: ExpressionType end -"Constraint container for linking product expressions and variables." -struct BilinearProductLinkingConstraint <: ConstraintType end - -# --- Bilinear approximation config hierarchy --- - -"Abstract supertype for bilinear approximation method configurations." -abstract type BilinearApproxConfig end - -""" -Config for Bin2 bilinear approximation using z = ½((x+y)² − x² − y²). - -# Fields -- `quad_config::QuadraticApproxConfig`: quadratic method used for x², y², and (x+y)² -- `add_mccormick::Bool`: whether to add reformulated McCormick cuts through separable variables (default true) -""" -struct Bin2Config <: BilinearApproxConfig - quad_config::QuadraticApproxConfig - add_mccormick::Bool -end -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_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, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_var, - y_var, - 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_bounds, meta * "_x", - ) - ysq = _add_quadratic_approx!( - config.quad_config, container, C, names, time_steps, - y_var, y_bounds, meta * "_y", - ) - return _add_bilinear_approx!( - config, container, C, names, time_steps, - xsq, ysq, x_var, y_var, - x_bounds, y_bounds, meta, - ) -end - -""" - _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, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - xsq, - ysq, - x_var, - y_var, - 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 (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" - - p_expr = add_expression_container!( - container, - VariableSumExpression, - C, - names, - time_steps; - meta = meta_plus, - ) - for name in names, t in time_steps - p = JuMP.AffExpr(0.0) - add_proportional_to_jump_expression!(p, x_var[name, t], 1.0) - add_proportional_to_jump_expression!(p, y_var[name, t], 1.0) - p_expr[name, t] = p - end - - # Approximate p² = (x+y)² using the provided quadratic config - psq = _add_quadratic_approx!( - config.quad_config, container, C, names, time_steps, - p_expr, p_bounds, meta_plus, - ) - - result_expr = add_expression_container!( - container, - BilinearProductExpression, - C, - names, - time_steps; - meta, - ) - - for name in names, t in time_steps - # z = (1/2) * (p² − x² − y²) - result = result_expr[name, t] = JuMP.AffExpr(0.0) - add_proportional_to_jump_expression!(result, psq[name, t], 0.5) - add_proportional_to_jump_expression!(result, xsq[name, t], -0.5) - add_proportional_to_jump_expression!(result, ysq[name, t], -0.5) - end - - # --- Reformulated McCormick cuts (optional) --- - if config.add_mccormick - _add_reformulated_mccormick!( - container, C, names, time_steps, - x_var, y_var, psq, xsq, ysq, - x_bounds, y_bounds, meta, - ) - end - - return result_expr -end diff --git a/src/bilinear_approximations/hybs.jl b/src/bilinear_approximations/hybs.jl deleted file mode 100644 index 8df74bff..00000000 --- a/src/bilinear_approximations/hybs.jl +++ /dev/null @@ -1,235 +0,0 @@ -# HybS (Hybrid Separable) MIP relaxation for bilinear products z = x·y. -# Combines Bin2 lower bound and Bin3 upper bound with shared sawtooth for x², y² -# and LP-only epigraph for (x+y)², (x−y)². Uses 2L binaries instead of 3L (Bin2). -# Reference: Beach, Burlacu, Bärmann, Hager, Hildebrand (2024), Definition 10. - -"Two-sided HybS bound constraints: Bin2 lower + Bin3 upper." -struct HybSBoundConstraint <: ConstraintType end - -""" -Config for HybS (Hybrid Separable) bilinear approximation. - -Combines Bin2 lower bound and Bin3 upper bound with shared quadratic for x², y² -and LP-only epigraph for (x+y)², (x−y)². - -# Fields -- `quad_config::QuadraticApproxConfig`: quadratic method used for the shared x² and y² terms -- `epigraph_depth::Int`: depth for the epigraph Q^{L1} LP-only approximation of cross-terms (x±y)² -- `add_mccormick::Bool`: whether to add standard McCormick envelope cuts on the product variable (default false) -""" -struct HybSConfig <: BilinearApproxConfig - quad_config::QuadraticApproxConfig - epigraph_depth::Int - add_mccormick::Bool -end -HybSConfig(quad_config::QuadraticApproxConfig, epigraph_depth::Int) = - HybSConfig(quad_config, epigraph_depth, false) - -# --- Unified HybS dispatch methods --- - -""" - _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, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_var, - y_var, - 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_bounds, meta * "_x", - ) - ysq = _add_quadratic_approx!( - config.quad_config, container, C, names, time_steps, - y_var, y_bounds, meta * "_y", - ) - return _add_bilinear_approx!( - config, container, C, names, time_steps, - xsq, ysq, x_var, y_var, - x_bounds, y_bounds, meta, - ) -end - -""" - _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². - -Combines Bin2 and Bin3 separable identities: -- Bin2 lower bound: z ≥ ½(z_p1 − z_x − z_y) where z_p1 lower-bounds (x+y)² -- 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, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - xsq, - ysq, - x_var, - y_var, - x_bounds::Vector{MinMax}, - y_bounds::Vector{MinMax}, - 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) - ] - - jump_model = get_jump_model(container) - - # Meta suffixes for cross-term expressions - meta_p1 = meta * "_plus" - meta_p2 = meta * "_diff" - - p1_expr = add_expression_container!( - container, - VariableSumExpression, - C, - names, - time_steps; - meta = meta_p1, - ) - p2_expr = add_expression_container!( - container, - VariableDifferenceExpression, - C, - names, - time_steps; - meta = meta_p2, - ) - - for name in names, t in time_steps - x = x_var[name, t] - y = y_var[name, t] - - # p1 = x + y - p1 = p1_expr[name, t] = JuMP.AffExpr(0.0) - add_proportional_to_jump_expression!(p1, x, 1.0) - add_proportional_to_jump_expression!(p1, y, 1.0) - - # p2 = x − y - p2 = p2_expr[name, t] = JuMP.AffExpr(0.0) - add_proportional_to_jump_expression!(p2, x, 1.0) - add_proportional_to_jump_expression!(p2, y, -1.0) - end - - # --- Epigraph Q^{L1} lower bound for (x+y)² and (x−y)² (no binaries) --- - epi_cfg = EpigraphQuadConfig(config.epigraph_depth) - zp1_expr = _add_quadratic_approx!( - epi_cfg, - container, C, names, time_steps, - p1_expr, p1_bounds, meta_p1, - ) - zp2_expr = _add_quadratic_approx!( - epi_cfg, - container, C, names, time_steps, - p2_expr, p2_bounds, meta_p2, - ) - - # --- Create z variable and two-sided HybS bounds --- - z_var = add_variable_container!( - container, - BilinearProductVariable, - C, - names, - time_steps; - meta, - ) - hybrid_cons = add_constraints_container!( - container, - HybSBoundConstraint, - C, - names, - 1:2, - time_steps; - sparse = true, - meta, - ) - result_expr = add_expression_container!( - container, - BilinearProductExpression, - C, - names, - time_steps; - 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 - - # 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) - - z = - z_var[name, t] = JuMP.@variable( - jump_model, - base_name = "HybSProduct_$(C)_{$(name), $(t)}", - lower_bound = z_lo, - upper_bound = z_hi, - ) - - zx = xsq[name, t] - zy = ysq[name, t] - zp1 = zp1_expr[name, t] - zp2 = zp2_expr[name, t] - - # Bin2 lower bound: z ≥ ½(z_p1 − z_x − z_y) - hybrid_cons[(name, 1, t)] = JuMP.@constraint( - jump_model, - z >= 0.5 * (zp1 - zx - zy), - ) - # Bin3 upper bound: z ≤ ½(z_x + z_y − z_p2) - hybrid_cons[(name, 2, t)] = JuMP.@constraint( - jump_model, - z <= 0.5 * (zx + zy - zp2), - ) - - result_expr[name, t] = JuMP.AffExpr(0.0, z => 1.0) - end - - # --- Standard McCormick envelope cuts on the product variable --- - if config.add_mccormick - _add_mccormick_envelope!( - container, C, names, time_steps, - x_var, y_var, z_var, - x_bounds, y_bounds, meta, - ) - end - - return result_expr -end diff --git a/src/bilinear_approximations/mccormick.jl b/src/bilinear_approximations/mccormick.jl deleted file mode 100644 index 8ee2e3d2..00000000 --- a/src/bilinear_approximations/mccormick.jl +++ /dev/null @@ -1,404 +0,0 @@ -# McCormick envelope for bilinear products z = x·y. -# Adds 4 linear inequalities that bound z given variable bounds on x and y. - -"Standard McCormick envelope constraints bounding the bilinear product z = x·y." -struct McCormickConstraint <: ConstraintType end - -"Reformulated McCormick constraints on Bin2 separable variables." -struct ReformulatedMcCormickConstraint <: ConstraintType end - -""" - _mc_setindex!(cons, index, n, constraint) - -Helper function for setting constraints by-index in a McCormick constraint container. - -Supports 2- and 3-length tuples. -""" -@inline function _mc_setindex!(cons, index::Tuple{A, B}, n::Int, constraint) where {A, B} - cons[index[1], n, index[2]] = constraint -end - -@inline function _mc_setindex!( - cons, - index::Tuple{A, B, C}, - n::Int, - constraint, -) where {A, B, C} - cons[index[1], index[2], n, index[3]] = constraint -end - -function _add_mccormick_envelope!( - jump_model::JuMP.Model, - cons, - index, - x::JuMP.AbstractJuMPScalar, - y::JuMP.AbstractJuMPScalar, - z::JuMP.AbstractJuMPScalar, - x_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64; - lower_bounds::Bool = true, -) - if lower_bounds - _mc_setindex!( - cons, - index, - 1, - JuMP.@constraint( - jump_model, - z >= x_min * y + x * y_min - x_min * y_min, - ) - ) - _mc_setindex!( - cons, - index, - 2, - JuMP.@constraint( - jump_model, - z >= x_max * y + x * y_max - x_max * y_max, - ) - ) - end - _mc_setindex!( - cons, - index, - 3, - JuMP.@constraint( - jump_model, - z <= x_max * y + x * y_min - x_max * y_min, - ) - ) - _mc_setindex!( - cons, - index, - 4, - JuMP.@constraint( - jump_model, - z <= x_min * y + x * y_max - x_min * y_max, - ) - ) -end - -""" - _add_mccormick_envelope!(container, C, names, time_steps, x_var, y_var, z_var, x_min, x_max, y_min, y_max, meta) - -Add McCormick envelope constraints for the bilinear product z ≈ x·y. - -For each (name, t), adds 4 linear inequalities: -``` -z ≥ x_min·y + x·y_min − x_min·y_min -z ≥ x_max·y + x·y_max − x_max·y_max -z ≤ x_max·y + x·y_min − x_max·y_min -z ≤ x_min·y + x·y_max − x_min·y_max -``` - -# 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_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 -- `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_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64, - meta::String; - lower_bounds::Bool = true, -) 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!( - container, - McCormickConstraint, - C, - names, - 1:4, - time_steps; - sparse = true, - meta, - ) - - for name in names, t in time_steps - _add_mccormick_envelope!( - jump_model, mc_cons, (name, t), - x_var[name, t], y_var[name, t], z_var[name, t], - x_min, x_max, y_min, y_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, - x_min::Float64, - x_max::Float64, - meta::String; - lower_bounds::Bool = true, -) where {C <: IS.InfrastructureSystemsComponent} - _add_mccormick_envelope!( - container, C, names, time_steps, - x_var, x_var, z_var, - x_min, x_max, x_min, x_max, - meta; lower_bounds, - ) - 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, - index, - x::JuMP.VariableRef, - z::JuMP.VariableRef, - x_min::Float64, - x_max::Float64; - lower_bounds::Bool = true, -) - _add_mccormick_envelope!( - jump_model, cons, index, - x, x, z, - x_min, x_max, x_min, x_max; - lower_bounds, - ) -end - -# Lower McCormick bounds on (z_p1 − z_x − z_y) for the Bin2 reformulation. -function _add_reformulated_lower_mccormick!( - jump_model::JuMP.Model, - cons, - index, - x::JuMP.AbstractJuMPScalar, - y::JuMP.AbstractJuMPScalar, - zp1::JuMP.AbstractJuMPScalar, - zx::JuMP.AbstractJuMPScalar, - zy::JuMP.AbstractJuMPScalar, - x_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64, -) - _mc_setindex!( - cons, - index, - 1, - JuMP.@constraint( - jump_model, - zp1 - zx - zy >= 2.0 * (x_min * y + x * y_min - x_min * y_min), - ) - ) - _mc_setindex!( - cons, - index, - 2, - JuMP.@constraint( - jump_model, - zp1 - zx - zy >= 2.0 * (x_max * y + x * y_max - x_max * y_max), - ) - ) -end - -function _add_reformulated_mccormick_bin2!( - jump_model::JuMP.Model, - cons, - index, - x::JuMP.AbstractJuMPScalar, - y::JuMP.AbstractJuMPScalar, - zp1::JuMP.AbstractJuMPScalar, - zx::JuMP.AbstractJuMPScalar, - zy::JuMP.AbstractJuMPScalar, - x_min::Float64, - x_max::Float64, - y_min::Float64, - y_max::Float64, -) - _add_reformulated_lower_mccormick!( - jump_model, cons, index, x, y, zp1, zx, zy, x_min, x_max, y_min, y_max, - ) - # Upper bounds also on (z_p1 − z_x − z_y) since Bin2 has no z_p2 - _mc_setindex!( - cons, - index, - 3, - JuMP.@constraint( - jump_model, - zp1 - zx - zy <= 2.0 * (x_max * y + x * y_min - x_max * y_min), - ) - ) - _mc_setindex!( - cons, - index, - 4, - JuMP.@constraint( - jump_model, - zp1 - zx - zy <= 2.0 * (x_min * y + x * y_max - x_min * y_max), - ) - ) -end - -""" - _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, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_var, - y_var, - psq, - xsq, - ysq, - x_bounds::Vector{MinMax}, - y_bounds::Vector{MinMax}, - meta::String, -) where {C <: IS.InfrastructureSystemsComponent} - jump_model = get_jump_model(container) - - mc_cons = add_constraints_container!( - container, - ReformulatedMcCormickConstraint, - 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_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], - xb.min, xb.max, yb.min, yb.max, - ) - end - - return -end diff --git a/src/bilinear_approximations/nmdt.jl b/src/bilinear_approximations/nmdt.jl deleted file mode 100644 index 247ca71f..00000000 --- a/src/bilinear_approximations/nmdt.jl +++ /dev/null @@ -1,267 +0,0 @@ -# DNMDT (Double Normalized Multiparametric Disaggregation Technique) bilinear approximation of x·y. -# Independently discretizes both x and y, forms four cross binary-continuous products, then -# combines two NMDT estimates with a convex weighting λ (default 0.5). Reduces to the NMDT -# formulation when applied to x·x (quadratic case). -# Reference: Teles, Castro, Matos (2013), Multiparametric disaggregation technique for global -# optimization of polynomial programming problems. - -""" -Config for double-NMDT bilinear approximation (discretizes both x and y). - -# Fields -- `depth::Int`: number of binary discretization levels L for both x and y -""" -struct DNMDTBilinearConfig <: BilinearApproxConfig - depth::Int -end - -""" -Config for single-NMDT bilinear approximation (discretizes x only). - -# Fields -- `depth::Int`: number of binary discretization levels L for x -""" -struct NMDTBilinearConfig <: BilinearApproxConfig - depth::Int -end - -# --- DNMDT bilinear approximation --- - -""" - _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. - -Constructs all four cross binary-continuous products (β_x·yh, β_y·δx, β_y·xh, β_x·δy) -then delegates to the core DNMDT assembler. Stores results in a `BilinearProductExpression` -container. - -# Arguments -- `config::DNMDTBilinearConfig`: configuration for the DNMDT bilinear approximation -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `names::Vector{String}`: component names -- `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!( - config::DNMDTBilinearConfig, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_disc::NMDTDiscretization, - y_disc::NMDTDiscretization, - x_bounds::Vector{MinMax}, - y_bounds::Vector{MinMax}, - meta::String, -) where {C <: IS.InfrastructureSystemsComponent} - bx_yh_expr = _binary_continuous_product!( - container, C, names, time_steps, - x_disc, y_disc.norm_expr, 0.0, 1.0, - config.depth, meta * "_bx_yh", - ) - by_dx_expr = _binary_continuous_product!( - container, C, names, time_steps, - y_disc, x_disc.delta_var, 0.0, 2.0^(-config.depth), - config.depth, meta * "_by_dx", - ) - by_xh_expr = _binary_continuous_product!( - container, C, names, time_steps, - y_disc, x_disc.norm_expr, 0.0, 1.0, - config.depth, meta * "_by_xh", - ) - bx_dy_expr = _binary_continuous_product!( - container, C, names, time_steps, - x_disc, y_disc.delta_var, 0.0, 2.0^(-config.depth), - config.depth, meta * "_bx_dy", - ) - - 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_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_bounds, y_bounds, meta) - -Approximate x·y using the DNMDT method from raw variable inputs. - -Discretizes both x and y independently via `_discretize!` then delegates to the -pre-discretized overload. - -# Arguments -- `config::DNMDTBilinearConfig`: configuration for the DNMDT bilinear approximation -- `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) -- `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!( - config::DNMDTBilinearConfig, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_var, - y_var, - x_bounds::Vector{MinMax}, - y_bounds::Vector{MinMax}, - meta::String, -) where {C <: IS.InfrastructureSystemsComponent} - x_disc = _discretize!( - container, - C, - names, - time_steps, - x_var, - x_bounds, - config.depth, - meta * "_x", - ) - y_disc = _discretize!( - container, - C, - names, - time_steps, - y_var, - y_bounds, - config.depth, - meta * "_y", - ) - return _add_bilinear_approx!( - config, - container, - C, - names, - time_steps, - x_disc, - y_disc, - x_bounds, - y_bounds, - meta, - ) -end - -# --- NMDT bilinear approximation --- - -""" - _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. - -Discretizes only x (using `x_disc`) while y is already normalized to yh ∈ [0,1]. -Computes binary-continuous product β_x·yh and residual product δ_x·yh, then assembles -x·y via `_assemble_product!`. Stores results in a `BilinearProductExpression` container. - -# Arguments -- `config::NMDTBilinearConfig`: configuration for the NMDT bilinear approximation -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `names::Vector{String}`: component names -- `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) -- `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!( - config::NMDTBilinearConfig, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_disc::NMDTDiscretization, - yh_expr, - x_bounds::Vector{MinMax}, - y_bounds::Vector{MinMax}, - meta::String; -) where {C <: IS.InfrastructureSystemsComponent} - bx_y_expr = _binary_continuous_product!( - container, C, names, time_steps, - x_disc, yh_expr, 0.0, 1.0, - config.depth, meta, - ) - dz = _residual_product!( - container, C, names, time_steps, - x_disc, yh_expr, 1.0, config.depth, meta; - ) - - return _assemble_product!( - container, C, names, time_steps, - [bx_y_expr], dz, - 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_bounds, y_bounds, meta) - -Approximate x·y using the NMDT method from raw variable inputs. - -Discretizes x via `_discretize!` and normalizes y via `_normed_variable!`, then -delegates to the pre-discretized overload. - -# Arguments -- `config::NMDTBilinearConfig`: configuration for the NMDT bilinear approximation -- `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) -- `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!( - config::NMDTBilinearConfig, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_var, - y_var, - x_bounds::Vector{MinMax}, - y_bounds::Vector{MinMax}, - meta::String, -) where {C <: IS.InfrastructureSystemsComponent} - x_disc = _discretize!( - container, - C, - names, - time_steps, - x_var, - x_bounds, - config.depth, - meta * "_x", - ) - yh_expr = - _normed_variable!(container, C, names, time_steps, y_var, y_bounds, meta * "_y") - return _add_bilinear_approx!( - config, - container, - C, - names, - time_steps, - x_disc, - yh_expr, - x_bounds, - y_bounds, - meta, - ) -end diff --git a/src/bilinear_approximations/no_approx.jl b/src/bilinear_approximations/no_approx.jl deleted file mode 100644 index a6fd81fe..00000000 --- a/src/bilinear_approximations/no_approx.jl +++ /dev/null @@ -1,100 +0,0 @@ -# No-op bilinear approximation: returns exact x·y as a QuadExpr. -# For NLP-capable solvers or testing purposes. - -"No-op bilinear config: returns exact x·y as a QuadExpr." -struct NoBilinearApproxConfig <: BilinearApproxConfig end - -""" - _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. - -# Arguments -- `::NoBilinearApproxConfig`: no-op configuration (no fields) -- `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) -- `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!( - ::NoBilinearApproxConfig, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_var, - y_var, - x_bounds::Vector{MinMax}, - y_bounds::Vector{MinMax}, - meta::String, -) where {C <: IS.InfrastructureSystemsComponent} - result_expr = add_expression_container!( - container, - BilinearProductExpression, - C, - names, - time_steps; - meta, - expr_type = JuMP.QuadExpr, - ) - for name in names, t in time_steps - result_expr[name, t] = x_var[name, t] * y_var[name, t] - end - return result_expr -end - -""" - _add_bilinear_approx!(::NoBilinearApproxConfig, container, C, names, time_steps, xsq, ysq, x_var, y_var, x_bounds, y_bounds, meta) - -Precomputed-form no-op bilinear approximation: returns exact x·y as a QuadExpr. -`xsq` and `ysq` are accepted for signature parity with the precomputed-form -dispatches of `Bin2Config` and `HybSConfig` (so a caller can swap configs -without changing the call site) but are unused here. - -# Arguments -- `::NoBilinearApproxConfig`: no-op configuration (no fields) -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `names::Vector{String}`: component names -- `time_steps::UnitRange{Int}`: time periods -- `xsq`: precomputed x² container (ignored) -- `ysq`: precomputed y² container (ignored) -- `x_var`: container of x variables indexed by (name, t) -- `y_var`: container of y variables indexed by (name, t) -- `x_bounds::Vector{MinMax}`: per-component bounds on x -- `y_bounds::Vector{MinMax}`: per-component bounds on y -- `meta::String`: variable type identifier for the approximation -""" -function _add_bilinear_approx!( - ::NoBilinearApproxConfig, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - xsq, - ysq, - x_var, - y_var, - x_bounds::Vector{MinMax}, - y_bounds::Vector{MinMax}, - meta::String, -) where {C <: IS.InfrastructureSystemsComponent} - result_expr = add_expression_container!( - container, - BilinearProductExpression, - C, - names, - time_steps; - meta, - expr_type = JuMP.QuadExpr, - ) - for name in names, t in time_steps - result_expr[name, t] = x_var[name, t] * y_var[name, t] - end - return result_expr -end diff --git a/src/quadratic_approximations/common.jl b/src/quadratic_approximations/common.jl deleted file mode 100644 index be37adfa..00000000 --- a/src/quadratic_approximations/common.jl +++ /dev/null @@ -1,54 +0,0 @@ -"Expression container for the normalized variable xh = (x − x_min) / (x_max − x_min) ∈ [0,1]." -struct NormedVariableExpression <: ExpressionType end - -"Expression container for quadratic (x²) approximation results." -struct QuadraticExpression <: ExpressionType end - -# --- Quadratic approximation config hierarchy --- - -"Abstract supertype for quadratic approximation method configurations." -abstract type QuadraticApproxConfig end - -""" - _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]. - -Stores results in a `NormedVariableExpression` expression container. - -# 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 variables indexed by (name, t) -- `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!( - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_var, - bounds::Vector{MinMax}, - meta::String, -) where {C <: IS.InfrastructureSystemsComponent} - result_expr = add_expression_container!( - container, - NormedVariableExpression, - C, - names, - time_steps; - meta, - ) - - 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, -b.min / lx) - end - return result_expr -end diff --git a/src/quadratic_approximations/epigraph.jl b/src/quadratic_approximations/epigraph.jl deleted file mode 100644 index 2b180b6c..00000000 --- a/src/quadratic_approximations/epigraph.jl +++ /dev/null @@ -1,193 +0,0 @@ -# Epigraph (Q^{L1}) LP-only lower bound for x² using tangent-line cuts. -# Pure LP — zero binary variables. Creates a variable z ≥ x² (approximately) -# bounded from below by supporting hyperplanes of the parabola. -# Reference: Beach, Burlacu, Hager, Hildebrand (2024), Q^{L1} relaxation. - -"Expression container for epigraph quadratic approximation results." -struct EpigraphExpression <: ExpressionType end - -"Variable representing a lower-bounded approximation of x² in epigraph relaxation." -struct EpigraphVariable <: VariableType end -"Tangent-line lower-bound constraints in epigraph relaxation." -struct EpigraphTangentConstraint <: ConstraintType end -"Tangent-line lower-bound expression fL used in the epigraph formulation." -struct EpigraphTangentExpression <: ExpressionType end - -""" -Config for epigraph (Q^{L1}) LP-only lower-bound quadratic approximation. - -# Fields -- `depth::Int`: number of tangent-line breakpoints (2^depth + 1 tangent lines); pure LP, zero binary variables -""" -struct EpigraphQuadConfig <: QuadraticApproxConfig - depth::Int -end - -""" - _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). - -For each (name, t), creates a variable z and adds 2^depth + 1 tangent-line -constraints of the form `z ≥ 2·aₖ·x − aₖ²` at uniformly spaced breakpoints -aₖ = x_min + k·Δ/2^depth for k = 0,…,2^depth. Pure LP — zero binary variables. - -Stores affine expressions that lower-bound x² in an `EpigraphExpression` expression container. - -The maximum underestimation gap between the tangent envelope and x² is -Δ²·2^{−2·depth−2} where Δ = x_max − x_min. - -# Arguments -- `config::EpigraphQuadConfig`: configuration with `depth` field controlling the number of tangent-line breakpoints (2^depth + 1 tangent lines) -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `names::Vector{String}`: component names -- `time_steps::UnitRange{Int}`: time periods -- `x_var`: container of variables indexed by (name, t) -- `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!( - config::EpigraphQuadConfig, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_var, - bounds::Vector{MinMax}, - meta::String, -) where {C <: IS.InfrastructureSystemsComponent} - IS.@assert_op config.depth >= 1 - jump_model = get_jump_model(container) - g_levels = 0:(config.depth) - - z_var = add_variable_container!( - container, - EpigraphVariable, - C, - names, - time_steps; - meta, - ) - g_var = add_variable_container!( - container, - SawtoothAuxVariable, - C, - names, - g_levels, - time_steps; - meta, - ) - lp_cons = add_constraints_container!( - container, - SawtoothLPConstraint, - C, - names, - 1:2, - time_steps; - meta, - ) - link_cons = add_constraints_container!( - container, - SawtoothLinkingConstraint, - C, - names, - time_steps; - meta, - ) - fL_expr = add_expression_container!( - container, - EpigraphTangentExpression, - C, - names, - time_steps; - meta, - ) - tangent_cons = add_constraints_container!( - container, - EpigraphTangentConstraint, - C, - names, - 1:(config.depth + 2), - time_steps; - sparse = true, - meta, - ) - result_expr = add_expression_container!( - container, - EpigraphExpression, - C, - names, - time_steps; - meta, - ) - - 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] - for j in g_levels - g_var[name, j, t] = JuMP.@variable( - jump_model, - base_name = "SawtoothAux_$(C)_{$(name), $(j), $(t)}", - lower_bound = 0.0, - upper_bound = 1.0, - ) - end - g0 = g_var[name, 0, t] - - # Linking constraint: g_0 = (x - x_min) / Δ - link_cons[name, t] = JuMP.@constraint( - jump_model, - g0 == (x - b.min) / delta, - ) - - # T^L constraints for j = 1,...,L - for j in 1:(config.depth) - g_prev = g_var[name, j - 1, t] - g_curr = g_var[name, j, t] - - # g_j ≤ 2 g_{j-1} - lp_cons[name, 1, t] = JuMP.@constraint(jump_model, g_curr <= 2.0 * g_prev) - # g_j ≤ 2(1 - g_{j-1}) - lp_cons[name, 2, t] = - JuMP.@constraint(jump_model, g_curr <= 2.0 * (1.0 - g_prev)) - end - - # Create the epigraph variable (bounded from below by tangent cuts) - z = - z_var[name, t] = JuMP.@variable( - jump_model, - base_name = "EpigraphVar_$(C)_{$(name), $(t)}", - lower_bound = 0.0, - upper_bound = z_ub, - ) - - fL = fL_expr[name, t] = JuMP.AffExpr(0.0) - for j in 1:(config.depth) - add_proportional_to_jump_expression!( - fL, - g_var[name, j, t], - delta * delta * 2.0^(-2j), - ) - tangent_cons[(name, j + 1, t)] = JuMP.@constraint( - jump_model, - z >= - 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 * b.min - 1.0 + 2.0 * delta * g0 - ) - - result_expr[name, t] = JuMP.AffExpr(0.0, z => 1.0) - end - - return result_expr -end diff --git a/src/quadratic_approximations/manual_sos2.jl b/src/quadratic_approximations/manual_sos2.jl deleted file mode 100644 index 0275424a..00000000 --- a/src/quadratic_approximations/manual_sos2.jl +++ /dev/null @@ -1,235 +0,0 @@ -# SOS2-based piecewise linear approximation of x² for use in constraints. -# Uses manually-implemented SOS2 adjacency via binary variables and linear constraints. - -"Binary segment-selection variables (z) for manual SOS2 quadratic approximation." -struct ManualSOS2BinaryVariable <: SparseVariableType end -"Ensures exactly one segment is active (∑zⱼ = 1) in manual SOS2 quadratic approximation." -struct ManualSOS2SegmentSelectionConstraint <: ConstraintType end -"Expression for the segment selection sum Σ z_j in manual SOS2 quadratic approximation." -struct ManualSOS2SegmentSelectionExpression <: ExpressionType end -"Links active segment to lambda variables." -struct ManualSOS2AdjacencyConstraint <: ConstraintType end - -""" -Config for manual binary-variable SOS2 quadratic approximation. - -# Fields -- `depth::Int`: number of PWL segments (breakpoints = depth + 1) -- `pwmcc_segments::Int`: number of piecewise McCormick cut partitions; 0 to disable (default 4) -""" -struct ManualSOS2QuadConfig <: QuadraticApproxConfig - depth::Int - pwmcc_segments::Int -end -ManualSOS2QuadConfig(depth::Int) = ManualSOS2QuadConfig(depth, 4) - -""" - _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. - -Creates lambda (λ) variables representing convex combination weights over breakpoints, -adds linking, normalization, and manual adjacency constraints using binary variables, -and stores affine expressions approximating x² in a `QuadraticExpression` -expression container. - -# Arguments -- `config::ManualSOS2QuadConfig`: configuration with `depth` (number of PWL segments) and `pwmcc_segments` (PWMCC cut partitions; 0 to disable, default 4) -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `names::Vector{String}`: component names -- `time_steps::UnitRange{Int}`: time periods -- `x_var`: container of variables indexed by (name, t) -- `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!( - config::ManualSOS2QuadConfig, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_var, - bounds::Vector{MinMax}, - meta::String, -) where {C <: IS.InfrastructureSystemsComponent} - x_bkpts, x_sq_bkpts = - _get_breakpoints_for_pwl_function( - 0.0, - 1.0, - _square; - num_segments = config.depth, - ) - n_points = config.depth + 1 - n_bins = n_points - 1 - jump_model = get_jump_model(container) - - # Create all containers upfront - 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, - C, - names, - time_steps; - meta, - ) - link_expr = add_expression_container!( - container, - SOS2LinkingExpression, - C, - names, - time_steps; - meta, - ) - norm_cons = add_constraints_container!( - container, - SOS2NormConstraint, - C, - names, - time_steps; - meta, - ) - norm_expr = add_expression_container!( - container, - SOS2NormExpression, - C, - names, - time_steps; - meta, - ) - seg_cons = add_constraints_container!( - container, - ManualSOS2SegmentSelectionConstraint, - C, - names, - time_steps; - meta, - ) - seg_expr = add_expression_container!( - container, - ManualSOS2SegmentSelectionExpression, - C, - names, - time_steps; - meta, - ) - adj_cons = add_constraints_container!( - container, - ManualSOS2AdjacencyConstraint, - C, - names, - 1:n_points, - time_steps; - meta, - ) - result_expr = add_expression_container!( - container, - QuadraticExpression, - C, - names, - time_steps; - meta, - ) - - 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 - x = x_var[name, t] - - # Create lambda variables: λ_i ∈ [0, 1] - lambda = Vector{JuMP.VariableRef}(undef, n_points) - for i in 1:n_points - lambda[i] = - lambda_container[name, i, t] = JuMP.@variable( - jump_model, - base_name = "QuadraticVariable_$(C)_{$(name), pwl_$(i), $(t)}", - lower_bound = 0.0, - upper_bound = 1.0, - ) - end - - # x = Σ λ_i * x_i - link = link_expr[name, t] = JuMP.AffExpr(0.0) - 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 - b.min) / lx == link) - - # Σ λ_i = 1 - norm = norm_expr[name, t] = JuMP.AffExpr(0.0) - for l in lambda - add_proportional_to_jump_expression!(norm, l, 1.0) - end - norm_cons[name, t] = JuMP.@constraint(jump_model, norm == 1.0) - - # Create binary segment-selection variables z_j - z_vars = Vector{JuMP.VariableRef}(undef, n_bins) - for j in 1:n_bins - z_vars[j] = - z_container[name, j, t] = JuMP.@variable( - jump_model, - base_name = "ManualSOS2Binary_$(C)_{$(name), $(j), $(t)}", - binary = true, - ) - end - - # Σ z_j = 1 (segment selection) - seg = seg_expr[name, t] = JuMP.AffExpr(0.0) - for z in z_vars - add_proportional_to_jump_expression!(seg, z, 1.0) - end - seg_cons[name, t] = JuMP.@constraint(jump_model, seg == 1) - - # Adjacency constraints: λ_i ≤ z_{i-1} + z_i (with boundary cases) - # λ_1 ≤ z_1 - adj_cons[name, 1, t] = JuMP.@constraint(jump_model, lambda[1] <= z_vars[1]) - # λ_i ≤ z_{i-1} + z_i for i = 2..n-1 - for i in 2:(n_points - 1) - adj_cons[name, i + 1, t] = - JuMP.@constraint(jump_model, lambda[i] <= z_vars[i - 1] + z_vars[i]) - end - # λ_n ≤ z_{n-1} - adj_cons[name, n_points, t] = - JuMP.@constraint(jump_model, lambda[n_points] <= z_vars[n_bins]) - - # Build x̂² = Σ λ_i * x_i² as an affine expression - x_hat_sq = JuMP.AffExpr(0.0) - for i in 1:n_points - add_proportional_to_jump_expression!(x_hat_sq, lambda[i], x_sq_bkpts[i]) - 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 * 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, bounds, - config.pwmcc_segments, meta * "_pwmcc", - ) - end - - return result_expr -end diff --git a/src/quadratic_approximations/nmdt.jl b/src/quadratic_approximations/nmdt.jl deleted file mode 100644 index 343dbaec..00000000 --- a/src/quadratic_approximations/nmdt.jl +++ /dev/null @@ -1,229 +0,0 @@ -# NMDT (Normalized Multiparametric Disaggregation Technique) quadratic approximation of x². -# Normalizes x to [0,1], discretizes using L binary variables β₁,…,β_L plus a -# residual δ ∈ [0, 2^{−L}], then replaces each binary-continuous product β_i·xh -# with a McCormick-linearized auxiliary variable. Assembles the result via the -# separable identity x² = (lx·xh + x_min)². Optionally tightens with an epigraph -# lower bound on xh². -# NMDT Reference: Teles, Castro, Matos (2013), Multiparametric disaggregation -# technique for global optimization of polynomial programming problems. - -""" -Config for double-NMDT quadratic approximation. - -# Fields -- `depth::Int`: number of binary discretization levels L -- `epigraph_depth::Int`: LP tightening depth via epigraph Q^{L1} lower bound; 0 to disable (default 3×depth) -""" -struct DNMDTQuadConfig <: QuadraticApproxConfig - depth::Int - epigraph_depth::Int -end -DNMDTQuadConfig(depth::Int) = DNMDTQuadConfig(depth, 3 * depth) - -""" -Config for single-NMDT quadratic approximation. - -# Fields -- `depth::Int`: number of binary discretization levels L -- `epigraph_depth::Int`: LP tightening depth via epigraph Q^{L1} lower bound; 0 to disable (default 3×depth) -""" -struct NMDTQuadConfig <: QuadraticApproxConfig - depth::Int - epigraph_depth::Int -end -NMDTQuadConfig(depth::Int) = NMDTQuadConfig(depth, 3 * depth) - -""" - _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. - -Constructs two binary-continuous products (β·xh and β·δ) and delegates to the core -DNMDT assembler, storing results in a `QuadraticExpression` container. Optionally -tightens lower bounds with an epigraph relaxation via `_tighten_lower_bounds!`. - -# Arguments -- `config::DNMDTQuadConfig`: configuration with `depth` (binary discretization levels) and `epigraph_depth` (LP tightening depth; 0 to disable, default 3×depth) -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `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!( - config::DNMDTQuadConfig, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_disc::NMDTDiscretization, - bounds::Vector{MinMax}, - meta::String, -) where {C <: IS.InfrastructureSystemsComponent} - tighten = config.epigraph_depth > 0 - bx_xh_expr = _binary_continuous_product!( - container, C, names, time_steps, - x_disc, x_disc.norm_expr, 0.0, 1.0, - config.depth, meta * "_bx_xh"; tighten, - ) - bx_dx_expr = _binary_continuous_product!( - container, C, names, time_steps, - x_disc, x_disc.delta_var, 0.0, 2.0^(-config.depth), - config.depth, meta * "_bx_dx"; tighten, - ) - - 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, bounds, bounds, - config.depth, meta; tighten, - result_type = QuadraticExpression, - ) - - if config.epigraph_depth > 0 - _tighten_lower_bounds!( - container, C, names, time_steps, - result_expr, x_disc, config.epigraph_depth, meta, - ) - end - - return result_expr -end - -""" - _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. - -Discretizes x via `_discretize!` then delegates to the `NMDTDiscretization` overload. -Stores results in a `QuadraticExpression` container. - -# Arguments -- `config::DNMDTQuadConfig`: configuration with `depth` (binary discretization levels) and `epigraph_depth` (LP tightening depth; 0 to disable, default 3×depth) -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `names::Vector{String}`: component names -- `time_steps::UnitRange{Int}`: time periods -- `x_var`: container of variables indexed by (name, t) -- `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!( - config::DNMDTQuadConfig, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_var, - bounds::Vector{MinMax}, - meta::String, -) where {C <: IS.InfrastructureSystemsComponent} - x_disc = _discretize!( - container, C, names, time_steps, - x_var, bounds, config.depth, meta, - ) - - return _add_quadratic_approx!( - config, container, C, names, time_steps, - x_disc, bounds, meta, - ) -end - -""" - _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. - -Computes the binary-continuous product β·xh and residual product δ·xh, then -assembles x² via `_assemble_product!`. Stores results in a `QuadraticExpression` -container. Optionally tightens lower bounds with an epigraph relaxation. - -# Arguments -- `config::NMDTQuadConfig`: configuration with `depth` (binary discretization levels) and `epigraph_depth` (LP tightening depth; 0 to disable, default 3×depth) -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `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!( - config::NMDTQuadConfig, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_disc::NMDTDiscretization, - bounds::Vector{MinMax}, - meta::String, -) where {C <: IS.InfrastructureSystemsComponent} - tighten = config.epigraph_depth > 0 - bx_y_expr = _binary_continuous_product!( - container, C, names, time_steps, - x_disc, x_disc.norm_expr, 0.0, 1.0, - config.depth, meta; tighten, - ) - dz = _residual_product!( - container, C, names, time_steps, - x_disc, x_disc.norm_expr, 1.0, - config.depth, meta; tighten, - ) - - result_expr = _assemble_product!( - container, C, names, time_steps, - [bx_y_expr], dz, - x_disc, x_disc, bounds, bounds, - meta; result_type = QuadraticExpression, - ) - - if config.epigraph_depth > 0 - _tighten_lower_bounds!( - container, C, names, time_steps, - result_expr, x_disc, config.epigraph_depth, meta, - ) - end - - return result_expr -end - -""" - _add_quadratic_approx!(config::NMDTQuadConfig, container, C, names, time_steps, x_var, bounds, meta) - -Approximate x² using the NMDT method from raw variable inputs. - -Discretizes x via `_discretize!` then delegates to the `NMDTDiscretization` overload. -Stores results in a `QuadraticExpression` container. - -# Arguments -- `config::NMDTQuadConfig`: configuration with `depth` (binary discretization levels) and `epigraph_depth` (LP tightening depth; 0 to disable, default 3×depth) -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `names::Vector{String}`: component names -- `time_steps::UnitRange{Int}`: time periods -- `x_var`: container of variables indexed by (name, t) -- `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!( - config::NMDTQuadConfig, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_var, - bounds::Vector{MinMax}, - meta::String, -) where {C <: IS.InfrastructureSystemsComponent} - x_disc = _discretize!( - container, C, names, time_steps, - x_var, bounds, config.depth, meta, - ) - - return _add_quadratic_approx!( - config, container, C, names, time_steps, - x_disc, bounds, meta, - ) -end diff --git a/src/quadratic_approximations/nmdt_common.jl b/src/quadratic_approximations/nmdt_common.jl deleted file mode 100644 index af35a90c..00000000 --- a/src/quadratic_approximations/nmdt_common.jl +++ /dev/null @@ -1,554 +0,0 @@ -"Binary discretization variables β_i ∈ {0,1} in the NMDT decomposition of xh." -struct NMDTBinaryVariable <: VariableType end -"Residual variable δ ∈ [0, 2^{−L}] capturing the NMDT discretization error." -struct NMDTResidualVariable <: VariableType end -"McCormick linearization variables u_i ≈ β_i · y in NMDT binary-continuous products." -struct NMDTBinaryContinuousProductVariable <: VariableType end -"Variable z ≈ δ · y linearizing the residual-continuous product in NMDT." -struct NMDTResidualProductVariable <: VariableType end - -"Expression container for the NMDT binary discretization: Σ 2^{−i}·β_i + δ ≈ xh." -struct NMDTDiscretizationExpression <: ExpressionType end -"Expression container for the NMDT binary-continuous product: Σ 2^{−i}·u_i ≈ β·y." -struct NMDTBinaryContinuousProductExpression <: ExpressionType end -"Expression container for the final NMDT quadratic approximation result." -struct NMDTResultExpression <: ExpressionType end - -"Constraint enforcing xh = Σ 2^{−i}·β_i + δ in the NMDT discretization." -struct NMDTEDiscretizationConstraint <: ConstraintType end -"McCormick envelope constraints for binary-continuous products u_i ≈ β_i·y in NMDT." -struct NMDTBinaryContinuousProductConstraint <: ConstraintType end -"Epigraph lower-bound tightening constraint on the NMDT quadratic result." -struct NMDTTightenConstraint <: ConstraintType end - -""" -Stores the result of discretizing a normalized variable for use in NMDT products. - -Fields: -- `norm_expr`: affine expression for xh = (x − x_min)/(x_max − x_min) ∈ [0,1] -- `beta_var`: binary variables β_i ∈ {0,1} indexed by (name, i, t) -- `delta_var`: residual variables δ ∈ [0, 2^{−depth}] indexed by (name, t) -""" -struct NMDTDiscretization{NE, BV, DV} - norm_expr::NE - beta_var::BV - delta_var::DV -end - -""" - _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. - -Creates L binary variables β₁,…,β_L and one residual δ ∈ [0, 2^{−L}] such that -xh = Σᵢ 2^{−i}·β_i + δ. Enforces this via a `NMDTEDiscretizationConstraint` and -returns an `NMDTDiscretization` struct holding all components. - -# 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 variables indexed by (name, t) -- `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 -""" -function _discretize!( - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_var, - bounds::Vector{MinMax}, - depth::Int, - meta::String; -) where {C <: IS.InfrastructureSystemsComponent} - IS.@assert_op depth >= 1 - jump_model = get_jump_model(container) - - beta_var = add_variable_container!( - container, - NMDTBinaryVariable, - C, - names, - 1:depth, - time_steps; - meta, - ) - delta_var = add_variable_container!( - container, - NMDTResidualVariable, - C, - names, - time_steps; - meta, - ) - disc_expr = add_expression_container!( - container, - NMDTDiscretizationExpression, - C, - names, - time_steps; - meta, - ) - disc_cons = add_constraints_container!( - container, - NMDTEDiscretizationConstraint, - C, - names, - time_steps; - meta, - ) - - xh_expr = _normed_variable!( - container, C, names, time_steps, - x_var, bounds, meta, - ) - - for name in names, t in time_steps - disc = disc_expr[name, t] = JuMP.AffExpr(0.0) - for i in 1:depth - beta = - beta_var[name, i, t] = JuMP.@variable( - jump_model, - base_name = "NMDTBinary_$(C)_{$(name), $(i), $(t)}", - binary = true - ) - add_proportional_to_jump_expression!(disc, beta, 2.0^(-i)) - end - delta = - delta_var[name, t] = JuMP.@variable( - jump_model, - base_name = "NMDTResidual_$(C)_{$(name), $(t)}", - lower_bound = 0.0, - upper_bound = 2.0^(-depth) - ) - add_proportional_to_jump_expression!(disc, delta, 1.0) - disc_cons[name, t] = JuMP.@constraint( - jump_model, - xh_expr[name, t] == disc - ) - end - - return NMDTDiscretization(xh_expr, beta_var, delta_var) -end - -""" - _binary_continuous_product!(container, C, names, time_steps, bin_disc, cont_var, cont_min, cont_max, depth, meta; tighten) - -Linearize each binary-continuous product β_i·y using McCormick envelopes. - -For each depth level i, creates a variable u_i ≈ β_i·y with bounds [cont_min, cont_max] -and adds 4 McCormick constraints via `_add_mccormick_envelope!`. Assembles the weighted sum -Σᵢ 2^{−i}·u_i into a `NMDTBinaryContinuousProductExpression`. - -# Arguments -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `names::Vector{String}`: component names -- `time_steps::UnitRange{Int}`: time periods -- `bin_disc`: `NMDTDiscretization` providing β_i variables and depth -- `cont_var`: container of continuous variables y indexed by (name, t) -- `cont_min::Float64`: lower bound of y -- `cont_max::Float64`: upper bound of y -- `depth::Int`: number of binary discretization levels -- `meta::String`: identifier encoding the original variable type being approximated -- `tighten::Bool`: if true, omit McCormick lower bounds (for use when a tighter bound is applied elsewhere) -""" -function _binary_continuous_product!( - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - bin_disc, - cont_var, - cont_min::Float64, - cont_max::Float64, - depth::Int, - meta::String; - tighten::Bool = false, -) where {C <: IS.InfrastructureSystemsComponent} - jump_model = get_jump_model(container) - - u_var = add_variable_container!( - container, - NMDTBinaryContinuousProductVariable, - C, - names, - 1:depth, - time_steps; - meta, - ) - u_cons = add_constraints_container!( - container, - NMDTBinaryContinuousProductConstraint, - C, - names, - 1:depth, - 1:4, - time_steps; - meta, - ) - result_expr = add_expression_container!( - container, - NMDTBinaryContinuousProductExpression, - C, - names, - time_steps; - meta, - ) - - for name in names, t in time_steps - result = result_expr[name, t] = JuMP.AffExpr(0.0) - for i in 1:depth - u_i = - u_var[name, i, t] = JuMP.@variable( - jump_model, - base_name = "NMDTBinContProd_$(C)_{$(name), $(i), $(t)}", - lower_bound = cont_min, - upper_bound = cont_max - ) - _add_mccormick_envelope!( - jump_model, u_cons, (name, i, t), - cont_var[name, t], bin_disc.beta_var[name, i, t], u_i, - cont_min, cont_max, 0.0, 1.0; - lower_bounds = !tighten, - ) - add_proportional_to_jump_expression!(result, u_i, 2.0^(-i)) - end - end - - return result_expr -end - -""" - _tighten_lower_bounds!(container, C, names, time_steps, result_expr, x_disc, epigraph_depth, meta) - -Add epigraph lower-bound constraints to tighten an NMDT quadratic approximation. - -Computes an epigraph Q^{L1} lower bound on xh² -and adds a `NMDTTightenConstraint` enforcing `result_expr[name,t] ≥ epi_expr[name,t]` -for each (name, t). This improves the lower bound quality of NMDT without adding binaries. - -# Arguments -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `names::Vector{String}`: component names -- `time_steps::UnitRange{Int}`: time periods -- `result_expr`: expression container for the NMDT quadratic result to be tightened -- `x_disc`: `NMDTDiscretization` for x, providing `norm_expr` and `depth` -- `epigraph_depth::Int`: depth for the epigraph Q^{L1} lower-bound approximation -- `meta::String`: identifier encoding the original variable type being approximated -""" -function _tighten_lower_bounds!( - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - result_expr, - x_disc, - epigraph_depth::Int, - meta::String; -) where {C <: IS.InfrastructureSystemsComponent} - jump_model = get_jump_model(container) - - epi_expr = _add_quadratic_approx!( - EpigraphQuadConfig(epigraph_depth), - container, C, names, time_steps, - x_disc.norm_expr, fill(MinMax((min = 0.0, max = 1.0)), length(names)), - meta * "_epi", - ) - epi_cons = add_constraints_container!( - container, - NMDTTightenConstraint, - C, - names, - time_steps; - meta, - ) - for name in names, t in time_steps - epi_cons[name, t] = JuMP.@constraint( - jump_model, - result_expr[name, t] >= epi_expr[name, t], - ) - end -end - -""" - _residual_product!(container, C, names, time_steps, x_disc, y_var, y_max, meta; tighten) - -Linearize the residual-continuous product z ≈ δ·y using McCormick envelopes. - -Creates a variable z ∈ [0, 2^{−L}·y_max] for each (name, t) and bounds it with -McCormick constraints on (δ, y) where δ ∈ [0, 2^{−L}]. Stores results in a -`NMDTResidualProductVariable` container. - -# Arguments -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `names::Vector{String}`: component names -- `time_steps::UnitRange{Int}`: time periods -- `x_disc`: `NMDTDiscretization` for x, providing `delta_var` and `depth` -- `y_var`: container of continuous variables y indexed by (name, t) -- `y_max::Float64`: upper bound of y (lower bound assumed 0) -- `meta::String`: identifier encoding the original variable type being approximated -- `tighten::Bool`: if true, omit McCormick lower bounds (default: false) -""" -function _residual_product!( - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_disc, - y_var, - y_max::Float64, - depth::Int, - meta::String; - tighten::Bool = false, -) where {C <: IS.InfrastructureSystemsComponent} - x_max = 2.0^(-depth) - jump_model = get_jump_model(container) - - z_var = add_variable_container!( - container, - NMDTResidualProductVariable, - C, - names, - time_steps; - meta, - ) - - for name in names, t in time_steps - z_var[name, t] = JuMP.@variable( - jump_model, - base_name = "NMDTResidualProduct_$(C)_{$(name), $(t)}", - lower_bound = 0.0, - upper_bound = x_max * y_max, - ) - end - - _add_mccormick_envelope!( - container, C, names, time_steps, - x_disc.delta_var, y_var, z_var, - 0.0, x_max, 0.0, y_max, - meta; lower_bounds = !tighten, - ) - - return z_var -end - -""" - _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. - -Applies the affine rescaling: -``` -x·y = lx·ly·zh + lx·y_min·xh + ly·x_min·yh + x_min·y_min -``` -where `zh = Σ terms[name,t] + dz_var[name,t]` collects the binary-continuous and -residual product contributions, lx = x_max − x_min, ly = y_max − y_min. - -Stores results in an expression container of type `result_type`. - -# Arguments -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `names::Vector{String}`: component names -- `time_steps::UnitRange{Int}`: time periods -- `terms`: iterable of expression containers indexed by (name, t) for the binary-continuous products -- `dz_var`: variable container for the residual product δ·y -- `xh_norm`: normed expression container for x -- `yh_norm`: normed expression container for y -- `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`) -""" -function _assemble_product!( - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - terms, - dz_var, - xh_expr, - yh_expr, - x_bounds::Vector{MinMax}, - y_bounds::Vector{MinMax}, - meta::String; - result_type = NMDTResultExpression, -) where {C <: IS.InfrastructureSystemsComponent} - result_expr = add_expression_container!( - container, - result_type, - C, - names, - time_steps; - 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 - 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 - add_proportional_to_jump_expression!(zh, term[name, t], 1.0) - end - 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 * 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_bounds, y_bounds, meta; result_type) - -Convenience overload: extracts `norm_expr` from both discretizations and delegates to the core `_assemble_product!`. -""" -function _assemble_product!( - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - terms, - dz_var, - x_disc::NMDTDiscretization, - y_disc::NMDTDiscretization, - 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_bounds, y_bounds, - meta; result_type, - ) -end - -""" - _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!`. -""" -function _assemble_product!( - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - terms, - dz_var, - x_disc::NMDTDiscretization, - yh_expr, - 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_bounds, y_bounds, - meta; result_type, - ) -end - -""" - _assemble_dnmdt!(container, C, names, time_steps, bx_yh_expr, by_dx_expr, by_xh_expr, bx_dy_expr, x_disc, y_disc, meta; lambda, result_type) - -Core assembler for the DNMDT bilinear approximation of x·y from pre-computed cross products. - -Builds two NMDT product estimates from opposite discretization pairings and combines them: -- z₁ = assemble(bx·yh + by·δx + δx·δy, x_disc, y_disc) -- z₂ = assemble(by·xh + bx·δy + δx·δy, y_disc, x_disc) -- result = λ·z₁ + (1−λ)·z₂ - -The shared residual product δx·δy is computed internally. Stores results in an expression -container of type `result_type`. - -# Arguments -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `names::Vector{String}`: component names -- `time_steps::UnitRange{Int}`: time periods -- `bx_yh_expr`: expression for β_x·yh binary-continuous products -- `by_dx_expr`: expression for β_y·δx binary-continuous products -- `by_xh_expr`: expression for β_y·xh binary-continuous products -- `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`) -""" -function _assemble_dnmdt!( - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - bx_yh_expr, - by_dx_expr, - by_xh_expr, - bx_dy_expr, - x_disc::NMDTDiscretization, - y_disc::NMDTDiscretization, - x_bounds::Vector{MinMax}, - y_bounds::Vector{MinMax}, - depth::Int, - meta::String; - lambda::Float64 = DNMDT_LAMBDA, - result_type::Type = NMDTResultExpression, - tighten::Bool = false, -) where {C <: IS.InfrastructureSystemsComponent} - result_expr = add_expression_container!( - container, - result_type, - C, - names, - time_steps; - meta, - ) - - dz = _residual_product!( - container, C, names, time_steps, - x_disc, y_disc.delta_var, 2.0^(-depth), - depth, meta; tighten, - ) - z1_expr = _assemble_product!( - container, C, names, time_steps, - [bx_yh_expr, by_dx_expr], dz, - 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_bounds, x_bounds, - meta * "_nmdt2", - ) - - for name in names, t in time_steps - result = result_expr[name, t] = JuMP.AffExpr(0.0) - add_proportional_to_jump_expression!(result, z1_expr[name, t], lambda) - add_proportional_to_jump_expression!(result, z2_expr[name, t], 1.0 - lambda) - end - - return result_expr -end diff --git a/src/quadratic_approximations/no_approx.jl b/src/quadratic_approximations/no_approx.jl deleted file mode 100644 index e2d88498..00000000 --- a/src/quadratic_approximations/no_approx.jl +++ /dev/null @@ -1,45 +0,0 @@ -# No-op quadratic approximation: returns exact x² as a QuadExpr. -# For NLP-capable solvers or testing purposes. - -"No-op config: returns exact x² as a QuadExpr (for NLP-capable solvers or testing)." -struct NoQuadApproxConfig <: QuadraticApproxConfig end - -""" - _add_quadratic_approx!(::NoQuadApproxConfig, container, C, names, time_steps, x_var, bounds, meta) - -No-op quadratic approximation: returns exact x² as a QuadExpr. - -# Arguments -- `::NoQuadApproxConfig`: no-op configuration (no fields) -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `names::Vector{String}`: component names -- `time_steps::UnitRange{Int}`: time periods -- `x_var`: container of variables indexed by (name, t) -- `bounds::Vector{MinMax}`: per-name lower and upper bounds of x domain -- `meta::String`: variable type identifier for the approximation -""" -function _add_quadratic_approx!( - ::NoQuadApproxConfig, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_var, - bounds::Vector{MinMax}, - meta::String, -) where {C <: IS.InfrastructureSystemsComponent} - result_expr = add_expression_container!( - container, - QuadraticExpression, - C, - names, - time_steps; - meta, - expr_type = JuMP.QuadExpr, - ) - for name in names, t in time_steps - result_expr[name, t] = x_var[name, t] * x_var[name, t] - end - return result_expr -end diff --git a/src/quadratic_approximations/pwl_utils.jl b/src/quadratic_approximations/pwl_utils.jl deleted file mode 100644 index ab68812d..00000000 --- a/src/quadratic_approximations/pwl_utils.jl +++ /dev/null @@ -1,60 +0,0 @@ -# Shared utilities for piecewise linear approximation methods. - -""" - _get_breakpoints_for_pwl_function(min_val, max_val, f; num_segments = DEFAULT_INTERPOLATION_LENGTH) - -Generate breakpoints for piecewise linear (PWL) approximation of a nonlinear function. - -This function creates equally-spaced breakpoints over the specified domain [min_val, max_val] -and evaluates the given function at each breakpoint to construct a piecewise linear approximation. -The breakpoints are used in optimization problems to linearize nonlinear constraints or objectives. - -# Arguments -- `min_val::Float64`: Minimum value of the domain for the PWL approximation -- `max_val::Float64`: Maximum value of the domain for the PWL approximation -- `f`: Function to be approximated (must be callable with Float64 input) -- `num_segments::Int`: Number of linear segments in the PWL approximation (default: DEFAULT_INTERPOLATION_LENGTH) - -# Returns -- `Tuple{Vector{Float64}, Vector{Float64}}`: A tuple containing: - - `x_bkpts`: Vector of x-coordinates (breakpoints) in the domain - - `y_bkpts`: Vector of y-coordinates (function values at breakpoints) - -# Notes -- The number of breakpoints is `num_segments + 1` -- Breakpoints are equally spaced across the domain -- The first breakpoint is always at `min_val` and the last at `max_val` -""" -function _get_breakpoints_for_pwl_function( - min_val::Float64, - max_val::Float64, - f; - num_segments = DEFAULT_INTERPOLATION_LENGTH, -) - # Calculate total number of breakpoints (one more than segments) - # num_segments is the number of linear segments in the PWL approximation - # num_bkpts is the total number of breakpoints needed for the segments - num_bkpts = num_segments + 1 - - # Calculate step size for equally-spaced breakpoints - step = (max_val - min_val) / num_segments - - # Pre-allocate vectors for breakpoint coordinates - x_bkpts = Vector{Float64}(undef, num_bkpts) # Domain values (x-coordinates) - y_bkpts = Vector{Float64}(undef, num_bkpts) # Function values (y-coordinates) - - # Set the first breakpoint at the minimum domain value - x_bkpts[1] = min_val - y_bkpts[1] = f(min_val) - - # Generate remaining breakpoints by stepping through the domain - for i in 1:num_segments - x_val = min_val + step * i # Calculate x-coordinate of current breakpoint - x_bkpts[i + 1] = x_val - y_bkpts[i + 1] = f(x_val) # Evaluate function at current breakpoint - end - return x_bkpts, y_bkpts -end - -"Helper: returns x² (used as the default function for PWL breakpoint generation)." -_square(x::Float64) = x * x diff --git a/src/quadratic_approximations/pwmcc_cuts.jl b/src/quadratic_approximations/pwmcc_cuts.jl deleted file mode 100644 index a636e6f8..00000000 --- a/src/quadratic_approximations/pwmcc_cuts.jl +++ /dev/null @@ -1,231 +0,0 @@ -# Piecewise McCormick (PWMCC) cuts for concave terms in Bin2 bilinear approximation. -# Adds K local chord upper bounds on the v^2 SoS2 approximation by partitioning -# each concave term's domain into K sub-intervals. -# LP gap shrinks from Delta^2/4 to Delta^2/(4K^2). -# These cuts supplement (do not replace) existing SoS2 constraints. - -"Binary interval selector for piecewise McCormick cuts." -struct PiecewiseMcCormickBinary <: SparseVariableType end - -"Disaggregated variable for piecewise McCormick cuts." -struct PiecewiseMcCormickDisaggregated <: SparseVariableType end - -"Selector sum constraint: sum_k delta_k = 1." -struct PiecewiseMcCormickSelectorSum <: ConstraintType end - -"Disaggregation linking constraint: v = sum_k v^d_k." -struct PiecewiseMcCormickLinking <: ConstraintType end - -"Interval activation lower bound: t_{k-1} * delta_k <= v^d_k." -struct PiecewiseMcCormickIntervalLB <: ConstraintType end - -"Interval activation upper bound: v^d_k <= t_k * delta_k." -struct PiecewiseMcCormickIntervalUB <: ConstraintType end - -"Piecewise McCormick chord upper-bound constraint on v^2 approximation." -struct PiecewiseMcCormickChordUB <: ConstraintType end - -"Piecewise McCormick tangent lower-bound constraint (left endpoint)." -struct PiecewiseMcCormickTangentLBL <: ConstraintType end - -"Piecewise McCormick tangent lower-bound constraint (right endpoint)." -struct PiecewiseMcCormickTangentLBR <: ConstraintType end - -""" - _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 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. - -# Arguments -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `names::Vector{String}`: component names -- `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)) -- `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" -""" -function _add_pwmcc_concave_cuts!( - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - v_var, - q_expr, - bounds::Vector{MinMax}, - K::Int, - meta::String, -) where {C <: IS.InfrastructureSystemsComponent} - IS.@assert_op K >= 1 - - jump_model = get_jump_model(container) - - # Create containers - 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, - C, - names, - time_steps; - meta, - ) - linking_cons = add_constraints_container!( - container, - PiecewiseMcCormickLinking, - C, - names, - time_steps; - meta, - ) - interval_lb_cons = add_constraints_container!( - container, - PiecewiseMcCormickIntervalLB, - C, - names, - 1:K, - time_steps; - meta, - ) - interval_ub_cons = add_constraints_container!( - container, - PiecewiseMcCormickIntervalUB, - C, - names, - 1:K, - time_steps; - meta, - ) - chord_ub_cons = add_constraints_container!( - container, - PiecewiseMcCormickChordUB, - C, - names, - time_steps; - meta, - ) - tangent_lb_l_cons = add_constraints_container!( - container, - PiecewiseMcCormickTangentLBL, - C, - names, - time_steps; - meta, - ) - tangent_lb_r_cons = add_constraints_container!( - container, - PiecewiseMcCormickTangentLBR, - C, - names, - time_steps; - meta, - ) - - delta = Vector{JuMP.VariableRef}(undef, K) - vd = Vector{JuMP.VariableRef}(undef, K) - - 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] - - for k in 1:K - delta[k] = - 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( - jump_model, - base_name = "PwMcCDis_$(C)_{$(name), $(k), $(t)}", - ) - end - - sel_expr = JuMP.AffExpr(0.0) - for k in 1:K - JuMP.add_to_expression!(sel_expr, delta[k]) - end - selector_cons[name, t] = JuMP.@constraint(jump_model, sel_expr == 1.0) - - link_expr = JuMP.AffExpr(0.0) - for k in 1:K - JuMP.add_to_expression!(link_expr, vd[k]) - end - linking_cons[name, t] = JuMP.@constraint(jump_model, link_expr == v) - - for k in 1:K - interval_lb_cons[name, k, t] = JuMP.@constraint( - jump_model, - brk[k] * delta[k] <= vd[k] - ) - interval_ub_cons[name, k, t] = JuMP.@constraint( - jump_model, - vd[k] <= brk[k + 1] * delta[k] - ) - end - - # Chord upper bound: prevents q from exceeding the local piecewise chord - # of v^2 in the LP relaxation (tightens from global chord to piecewise). - chord_rhs = JuMP.AffExpr(0.0) - for k in 1:K - JuMP.add_to_expression!(chord_rhs, sum_brk[k], vd[k]) - JuMP.add_to_expression!(chord_rhs, -prod_brk[k], delta[k]) - end - chord_ub_cons[name, t] = JuMP.@constraint(jump_model, q <= chord_rhs) - - # Tangent lower bounds from convexity of v^2 at interval endpoints. - tang_l_rhs = JuMP.AffExpr(0.0) - for k in 1:K - JuMP.add_to_expression!(tang_l_rhs, two_brk_l[k], vd[k]) - JuMP.add_to_expression!(tang_l_rhs, -sq_brk_l[k], delta[k]) - end - tangent_lb_l_cons[name, t] = JuMP.@constraint(jump_model, q >= tang_l_rhs) - - tang_r_rhs = JuMP.AffExpr(0.0) - for k in 1:K - JuMP.add_to_expression!(tang_r_rhs, two_brk_r[k], vd[k]) - JuMP.add_to_expression!(tang_r_rhs, -sq_brk_r[k], delta[k]) - end - tangent_lb_r_cons[name, t] = JuMP.@constraint(jump_model, q >= tang_r_rhs) - end - - return -end diff --git a/src/quadratic_approximations/sawtooth.jl b/src/quadratic_approximations/sawtooth.jl deleted file mode 100644 index 0a0d0d2e..00000000 --- a/src/quadratic_approximations/sawtooth.jl +++ /dev/null @@ -1,227 +0,0 @@ -# Sawtooth MIP approximation of x² for use in constraints. -# Uses recursive tooth function compositions with O(log(1/ε)) binary variables. -# Reference: Beach, Burlacu, Hager, Hildebrand (2024). - -"Auxiliary continuous variables (g₀, …, g_L) for sawtooth quadratic approximation." -struct SawtoothAuxVariable <: VariableType end -"Binary variables (α₁, …, α_L) for sawtooth quadratic approximation." -struct SawtoothBinaryVariable <: VariableType end -"Variable result in tightened version." -struct SawtoothTightenedVariable <: VariableType end -"Links g₀ to the normalized x value in sawtooth quadratic approximation." -struct SawtoothLinkingConstraint <: ConstraintType end -"Constrains g_j based on g_{j-1}." -struct SawtoothMIPConstraint <: ConstraintType end -"LP relaxation constraints (g_j ≤ 2g_{j-1}, g_j ≤ 2(1−g_{j-1})) used in epigraph tightening." -struct SawtoothLPConstraint <: ConstraintType end -"Bounds tightened variable." -struct SawtoothTightenedConstraint <: ConstraintType end - -""" -Config for sawtooth MIP quadratic approximation. - -# Fields -- `depth::Int`: recursion depth L; uses L binary variables for 2^L + 1 breakpoints -- `epigraph_depth::Int`: LP tightening depth via epigraph Q^{L1} lower bound; 0 to disable (default 0) -""" -struct SawtoothQuadConfig <: QuadraticApproxConfig - depth::Int - epigraph_depth::Int -end -SawtoothQuadConfig(depth::Int) = SawtoothQuadConfig(depth, 0) - -""" - _add_quadratic_approx!(config::SawtoothQuadConfig, container, C, names, time_steps, x_var, bounds, meta) - -Approximate x² using the sawtooth MIP formulation. - -Creates auxiliary continuous variables g_0,...,g_L and binary variables α_1,...,α_L, -adds S^L constraints (4 per level) and a linking constraint for each component and -time step, and stores affine expressions approximating x² in a -`QuadraticExpression` expression container. - -For depth L, the approximation interpolates x² at 2^L + 1 uniformly spaced breakpoints -with maximum overestimation error Δ² · 2^{-2L-2} where Δ = x_max - x_min. - -# Arguments -- `config::SawtoothQuadConfig`: configuration with `depth` (recursion depth L; uses L binary variables) and `epigraph_depth` (LP tightening depth; 0 to disable) -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `names::Vector{String}`: component names -- `time_steps::UnitRange{Int}`: time periods -- `x_var`: container of variables indexed by (name, t) -- `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!( - config::SawtoothQuadConfig, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_var, - bounds::Vector{MinMax}, - meta::String, -) where {C <: IS.InfrastructureSystemsComponent} - IS.@assert_op config.depth >= 1 - jump_model = get_jump_model(container) - - # Create containers with known dimensions - g_levels = 0:(config.depth) - alpha_levels = 1:(config.depth) - g_var = add_variable_container!( - container, - SawtoothAuxVariable, - C, - names, - g_levels, - time_steps; - meta, - ) - alpha_var = add_variable_container!( - container, - SawtoothBinaryVariable, - C, - names, - alpha_levels, - time_steps; - meta, - ) - mip_cons = add_constraints_container!( - container, - SawtoothMIPConstraint, - C, - names, - 1:4, - time_steps; - sparse = true, - meta, - ) - link_cons = add_constraints_container!( - container, - SawtoothLinkingConstraint, - C, - names, - time_steps; - meta, - ) - result_expr = add_expression_container!( - container, - QuadraticExpression, - C, - names, - time_steps; - meta, - ) - - if config.epigraph_depth > 0 - lp_expr = _add_quadratic_approx!( - EpigraphQuadConfig(config.epigraph_depth), - container, C, names, time_steps, - x_var, bounds, meta * "_lb", - ) - z_var = add_variable_container!( - container, - SawtoothTightenedVariable, - C, - names, - time_steps; - meta, - ) - tight_cons = add_constraints_container!( - container, - SawtoothTightenedConstraint, - C, - names, - 1:2, - time_steps; - meta, - ) - end - - 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] - for j in g_levels - g_var[name, j, t] = JuMP.@variable( - jump_model, - base_name = "SawtoothAux_$(C)_{$(name), $(j), $(t)}", - lower_bound = 0.0, - upper_bound = 1.0, - ) - end - - # Binary variables α_1,...,α_L - for j in alpha_levels - alpha_var[name, j, t] = JuMP.@variable( - jump_model, - base_name = "SawtoothBin_$(C)_{$(name), $(j), $(t)}", - binary = true, - ) - end - - # Linking constraint: g_0 = (x - x_min) / Δ - link_cons[name, t] = JuMP.@constraint( - jump_model, - g_var[name, 0, t] == (x - b.min) / delta, - ) - - # S^L constraints for j = 1,...,L - for j in alpha_levels - g_prev = g_var[name, j - 1, t] - g_curr = g_var[name, j, t] - alpha_j = alpha_var[name, j, t] - - # g_j ≤ 2 g_{j-1} - mip_cons[name, 1, t] = JuMP.@constraint(jump_model, g_curr <= 2.0 * g_prev) - # g_j ≤ 2(1 - g_{j-1}) - mip_cons[name, 2, t] = - JuMP.@constraint(jump_model, g_curr <= 2.0 * (1.0 - g_prev)) - # g_j ≥ 2(g_{j-1} - α_j) - mip_cons[name, 3, t] = - JuMP.@constraint(jump_model, g_curr >= 2.0 * (g_prev - alpha_j)) - # g_j ≥ 2(α_j - g_{j-1}) - mip_cons[name, 4, t] = - JuMP.@constraint(jump_model, g_curr >= 2.0 * (alpha_j - g_prev)) - end - - # Build x² ≈ x_min² + (2 x_min Δ + Δ²) g_0 - Σ_{j=1}^L Δ² 2^{-2j} g_j - x_sq_approx = JuMP.AffExpr(b.min * b.min) - add_proportional_to_jump_expression!( - x_sq_approx, - g_var[name, 0, t], - 2.0 * b.min * delta + delta * delta, - ) - for j in alpha_levels - add_proportional_to_jump_expression!( - x_sq_approx, - g_var[name, j, t], - -saw_coeffs[j], - ) - end - - if config.epigraph_depth > 0 - z = - z_var[name, t] = JuMP.@variable( - jump_model, - base_name = "TightenedSawtooth_$(C)_{$(name), $(t)}", - lower_bound = z_min, - upper_bound = z_max - ) - tight_cons[name, 1, t] = JuMP.@constraint(jump_model, z <= x_sq_approx) - tight_cons[name, 2, t] = JuMP.@constraint(jump_model, z >= lp_expr[name, t]) - result_expr[name, t] = JuMP.AffExpr(0.0, z => 1.0) - else - result_expr[name, t] = x_sq_approx - end - end - - return result_expr -end diff --git a/src/quadratic_approximations/solver_sos2.jl b/src/quadratic_approximations/solver_sos2.jl deleted file mode 100644 index 110428fe..00000000 --- a/src/quadratic_approximations/solver_sos2.jl +++ /dev/null @@ -1,186 +0,0 @@ -# SOS2-based piecewise linear approximation of x² for use in constraints. -# Uses solver-native MOI.SOS2 constraints for adjacency enforcement. - -"lambda_var (λ) convex combination weight variables for SOS2 quadratic approximation." -struct QuadraticVariable <: SparseVariableType end -"Links x to the weighted sum of breakpoints in SOS2 quadratic approximation." -struct SOS2LinkingConstraint <: ConstraintType end -"Expression for the weighted sum of breakpoints Σ λ_i * x_i linking x to lambda variables." -struct SOS2LinkingExpression <: ExpressionType end -"Ensures the sum of λ weights equals 1 in SOS2 quadratic approximation." -struct SOS2NormConstraint <: ConstraintType end -"Expression for the normalization sum Σ λ_i in SOS2 quadratic approximation." -struct SOS2NormExpression <: ExpressionType end - -"Solver-native MOI.SOS2 adjacency constraint on lambda variables." -struct SolverSOS2Constraint <: ConstraintType end - -""" -Config for solver-native SOS2 quadratic approximation (MOI.SOS2 adjacency). - -# Fields -- `depth::Int`: number of PWL segments (breakpoints = depth + 1) -- `pwmcc_segments::Int`: number of piecewise McCormick cut partitions; 0 to disable (default 4) -""" -struct SolverSOS2QuadConfig <: QuadraticApproxConfig - depth::Int - pwmcc_segments::Int -end -SolverSOS2QuadConfig(depth::Int) = SolverSOS2QuadConfig(depth, 4) - -""" - _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. - -Creates lambda_var (λ) variables representing convex combination weights over breakpoints, -adds linking, normalization, and MOI.SOS2 constraints, and stores affine expressions -approximating x² in a `QuadraticExpression` expression container. - -# Arguments -- `config::SolverSOS2QuadConfig`: configuration with `depth` (number of PWL segments) and `pwmcc_segments` (PWMCC cut partitions; 0 to disable, default 4) -- `container::OptimizationContainer`: the optimization container -- `::Type{C}`: component type -- `names::Vector{String}`: component names -- `time_steps::UnitRange{Int}`: time periods -- `x_var`: container of variables indexed by (name, t) -- `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!( - config::SolverSOS2QuadConfig, - container::OptimizationContainer, - ::Type{C}, - names::Vector{String}, - time_steps::UnitRange{Int}, - x_var, - bounds::Vector{MinMax}, - meta::String, -) where {C <: IS.InfrastructureSystemsComponent} - x_bkpts, x_sq_bkpts = - _get_breakpoints_for_pwl_function( - 0.0, - 1.0, - _square; - num_segments = config.depth, - ) - n_points = config.depth + 1 - jump_model = get_jump_model(container) - - # Create all containers upfront - lambda_var = add_variable_container!( - container, - QuadraticVariable, - C, - names, - 1:n_points, - time_steps; - meta, - ) - link_cons = add_constraints_container!( - container, - SOS2LinkingConstraint, - C, - names, - time_steps; - meta, - ) - link_expr = add_expression_container!( - container, - SOS2LinkingExpression, - C, - names, - time_steps; - meta, - ) - norm_cons = add_constraints_container!( - container, - SOS2NormConstraint, - C, - names, - time_steps; - meta, - ) - norm_expr = add_expression_container!( - container, - SOS2NormExpression, - C, - names, - time_steps; - meta, - ) - sos_cons = add_constraints_container!( - container, - SolverSOS2Constraint, - C, - names, - time_steps; - meta, - ) - result_expr = add_expression_container!( - container, - QuadraticExpression, - C, - names, - time_steps; - meta, - ) - - 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 - x = x_var[name, t] - - # Create lambda_var variables: λ_i ∈ [0, 1] - lambda = Vector{JuMP.VariableRef}(undef, n_points) - for i in 1:n_points - lambda[i] = - lambda_var[name, i, t] = JuMP.@variable( - jump_model, - base_name = "QuadraticVariable_$(C)_{$(name), pwl_$(i), $(t)}", - lower_bound = 0.0, - upper_bound = 1.0, - ) - end - - # x = Σ λ_i * x_i - link = link_expr[name, t] = JuMP.AffExpr(0.0) - 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 - b.min) / lx == link) - - # Σ λ_i = 1 - norm = norm_expr[name, t] = JuMP.AffExpr(0.0) - for l in lambda - add_proportional_to_jump_expression!(norm, l, 1.0) - end - norm_cons[name, t] = JuMP.@constraint(jump_model, norm == 1.0) - - # λ ∈ SOS2 (solver-native) - sos_cons[name, t] = - JuMP.@constraint(jump_model, lambda in MOI.SOS2(collect(1:n_points))) - - # Build x̂² = Σ λ_i * x_i² as an affine expression - x_hat_sq = JuMP.AffExpr(0.0) - for i in 1:n_points - add_proportional_to_jump_expression!(x_hat_sq, lambda[i], x_sq_bkpts[i]) - 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 * 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, bounds, - config.pwmcc_segments, meta * "_pwmcc", - ) - end - - return result_expr -end diff --git a/test/InfrastructureOptimizationModelsTests.jl b/test/InfrastructureOptimizationModelsTests.jl index 97968b5e..7b761aa9 100644 --- a/test/InfrastructureOptimizationModelsTests.jl +++ b/test/InfrastructureOptimizationModelsTests.jl @@ -133,7 +133,8 @@ function run_tests() include(joinpath(TEST_DIR, "test_jump_utils.jl")) include(joinpath(TEST_DIR, "test_pwl_methods.jl")) - # --- quadratic_approximations/ subfolder --- + # --- approximations/ subfolder --- + # IOM-wrapper regression tests (exercise add_*_approx! end-to-end). include(joinpath(TEST_DIR, "test_quadratic_approximations.jl")) include(joinpath(TEST_DIR, "test_bilinear_approximations.jl")) include(joinpath(TEST_DIR, "test_hybs_approximations.jl")) diff --git a/test/performance/bilinear_delta_benchmark.jl b/test/performance/bilinear_delta_benchmark.jl index d51f6f84..d96439b5 100644 --- a/test/performance/bilinear_delta_benchmark.jl +++ b/test/performance/bilinear_delta_benchmark.jl @@ -234,7 +234,7 @@ function build_gen_bilinear( container, net::MockNetworkProblem, V_container, I_container, time_steps, bilinear_config::Union{IOM.Bin2Config, IOM.HybSConfig}, quad_config, ) - V_sq = IOM._add_quadratic_approx!( + V_sq = IOM.add_quadratic_approx!( quad_config, container, MockNetworkNode, @@ -245,7 +245,7 @@ function build_gen_bilinear( V_MAX, "gen_V_sq", ) - I_sq = IOM._add_quadratic_approx!( + I_sq = IOM.add_quadratic_approx!( quad_config, container, MockNetworkNode, @@ -256,7 +256,7 @@ function build_gen_bilinear( I_GEN_MAX, "gen_I_sq", ) - z_gen = IOM._add_bilinear_approx!( + z_gen = IOM.add_bilinear_approx!( bilinear_config, container, MockNetworkNode, @@ -291,11 +291,11 @@ function build_gen_bilinear( container, MockNetworkNode, net.gen_nodes, time_steps, I_container, I_GEN_MIN, I_GEN_MAX, quad_config.depth, "gen_I", ) - I_sq = IOM._add_quadratic_approx!( + I_sq = IOM.add_quadratic_approx!( quad_config, container, MockNetworkNode, net.gen_nodes, time_steps, I_disc, I_GEN_MIN, I_GEN_MAX, "gen_I_sq", ) - z_gen = IOM._add_bilinear_approx!( + z_gen = IOM.add_bilinear_approx!( bilinear_config, container, MockNetworkNode, net.gen_nodes, time_steps, V_disc, I_disc, V_MIN, V_MAX, I_GEN_MIN, I_GEN_MAX, "gen", ) @@ -309,7 +309,7 @@ function build_gen_bilinear( container, net::MockNetworkProblem, V_container, I_container, time_steps, bilinear_config::IOM.NoBilinearApproxConfig, quad_config::IOM.NoQuadApproxConfig, ) - z_gen = IOM._add_bilinear_approx!( + z_gen = IOM.add_bilinear_approx!( bilinear_config, container, MockNetworkNode, @@ -323,7 +323,7 @@ function build_gen_bilinear( I_GEN_MAX, "gen", ) - I_sq = IOM._add_quadratic_approx!( + I_sq = IOM.add_quadratic_approx!( quad_config, container, MockNetworkNode, @@ -393,7 +393,7 @@ function build_mip_model( ) # --- Bilinear dem: always uses the config-based dispatch --- - z_dem = IOM._add_bilinear_approx!( + z_dem = IOM.add_bilinear_approx!( bilinear_config, container, MockNetworkNode, net.dem_nodes, time_steps, V_container, I_container, V_MIN, V_MAX, I_DEM_MIN, I_DEM_MAX, "dem", diff --git a/test/test_bilinear_approximations.jl b/test/test_bilinear_approximations.jl index 8619eb5d..88081904 100644 --- a/test/test_bilinear_approximations.jl +++ b/test/test_bilinear_approximations.jl @@ -12,12 +12,10 @@ const BILINEAR_META = "BilinearTest" JuMP.fix(setup.x_var_container["dev1", 1], 2.5; force = true) JuMP.fix(setup.y_var_container["dev1", 1], 1.5; force = true) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.Bin2Config(IOM.SolverSOS2QuadConfig(4, 0), add_mc), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 4.0)], @@ -53,12 +51,10 @@ const BILINEAR_META = "BilinearTest" JuMP.set_lower_bound(y_var, 0.0) JuMP.set_upper_bound(y_var, 4.0) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.Bin2Config(IOM.SolverSOS2QuadConfig(8, 0)), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 4.0)], @@ -86,12 +82,10 @@ const BILINEAR_META = "BilinearTest" JuMP.fix(setup2.x_var_container["dev1", 1], 2.0; force = true) JuMP.fix(setup2.y_var_container["dev1", 1], 3.0; force = true) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.Bin2Config(IOM.SolverSOS2QuadConfig(8, 0)), setup2.container, MockThermalGen, - ["dev1"], - 1:1, setup2.x_var_container, setup2.y_var_container, [(min = 0.0, max = 4.0)], @@ -124,12 +118,10 @@ const BILINEAR_META = "BilinearTest" w = JuMP.@variable(setup.jump_model, base_name = "w") - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.Bin2Config(IOM.SolverSOS2QuadConfig(8, 0)), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 4.0)], @@ -165,12 +157,10 @@ const BILINEAR_META = "BilinearTest" JuMP.set_lower_bound(y_var, 0.0) JuMP.set_upper_bound(y_var, 4.0) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.Bin2Config(IOM.SolverSOS2QuadConfig(8, 0)), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 4.0)], @@ -206,12 +196,10 @@ const BILINEAR_META = "BilinearTest" JuMP.fix(x_var, 2.5; force = true) JuMP.fix(y_var, 1.5; force = true) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.Bin2Config(IOM.SolverSOS2QuadConfig(num_segments, 0)), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 4.0)], @@ -246,12 +234,10 @@ const BILINEAR_META = "BilinearTest" 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.add_bilinear_approx!( IOM.Bin2Config(IOM.ManualSOS2QuadConfig(8, 0)), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 4.0)], @@ -282,12 +268,10 @@ const BILINEAR_META = "BilinearTest" w = JuMP.@variable(setup.jump_model, base_name = "w") - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.Bin2Config(IOM.ManualSOS2QuadConfig(8, 0)), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 4.0)], @@ -321,12 +305,10 @@ const BILINEAR_META = "BilinearTest" JuMP.set_lower_bound(y_var, 0.0) JuMP.set_upper_bound(y_var, 4.0) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.Bin2Config(IOM.ManualSOS2QuadConfig(8, 0)), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 4.0)], @@ -358,12 +340,10 @@ const BILINEAR_META = "BilinearTest" JuMP.fix(setup.x_var_container["dev1", 1], 2.5; force = true) JuMP.fix(setup.y_var_container["dev1", 1], 1.5; force = true) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.Bin2Config(IOM.ManualSOS2QuadConfig(num_segments, 0)), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 4.0)], @@ -398,12 +378,10 @@ const BILINEAR_META = "BilinearTest" 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.add_bilinear_approx!( IOM.Bin2Config(IOM.SawtoothQuadConfig(3)), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 4.0)], @@ -434,12 +412,10 @@ const BILINEAR_META = "BilinearTest" w = JuMP.@variable(setup.jump_model, base_name = "w") - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.Bin2Config(IOM.SawtoothQuadConfig(3)), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 4.0)], @@ -473,12 +449,10 @@ const BILINEAR_META = "BilinearTest" JuMP.set_lower_bound(y_var, 0.0) JuMP.set_upper_bound(y_var, 4.0) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.Bin2Config(IOM.SawtoothQuadConfig(3)), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 4.0)], @@ -510,12 +484,10 @@ const BILINEAR_META = "BilinearTest" JuMP.fix(setup.x_var_container["dev1", 1], 2.5; force = true) JuMP.fix(setup.y_var_container["dev1", 1], 1.5; force = true) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.Bin2Config(IOM.SawtoothQuadConfig(depth)), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 4.0)], diff --git a/test/test_hybs_approximations.jl b/test/test_hybs_approximations.jl index da2ab025..aee3aa25 100644 --- a/test/test_hybs_approximations.jl +++ b/test/test_hybs_approximations.jl @@ -7,12 +7,10 @@ const HYBS_BILINEAR_META = "BilinearTest" x_var = setup.var_container["dev1", 1] JuMP.fix(x_var, 0.35; force = true) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.EpigraphQuadConfig(4), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.var_container, [(min = 0.0, max = 1.0)], HYBS_META, @@ -40,12 +38,10 @@ const HYBS_BILINEAR_META = "BilinearTest" x_var = setup.var_container["dev1", 1] JuMP.fix(x_var, 1.3; force = true) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.EpigraphQuadConfig(4), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.var_container, [(min = 0.0, max = 2.0)], HYBS_META, @@ -76,12 +72,10 @@ const HYBS_BILINEAR_META = "BilinearTest" x_var = setup.var_container["dev1", 1] JuMP.fix(x_var, 0.35; force = true) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.EpigraphQuadConfig(depth), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.var_container, [(min = 0.0, max = 1.0)], HYBS_META, @@ -118,12 +112,10 @@ end 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.add_bilinear_approx!( IOM.HybSConfig(IOM.SawtoothQuadConfig(2), 2), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 1.0)], @@ -159,12 +151,10 @@ end JuMP.fix(x_var, 2.0; force = true) JuMP.fix(y_var, 3.0; force = true) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.HybSConfig(IOM.SawtoothQuadConfig(3), 3), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 4.0)], @@ -197,12 +187,10 @@ end w = JuMP.@variable(setup.jump_model, base_name = "w") - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.HybSConfig(IOM.SawtoothQuadConfig(3), 3), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 4.0)], @@ -235,12 +223,10 @@ end JuMP.fix(setup.x_var_container["dev1", 1], 0.4; force = true) JuMP.fix(setup.y_var_container["dev1", 1], 0.7; force = true) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.HybSConfig(IOM.SawtoothQuadConfig(depth), depth), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 1.0)], @@ -276,12 +262,10 @@ end JuMP.fix(setup.x_var_container["dev1", 1], 3.5; force = true) JuMP.fix(setup.y_var_container["dev1", 1], 2.1; force = true) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.HybSConfig(IOM.SawtoothQuadConfig(3), 3), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 2.0, max = 5.0)], @@ -317,12 +301,10 @@ end JuMP.set_lower_bound(y_var, 0.0) JuMP.set_upper_bound(y_var, 4.0) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.HybSConfig(IOM.SawtoothQuadConfig(2), 2), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 4.0)], @@ -351,12 +333,10 @@ end for depth in [1, 2, 4] # HybS setup_h = _setup_bilinear_test(["dev1"], 1:1) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.HybSConfig(IOM.SawtoothQuadConfig(depth), depth), setup_h.container, MockThermalGen, - ["dev1"], - 1:1, setup_h.x_var_container, setup_h.y_var_container, [(min = 0.0, max = 1.0)], @@ -368,12 +348,10 @@ end # Bin2 (sawtooth) setup_b = _setup_bilinear_test(["dev1"], 1:1) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.Bin2Config(IOM.SawtoothQuadConfig(depth)), setup_b.container, MockThermalGen, - ["dev1"], - 1:1, setup_b.x_var_container, setup_b.y_var_container, [(min = 0.0, max = 1.0)], diff --git a/test/test_nmdt_approximations.jl b/test/test_nmdt_approximations.jl index a5c2ecfa..05f233e6 100644 --- a/test/test_nmdt_approximations.jl +++ b/test/test_nmdt_approximations.jl @@ -12,9 +12,9 @@ const NMDT_BILINEAR_META = "NMDTBilinearTest" JuMP.set_upper_bound(setup.var_container["gen1", 1], 1.0) JuMP.fix(setup.var_container["gen1", 1], 0.6; force = true) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.DNMDTQuadConfig(4, 0), - setup.container, MockThermalGen, names, ts, + setup.container, MockThermalGen, setup.var_container, [(min = 0.0, max = 1.0)], DNMDT_META, ) @@ -46,9 +46,9 @@ const NMDT_BILINEAR_META = "NMDTBilinearTest" setup = _setup_qa_test(["gen1"], 1:1) JuMP.fix(setup.var_container["gen1", 1], x0; force = true) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.DNMDTQuadConfig(3, 0), - setup.container, MockThermalGen, ["gen1"], 1:1, + setup.container, MockThermalGen, setup.var_container, [(min = 0.0, max = 1.0)], DNMDT_META, ) expr = IOM.get_expression( @@ -77,9 +77,9 @@ const NMDT_BILINEAR_META = "NMDTBilinearTest" setup = _setup_qa_test(["gen1"], 1:1) JuMP.fix(setup.var_container["gen1", 1], x0; force = true) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.DNMDTQuadConfig(2 * L, 0), - setup.container, MockThermalGen, ["gen1"], 1:1, + setup.container, MockThermalGen, setup.var_container, [(min = 0.0, max = 1.0)], DNMDT_META, ) expr = IOM.get_expression( @@ -109,9 +109,9 @@ end setup = _setup_qa_test(["gen1"], 1:1) JuMP.fix(setup.var_container["gen1", 1], x0; force = true) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( (tighten ? IOM.DNMDTQuadConfig(2) : IOM.DNMDTQuadConfig(2, 0)), - setup.container, MockThermalGen, ["gen1"], 1:1, + setup.container, MockThermalGen, setup.var_container, [(min = 0.0, max = 1.0)], DNMDT_META, ) expr = IOM.get_expression( @@ -146,9 +146,9 @@ end setup = _setup_qa_test(["gen1"], 1:1) JuMP.fix(setup.var_container["gen1", 1], 0.35; force = true) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.DNMDTQuadConfig(L), - setup.container, MockThermalGen, ["gen1"], 1:1, + setup.container, MockThermalGen, setup.var_container, [(min = 0.0, max = 1.0)], DNMDT_META, ) expr = IOM.get_expression( @@ -180,9 +180,9 @@ end 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.add_bilinear_approx!( IOM.DNMDTBilinearConfig(2), - setup.container, MockThermalGen, ["dev1"], 1:1, + setup.container, MockThermalGen, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], DNMDT_META, ) @@ -214,9 +214,9 @@ end 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.add_bilinear_approx!( IOM.DNMDTBilinearConfig(2 * L), - setup.container, MockThermalGen, ["dev1"], 1:1, + setup.container, MockThermalGen, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], DNMDT_META, ) @@ -249,9 +249,9 @@ end 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.add_bilinear_approx!( IOM.DNMDTBilinearConfig(8), - setup.container, MockThermalGen, ["dev1"], 1:1, + setup.container, MockThermalGen, setup.x_var_container, setup.y_var_container, [(min = x_min, max = x_max)], [(min = y_min, max = y_max)], DNMDT_META, ) @@ -276,8 +276,8 @@ end # @testset "McCormick toggle" begin # setup = _setup_bilinear_test(["dev1"], 1:1) - # IOM._add_bilinear_approx!( - # setup.container, MockThermalGen, ["dev1"], 1:1, + # IOM.add_bilinear_approx!( + # setup.container, MockThermalGen, # setup.x_var_container, setup.y_var_container, # 0.0, 1.0, 0.0, 1.0, 2, DNMDT_META; # add_mccormick = false, @@ -293,9 +293,9 @@ end 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.add_bilinear_approx!( IOM.DNMDTBilinearConfig(3), - setup.container, MockThermalGen, ["dev1"], 1:1, + setup.container, MockThermalGen, setup.y_var_container, setup.x_var_container, [(min = 0.0, max = 4.0)], [(min = 0.0, max = 4.0)], DNMDT_META, ) @@ -321,9 +321,9 @@ end JuMP.fix(setup_d.x_var_container["dev1", 1], 0.4; force = true) JuMP.fix(setup_d.y_var_container["dev1", 1], 0.7; force = true) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.DNMDTBilinearConfig(depth), - setup_d.container, MockThermalGen, ["dev1"], 1:1, + setup_d.container, MockThermalGen, setup_d.x_var_container, setup_d.y_var_container, [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], DNMDT_META, ) @@ -343,9 +343,9 @@ end JuMP.fix(setup_h.x_var_container["dev1", 1], 0.4; force = true) JuMP.fix(setup_h.y_var_container["dev1", 1], 0.7; force = true) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.HybSConfig(IOM.SawtoothQuadConfig(depth), depth), - setup_h.container, MockThermalGen, ["dev1"], 1:1, + setup_h.container, MockThermalGen, setup_h.x_var_container, setup_h.y_var_container, [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], DNMDT_HYBS_META, ) @@ -380,9 +380,9 @@ end JuMP.set_upper_bound(setup.var_container["gen1", 1], 1.0) JuMP.fix(setup.var_container["gen1", 1], 0.6; force = true) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.NMDTQuadConfig(4, 0), - setup.container, MockThermalGen, names, ts, + setup.container, MockThermalGen, setup.var_container, [(min = 0.0, max = 1.0)], NMDT_META, ) @@ -414,9 +414,9 @@ end setup = _setup_qa_test(["gen1"], 1:1) JuMP.fix(setup.var_container["gen1", 1], x0; force = true) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.NMDTQuadConfig(3, 0), - setup.container, MockThermalGen, ["gen1"], 1:1, + setup.container, MockThermalGen, setup.var_container, [(min = 0.0, max = 1.0)], NMDT_META, ) expr = IOM.get_expression( @@ -447,9 +447,9 @@ end setup = _setup_qa_test(["gen1"], 1:1) JuMP.fix(setup.var_container["gen1", 1], x0; force = true) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.NMDTQuadConfig(L, 0), - setup.container, MockThermalGen, ["gen1"], 1:1, + setup.container, MockThermalGen, setup.var_container, [(min = 0.0, max = 1.0)], NMDT_META, ) expr = IOM.get_expression( @@ -477,9 +477,9 @@ end setup = _setup_qa_test(["gen1"], 1:1) JuMP.fix(setup.var_container["gen1", 1], x0; force = true) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( (tighten ? IOM.NMDTQuadConfig(2) : IOM.NMDTQuadConfig(2, 0)), - setup.container, MockThermalGen, ["gen1"], 1:1, + setup.container, MockThermalGen, setup.var_container, [(min = 0.0, max = 1.0)], NMDT_META, ) expr = IOM.get_expression( @@ -523,9 +523,9 @@ end setup = _setup_qa_test(["gen1"], 1:1) JuMP.fix(setup.var_container["gen1", 1], x0; force = true) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( config_fn(L), - setup.container, MockThermalGen, ["gen1"], 1:1, + setup.container, MockThermalGen, setup.var_container, [(min = 0.0, max = 1.0)], NMDT_META, ) expr = IOM.get_expression( @@ -558,17 +558,17 @@ end # Both D-NMDT and NMDT univariate use L binary variables for depth=L depth = 4 setup_n = _setup_qa_test(["gen1"], 1:1) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.NMDTQuadConfig(depth, 0), - setup_n.container, MockThermalGen, ["gen1"], 1:1, + setup_n.container, MockThermalGen, 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)) setup_d = _setup_qa_test(["gen1"], 1:1) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.DNMDTQuadConfig(depth, 0), - setup_d.container, MockThermalGen, ["gen1"], 1:1, + setup_d.container, MockThermalGen, 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)) @@ -592,9 +592,9 @@ end 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.add_bilinear_approx!( IOM.NMDTBilinearConfig(3), - setup.container, MockThermalGen, ["dev1"], 1:1, + setup.container, MockThermalGen, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], NMDT_BILINEAR_META, ) @@ -627,9 +627,9 @@ end 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.add_bilinear_approx!( IOM.NMDTBilinearConfig(L), - setup.container, MockThermalGen, ["dev1"], 1:1, + setup.container, MockThermalGen, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], NMDT_BILINEAR_META, @@ -658,9 +658,9 @@ end JuMP.fix(setup.x_var_container["dev1", 1], 0.75; force = true) JuMP.fix(setup.y_var_container["dev1", 1], 0.5; force = true) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.NMDTBilinearConfig(4), - setup.container, MockThermalGen, ["dev1"], 1:1, + setup.container, MockThermalGen, setup.x_var_container, setup.y_var_container, [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], NMDT_BILINEAR_META, ) @@ -684,18 +684,18 @@ end @testset "NMDT uses L binaries, D-NMDT uses 2L" begin L = 3 setup_n = _setup_bilinear_test(["dev1"], 1:1) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.NMDTBilinearConfig(L), - setup_n.container, MockThermalGen, ["dev1"], 1:1, + setup_n.container, MockThermalGen, setup_n.x_var_container, setup_n.y_var_container, [(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)) setup_d = _setup_bilinear_test(["dev1"], 1:1) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.DNMDTBilinearConfig(L), - setup_d.container, MockThermalGen, ["dev1"], 1:1, + setup_d.container, MockThermalGen, setup_d.x_var_container, setup_d.y_var_container, [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], DNMDT_META, ) @@ -714,9 +714,9 @@ end JuMP.fix(setup_n.x_var_container["dev1", 1], 0.4; force = true) JuMP.fix(setup_n.y_var_container["dev1", 1], 0.7; force = true) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.NMDTBilinearConfig(L), - setup_n.container, MockThermalGen, ["dev1"], 1:1, + setup_n.container, MockThermalGen, setup_n.x_var_container, setup_n.y_var_container, [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], NMDT_BILINEAR_META, ) @@ -735,9 +735,9 @@ end JuMP.fix(setup_d.x_var_container["dev1", 1], 0.4; force = true) JuMP.fix(setup_d.y_var_container["dev1", 1], 0.7; force = true) - IOM._add_bilinear_approx!( + IOM.add_bilinear_approx!( IOM.DNMDTBilinearConfig(L), - setup_d.container, MockThermalGen, ["dev1"], 1:1, + setup_d.container, MockThermalGen, setup_d.x_var_container, setup_d.y_var_container, [(min = 0.0, max = 1.0)], [(min = 0.0, max = 1.0)], DNMDT_META, ) diff --git a/test/test_quadratic_approximations.jl b/test/test_quadratic_approximations.jl index 4bc8f4d0..b0e68ee5 100644 --- a/test/test_quadratic_approximations.jl +++ b/test/test_quadratic_approximations.jl @@ -11,12 +11,10 @@ const TEST_META = "TestVar" JuMP.set_lower_bound(x_var, 0.0) JuMP.set_upper_bound(x_var, 4.0) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.SolverSOS2QuadConfig(4, 0), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.var_container, [(min = 0.0, max = 4.0)], TEST_META, @@ -47,12 +45,10 @@ const TEST_META = "TestVar" y = JuMP.@variable(setup.jump_model, base_name = "y") - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.SolverSOS2QuadConfig(4, 0), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.var_container, [(min = 0.0, max = 4.0)], TEST_META, @@ -86,12 +82,10 @@ const TEST_META = "TestVar" JuMP.set_lower_bound(x_var, 0.0) JuMP.set_upper_bound(x_var, 6.0) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.SolverSOS2QuadConfig(num_segments, 0), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.var_container, [(min = 0.0, max = 6.0)], TEST_META, @@ -126,12 +120,10 @@ const TEST_META = "TestVar" JuMP.set_lower_bound(x_var, 0.0) JuMP.set_upper_bound(x_var, 4.0) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.ManualSOS2QuadConfig(4, 0), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.var_container, [(min = 0.0, max = 4.0)], TEST_META, @@ -161,12 +153,10 @@ const TEST_META = "TestVar" y = JuMP.@variable(setup.jump_model, base_name = "y") - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.ManualSOS2QuadConfig(4, 0), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.var_container, [(min = 0.0, max = 4.0)], TEST_META, @@ -198,12 +188,10 @@ const TEST_META = "TestVar" JuMP.set_lower_bound(x_var, 0.0) JuMP.set_upper_bound(x_var, 4.0) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.SawtoothQuadConfig(2), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.var_container, [(min = 0.0, max = 4.0)], TEST_META, @@ -233,12 +221,10 @@ const TEST_META = "TestVar" y = JuMP.@variable(setup.jump_model, base_name = "y") - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.SawtoothQuadConfig(2), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.var_container, [(min = 0.0, max = 4.0)], TEST_META, @@ -271,12 +257,10 @@ const TEST_META = "TestVar" JuMP.set_lower_bound(x_var, 0.0) JuMP.set_upper_bound(x_var, 6.0) - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.SawtoothQuadConfig(depth), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.var_container, [(min = 0.0, max = 6.0)], TEST_META, @@ -313,23 +297,19 @@ const TEST_META = "TestVar" JuMP.set_upper_bound(x_var, 4.0) if method == :sos2 - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.SolverSOS2QuadConfig(2^depth, 0), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.var_container, [(min = 0.0, max = 4.0)], TEST_META, ) else - IOM._add_quadratic_approx!( + IOM.add_quadratic_approx!( IOM.SawtoothQuadConfig(depth), setup.container, MockThermalGen, - ["dev1"], - 1:1, setup.var_container, [(min = 0.0, max = 4.0)], TEST_META,