Refactor approximations into a pure-JuMP layer

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)

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

src/approximations/common.jl

@@ -0,0 +1,34 @@
+# Shared infrastructure for quadratic and bilinear approximations.
+#
+# Each method ships a scalar `build_<method>` function (pure-JuMP, no IOM
+# dependencies) and an `add_<method>_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