From b914f6e68b3e9daace3eda5973362c290b4d5c51 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Wed, 5 Aug 2026 12:06:07 -0400 Subject: [PATCH] Initial --- GDPOptimizer.jl/Project.toml | 20 + GDPOptimizer.jl/src/GDPOptimizer.jl | 13 + GDPOptimizer.jl/src/cuts.jl | 202 ++++++++ GDPOptimizer.jl/src/loa.jl | 257 +++++++++ GDPOptimizer.jl/src/master.jl | 181 +++++++ GDPOptimizer.jl/src/nlp.jl | 168 ++++++ GDPOptimizer.jl/src/optimizer.jl | 370 +++++++++++++ GDPOptimizer.jl/src/problem.jl | 159 ++++++ GDPOptimizer.jl/src/sets.jl | 96 ++++ GDPOptimizer.jl/test/loa.jl | 489 ++++++++++++++++++ GDPOptimizer.jl/test/moi.jl | 35 ++ GDPOptimizer.jl/test/optimizer.jl | 92 ++++ GDPOptimizer.jl/test/runtests.jl | 9 + GDPOptimizer.jl/test/sets.jl | 48 ++ Project.toml | 14 +- ext/GDPOptimizerDisjunctiveProgramming.jl | 76 +++ src/extension_api.jl | 23 + test/aqua.jl | 5 +- .../GDPOptimizerDisjunctiveProgramming.jl | 170 ++++++ test/runtests.jl | 1 + 20 files changed, 2423 insertions(+), 5 deletions(-) create mode 100644 GDPOptimizer.jl/Project.toml create mode 100644 GDPOptimizer.jl/src/GDPOptimizer.jl create mode 100644 GDPOptimizer.jl/src/cuts.jl create mode 100644 GDPOptimizer.jl/src/loa.jl create mode 100644 GDPOptimizer.jl/src/master.jl create mode 100644 GDPOptimizer.jl/src/nlp.jl create mode 100644 GDPOptimizer.jl/src/optimizer.jl create mode 100644 GDPOptimizer.jl/src/problem.jl create mode 100644 GDPOptimizer.jl/src/sets.jl create mode 100644 GDPOptimizer.jl/test/loa.jl create mode 100644 GDPOptimizer.jl/test/moi.jl create mode 100644 GDPOptimizer.jl/test/optimizer.jl create mode 100644 GDPOptimizer.jl/test/runtests.jl create mode 100644 GDPOptimizer.jl/test/sets.jl create mode 100644 ext/GDPOptimizerDisjunctiveProgramming.jl create mode 100644 test/extensions/GDPOptimizerDisjunctiveProgramming.jl diff --git a/GDPOptimizer.jl/Project.toml b/GDPOptimizer.jl/Project.toml new file mode 100644 index 00000000..d3440858 --- /dev/null +++ b/GDPOptimizer.jl/Project.toml @@ -0,0 +1,20 @@ +name = "GDPOptimizer" +uuid = "7da86d0c-51ea-4770-ae9e-8cf6d5933256" +authors = ["Daniel Nguyen"] +version = "0.1.0" + +[deps] +MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" + +[compat] +MathOptInterface = "1" +julia = "1.10" + +[extras] +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["HiGHS", "Ipopt", "JuMP", "Test"] diff --git a/GDPOptimizer.jl/src/GDPOptimizer.jl b/GDPOptimizer.jl/src/GDPOptimizer.jl new file mode 100644 index 00000000..12fb3253 --- /dev/null +++ b/GDPOptimizer.jl/src/GDPOptimizer.jl @@ -0,0 +1,13 @@ +module GDPOptimizer + +import MathOptInterface as MOI + +include("sets.jl") +include("optimizer.jl") +include("problem.jl") +include("master.jl") +include("nlp.jl") +include("cuts.jl") +include("loa.jl") + +end diff --git a/GDPOptimizer.jl/src/cuts.jl b/GDPOptimizer.jl/src/cuts.jl new file mode 100644 index 00000000..079cc2e3 --- /dev/null +++ b/GDPOptimizer.jl/src/cuts.jl @@ -0,0 +1,202 @@ +################################################################################ +# LINEARIZATION +################################################################################ +# First-order Taylor expansions of the nonlinear rows, sharing one +# `MOI.Nonlinear` reverse-mode evaluator per row across iterations +# (only the evaluation point changes). +struct _Linearizer + evaluators::Dict{UInt64, + Tuple{MOI.Nonlinear.Evaluator, Vector{MOI.VariableIndex}}} +end + +_Linearizer() = _Linearizer(Dict{UInt64, + Tuple{MOI.Nonlinear.Evaluator, Vector{MOI.VariableIndex}}}()) + +function _append_variables( + variables::Vector{MOI.VariableIndex}, + func::MOI.VariableIndex + ) + return push!(variables, func) +end +function _append_variables( + variables::Vector{MOI.VariableIndex}, + func::MOI.ScalarAffineFunction{Float64} + ) + return append!(variables, term.variable for term in func.terms) +end +function _append_variables( + variables::Vector{MOI.VariableIndex}, + func::MOI.ScalarQuadraticFunction{Float64} + ) + append!(variables, term.variable for term in func.affine_terms) + for term in func.quadratic_terms + push!(variables, term.variable_1, term.variable_2) + end + return variables +end +function _append_variables( + variables::Vector{MOI.VariableIndex}, + func::MOI.ScalarNonlinearFunction + ) + for arg in func.args + arg isa Real || _append_variables(variables, arg) + end + return variables +end + +function _evaluator(linearizer::_Linearizer, func) + return get!(linearizer.evaluators, objectid(func)) do + variables = _append_variables(MOI.VariableIndex[], func) + unique!(variables) + nonlinear = MOI.Nonlinear.Model() + MOI.Nonlinear.set_objective(nonlinear, func) + evaluator = MOI.Nonlinear.Evaluator(nonlinear, + MOI.Nonlinear.SparseReverseMode(), variables) + MOI.initialize(evaluator, [:Grad]) + return (evaluator, variables) + end +end + +# Exact for affine functions, first-order Taylor at `point` otherwise. +# The result stays in cache space; callers remap when adding cuts. +_linearize(::_Linearizer, func::MOI.ScalarAffineFunction{Float64}, point) = func +_linearize(::_Linearizer, func::MOI.VariableIndex, point) = _to_affine(func) +function _linearize( + linearizer::_Linearizer, + func::Union{MOI.ScalarQuadraticFunction{Float64}, + MOI.ScalarNonlinearFunction}, + point::AbstractDict + ) + evaluator, variables = _evaluator(linearizer, func) + x = [point[vi] for vi in variables] + value = MOI.eval_objective(evaluator, x) + gradient = zeros(length(x)) + MOI.eval_objective_gradient(evaluator, gradient, x) + constant = value - sum(gradient[i] * x[i] for i in eachindex(x); + init = 0.0) + terms = [MOI.ScalarAffineTerm(gradient[i], variables[i]) + for i in eachindex(variables) if gradient[i] != 0.0] + return MOI.ScalarAffineFunction(terms, constant) +end + +################################################################################ +# OA CUT EMISSION +################################################################################ +_penalty_sign(sense::MOI.OptimizationSense) = sense == MOI.MAX_SENSE ? -1 : 1 + +# The `<= 0` directions of an OA cut for `set`: `lin - rhs` for +# LessThan, `rhs - lin` for GreaterThan, both for EqualTo / Interval. +_oa_cut_terms(set::MOI.LessThan{Float64}, lin) = + (MOI.Utilities.operate(-, Float64, lin, set.upper),) +_oa_cut_terms(set::MOI.GreaterThan{Float64}, lin) = + (MOI.Utilities.operate(-, Float64, set.lower, lin),) +_oa_cut_terms(set::MOI.EqualTo{Float64}, lin) = + (MOI.Utilities.operate(-, Float64, lin, set.value), + MOI.Utilities.operate(-, Float64, set.value, lin)) +_oa_cut_terms(set::MOI.Interval{Float64}, lin) = + (MOI.Utilities.operate(-, Float64, lin, set.upper), + MOI.Utilities.operate(-, Float64, set.lower, lin)) + +# Emit all OA cuts for one NLP result: the objective cut, a slacked row +# per nonlinear global, and a gated cut per active nonlinear disjunct +# row. Every slack keeps a nonconvex linearization from making the +# master infeasible. +function _add_oa_cuts( + model::Optimizer, + problem::_Problem, + master::_Master, + linearizer::_Linearizer, + result::NamedTuple + ) + result.point === nothing && return + sign = _penalty_sign(master.sense) + _add_objective_cut(model, master, linearizer, result.point, sign) + for (func, set) in problem.nonlinear_rows + _is_linear(func) && continue + lin = _master_linearization(master, linearizer, func, result.point) + _add_global_oa_row(model, master, lin, set, sign) + end + for disjunction in problem.disjunctions, disjunct in disjunction.disjuncts + _disjunct_active(result.combination, disjunct) || continue + for (func, set) in zip(disjunct.functions, disjunct.sets) + _is_linear(func) && continue + lin = _master_linearization(master, linearizer, func, result.point) + _add_disjunct_oa_cut(model, master, disjunct, lin, set, sign) + end + end + return +end + +function _master_linearization( + master::_Master, + linearizer::_Linearizer, + func, + point + ) + return _map_to(master.variable_map, _linearize(linearizer, func, point)) +end + +# Slacked objective cut. MIN: `lin <= alpha_oa + slack`; MAX symmetric. +function _add_objective_cut( + model::Optimizer, + master::_Master, + linearizer::_Linearizer, + point, + sign::Int + ) + lin = _master_linearization(master, linearizer, master.objective, point) + slack = _add_penalized_slack(master, model.options, sign) + alpha = _to_affine(master.alpha_oa) + if master.sense == MOI.MAX_SENSE + body = MOI.Utilities.operate(-, Float64, + MOI.Utilities.operate(-, Float64, alpha, lin), _to_affine(slack)) + else + body = MOI.Utilities.operate(-, Float64, + MOI.Utilities.operate(-, Float64, lin, alpha), _to_affine(slack)) + end + MOI.Utilities.normalize_and_add_constraint(master.model, body, + MOI.LessThan(0.0)) + return +end + +# Slacked global OA row(s): each direction `term - slack <= 0`. +function _add_global_oa_row( + model::Optimizer, + master::_Master, + lin::MOI.ScalarAffineFunction{Float64}, + set::MOI.AbstractScalarSet, + sign::Int + ) + slack = _add_penalized_slack(master, model.options, sign) + for term in _oa_cut_terms(set, lin) + body = MOI.Utilities.operate(-, Float64, term, _to_affine(slack)) + MOI.Utilities.normalize_and_add_constraint(master.model, body, + MOI.LessThan(0.0)) + end + return +end + +# Gated disjunct cut: each direction `term - slack <= M * (1 - z)` for +# an activation `z` (its complement gates on `M * z`). +function _add_disjunct_oa_cut( + model::Optimizer, + master::_Master, + disjunct::_Disjunct, + lin::MOI.ScalarAffineFunction{Float64}, + set::MOI.AbstractScalarSet, + sign::Int + ) + M = Float64(_option(model, "M_value")) + activation = _map_to(master.variable_map, disjunct.activation) + # M * (1 - activation) moved left: term - slack + M * activation - M + gate = MOI.Utilities.operate(-, Float64, + MOI.Utilities.operate(*, Float64, M, activation), M) + slack = _add_penalized_slack(master, model.options, sign) + for term in _oa_cut_terms(set, lin) + body = MOI.Utilities.operate(+, Float64, + MOI.Utilities.operate(-, Float64, term, _to_affine(slack)), gate) + MOI.Utilities.normalize_and_add_constraint(master.model, body, + MOI.LessThan(0.0)) + end + return +end diff --git a/GDPOptimizer.jl/src/loa.jl b/GDPOptimizer.jl/src/loa.jl new file mode 100644 index 00000000..d68fdd6a --- /dev/null +++ b/GDPOptimizer.jl/src/loa.jl @@ -0,0 +1,257 @@ +################################################################################ +# LOGIC-BASED OUTER APPROXIMATION +################################################################################ +_worst_objective(sense::MOI.OptimizationSense) = + sense == MOI.MAX_SENSE ? -Inf : Inf +_is_better(sense::MOI.OptimizationSense, new, best) = + sense == MOI.MAX_SENSE ? new > best : new < best +_gap(sense::MOI.OptimizationSense, best, bound) = + sense == MOI.MAX_SENSE ? bound - best : best - bound + +# one record per NLP solve, for convergence traces +function _log_progress( + model::Optimizer, + t_start::Float64, + best_objective, + master_bound + ) + model.silent && return + bound = master_bound === nothing ? NaN : master_bound + @info "LOA progress: elapsed=$(time() - t_start) " * + "incumbent=$best_objective bound=$bound" + return +end + +# one target per (binary, value); linear disjuncts need no cover +# since the master already carries them exactly +function _cover_disjuncts(problem::_Problem) + seen = Set{Tuple{MOI.VariableIndex, Bool}}() + cover = _Disjunct[] + for disjunction in problem.disjunctions, disjunct in disjunction.disjuncts + _is_nonlinear_disjunct(disjunct) || continue + key = (disjunct.binary, disjunct.active_value) + key in seen && continue + push!(seen, key) + push!(cover, disjunct) + end + return cover +end + +# GDPopt's covering weights: an uncovered disjunct outweighs all +# covered ones +function _cover_objective( + master::_Master, + cover::Vector{_Disjunct}, + needs_cover, + num_covered::Int + ) + objective = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0) + for i in eachindex(cover) + weight = Float64(needs_cover[i] ? num_covered + 1 : 1) + activation = _map_to(master.variable_map, cover[i].activation) + objective = MOI.Utilities.operate(+, Float64, objective, + MOI.Utilities.operate(*, Float64, weight, activation)) + end + return objective +end + +# variable starts warm start the first NLP; later iterations use +# the last feasible primal +function _user_start_values(model::Optimizer, problem::_Problem) + point = Dict{MOI.VariableIndex, Float64}() + for vi in problem.variables + start = MOI.get(model.cache, MOI.VariablePrimalStart(), vi) + start === nothing || (point[vi] = start) + end + return isempty(point) ? nothing : (point = point,) +end + +function _set_master_objective(master::_Master, sense, objective) + MOI.set(master.model, MOI.ObjectiveSense(), sense) + MOI.set(master.model, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + objective) + return +end + +################################################################################ +# MAIN LOOP +################################################################################ +function MOI.optimize!(model::Optimizer) + t_start = time() + _reset_results(model) + problem = _build_problem(model) + master = _build_master(model, problem) + subproblem = _build_subproblem(model, problem) + linearizer = _Linearizer() + sense = problem.sense + overall_deadline = t_start + Float64(_option(model, "time_limit")) + loop_deadline = min(overall_deadline, + t_start + Float64(_option(model, "iteration_time_limit"))) + + best_objective = _worst_objective(sense) + best_result = nothing + previous_result = _user_start_values(model, problem) + master_bound = nothing + master_status = nothing + converged = false + + # Shared iteration tail: no-good cut, OA cuts, incumbent update. + process_result = result -> begin + _avoid_combination(master, result.combination) + _add_oa_cuts(model, problem, master, linearizer, result) + if result.feasible && + _is_better(sense, result.objective, best_objective) + best_objective = result.objective + best_result = result + end + result.feasible && (previous_result = result) + _log_progress(model, t_start, best_objective, master_bound) + return + end + warm_start = () -> + previous_result === nothing ? nothing : previous_result.point + + # set covering: reuse the master with a coverage objective so + # every nonlinear disjunct gets visited once + cover = _cover_disjuncts(problem) + needs_cover = trues(length(cover)) + num_covered = 0 + for iteration in 1:_option(model, "set_cover_max_iter") + (iteration == 1 || any(needs_cover)) || break + time() < loop_deadline || break + _set_master_objective(master, MOI.MAX_SENSE, + _cover_objective(master, cover, needs_cover, num_covered)) + _cap_remaining_time(master.model, loop_deadline) + MOI.optimize!(master.model) + solved = _solved_and_feasible(master.model) + # capture the status before the objective restore invalidates it + status = MOI.get(master.model, MOI.TerminationStatus()) + combination = solved ? _extract_combination(problem, master) : nothing + _set_master_objective(master, master.sense, master.oa_objective) + if !solved + master_status = status + break + end + result = _solve_nlp(model, problem, subproblem, combination, + warm_start(); deadline = loop_deadline) + process_result(result) + # covered only once active in a feasible NLP; infeasible + # combinations just leave their no-good cut + if result.feasible + for i in eachindex(cover) + needs_cover[i] || continue + _disjunct_active(result.combination, cover[i]) && + (needs_cover[i] = false) + end + num_covered = count(!, needs_cover) + end + end + + # main loop: alpha_oa gives the bound, the NLP the incumbent + if master_status === nothing + for _ in 1:_option(model, "max_iter") + time() < loop_deadline || break + _cap_remaining_time(master.model, loop_deadline) + MOI.optimize!(master.model) + if !_solved_and_feasible(master.model) + master_status = MOI.get(master.model, MOI.TerminationStatus()) + break + end + master_bound = MOI.get(master.model, MOI.VariablePrimal(), + master.alpha_oa) + if best_result !== nothing + gap = _gap(sense, best_objective, master_bound) + total_slack = abs(MOI.get(master.model, + MOI.ObjectiveValue()) - master_bound) / + Float64(_option(model, "oa_penalty")) + tol = Float64(_option(model, "convergence_tol")) * + max(abs(best_objective), 1.0) + if gap <= tol && + total_slack <= Float64(_option(model, "slack_tol")) + converged = true + break + end + end + combination = _extract_combination(problem, master) + result = _solve_nlp(model, problem, subproblem, combination, + warm_start(); deadline = loop_deadline) + process_result(result) + end + end + + _store_results(model, sense, best_objective, best_result, master_bound, + master_status, converged, loop_deadline) + model.solve_time = time() - t_start + return +end + +################################################################################ +# RESULT SYNTHESIS +################################################################################ +# the OA bound is valid only for convex problems, so report +# LOCALLY_SOLVED, never OPTIMAL +function _store_results( + model::Optimizer, + sense::MOI.OptimizationSense, + best_objective::Float64, + best_result, + master_bound, + master_status, + converged::Bool, + loop_deadline::Float64 + ) + timed_out = time() >= loop_deadline + if best_result === nothing + model.primal_status = MOI.NO_SOLUTION + model.objective_value = NaN + model.raw_status = "No feasible incumbent found." + if master_status == MOI.INFEASIBLE + model.termination_status = MOI.INFEASIBLE + model.raw_status = "No feasible incumbent: the master " * + "problem is infeasible." + elseif timed_out + model.termination_status = MOI.TIME_LIMIT + elseif master_status === nothing + model.termination_status = MOI.ITERATION_LIMIT + else + model.termination_status = MOI.OTHER_LIMIT + model.raw_status = "No feasible incumbent: the master " * + "solve finished with status $master_status." + end + return + end + model.primal_status = MOI.FEASIBLE_POINT + model.incumbent = best_result.point + model.objective_value = best_objective + model.objective_bound = master_bound === nothing ? nothing : + Float64(master_bound) + if converged + model.termination_status = MOI.LOCALLY_SOLVED + elseif timed_out + model.termination_status = MOI.TIME_LIMIT + elseif master_status == MOI.INFEASIBLE + # all combinations visited; the incumbent is best over all + # of them, but the OA bound is gone + model.termination_status = MOI.LOCALLY_SOLVED + elseif master_status !== nothing + model.termination_status = MOI.OTHER_LIMIT + else + model.termination_status = MOI.ITERATION_LIMIT + end + label = converged ? "converged" : (master_status == MOI.INFEASIBLE ? + "combinations exhausted" : "limit hit") + if master_bound === nothing + model.raw_status = "LOA finished [$label]: incumbent " * + "$best_objective (master produced no bound)." + else + gap = _gap(sense, best_objective, master_bound) + relative = abs(best_objective) > 1e-10 ? + gap / abs(best_objective) : gap + model.relative_gap = relative + model.raw_status = "LOA finished [$label]: incumbent " * + "$best_objective, master bound $master_bound, gap $gap " * + "(relative $relative)." + end + return +end diff --git a/GDPOptimizer.jl/src/master.jl b/GDPOptimizer.jl/src/master.jl new file mode 100644 index 00000000..8d3dd2d7 --- /dev/null +++ b/GDPOptimizer.jl/src/master.jl @@ -0,0 +1,181 @@ +################################################################################ +# MASTER CONSTRUCTION +################################################################################ +# `alpha_oa` carries the objective; `oa_objective` also tracks the +# slack penalties so set covering can swap objectives out and back +mutable struct _Master + model::MOI.ModelLike + variable_map::Dict{MOI.VariableIndex, MOI.VariableIndex} + sense::MOI.OptimizationSense + objective::MOI.AbstractScalarFunction + alpha_oa::MOI.VariableIndex + oa_objective::MOI.ScalarAffineFunction{Float64} +end + +function _instantiate(factory, name::String) + factory === nothing && + error("GDPOptimizer requires the `$name` optimizer factory.") + solver = MOI.instantiate(factory; + with_cache_type = Float64, with_bridge_type = Float64) + MOI.set(solver, MOI.Silent(), true) + return solver +end + +_map_to(variable_map::AbstractDict, func) = + MOI.Utilities.map_indices(vi -> variable_map[vi], func) + +function _solved_and_feasible(solver::MOI.ModelLike) + status = MOI.get(solver, MOI.TerminationStatus()) + return status in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) && + MOI.get(solver, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT +end + +function _cap_remaining_time(solver::MOI.ModelLike, deadline::Float64) + isfinite(deadline) || return + MOI.set(solver, MOI.TimeLimitSec(), max(0.0, deadline - time())) + return +end + +function _build_master(model::Optimizer, problem::_Problem) + mip = _instantiate(model.mip_solver, "mip_solver") + variable_map = Dict{MOI.VariableIndex, MOI.VariableIndex}( + vi => MOI.add_variable(mip) for vi in problem.variables) + for ci in problem.variable_cis + vi = MOI.get(model.cache, MOI.ConstraintFunction(), ci) + MOI.add_constraint(mip, variable_map[vi], + MOI.get(model.cache, MOI.ConstraintSet(), ci)) + end + for ci in problem.linear_cis + func = MOI.get(model.cache, MOI.ConstraintFunction(), ci) + MOI.add_constraint(mip, _map_to(variable_map, func), + MOI.get(model.cache, MOI.ConstraintSet(), ci)) + end + # nonlinear-typed rows that demoted to affine stay in the master + for (func, set) in problem.nonlinear_rows + _is_linear(func) || continue + MOI.add_constraint(mip, _map_to(variable_map, _to_affine(func)), set) + end + for disjunction in problem.disjunctions + _add_exactly_one(mip, variable_map, disjunction) + for disjunct in disjunction.disjuncts + _add_gated_rows(mip, variable_map, disjunct, model.options) + end + end + alpha_oa = MOI.add_variable(mip) + sense = problem.sense + oa_objective = MOI.ScalarAffineFunction( + [MOI.ScalarAffineTerm(1.0, alpha_oa)], 0.0) + MOI.set(mip, MOI.ObjectiveSense(), sense) + MOI.set(mip, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + oa_objective) + return _Master(mip, variable_map, sense, problem.objective, alpha_oa, + oa_objective) +end + +# Indicators sum to the activation (1 top-level, parent indicator +# nested); a complement pair normalizes to a trivial `1 == 1` row +function _add_exactly_one( + mip::MOI.ModelLike, + variable_map::AbstractDict, + disjunction::_Disjunction + ) + total = MOI.Utilities.operate(+, Float64, + (_map_to(variable_map, disjunct.activation) + for disjunct in disjunction.disjuncts)...) + body = MOI.Utilities.operate(-, Float64, total, + _map_to(variable_map, disjunction.activation)) + MOI.Utilities.normalize_and_add_constraint(mip, body, MOI.EqualTo(0.0)) + return +end + +# split Interval rows; the indicator bridge takes one-sided sets only +_indicator_sets(set::MOI.Interval{Float64}) = + (MOI.LessThan(set.upper), MOI.GreaterThan(set.lower)) +_indicator_sets(set::MOI.AbstractScalarSet) = (set,) + +# Gate each linear row with an indicator constraint or big-M per +# `master_gating`; nonlinear rows enter the master only as OA cuts +function _add_gated_rows( + mip::MOI.ModelLike, + variable_map::AbstractDict, + disjunct::_Disjunct, + options::Dict{String, Any} + ) + gating = options["master_gating"] + gating in ("indicator", "bigm") || + error("Unknown `master_gating` value `$gating`.") + activate = disjunct.active_value ? MOI.ACTIVATE_ON_ONE : + MOI.ACTIVATE_ON_ZERO + binary = variable_map[disjunct.binary] + activation = _map_to(variable_map, disjunct.activation) + M = Float64(options["M_value"]) + # M * (1 - activation) moved left, as in the disjunct OA cuts + gate = MOI.Utilities.operate(-, Float64, + MOI.Utilities.operate(*, Float64, M, activation), M) + for (func, set) in zip(disjunct.functions, disjunct.sets) + _is_linear(func) || continue + row = _map_to(variable_map, _to_affine(func)) + if gating == "bigm" + for term in _oa_cut_terms(set, row) + body = MOI.Utilities.operate(+, Float64, term, gate) + MOI.Utilities.normalize_and_add_constraint(mip, body, + MOI.LessThan(0.0)) + end + else + for inner in _indicator_sets(set) + gated = MOI.Utilities.operate(vcat, Float64, binary, row) + MOI.add_constraint(mip, gated, + MOI.Indicator{activate}(inner)) + end + end + end + return +end + +# bounded penalized slack so an invalid cut cannot blow up the master +function _add_penalized_slack( + master::_Master, + options::Dict{String, Any}, + penalty_sign::Int + ) + slack = MOI.add_variable(master.model) + MOI.add_constraint(master.model, slack, MOI.GreaterThan(0.0)) + MOI.add_constraint(master.model, slack, + MOI.LessThan(Float64(options["max_slack"]))) + penalty = penalty_sign * Float64(options["oa_penalty"]) + master.oa_objective = MOI.Utilities.operate(+, Float64, + master.oa_objective, MOI.ScalarAffineFunction( + [MOI.ScalarAffineTerm(penalty, slack)], 0.0)) + MOI.set(master.model, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + master.oa_objective) + return slack +end + +# round to Bool; MILP values are only within integer tolerance +function _extract_combination(problem::_Problem, master::_Master) + return Dict{MOI.VariableIndex, Bool}( + binary => round(Bool, MOI.get(master.model, MOI.VariablePrimal(), + master.variable_map[binary])) + for binary in problem.binaries) +end + +# No-good cut: active `1 - z` plus inactive `z` terms must reach 1 +function _avoid_combination(master::_Master, combination::AbstractDict) + terms = MOI.ScalarAffineTerm{Float64}[] + constant = 0.0 + for (binary, value) in combination + mapped = master.variable_map[binary] + if value + push!(terms, MOI.ScalarAffineTerm(-1.0, mapped)) + constant += 1.0 + else + push!(terms, MOI.ScalarAffineTerm(1.0, mapped)) + end + end + cut = MOI.ScalarAffineFunction(terms, constant) + MOI.Utilities.normalize_and_add_constraint(master.model, cut, + MOI.GreaterThan(1.0)) + return +end diff --git a/GDPOptimizer.jl/src/nlp.jl b/GDPOptimizer.jl/src/nlp.jl new file mode 100644 index 00000000..eaff8575 --- /dev/null +++ b/GDPOptimizer.jl/src/nlp.jl @@ -0,0 +1,168 @@ +################################################################################ +# NLP SUBPROBLEM +################################################################################ +# Built once; each iteration overwrites the binary fixes in place +# and swaps the active disjuncts' rows. No big-M anywhere. +struct _Subproblem + model::MOI.ModelLike + variable_map::Dict{MOI.VariableIndex, MOI.VariableIndex} + fixes::Dict{MOI.VariableIndex, + MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{Float64}}} + rows::Vector{MOI.ConstraintIndex} +end + +function _build_subproblem(model::Optimizer, problem::_Problem) + nlp = _instantiate(model.nlp_solver, "nlp_solver") + variable_map = Dict{MOI.VariableIndex, MOI.VariableIndex}( + vi => MOI.add_variable(nlp) for vi in problem.variables) + indicators = Set(problem.binaries) + for ci in problem.variable_cis + vi = MOI.get(model.cache, MOI.ConstraintFunction(), ci) + vi in indicators && continue + MOI.add_constraint(nlp, variable_map[vi], + MOI.get(model.cache, MOI.ConstraintSet(), ci)) + end + fixes = Dict(binary => MOI.add_constraint(nlp, variable_map[binary], + MOI.EqualTo(0.0)) for binary in problem.binaries) + for ci in problem.linear_cis + func = MOI.get(model.cache, MOI.ConstraintFunction(), ci) + MOI.add_constraint(nlp, _map_to(variable_map, func), + MOI.get(model.cache, MOI.ConstraintSet(), ci)) + end + for (func, set) in problem.nonlinear_rows + MOI.add_constraint(nlp, _map_to(variable_map, func), set) + end + MOI.set(nlp, MOI.ObjectiveSense(), problem.sense) + objective = _map_to(variable_map, problem.objective) + MOI.set(nlp, MOI.ObjectiveFunction{typeof(objective)}(), objective) + return _Subproblem(nlp, variable_map, fixes, MOI.ConstraintIndex[]) +end + +function _set_warm_start(nlp::MOI.ModelLike, variable_map::AbstractDict, point) + point === nothing && return + for (vi, value) in point + MOI.set(nlp, MOI.VariablePrimalStart(), variable_map[vi], value) + end + return +end + +function _extract_point( + nlp::MOI.ModelLike, + problem::_Problem, + variable_map::AbstractDict + ) + return Dict{MOI.VariableIndex, Float64}( + vi => MOI.get(nlp, MOI.VariablePrimal(), variable_map[vi]) + for vi in problem.variables) +end + +# Solve the NLP at a fixed combination: overwrite the binary fixes, +# swap the active disjuncts' rows, and optimize. If infeasible, fall +# through to NLPF (a slacked version that always solves) so the master +# still gets a linearization site, not just a no-good cut. +function _solve_nlp( + model::Optimizer, + problem::_Problem, + sub::_Subproblem, + combination::AbstractDict, + warm_start; + deadline::Float64 = Inf + ) + for (binary, value) in combination + MOI.set(sub.model, MOI.ConstraintSet(), sub.fixes[binary], + MOI.EqualTo(value ? 1.0 : 0.0)) + end + for ci in sub.rows + MOI.delete(sub.model, ci) + end + empty!(sub.rows) + for disjunction in problem.disjunctions, disjunct in disjunction.disjuncts + _disjunct_active(combination, disjunct) || continue + for (func, set) in zip(disjunct.functions, disjunct.sets) + push!(sub.rows, MOI.add_constraint(sub.model, + _map_to(sub.variable_map, func), set)) + end + end + _set_warm_start(sub.model, sub.variable_map, warm_start) + _cap_remaining_time(sub.model, deadline) + MOI.optimize!(sub.model) + if _solved_and_feasible(sub.model) + return (combination = combination, + point = _extract_point(sub.model, problem, sub.variable_map), + objective = MOI.get(sub.model, MOI.ObjectiveValue()), + feasible = true) + end + if Bool(_option(model, "use_nlpf")) + result = _solve_nlpf(model, problem, combination, warm_start; + deadline = deadline) + result === nothing || return result + end + return (combination = combination, + point = nothing, objective = Inf, feasible = false) +end + +################################################################################ +# NLPF (FEASIBILITY SUBPROBLEM) +################################################################################ +_nlpf_slacked(func, u, ::MOI.LessThan{Float64}) = + MOI.Utilities.operate(-, Float64, func, u) +_nlpf_slacked(func, u, ::MOI.GreaterThan{Float64}) = + MOI.Utilities.operate(+, Float64, func, u) +_nlpf_slacked(func, u, ::MOI.AbstractScalarSet) = nothing + +# The slacked feasibility NLP: one nonnegative `u` relaxes every scalar +# inequality row (bounds and equalities stay exact) and is minimized. +# Its solution is a linearization site for an infeasible combination. +function _solve_nlpf( + model::Optimizer, + problem::_Problem, + combination::AbstractDict, + warm_start; + deadline::Float64 = Inf + ) + nlp = _instantiate(model.nlp_solver, "nlp_solver") + variable_map = Dict{MOI.VariableIndex, MOI.VariableIndex}( + vi => MOI.add_variable(nlp) for vi in problem.variables) + u = MOI.add_variable(nlp) + MOI.add_constraint(nlp, u, MOI.GreaterThan(0.0)) + for ci in problem.variable_cis + vi = MOI.get(model.cache, MOI.ConstraintFunction(), ci) + haskey(combination, vi) && continue + MOI.add_constraint(nlp, variable_map[vi], + MOI.get(model.cache, MOI.ConstraintSet(), ci)) + end + for (binary, value) in combination + MOI.add_constraint(nlp, variable_map[binary], + MOI.EqualTo(value ? 1.0 : 0.0)) + end + rows = Tuple{MOI.AbstractScalarFunction, MOI.AbstractScalarSet}[] + for ci in problem.linear_cis + push!(rows, (MOI.get(model.cache, MOI.ConstraintFunction(), ci), + MOI.get(model.cache, MOI.ConstraintSet(), ci))) + end + append!(rows, problem.nonlinear_rows) + for disjunction in problem.disjunctions, disjunct in disjunction.disjuncts + _disjunct_active(combination, disjunct) || continue + append!(rows, zip(disjunct.functions, disjunct.sets)) + end + for (func, set) in rows + mapped = _map_to(variable_map, func) + slacked = _nlpf_slacked(mapped, u, set) + if slacked === nothing + MOI.add_constraint(nlp, mapped, set) + else + MOI.add_constraint(nlp, slacked, set) + end + end + MOI.set(nlp, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.set(nlp, MOI.ObjectiveFunction{MOI.VariableIndex}(), u) + _set_warm_start(nlp, variable_map, warm_start) + _cap_remaining_time(nlp, deadline) + MOI.optimize!(nlp) + # Use the primal only at a genuine feasible point; a solver can + # report values at a nonfeasible/NaN primal that poisons the cut. + _solved_and_feasible(nlp) || return nothing + return (combination = combination, + point = _extract_point(nlp, problem, variable_map), + objective = Inf, feasible = false) +end diff --git a/GDPOptimizer.jl/src/optimizer.jl b/GDPOptimizer.jl/src/optimizer.jl new file mode 100644 index 00000000..1bfc110f --- /dev/null +++ b/GDPOptimizer.jl/src/optimizer.jl @@ -0,0 +1,370 @@ +################################################################################ +# OPTIMIZER +################################################################################ +const _Cache = MOI.Utilities.UniversalFallback{MOI.Utilities.Model{Float64}} + +# Raw options, mirroring DisjunctiveProgramming.jl's LOA defaults. +const _DEFAULT_OPTIONS = Dict{String, Any}( + "max_iter" => 10, + "set_cover_max_iter" => 8, + "M_value" => 1e9, + "master_gating" => "indicator", + "max_slack" => 1e3, + "oa_penalty" => 1e3, + "use_nlpf" => true, + "convergence_tol" => 1e-6, + "slack_tol" => 1e-4, + "iteration_time_limit" => Inf, + "time_limit" => 3600.0, +) + +""" + Optimizer(; nlp_solver, mip_solver = nlp_solver, kwargs...) + +Logic-based outer approximation solver for models containing +[`DisjunctionSet`](@ref) constraints. `nlp_solver` and `mip_solver` +are optimizer factories as accepted by `MOI.instantiate`. The +remaining keyword arguments set raw options (also reachable through +`MOI.RawOptimizerAttribute`): + +- `max_iter = 10`: master/NLP iterations after the set-covering seed. +- `set_cover_max_iter = 8`: set-covering initialization iterations. +- `M_value = 1e9`: big-M gating the disjunct OA cuts in the master. +- `master_gating = "indicator"`: how linear disjunct rows enter the + master, `"indicator"` (constraint gated by the binary) or `"bigm"` + (rows relaxed by `M_value * (1 - z)`, tighter for solvers that + cannot strengthen indicators, e.g. with presolve disabled). +- `max_slack = 1e3`: upper bound of each OA cut slack. +- `oa_penalty = 1e3`: objective penalty per unit of cut slack. +- `use_nlpf = true`: solve a slacked feasibility NLP when the primary + NLP is infeasible, so its point still seeds OA cuts. +- `convergence_tol = 1e-6`: relative incumbent/bound gap tolerance. +- `slack_tol = 1e-4`: total cut slack tolerance for convergence. +- `iteration_time_limit = Inf`: seconds allotted to the LOA loop. +- `time_limit = 3600.0`: overall seconds budget. +""" +mutable struct Optimizer <: MOI.AbstractOptimizer + nlp_solver::Any + mip_solver::Any + cache::_Cache + options::Dict{String, Any} + silent::Bool + # results (filled by MOI.optimize!) + termination_status::MOI.TerminationStatusCode + primal_status::MOI.ResultStatusCode + incumbent::Dict{MOI.VariableIndex, Float64} + objective_value::Float64 + objective_bound::Union{Nothing, Float64} + relative_gap::Float64 + raw_status::String + solve_time::Float64 +end + +function Optimizer(; + nlp_solver = nothing, + mip_solver = nlp_solver, + kwargs... + ) + options = copy(_DEFAULT_OPTIONS) + for (key, value) in kwargs + haskey(options, string(key)) || throw(ArgumentError( + "Unknown option `$key`.")) + options[string(key)] = value + end + return Optimizer(nlp_solver, mip_solver, + MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), + options, false, MOI.OPTIMIZE_NOT_CALLED, MOI.NO_SOLUTION, + Dict{MOI.VariableIndex, Float64}(), NaN, nothing, NaN, "", NaN) +end + +_option(model::Optimizer, name::String) = model.options[name] + +MOI.get(::Optimizer, ::MOI.SolverName) = "GDPOptimizer" +MOI.get(::Optimizer, ::MOI.SolverVersion) = "0.1.0" + +MOI.is_empty(model::Optimizer) = MOI.is_empty(model.cache) + +function MOI.empty!(model::Optimizer) + MOI.empty!(model.cache) + _reset_results(model) + return +end + +# The non-cache part of MOI.empty!, also run at the top of every solve +# so a re-optimize cannot leak the previous solve's bound, gap, or +# point. +function _reset_results(model::Optimizer) + model.termination_status = MOI.OPTIMIZE_NOT_CALLED + model.primal_status = MOI.NO_SOLUTION + empty!(model.incumbent) + model.objective_value = NaN + model.objective_bound = nothing + model.relative_gap = NaN + model.raw_status = "" + model.solve_time = NaN + return +end + +MOI.supports_incremental_interface(::Optimizer) = true + +function MOI.copy_to(model::Optimizer, src::MOI.ModelLike) + return MOI.Utilities.default_copy_to(model, src) +end + +################################################################################ +# MODEL-BUILDING FORWARDING +################################################################################ +# Everything model-building lands in the cache; the solve partitions it +# in MOI.optimize!, so incremental additions need no invalidation. +const _Index = Union{MOI.VariableIndex, MOI.ConstraintIndex} + +MOI.add_variable(model::Optimizer) = MOI.add_variable(model.cache) + +MOI.add_variables(model::Optimizer, n::Int) = MOI.add_variables(model.cache, n) + +function MOI.add_constrained_variable( + model::Optimizer, + set::MOI.AbstractScalarSet + ) + return MOI.add_constrained_variable(model.cache, set) +end + +function MOI.add_constrained_variables( + model::Optimizer, + set::MOI.AbstractVectorSet + ) + return MOI.add_constrained_variables(model.cache, set) +end + +function MOI.add_constraint( + model::Optimizer, + func::MOI.AbstractFunction, + set::MOI.AbstractSet + ) + return MOI.add_constraint(model.cache, func, set) +end + +# Honest constraint support (the cache would claim everything, which +# disables the bridges that rewrite e.g. vector cones into supported +# scalar rows): only what `_build_problem` actually partitions. +const _ScalarFunction = Union{MOI.ScalarAffineFunction{Float64}, + MOI.ScalarQuadraticFunction{Float64}, MOI.ScalarNonlinearFunction} +const _VectorFunction = Union{MOI.VectorOfVariables, + MOI.VectorAffineFunction{Float64}, + MOI.VectorQuadraticFunction{Float64}, MOI.VectorNonlinearFunction} + +function MOI.supports_constraint( + ::Optimizer, + ::Type{<:MOI.AbstractFunction}, + ::Type{<:MOI.AbstractSet} + ) + return false +end + +# Non-indicator discrete variables pass through to the subproblems +# with their integrality intact, so the nlp_solver must handle them. +function MOI.supports_constraint( + ::Optimizer, + ::Type{MOI.VariableIndex}, + ::Type{<:Union{_SupportedInnerSet, MOI.ZeroOne, MOI.Integer}} + ) + return true +end + +function MOI.supports_constraint( + ::Optimizer, + ::Type{<:_ScalarFunction}, + ::Type{<:_SupportedInnerSet} + ) + return true +end + +function MOI.supports_constraint( + ::Optimizer, + ::Type{<:_VectorFunction}, + ::Type{DisjunctionSet} + ) + return true +end + +# A `DisjunctionSet` never constrains variables on creation: without +# this, `copy_to` turns a `VectorOfVariables` disjunction (which may +# repeat a variable across rows) into `add_constrained_variables`. +function MOI.supports_add_constrained_variables( + ::Optimizer, + ::Type{DisjunctionSet} + ) + return false +end + +MOI.is_valid(model::Optimizer, index::_Index) = + MOI.is_valid(model.cache, index) + +MOI.delete(model::Optimizer, index::_Index) = MOI.delete(model.cache, index) + +function MOI.get(model::Optimizer, T::Type{<:_Index}, name::String) + return MOI.get(model.cache, T, name) +end + +# Non-result attributes forward to the cache; the result attributes +# below override these generic methods, and any solve-set attribute +# without an override (duals, ConstraintPrimal) is refused rather than +# forwarded to the cache, which never holds results. +# Forward to the wrapped inner model (honest), not the UniversalFallback, +# which claims support for every attribute and so never rejects an +# unsupported one on copy_to. +MOI.supports(model::Optimizer, attr::MOI.AbstractModelAttribute) = + MOI.supports(model.cache.model, attr) + +function MOI.set(model::Optimizer, attr::MOI.AbstractModelAttribute, value) + MOI.supports(model, attr) || throw(MOI.UnsupportedAttribute(attr)) + return MOI.set(model.cache, attr, value) +end + +function MOI.get(model::Optimizer, attr::MOI.AbstractModelAttribute) + MOI.is_set_by_optimize(attr) && throw(MOI.GetAttributeNotAllowed(attr)) + return MOI.get(model.cache, attr) +end + +function MOI.supports( + model::Optimizer, + attr::MOI.AbstractVariableAttribute, + ::Type{MOI.VariableIndex} + ) + return MOI.supports(model.cache.model, attr, MOI.VariableIndex) +end + +function MOI.set( + model::Optimizer, + attr::MOI.AbstractVariableAttribute, + vi::MOI.VariableIndex, + value + ) + MOI.supports(model, attr, MOI.VariableIndex) || + throw(MOI.UnsupportedAttribute(attr)) + return MOI.set(model.cache, attr, vi, value) +end + +function MOI.get( + model::Optimizer, + attr::MOI.AbstractVariableAttribute, + vi::MOI.VariableIndex + ) + MOI.is_set_by_optimize(attr) && throw(MOI.GetAttributeNotAllowed(attr)) + return MOI.get(model.cache, attr, vi) +end + +function MOI.supports( + model::Optimizer, + attr::MOI.AbstractConstraintAttribute, + C::Type{<:MOI.ConstraintIndex} + ) + return MOI.supports(model.cache.model, attr, C) +end + +function MOI.set( + model::Optimizer, + attr::MOI.AbstractConstraintAttribute, + ci::MOI.ConstraintIndex, + value + ) + MOI.supports(model, attr, typeof(ci)) || + throw(MOI.UnsupportedAttribute(attr)) + return MOI.set(model.cache, attr, ci, value) +end + +function MOI.get( + model::Optimizer, + attr::MOI.AbstractConstraintAttribute, + ci::MOI.ConstraintIndex + ) + MOI.is_set_by_optimize(attr) && throw(MOI.GetAttributeNotAllowed(attr)) + return MOI.get(model.cache, attr, ci) +end + +################################################################################ +# OPTIMIZER ATTRIBUTES +################################################################################ +MOI.supports(::Optimizer, ::MOI.Silent) = true + +function MOI.set(model::Optimizer, ::MOI.Silent, value::Bool) + model.silent = value + return +end + +MOI.get(model::Optimizer, ::MOI.Silent) = model.silent + +MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true + +function MOI.set( + model::Optimizer, + ::MOI.TimeLimitSec, + value::Union{Nothing, Real} + ) + model.options["time_limit"] = value === nothing ? Inf : Float64(value) + return +end + +function MOI.get(model::Optimizer, ::MOI.TimeLimitSec) + limit = model.options["time_limit"] + return isfinite(limit) ? limit : nothing +end + +function MOI.supports(model::Optimizer, attr::MOI.RawOptimizerAttribute) + return haskey(model.options, attr.name) +end + +function MOI.set(model::Optimizer, attr::MOI.RawOptimizerAttribute, value) + MOI.supports(model, attr) || throw(MOI.UnsupportedAttribute(attr)) + model.options[attr.name] = value + return +end + +function MOI.get(model::Optimizer, attr::MOI.RawOptimizerAttribute) + MOI.supports(model, attr) || throw(MOI.UnsupportedAttribute(attr)) + return model.options[attr.name] +end + +################################################################################ +# RESULT ATTRIBUTES +################################################################################ +MOI.get(model::Optimizer, ::MOI.TerminationStatus) = model.termination_status + +MOI.get(model::Optimizer, ::MOI.RawStatusString) = model.raw_status + +function MOI.get(model::Optimizer, ::MOI.ResultCount) + return model.primal_status == MOI.NO_SOLUTION ? 0 : 1 +end + +function MOI.get(model::Optimizer, attr::MOI.PrimalStatus) + return attr.result_index == 1 ? model.primal_status : MOI.NO_SOLUTION +end + +MOI.get(model::Optimizer, ::MOI.DualStatus) = MOI.NO_SOLUTION + +function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue) + MOI.check_result_index_bounds(model, attr) + return model.objective_value +end + +function MOI.get(model::Optimizer, ::MOI.ObjectiveBound) + bound = model.objective_bound + if bound === nothing + sense = MOI.get(model.cache, MOI.ObjectiveSense()) + return sense == MOI.MAX_SENSE ? Inf : -Inf + end + return bound +end + +function MOI.get( + model::Optimizer, + attr::MOI.VariablePrimal, + vi::MOI.VariableIndex + ) + MOI.check_result_index_bounds(model, attr) + return model.incumbent[vi] +end + +MOI.get(model::Optimizer, ::MOI.RelativeGap) = model.relative_gap + +MOI.get(model::Optimizer, ::MOI.SolveTimeSec) = model.solve_time diff --git a/GDPOptimizer.jl/src/problem.jl b/GDPOptimizer.jl/src/problem.jl new file mode 100644 index 00000000..121f81be --- /dev/null +++ b/GDPOptimizer.jl/src/problem.jl @@ -0,0 +1,159 @@ +################################################################################ +# PROBLEM PARTITION +################################################################################ +# activation is `z` or `1 - z`; rows stay in cache space +struct _Disjunct + activation::MOI.ScalarAffineFunction{Float64} + binary::MOI.VariableIndex + active_value::Bool + functions::Vector{MOI.AbstractScalarFunction} + sets::Vector{MOI.AbstractScalarSet} +end + +# activation: 1 at top level, the parent indicator when nested +struct _Disjunction + activation::MOI.ScalarAffineFunction{Float64} + disjuncts::Vector{_Disjunct} +end + +# Constraint indices stay in cache space; the builders remap them +struct _Problem + variables::Vector{MOI.VariableIndex} + binaries::Vector{MOI.VariableIndex} + disjunctions::Vector{_Disjunction} + variable_cis::Vector{MOI.ConstraintIndex} + linear_cis::Vector{MOI.ConstraintIndex} + # stable function objects so the linearizer can cache evaluators + nonlinear_rows::Vector{Tuple{MOI.AbstractScalarFunction, + MOI.AbstractScalarSet}} + sense::MOI.OptimizationSense + objective::MOI.AbstractScalarFunction +end + +_is_linear(::Union{MOI.VariableIndex, MOI.ScalarAffineFunction}) = true +_is_linear(::MOI.AbstractScalarFunction) = false + +_to_affine(func::MOI.ScalarAffineFunction{Float64}) = func +function _to_affine(func::MOI.AbstractScalarFunction) + return convert(MOI.ScalarAffineFunction{Float64}, func) +end + +# a raw `VariableIndex` row becomes a bound and collides with the +# variable's own bounds in the subproblem +_as_row(func::MOI.VariableIndex) = _to_affine(func) +_as_row(func::MOI.AbstractScalarFunction) = func + +# demote rows the enclosing vector function promoted +_demote(func::MOI.AbstractScalarFunction) = func +function _demote(func::MOI.ScalarQuadraticFunction{Float64}) + return _try_convert(MOI.ScalarAffineFunction{Float64}, func) +end +function _demote(func::MOI.ScalarNonlinearFunction) + return _try_convert(MOI.ScalarAffineFunction{Float64}, + _try_convert(MOI.ScalarQuadraticFunction{Float64}, func)) +end + +function _try_convert(T::Type, func) + return try + convert(T, func) + catch + func + end +end + +_scalarize(func::MOI.AbstractVectorFunction) = + collect(MOI.Utilities.eachscalar(func)) + +# `z` -> (z, true), `1 - z` -> (z, false) +function _activation_binary(activation::MOI.ScalarAffineFunction{Float64}) + canonical = MOI.Utilities.canonical(activation) + if length(canonical.terms) == 1 + term = only(canonical.terms) + if term.coefficient == 1.0 && canonical.constant == 0.0 + return term.variable, true + elseif term.coefficient == -1.0 && canonical.constant == 1.0 + return term.variable, false + end + end + return error("Unsupported indicator expression `$activation`: each " * + "`DisjunctionSet` indicator must be a binary variable `z` or " * + "its complement `1 - z`.") +end + +function _parse_disjunction( + cache::_Cache, + ci::MOI.ConstraintIndex{F, DisjunctionSet} + ) where {F} + set = MOI.get(cache, MOI.ConstraintSet(), ci) + rows = _scalarize(MOI.get(cache, MOI.ConstraintFunction(), ci)) + disjunction_activation = _to_affine( + _demote(rows[activation_index(set)])) + disjuncts = _Disjunct[] + for (i, j) in enumerate(indicator_indices(set)) + activation = _to_affine(_demote(rows[j])) + binary, active_value = _activation_binary(activation) + zero_one = MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}( + binary.value) + MOI.is_valid(cache, zero_one) || error("The indicator variable " * + "of a `DisjunctionSet` disjunct must be `MOI.ZeroOne`.") + functions = MOI.AbstractScalarFunction[ + _as_row(_demote(rows[k])) for k in row_indices(set, i)] + push!(disjuncts, _Disjunct(activation, binary, active_value, + functions, set.inner_sets[i])) + end + return _Disjunction(disjunction_activation, disjuncts) +end + +# partition the cache for the LOA loop +function _build_problem(model::Optimizer) + cache = model.cache + sense = MOI.get(cache, MOI.ObjectiveSense()) + if sense == MOI.FEASIBILITY_SENSE + # feasibility model: minimize a constant zero objective + sense = MOI.MIN_SENSE + objective = MOI.ScalarAffineFunction( + MOI.ScalarAffineTerm{Float64}[], 0.0) + else + F = MOI.get(cache, MOI.ObjectiveFunctionType()) + objective = _demote(MOI.get(cache, MOI.ObjectiveFunction{F}())) + end + variables = MOI.get(cache, MOI.ListOfVariableIndices()) + disjunctions = _Disjunction[] + variable_cis = MOI.ConstraintIndex[] + linear_cis = MOI.ConstraintIndex[] + nonlinear_rows = Tuple{MOI.AbstractScalarFunction, + MOI.AbstractScalarSet}[] + for (FC, S) in MOI.get(cache, MOI.ListOfConstraintTypesPresent()) + cis = MOI.get(cache, MOI.ListOfConstraintIndices{FC, S}()) + if S === DisjunctionSet + append!(disjunctions, + _parse_disjunction(cache, ci) for ci in cis) + elseif FC === MOI.VariableIndex && S <: MOI.AbstractScalarSet + append!(variable_cis, cis) + elseif FC === MOI.ScalarAffineFunction{Float64} && + S <: MOI.AbstractScalarSet + append!(linear_cis, cis) + elseif FC <: Union{MOI.ScalarQuadraticFunction{Float64}, + MOI.ScalarNonlinearFunction} && S <: MOI.AbstractScalarSet + append!(nonlinear_rows, + (_demote(MOI.get(cache, MOI.ConstraintFunction(), ci)), + MOI.get(cache, MOI.ConstraintSet(), ci)) + for ci in cis) + else + error("GDPOptimizer does not support `$FC`-in-`$S` constraints.") + end + end + binaries = unique!([disjunct.binary for disjunction in disjunctions + for disjunct in disjunction.disjuncts]) + # non-indicator discrete variables keep their integrality; the + # nlp_solver must handle whatever remains + return _Problem(variables, binaries, disjunctions, variable_cis, + linear_cis, nonlinear_rows, sense, objective) +end + +function _disjunct_active(combination::AbstractDict, disjunct::_Disjunct) + return combination[disjunct.binary] == disjunct.active_value +end + +_is_nonlinear_disjunct(disjunct::_Disjunct) = + any(!_is_linear(func) for func in disjunct.functions) diff --git a/GDPOptimizer.jl/src/sets.jl b/GDPOptimizer.jl/src/sets.jl new file mode 100644 index 00000000..d9f95d46 --- /dev/null +++ b/GDPOptimizer.jl/src/sets.jl @@ -0,0 +1,96 @@ +################################################################################ +# DISJUNCTION SET +################################################################################ +const _SupportedInnerSet = Union{ + MOI.LessThan{Float64}, + MOI.GreaterThan{Float64}, + MOI.EqualTo{Float64}, + MOI.Interval{Float64} +} + +""" + DisjunctionSet(inner_sets::Vector{<:Vector}) + +The vector set of a disjunction with `length(inner_sets)` disjuncts. +A constraint function in this set stacks the disjunction's activation +expression, one indicator component per disjunct, and the constraint +rows of each disjunct in order: +`[p, z_1, ..., z_k, rows of disjunct 1, ..., rows of disjunct k]`, +where disjunct `i` owns `length(inner_sets[i])` rows. A point is in +the set when the indicators sum to the activation value and every +disjunct whose indicator equals 1 has its rows in their scalar sets. + +A top-level disjunction has the constant activation `1`, so exactly +one disjunct is active. A nested disjunction's activation is the +indicator expression of its parent disjunct, so it selects exactly one +disjunct while the parent is active and is vacuous otherwise. +Indicator components must be affine in binary (`MOI.ZeroOne`) +variables; a plain binary variable and its `1 - z` complement are both +valid. + +Use [`num_disjuncts`](@ref), [`activation_index`](@ref), +[`indicator_indices`](@ref), and [`row_indices`](@ref) to locate the +components of a constraint function in this set. + +Supported inner sets: `MOI.LessThan{Float64}`, +`MOI.GreaterThan{Float64}`, `MOI.EqualTo{Float64}`, +`MOI.Interval{Float64}`. +""" +struct DisjunctionSet <: MOI.AbstractVectorSet + inner_sets::Vector{Vector{MOI.AbstractScalarSet}} + function DisjunctionSet(inner_sets::Vector{<:Vector}) + isempty(inner_sets) && throw(ArgumentError( + "A `DisjunctionSet` requires at least one disjunct.")) + for sets in inner_sets, set in sets + set isa _SupportedInnerSet || throw(ArgumentError( + "Unsupported inner set `$set` in `DisjunctionSet`.")) + end + return new([MOI.AbstractScalarSet[set for set in sets] + for sets in inner_sets]) + end +end + +""" + num_disjuncts(set::DisjunctionSet) + +Return the number of disjuncts of the disjunction. +""" +num_disjuncts(set::DisjunctionSet) = length(set.inner_sets) + +""" + activation_index(set::DisjunctionSet) + +Return the component index of the disjunction's activation expression +within a constraint function in `set`. +""" +activation_index(::DisjunctionSet) = 1 + +""" + indicator_indices(set::DisjunctionSet) + +Return the component range of the disjunct indicator expressions +within a constraint function in `set`. +""" +indicator_indices(set::DisjunctionSet) = 1 .+ (1:num_disjuncts(set)) + +""" + row_indices(set::DisjunctionSet, i::Int) + +Return the component range of disjunct `i`'s constraint rows within a +constraint function in `set`. +""" +function row_indices(set::DisjunctionSet, i::Int) + offset = 1 + num_disjuncts(set) + + sum(length(set.inner_sets[j]) for j in 1:(i - 1); init = 0) + return offset .+ (1:length(set.inner_sets[i])) +end + +function MOI.dimension(set::DisjunctionSet) + return 1 + num_disjuncts(set) + sum(length, set.inner_sets) +end + +Base.copy(set::DisjunctionSet) = DisjunctionSet(set.inner_sets) + +function Base.:(==)(a::DisjunctionSet, b::DisjunctionSet) + return a.inner_sets == b.inner_sets +end diff --git a/GDPOptimizer.jl/test/loa.jl b/GDPOptimizer.jl/test/loa.jl new file mode 100644 index 00000000..bc23fa87 --- /dev/null +++ b/GDPOptimizer.jl/test/loa.jl @@ -0,0 +1,489 @@ +using JuMP +import HiGHS, Ipopt + +function _loa_optimizer(; kwargs...) + return () -> GDPO.Optimizer(; nlp_solver = Ipopt.Optimizer, + mip_solver = HiGHS.Optimizer, kwargs...) +end + +# min x with x >= 2 (disjunct 1) or x >= 5 (disjunct 2). The loop +# enumerates both combinations and keeps the incumbent at 2. +function test_linear_disjunction() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], x, x] in GDPO.DisjunctionSet([ + [MOI.GreaterThan(2.0)], [MOI.GreaterThan(5.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test primal_status(model) == MOI.FEASIBLE_POINT + @test objective_value(model) ≈ 2.0 atol = 1e-5 + @test value(x) ≈ 2.0 atol = 1e-5 + @test value(z[1]) ≈ 1.0 atol = 1e-5 + @test value(z[2]) ≈ 0.0 atol = 1e-5 + @test occursin("LOA finished", raw_status(model)) + @test solve_time(model) > 0.0 + @test dual_status(model) == MOI.NO_SOLUTION +end + +# Convex quadratic objective over linear disjuncts: y >= x or +# y >= 2 - x. The optimum sits at x = 3, y = 0 in disjunct 2. +function test_quadratic_objective() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 4) + @variable(model, 0 <= y <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], y - x, y - (2 - x)] in + GDPO.DisjunctionSet([ + [MOI.GreaterThan(0.0)], [MOI.GreaterThan(0.0)]])) + @objective(model, Min, (x - 3)^2 + y) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 0.0 atol = 1e-5 + @test value(z[2]) ≈ 1.0 atol = 1e-5 + @test value(x) ≈ 3.0 atol = 1e-4 + @test value(y) ≈ 0.0 atol = 1e-5 + # exhaustion path: the bound is the last master's, still valid + @test objective_bound(model) <= objective_value(model) + 1e-6 + @test !isnan(relative_gap(model)) +end + +# max x with (x <= 3) or (x <= 7): port of DP.jl test_loa_solve_simple. +function test_max_sense_linear() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], x, x] in GDPO.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(7.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 7.0 atol = 1e-4 + @test value(z[2]) ≈ 1.0 atol = 1e-5 +end + +# Two disjunctions: port of DP.jl test_loa_solve_two_disjunctions. +function test_two_disjunctions() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, 0 <= w <= 10) + @variable(model, zx[1:2], Bin) + @variable(model, zw[1:2], Bin) + @constraint(model, [1, zx[1], zx[2], x, x] in GDPO.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(7.0)]])) + @constraint(model, [1, zw[1], zw[2], w, w] in GDPO.DisjunctionSet([ + [MOI.LessThan(2.0)], [MOI.LessThan(5.0)]])) + @objective(model, Max, x + w) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 12.0 atol = 1e-4 +end + +# Nonlinear global x^2 <= 25 binds before the chosen disjunct: port of +# DP.jl test_loa_nonlinear_global (no Juniper needed - the layer's NLP +# has its binaries fixed). +function test_nonlinear_global() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, x^2 <= 25) + @constraint(model, [1, z[1], z[2], x, x] in GDPO.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(8.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 5.0 atol = 1e-3 + @test value(z[2]) ≈ 1.0 atol = 1e-5 + @test isfinite(relative_gap(model)) +end + +# Nonlinear global equality: one seed combination is NLP-infeasible, so +# the NLPF path supplies the linearization site. Port of DP.jl +# test_loa_nonlinear_equality_global. +function test_nonlinear_equality_global() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, x^2 == 25) + @constraint(model, [1, z[1], z[2], x, x] in GDPO.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(8.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 5.0 atol = 1e-3 + @test value(x) ≈ 5.0 atol = 1e-3 + @test value(z[2]) ≈ 1.0 atol = 1e-5 +end + +# Nonlinear equality inside a disjunct: set covering must activate it +# and the cut emits both gated directions. Port of DP.jl +# test_loa_nonlinear_equality_disjunct. +function test_nonlinear_equality_disjunct() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, 1.0z[1], 1.0z[2], 1.0x, x^2] in + GDPO.DisjunctionSet([[MOI.LessThan(3.0)], [MOI.EqualTo(64.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 8.0 atol = 1e-3 + @test value(z[2]) ≈ 1.0 atol = 1e-5 + # deterministic gap convergence: set covering visits the nonlinear + # disjunct first, then the master proves the other is worse + @test occursin("converged", raw_status(model)) + @test relative_gap(model) <= 1e-5 +end + +# Nonlinear Interval row inside a disjunct: port of DP.jl +# test_loa_nonlinear_interval_disjunct. +function test_nonlinear_interval_disjunct() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, 1.0z[1], 1.0z[2], 1.0x, x^2] in + GDPO.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.Interval(36.0, 64.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 8.0 atol = 1e-3 + @test value(z[2]) ≈ 1.0 atol = 1e-5 +end + +# A single binary drives both disjuncts through its complement: port of +# DP.jl test_loa_complement_indicator_nonlinear_disjunct. +function test_complement_indicator() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z, Bin) + @constraint(model, [1, 1.0z, 1 - z, 1.0x, x^2] in + GDPO.DisjunctionSet([[MOI.LessThan(3.0)], [MOI.EqualTo(64.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 8.0 atol = 1e-3 + @test value(z) ≈ 0.0 atol = 1e-5 +end + +# `use_nlpf = false` still solves by enumeration when the seeds are +# feasible. +function test_nlpf_disabled() + model = Model(_loa_optimizer(use_nlpf = false)) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, 1.0z[1], 1.0z[2], 1.0x, x^2] in + GDPO.DisjunctionSet([[MOI.LessThan(3.0)], [MOI.EqualTo(64.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 8.0 atol = 1e-3 +end + +# Big-M master gating solves the same instances as indicator gating +# (interval split included). +function test_bigm_master_gating() + model = Model(_loa_optimizer(master_gating = "bigm", M_value = 100.0)) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, 1.0z[1], 1.0z[2], 1.0x, 1.0x] in + GDPO.DisjunctionSet([ + [MOI.LessThan(2.0)], [MOI.Interval(3.0, 7.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 7.0 atol = 1e-4 +end + +# Integer-typed options convert at their use sites, including the +# Bool read of use_nlpf on the infeasible-seed path. +function test_integer_options() + model = Model(_loa_optimizer(use_nlpf = 0, M_value = 10^9, + time_limit = 3600)) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, 1.0z[1], 1.0z[2], 1.0x, x^2] in + GDPO.DisjunctionSet([[MOI.LessThan(3.0)], [MOI.EqualTo(64.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 8.0 atol = 1e-3 +end + +# A binary that is not a disjunction indicator keeps its integrality +# in the subproblem, which the nlp_solver must then handle (HiGHS +# both roles here since the model is linear). +function test_non_indicator_binary() + factory = () -> GDPO.Optimizer(nlp_solver = HiGHS.Optimizer, + mip_solver = HiGHS.Optimizer) + model = Model(factory) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @variable(model, w, Bin) + @constraint(model, 1.0x + w <= 8) + @constraint(model, [1, z[1], z[2], x, x] in GDPO.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(7.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 7.0 atol = 1e-5 +end + +# Globally infeasible model: the master is infeasible before any +# incumbent exists. +function test_infeasible_no_incumbent() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, 1.0x >= 20) + @constraint(model, [1, z[1], z[2], x, x] in GDPO.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(7.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.INFEASIBLE + @test primal_status(model) == MOI.NO_SOLUTION + @test result_count(model) == 0 +end + +# A zero time limit exits before any solve. +function test_time_limit() + model = Model(_loa_optimizer(time_limit = 0.0)) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], x, x] in GDPO.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(7.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.TIME_LIMIT + @test result_count(model) == 0 +end + +# Nonconvex objective: the OA bound is no certificate, so the layer +# must never report OPTIMAL and the raw status carries the gap record. +function test_nonconvex_never_optimal() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 2) + @variable(model, z[1:2], Bin) + @constraint(model, [1, 1.0z[1], 1.0z[2], 1.0x, 1.0x] in + GDPO.DisjunctionSet([[MOI.LessThan(1.0)], [MOI.LessThan(2.0)]])) + @objective(model, Min, -x^2) + optimize!(model) + @test termination_status(model) != MOI.OPTIMAL + @test termination_status(model) in + (MOI.LOCALLY_SOLVED, MOI.ITERATION_LIMIT, MOI.TIME_LIMIT) + @test primal_status(model) == MOI.FEASIBLE_POINT + @test occursin("LOA finished", raw_status(model)) +end + +# A feasibility model (no objective) minimizes a constant zero. +function test_feasibility_sense() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + set_start_value(x, 6.0) + @constraint(model, [1, z[1], z[2], x, x] in GDPO.DisjunctionSet([ + [MOI.GreaterThan(5.0)], [MOI.LessThan(1.0)]])) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test primal_status(model) == MOI.FEASIBLE_POINT + @test value(x) >= 5.0 - 1e-6 || value(x) <= 1.0 + 1e-6 +end + +# A linear Interval disjunct row reaches the master as two one-sided +# indicator constraints (Indicator{A}(Interval) has no MILP bridge). +function test_linear_interval_disjunct() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, 1.0z[1], 1.0z[2], 1.0x, 1.0x] in + GDPO.DisjunctionSet([ + [MOI.LessThan(2.0)], [MOI.Interval(3.0, 7.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 7.0 atol = 1e-4 + @test value(z[2]) ≈ 1.0 atol = 1e-5 +end + +# A second solve on the same optimizer must not leak the first solve's +# bound, gap, or incumbent. +function test_reoptimize_resets_results() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], x, x] in GDPO.DisjunctionSet([ + [MOI.GreaterThan(2.0)], [MOI.GreaterThan(5.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 2.0 atol = 1e-5 + @constraint(model, 1.0x >= 20) + optimize!(model) + @test termination_status(model) == MOI.INFEASIBLE + @test result_count(model) == 0 + @test objective_bound(model) == -Inf + @test isnan(relative_gap(model)) +end + +# An unsupported constraint type is rewritten by the JuMP bridge layer +# into supported rows instead of erroring inside the solve. +function test_bridged_vector_constraint() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1.0x - 5.0] in MOI.Nonnegatives(1)) + @constraint(model, [1, z[1], z[2], x, x] in GDPO.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(7.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 5.0 atol = 1e-4 + @test value(z[2]) ≈ 1.0 atol = 1e-5 +end + +################################################################################ +# UNIT TESTS +################################################################################ +function test_cut_term_directions() + x = MOI.VariableIndex(1) + lin = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(2.0, x)], 1.0) + @test length(GDPO._oa_cut_terms(MOI.LessThan(5.0), lin)) == 1 + @test length(GDPO._oa_cut_terms(MOI.GreaterThan(5.0), lin)) == 1 + @test length(GDPO._oa_cut_terms(MOI.EqualTo(5.0), lin)) == 2 + @test length(GDPO._oa_cut_terms(MOI.Interval(0.0, 5.0), lin)) == 2 + less = only(GDPO._oa_cut_terms(MOI.LessThan(5.0), lin)) + @test less.constant == -4.0 + greater = only(GDPO._oa_cut_terms(MOI.GreaterThan(5.0), lin)) + @test greater.constant == 4.0 + @test only(greater.terms).coefficient == -2.0 +end + +function test_activation_binary() + z = MOI.VariableIndex(1) + saf(a, c) = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(a, z)], c) + @test GDPO._activation_binary(saf(1.0, 0.0)) == (z, true) + @test GDPO._activation_binary(saf(-1.0, 1.0)) == (z, false) + @test_throws ErrorException GDPO._activation_binary(saf(2.0, 0.0)) + @test_throws ErrorException GDPO._activation_binary(saf(1.0, 3.0)) +end + +function test_sense_primitives() + @test GDPO._penalty_sign(MOI.MIN_SENSE) == 1 + @test GDPO._penalty_sign(MOI.MAX_SENSE) == -1 + @test GDPO._worst_objective(MOI.MIN_SENSE) == Inf + @test GDPO._worst_objective(MOI.MAX_SENSE) == -Inf + @test GDPO._is_better(MOI.MIN_SENSE, 1.0, 2.0) + @test GDPO._is_better(MOI.MAX_SENSE, 2.0, 1.0) + @test GDPO._gap(MOI.MIN_SENSE, 5.0, 3.0) == 2.0 + @test GDPO._gap(MOI.MAX_SENSE, 3.0, 5.0) == 2.0 +end + +function test_linearize_quadratic() + x = MOI.VariableIndex(1) + func = MOI.ScalarQuadraticFunction( + [MOI.ScalarQuadraticTerm(2.0, x, x)], + MOI.ScalarAffineTerm{Float64}[], 0.0) # x^2 + linearizer = GDPO._Linearizer() + lin = GDPO._linearize(linearizer, func, Dict(x => 3.0)) + # x^2 at x = 3: 9 + 6 (x - 3) = 6 x - 9 + @test lin.constant ≈ -9.0 + @test only(lin.terms).coefficient ≈ 6.0 + @test length(linearizer.evaluators) == 1 + GDPO._linearize(linearizer, func, Dict(x => 4.0)) + @test length(linearizer.evaluators) == 1 +end + +# Nested disjunction: the inner disjunction's activation component is +# the outer indicator, so it selects a mode only while the outer +# disjunct is active and is vacuous otherwise. +function test_nested_disjunction() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, zout[1:2], Bin) + @variable(model, zin[1:2], Bin) + @constraint(model, [1, zout[1], zout[2], x] in GDPO.DisjunctionSet([ + [MOI.LessThan(2.0)], MOI.AbstractScalarSet[]])) + @constraint(model, [1.0zout[2], zin[1], zin[2], x^2, x] in + GDPO.DisjunctionSet([[MOI.LessThan(25.0)], [MOI.LessThan(8.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 8.0 atol = 1e-3 + @test value(zout[2]) ≈ 1.0 atol = 1e-5 + @test value(zin[2]) ≈ 1.0 atol = 1e-5 +end + +# When the parent disjunct is off, the inner indicators sum to zero +# and the inner rows impose nothing. +function test_nested_disjunction_vacuous() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, zout[1:2], Bin) + @variable(model, zin[1:2], Bin) + @constraint(model, [1, zout[1], zout[2], x] in GDPO.DisjunctionSet([ + [MOI.LessThan(2.0)], MOI.AbstractScalarSet[]])) + @constraint(model, [1.0zout[2], zin[1], zin[2], x, x^2] in + GDPO.DisjunctionSet([ + [MOI.GreaterThan(5.0)], [MOI.GreaterThan(36.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 0.0 atol = 1e-4 + @test value(zout[1]) ≈ 1.0 atol = 1e-5 + @test value(zin[1]) + value(zin[2]) ≈ 0.0 atol = 1e-5 +end + +@testset "LOA loop" begin + test_linear_disjunction() + test_nested_disjunction() + test_nested_disjunction_vacuous() + test_quadratic_objective() + test_max_sense_linear() + test_two_disjunctions() + test_nonlinear_global() + test_nonlinear_equality_global() + test_nonlinear_equality_disjunct() + test_nonlinear_interval_disjunct() + test_complement_indicator() + test_nlpf_disabled() + test_bigm_master_gating() + test_integer_options() + test_non_indicator_binary() + test_infeasible_no_incumbent() + test_time_limit() + test_nonconvex_never_optimal() + test_feasibility_sense() + test_linear_interval_disjunct() + test_reoptimize_resets_results() + test_bridged_vector_constraint() +end + +@testset "LOA units" begin + test_cut_term_directions() + test_activation_binary() + test_sense_primitives() + test_linearize_quadratic() +end diff --git a/GDPOptimizer.jl/test/moi.jl b/GDPOptimizer.jl/test/moi.jl new file mode 100644 index 00000000..0dbb762e --- /dev/null +++ b/GDPOptimizer.jl/test/moi.jl @@ -0,0 +1,35 @@ +# MOI contract conformance. MOI.Test cannot build a `DisjunctionSet`, so it +# only exercises the API-contract surface (model API + optimizer attributes), +# not the LOA solving logic -- that is covered by the other test files. The +# solving families and the dual/basis attributes are out of scope by design. +import HiGHS +import Ipopt + +@testset "MOI contract (model API + attributes)" begin + optimizer = MOI.instantiate( + () -> GDPO.Optimizer( + nlp_solver = Ipopt.Optimizer, + mip_solver = HiGHS.Optimizer, + time_limit = 20.0, + ); + with_cache_type = Float64, + with_bridge_type = Float64, + ) + config = MOI.Test.Config( + atol = 1e-3, + rtol = 1e-3, + optimal_status = MOI.LOCALLY_SOLVED, # never reports OPTIMAL + exclude = Any[ + MOI.ConstraintDual, + MOI.DualObjectiveValue, + MOI.ConstraintBasisStatus, + MOI.VariableBasisStatus, + ], + ) + MOI.Test.runtests( + optimizer, + config; + include = Regex[r"test_model_", r"test_attribute_"], + warn_unsupported = false, + ) +end diff --git a/GDPOptimizer.jl/test/optimizer.jl b/GDPOptimizer.jl/test/optimizer.jl new file mode 100644 index 00000000..b01b9161 --- /dev/null +++ b/GDPOptimizer.jl/test/optimizer.jl @@ -0,0 +1,92 @@ +# A two-disjunct toy stored straight into a model: binaries z1, z2 and +# continuous x with x <= 0 (disjunct 1) or x >= 1 (disjunct 2). +function _build_toy_disjunction(model) + x = MOI.add_variable(model) + z = MOI.add_variables(model, 2) + for zi in z + MOI.add_constraint(model, zi, MOI.ZeroOne()) + end + set = GDPO.DisjunctionSet([[MOI.LessThan(0.0)], [MOI.GreaterThan(1.0)]]) + func = MOI.VectorAffineFunction( + [MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, z[1])), + MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(1.0, z[2])), + MOI.VectorAffineTerm(4, MOI.ScalarAffineTerm(1.0, x)), + MOI.VectorAffineTerm(5, MOI.ScalarAffineTerm(1.0, x))], + [1.0, 0.0, 0.0, 0.0, 0.0]) # constant activation in component 1 + ci = MOI.add_constraint(model, func, set) + return x, z, set, ci +end + +function test_optimizer_scaffold() + optimizer = GDPO.Optimizer() + @test MOI.is_empty(optimizer) + @test MOI.get(optimizer, MOI.SolverName()) == "GDPOptimizer" + @test MOI.get(optimizer, MOI.TerminationStatus()) == + MOI.OPTIMIZE_NOT_CALLED + @test MOI.get(optimizer, MOI.ResultCount()) == 0 + @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.NO_SOLUTION + @test MOI.get(optimizer, MOI.DualStatus()) == MOI.NO_SOLUTION +end + +function test_optimizer_options() + optimizer = GDPO.Optimizer(max_iter = 25) + @test MOI.get(optimizer, MOI.RawOptimizerAttribute("max_iter")) == 25 + MOI.set(optimizer, MOI.RawOptimizerAttribute("max_iter"), 5) + @test MOI.get(optimizer, MOI.RawOptimizerAttribute("max_iter")) == 5 + @test MOI.get(optimizer, MOI.RawOptimizerAttribute("M_value")) == 1e9 + bogus = MOI.RawOptimizerAttribute("bogus") + @test !MOI.supports(optimizer, bogus) + @test_throws MOI.UnsupportedAttribute MOI.get(optimizer, bogus) + @test_throws MOI.UnsupportedAttribute MOI.set(optimizer, bogus, 1) + @test_throws ArgumentError GDPO.Optimizer(bogus = 1) + MOI.set(optimizer, MOI.Silent(), true) + @test MOI.get(optimizer, MOI.Silent()) + MOI.set(optimizer, MOI.TimeLimitSec(), 100) + @test MOI.get(optimizer, MOI.TimeLimitSec()) == 100.0 + @test MOI.get(optimizer, MOI.RawOptimizerAttribute("time_limit")) == 100.0 + MOI.set(optimizer, MOI.TimeLimitSec(), nothing) + @test MOI.get(optimizer, MOI.TimeLimitSec()) === nothing +end + +function test_model_building() + optimizer = GDPO.Optimizer() + x, z, set, ci = _build_toy_disjunction(optimizer) + @test MOI.is_valid(optimizer, ci) + @test MOI.get(optimizer, MOI.NumberOfVariables()) == 3 + @test (MOI.VectorAffineFunction{Float64}, GDPO.DisjunctionSet) in + MOI.get(optimizer, MOI.ListOfConstraintTypesPresent()) + @test MOI.get(optimizer, MOI.ConstraintSet(), ci) == set + objective = MOI.ScalarAffineFunction( + [MOI.ScalarAffineTerm(1.0, x)], 0.0) + MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.set(optimizer, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + objective) + @test MOI.get(optimizer, MOI.ObjectiveSense()) == MOI.MIN_SENSE + @test !MOI.is_empty(optimizer) + MOI.empty!(optimizer) + @test MOI.is_empty(optimizer) +end + +function test_copy_to() + src = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) + x, z, set, ci = _build_toy_disjunction(src) + MOI.set(src, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.set(src, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0)) + dest = GDPO.Optimizer() + index_map = MOI.copy_to(dest, src) + @test MOI.get(dest, MOI.NumberOfVariables()) == 3 + @test (MOI.VectorAffineFunction{Float64}, GDPO.DisjunctionSet) in + MOI.get(dest, MOI.ListOfConstraintTypesPresent()) + @test MOI.get(dest, MOI.ConstraintSet(), index_map[ci]) == set + dest_func = MOI.get(dest, MOI.ConstraintFunction(), index_map[ci]) + @test dest_func.terms[1].scalar_term.variable == index_map[z[1]] +end + +@testset "Optimizer scaffold" begin + test_optimizer_scaffold() + test_optimizer_options() + test_model_building() + test_copy_to() +end diff --git a/GDPOptimizer.jl/test/runtests.jl b/GDPOptimizer.jl/test/runtests.jl new file mode 100644 index 00000000..6481e42f --- /dev/null +++ b/GDPOptimizer.jl/test/runtests.jl @@ -0,0 +1,9 @@ +using Test +import MathOptInterface as MOI +using GDPOptimizer +const GDPO = GDPOptimizer + +include("sets.jl") +include("optimizer.jl") +include("loa.jl") +include("moi.jl") diff --git a/GDPOptimizer.jl/test/sets.jl b/GDPOptimizer.jl/test/sets.jl new file mode 100644 index 00000000..307a13ad --- /dev/null +++ b/GDPOptimizer.jl/test/sets.jl @@ -0,0 +1,48 @@ +function test_set_construction() + set = GDPO.DisjunctionSet([ + [MOI.LessThan(1.0), MOI.GreaterThan(0.0)], + [MOI.EqualTo(2.0)], + ]) + @test GDPO.num_disjuncts(set) == 2 + @test MOI.dimension(set) == 1 + 2 + 3 + @test GDPO.activation_index(set) == 1 + @test GDPO.indicator_indices(set) == 2:3 + @test GDPO.row_indices(set, 1) == 4:5 + @test GDPO.row_indices(set, 2) == 6:6 + copied = copy(set) + @test copied == set + @test copied.inner_sets !== set.inner_sets +end + +function test_set_interval_inner() + set = GDPO.DisjunctionSet([[MOI.Interval(0.0, 1.0)], [MOI.LessThan(0.0)]]) + @test MOI.dimension(set) == 5 +end + +function test_set_empty_disjunct() + set = GDPO.DisjunctionSet([ + MOI.AbstractScalarSet[], + [MOI.LessThan(0.0)], + ]) + @test MOI.dimension(set) == 4 + @test GDPO.row_indices(set, 1) == 4:3 + @test isempty(GDPO.row_indices(set, 1)) + @test GDPO.row_indices(set, 2) == 4:4 +end + +function test_set_validation() + @test_throws ArgumentError GDPO.DisjunctionSet( + Vector{MOI.AbstractScalarSet}[]) + @test_throws ArgumentError GDPO.DisjunctionSet([[MOI.ZeroOne()]]) + @test_throws ArgumentError GDPO.DisjunctionSet([ + [MOI.LessThan(0.0)], + [MOI.Nonnegatives(2)], + ]) +end + +@testset "DisjunctionSet" begin + test_set_construction() + test_set_interval_inner() + test_set_empty_disjunct() + test_set_validation() +end diff --git a/Project.toml b/Project.toml index bdec9be4..e2e9354b 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,10 @@ authors = ["hdavid16 "] version = "0.6.1" [deps] +# GDPOptimizer is unregistered: it lives in [deps] (dev'ed by path) so +# the extension loads through its strong dependency. Move it to +# [weakdeps] once registered. +GDPOptimizer = "7da86d0c-51ea-4770-ae9e-8cf6d5933256" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" @@ -11,23 +15,25 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" [extensions] +GDPOptimizerDisjunctiveProgramming = "GDPOptimizer" InfiniteDisjunctiveProgramming = "InfiniteOpt" [compat] Aqua = "0.8" +GDPOptimizer = "0.1" +InfiniteOpt = "0.6.3" +Ipopt = "1.9.0" JuMP = "1.18" +Juniper = "0.9.3" Reexport = "1" julia = "1.10" -Juniper = "0.9.3" -Ipopt = "1.9.0" -InfiniteOpt = "0.6.3" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" Juniper = "2ddba703-00a4-53a7-87a5-e8b9971dde84" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["Aqua", "HiGHS", "Test", "Juniper", "Ipopt", "InfiniteOpt"] diff --git a/ext/GDPOptimizerDisjunctiveProgramming.jl b/ext/GDPOptimizerDisjunctiveProgramming.jl new file mode 100644 index 00000000..c1e776e1 --- /dev/null +++ b/ext/GDPOptimizerDisjunctiveProgramming.jl @@ -0,0 +1,76 @@ +module GDPOptimizerDisjunctiveProgramming + +import DisjunctiveProgramming as DP +import GDPOptimizer +import JuMP +import JuMP.MOI as _MOI + +################################################################################ +# DISJUNCTION SET LOWERING +################################################################################ +# One flat constraint per disjunction: [activation, indicators, +# rows]. Nested disjunctions get their own constraint with the parent +# indicator as the activation. InfiniteOpt transcribes per support. +function DP.reformulate_disjunction( + model::JuMP.AbstractModel, + disj::DP.Disjunction, + method::DP.MOIDisjunction + ) + constraints = JuMP.AbstractConstraint[] + _lower_disjunction(constraints, model, disj, 1.0) + return constraints +end + +function _lower_disjunction( + constraints::Vector{JuMP.AbstractConstraint}, + model::JuMP.AbstractModel, + disj::DP.Disjunction, + activation + ) + indicators = JuMP.AbstractJuMPScalar[] + rows = JuMP.AbstractJuMPScalar[] + inner_sets = Vector{_MOI.AbstractScalarSet}[] + for lvref in disj.indicators + indicator = 1.0 * DP.binary_variable(lvref) + push!(indicators, indicator) + sets = _MOI.AbstractScalarSet[] + for cref in get(DP._indicator_to_constraints(model), lvref, []) + constraint = JuMP.constraint_object(cref) + if constraint isa DP.Disjunction + haskey(DP._exactly1_constraints(model), cref) || error( + "`MOIDisjunction` requires nested " * + "disjunctions created with `exactly1 = true`.") + _lower_disjunction(constraints, model, constraint, + indicator) + else + push!(rows, constraint.func) + push!(sets, constraint.set) + end + end + push!(inner_sets, sets) + end + push!(constraints, JuMP.VectorConstraint( + collect(promote(1.0 * activation, indicators..., rows...)), + GDPOptimizer.DisjunctionSet(inner_sets))) + return +end + +# transcription may collapse the constant activation to a number; +# promote back to a uniform expression vector +function JuMP.build_constraint( + ::Function, + func::AbstractVector{<:JuMP.AbstractJuMPScalar}, + set::GDPOptimizer.DisjunctionSet + ) + return JuMP.VectorConstraint(func, set) +end + +function JuMP.build_constraint( + ::Function, + func::AbstractVector, + set::GDPOptimizer.DisjunctionSet + ) + return JuMP.VectorConstraint(collect(promote(func...)), set) +end + +end diff --git a/src/extension_api.jl b/src/extension_api.jl index ad3566f6..c06670fe 100644 --- a/src/extension_api.jl +++ b/src/extension_api.jl @@ -15,6 +15,29 @@ julia> InfiniteGDPModel() """ function InfiniteGDPModel end +""" + MOIDisjunction() + +Reformulation method that lowers each disjunction to a single vector +constraint in `GDPOptimizer.DisjunctionSet` so the model can be solved +by the GDPOptimizer.jl MOI solver layer (logic-based outer +approximation). Requires GDPOptimizer to be imported first and the +model's optimizer to be a `GDPOptimizer.Optimizer`. + +**Example** +```julia +julia> using DisjunctiveProgramming, GDPOptimizer, HiGHS, Ipopt + +julia> model = GDPModel(() -> GDPOptimizer.Optimizer( + nlp_solver = Ipopt.Optimizer, mip_solver = HiGHS.Optimizer)); + +julia> optimize!(model, gdp_method = MOIDisjunction()) +``` +""" +struct MOIDisjunction <: AbstractReformulationMethod end + +requires_exactly1(::MOIDisjunction) = true + """ InfiniteLogical(prefs...) diff --git a/test/aqua.jl b/test/aqua.jl index b7da33c8..ac69ab5d 100644 --- a/test/aqua.jl +++ b/test/aqua.jl @@ -1,6 +1,9 @@ using Aqua using DisjunctiveProgramming -Aqua.test_all(DisjunctiveProgramming, deps_compat = false, ambiguities = false) +# GDPOptimizer is a strong dep only because it is unregistered (its +# extension loads through it); DP.jl src never imports it directly. +Aqua.test_all(DisjunctiveProgramming, deps_compat = false, + ambiguities = false, stale_deps = (ignore = [:GDPOptimizer],)) Aqua.test_deps_compat(DisjunctiveProgramming, check_extras = (ignore=[:Test, :HiGHS],)) Aqua.test_ambiguities(DisjunctiveProgramming) \ No newline at end of file diff --git a/test/extensions/GDPOptimizerDisjunctiveProgramming.jl b/test/extensions/GDPOptimizerDisjunctiveProgramming.jl new file mode 100644 index 00000000..14eb9a7f --- /dev/null +++ b/test/extensions/GDPOptimizerDisjunctiveProgramming.jl @@ -0,0 +1,170 @@ +using HiGHS, Ipopt, Juniper, InfiniteOpt +import GDPOptimizer +import DisjunctiveProgramming as DP + +function _gdp_optimizer_factory() + return () -> GDPOptimizer.Optimizer(nlp_solver = Ipopt.Optimizer, + mip_solver = HiGHS.Optimizer) +end + +# The lowering emits one vector constraint per disjunction with the +# activation first, then the indicators, then the rows grouped by +# disjunct. +function test_disjunction_set_lowering() + model = GDPModel() + @variable(model, 0 <= x <= 10) + @variable(model, Y[1:2], Logical) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x^2 == 64, Disjunct(Y[2])) + @disjunction(model, Y) + reformulate_model(model, MOIDisjunction()) + crefs = DP._reformulation_constraints(model) + vector_crefs = filter(crefs) do cref + constraint_object(cref).set isa GDPOptimizer.DisjunctionSet + end + @test length(vector_crefs) == 1 + constraint = constraint_object(only(vector_crefs)) + set = constraint.set + @test set.inner_sets == + [[MOI.LessThan(3.0)], [MOI.EqualTo(64.0)]] + @test length(constraint.func) == MOI.dimension(set) == 5 + @test isequal_canonical(constraint.func[1], + one(constraint.func[1])) # top-level activation is the constant 1 + @test coefficient(constraint.func[2], binary_variable(Y[1])) == 1.0 +end + +# A nested disjunction lowers to its own flat constraint whose +# activation is the parent indicator; the parent carries no rows for it. +function test_lowering_nested() + model = GDPModel() + @variable(model, 0 <= x <= 10) + @variable(model, Y[1:2], Logical) + @variable(model, W[1:2], Logical) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x <= 5, Disjunct(W[1])) + @constraint(model, x <= 6, Disjunct(W[2])) + @disjunction(model, W, Disjunct(Y[2])) + @disjunction(model, Y) + reformulate_model(model, MOIDisjunction()) + objs = [constraint_object(cref) + for cref in DP._reformulation_constraints(model) + if constraint_object(cref).set isa GDPOptimizer.DisjunctionSet] + @test length(objs) == 2 + inner = only(filter(c -> c.set.inner_sets == + [[MOI.LessThan(5.0)], [MOI.LessThan(6.0)]], objs)) + outer = only(filter(c -> c.set.inner_sets == + [[MOI.LessThan(3.0)], MOI.AbstractScalarSet[]], objs)) + @test isequal_canonical(outer.func[1], one(outer.func[1])) + @test coefficient(inner.func[1], binary_variable(Y[2])) == 1.0 + @test length(outer.func) == MOI.dimension(outer.set) == 4 +end + +function test_lowering_nested_requires_exactly1() + model = GDPModel() + @variable(model, 0 <= x <= 10) + @variable(model, Y[1:2], Logical) + @variable(model, W[1:2], Logical) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x <= 5, Disjunct(W[1])) + @constraint(model, x <= 6, Disjunct(W[2])) + @disjunction(model, W, Disjunct(Y[2]), exactly1 = false) + @disjunction(model, Y) + @test_throws ErrorException reformulate_model(model, + MOIDisjunction()) +end + +# Nested solve: max x with x <= 2, or a nested mode choice between +# x^2 <= 25 and x <= 8. +function test_lowering_solve_nested() + model = GDPModel(_gdp_optimizer_factory()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, Y[1:2], Logical) + @variable(model, W[1:2], Logical) + @constraint(model, x <= 2, Disjunct(Y[1])) + @constraint(model, x^2 <= 25, Disjunct(W[1])) + @constraint(model, x <= 8, Disjunct(W[2])) + @disjunction(model, W, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, x) + optimize!(model, gdp_method = MOIDisjunction()) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 8.0 atol = 1e-3 + @test value(Y[2]) + @test value(W[2]) +end + +# The same GDP solved through a BigM reformulation and through the +# lowering into GDPOptimizer must agree: max x with (x <= 3) or +# (x <= 7). +function test_lowering_solve_linear() + model = GDPModel(_gdp_optimizer_factory()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, Y[1:2], Logical) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x <= 7, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, x) + optimize!(model, gdp_method = MOIDisjunction()) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 7.0 atol = 1e-4 + @test value(x) ≈ 7.0 atol = 1e-4 + @test value(Y[2]) + + reference = GDPModel(HiGHS.Optimizer) + set_silent(reference) + @variable(reference, 0 <= x2 <= 10) + @variable(reference, Y2[1:2], Logical) + @constraint(reference, x2 <= 3, Disjunct(Y2[1])) + @constraint(reference, x2 <= 7, Disjunct(Y2[2])) + @disjunction(reference, Y2) + @objective(reference, Max, x2) + optimize!(reference, gdp_method = BigM()) + @test objective_value(model) ≈ objective_value(reference) atol = 1e-6 +end + +# Nonlinear disjunct row: max x with (x <= 3) or (x^2 == 64); the +# unique optimum 8 is checked directly +function test_lowering_solve_nonlinear() + model = GDPModel(_gdp_optimizer_factory()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, Y[1:2], Logical) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x^2 == 64, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, x) + optimize!(model, gdp_method = MOIDisjunction()) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 8.0 atol = 1e-3 + @test value(Y[2]) +end + +# An InfiniteGDPModel lowers through the same method: the disjunction +# transcribes to one DisjunctionSet constraint per support. +function test_lowering_infinite() + model = InfiniteGDPModel(_gdp_optimizer_factory()) + set_silent(model) + @infinite_parameter(model, t in [0, 1], num_supports = 3) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y[1:2], InfiniteLogical(t)) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x >= 5, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, integral(x, t)) + optimize!(model, gdp_method = MOIDisjunction()) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 10.0 atol = 1e-4 + @test all(value(x) .>= 5.0 .- 1e-4) +end + +@testset "GDPOptimizer lowering" begin + test_disjunction_set_lowering() + test_lowering_nested() + test_lowering_nested_requires_exactly1() + test_lowering_solve_linear() + test_lowering_solve_nonlinear() + test_lowering_solve_nested() + test_lowering_infinite() +end diff --git a/test/runtests.jl b/test/runtests.jl index 06e8813d..987c7a3f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -24,4 +24,5 @@ include("constraints/disjunction.jl") include("print.jl") include("solve.jl") include("extensions/InfiniteDisjunctiveProgramming.jl") +include("extensions/GDPOptimizerDisjunctiveProgramming.jl")