diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000..5f06744 --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,24 @@ +# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates +version: 2 +updates: + - package-ecosystem: "github-actions" + directory: "/" + schedule: + interval: "weekly" + # To group all GitHub Actions updates into a single PR, uncomment the following: + # groups: + # github-actions: + # patterns: + # - "*" + - package-ecosystem: "julia" + directories: + - "/" + - "/docs" + - "/test" + schedule: + interval: "weekly" + # To group all Julia dependency updates into a single PR, uncomment the following: + # groups: + # julia-dependencies: + # patterns: + # - "*" diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml new file mode 100644 index 0000000..2e521a4 --- /dev/null +++ b/.github/workflows/CI.yml @@ -0,0 +1,72 @@ +name: CI +on: + push: + branches: + - main + tags: ['*'] + pull_request: + workflow_dispatch: +concurrency: + # Skip intermediate builds: always. + # Cancel intermediate builds: only if it is a pull request build. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} +jobs: + test: + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} + runs-on: ${{ matrix.os }} + timeout-minutes: 60 + permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created + actions: write + contents: read + strategy: + fail-fast: false + matrix: + version: + - '1.10' + - '1.11' + - 'pre' + os: + - ubuntu-latest + arch: + - x64 + steps: + - uses: actions/checkout@v6 + - uses: julia-actions/setup-julia@v3 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v2 + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-runtest@v1 + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v6 + with: + files: lcov.info + token: ${{ secrets.CODECOV_TOKEN }} + fail_ci_if_error: false + docs: + name: Documentation + runs-on: ubuntu-latest + permissions: + actions: write # needed to allow julia-actions/cache to proactively delete old caches that it has created + contents: write + statuses: write + steps: + - uses: actions/checkout@v6 + - uses: julia-actions/setup-julia@v3 + with: + version: '1' + - uses: julia-actions/cache@v2 + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-docdeploy@v1 + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} + - name: Run doctests + shell: julia --project=docs --color=yes {0} + run: | + using Documenter: DocMeta, doctest + using DisjunctiveAlgorithms + DocMeta.setdocmeta!(DisjunctiveAlgorithms, :DocTestSetup, :(using DisjunctiveAlgorithms); recursive=true) + doctest(DisjunctiveAlgorithms) diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml new file mode 100644 index 0000000..ff25ebc --- /dev/null +++ b/.github/workflows/CompatHelper.yml @@ -0,0 +1,18 @@ +name: CompatHelper +on: + schedule: + - cron: 0 0 * * * + workflow_dispatch: +jobs: + CompatHelper: + runs-on: ubuntu-latest + steps: + - name: Install CompatHelper + run: using Pkg; Pkg.add("CompatHelper") + shell: julia --color=yes {0} + - name: Run CompatHelper + run: using CompatHelper; CompatHelper.main() + shell: julia --color=yes {0} + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }} diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml new file mode 100644 index 0000000..3f42eb1 --- /dev/null +++ b/.github/workflows/TagBot.yml @@ -0,0 +1,19 @@ +name: TagBot +on: + issue_comment: + types: + - created + workflow_dispatch: +jobs: + TagBot: + if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot' + runs-on: ubuntu-latest + steps: + - uses: JuliaRegistries/TagBot@v1 + with: + token: ${{ secrets.GITHUB_TOKEN }} + # For commits that modify workflow files: SSH key enables tagging, but + # releases require manual creation. For full automation of such commits, + # use a PAT with `workflow` scope instead of GITHUB_TOKEN. + # See: https://github.com/JuliaRegistries/TagBot#commits-that-modify-workflow-files + ssh: ${{ secrets.DOCUMENTER_KEY }} diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..ff640ed --- /dev/null +++ b/Project.toml @@ -0,0 +1,23 @@ +name = "DisjunctiveAlgorithms" +uuid = "af3b8e3a-4098-461f-b6c0-2c9706adb26e" +authors = ["Daniel Nguyen"] +version = "0.1.0" + +[deps] +DisjunctiveProgramming = "0d27d021-0159-4c7d-b4a7-9ccb5d9366cf" +MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" + +[compat] +DisjunctiveProgramming = "0.6" +MathOptInterface = "1" +julia = "1.10" + +[extras] +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" +InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["HiGHS", "InfiniteOpt", "Ipopt", "JuMP", "Test"] diff --git a/README.md b/README.md index d6ac54f..1d82799 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,40 @@ # DisjunctiveAlgorithms.jl -An optimizer suite for generalized disjunctive programming. + +[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://infiniteopt.github.io/DisjunctiveAlgorithms.jl/stable/) +[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://infiniteopt.github.io/DisjunctiveAlgorithms.jl/dev/) +[![Build Status](https://github.com/infiniteopt/DisjunctiveAlgorithms.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/infiniteopt/DisjunctiveAlgorithms.jl/actions/workflows/CI.yml?query=branch%3Amain) +[![Coverage](https://codecov.io/gh/infiniteopt/DisjunctiveAlgorithms.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/infiniteopt/DisjunctiveAlgorithms.jl) + +An optimizer suite for generalized disjunctive programming (GDP). + +DisjunctiveAlgorithms.jl is an MOI-layer solver for models that contain +disjunctions encoded as vector constraints in +`DisjunctiveProgramming.DisjunctionSet`. Disjunction-aware algorithms +(currently logic-based outer approximation) solve the model by +dispatching subproblems to user-provided MIP and NLP solvers. + +The design follows +[MultiObjectiveAlgorithms.jl](https://github.com/jump-dev/MultiObjectiveAlgorithms.jl): +one `Optimizer` that wraps inner solvers, with the algorithm and its +options selected through optimizer attributes. + +## Usage with DisjunctiveProgramming.jl + +```julia +using DisjunctiveProgramming, DisjunctiveAlgorithms, HiGHS, Ipopt +import DisjunctiveAlgorithms as DA + +model = GDPModel(() -> DA.Optimizer(nlp_solver = Ipopt.Optimizer, + mip_solver = HiGHS.Optimizer)) +@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()) +``` + +`MOIDisjunction()` lowers each disjunction to a single +`DisjunctionSet` constraint that this package consumes directly; no +Big-M or Hull reformulation is performed on the modeling side. diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..abda110 --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,6 @@ +[deps] +DisjunctiveAlgorithms = "af3b8e3a-4098-461f-b6c0-2c9706adb26e" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" + +[compat] +Documenter = "1.17.0" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..20dc440 --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,22 @@ +using DisjunctiveAlgorithms +using Documenter + +DocMeta.setdocmeta!(DisjunctiveAlgorithms, :DocTestSetup, :(using DisjunctiveAlgorithms); recursive=true) + +makedocs(; + modules=[DisjunctiveAlgorithms], + authors="Daniel Nguyen", + sitename="DisjunctiveAlgorithms.jl", + format=Documenter.HTML(; + canonical="https://infiniteopt.github.io/DisjunctiveAlgorithms.jl", + edit_link="main", + assets=String[], + ), + pages=[ + "Home" => "index.md", + ], +) + +deploydocs(; + repo="github.com/infiniteopt/DisjunctiveAlgorithms.jl", +) diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..4c8e161 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,14 @@ +```@meta +CurrentModule = DisjunctiveAlgorithms +``` + +# DisjunctiveAlgorithms + +Documentation for [DisjunctiveAlgorithms](https://github.com/infiniteopt/DisjunctiveAlgorithms.jl). + +```@index +``` + +```@autodocs +Modules = [DisjunctiveAlgorithms] +``` diff --git a/src/DisjunctiveAlgorithms.jl b/src/DisjunctiveAlgorithms.jl new file mode 100644 index 0000000..6831d27 --- /dev/null +++ b/src/DisjunctiveAlgorithms.jl @@ -0,0 +1,15 @@ +module DisjunctiveAlgorithms + +import MathOptInterface as MOI +import DisjunctiveProgramming: DisjunctionSet, num_disjuncts, + activation_index, indicator_indices, row_indices, + _SupportedInnerSet + +include("optimizer.jl") +include("problem.jl") +include("master.jl") +include("nlp.jl") +include("cuts.jl") +include("loa.jl") + +end diff --git a/src/cuts.jl b/src/cuts.jl new file mode 100644 index 0000000..079cc2e --- /dev/null +++ b/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/src/loa.jl b/src/loa.jl new file mode 100644 index 0000000..d68fdd6 --- /dev/null +++ b/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/src/master.jl b/src/master.jl new file mode 100644 index 0000000..7e35d78 --- /dev/null +++ b/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("DisjunctiveAlgorithms 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/src/nlp.jl b/src/nlp.jl new file mode 100644 index 0000000..eaff857 --- /dev/null +++ b/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/src/optimizer.jl b/src/optimizer.jl new file mode 100644 index 0000000..b55dd60 --- /dev/null +++ b/src/optimizer.jl @@ -0,0 +1,380 @@ +################################################################################ +# 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) = "DisjunctiveAlgorithms" +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 + +# The LOA loop consumes warm starts, so accept them even though the +# honest inner model does not store them (the fallback does). +function MOI.supports( + ::Optimizer, + ::MOI.VariablePrimalStart, + ::Type{MOI.VariableIndex} + ) + return true +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/src/problem.jl b/src/problem.jl new file mode 100644 index 0000000..ebce607 --- /dev/null +++ b/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("DisjunctiveAlgorithms 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/test/integration.jl b/test/integration.jl new file mode 100644 index 0000000..41067e0 --- /dev/null +++ b/test/integration.jl @@ -0,0 +1,121 @@ +using DisjunctiveProgramming, HiGHS, Ipopt, InfiniteOpt + +function _optimizer_factory() + return () -> DA.Optimizer(nlp_solver = Ipopt.Optimizer, + mip_solver = HiGHS.Optimizer) +end + +# The same GDP solved through a BigM reformulation and through the +# lowering into DisjunctiveAlgorithms must agree: max x with +# (x <= 3) or (x <= 7). +function test_lowering_solve_linear() + model = GDPModel(_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(_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 + +# 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(_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 + +# An InfiniteGDPModel lowers through the same method: the disjunction +# transcribes to one DisjunctionSet constraint per support. +function test_lowering_infinite() + model = InfiniteGDPModel(_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 + +# A nested infinite disjunction transcribes with the parent indicator +# as its activation, so the rebuilt constraint function is a uniform +# vector of scalar expressions (no constant to promote). +function test_lowering_infinite_nested() + model = InfiniteGDPModel(_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)) + @variable(model, W[1:2], InfiniteLogical(t)) + @constraint(model, x <= 2, Disjunct(Y[1])) + @constraint(model, x >= 5, Disjunct(W[1])) + @constraint(model, x >= 3, Disjunct(W[2])) + @disjunction(model, W, 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 +end + +@testset "DisjunctiveProgramming integration" begin + test_lowering_solve_linear() + test_lowering_solve_nonlinear() + test_lowering_solve_nested() + test_lowering_infinite() + test_lowering_infinite_nested() +end diff --git a/test/loa.jl b/test/loa.jl new file mode 100644 index 0000000..5102e92 --- /dev/null +++ b/test/loa.jl @@ -0,0 +1,791 @@ +using JuMP +import HiGHS, Ipopt + +function _loa_optimizer(; kwargs...) + return () -> DA.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 DA.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 + DA.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 DA.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 DA.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(7.0)]])) + @constraint(model, [1, zw[1], zw[2], w, w] in DA.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 DA.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 DA.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 + DA.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 + DA.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 + DA.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 + DA.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 + DA.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 + DA.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 = () -> DA.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 DA.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 DA.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 DA.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 + DA.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 DA.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 + DA.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 DA.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 DA.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 + +# Genuinely nonlinear (non-quadratic) global rows stay +# `ScalarNonlinearFunction`s through demotion and the linearizer walks +# their argument trees (variable and affine arguments). +function test_nonlinear_exp_global() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, exp(x) <= exp(5.0)) + @constraint(model, exp(x + 1.0) <= exp(6.5)) + @constraint(model, [1, z[1], z[2], x, x] in DA.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 +end + +# `use_nlpf = false` with an infeasible combination: the solve keeps +# only the no-good cut and still finishes from the other disjunct. +function test_nlpf_disabled_infeasible_combination() + model = Model(_loa_optimizer(use_nlpf = false)) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, 1.0x >= 5) + @constraint(model, [1, 1.0z[1], 1.0z[2], x^2, x^2] in + DA.DisjunctionSet([[MOI.LessThan(9.0)], [MOI.LessThan(64.0)]])) + @objective(model, Min, 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 +end + +# The NLPF slacks GreaterThan rows and copies the global linear rows +# when restoring an infeasible combination. +function test_nlpf_slacked_rows() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, 1.0x >= 4) + @constraint(model, x^2 >= 25) + @constraint(model, [1, 1.0z[1], 1.0z[2], x^2, x^2] in + DA.DisjunctionSet([[MOI.LessThan(9.0)], [MOI.LessThan(64.0)]])) + @objective(model, Min, 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 +end + +# A quadratic-typed global row with no quadratic terms demotes to +# affine and lands in the master directly. +function test_demoted_affine_global() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, zero(QuadExpr) + 1.0x <= 6) + @constraint(model, [1, z[1], z[2], x, x] in DA.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) ≈ 6.0 atol = 1e-4 +end + +# An all-variable disjunction arrives as `VectorOfVariables`; its raw +# variable rows are promoted to affine rows rather than bounds. +function test_vector_of_variables_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 DA.DisjunctionSet([ + [MOI.LessThan(2.0)], MOI.AbstractScalarSet[]])) + @constraint(model, [zout[2], zin[1], zin[2], x, x] in + DA.DisjunctionSet([[MOI.LessThan(5.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-4 + @test value(zout[2]) ≈ 1.0 atol = 1e-5 + @test value(zin[2]) ≈ 1.0 atol = 1e-5 +end + +# Zero iteration budgets exit before any solve, without an incumbent. +function test_iteration_limit_no_incumbent() + model = Model(_loa_optimizer(set_cover_max_iter = 0, max_iter = 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 DA.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(7.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.ITERATION_LIMIT + @test result_count(model) == 0 +end + +# `max_iter = 0` keeps the set-covering incumbent but produces no +# master bound. +function test_iteration_limit_with_incumbent() + model = Model(_loa_optimizer(max_iter = 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 DA.DisjunctionSet([ + [MOI.GreaterThan(2.0)], [MOI.GreaterThan(5.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.ITERATION_LIMIT + @test primal_status(model) == MOI.FEASIBLE_POINT + @test occursin("limit hit", raw_status(model)) + @test occursin("no bound", raw_status(model)) + @test objective_bound(model) == -Inf +end + +# An unbounded master (no OA cuts yet bound alpha_oa) surfaces its +# status instead of looping. +function test_master_abnormal_status() + model = Model(_loa_optimizer(set_cover_max_iter = 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 DA.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(7.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.OTHER_LIMIT + @test result_count(model) == 0 + @test occursin("master solve finished", raw_status(model)) +end + +################################################################################ +# MOCK SOLVER +################################################################################ +# Delegates every MOI call to a wrapped optimizer, but sleeps +# `sleep_time` seconds in each solve and, from solve `fail_from` on, +# skips the inner solve and reports `fail_status` with no solution. +# Deterministic triggers for the deadline and abnormal-master paths. +mutable struct MockSolver <: MOI.AbstractOptimizer + inner::MOI.AbstractOptimizer + sleep_time::Float64 + fail_from::Int + fail_status::MOI.TerminationStatusCode + solves::Int + failing::Bool +end + +function MockSolver( + factory; + sleep_time::Float64 = 0.0, + fail_from::Int = typemax(Int), + fail_status::MOI.TerminationStatusCode = MOI.NODE_LIMIT + ) + return MockSolver(MOI.instantiate(factory), sleep_time, fail_from, + fail_status, 0, false) +end + +function MOI.optimize!(model::MockSolver) + model.solves += 1 + model.sleep_time > 0 && sleep(model.sleep_time) + model.failing = model.solves >= model.fail_from + model.failing || MOI.optimize!(model.inner) + return +end + +const _WrappedAttr = Union{MOI.AbstractModelAttribute, + MOI.AbstractOptimizerAttribute} +const _WrappedIndexAttr = Union{MOI.AbstractVariableAttribute, + MOI.AbstractConstraintAttribute} +const _WrappedIndex = Union{MOI.VariableIndex, MOI.ConstraintIndex} + +function MOI.get(model::MockSolver, attr::_WrappedAttr) + if model.failing + attr isa MOI.TerminationStatus && return model.fail_status + attr isa MOI.PrimalStatus && return MOI.NO_SOLUTION + end + return MOI.get(model.inner, attr) +end + +MOI.is_empty(model::MockSolver) = MOI.is_empty(model.inner) +MOI.empty!(model::MockSolver) = MOI.empty!(model.inner) +MOI.supports_incremental_interface(::MockSolver) = true +MOI.copy_to(model::MockSolver, src::MOI.ModelLike) = + MOI.copy_to(model.inner, src) +MOI.add_variable(model::MockSolver) = MOI.add_variable(model.inner) +MOI.delete(model::MockSolver, index) = MOI.delete(model.inner, index) +MOI.is_valid(model::MockSolver, index) = MOI.is_valid(model.inner, index) + +function MOI.add_constraint( + model::MockSolver, + func::MOI.AbstractFunction, + set::MOI.AbstractSet + ) + return MOI.add_constraint(model.inner, func, set) +end + +function MOI.supports_constraint( + model::MockSolver, + F::Type{<:MOI.AbstractFunction}, + S::Type{<:MOI.AbstractSet} + ) + return MOI.supports_constraint(model.inner, F, S) +end + +MOI.supports(model::MockSolver, attr::_WrappedAttr) = + MOI.supports(model.inner, attr) + +MOI.set(model::MockSolver, attr::_WrappedAttr, value) = + MOI.set(model.inner, attr, value) + +function MOI.supports( + model::MockSolver, + attr::_WrappedIndexAttr, + I::Type{<:_WrappedIndex} + ) + return MOI.supports(model.inner, attr, I) +end + +function MOI.get( + model::MockSolver, + attr::_WrappedIndexAttr, + index::_WrappedIndex + ) + return MOI.get(model.inner, attr, index) +end + +function MOI.set( + model::MockSolver, + attr::_WrappedIndexAttr, + index::_WrappedIndex, + value + ) + return MOI.set(model.inner, attr, index, value) +end + +function _mock_time_limit_model(limit::Float64; kwargs...) + factory = () -> DA.Optimizer( + nlp_solver = () -> MockSolver(Ipopt.Optimizer; kwargs...), + mip_solver = HiGHS.Optimizer, + iteration_time_limit = limit) + model = Model(factory) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.GreaterThan(2.0)], [MOI.GreaterThan(5.0)]])) + @objective(model, Min, x) + return model +end + +# The loop deadline passes right after the covering incumbent: the +# mock NLP sleeps past `iteration_time_limit`, so the main loop never +# starts and the incumbent is reported against the time limit. The +# first solve warms up the mock stack so compilation latency cannot +# eat the deadline before the covering pass runs. +function test_time_limit_with_incumbent() + warmup = _mock_time_limit_model(Inf) + optimize!(warmup) + @test objective_value(warmup) ≈ 2.0 atol = 1e-5 + model = _mock_time_limit_model(0.5, sleep_time = 1.5) + optimize!(model) + @test termination_status(model) == MOI.TIME_LIMIT + @test primal_status(model) == MOI.FEASIBLE_POINT + @test occursin("limit hit", raw_status(model)) +end + +# The master finishes abnormally after an incumbent exists: the mock +# master reports a node limit on its second solve. +function test_master_abnormal_status_with_incumbent() + factory = () -> DA.Optimizer( + nlp_solver = Ipopt.Optimizer, + mip_solver = () -> MockSolver(HiGHS.Optimizer, fail_from = 2)) + model = Model(factory) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.GreaterThan(2.0)], [MOI.GreaterThan(5.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.OTHER_LIMIT + @test primal_status(model) == MOI.FEASIBLE_POINT + @test occursin("limit hit", raw_status(model)) +end + +################################################################################ +# UNIT TESTS +################################################################################ +function test_cut_term_directions() + x = MOI.VariableIndex(1) + lin = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(2.0, x)], 1.0) + @test length(DA._oa_cut_terms(MOI.LessThan(5.0), lin)) == 1 + @test length(DA._oa_cut_terms(MOI.GreaterThan(5.0), lin)) == 1 + @test length(DA._oa_cut_terms(MOI.EqualTo(5.0), lin)) == 2 + @test length(DA._oa_cut_terms(MOI.Interval(0.0, 5.0), lin)) == 2 + less = only(DA._oa_cut_terms(MOI.LessThan(5.0), lin)) + @test less.constant == -4.0 + greater = only(DA._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 DA._activation_binary(saf(1.0, 0.0)) == (z, true) + @test DA._activation_binary(saf(-1.0, 1.0)) == (z, false) + @test_throws ErrorException DA._activation_binary(saf(2.0, 0.0)) + @test_throws ErrorException DA._activation_binary(saf(1.0, 3.0)) +end + +function test_sense_primitives() + @test DA._penalty_sign(MOI.MIN_SENSE) == 1 + @test DA._penalty_sign(MOI.MAX_SENSE) == -1 + @test DA._worst_objective(MOI.MIN_SENSE) == Inf + @test DA._worst_objective(MOI.MAX_SENSE) == -Inf + @test DA._is_better(MOI.MIN_SENSE, 1.0, 2.0) + @test DA._is_better(MOI.MAX_SENSE, 2.0, 1.0) + @test DA._gap(MOI.MIN_SENSE, 5.0, 3.0) == 2.0 + @test DA._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 = DA._Linearizer() + lin = DA._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 + DA._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 DA.DisjunctionSet([ + [MOI.LessThan(2.0)], MOI.AbstractScalarSet[]])) + @constraint(model, [1.0zout[2], zin[1], zin[2], x^2, x] in + DA.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 DA.DisjunctionSet([ + [MOI.LessThan(2.0)], MOI.AbstractScalarSet[]])) + @constraint(model, [1.0zout[2], zin[1], zin[2], x, x^2] in + DA.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_exp_global() + test_nonlinear_equality_global() + test_nonlinear_equality_disjunct() + test_nonlinear_interval_disjunct() + test_complement_indicator() + test_nlpf_disabled() + test_nlpf_disabled_infeasible_combination() + test_nlpf_slacked_rows() + test_demoted_affine_global() + test_vector_of_variables_disjunction() + test_bigm_master_gating() + test_integer_options() + test_non_indicator_binary() + test_infeasible_no_incumbent() + test_time_limit() + test_iteration_limit_no_incumbent() + test_iteration_limit_with_incumbent() + test_master_abnormal_status() + test_time_limit_with_incumbent() + test_master_abnormal_status_with_incumbent() + 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/test/moi.jl b/test/moi.jl new file mode 100644 index 0000000..25fd8b1 --- /dev/null +++ b/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( + () -> DA.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/test/optimizer.jl b/test/optimizer.jl new file mode 100644 index 0000000..e5f004c --- /dev/null +++ b/test/optimizer.jl @@ -0,0 +1,127 @@ +# 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 = DA.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 = DA.Optimizer() + @test MOI.is_empty(optimizer) + @test MOI.get(optimizer, MOI.SolverName()) == "DisjunctiveAlgorithms" + @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 = DA.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 DA.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 = DA.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}, DA.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 = DA.Optimizer() + index_map = MOI.copy_to(dest, src) + @test MOI.get(dest, MOI.NumberOfVariables()) == 3 + @test (MOI.VectorAffineFunction{Float64}, DA.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 + +# The incremental MOI surface forwards to the cache: names, starts, +# constrained variables, and deletion. +function test_moi_forwarding() + optimizer = DA.Optimizer() + @test MOI.supports_incremental_interface(optimizer) + x = MOI.add_variable(optimizer) + MOI.set(optimizer, MOI.VariableName(), x, "x") + @test MOI.get(optimizer, MOI.VariableIndex, "x") == x + MOI.set(optimizer, MOI.VariablePrimalStart(), x, 2.5) + @test MOI.get(optimizer, MOI.VariablePrimalStart(), x) == 2.5 + y, ci_y = MOI.add_constrained_variables(optimizer, MOI.Nonnegatives(2)) + @test length(y) == 2 + @test MOI.is_valid(optimizer, ci_y) + func = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0) + ci = MOI.add_constraint(optimizer, func, MOI.LessThan(1.0)) + MOI.set(optimizer, MOI.ConstraintName(), ci, "c") + @test MOI.get(optimizer, typeof(ci), "c") == ci + MOI.delete(optimizer, ci) + @test !MOI.is_valid(optimizer, ci) + @test !MOI.supports_add_constrained_variables(optimizer, + DA.DisjunctionSet) +end + +# A constraint type the partition cannot place errors at solve time. +function test_unsupported_constraint_type() + optimizer = DA.Optimizer() + x = MOI.add_variable(optimizer) + func = MOI.VectorAffineFunction( + [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x))], [0.0]) + MOI.add_constraint(optimizer, func, MOI.Nonnegatives(1)) + @test_throws ErrorException MOI.optimize!(optimizer) +end + +@testset "Optimizer scaffold" begin + test_optimizer_scaffold() + test_optimizer_options() + test_model_building() + test_copy_to() + test_moi_forwarding() + test_unsupported_constraint_type() +end diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..f659abc --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,9 @@ +using Test +import MathOptInterface as MOI +using DisjunctiveAlgorithms +const DA = DisjunctiveAlgorithms + +include("optimizer.jl") +include("loa.jl") +include("moi.jl") +include("integration.jl")