From 375cbf8e36b8e7a29fbc8156a239f1e1b5b761af Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Tue, 3 Mar 2026 17:32:09 -0700 Subject: [PATCH 1/7] add bilinear models implementation --- src/quadratic_approximations/bilinear.jl | 278 +++++++++++++++++++++++ 1 file changed, 278 insertions(+) create mode 100644 src/quadratic_approximations/bilinear.jl diff --git a/src/quadratic_approximations/bilinear.jl b/src/quadratic_approximations/bilinear.jl new file mode 100644 index 0000000..273a2ba --- /dev/null +++ b/src/quadratic_approximations/bilinear.jl @@ -0,0 +1,278 @@ +# Bin2 separable approximation of bilinear products z = x·y. +# Uses the difference-of-squares identity: x·y = ¼((x+y)² − (x−y)²). +# Calls existing quadratic approximation functions twice for u² and v². + +struct BilinearApproxSumVariable <: VariableType end # u = x + y +struct BilinearApproxDiffVariable <: VariableType end # v = x − y +struct BilinearApproxSumLinkingConstraint <: ConstraintType end # u == x + y +struct BilinearApproxDiffLinkingConstraint <: ConstraintType end # v == x − y +struct BilinearProductVariable <: VariableType end # z ≈ x·y + +""" + _add_bilinear_approx_impl!(container, C, names, time_steps, x_var_container, y_var_container, x_min, x_max, y_min, y_max, quad_approx_fn, meta; add_mccormick) + +Internal implementation for Bin2 bilinear approximation using z = ¼((x+y)² − (x−y)²). + +Creates auxiliary variables u = x+y and v = x−y, calls `quad_approx_fn` twice to +approximate u² and v², then combines via the difference-of-squares identity. + +# Arguments +- `container::OptimizationContainer`: the optimization container +- `::Type{C}`: component type +- `names::Vector{String}`: component names +- `time_steps::UnitRange{Int}`: time periods +- `x_var_container`: container of x variables indexed by (name, t) +- `y_var_container`: container of y variables indexed by (name, t) +- `x_min::Float64`: lower bound of x +- `x_max::Float64`: upper bound of x +- `y_min::Float64`: lower bound of y +- `y_max::Float64`: upper bound of y +- `quad_approx_fn`: callable with signature (container, C, names, ts, var_cont, lo, hi, meta) → Dict +- `meta::String`: identifier for container keys +- `add_mccormick::Bool`: whether to add McCormick envelope constraints (default: false) + +# Returns +- `Dict{Tuple{String, Int}, JuMP.AffExpr}`: maps (name, t) to affine expression approximating x·y +""" +function _add_bilinear_approx_impl!( + container::OptimizationContainer, + ::Type{C}, + names::Vector{String}, + time_steps::UnitRange{Int}, + x_var_container, + y_var_container, + x_min::Float64, + x_max::Float64, + y_min::Float64, + y_max::Float64, + quad_approx_fn, + meta::String; + add_mccormick::Bool = false, +) where {C <: IS.InfrastructureSystemsComponent} + # Bounds for u = x + y and v = x − y + u_min = x_min + y_min + u_max = x_max + y_max + v_min = x_min - y_max + v_max = x_max - y_min + IS.@assert_op u_max > u_min + IS.@assert_op v_max > v_min + + jump_model = get_jump_model(container) + meta_plus = meta * "_plus" + meta_minus = meta * "_minus" + + # Create u and v variable containers + u_container = add_variable_container!( + container, + BilinearApproxSumVariable(), + C, + names, + time_steps; + meta = meta_plus, + ) + v_container = add_variable_container!( + container, + BilinearApproxDiffVariable(), + C, + names, + time_steps; + meta = meta_minus, + ) + + # Create linking constraint containers + u_link_container = add_constraints_container!( + container, + BilinearApproxSumLinkingConstraint(), + C, + names, + time_steps; + meta = meta_plus, + ) + v_link_container = add_constraints_container!( + container, + BilinearApproxDiffLinkingConstraint(), + C, + names, + time_steps; + meta = meta_minus, + ) + + # Create u, v variables and linking constraints + for name in names, t in time_steps + x = x_var_container[name, t] + y = y_var_container[name, t] + + u_container[name, t] = JuMP.@variable( + jump_model, + base_name = "BilinearSum_$(C)_{$(name), $(t)}", + lower_bound = u_min, + upper_bound = u_max, + ) + v_container[name, t] = JuMP.@variable( + jump_model, + base_name = "BilinearDiff_$(C)_{$(name), $(t)}", + lower_bound = v_min, + upper_bound = v_max, + ) + + u_link_container[name, t] = + JuMP.@constraint(jump_model, u_container[name, t] == x + y) + v_link_container[name, t] = + JuMP.@constraint(jump_model, v_container[name, t] == x - y) + end + + # Approximate u² and v² using the provided quadratic approximation function + zu_dict = quad_approx_fn(container, C, names, time_steps, u_container, u_min, u_max, meta_plus) + zv_dict = quad_approx_fn(container, C, names, time_steps, v_container, v_min, v_max, meta_minus) + + # Create z variable container for the bilinear product + z_container = add_variable_container!( + container, + BilinearProductVariable(), + C, + names, + time_steps; + meta, + ) + + result = Dict{Tuple{String, Int}, JuMP.AffExpr}() + + for name in names, t in time_steps + z_var = JuMP.@variable( + jump_model, + base_name = "BilinearProduct_$(C)_{$(name), $(t)}", + ) + z_container[name, t] = z_var + + # z = 0.25 * (u² − v²) + z_expr = JuMP.AffExpr(0.0) + JuMP.add_to_expression!(z_expr, 0.25, zu_dict[(name, t)]) + JuMP.add_to_expression!(z_expr, -0.25, zv_dict[(name, t)]) + + JuMP.@constraint(jump_model, z_var == z_expr) + + result[(name, t)] = JuMP.AffExpr(0.0, z_var => 1.0) + end + + # Optional McCormick envelope + if add_mccormick + _add_mccormick_envelope!( + container, C, names, time_steps, + x_var_container, y_var_container, z_container, + x_min, x_max, y_min, y_max, meta, + ) + end + + return result +end + +""" + _add_sos2_bilinear_approx!(container, C, names, time_steps, x_var_container, y_var_container, x_min, x_max, y_min, y_max, num_segments, meta; add_mccormick) + +Approximate x·y using Bin2 decomposition with solver-native SOS2 quadratic approximations. + +# Arguments +Same as `_add_bilinear_approx_impl!` plus: +- `num_segments::Int`: number of PWL segments for each quadratic approximation + +# Returns +- `Dict{Tuple{String, Int}, JuMP.AffExpr}`: maps (name, t) to affine expression approximating x·y +""" +function _add_sos2_bilinear_approx!( + container::OptimizationContainer, + ::Type{C}, + names::Vector{String}, + time_steps::UnitRange{Int}, + x_var_container, + y_var_container, + x_min::Float64, + x_max::Float64, + y_min::Float64, + y_max::Float64, + num_segments::Int, + meta::String; + add_mccormick::Bool = false, +) where {C <: IS.InfrastructureSystemsComponent} + quad_fn = (cont, CT, nms, ts, vc, lo, hi, m) -> + _add_sos2_quadratic_approx!(cont, CT, nms, ts, vc, lo, hi, num_segments, m) + return _add_bilinear_approx_impl!( + container, C, names, time_steps, + x_var_container, y_var_container, + x_min, x_max, y_min, y_max, quad_fn, meta; + add_mccormick, + ) +end + +""" + _add_manual_sos2_bilinear_approx!(container, C, names, time_steps, x_var_container, y_var_container, x_min, x_max, y_min, y_max, num_segments, meta; add_mccormick) + +Approximate x·y using Bin2 decomposition with manual SOS2 quadratic approximations. + +# Arguments +Same as `_add_bilinear_approx_impl!` plus: +- `num_segments::Int`: number of PWL segments for each quadratic approximation + +# Returns +- `Dict{Tuple{String, Int}, JuMP.AffExpr}`: maps (name, t) to affine expression approximating x·y +""" +function _add_manual_sos2_bilinear_approx!( + container::OptimizationContainer, + ::Type{C}, + names::Vector{String}, + time_steps::UnitRange{Int}, + x_var_container, + y_var_container, + x_min::Float64, + x_max::Float64, + y_min::Float64, + y_max::Float64, + num_segments::Int, + meta::String; + add_mccormick::Bool = false, +) where {C <: IS.InfrastructureSystemsComponent} + quad_fn = (cont, CT, nms, ts, vc, lo, hi, m) -> + _add_manual_sos2_quadratic_approx!(cont, CT, nms, ts, vc, lo, hi, num_segments, m) + return _add_bilinear_approx_impl!( + container, C, names, time_steps, + x_var_container, y_var_container, + x_min, x_max, y_min, y_max, quad_fn, meta; + add_mccormick, + ) +end + +""" + _add_sawtooth_bilinear_approx!(container, C, names, time_steps, x_var_container, y_var_container, x_min, x_max, y_min, y_max, depth, meta; add_mccormick) + +Approximate x·y using Bin2 decomposition with sawtooth quadratic approximations. + +# Arguments +Same as `_add_bilinear_approx_impl!` plus: +- `depth::Int`: sawtooth depth (number of binary variables per quadratic approximation) + +# Returns +- `Dict{Tuple{String, Int}, JuMP.AffExpr}`: maps (name, t) to affine expression approximating x·y +""" +function _add_sawtooth_bilinear_approx!( + container::OptimizationContainer, + ::Type{C}, + names::Vector{String}, + time_steps::UnitRange{Int}, + x_var_container, + y_var_container, + x_min::Float64, + x_max::Float64, + y_min::Float64, + y_max::Float64, + depth::Int, + meta::String; + add_mccormick::Bool = false, +) where {C <: IS.InfrastructureSystemsComponent} + quad_fn = (cont, CT, nms, ts, vc, lo, hi, m) -> + _add_sawtooth_quadratic_approx!(cont, CT, nms, ts, vc, lo, hi, depth, m) + return _add_bilinear_approx_impl!( + container, C, names, time_steps, + x_var_container, y_var_container, + x_min, x_max, y_min, y_max, quad_fn, meta; + add_mccormick, + ) +end From c171cf72d9cf794255f045d242d7178fa6cd2315 Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Tue, 3 Mar 2026 17:32:24 -0700 Subject: [PATCH 2/7] add a basic mccormick --- src/quadratic_approximations/mccormick.jl | 93 +++++++++++++++++++++++ 1 file changed, 93 insertions(+) create mode 100644 src/quadratic_approximations/mccormick.jl diff --git a/src/quadratic_approximations/mccormick.jl b/src/quadratic_approximations/mccormick.jl new file mode 100644 index 0000000..f1bdb96 --- /dev/null +++ b/src/quadratic_approximations/mccormick.jl @@ -0,0 +1,93 @@ +# McCormick envelope for bilinear products z = x·y. +# Adds 4 linear inequalities that bound z given variable bounds on x and y. + +struct McCormickConstraint <: ConstraintType end + +""" + _add_mccormick_envelope!(container, C, names, time_steps, x_var_container, y_var_container, z_var_container, 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`: container of x variables indexed by (name, t) +- `y_var_container`: container of y variables indexed by (name, t) +- `z_var_container`: 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_container, + y_var_container, + z_var_container, + x_min::Float64, + x_max::Float64, + y_min::Float64, + y_max::Float64, + meta::String, +) where {C <: IS.InfrastructureSystemsComponent} + IS.@assert_op x_max > x_min + IS.@assert_op y_max > y_min + jump_model = get_jump_model(container) + + mc_container = add_constraints_container!( + container, + McCormickConstraint(), + C, + names, + 1:4, + time_steps; + meta, + sparse = true, + ) + + for name in names, t in time_steps + x = x_var_container[name, t] + y = y_var_container[name, t] + z = z_var_container[name, t] + + # z ≥ x_min·y + x·y_min − x_min·y_min + mc_container[(name, 1, t)] = JuMP.@constraint( + jump_model, + z >= x_min * y + x * y_min - x_min * y_min, + ) + # z ≥ x_max·y + x·y_max − x_max·y_max + mc_container[(name, 2, t)] = JuMP.@constraint( + jump_model, + z >= x_max * y + x * y_max - x_max * y_max, + ) + # z ≤ x_max·y + x·y_min − x_max·y_min + mc_container[(name, 3, t)] = JuMP.@constraint( + jump_model, + z <= x_max * y + x * y_min - x_max * y_min, + ) + # z ≤ x_min·y + x·y_max − x_min·y_max + mc_container[(name, 4, t)] = JuMP.@constraint( + jump_model, + z <= x_min * y + x * y_max - x_min * y_max, + ) + end + + return nothing +end From d78befea4d113b4e8e2585da5231c0d8a338a05a Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Tue, 3 Mar 2026 17:32:33 -0700 Subject: [PATCH 3/7] add testing --- src/InfrastructureOptimizationModels.jl | 2 + test/InfrastructureOptimizationModelsTests.jl | 1 + test/test_bilinear_approximations.jl | 714 ++++++++++++++++++ 3 files changed, 717 insertions(+) create mode 100644 test/test_bilinear_approximations.jl diff --git a/src/InfrastructureOptimizationModels.jl b/src/InfrastructureOptimizationModels.jl index 27bb181..437b7ee 100644 --- a/src/InfrastructureOptimizationModels.jl +++ b/src/InfrastructureOptimizationModels.jl @@ -597,6 +597,8 @@ include("objective_function/market_bid.jl") include("quadratic_approximations/solver_sos2.jl") include("quadratic_approximations/manual_sos2.jl") include("quadratic_approximations/sawtooth.jl") +include("quadratic_approximations/mccormick.jl") +include("quadratic_approximations/bilinear.jl") # add_param_container! wrappers — must come after piecewise_linear.jl # (which defines VariableValueParameter and FixValueParameter) diff --git a/test/InfrastructureOptimizationModelsTests.jl b/test/InfrastructureOptimizationModelsTests.jl index 7c7e829..a93b161 100644 --- a/test/InfrastructureOptimizationModelsTests.jl +++ b/test/InfrastructureOptimizationModelsTests.jl @@ -142,6 +142,7 @@ function run_tests() # --- quadratic_approximations/ subfolder --- include(joinpath(TEST_DIR, "test_quadratic_approximations.jl")) + include(joinpath(TEST_DIR, "test_bilinear_approximations.jl")) end #= diff --git a/test/test_bilinear_approximations.jl b/test/test_bilinear_approximations.jl new file mode 100644 index 0000000..9ac6f61 --- /dev/null +++ b/test/test_bilinear_approximations.jl @@ -0,0 +1,714 @@ +const BILINEAR_META = "BilinearTest" + +function _setup_bilinear_test(device_names::Vector{String}, time_steps::UnitRange{Int}) + container = _setup_qa_container(time_steps) + x_var_container = IOM.add_variable_container!( + container, + TestOriginalVariable(), + MockThermalGen, + device_names, + time_steps, + ) + y_var_container = IOM.add_variable_container!( + container, + TestApproximatedVariable(), + MockThermalGen, + device_names, + time_steps, + ) + jump_model = IOM.get_jump_model(container) + for name in device_names, t in time_steps + x_var_container[name, t] = + JuMP.@variable(jump_model, base_name = "x_$(name)_$(t)") + y_var_container[name, t] = + JuMP.@variable(jump_model, base_name = "y_$(name)_$(t)") + end + return (; container, x_var_container, y_var_container, jump_model) +end + +@testset "Bilinear Approximations" begin + @testset "Solver SOS2 Bilinear" begin + @testset "Constraint structure" begin + setup = _setup_bilinear_test(["dev1"], 1:1) + + result = IOM._add_sos2_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:1, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + 4, + BILINEAR_META, + ) + + @test haskey(result, ("dev1", 1)) + @test result[("dev1", 1)] isa JuMP.AffExpr + + # u = x + y variable container should exist + @test IOM.has_container_key( + setup.container, + IOM.BilinearApproxSumVariable, + MockThermalGen, + BILINEAR_META * "_plus", + ) + # v = x - y variable container should exist + @test IOM.has_container_key( + setup.container, + IOM.BilinearApproxDiffVariable, + MockThermalGen, + BILINEAR_META * "_minus", + ) + # z ≈ x·y variable container should exist + @test IOM.has_container_key( + setup.container, + IOM.BilinearProductVariable, + MockThermalGen, + BILINEAR_META, + ) + # Linking constraints should exist + @test IOM.has_container_key( + setup.container, + IOM.BilinearApproxSumLinkingConstraint, + MockThermalGen, + BILINEAR_META * "_plus", + ) + @test IOM.has_container_key( + setup.container, + IOM.BilinearApproxDiffLinkingConstraint, + MockThermalGen, + BILINEAR_META * "_minus", + ) + # Inner quadratic approx containers should exist with _plus/_minus meta + @test IOM.has_container_key( + setup.container, + IOM.QuadraticApproxVariable, + MockThermalGen, + BILINEAR_META * "_plus", + ) + @test IOM.has_container_key( + setup.container, + IOM.QuadraticApproxVariable, + MockThermalGen, + BILINEAR_META * "_minus", + ) + + # u bounds should be [0+0, 4+4] = [0, 8] + u_container = IOM.get_variable( + setup.container, + IOM.BilinearApproxSumVariable(), + MockThermalGen, + BILINEAR_META * "_plus", + ) + @test JuMP.lower_bound(u_container["dev1", 1]) == 0.0 + @test JuMP.upper_bound(u_container["dev1", 1]) == 8.0 + + # v bounds should be [0-4, 4-0] = [-4, 4] + v_container = IOM.get_variable( + setup.container, + IOM.BilinearApproxDiffVariable(), + MockThermalGen, + BILINEAR_META * "_minus", + ) + @test JuMP.lower_bound(v_container["dev1", 1]) == -4.0 + @test JuMP.upper_bound(v_container["dev1", 1]) == 4.0 + end + + @testset "Constraint structure with McCormick" begin + setup = _setup_bilinear_test(["dev1"], 1:1) + + IOM._add_sos2_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:1, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + 4, + BILINEAR_META; + add_mccormick = true, + ) + + @test IOM.has_container_key( + setup.container, + IOM.McCormickConstraint, + MockThermalGen, + BILINEAR_META, + ) + end + + @testset "Fixed-variable correctness" begin + # Fix x=3, y ∈ [0,4]: min xy should give z≈0 at y=0 + setup = _setup_bilinear_test(["dev1"], 1:1) + x_var = setup.x_var_container["dev1", 1] + y_var = setup.y_var_container["dev1", 1] + JuMP.fix(x_var, 3.0; force = true) + JuMP.set_lower_bound(y_var, 0.0) + JuMP.set_upper_bound(y_var, 4.0) + + result = IOM._add_sos2_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:1, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + 8, + BILINEAR_META, + ) + z_expr = result[("dev1", 1)] + + JuMP.@objective(setup.jump_model, Min, z_expr) + JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) + JuMP.set_silent(setup.jump_model) + JuMP.optimize!(setup.jump_model) + + @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL + @test JuMP.objective_value(setup.jump_model) ≈ 0.0 atol = 1e-4 + + # Fix x=0, y=0: z should be exactly 0 + setup2 = _setup_bilinear_test(["dev1"], 1:1) + JuMP.fix(setup2.x_var_container["dev1", 1], 0.0; force = true) + JuMP.fix(setup2.y_var_container["dev1", 1], 0.0; force = true) + + result2 = IOM._add_sos2_bilinear_approx!( + setup2.container, + MockThermalGen, + ["dev1"], + 1:1, + setup2.x_var_container, + setup2.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + 4, + BILINEAR_META, + ) + z_expr2 = result2[("dev1", 1)] + + JuMP.@objective(setup2.jump_model, Max, z_expr2) + JuMP.set_optimizer(setup2.jump_model, HiGHS.Optimizer) + JuMP.set_silent(setup2.jump_model) + JuMP.optimize!(setup2.jump_model) + + @test JuMP.termination_status(setup2.jump_model) == JuMP.OPTIMAL + @test JuMP.objective_value(setup2.jump_model) ≈ 0.0 atol = 1e-6 + end + + @testset "Constraint usage: x·y + w = 10 with x=2" begin + setup = _setup_bilinear_test(["dev1"], 1:1) + x_var = setup.x_var_container["dev1", 1] + y_var = setup.y_var_container["dev1", 1] + JuMP.fix(x_var, 2.0; force = true) + JuMP.fix(y_var, 3.0; force = true) + + w = JuMP.@variable(setup.jump_model, base_name = "w") + + result = IOM._add_sos2_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:1, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + 8, + BILINEAR_META, + ) + z_expr = result[("dev1", 1)] + + # x·y + w = 10 → 2·3 + w = 10 → w = 4 + JuMP.@constraint(setup.jump_model, z_expr + w == 10.0) + JuMP.@objective(setup.jump_model, Min, w) + JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) + JuMP.set_silent(setup.jump_model) + JuMP.optimize!(setup.jump_model) + + @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL + @test JuMP.value(w) ≈ 4.0 atol = 1e-4 + end + + @testset "Vertex optimum" begin + # min x·y on [0,4]×[0,4]; minimum is 0 at a corner + setup = _setup_bilinear_test(["dev1"], 1:1) + x_var = setup.x_var_container["dev1", 1] + y_var = setup.y_var_container["dev1", 1] + JuMP.set_lower_bound(x_var, 0.0) + JuMP.set_upper_bound(x_var, 4.0) + JuMP.set_lower_bound(y_var, 0.0) + JuMP.set_upper_bound(y_var, 4.0) + + result = IOM._add_sos2_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:1, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + 8, + BILINEAR_META, + ) + z_expr = result[("dev1", 1)] + + JuMP.@objective(setup.jump_model, Min, z_expr) + JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) + JuMP.set_silent(setup.jump_model) + JuMP.optimize!(setup.jump_model) + + @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL + @test JuMP.objective_value(setup.jump_model) ≈ 0.0 atol = 1e-4 + end + + @testset "Multiple time steps" begin + setup = _setup_bilinear_test(["dev1"], 1:3) + result = IOM._add_sos2_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:3, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + 4, + BILINEAR_META, + ) + + for t in 1:3 + @test haskey(result, ("dev1", t)) + @test result[("dev1", t)] isa JuMP.AffExpr + end + end + + @testset "Approximation quality improves with segments" begin + # Fix x=2.5, y=1.5: true product = 3.75 + # Sweep segments, verify gap shrinks + true_product = 2.5 * 1.5 + errors = Float64[] + for num_segments in 2 .^ (1:5) + setup = _setup_bilinear_test(["dev1"], 1:1) + x_var = setup.x_var_container["dev1", 1] + y_var = setup.y_var_container["dev1", 1] + JuMP.fix(x_var, 2.5; force = true) + JuMP.fix(y_var, 1.5; force = true) + + result = IOM._add_sos2_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:1, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + num_segments, + BILINEAR_META, + ) + z_expr = result[("dev1", 1)] + + JuMP.@objective(setup.jump_model, Max, z_expr) + JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) + JuMP.set_silent(setup.jump_model) + JuMP.optimize!(setup.jump_model) + + obj_val = JuMP.objective_value(setup.jump_model) + push!(errors, abs(obj_val - true_product)) + end + for i in 2:length(errors) + @test errors[i] <= errors[i - 1] + 1e-10 + end + end + end + + @testset "Manual SOS2 Bilinear" begin + @testset "Constraint structure" begin + setup = _setup_bilinear_test(["dev1"], 1:1) + + result = IOM._add_manual_sos2_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:1, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + 4, + BILINEAR_META, + ) + + @test haskey(result, ("dev1", 1)) + @test result[("dev1", 1)] isa JuMP.AffExpr + + # Binary variables should exist for both u² and v² paths + @test IOM.has_container_key( + setup.container, + IOM.ManualSOS2BinaryVariable, + MockThermalGen, + BILINEAR_META * "_plus", + ) + @test IOM.has_container_key( + setup.container, + IOM.ManualSOS2BinaryVariable, + MockThermalGen, + BILINEAR_META * "_minus", + ) + + # NO solver SOS2 constraints + sos2_count = JuMP.num_constraints( + setup.jump_model, + Vector{JuMP.VariableRef}, + MOI.SOS2{Float64}, + ) + @test sos2_count == 0 + end + + @testset "Fixed-variable correctness" begin + setup = _setup_bilinear_test(["dev1"], 1:1) + JuMP.fix(setup.x_var_container["dev1", 1], 2.0; force = true) + JuMP.fix(setup.y_var_container["dev1", 1], 3.0; force = true) + + result = IOM._add_manual_sos2_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:1, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + 8, + BILINEAR_META, + ) + z_expr = result[("dev1", 1)] + + JuMP.@objective(setup.jump_model, Max, z_expr) + JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) + JuMP.set_silent(setup.jump_model) + JuMP.optimize!(setup.jump_model) + + @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL + @test JuMP.objective_value(setup.jump_model) ≈ 6.0 atol = 1e-4 + end + + @testset "Constraint usage: x·y + w = 10 with x=2, y=3" begin + setup = _setup_bilinear_test(["dev1"], 1:1) + JuMP.fix(setup.x_var_container["dev1", 1], 2.0; force = true) + JuMP.fix(setup.y_var_container["dev1", 1], 3.0; force = true) + + w = JuMP.@variable(setup.jump_model, base_name = "w") + + result = IOM._add_manual_sos2_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:1, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + 8, + BILINEAR_META, + ) + z_expr = result[("dev1", 1)] + + JuMP.@constraint(setup.jump_model, z_expr + w == 10.0) + JuMP.@objective(setup.jump_model, Min, w) + JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) + JuMP.set_silent(setup.jump_model) + JuMP.optimize!(setup.jump_model) + + @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL + @test JuMP.value(w) ≈ 4.0 atol = 1e-4 + end + + @testset "Vertex optimum" begin + setup = _setup_bilinear_test(["dev1"], 1:1) + x_var = setup.x_var_container["dev1", 1] + y_var = setup.y_var_container["dev1", 1] + JuMP.set_lower_bound(x_var, 0.0) + JuMP.set_upper_bound(x_var, 4.0) + JuMP.set_lower_bound(y_var, 0.0) + JuMP.set_upper_bound(y_var, 4.0) + + result = IOM._add_manual_sos2_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:1, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + 8, + BILINEAR_META, + ) + z_expr = result[("dev1", 1)] + + JuMP.@objective(setup.jump_model, Min, z_expr) + JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) + JuMP.set_silent(setup.jump_model) + JuMP.optimize!(setup.jump_model) + + @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL + @test JuMP.objective_value(setup.jump_model) ≈ 0.0 atol = 1e-4 + end + + @testset "Multiple time steps" begin + setup = _setup_bilinear_test(["dev1"], 1:3) + result = IOM._add_manual_sos2_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:3, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + 4, + BILINEAR_META, + ) + + for t in 1:3 + @test haskey(result, ("dev1", t)) + end + end + + @testset "Approximation quality improves with segments" begin + true_product = 2.5 * 1.5 + errors = Float64[] + for num_segments in 2 .^ (1:5) + setup = _setup_bilinear_test(["dev1"], 1:1) + JuMP.fix(setup.x_var_container["dev1", 1], 2.5; force = true) + JuMP.fix(setup.y_var_container["dev1", 1], 1.5; force = true) + + result = IOM._add_manual_sos2_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:1, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + num_segments, + BILINEAR_META, + ) + z_expr = result[("dev1", 1)] + + JuMP.@objective(setup.jump_model, Max, z_expr) + JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) + JuMP.set_silent(setup.jump_model) + JuMP.optimize!(setup.jump_model) + + obj_val = JuMP.objective_value(setup.jump_model) + push!(errors, abs(obj_val - true_product)) + end + for i in 2:length(errors) + @test errors[i] <= errors[i - 1] + 1e-10 + end + end + end + + @testset "Sawtooth Bilinear" begin + @testset "Constraint structure" begin + setup = _setup_bilinear_test(["dev1"], 1:1) + + result = IOM._add_sawtooth_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:1, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + 2, + BILINEAR_META, + ) + + @test haskey(result, ("dev1", 1)) + @test result[("dev1", 1)] isa JuMP.AffExpr + + # Sawtooth aux/binary variables for both u² and v² paths + @test IOM.has_container_key( + setup.container, + IOM.SawtoothAuxVariable, + MockThermalGen, + BILINEAR_META * "_plus", + ) + @test IOM.has_container_key( + setup.container, + IOM.SawtoothBinaryVariable, + MockThermalGen, + BILINEAR_META * "_plus", + ) + @test IOM.has_container_key( + setup.container, + IOM.SawtoothAuxVariable, + MockThermalGen, + BILINEAR_META * "_minus", + ) + @test IOM.has_container_key( + setup.container, + IOM.SawtoothBinaryVariable, + MockThermalGen, + BILINEAR_META * "_minus", + ) + end + + @testset "Fixed-variable correctness" begin + setup = _setup_bilinear_test(["dev1"], 1:1) + JuMP.fix(setup.x_var_container["dev1", 1], 2.0; force = true) + JuMP.fix(setup.y_var_container["dev1", 1], 3.0; force = true) + + result = IOM._add_sawtooth_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:1, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + 3, + BILINEAR_META, + ) + z_expr = result[("dev1", 1)] + + JuMP.@objective(setup.jump_model, Max, z_expr) + JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) + JuMP.set_silent(setup.jump_model) + JuMP.optimize!(setup.jump_model) + + @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL + @test JuMP.objective_value(setup.jump_model) ≈ 6.0 atol = 1e-3 + end + + @testset "Constraint usage: x·y + w = 10 with x=2, y=3" begin + setup = _setup_bilinear_test(["dev1"], 1:1) + JuMP.fix(setup.x_var_container["dev1", 1], 2.0; force = true) + JuMP.fix(setup.y_var_container["dev1", 1], 3.0; force = true) + + w = JuMP.@variable(setup.jump_model, base_name = "w") + + result = IOM._add_sawtooth_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:1, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + 3, + BILINEAR_META, + ) + z_expr = result[("dev1", 1)] + + JuMP.@constraint(setup.jump_model, z_expr + w == 10.0) + JuMP.@objective(setup.jump_model, Min, w) + JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) + JuMP.set_silent(setup.jump_model) + JuMP.optimize!(setup.jump_model) + + @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL + @test JuMP.value(w) ≈ 4.0 atol = 1e-3 + end + + @testset "Vertex optimum" begin + setup = _setup_bilinear_test(["dev1"], 1:1) + x_var = setup.x_var_container["dev1", 1] + y_var = setup.y_var_container["dev1", 1] + JuMP.set_lower_bound(x_var, 0.0) + JuMP.set_upper_bound(x_var, 4.0) + JuMP.set_lower_bound(y_var, 0.0) + JuMP.set_upper_bound(y_var, 4.0) + + result = IOM._add_sawtooth_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:1, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + 3, + BILINEAR_META, + ) + z_expr = result[("dev1", 1)] + + JuMP.@objective(setup.jump_model, Min, z_expr) + JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) + JuMP.set_silent(setup.jump_model) + JuMP.optimize!(setup.jump_model) + + @test JuMP.termination_status(setup.jump_model) == JuMP.OPTIMAL + @test JuMP.objective_value(setup.jump_model) ≈ 0.0 atol = 1e-3 + end + + @testset "Multiple time steps" begin + setup = _setup_bilinear_test(["dev1"], 1:3) + result = IOM._add_sawtooth_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:3, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + 2, + BILINEAR_META, + ) + + for t in 1:3 + @test haskey(result, ("dev1", t)) + end + end + + @testset "Approximation quality improves with depth" begin + true_product = 2.5 * 1.5 + errors = Float64[] + for depth in 1:5 + setup = _setup_bilinear_test(["dev1"], 1:1) + JuMP.fix(setup.x_var_container["dev1", 1], 2.5; force = true) + JuMP.fix(setup.y_var_container["dev1", 1], 1.5; force = true) + + result = IOM._add_sawtooth_bilinear_approx!( + setup.container, + MockThermalGen, + ["dev1"], + 1:1, + setup.x_var_container, + setup.y_var_container, + 0.0, 4.0, + 0.0, 4.0, + depth, + BILINEAR_META, + ) + z_expr = result[("dev1", 1)] + + JuMP.@objective(setup.jump_model, Max, z_expr) + JuMP.set_optimizer(setup.jump_model, HiGHS.Optimizer) + JuMP.set_silent(setup.jump_model) + JuMP.optimize!(setup.jump_model) + + obj_val = JuMP.objective_value(setup.jump_model) + push!(errors, abs(obj_val - true_product)) + end + for i in 2:length(errors) + @test errors[i] <= errors[i - 1] + 1e-10 + end + end + end +end From 60d6ec457fdfa05e8d943b0afe8df5f6ac160a03 Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Wed, 4 Mar 2026 08:27:37 -0700 Subject: [PATCH 4/7] add benchmark auto generated code --- scripts/bilinear_delta_benchmark.jl | 476 ++++++++++++++++++++++++++++ 1 file changed, 476 insertions(+) create mode 100644 scripts/bilinear_delta_benchmark.jl diff --git a/scripts/bilinear_delta_benchmark.jl b/scripts/bilinear_delta_benchmark.jl new file mode 100644 index 0000000..05022c4 --- /dev/null +++ b/scripts/bilinear_delta_benchmark.jl @@ -0,0 +1,476 @@ +""" +Benchmark script for the bilinear delta model described in `bilinear_delta_model.tex`. + +Builds a network OPF with delta (incremental) PWL costs and bilinear V·I power +balance constraints. Compares three Bin2 bilinear approximation methods from +InfrastructureOptimizationModels: + + 1. Solver-native SOS2 + 2. Manual SOS2 (binary-variable adjacency) + 3. Sawtooth MIP + +Usage: + julia --project=test scripts/bilinear_delta_benchmark.jl [N] [K] [seed] + + N = number of nodes (default 10, must be even) + K = number of PWL cost segments per generator (default 3) + seed = random seed for network generation (default 42) +""" + +using InfrastructureOptimizationModels +using JuMP +using HiGHS +using Dates +using Random +using Printf + +const IOM = InfrastructureOptimizationModels +const IS = IOM.IS + +# Try to load Ipopt for NLP reference; skip if unavailable +const HAS_IPOPT = try + @eval using Ipopt + true +catch + false +end + +# ─── Minimal types for IOM container infrastructure ────────────────────────── + +mutable struct BenchmarkSystem <: IS.InfrastructureSystemsContainer + base_power::Float64 +end + +IOM.get_base_power(sys::BenchmarkSystem) = sys.base_power +IOM.stores_time_series_in_memory(::BenchmarkSystem) = false + +struct NetworkNode <: IS.InfrastructureSystemsComponent end + +struct VoltageVariable <: IOM.VariableType end +struct CurrentVariable <: IOM.VariableType end + +# ─── Network data ──────────────────────────────────────────────────────────── + +struct NetworkData + N::Int + K::Int + gen_nodes::Vector{String} + dem_nodes::Vector{String} + all_nodes::Vector{String} + edges::Vector{Tuple{String, String}} + conductances::Dict{Tuple{String, String}, Float64} + demands::Dict{String, Float64} + marginal_costs::Dict{String, Vector{Float64}} + segment_widths::Dict{String, Vector{Float64}} +end + +""" +Generate a random connected network matching the model in bilinear_delta_model.tex. + +- `N` nodes split evenly into generators and demands. +- `K` PWL cost segments per generator with nondecreasing marginal costs. +- Connected graph via random spanning tree plus extra edges. +- Conductances scaled for feasibility given voltage/current bounds. +- Demands d_i ∈ [0.05, 0.15]. +""" +function generate_network(; N::Int = 10, K::Int = 3, seed::Int = 42) + @assert iseven(N) "N must be even" + rng = MersenneTwister(seed) + + all_nodes = ["n$(i)" for i in 1:N] + gen_nodes = all_nodes[1:(N ÷ 2)] + dem_nodes = all_nodes[(N ÷ 2 + 1):N] + + edges = Tuple{String, String}[] + conductances = Dict{Tuple{String, String}, Float64}() + + # Random spanning tree: shuffle nodes, connect consecutive pairs + perm = shuffle(rng, 1:N) + for idx in 1:(N - 1) + a, b = all_nodes[perm[idx]], all_nodes[perm[idx + 1]] + e = a < b ? (a, b) : (b, a) + if e ∉ edges + push!(edges, e) + conductances[e] = 1.0 + 4.0 * rand(rng) + end + end + + # Extra edges for density + for _ in 1:(N ÷ 2) + i, j = rand(rng, 1:N), rand(rng, 1:N) + i == j && continue + a, b = all_nodes[i], all_nodes[j] + e = a < b ? (a, b) : (b, a) + if e ∉ edges + push!(edges, e) + conductances[e] = 1.0 + 4.0 * rand(rng) + end + end + + # Demands at demand nodes (small enough for network to transport) + demands = Dict(d => 0.05 + 0.1 * rand(rng) for d in dem_nodes) + + # PWL cost data: K segments per generator with nondecreasing marginal costs + marginal_costs = Dict{String, Vector{Float64}}() + segment_widths = Dict{String, Vector{Float64}}() + for g in gen_nodes + mc = sort(rand(rng, K) .* 10.0 .+ 1.0) + marginal_costs[g] = mc + widths = rand(rng, K) .+ 0.1 + widths .*= 1.5 / sum(widths) + segment_widths[g] = widths + end + + return NetworkData( + N, K, gen_nodes, dem_nodes, all_nodes, + edges, conductances, demands, + marginal_costs, segment_widths, + ) +end + +"""Build adjacency list from edge set.""" +function adjacency_list(net::NetworkData) + adj = Dict{String, Vector{Tuple{String, Float64}}}() + for n in net.all_nodes + adj[n] = Tuple{String, Float64}[] + end + for (a, b) in net.edges + g = net.conductances[(a, b)] + push!(adj[a], (b, g)) + push!(adj[b], (a, g)) + end + return adj +end + +# ─── Model constants ───────────────────────────────────────────────────────── + +const V_MIN = 0.8 +const V_MAX = 1.2 +const I_GEN_MIN = 0.0 +const I_GEN_MAX = 1.0 +const I_DEM_MIN = -1.0 +const I_DEM_MAX = 0.0 +const P_MAX = 1.5 + +# ─── NLP reference model (direct bilinear constraints via Ipopt) ───────────── + +function build_nlp_model(net::NetworkData) + model = JuMP.Model() + adj = adjacency_list(net) + + # Variables + @variable(model, V_MIN <= V[n in net.all_nodes] <= V_MAX) + @variable(model, I_GEN_MIN <= Ig[g in net.gen_nodes] <= I_GEN_MAX) + @variable(model, I_DEM_MIN <= Id[d in net.dem_nodes] <= I_DEM_MAX) + @variable(model, 0.0 <= Pg[g in net.gen_nodes] <= P_MAX) + @variable(model, + 0.0 <= delta[g in net.gen_nodes, k in 1:(net.K)] <= net.segment_widths[g][k]) + + # Objective: min Σ m_{i,k} · δ_{i,k} + @objective(model, Min, sum( + net.marginal_costs[g][k] * delta[g, k] + for g in net.gen_nodes, k in 1:(net.K) + )) + + # Bilinear: Pg = V · I (generators) + @constraint(model, [g in net.gen_nodes], Pg[g] == V[g] * Ig[g]) + + # Bilinear: V · I = -d (demands) + @constraint(model, [d in net.dem_nodes], V[d] * Id[d] == -net.demands[d]) + + # KCL + for g in net.gen_nodes + @constraint(model, Ig[g] == sum(c * (V[g] - V[j]) for (j, c) in adj[g])) + end + for d in net.dem_nodes + @constraint(model, Id[d] == sum(c * (V[d] - V[j]) for (j, c) in adj[d])) + end + + # Delta link: Pg = Σ δ_{g,k} + @constraint(model, [g in net.gen_nodes], + Pg[g] == sum(delta[g, k] for k in 1:(net.K))) + + return model +end + +# ─── IOM container setup ───────────────────────────────────────────────────── + +function make_container() + sys = BenchmarkSystem(100.0) + settings = IOM.Settings(sys; horizon = Dates.Hour(1), resolution = Dates.Hour(1)) + container = IOM.OptimizationContainer(sys, settings, JuMP.Model(), IS.Deterministic) + IOM.set_time_steps!(container, 1:1) + return container +end + +# ─── MIP model using IOM bilinear approximations ───────────────────────────── + +""" + build_mip_model(net, method, refinement) -> NamedTuple + +Build the MIP approximation of the bilinear delta model. + +# Arguments +- `net::NetworkData`: network data +- `method::Symbol`: `:sos2`, `:manual_sos2`, or `:sawtooth` +- `refinement::Int`: number of PWL segments (SOS2) or sawtooth depth +""" +function build_mip_model(net::NetworkData, method::Symbol, refinement::Int) + container = make_container() + jump_model = IOM.get_jump_model(container) + time_steps = 1:1 + adj = adjacency_list(net) + + # --- V and I variable containers via IOM --- + V_container = IOM.add_variable_container!( + container, VoltageVariable(), NetworkNode, + net.all_nodes, time_steps, + ) + I_container = IOM.add_variable_container!( + container, CurrentVariable(), NetworkNode, + net.all_nodes, time_steps, + ) + + for n in net.all_nodes + V_container[n, 1] = JuMP.@variable( + jump_model, base_name = "V_$(n)", + lower_bound = V_MIN, upper_bound = V_MAX, + ) + end + for g in net.gen_nodes + I_container[g, 1] = JuMP.@variable( + jump_model, base_name = "I_$(g)", + lower_bound = I_GEN_MIN, upper_bound = I_GEN_MAX, + ) + end + for d in net.dem_nodes + I_container[d, 1] = JuMP.@variable( + jump_model, base_name = "I_$(d)", + lower_bound = I_DEM_MIN, upper_bound = I_DEM_MAX, + ) + end + + # --- Dispatch to the chosen bilinear approximation method --- + bilinear_fn! = if method === :sos2 + (cont, names, xc, yc, ylo, yhi, meta) -> + IOM._add_sos2_bilinear_approx!( + cont, NetworkNode, names, time_steps, xc, yc, + V_MIN, V_MAX, ylo, yhi, refinement, meta, + ) + elseif method === :manual_sos2 + (cont, names, xc, yc, ylo, yhi, meta) -> + IOM._add_manual_sos2_bilinear_approx!( + cont, NetworkNode, names, time_steps, xc, yc, + V_MIN, V_MAX, ylo, yhi, refinement, meta, + ) + elseif method === :sawtooth + (cont, names, xc, yc, ylo, yhi, meta) -> + IOM._add_sawtooth_bilinear_approx!( + cont, NetworkNode, names, time_steps, xc, yc, + V_MIN, V_MAX, ylo, yhi, refinement, meta, + ) + else + error("Unknown method: $method. Use :sos2, :manual_sos2, or :sawtooth.") + end + + # Generator bilinear: z_gen ≈ V · I_gen + z_gen = bilinear_fn!( + container, net.gen_nodes, V_container, I_container, + I_GEN_MIN, I_GEN_MAX, "gen", + ) + + # Demand bilinear: z_dem ≈ V · I_dem + z_dem = bilinear_fn!( + container, net.dem_nodes, V_container, I_container, + I_DEM_MIN, I_DEM_MAX, "dem", + ) + + # --- Remaining linear model components (directly in JuMP) --- + Pg = Dict{String, JuMP.VariableRef}() + delta = Dict{Tuple{String, Int}, JuMP.VariableRef}() + + for g in net.gen_nodes + Pg[g] = JuMP.@variable(jump_model, base_name = "Pg_$(g)", + lower_bound = 0.0, upper_bound = P_MAX) + for k in 1:(net.K) + delta[(g, k)] = JuMP.@variable(jump_model, base_name = "delta_$(g)_$(k)", + lower_bound = 0.0, upper_bound = net.segment_widths[g][k]) + end + end + + # Objective: min Σ m_{i,k} · δ_{i,k} + JuMP.@objective(jump_model, Min, sum( + net.marginal_costs[g][k] * delta[(g, k)] + for g in net.gen_nodes, k in 1:(net.K) + )) + + # Pg == bilinear approx of V·I (generators) + for g in net.gen_nodes + JuMP.@constraint(jump_model, Pg[g] == z_gen[(g, 1)]) + end + + # V·I == -d (demands) + for d in net.dem_nodes + JuMP.@constraint(jump_model, z_dem[(d, 1)] == -net.demands[d]) + end + + # KCL: I_i = Σ g_{ij}(V_i - V_j) + for n in net.all_nodes + JuMP.@constraint(jump_model, + I_container[n, 1] == sum( + c * (V_container[n, 1] - V_container[j, 1]) + for (j, c) in adj[n] + )) + end + + # Delta link: Pg = Σ δ_{g,k} + for g in net.gen_nodes + JuMP.@constraint(jump_model, + Pg[g] == sum(delta[(g, k)] for k in 1:(net.K))) + end + + return (; container, jump_model, V_container, I_container, z_gen, z_dem) +end + +# ─── Metrics ────────────────────────────────────────────────────────────────── + +function compute_bilinear_residuals(result, net) + V = result.V_container + I = result.I_container + residuals = Float64[] + for g in net.gen_nodes + v_val = JuMP.value(V[g, 1]) + i_val = JuMP.value(I[g, 1]) + z_val = JuMP.value(result.z_gen[(g, 1)]) + push!(residuals, abs(v_val * i_val - z_val)) + end + for d in net.dem_nodes + v_val = JuMP.value(V[d, 1]) + i_val = JuMP.value(I[d, 1]) + z_val = JuMP.value(result.z_dem[(d, 1)]) + push!(residuals, abs(v_val * i_val - z_val)) + end + return residuals +end + +function model_size(jump_model) + nv = JuMP.num_variables(jump_model) + nc = sum( + JuMP.num_constraints(jump_model, f, s) + for (f, s) in JuMP.list_of_constraint_types(jump_model) + ) + nb = JuMP.num_constraints(jump_model, JuMP.VariableRef, JuMP.MOI.ZeroOne) + return (; variables = nv, constraints = nc, binaries = nb) +end + +# ─── Main benchmark ────────────────────────────────────────────────────────── + +""" + run_benchmark(; N, K, seed, segment_counts, sawtooth_depths) + +Run the full benchmark comparing the three Bin2 bilinear approximation methods. + +# Keyword arguments +- `N::Int = 10`: number of nodes (must be even) +- `K::Int = 3`: number of PWL cost segments per generator +- `seed::Int = 42`: random seed for network generation +- `segment_counts::Vector{Int} = [2, 4, 8, 16]`: segment counts for SOS2 methods +- `sawtooth_depths::Vector{Int} = [2, 4, 8]`: depths for sawtooth method +""" +function run_benchmark(; + N::Int = 10, + K::Int = 3, + seed::Int = 42, + segment_counts::Vector{Int} = [2, 4, 8, 16], + sawtooth_depths::Vector{Int} = [2, 4, 8], +) + net = generate_network(; N, K, seed) + println("Network: $(net.N) nodes, $(length(net.edges)) edges, $(net.K) cost segments") + println("Generators: $(length(net.gen_nodes)), Demands: $(length(net.dem_nodes))") + println() + + # --- NLP reference --- + nlp_obj = NaN + if HAS_IPOPT + println("=" ^ 90) + println("NLP Reference (Ipopt, bilinear V·I constraints)") + println("=" ^ 90) + nlp_model = build_nlp_model(net) + JuMP.set_optimizer(nlp_model, Ipopt.Optimizer) + JuMP.set_optimizer_attribute(nlp_model, "print_level", 0) + nlp_t = @elapsed JuMP.optimize!(nlp_model) + nlp_status = JuMP.termination_status(nlp_model) + sz = model_size(nlp_model) + if nlp_status == JuMP.LOCALLY_SOLVED + nlp_obj = JuMP.objective_value(nlp_model) + end + @printf(" Status: %s\n", nlp_status) + @printf(" Objective: %.6f\n", nlp_obj) + @printf(" Variables: %d\n", sz.variables) + @printf(" Constraints: %d\n", sz.constraints) + @printf(" Solve time: %.4f s\n", nlp_t) + else + println("Ipopt not available — skipping NLP reference.") + println("Install Ipopt.jl in the test environment for NLP comparison.") + end + println() + + # --- MIP benchmarks --- + method_configs = [ + (:sos2, "Solver SOS2", segment_counts), + (:manual_sos2, "Manual SOS2", segment_counts), + (:sawtooth, "Sawtooth", sawtooth_depths), + ] + + println("=" ^ 90) + println("MIP Bilinear Approximations (HiGHS)") + println(" Refinement = num_segments for SOS2 methods, depth for Sawtooth") + println("=" ^ 90) + @printf("%-14s %4s %6s %7s %6s %12s %9s %10s %8s\n", + "Method", "Ref", "Vars", "Constrs", "Bins", "Objective", "Gap(%)", "Max Resid", "Time(s)") + println("-" ^ 90) + + for (method, label, refs) in method_configs + for ref in refs + build_t = @elapsed begin + result = build_mip_model(net, method, ref) + end + + JuMP.set_optimizer(result.jump_model, HiGHS.Optimizer) + JuMP.set_optimizer_attribute(result.jump_model, "log_to_console", false) + JuMP.set_optimizer_attribute(result.jump_model, "time_limit", 300.0) + + solve_t = @elapsed JuMP.optimize!(result.jump_model) + status = JuMP.termination_status(result.jump_model) + sz = model_size(result.jump_model) + + if status in (JuMP.OPTIMAL, JuMP.TIME_LIMIT) && + JuMP.has_values(result.jump_model) + obj = JuMP.objective_value(result.jump_model) + gap = isnan(nlp_obj) ? NaN : abs(nlp_obj - obj) / max(abs(nlp_obj), 1e-10) * 100.0 + resids = compute_bilinear_residuals(result, net) + gap_str = isnan(gap) ? " -" : @sprintf("%8.4f", gap) + @printf("%-14s %4d %6d %7d %6d %12.6f %s %10.2e %8.4f\n", + label, ref, + sz.variables, sz.constraints, sz.binaries, + obj, gap_str, maximum(resids), solve_t) + else + @printf("%-14s %4d %6d %7d %6d %12s %9s %10s %8.4f\n", + label, ref, + sz.variables, sz.constraints, sz.binaries, + string(status), "-", "-", solve_t) + end + end + println() + end + println("=" ^ 90) +end + +# ─── Entry point ────────────────────────────────────────────────────────────── + +if abspath(PROGRAM_FILE) == @__FILE__ + N = get(ARGS, 1, "10") |> x -> parse(Int, x) + K = get(ARGS, 2, "3") |> x -> parse(Int, x) + seed = get(ARGS, 3, "42") |> x -> parse(Int, x) + run_benchmark(; N, K, seed) +end From 1d6f98d22c12da47925a3017b5aa72abfe8b6005 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Wed, 4 Mar 2026 13:35:20 -0500 Subject: [PATCH 5/7] replaced bin1 implementation with bin2, renamed u -> p aligned with literature --- src/quadratic_approximations/bilinear.jl | 110 +++++++++++------------ 1 file changed, 54 insertions(+), 56 deletions(-) diff --git a/src/quadratic_approximations/bilinear.jl b/src/quadratic_approximations/bilinear.jl index 273a2ba..b5ed838 100644 --- a/src/quadratic_approximations/bilinear.jl +++ b/src/quadratic_approximations/bilinear.jl @@ -11,10 +11,10 @@ struct BilinearProductVariable <: VariableType end # z ≈ x·y """ _add_bilinear_approx_impl!(container, C, names, time_steps, x_var_container, y_var_container, x_min, x_max, y_min, y_max, quad_approx_fn, meta; add_mccormick) -Internal implementation for Bin2 bilinear approximation using z = ¼((x+y)² − (x−y)²). +Internal implementation for Bin2 bilinear approximation using z = (1/2)((x+y)² − x² - y²). -Creates auxiliary variables u = x+y and v = x−y, calls `quad_approx_fn` twice to -approximate u² and v², then combines via the difference-of-squares identity. +Creates auxiliary variables p = x+y, calls `quad_approx_fn` to +approximate p², then combines via multiplicative identity. # Arguments - `container::OptimizationContainer`: the optimization container @@ -49,20 +49,16 @@ function _add_bilinear_approx_impl!( meta::String; add_mccormick::Bool = false, ) where {C <: IS.InfrastructureSystemsComponent} - # Bounds for u = x + y and v = x − y - u_min = x_min + y_min - u_max = x_max + y_max - v_min = x_min - y_max - v_max = x_max - y_min - IS.@assert_op u_max > u_min - IS.@assert_op v_max > v_min + # Bounds for p = x + y + p_min = x_min + y_min + p_max = x_max + y_max + IS.@assert_op p_min <= p_max jump_model = get_jump_model(container) meta_plus = meta * "_plus" - meta_minus = meta * "_minus" - # Create u and v variable containers - u_container = add_variable_container!( + # Create p variable container + p_container = add_variable_container!( container, BilinearApproxSumVariable(), C, @@ -70,17 +66,9 @@ function _add_bilinear_approx_impl!( time_steps; meta = meta_plus, ) - v_container = add_variable_container!( - container, - BilinearApproxDiffVariable(), - C, - names, - time_steps; - meta = meta_minus, - ) # Create linking constraint containers - u_link_container = add_constraints_container!( + p_link_container = add_constraints_container!( container, BilinearApproxSumLinkingConstraint(), C, @@ -88,42 +76,38 @@ function _add_bilinear_approx_impl!( time_steps; meta = meta_plus, ) - v_link_container = add_constraints_container!( - container, - BilinearApproxDiffLinkingConstraint(), - C, - names, - time_steps; - meta = meta_minus, - ) - # Create u, v variables and linking constraints + # Create p variable and linking constraint for name in names, t in time_steps x = x_var_container[name, t] y = y_var_container[name, t] - u_container[name, t] = JuMP.@variable( + p_container[name, t] = JuMP.@variable( jump_model, base_name = "BilinearSum_$(C)_{$(name), $(t)}", - lower_bound = u_min, - upper_bound = u_max, - ) - v_container[name, t] = JuMP.@variable( - jump_model, - base_name = "BilinearDiff_$(C)_{$(name), $(t)}", - lower_bound = v_min, - upper_bound = v_max, + lower_bound = p_min, + upper_bound = p_max, ) - u_link_container[name, t] = - JuMP.@constraint(jump_model, u_container[name, t] == x + y) - v_link_container[name, t] = - JuMP.@constraint(jump_model, v_container[name, t] == x - y) + p_link_container[name, t] = + JuMP.@constraint(jump_model, p_container[name, t] == x + y) end - # Approximate u² and v² using the provided quadratic approximation function - zu_dict = quad_approx_fn(container, C, names, time_steps, u_container, u_min, u_max, meta_plus) - zv_dict = quad_approx_fn(container, C, names, time_steps, v_container, v_min, v_max, meta_minus) + # Approximate p² using the provided quadratic approximation function + zp_dict = quad_approx_fn( + container, + C, + names, + time_steps, + p_container, + p_min, + p_max, + meta_plus, + ) + zx_dict = + quad_approx_fn(container, C, names, time_steps, x_var_container, x_min, x_max, meta) + zy_dict = + quad_approx_fn(container, C, names, time_steps, y_var_container, y_min, y_max, meta) # Create z variable container for the bilinear product z_container = add_variable_container!( @@ -144,10 +128,11 @@ function _add_bilinear_approx_impl!( ) z_container[name, t] = z_var - # z = 0.25 * (u² − v²) + # z = (1/2) * (p² − x² - y²) z_expr = JuMP.AffExpr(0.0) - JuMP.add_to_expression!(z_expr, 0.25, zu_dict[(name, t)]) - JuMP.add_to_expression!(z_expr, -0.25, zv_dict[(name, t)]) + JuMP.add_to_expression!(z_expr, 0.5, zp_dict[(name, t)]) + JuMP.add_to_expression!(z_expr, -0.5, zx_dict[(name, t)]) + JuMP.add_to_expression!(z_expr, -0.5, zy_dict[(name, t)]) JuMP.@constraint(jump_model, z_var == z_expr) @@ -193,8 +178,9 @@ function _add_sos2_bilinear_approx!( meta::String; add_mccormick::Bool = false, ) where {C <: IS.InfrastructureSystemsComponent} - quad_fn = (cont, CT, nms, ts, vc, lo, hi, m) -> - _add_sos2_quadratic_approx!(cont, CT, nms, ts, vc, lo, hi, num_segments, m) + quad_fn = + (cont, CT, nms, ts, vc, lo, hi, m) -> + _add_sos2_quadratic_approx!(cont, CT, nms, ts, vc, lo, hi, num_segments, m) return _add_bilinear_approx_impl!( container, C, names, time_steps, x_var_container, y_var_container, @@ -230,8 +216,19 @@ function _add_manual_sos2_bilinear_approx!( meta::String; add_mccormick::Bool = false, ) where {C <: IS.InfrastructureSystemsComponent} - quad_fn = (cont, CT, nms, ts, vc, lo, hi, m) -> - _add_manual_sos2_quadratic_approx!(cont, CT, nms, ts, vc, lo, hi, num_segments, m) + quad_fn = + (cont, CT, nms, ts, vc, lo, hi, m) -> + _add_manual_sos2_quadratic_approx!( + cont, + CT, + nms, + ts, + vc, + lo, + hi, + num_segments, + m, + ) return _add_bilinear_approx_impl!( container, C, names, time_steps, x_var_container, y_var_container, @@ -267,8 +264,9 @@ function _add_sawtooth_bilinear_approx!( meta::String; add_mccormick::Bool = false, ) where {C <: IS.InfrastructureSystemsComponent} - quad_fn = (cont, CT, nms, ts, vc, lo, hi, m) -> - _add_sawtooth_quadratic_approx!(cont, CT, nms, ts, vc, lo, hi, depth, m) + quad_fn = + (cont, CT, nms, ts, vc, lo, hi, m) -> + _add_sawtooth_quadratic_approx!(cont, CT, nms, ts, vc, lo, hi, depth, m) return _add_bilinear_approx_impl!( container, C, names, time_steps, x_var_container, y_var_container, From 1ef8cc6011199af9fd5ee8414de5490b54944837 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Thu, 5 Mar 2026 13:04:14 -0500 Subject: [PATCH 6/7] added unique meta keys for x and y variable containers; removed _minus from tests; tests now pass --- src/quadratic_approximations/bilinear.jl | 34 ++++++++++++---- test/test_bilinear_approximations.jl | 51 +----------------------- 2 files changed, 28 insertions(+), 57 deletions(-) diff --git a/src/quadratic_approximations/bilinear.jl b/src/quadratic_approximations/bilinear.jl index b5ed838..0d2866d 100644 --- a/src/quadratic_approximations/bilinear.jl +++ b/src/quadratic_approximations/bilinear.jl @@ -1,11 +1,9 @@ # Bin2 separable approximation of bilinear products z = x·y. -# Uses the difference-of-squares identity: x·y = ¼((x+y)² − (x−y)²). -# Calls existing quadratic approximation functions twice for u² and v². +# Uses the identity: x·y = (1/2)*((x+y)² − x² - y²). +# Calls existing quadratic approximation functions for p²=(x+y)² -struct BilinearApproxSumVariable <: VariableType end # u = x + y -struct BilinearApproxDiffVariable <: VariableType end # v = x − y -struct BilinearApproxSumLinkingConstraint <: ConstraintType end # u == x + y -struct BilinearApproxDiffLinkingConstraint <: ConstraintType end # v == x − y +struct BilinearApproxSumVariable <: VariableType end # p = x + y +struct BilinearApproxSumLinkingConstraint <: ConstraintType end # p == x + y struct BilinearProductVariable <: VariableType end # z ≈ x·y """ @@ -56,6 +54,8 @@ function _add_bilinear_approx_impl!( jump_model = get_jump_model(container) meta_plus = meta * "_plus" + meta_x = meta * "_x" + meta_y = meta * "_y" # Create p variable container p_container = add_variable_container!( @@ -105,9 +105,27 @@ function _add_bilinear_approx_impl!( meta_plus, ) zx_dict = - quad_approx_fn(container, C, names, time_steps, x_var_container, x_min, x_max, meta) + quad_approx_fn( + container, + C, + names, + time_steps, + x_var_container, + x_min, + x_max, + meta_x, + ) zy_dict = - quad_approx_fn(container, C, names, time_steps, y_var_container, y_min, y_max, meta) + quad_approx_fn( + container, + C, + names, + time_steps, + y_var_container, + y_min, + y_max, + meta_y, + ) # Create z variable container for the bilinear product z_container = add_variable_container!( diff --git a/test/test_bilinear_approximations.jl b/test/test_bilinear_approximations.jl index 9ac6f61..6940ec6 100644 --- a/test/test_bilinear_approximations.jl +++ b/test/test_bilinear_approximations.jl @@ -47,20 +47,13 @@ end @test haskey(result, ("dev1", 1)) @test result[("dev1", 1)] isa JuMP.AffExpr - # u = x + y variable container should exist + # p = x + y variable container should exist @test IOM.has_container_key( setup.container, IOM.BilinearApproxSumVariable, MockThermalGen, BILINEAR_META * "_plus", ) - # v = x - y variable container should exist - @test IOM.has_container_key( - setup.container, - IOM.BilinearApproxDiffVariable, - MockThermalGen, - BILINEAR_META * "_minus", - ) # z ≈ x·y variable container should exist @test IOM.has_container_key( setup.container, @@ -75,12 +68,6 @@ end MockThermalGen, BILINEAR_META * "_plus", ) - @test IOM.has_container_key( - setup.container, - IOM.BilinearApproxDiffLinkingConstraint, - MockThermalGen, - BILINEAR_META * "_minus", - ) # Inner quadratic approx containers should exist with _plus/_minus meta @test IOM.has_container_key( setup.container, @@ -88,14 +75,8 @@ end MockThermalGen, BILINEAR_META * "_plus", ) - @test IOM.has_container_key( - setup.container, - IOM.QuadraticApproxVariable, - MockThermalGen, - BILINEAR_META * "_minus", - ) - # u bounds should be [0+0, 4+4] = [0, 8] + # p bounds should be [0+0, 4+4] = [0, 8] u_container = IOM.get_variable( setup.container, IOM.BilinearApproxSumVariable(), @@ -104,16 +85,6 @@ end ) @test JuMP.lower_bound(u_container["dev1", 1]) == 0.0 @test JuMP.upper_bound(u_container["dev1", 1]) == 8.0 - - # v bounds should be [0-4, 4-0] = [-4, 4] - v_container = IOM.get_variable( - setup.container, - IOM.BilinearApproxDiffVariable(), - MockThermalGen, - BILINEAR_META * "_minus", - ) - @test JuMP.lower_bound(v_container["dev1", 1]) == -4.0 - @test JuMP.upper_bound(v_container["dev1", 1]) == 4.0 end @testset "Constraint structure with McCormick" begin @@ -355,12 +326,6 @@ end MockThermalGen, BILINEAR_META * "_plus", ) - @test IOM.has_container_key( - setup.container, - IOM.ManualSOS2BinaryVariable, - MockThermalGen, - BILINEAR_META * "_minus", - ) # NO solver SOS2 constraints sos2_count = JuMP.num_constraints( @@ -551,18 +516,6 @@ end MockThermalGen, BILINEAR_META * "_plus", ) - @test IOM.has_container_key( - setup.container, - IOM.SawtoothAuxVariable, - MockThermalGen, - BILINEAR_META * "_minus", - ) - @test IOM.has_container_key( - setup.container, - IOM.SawtoothBinaryVariable, - MockThermalGen, - BILINEAR_META * "_minus", - ) end @testset "Fixed-variable correctness" begin From 0e57ff47340202319cccc060b9c760b264bb1f26 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Thu, 5 Mar 2026 13:19:25 -0500 Subject: [PATCH 7/7] formatting --- src/quadratic_approximations/bilinear.jl | 42 +++++++++++------------- test/test_bilinear_approximations.jl | 2 +- 2 files changed, 21 insertions(+), 23 deletions(-) diff --git a/src/quadratic_approximations/bilinear.jl b/src/quadratic_approximations/bilinear.jl index 0d2866d..ba9efcc 100644 --- a/src/quadratic_approximations/bilinear.jl +++ b/src/quadratic_approximations/bilinear.jl @@ -104,28 +104,26 @@ function _add_bilinear_approx_impl!( p_max, meta_plus, ) - zx_dict = - quad_approx_fn( - container, - C, - names, - time_steps, - x_var_container, - x_min, - x_max, - meta_x, - ) - zy_dict = - quad_approx_fn( - container, - C, - names, - time_steps, - y_var_container, - y_min, - y_max, - meta_y, - ) + zx_dict = quad_approx_fn( + container, + C, + names, + time_steps, + x_var_container, + x_min, + x_max, + meta_x, + ) + zy_dict = quad_approx_fn( + container, + C, + names, + time_steps, + y_var_container, + y_min, + y_max, + meta_y, + ) # Create z variable container for the bilinear product z_container = add_variable_container!( diff --git a/test/test_bilinear_approximations.jl b/test/test_bilinear_approximations.jl index 6940ec6..e02e424 100644 --- a/test/test_bilinear_approximations.jl +++ b/test/test_bilinear_approximations.jl @@ -68,7 +68,7 @@ end MockThermalGen, BILINEAR_META * "_plus", ) - # Inner quadratic approx containers should exist with _plus/_minus meta + # Inner quadratic approx containers should exist with _plus meta @test IOM.has_container_key( setup.container, IOM.QuadraticApproxVariable,