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 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/src/quadratic_approximations/bilinear.jl b/src/quadratic_approximations/bilinear.jl new file mode 100644 index 0000000..ba9efcc --- /dev/null +++ b/src/quadratic_approximations/bilinear.jl @@ -0,0 +1,292 @@ +# 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)² + +struct BilinearApproxSumVariable <: VariableType end # p = x + y +struct BilinearApproxSumLinkingConstraint <: ConstraintType end # p == 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 = (1/2)((x+y)² − x² - y²). + +Creates auxiliary variables p = x+y, calls `quad_approx_fn` to +approximate p², then combines via multiplicative 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 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_x = meta * "_x" + meta_y = meta * "_y" + + # Create p variable container + p_container = add_variable_container!( + container, + BilinearApproxSumVariable(), + C, + names, + time_steps; + meta = meta_plus, + ) + + # Create linking constraint containers + p_link_container = add_constraints_container!( + container, + BilinearApproxSumLinkingConstraint(), + C, + names, + time_steps; + meta = meta_plus, + ) + + # 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] + + p_container[name, t] = JuMP.@variable( + jump_model, + base_name = "BilinearSum_$(C)_{$(name), $(t)}", + lower_bound = p_min, + upper_bound = p_max, + ) + + p_link_container[name, t] = + JuMP.@constraint(jump_model, p_container[name, t] == x + y) + end + + # 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_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!( + 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 = (1/2) * (p² − x² - y²) + z_expr = JuMP.AffExpr(0.0) + 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) + + 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 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 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..e02e424 --- /dev/null +++ b/test/test_bilinear_approximations.jl @@ -0,0 +1,667 @@ +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 + + # p = x + y variable container should exist + @test IOM.has_container_key( + setup.container, + IOM.BilinearApproxSumVariable, + MockThermalGen, + BILINEAR_META * "_plus", + ) + # 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", + ) + # Inner quadratic approx containers should exist with _plus meta + @test IOM.has_container_key( + setup.container, + IOM.QuadraticApproxVariable, + MockThermalGen, + BILINEAR_META * "_plus", + ) + + # p 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 + 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", + ) + + # 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", + ) + 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