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..993afe0 --- /dev/null +++ b/.github/workflows/CI.yml @@ -0,0 +1,84 @@ +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 + # DisjunctionSet ships in an unreleased DisjunctiveProgramming: + # add it from the gdp_optimizer branch until a DP release + - name: Add unreleased DisjunctiveProgramming + run: julia --project=. -e 'using Pkg; Pkg.add(PackageSpec(url="https://github.com/dnguyen227/DisjunctiveProgramming.jl", rev="gdp_optimizer"))' + shell: bash + - 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 + # DisjunctionSet ships in an unreleased DisjunctiveProgramming: + # add it from the gdp_optimizer branch until a DP release + - name: Add unreleased DisjunctiveProgramming + run: | + julia --project=. -e 'using Pkg; Pkg.add(PackageSpec(url="https://github.com/dnguyen227/DisjunctiveProgramming.jl", rev="gdp_optimizer"))' + julia --project=docs -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.add(PackageSpec(url="https://github.com/dnguyen227/DisjunctiveProgramming.jl", rev="gdp_optimizer")); Pkg.instantiate()' + shell: bash + - 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..e922c08 --- /dev/null +++ b/Project.toml @@ -0,0 +1,24 @@ +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" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[compat] +DisjunctiveProgramming = "0.6, 0.7" +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..143013b 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,53 @@ # 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. + +## Installation + +DisjunctiveAlgorithms needs the `DisjunctionSet` set and the `Direct()` +reformulation, which live on DisjunctiveProgramming's `gdp_optimizer` +branch and are not in a registered release yet. Until they ship, +install the branch first: + +```julia +import Pkg +Pkg.add(url = "https://github.com/dnguyen227/DisjunctiveProgramming.jl", + rev = "gdp_optimizer") +Pkg.add(url = "https://github.com/infiniteopt/DisjunctiveAlgorithms.jl") +``` + +## Usage with DisjunctiveProgramming.jl + +```julia +using DisjunctiveProgramming, DisjunctiveAlgorithms, HiGHS, Ipopt +import DisjunctiveAlgorithms as DA + +model = GDPModel(() -> DA.Optimizer(Ipopt.Optimizer, 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 = Direct()) +``` + +`Direct()` 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..3fbe4b2 --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,3 @@ +[deps] +DisjunctiveAlgorithms = "af3b8e3a-4098-461f-b6c0-2c9706adb26e" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" 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..60155c3 --- /dev/null +++ b/src/DisjunctiveAlgorithms.jl @@ -0,0 +1,16 @@ +module DisjunctiveAlgorithms + +import MathOptInterface as MOI +import Random +import DisjunctiveProgramming: DisjunctionSet, activation_index, + indicator_indices, row_indices, SupportedInnerSet + +include("optimizer.jl") +include("problem.jl") +include("algorithms/LOA/master.jl") +include("combination_sources.jl") +include("algorithms/LOA/nlp.jl") +include("algorithms/LOA/cuts.jl") +include("algorithms/LOA.jl") + +end diff --git a/src/algorithms/LOA.jl b/src/algorithms/LOA.jl new file mode 100644 index 0000000..93b4867 --- /dev/null +++ b/src/algorithms/LOA.jl @@ -0,0 +1,502 @@ +################################################################################ +# LOGIC-BASED OUTER APPROXIMATION +################################################################################ +""" + LOA() + +`LOA` implements logic-based outer approximation (Turkay and +Grossmann 1996): a MILP master keeps the linear rows exactly, gates +each disjunct's rows on its indicator, and accumulates OA cuts from +NLP subproblems solved at fixed indicator combinations. The +set-covering pass of Turkay and Grossmann (1996, Appendix D), +restricted to the nonlinear disjuncts as in GDPopt (Chen, Johnson, +Siirola and Grossmann 2018), initializes it with combinations that +activate every nonlinear disjunct once; [`SetCoverProblem`](@ref) +picks the covering problem. The OA bound is only +valid on convex problems, so convergence reports `LOCALLY_SOLVED`, +never `OPTIMAL`. + +## Supported optimizer attributes + +- [`NumIterationLimit`](@ref) +- [`SetCoverIterationLimit`](@ref) +- [`SetCoverProblem`](@ref) +- [`M_Value`](@ref) +- [`MasterReformulation`](@ref) +- [`MaxSlack`](@ref) +- [`OASlack`](@ref) +- [`UseNLPF`](@ref) +- [`ConvergenceTolerance`](@ref) +- [`SlackTolerance`](@ref) +- [`IterationTimeLimit`](@ref) +- [`MultiGenerationSize`](@ref) +- [`CombinationSource`](@ref) +- [`SubproblemMethod`](@ref) +""" +mutable struct LOA <: AbstractAlgorithm + num_iteration_limit::Union{Nothing, Int} + set_cover_iteration_limit::Union{Nothing, Int} + set_cover_problem::Union{Nothing, AbstractSetCover} + m_value::Union{Nothing, Float64} + master_reformulation::Union{Nothing, String} + max_slack::Union{Nothing, Float64} + oa_slack::Union{Nothing, Float64} + use_nlpf::Union{Nothing, Bool} + convergence_tolerance::Union{Nothing, Float64} + slack_tolerance::Union{Nothing, Float64} + iteration_time_limit::Union{Nothing, Float64} + multi_generation_size::Union{Nothing, Int} + combination_source::Any + subproblem_method::Any + + LOA() = new(nothing, nothing, nothing, nothing, nothing, nothing, + nothing, nothing, nothing, nothing, nothing, nothing, nothing, + nothing) +end + +_default(::Algorithm) = LOA() + +# Attribute values that would explode mid-solve fail at set time +# instead; attributes without a method to extend are unvalidated. +_validate(::AbstractAlgorithmAttribute, value) = nothing + +function _validate(::MasterReformulation, value) + value in ("indicator", "bigm") || error( + "`MasterReformulation` must be \"indicator\" or \"bigm\" " * + "(got `$value`).") + return nothing +end + +function _validate(::SetCoverProblem, value) + value isa AbstractSetCover || error("`SetCoverProblem` must be an " * + "`AbstractSetCover`, e.g. `LogicOnlyCover()` (got `$value`).") + return nothing +end +function _validate(::OASlack, value) + value > 0 || + error("`OASlack` must be positive (got `$value`).") + return nothing +end + +function _validate(::MultiGenerationSize, value) + value isa Integer && value >= 1 || error("`MultiGenerationSize` " * + "must be a positive integer (got `$value`).") + return nothing +end + +# An unset field (`nothing`) reads back as the attribute default. +for (attr, field) in ( + (NumIterationLimit, :num_iteration_limit), + (SetCoverIterationLimit, :set_cover_iteration_limit), + (SetCoverProblem, :set_cover_problem), + (M_Value, :m_value), + (MasterReformulation, :master_reformulation), + (MaxSlack, :max_slack), + (OASlack, :oa_slack), + (UseNLPF, :use_nlpf), + (ConvergenceTolerance, :convergence_tolerance), + (SlackTolerance, :slack_tolerance), + (IterationTimeLimit, :iteration_time_limit), + (MultiGenerationSize, :multi_generation_size), + ) + @eval begin + MOI.supports(::LOA, ::$attr) = true + function MOI.set(algorithm::LOA, attr::$attr, value) + _validate(attr, value) + algorithm.$field = value + return + end + function MOI.get(algorithm::LOA, attr::$attr) + return something(algorithm.$field, _default(algorithm, attr)) + end + end +end + +# `nothing` is these attributes' default, so they skip the +# `something`-based fallback the loop above generates +MOI.supports(::LOA, ::SubproblemMethod) = true +function MOI.set(algorithm::LOA, ::SubproblemMethod, value) + algorithm.subproblem_method = value + return +end +MOI.get(algorithm::LOA, ::SubproblemMethod) = algorithm.subproblem_method + +MOI.supports(::LOA, ::CombinationSource) = true +function MOI.set(algorithm::LOA, ::CombinationSource, value) + algorithm.combination_source = value + return +end +MOI.get(algorithm::LOA, ::CombinationSource) = algorithm.combination_source + +_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 + +# The NLP fixes a feasible indicator combination, so an unbounded NLP +# proves the whole model unbounded; only a genuine infeasibility +# certificate lets a combination count toward exhaustion. +_nlp_unbounded(status::MOI.TerminationStatusCode) = + status in (MOI.DUAL_INFEASIBLE, MOI.NORM_LIMIT) +_nlp_infeasible(status::MOI.TerminationStatusCode) = + status in (MOI.INFEASIBLE, MOI.LOCALLY_INFEASIBLE) + +# 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 + +# The covering problem: a fresh logic-only model, or the master itself +_set_cover_model(::LogicOnlyCover, model, problem, master) = + _build_set_cover(model, problem) +_set_cover_model(::MasterCover, model, problem, master) = master + +# A covering solve on the master leaves the covering objective behind; +# the OA objective goes back before any cut is added +_restore_objective(::LogicOnlyCover, master::_Master) = nothing +_restore_objective(::MasterCover, master::_Master) = + _set_master_objective(master, master.sense, master.oa_objective) + +# The logic model needs its own no-good cut to move on; the master +# receives the combination's no-good cut from `process_result` +_exclude_combination(::LogicOnlyCover, set_cover, combination) = + _avoid_combination(set_cover, combination) +_exclude_combination(::MasterCover, set_cover, combination) = nothing + +# Turkay and Grossmann's covering weights, as GDPopt uses them: an +# uncovered disjunct outweighs all covered ones +function _cover_objective( + set_cover::Union{_Master, _SetCoverModel}, + 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(set_cover.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 _solve_master( + model::Optimizer, + master::Union{_Master, _SetCoverModel}, + deadline::Float64 + ) + _cap_remaining_time(master.model, deadline) + model.master_time += @elapsed MOI.optimize!(master.model) + model.num_master_solves += 1 + return _solved_and_feasible(master.model) +end + +# The master's proposal plus up to `count - 1` more from the +# `CombinationSource`. The second return value says whether the source +# already added the no-good cuts (re-solving sources must), so the +# caller skips them. +function _extract_combinations( + source, + model::Optimizer, + problem::_Problem, + master::_Master, + incumbent, + count::Int, + deadline::Float64 + ) + proposal = _extract_combination(problem, master) + count > 1 || return [proposal], false + extra, excluded = _candidate_combinations(source, model, problem, + master, proposal, incumbent, count - 1, deadline) + return [proposal; extra], excluded +end + +function _set_master_objective( + master::Union{_Master, _SetCoverModel}, + sense, + objective + ) + MOI.set(master.model, MOI.ObjectiveSense(), sense) + MOI.set(master.model, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + objective) + return +end + +################################################################################ +# MAIN LOOP +################################################################################ +function _optimize!(algorithm::LOA, model::Optimizer) + t_start = time() + _reset_results(model) + problem = _build_problem(model) + method = MOI.get(algorithm, SubproblemMethod()) + _check_inner_support(method, model, problem) + master = _build_master(model, problem) + subproblem = _build_subproblem(method, model, problem) + linearizer = _Linearizer() + sense = problem.sense + overall_deadline = t_start + something(model.time_limit_sec, Inf) + loop_deadline = min(overall_deadline, + t_start + Float64(MOI.get(algorithm, IterationTimeLimit()))) + + best_objective = _worst_objective(sense) + best_result = nothing + previous_result = _user_start_values(model, problem) + master_bound = nothing + master_status = nothing + converged = false + unbounded = false + num_unresolved = 0 + + # Shared iteration tail: no-good cut, OA cuts, incumbent update. + # An unbounded NLP aborts the loops; a failed (neither feasible + # nor proven-infeasible) NLP still gets its no-good cut so the + # loop progresses, but is counted against the exhaustion claim. + process_result = (result; nogood = true) -> begin + if !result.feasible && _nlp_unbounded(result.status) + unbounded = true + return + end + result.feasible || _nlp_infeasible(result.status) || + (num_unresolved += 1) + result.feasible || (model.num_nlp_infeasible += 1) + nogood && _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: the covering problem (`SetCoverProblem`) picks + # combinations that activate every nonlinear disjunct once; the + # master receives their NLPs' cuts + cover = _cover_disjuncts(problem) + needs_cover = trues(length(cover)) + num_covered = 0 + covering_algorithm = MOI.get(algorithm, SetCoverProblem()) + set_cover = _set_cover_model(covering_algorithm, model, problem, master) + for iteration in 1:MOI.get(algorithm, SetCoverIterationLimit()) + (iteration == 1 || any(needs_cover)) || break + time() < loop_deadline || break + _set_master_objective(set_cover, MOI.MAX_SENSE, + _cover_objective(set_cover, cover, needs_cover, num_covered)) + solved = _solve_master(model, set_cover, loop_deadline) + combination = solved ? _extract_combination(problem, set_cover) : + nothing + _restore_objective(covering_algorithm, master) + solved || break + # the remaining targets are unreachable once a solve covers none + iteration == 1 || any(needs_cover[i] && + _disjunct_active(combination, cover[i]) + for i in eachindex(cover)) || break + _exclude_combination(covering_algorithm, set_cover, combination) + t_nlp = time() + result = _solve_nlp(method, model, problem, subproblem, + combination, warm_start(); deadline = loop_deadline) + model.nlp_time += time() - t_nlp + model.num_nlp_solves += 1 + process_result(result) + unbounded && break + # 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 && !unbounded + for _ in 1:MOI.get(algorithm, NumIterationLimit()) + time() < loop_deadline || break + if !_solve_master(model, master, loop_deadline) + 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(MOI.get(algorithm, OASlack())) + tol = Float64(MOI.get(algorithm, ConvergenceTolerance())) * + max(abs(best_objective), 1.0) + if gap <= tol && total_slack <= + Float64(MOI.get(algorithm, SlackTolerance())) + converged = true + break + end + end + combinations, excluded = _extract_combinations( + MOI.get(algorithm, CombinationSource()), model, problem, + master, best_result === nothing ? nothing : + best_result.combination, + MOI.get(algorithm, MultiGenerationSize()), loop_deadline) + t_nlp = time() + results = _solve_nlps(method, model, problem, subproblem, + combinations, warm_start(); deadline = loop_deadline) + model.nlp_time += time() - t_nlp + model.num_nlp_solves += length(results) + for result in results + process_result(result; nogood = !excluded) + unbounded && break + end + unbounded && break + end + end + + _store_results(model, sense, best_objective, best_result, master_bound, + master_status, converged, unbounded, num_unresolved, 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, + unbounded::Bool, + num_unresolved::Int, + loop_deadline::Float64 + ) + timed_out = time() >= loop_deadline + if unbounded + model.termination_status = MOI.DUAL_INFEASIBLE + model.primal_status = MOI.NO_SOLUTION + model.objective_value = NaN + model.raw_status = "An NLP subproblem at a fixed indicator " * + "combination is unbounded, so the model is unbounded." + return + end + 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 && num_unresolved == 0 + model.termination_status = MOI.INFEASIBLE + model.raw_status = "No feasible incumbent: the master " * + "problem is infeasible." + elseif master_status == MOI.INFEASIBLE + # exhaustion is not an infeasibility proof while some + # combination never solved to a certificate + model.termination_status = MOI.OTHER_LIMIT + model.raw_status = "No feasible incumbent: the master " * + "is infeasible, but $num_unresolved combination(s) " * + "failed without an infeasibility certificate." + elseif timed_out + model.termination_status = MOI.TIME_LIMIT + elseif master_status === nothing + model.termination_status = MOI.ITERATION_LIMIT + elseif master_status == MOI.DUAL_INFEASIBLE + model.termination_status = MOI.OTHER_LIMIT + model.raw_status = "No feasible incumbent: the master is " * + "unbounded (no OA cut bounds the objective yet); the " * + "model itself may be unbounded." + 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) + exhausted = master_status == MOI.INFEASIBLE && num_unresolved == 0 + if converged + model.termination_status = MOI.LOCALLY_SOLVED + elseif timed_out + model.termination_status = MOI.TIME_LIMIT + elseif exhausted + # 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 == MOI.INFEASIBLE + # visited, but some combination failed without a certificate, + # so the incumbent is not best over all of them + model.termination_status = MOI.OTHER_LIMIT + elseif master_status !== nothing + model.termination_status = MOI.OTHER_LIMIT + else + model.termination_status = MOI.ITERATION_LIMIT + end + label = converged ? "converged" : + exhausted ? "combinations exhausted" : + master_status == MOI.INFEASIBLE ? "combinations exhausted, " * + "$num_unresolved unresolved" : "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/algorithms/LOA/cuts.jl b/src/algorithms/LOA/cuts.jl new file mode 100644 index 0000000..0c79bf8 --- /dev/null +++ b/src/algorithms/LOA/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, 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, 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(MOI.get(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, 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/algorithms/LOA/master.jl b/src/algorithms/LOA/master.jl new file mode 100644 index 0000000..a67462c --- /dev/null +++ b/src/algorithms/LOA/master.jl @@ -0,0 +1,216 @@ +################################################################################ +# 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 + +# Turkay and Grossmann's set-covering problem: the disjunct binaries +# under the exactly-one rows and the propositional rows only, with no +# continuous variables and no cuts +struct _SetCoverModel + model::MOI.ModelLike + variable_map::Dict{MOI.VariableIndex, MOI.VariableIndex} +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 + +# The solver keeps any limit of its own (a per-solve cap from its +# factory) and only ever gets the remaining budget when that is less +function _cap_remaining_time(solver::MOI.ModelLike, deadline::Float64) + isfinite(deadline) || return + MOI.supports(solver, MOI.TimeLimitSec()) || return + remaining = max(0.0, deadline - time()) + limit = MOI.get(solver, MOI.TimeLimitSec()) + MOI.set(solver, MOI.TimeLimitSec(), + limit === nothing ? remaining : min(limit, remaining)) + return +end + +function _build_master(model::Optimizer, problem::_Problem) + mip = _instantiate(model.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) + 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 + +function _build_set_cover(model::Optimizer, problem::_Problem) + mip = _instantiate(model.mip_solver) + binaries = Set(problem.binaries) + variable_map = Dict{MOI.VariableIndex, MOI.VariableIndex}( + vi => MOI.add_variable(mip) for vi in problem.binaries) + for ci in problem.variable_cis + vi = MOI.get(model.cache, MOI.ConstraintFunction(), ci) + vi in binaries || continue + MOI.add_constraint(mip, variable_map[vi], + MOI.get(model.cache, MOI.ConstraintSet(), ci)) + end + # propositional logic arrives as affine rows over the binaries + for ci in problem.linear_cis + func = MOI.get(model.cache, MOI.ConstraintFunction(), ci) + all(term.variable in binaries for term in func.terms) || continue + MOI.add_constraint(mip, _map_to(variable_map, func), + MOI.get(model.cache, MOI.ConstraintSet(), ci)) + end + for disjunction in problem.disjunctions + _add_exactly_one(mip, variable_map, disjunction) + end + return _SetCoverModel(mip, variable_map) +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 +# `MasterReformulation`; nonlinear rows enter the master only as OA +# cuts +function _add_gated_rows( + mip::MOI.ModelLike, + variable_map::AbstractDict, + disjunct::_Disjunct, + model::Optimizer + ) + gating = MOI.get(model, MasterReformulation()) + 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(MOI.get(model, 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, + model::Optimizer, + 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(MOI.get(model, MaxSlack())))) + penalty = penalty_sign * Float64(MOI.get(model, OASlack())) + 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::Union{_Master, _SetCoverModel}, + index = 1 + ) + return Dict{MOI.VariableIndex, Bool}( + binary => round(Bool, MOI.get(master.model, + MOI.VariablePrimal(index), 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::Union{_Master, _SetCoverModel}, + 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/algorithms/LOA/nlp.jl b/src/algorithms/LOA/nlp.jl new file mode 100644 index 0000000..0916505 --- /dev/null +++ b/src/algorithms/LOA/nlp.jl @@ -0,0 +1,369 @@ +################################################################################ +# NLP SUBPROBLEM +################################################################################ +# The MOI subproblem path dispatches on the `SubproblemMethod` value +# `nothing`; a method object routes `_build_subproblem`/`_solve_nlp` +# to its own construction (see the `SubproblemMethod` docstring). +# 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 + +# One copy of the variables, bounds, binary fixes (at 0 until a +# combination is set), and global rows; disjunct rows depend on the +# combination and are the caller's. +function _add_copy(nlp::MOI.ModelLike, model::Optimizer, problem::_Problem) + 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 + return variable_map, fixes +end + +# Point one copy at a combination: overwrite its binary fixes and add +# the active disjuncts' rows, returned so a reused copy can drop them. +function _fix_combination( + nlp::MOI.ModelLike, + variable_map::AbstractDict, + fixes::AbstractDict, + problem::_Problem, + combination::AbstractDict + ) + for (binary, value) in combination + MOI.set(nlp, MOI.ConstraintSet(), fixes[binary], + MOI.EqualTo(value ? 1.0 : 0.0)) + end + rows = MOI.ConstraintIndex[] + for disjunction in problem.disjunctions, disjunct in disjunction.disjuncts + _disjunct_active(combination, disjunct) || continue + for (func, set) in zip(disjunct.functions, disjunct.sets) + push!(rows, MOI.add_constraint(nlp, + _map_to(variable_map, func), set)) + end + end + return rows +end + +function _build_subproblem(::Nothing, model::Optimizer, problem::_Problem) + nlp = _instantiate(model.nlp_solver) + variable_map, fixes = _add_copy(nlp, model, problem) + 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( + ::Nothing, + model::Optimizer, + problem::_Problem, + sub::_Subproblem, + combination::AbstractDict, + warm_start; + deadline::Float64 = Inf + ) + for ci in sub.rows + MOI.delete(sub.model, ci) + end + append!(empty!(sub.rows), _fix_combination(sub.model, sub.variable_map, + sub.fixes, problem, combination)) + _set_warm_start(sub.model, sub.variable_map, warm_start) + _cap_remaining_time(sub.model, deadline) + MOI.optimize!(sub.model) + status = MOI.get(sub.model, MOI.TerminationStatus()) + 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, status = status) + end + if Bool(MOI.get(model, UseNLPF())) + result = _solve_nlpf(model, problem, combination, warm_start; + deadline = deadline) + result === nothing || return (; result..., status = status) + end + return (combination = combination, + point = nothing, objective = Inf, feasible = false, status = status) +end + +# Several combinations at once; the default method has no batch form +# and takes them one at a time. +function _solve_nlps( + method, + model::Optimizer, + problem::_Problem, + sub, + combinations::AbstractVector, + warm_start; + deadline::Float64 = Inf + ) + return [_solve_nlp(method, model, problem, sub, combination, warm_start; + deadline) for combination in combinations] +end + +################################################################################ +# BATCHED SUBPROBLEMS +################################################################################ +""" + BatchedSubproblems() + +`SubproblemMethod` value that solves a vector of indicator combinations +as one block-diagonal NLP through the `nlp_solver`: each combination +gets its own copy of the variables and rows, and the objective is the +sum of the copies' objectives. The copies do not interact, so the +stacked solution is the sequential one; what changes is that the +solver sees one large, regular problem, which is what GPU evaluators +and factorizations need to pay off. A single combination is a batch of +one. A batch whose joint solve fails (one infeasible copy makes the +whole stack infeasible) falls back to the sequential path copy by copy, +so the results match the default method exactly. +""" +struct BatchedSubproblems end + +# the sequential subproblem is the fallback for failed batches +struct _BatchedSubproblem + sequential::_Subproblem +end + +_check_nlp_support(::BatchedSubproblems, model::Optimizer, problem::_Problem) = + _check_nlp_support(nothing, model, problem) + +function _build_subproblem( + ::BatchedSubproblems, + model::Optimizer, + problem::_Problem + ) + return _BatchedSubproblem(_build_subproblem(nothing, model, problem)) +end + +function _solve_nlp( + method::BatchedSubproblems, + model::Optimizer, + problem::_Problem, + sub::_BatchedSubproblem, + combination::AbstractDict, + warm_start; + deadline::Float64 = Inf + ) + return only(_solve_nlps(method, model, problem, sub, [combination], + warm_start; deadline)) +end + +# One stacked solve of the combinations as given; `nothing` when the +# joint problem did not solve to a feasible point. +function _solve_stack( + model::Optimizer, + problem::_Problem, + combinations::AbstractVector, + warm_start, + deadline::Float64 + ) + nlp = _instantiate(model.nlp_solver) + copies = map(combinations) do combination + variable_map, fixes = _add_copy(nlp, model, problem) + _fix_combination(nlp, variable_map, fixes, problem, combination) + _set_warm_start(nlp, variable_map, warm_start) + return variable_map + end + objective = foldl((a, b) -> _combine(+, a, b), + (_map_to(variable_map, problem.objective) for variable_map in copies)) + MOI.set(nlp, MOI.ObjectiveSense(), problem.sense) + MOI.set(nlp, MOI.ObjectiveFunction{typeof(objective)}(), objective) + _cap_remaining_time(nlp, deadline) + MOI.optimize!(nlp) + _solved_and_feasible(nlp) || return nothing + status = MOI.get(nlp, MOI.TerminationStatus()) + return map(zip(combinations, copies)) do (combination, variable_map) + point = _extract_point(nlp, problem, variable_map) + objective = MOI.Utilities.eval_variables(vi -> point[vi], + model.cache, problem.objective) + return (combination = combination, point = point, + objective = objective, feasible = true, status = status) + end +end + +# Stacked feasibility pass: every copy slacked, sum of slacks minimized. +# Copies with zero slack are feasible; the others' points are their +# restoration points. `nothing` when even the slacked stack fails. +function _solve_slacked_stack( + model::Optimizer, + problem::_Problem, + combinations::AbstractVector, + warm_start, + deadline::Float64 + ) + nlp = _instantiate(model.nlp_solver) + copies = map(combinations) do combination + variable_map, u = _add_slacked_copy(nlp, model, problem, combination) + _set_warm_start(nlp, variable_map, warm_start) + return variable_map, u + end + MOI.set(nlp, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.set(nlp, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + MOI.ScalarAffineFunction( + [MOI.ScalarAffineTerm(1.0, u) for (_, u) in copies], 0.0)) + _cap_remaining_time(nlp, deadline) + MOI.optimize!(nlp) + _solved_and_feasible(nlp) || return nothing + tolerance = Float64(MOI.get(model, SlackTolerance())) + feasible = [MOI.get(nlp, MOI.VariablePrimal(), u) <= tolerance + for (_, u) in copies] + points = [_extract_point(nlp, problem, variable_map) + for (variable_map, _) in copies] + return feasible, points +end + +# Stack first; if one copy is infeasible the stack is, so classify the +# copies with a stacked feasibility pass and re-stack the feasible +# ones. Sequential solves remain only for a stack that fails on its +# own account. +function _solve_nlps( + ::BatchedSubproblems, + model::Optimizer, + problem::_Problem, + sub::_BatchedSubproblem, + combinations::AbstractVector, + warm_start; + deadline::Float64 = Inf + ) + sequential = () -> _solve_nlps(nothing, model, problem, sub.sequential, + combinations, warm_start; deadline) + results = _solve_stack(model, problem, combinations, warm_start, deadline) + results === nothing || return results + Bool(MOI.get(model, UseNLPF())) || return sequential() + restored = _solve_slacked_stack(model, problem, combinations, warm_start, + deadline) + restored === nothing && return sequential() + feasible, points = restored + results = Vector{Any}(undef, length(combinations)) + active = findall(feasible) + if !isempty(active) + solved = _solve_stack(model, problem, combinations[active], + warm_start, deadline) + solved === nothing && (solved = _solve_nlps(nothing, model, problem, + sub.sequential, combinations[active], warm_start; deadline)) + results[active] .= solved + end + for i in findall(!, feasible) + results[i] = (combination = combinations[i], point = points[i], + objective = Inf, feasible = false, + status = MOI.LOCALLY_INFEASIBLE) + end + return results +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. +# One slacked copy at a combination: bounds and binary fixes exact, +# every inequality row relaxed by the copy's nonnegative `u`. +function _add_slacked_copy( + nlp::MOI.ModelLike, + model::Optimizer, + problem::_Problem, + combination::AbstractDict + ) + 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) + MOI.add_constraint(nlp, something(slacked, mapped), set) + end + return variable_map, u +end + +function _solve_nlpf( + model::Optimizer, + problem::_Problem, + combination::AbstractDict, + warm_start; + deadline::Float64 = Inf + ) + nlp = _instantiate(model.nlp_solver) + variable_map, u = _add_slacked_copy(nlp, model, problem, combination) + 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 hurts 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/combination_sources.jl b/src/combination_sources.jl new file mode 100644 index 0000000..fb12ef2 --- /dev/null +++ b/src/combination_sources.jl @@ -0,0 +1,330 @@ +################################################################################ +# COMBINATION SOURCES +################################################################################ +# Where the extra combinations of a multi-generation iteration come +# from. `_candidate_combinations(source, model, problem, master, +# proposal, incumbent, count, deadline)` returns up to `count` +# combinations other than `proposal` and whether it already added their +# no-good cuts to the master (a re-solving source must, so the loop +# skips them). The default `nothing` source is the solver's pool, then +# re-solves; the rest are experimental and need nothing from the solver. + +const _Combination = Dict{MOI.VariableIndex, Bool} + +# one chosen disjunct per disjunction -> combination +function _combination(problem::_Problem, chosen) + combination = _Combination() + for (disjunction, active) in zip(problem.disjunctions, chosen), + disjunct in disjunction.disjuncts + combination[disjunct.binary] = disjunct === active ? + disjunct.active_value : !disjunct.active_value + end + return combination +end + +# the active disjunct of each disjunction under a combination +function _active_disjuncts(problem::_Problem, combination::AbstractDict) + return map(problem.disjunctions) do disjunction + index = findfirst(d -> _disjunct_active(combination, d), + disjunction.disjuncts) + return disjunction.disjuncts[something(index, 1)] + end +end + +_push_candidate!(candidates, combination, proposal) = + combination == proposal || combination in candidates || + push!(candidates, combination) + +################################################################################ +# DEFAULT: POOL, OR RE-SOLVES +################################################################################ +# Result indices past 1 are the solver's pool (Gurobi's PoolSolutions); +# a short pool means a short batch. Only a solver without a pool falls +# back to re-solving behind a no-good cut per combination taken. +function _candidate_combinations( + ::Nothing, + model::Optimizer, + problem::_Problem, + master::_Master, + proposal::AbstractDict, + incumbent, + count::Int, + deadline::Float64 + ) + candidates = _Combination[] + for index in 2:min(count + 1, MOI.get(master.model, MOI.ResultCount())) + _push_candidate!(candidates, + _extract_combination(problem, master, index), proposal) + end + MOI.get(master.model, MOI.ResultCount()) > 1 && return candidates, false + return _resolve_candidates(model, problem, master, proposal, candidates, + count, deadline) +end + +# Exclude the proposal and every candidate so far, then re-solve until +# `count` are in hand; the cuts stay, so the caller must not add them. +function _resolve_candidates( + model::Optimizer, + problem::_Problem, + master::_Master, + proposal::AbstractDict, + candidates::Vector{_Combination}, + count::Int, + deadline::Float64 + ) + _avoid_combination(master, proposal) + foreach(c -> _avoid_combination(master, c), candidates) + while length(candidates) < count && time() < deadline + _solve_master(model, master, deadline) || break + combination = _extract_combination(problem, master) + combination == proposal && break + combination in candidates && break + push!(candidates, combination) + _avoid_combination(master, combination) + end + return candidates, true +end + +################################################################################ +# CUTOFF RE-SOLVES +################################################################################ +""" + CutoffResolve(; tolerance = 0.05) + +[`CombinationSource`](@ref) that re-solves the master behind no-good +cuts like the default, but with the master objective capped at the +current optimum plus `tolerance` (relative) and the last solution as a +warm start, so each re-solve is a small warm MILP. Still one master +solve per extra combination. +""" +struct CutoffResolve + tolerance::Float64 + function CutoffResolve(; tolerance::Real = 0.05) + tolerance >= 0 || error("`CutoffResolve` tolerance must be " * + "nonnegative (got `$tolerance`).") + return new(Float64(tolerance)) + end +end + +function _candidate_combinations( + source::CutoffResolve, + model::Optimizer, + problem::_Problem, + master::_Master, + proposal::AbstractDict, + incumbent, + count::Int, + deadline::Float64 + ) + # read the solution before any modification invalidates it + value = MOI.get(master.model, MOI.ObjectiveValue()) + starts = _primal_starts(master) + slack = source.tolerance * max(1.0, abs(value)) + set = master.sense == MOI.MAX_SENSE ? MOI.GreaterThan(value - slack) : + MOI.LessThan(value + slack) + cutoff = MOI.add_constraint(master.model, master.oa_objective, set) + _set_primal_starts(master, starts) + candidates, excluded = _resolve_candidates(model, problem, master, + proposal, _Combination[], count, deadline) + MOI.delete(master.model, cutoff) + return candidates, excluded +end + +function _primal_starts(master::_Master) + MOI.supports(master.model, MOI.VariablePrimalStart(), + MOI.VariableIndex) || return nothing + return [vi => MOI.get(master.model, MOI.VariablePrimal(), vi) + for vi in MOI.get(master.model, MOI.ListOfVariableIndices())] +end + +function _set_primal_starts(master::_Master, starts) + starts === nothing && return + for (vi, value) in starts + MOI.set(master.model, MOI.VariablePrimalStart(), vi, value) + end + return +end + +################################################################################ +# NEIGHBORHOOD +################################################################################ +""" + Neighborhood(; around = :proposal) + +[`CombinationSource`](@ref) made of every combination that differs from +the master's proposal (or from the incumbent, `around = :incumbent`) in +exactly one disjunction's choice. Free: no solver call. A disjunction +of `m` disjuncts contributes `m - 1` neighbors; when there are more +than asked for, an evenly spaced subset over the disjunctions is taken. +""" +struct Neighborhood + around::Symbol + function Neighborhood(; around::Symbol = :proposal) + around in (:proposal, :incumbent) || error("`Neighborhood` " * + "`around` must be `:proposal` or `:incumbent` (got `$around`).") + return new(around) + end +end + +function _candidate_combinations( + source::Neighborhood, + model::Optimizer, + problem::_Problem, + master::_Master, + proposal::AbstractDict, + incumbent, + count::Int, + deadline::Float64 + ) + center = source.around == :incumbent && incumbent !== nothing ? + incumbent : proposal + chosen = _active_disjuncts(problem, center) + candidates = _Combination[] + for (i, disjunction) in enumerate(problem.disjunctions), + disjunct in disjunction.disjuncts + disjunct === chosen[i] && continue + flipped = copy(chosen) + flipped[i] = disjunct + _push_candidate!(candidates, _combination(problem, flipped), proposal) + end + length(candidates) <= count && return candidates, false + count == 1 && return candidates[1:1], false + picks = round.(Int, range(1, length(candidates); length = count)) + return candidates[unique(picks)], false +end + +################################################################################ +# LP RELAXATION ROUNDING +################################################################################ +""" + LPRounding(; seed = 0) + +[`CombinationSource`](@ref) that samples combinations from the master's +LP relaxation: integrality is dropped on a copy of the master, the LP +is solved, and each disjunction's active disjunct is drawn with +probability proportional to the relaxed value of its activation. One +LP per iteration, informed by every cut in the master. The LP solve is +not counted in `MasterSolveCount`. +""" +struct LPRounding + rng::Random.MersenneTwister + LPRounding(; seed::Integer = 0) = new(Random.MersenneTwister(seed)) +end + +function _candidate_combinations( + source::LPRounding, + model::Optimizer, + problem::_Problem, + master::_Master, + proposal::AbstractDict, + incumbent, + count::Int, + deadline::Float64 + ) + weights = _relaxed_activations(model, problem, master, deadline) + weights === nothing && return _Combination[], false + candidates = _Combination[] + for _ in 1:(20 * count) + chosen = [_sample(source.rng, disjunction.disjuncts, w) for + (disjunction, w) in zip(problem.disjunctions, weights)] + _push_candidate!(candidates, _combination(problem, chosen), proposal) + length(candidates) == count && break + end + return candidates, false +end + +# relaxed activation value of every disjunct, per disjunction +function _relaxed_activations( + model::Optimizer, + problem::_Problem, + master::_Master, + deadline::Float64 + ) + lp = _instantiate(model.mip_solver) + index_map = MOI.copy_to(lp, master.model) + # indicator-gated rows have no relaxation without their binary, so + # they leave with the integrality; the big-M cuts stay + for (F, S) in MOI.get(lp, MOI.ListOfConstraintTypesPresent()) + if S <: MOI.Indicator + foreach(ci -> MOI.delete(lp, ci), + collect(MOI.get(lp, MOI.ListOfConstraintIndices{F, S}()))) + continue + end + F === MOI.VariableIndex && S in (MOI.ZeroOne, MOI.Integer) || continue + for ci in collect(MOI.get(lp, MOI.ListOfConstraintIndices{F, S}())) + vi = MOI.get(lp, MOI.ConstraintFunction(), ci) + MOI.delete(lp, ci) + S === MOI.ZeroOne && _bound_unit_interval(lp, vi) + end + end + _cap_remaining_time(lp, deadline) + MOI.optimize!(lp) + _solved_and_feasible(lp) || return nothing + value = vi -> MOI.get(lp, MOI.VariablePrimal(), + index_map[master.variable_map[vi]]) + return [[clamp(MOI.Utilities.eval_variables(value, disjunct.activation), + 0.0, 1.0) for disjunct in disjunction.disjuncts] + for disjunction in problem.disjunctions] +end + +# [0, 1] bounds for a relaxed binary, unless bounds already exist +function _bound_unit_interval(lp::MOI.ModelLike, vi::MOI.VariableIndex) + interval = MOI.ConstraintIndex{MOI.VariableIndex, + MOI.Interval{Float64}}(vi.value) + MOI.is_valid(lp, interval) && return + lower = MOI.ConstraintIndex{MOI.VariableIndex, + MOI.GreaterThan{Float64}}(vi.value) + MOI.is_valid(lp, lower) || MOI.add_constraint(lp, vi, MOI.GreaterThan(0.0)) + upper = MOI.ConstraintIndex{MOI.VariableIndex, + MOI.LessThan{Float64}}(vi.value) + MOI.is_valid(lp, upper) || MOI.add_constraint(lp, vi, MOI.LessThan(1.0)) + return +end + +# weighted draw; uniform when the weights carry no information +function _sample(rng, disjuncts, weights) + total = sum(weights) + total > 0 || return rand(rng, disjuncts) + threshold = rand(rng) * total + running = 0.0 + for (disjunct, w) in zip(disjuncts, weights) + running += w + running >= threshold && return disjunct + end + return last(disjuncts) +end + +################################################################################ +# RANDOM COMBINATIONS +################################################################################ +""" + RandomCombinations(; seed = 0) + +[`CombinationSource`](@ref) drawing combinations uniformly, one active +disjunct per disjunction. The uninformed baseline for the sources +above. +""" +struct RandomCombinations + rng::Random.MersenneTwister + RandomCombinations(; seed::Integer = 0) = new(Random.MersenneTwister(seed)) +end + +function _candidate_combinations( + source::RandomCombinations, + model::Optimizer, + problem::_Problem, + master::_Master, + proposal::AbstractDict, + incumbent, + count::Int, + deadline::Float64 + ) + candidates = _Combination[] + for _ in 1:(20 * count) + chosen = [rand(source.rng, disjunction.disjuncts) + for disjunction in problem.disjunctions] + _push_candidate!(candidates, _combination(problem, chosen), proposal) + length(candidates) == count && break + end + return candidates, false +end diff --git a/src/optimizer.jl b/src/optimizer.jl new file mode 100644 index 0000000..e460a2c --- /dev/null +++ b/src/optimizer.jl @@ -0,0 +1,642 @@ +################################################################################ +# ALGORITHMS AND ATTRIBUTES +################################################################################ +""" + AbstractAlgorithm + +A super-type for the solution algorithms. Select one with the +[`Algorithm`](@ref) attribute. +""" +abstract type AbstractAlgorithm end + +""" + Algorithm() <: MOI.AbstractOptimizerAttribute + +The algorithm the optimizer runs. Defaults to [`LOA`](@ref). + +**Example** +```julia +julia> set_attribute(model, DisjunctiveAlgorithms.Algorithm(), + DisjunctiveAlgorithms.LOA()) +``` +""" +struct Algorithm <: MOI.AbstractOptimizerAttribute end + +""" + AbstractAlgorithmAttribute <: MOI.AbstractOptimizerAttribute + +A super-type for algorithm-specific optimizer attributes. Each +algorithm declares `MOI.supports` for the attributes it consumes and +documents them in its docstring. +""" +abstract type AbstractAlgorithmAttribute <: MOI.AbstractOptimizerAttribute end + +_default(::AbstractAlgorithm, attr::AbstractAlgorithmAttribute) = + _default(attr) + +""" + NumIterationLimit() <: AbstractAlgorithmAttribute -> Int + +Master/NLP iterations after the set-covering seed. Defaults to `10`. +""" +struct NumIterationLimit <: AbstractAlgorithmAttribute end +_default(::NumIterationLimit) = 10 + +""" + SetCoverIterationLimit() <: AbstractAlgorithmAttribute -> Int + +Set-covering initialization iterations: each solves the logic-only +cover MIP over the indicators and one NLP, until every nonlinear +disjunct has been active in a feasible NLP. Defaults to `8`. +""" +struct SetCoverIterationLimit <: AbstractAlgorithmAttribute end +_default(::SetCoverIterationLimit) = 8 + +""" + AbstractSetCover + +A super-type for the problems the set-covering pass can solve for its +initial combinations; see [`SetCoverProblem`](@ref). +""" +abstract type AbstractSetCover end + +""" + LogicOnlyCover() <: AbstractSetCover + +The covering problem of Turkay and Grossmann (1996, Appendix D): the +disjunct binaries under the exactly-one and propositional rows only. +""" +struct LogicOnlyCover <: AbstractSetCover end + +""" + MasterCover() <: AbstractSetCover + +The MILP master with the covering objective swapped in, as Pyomo's +GDPopt does (Chen, Johnson, Siirola and Grossmann 2018), so the +initial combinations also satisfy the linear rows and their NLPs are +far less often infeasible. +""" +struct MasterCover <: AbstractSetCover end + +""" + SetCoverProblem() <: AbstractAlgorithmAttribute -> AbstractSetCover + +Which problem the set-covering pass solves for its initial combinations, +[`LogicOnlyCover`](@ref) or [`MasterCover`](@ref). Defaults to +`LogicOnlyCover()`. +""" +struct SetCoverProblem <: AbstractAlgorithmAttribute end +_default(::SetCoverProblem) = LogicOnlyCover() + +""" + M_Value() <: AbstractAlgorithmAttribute -> Float64 + +Big-M gating the disjunct OA cuts in the master. Defaults to `1e9`. +""" +struct M_Value <: AbstractAlgorithmAttribute end +_default(::M_Value) = 1e9 + +""" + MasterReformulation() <: AbstractAlgorithmAttribute -> String + +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). Defaults to `"indicator"`. +""" +struct MasterReformulation <: AbstractAlgorithmAttribute end +_default(::MasterReformulation) = "indicator" + +""" + MaxSlack() <: AbstractAlgorithmAttribute -> Float64 + +Upper bound of each OA cut slack. Defaults to `1e3`. +""" +struct MaxSlack <: AbstractAlgorithmAttribute end +_default(::MaxSlack) = 1e3 + +""" + OASlack() <: AbstractAlgorithmAttribute -> Float64 + +Objective penalty per unit of cut slack. Defaults to `1e3`. +""" +struct OASlack <: AbstractAlgorithmAttribute end +_default(::OASlack) = 1e3 + +""" + UseNLPF() <: AbstractAlgorithmAttribute -> Bool + +Solve a slacked feasibility NLP when the primary NLP is infeasible, +so its point still seeds OA cuts. Defaults to `true`. +""" +struct UseNLPF <: AbstractAlgorithmAttribute end +_default(::UseNLPF) = true + +""" + ConvergenceTolerance() <: AbstractAlgorithmAttribute -> Float64 + +Relative incumbent/bound gap tolerance. Defaults to `1e-6`. +""" +struct ConvergenceTolerance <: AbstractAlgorithmAttribute end +_default(::ConvergenceTolerance) = 1e-6 + +""" + SlackTolerance() <: AbstractAlgorithmAttribute -> Float64 + +Total cut slack tolerance for convergence. Defaults to `1e-4`. +""" +struct SlackTolerance <: AbstractAlgorithmAttribute end +_default(::SlackTolerance) = 1e-4 + +""" + IterationTimeLimit() <: AbstractAlgorithmAttribute -> Float64 + +Seconds allotted to the solve loop. Defaults to `Inf`. The overall +budget is the standard `MOI.TimeLimitSec` attribute. +""" +struct IterationTimeLimit <: AbstractAlgorithmAttribute end +_default(::IterationTimeLimit) = Inf + +""" + MultiGenerationSize() <: AbstractAlgorithmAttribute -> Int + +Indicator combinations evaluated per master solve in the main loop +(multi-generation cuts). Defaults to `1`, plain LOA. Extra +combinations come from the master's solution pool when the +`mip_solver` exposes one through `MOI.ResultCount` (Gurobi with +`PoolSolutions`), as many as it holds up to this size, and otherwise +from re-solving the master behind a no-good cut per combination +already taken. All of them go through +`_solve_nlps`, one stacked NLP under [`BatchedSubproblems`](@ref), and +every result adds its cuts before the next master solve. +""" +struct MultiGenerationSize <: AbstractAlgorithmAttribute end +_default(::MultiGenerationSize) = 1 + +""" + CombinationSource() <: AbstractAlgorithmAttribute -> Any + +Where the extra combinations of a [`MultiGenerationSize`](@ref) +iteration come from. The default `nothing` reads the master's solution +pool (result indices past 1), and only without a pool re-solves the +master behind no-good cuts. The experimental sources in +`combination_sources.jl` need nothing from the solver: +[`Neighborhood`](@ref), [`LPRounding`](@ref), +[`RandomCombinations`](@ref), and [`CutoffResolve`](@ref). A source +extends `_candidate_combinations`. +""" +struct CombinationSource <: AbstractAlgorithmAttribute end +_default(::CombinationSource) = nothing + +""" + SubproblemMethod() <: AbstractAlgorithmAttribute -> Any + +How the NLP subproblems are built and solved. The default `nothing` +builds one MOI model from the `nlp_solver` and rewrites its binary +fixes and disjunct rows per combination. A method object routes the +subproblems elsewhere (e.g. a structure-preserving transcription of +the original model): it must extend `_build_subproblem` and +`_solve_nlp` on its type, and may extend `_check_nlp_support`. +`_solve_nlp` must return `(; combination, point, objective, feasible, +status)` with `point` indexed by this optimizer's variable indices. +""" +struct SubproblemMethod <: AbstractAlgorithmAttribute end +_default(::SubproblemMethod) = nothing + +################################################################################ +# OPTIMIZER +################################################################################ +const _Cache = MOI.Utilities.UniversalFallback{MOI.Utilities.Model{Float64}} + +""" + Optimizer(nlp_solver, mip_solver = nlp_solver) + +Solver for models containing +`DisjunctiveProgramming.DisjunctionSet` constraints. `nlp_solver` +and `mip_solver` are optimizer factories as accepted by +`MOI.instantiate`. The algorithm (default [`LOA`](@ref)) is selected +with the [`Algorithm`](@ref) attribute and configured through the +attributes it documents. A new optimizer starts with a 3600 second +`MOI.TimeLimitSec` safety limit; set it to `nothing` for no limit. + +**Example** +```julia +julia> model = GDPModel(() -> DisjunctiveAlgorithms.Optimizer( + Ipopt.Optimizer, HiGHS.Optimizer)); + +julia> optimize!(model, gdp_method = Direct()) +``` +""" +mutable struct Optimizer <: MOI.AbstractOptimizer + nlp_solver::Any + mip_solver::Any + cache::_Cache + algorithm::Union{Nothing, AbstractAlgorithm} + time_limit_sec::Union{Nothing, Float64} + 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 + num_master_solves::Int + num_nlp_solves::Int + num_nlp_infeasible::Int + master_time::Float64 + nlp_time::Float64 +end + +function Optimizer(nlp_solver, mip_solver = nlp_solver) + return Optimizer(nlp_solver, mip_solver, + MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), + nothing, 3600.0, false, MOI.OPTIMIZE_NOT_CALLED, MOI.NO_SOLUTION, + Dict{MOI.VariableIndex, Float64}(), NaN, nothing, NaN, "", NaN, 0, 0, + 0, 0.0, 0.0) +end + +_algorithm(model::Optimizer) = + something(model.algorithm, _default(Algorithm())) + +# Realize an inner solver factory, bridged and silenced. +function _instantiate(factory) + solver = MOI.instantiate(factory; + with_cache_type = Float64, with_bridge_type = Float64) + MOI.supports(solver, MOI.Silent()) && + MOI.set(solver, MOI.Silent(), true) + return solver +end + +MOI.get(::Optimizer, ::MOI.SolverName) = "DisjunctiveAlgorithms" +MOI.get(::Optimizer, ::MOI.SolverVersion) = string(pkgversion(@__MODULE__)) + +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 + model.num_master_solves = 0 + model.num_nlp_solves = 0 + model.num_nlp_infeasible = 0 + model.master_time = 0.0 + model.nlp_time = 0.0 + 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 + +MOI.optimize!(model::Optimizer) = _optimize!(_algorithm(model), model) + +################################################################################ +# 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) + +# Objective functions `_build_problem` can consume; the generic +# forwarding above would claim vector objectives the solve rejects. +const _ObjectiveFunction = Union{MOI.VariableIndex, + MOI.ScalarAffineFunction{Float64}, + MOI.ScalarQuadraticFunction{Float64}, MOI.ScalarNonlinearFunction} + +function MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{F}) where {F} + return F <: _ObjectiveFunction +end + +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 +# 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.time_limit_sec = value === nothing ? nothing : Float64(value) + return +end + +MOI.get(model::Optimizer, ::MOI.TimeLimitSec) = model.time_limit_sec + +MOI.supports(::Optimizer, ::Algorithm) = true + +MOI.get(model::Optimizer, ::Algorithm) = model.algorithm + +function MOI.set(model::Optimizer, ::Algorithm, algorithm::AbstractAlgorithm) + model.algorithm = algorithm + return +end + +function MOI.supports(model::Optimizer, attr::AbstractAlgorithmAttribute) + return MOI.supports(_algorithm(model), attr) +end + +function MOI.set(model::Optimizer, attr::AbstractAlgorithmAttribute, value) + model.algorithm === nothing && (model.algorithm = _default(Algorithm())) + MOI.set(model.algorithm, attr, value) + return +end + +function MOI.get(model::Optimizer, attr::AbstractAlgorithmAttribute) + return MOI.get(_algorithm(model), attr) +end + +################################################################################ +# RESULT ATTRIBUTES +################################################################################ +MOI.get(model::Optimizer, ::MOI.TerminationStatus) = model.termination_status + +MOI.get(model::Optimizer, ::MOI.RawStatusString) = model.raw_status + +""" + MasterSolveCount() <: MOI.AbstractModelAttribute -> Int + +Master MILP solves in the last `optimize!`, including the set-covering +pass and any pool-replacement re-solves. +""" +struct MasterSolveCount <: MOI.AbstractModelAttribute end + +""" + NLPSolveCount() <: MOI.AbstractModelAttribute -> Int + +Indicator combinations evaluated by NLP subproblems in the last +`optimize!`; a stacked batch of `k` counts `k`. +""" +struct NLPSolveCount <: MOI.AbstractModelAttribute end + +""" + NLPInfeasibleCount() <: MOI.AbstractModelAttribute -> Int + +Indicator combinations whose NLP subproblem returned no feasible point +in the last `optimize!` (infeasible or failed), out of +[`NLPSolveCount`](@ref). +""" +struct NLPInfeasibleCount <: MOI.AbstractModelAttribute end + +""" + MasterSolveTime() <: MOI.AbstractModelAttribute -> Float64 + +Seconds spent inside master MILP solves in the last `optimize!`. +""" +struct MasterSolveTime <: MOI.AbstractModelAttribute end + +""" + NLPSolveTime() <: MOI.AbstractModelAttribute -> Float64 + +Seconds spent inside NLP subproblem solves in the last `optimize!`, +feasibility restoration and stacked batches included. +""" +struct NLPSolveTime <: MOI.AbstractModelAttribute end + +const _SolveStatistic = Union{MasterSolveCount, NLPSolveCount, + NLPInfeasibleCount, MasterSolveTime, NLPSolveTime} +MOI.is_set_by_optimize(::_SolveStatistic) = true +MOI.get(model::Optimizer, ::MasterSolveCount) = model.num_master_solves +MOI.get(model::Optimizer, ::NLPSolveCount) = model.num_nlp_solves +MOI.get(model::Optimizer, ::NLPInfeasibleCount) = model.num_nlp_infeasible +MOI.get(model::Optimizer, ::MasterSolveTime) = model.master_time +MOI.get(model::Optimizer, ::NLPSolveTime) = model.nlp_time + +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..a745ae5 --- /dev/null +++ b/src/problem.jl @@ -0,0 +1,325 @@ +################################################################################ +# 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 +function _to_affine(constant::Real) + return MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], + Float64(constant)) +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 + +# Polynomial degree of an operator tree built from +, -, * and ^; +# `nothing` for any other operator. Promoted affine and quadratic rows +# only use these, so the degree says which MOI type they convert to. +_degree(::Real) = 0 +_degree(::MOI.VariableIndex) = 1 +_degree(::MOI.ScalarAffineFunction) = 1 +_degree(::MOI.ScalarQuadraticFunction) = 2 +function _degree(func::MOI.ScalarNonlinearFunction) + degrees = [_degree(arg) for arg in func.args] + any(isnothing, degrees) && return nothing + func.head in (:+, :-) && return maximum(degrees) + func.head == :* && return sum(degrees) + func.head == :^ && func.args[2] isa Real && + return degrees[1] * func.args[2] + return nothing +end + +# Evaluate a degree <= 2 operator tree into the matching typed +# function. MOI's own `convert` is documented as rough-and-ready: it +# reads a two-argument `*` as coefficient times variable and rejects +# `-`, so promoted rows such as `x * y` or `x * y - 3` throw there. +_polynomial(func::Real) = Float64(func) +_polynomial(func::MOI.VariableIndex) = _to_affine(func) +_polynomial(func::MOI.AbstractScalarFunction) = func +function _polynomial(func::MOI.ScalarNonlinearFunction) + args = [_polynomial(arg) for arg in func.args] + if func.head == :^ + base, exponent = args + base isa Real && return base^exponent + exponent == 0 && return 1.0 + exponent == 1 && return base + return MOI.Utilities.operate(*, Float64, base, base) + end + op = func.head == :+ ? (+) : func.head == :- ? (-) : (*) + length(args) == 1 && return _combine(op, args[1]) + return foldl((a, b) -> _combine(op, a, b), args) +end + +_combine(op, a::Real) = op(a) +_combine(op, a) = MOI.Utilities.operate(op, Float64, a) +_combine(op, a::Real, b::Real) = op(a, b) +_combine(op, a, b) = MOI.Utilities.operate(op, Float64, a, b) + +# demote rows the enclosing vector function promoted +_demote(func::MOI.AbstractScalarFunction) = func +function _demote(func::MOI.ScalarQuadraticFunction{Float64}) + return isempty(func.quadratic_terms) ? _to_affine(func) : func +end +function _demote(func::MOI.ScalarNonlinearFunction) + degree = _degree(func) + degree in (0, 1) && return _to_affine(_polynomial(func)) + degree == 2 && return _demote(_polynomial(func)) + return func +end + +_scalarize(func::MOI.AbstractVectorFunction) = + collect(MOI.Utilities.eachscalar(func)) + +function _indicator_error(func::MOI.AbstractScalarFunction) + return error("Unsupported indicator expression `$func`: each " * + "`DisjunctionSet` indicator must be a binary variable `z` or " * + "its complement `1 - z`.") +end + +# An activation/indicator row must demote to an affine expression; +# genuinely nonlinear rows get the curated error, not a raw convert +# failure. +function _activation_affine(row::MOI.AbstractScalarFunction) + demoted = _demote(row) + _is_linear(demoted) || _indicator_error(demoted) + return _to_affine(demoted) +end + +# `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 _indicator_error(activation) +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 = _activation_affine( + rows[activation_index(set)]) + disjuncts = _Disjunct[] + for (i, j) in enumerate(indicator_indices(set)) + activation = _activation_affine(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`.") + # transcribed rows can carry function constants: shift them + # into the scalar sets so consumers add constant-free rows + normalized = [MOI.Utilities.normalize_constant( + _as_row(_demote(rows[k])), inner) for (k, inner) in + zip(row_indices(set, i), set.inner_sets[i])] + push!(disjuncts, _Disjunct(activation, binary, active_value, + MOI.AbstractScalarFunction[f for (f, _) in normalized], + MOI.AbstractScalarSet[s for (_, s) in normalized])) + 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) + +################################################################################ +# INNER SOLVER SUPPORT CHECK +################################################################################ +# Constraint types the master receives: variable constraints, linear +# rows, and the gated linear disjunct rows per `MasterReformulation`. +# The exactly-one rows, OA cuts, no-good cuts, and slack bounds are +# always affine rows or variable bounds. +function _master_constraint_types(model::Optimizer, problem::_Problem) + affine = MOI.ScalarAffineFunction{Float64} + required = Set{Tuple{Type, Type}}([ + (affine, MOI.EqualTo{Float64}), + (affine, MOI.LessThan{Float64}), + (affine, MOI.GreaterThan{Float64}), + (MOI.VariableIndex, MOI.GreaterThan{Float64}), + (MOI.VariableIndex, MOI.LessThan{Float64}), + ]) + for ci in problem.variable_cis + push!(required, (MOI.VariableIndex, + typeof(MOI.get(model.cache, MOI.ConstraintSet(), ci)))) + end + for ci in problem.linear_cis + push!(required, (affine, + typeof(MOI.get(model.cache, MOI.ConstraintSet(), ci)))) + end + for (func, set) in problem.nonlinear_rows + _is_linear(func) && push!(required, (affine, typeof(set))) + end + MOI.get(model, MasterReformulation()) == "indicator" || return required + for disjunction in problem.disjunctions, + disjunct in disjunction.disjuncts + activate = disjunct.active_value ? MOI.ACTIVATE_ON_ONE : + MOI.ACTIVATE_ON_ZERO + for (func, set) in zip(disjunct.functions, disjunct.sets) + _is_linear(func) || continue + for inner in _indicator_sets(set) + push!(required, (MOI.VectorAffineFunction{Float64}, + MOI.Indicator{activate, typeof(inner)})) + end + end + end + return required +end + +# Constraint types the NLP subproblems receive: every global and +# disjunct row plus the binary fixes and the NLPF slack variable. +function _nlp_constraint_types(model::Optimizer, problem::_Problem) + required = Set{Tuple{Type, Type}}([ + (MOI.VariableIndex, MOI.EqualTo{Float64}), + (MOI.VariableIndex, MOI.GreaterThan{Float64}), + ]) + indicators = Set(problem.binaries) + for ci in problem.variable_cis + vi = MOI.get(model.cache, MOI.ConstraintFunction(), ci) + vi in indicators && continue + push!(required, (MOI.VariableIndex, + typeof(MOI.get(model.cache, MOI.ConstraintSet(), ci)))) + end + for ci in problem.linear_cis + push!(required, (MOI.ScalarAffineFunction{Float64}, + typeof(MOI.get(model.cache, MOI.ConstraintSet(), ci)))) + end + for (func, set) in problem.nonlinear_rows + push!(required, (typeof(func), typeof(set))) + end + for disjunction in problem.disjunctions, + disjunct in disjunction.disjuncts + for (func, set) in zip(disjunct.functions, disjunct.sets) + push!(required, (typeof(func), typeof(set))) + end + end + return required +end + +function _check_support(solver, name::String, destination::String, required) + for (F, S) in required + MOI.supports_constraint(solver, F, S) || error( + "The `$name` ($(MOI.get(solver, MOI.SolverName()))) does " * + "not support `$F`-in-`$S` constraints, which the " * + "$destination requires.") + end + return +end + +# Fail before any subproblem work when an inner solver cannot take the +# constraint types routed to it, naming the solver and the type. A +# subproblem method brings its own solver, so only the `nothing` +# (MOI-path) method checks the `nlp_solver`. +function _check_inner_support(method, model::Optimizer, problem::_Problem) + mip = _instantiate(model.mip_solver) + _check_support(mip, "mip_solver", "master problem", + _master_constraint_types(model, problem)) + _check_nlp_support(method, model, problem) + return +end + +_check_nlp_support(method, ::Optimizer, ::_Problem) = nothing + +function _check_nlp_support(::Nothing, model::Optimizer, problem::_Problem) + nlp = _instantiate(model.nlp_solver) + _check_support(nlp, "nlp_solver", "NLP subproblems", + _nlp_constraint_types(model, problem)) + F = typeof(problem.objective) + MOI.supports(nlp, MOI.ObjectiveFunction{F}()) || error( + "The `nlp_solver` ($(MOI.get(nlp, MOI.SolverName()))) does " * + "not support the `$F` objective the NLP subproblems require.") + return +end diff --git a/test/integration.jl b/test/integration.jl new file mode 100644 index 0000000..c628981 --- /dev/null +++ b/test/integration.jl @@ -0,0 +1,120 @@ +using DisjunctiveProgramming, HiGHS, Ipopt, InfiniteOpt + +function _optimizer_factory() + return () -> DA.Optimizer(Ipopt.Optimizer, 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 = Direct()) + @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 = Direct()) + @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 = Direct()) + @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 = Direct()) + @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 = Direct()) + @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..a4f4264 --- /dev/null +++ b/test/loa.jl @@ -0,0 +1,1124 @@ +using JuMP +import HiGHS, Ipopt + +function _loa_optimizer(attrs::Pair...) + return optimizer_with_attributes( + () -> DA.Optimizer(Ipopt.Optimizer, HiGHS.Optimizer), attrs...) +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 + +# Same instance with function constants left in the rows, as +# InfiniteOpt transcription produces: x - 7 >= -5 and x - 5 >= 0. +# The parse must shift the constants into the sets or direct solver +# wrappers reject the subproblem rows. +function test_row_function_constants() + 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 - 7, x - 5] in DA.DisjunctionSet([ + [MOI.GreaterThan(-5.0)], [MOI.GreaterThan(0.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 2.0 atol = 1e-5 + @test value(z[1]) ≈ 1.0 atol = 1e-5 +end + +# Convex quadratic objective over linear disjuncts: y >= x or +# y >= 2 - x. The optimum sits at x = 3, y = 0 in disjunct 2. +# Promoted rows reach DA as operator trees in the shapes JuMP emits: +# a bare `x * y` and a `-` node. MOI's own quadratic `convert` throws +# on both; `_demote` must not. +function test_quadratic_row_shapes() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10, start = 1) + @variable(model, 0 <= y <= 10, start = 1) + @variable(model, z[1:2], Bin) + product = NonlinearExpr(:*, Any[x, y]) + shifted = NonlinearExpr(:-, Any[NonlinearExpr(:*, Any[x, y]), 3.0]) + @constraint(model, [1, z[1], z[2], product, shifted] in + DA.DisjunctionSet([ + [MOI.GreaterThan(4.0)], [MOI.GreaterThan(6.0)]])) + @objective(model, Min, x + y) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 4.0 atol = 1e-4 + @test value(z[1]) ≈ 1.0 atol = 1e-5 + @test value(x) * value(y) >= 4.0 - 1e-5 +end + +# min x + y over [x >= 3] v [x <= 1] and [y >= x^2 + 1] v [y >= exp(x)] +# with 0 <= x <= 4, 0 <= y <= 10: optimum 1 at (0, 1); the combination +# (x >= 3, y >= exp(x)) is infeasible because exp(3) > 10. +function _batched_test_model(attrs::Pair...) + model = Model(_loa_optimizer(attrs...)) + set_silent(model) + @variable(model, 0 <= x <= 4) + @variable(model, 0 <= y <= 10) + @variable(model, z[1:2], Bin) + @variable(model, w[1:2], Bin) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.GreaterThan(3.0)], [MOI.LessThan(1.0)]])) + @constraint(model, [1, w[1], w[2], y - x^2, y - exp(x)] in + DA.DisjunctionSet([[MOI.GreaterThan(1.0)], [MOI.GreaterThan(0.0)]])) + @objective(model, Min, x + y) + return model, x, y +end + +# The set-cover problem holds the indicators, their integrality, the +# exactly-one rows and the propositional rows; a row mixing continuous +# variables stays out, and LOA still reaches the optimum from it. +function test_set_cover_model() + model, x, y = _batched_test_model() + z, w = model[:z], model[:w] + @constraint(model, z[1] + w[1] >= 1) + @constraint(model, x + z[1] <= 10) + optimize!(model) + @test objective_value(model) ≈ 1.0 atol = 1e-4 + optimizer = unsafe_backend(model) + problem = DA._build_problem(optimizer) + set_cover = DA._build_set_cover(optimizer, problem) + cover = set_cover.model + @test MOI.get(cover, MOI.NumberOfVariables()) == 4 + @test MOI.get(cover, MOI.NumberOfConstraints{MOI.VariableIndex, + MOI.ZeroOne}()) == 4 + affine = MOI.ScalarAffineFunction{Float64} + @test MOI.get(cover, MOI.NumberOfConstraints{affine, + MOI.EqualTo{Float64}}()) == 2 + @test MOI.get(cover, MOI.NumberOfConstraints{affine, + MOI.GreaterThan{Float64}}()) == 1 + @test MOI.get(cover, MOI.NumberOfConstraints{affine, + MOI.LessThan{Float64}}()) == 0 + combination = Dict(JuMP.index(z[1]) => true, JuMP.index(z[2]) => false, + JuMP.index(w[1]) => true, JuMP.index(w[2]) => false) + DA._avoid_combination(set_cover, combination) + @test MOI.get(cover, MOI.NumberOfConstraints{affine, + MOI.GreaterThan{Float64}}()) == 2 +end + +# every combination that activates exactly one disjunct per disjunction +function _all_combinations(problem::DA._Problem) + choices = Iterators.product((disjunction.disjuncts + for disjunction in problem.disjunctions)...) + return map(vec(collect(choices))) do active + combination = Dict{MOI.VariableIndex, Bool}() + for disjunction in problem.disjunctions, + disjunct in disjunction.disjuncts + combination[disjunct.binary] = disjunct in active ? + disjunct.active_value : !disjunct.active_value + end + return combination + end +end + +function test_batched_subproblems_loa() + model, x, y = _batched_test_model( + DA.SubproblemMethod() => DA.BatchedSubproblems()) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 1.0 atol = 1e-4 + @test value(x) ≈ 0.0 atol = 1e-4 + @test value(y) ≈ 1.0 atol = 1e-4 +end + +# The stacked solve must reproduce the sequential results, and a batch +# with an infeasible copy must fall back to them exactly. +function test_batched_matches_sequential() + model, _, _ = _batched_test_model() + optimize!(model) + optimizer = unsafe_backend(model) + problem = DA._build_problem(optimizer) + batched = DA.BatchedSubproblems() + sub = DA._build_subproblem(batched, optimizer, problem) + combinations = _all_combinations(problem) + @test length(combinations) == 4 + sequential = DA._solve_nlps(nothing, optimizer, problem, sub.sequential, + combinations, nothing) + @test count(r -> r.feasible, sequential) == 3 + feasible = [c for (c, r) in zip(combinations, sequential) if r.feasible] + stacked = DA._solve_nlps(batched, optimizer, problem, sub, feasible, + nothing) + expected = [r for r in sequential if r.feasible] + @test all(r.feasible for r in stacked) + @test all(r.status == MOI.LOCALLY_SOLVED for r in stacked) + for (r, e) in zip(stacked, expected) + @test r.combination === e.combination + @test r.objective ≈ e.objective atol = 1e-5 + for (vi, value) in e.point + @test r.point[vi] ≈ value atol = 1e-4 + end + end + # the stacked feasibility pass classifies the copies and gives the + # infeasible one a restoration point instead of dropping it + mixed = DA._solve_nlps(batched, optimizer, problem, sub, combinations, + nothing) + for (r, e) in zip(mixed, sequential) + @test r.feasible == e.feasible + @test r.point !== nothing + r.feasible && @test r.objective ≈ e.objective atol = 1e-5 + r.feasible || @test r.status == MOI.LOCALLY_INFEASIBLE + end +end + +# Two combinations per master solve through the no-good re-solve path +# (HiGHS has no solution pool); same optimum, fewer master solves than +# NLP solves. +function test_multi_generation_cuts() + model, x, y = _batched_test_model(DA.MultiGenerationSize() => 2) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 1.0 atol = 1e-4 + @test value(x) ≈ 0.0 atol = 1e-4 + optimizer = unsafe_backend(model) + @test MOI.get(optimizer, DA.NLPSolveCount()) >= 2 + @test MOI.get(optimizer, DA.MasterSolveCount()) >= 1 + @test MOI.get(model, DA.NLPSolveCount()) == + MOI.get(optimizer, DA.NLPSolveCount()) + @test 0 <= MOI.get(model, DA.NLPInfeasibleCount()) <= + MOI.get(model, DA.NLPSolveCount()) + master_time = MOI.get(model, DA.MasterSolveTime()) + nlp_time = MOI.get(model, DA.NLPSolveTime()) + @test master_time > 0 && nlp_time > 0 + @test master_time + nlp_time <= MOI.get(model, MOI.SolveTimeSec()) + stacked, _, _ = _batched_test_model(DA.MultiGenerationSize() => 2, + DA.SubproblemMethod() => DA.BatchedSubproblems()) + optimize!(stacked) + @test objective_value(stacked) ≈ 1.0 atol = 1e-4 + @test_throws ErrorException MOI.set(DA.LOA(), DA.MultiGenerationSize(), 0) + @test_throws ErrorException MOI.set(DA.LOA(), DA.MultiGenerationSize(), + 1.5) +end + +# Every source reaches the optimum through a three-per-master loop, and +# the neighborhood of a proposal on the 2 x 2 model is its two flips. +function test_combination_sources() + for source in (DA.Neighborhood(), DA.Neighborhood(around = :incumbent), + DA.LPRounding(seed = 1), DA.RandomCombinations(seed = 1), + DA.CutoffResolve(tolerance = 0.5)) + model, x, y = _batched_test_model(DA.MultiGenerationSize() => 3, + DA.CombinationSource() => source) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 1.0 atol = 1e-4 + @test value(x) ≈ 0.0 atol = 1e-4 + end + model, _, _ = _batched_test_model() + optimize!(model) + optimizer = unsafe_backend(model) + problem = DA._build_problem(optimizer) + master = DA._build_master(optimizer, problem) + proposal = _all_combinations(problem)[1] + neighbors, excluded = DA._candidate_combinations(DA.Neighborhood(), + optimizer, problem, master, proposal, nothing, 10, Inf) + @test !excluded + @test length(neighbors) == 2 + @test all(n != proposal for n in neighbors) + @test all(count(n[b] != proposal[b] for b in keys(proposal)) == 2 + for n in neighbors) + two, _ = DA._candidate_combinations(DA.Neighborhood(), optimizer, + problem, master, proposal, nothing, 1, Inf) + @test length(two) == 1 + @test_throws ErrorException DA.Neighborhood(around = :elsewhere) + @test_throws ErrorException DA.CutoffResolve(tolerance = -1) +end + +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 + +# `UseNLPF() => false` still solves by enumeration when the seeds +# are feasible. +function test_nlpf_disabled() + model = Model(_loa_optimizer(DA.UseNLPF() => 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_reformulation() + model = Model(_loa_optimizer(DA.MasterReformulation() => "bigm", + DA.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 attribute values convert on set, including the Bool +# `UseNLPF` on the infeasible-seed path. +function test_integer_options() + model = Model(_loa_optimizer(DA.UseNLPF() => 0, DA.M_Value() => 10^9, + MOI.TimeLimitSec() => 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(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(MOI.TimeLimitSec() => 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 + +# `UseNLPF() => 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(DA.UseNLPF() => 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 + +# `SetCoverProblem() => MasterCover()` covers from the full master: on +# a problem whose linear rows rule out most combinations (x + y >= 12 +# with x in a disjunct), the logic-only cover draws a combination whose +# NLP is infeasible, while the master cover lands on a feasible one and +# the run ends with the same incumbent. +function test_set_cover_problem_master() + for cover in (DA.LogicOnlyCover(), DA.MasterCover()) + model = Model(_loa_optimizer(DA.SetCoverProblem() => cover)) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, 0 <= y <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.LessThan(1.0)], [MOI.GreaterThan(4.0)]])) + @constraint(model, x + y >= 12) + @objective(model, Min, x^2 + y) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 24.0 atol = 1e-4 + @test value(z[2]) ≈ 1.0 atol = 1e-5 + infeasible = MOI.get(unsafe_backend(model), DA.NLPInfeasibleCount()) + cover isa DA.MasterCover && @test infeasible == 0 + end + @test_throws ErrorException Model(_loa_optimizer( + DA.SetCoverProblem() => "master")) +end + +# Zero iteration budgets exit before any solve, without an incumbent. +function test_iteration_limit_no_incumbent() + model = Model(_loa_optimizer(DA.SetCoverIterationLimit() => 0, + DA.NumIterationLimit() => 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 + +# `NumIterationLimit() => 0` keeps the set-covering incumbent but +# produces no master bound. +function test_iteration_limit_with_incumbent() + model = Model(_loa_optimizer(DA.NumIterationLimit() => 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(DA.SetCoverIterationLimit() => 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 + +# An unbounded NLP restriction proves the model unbounded: the solve +# reports `DUAL_INFEASIBLE` instead of an infeasibility or limit +# status. +function test_unbounded_model() + factory = optimizer_with_attributes( + () -> DA.Optimizer(HiGHS.Optimizer, HiGHS.Optimizer), + DA.MasterReformulation() => "bigm") + model = Model(factory) + set_silent(model) + @variable(model, x >= 0) + @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, Max, x) + optimize!(model) + @test termination_status(model) == MOI.DUAL_INFEASIBLE + @test result_count(model) == 0 + @test occursin("unbounded", raw_status(model)) +end + +# An NLP failure without an infeasibility certificate leaves its +# combination unresolved. Both disjuncts are nonlinear so the cover +# pass must visit both combinations: the first solves and gives the +# incumbent, the second fails with `NUMERICAL_ERROR`, and exhaustion +# must not claim local optimality over all combinations. +function test_nlp_failure_unresolved() + factory = () -> DA.Optimizer( + () -> MockSolver(Ipopt.Optimizer; + fail_from = 2, fail_status = MOI.NUMERICAL_ERROR), + HiGHS.Optimizer) + 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^2, x^2] in DA.DisjunctionSet([ + [MOI.GreaterThan(4.0)], [MOI.GreaterThan(25.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.OTHER_LIMIT + @test primal_status(model) == MOI.FEASIBLE_POINT + @test occursin("1 unresolved", raw_status(model)) +end + +# Bad attribute values fail at set time, not mid-solve. +function test_attribute_validation() + algorithm = DA.LOA() + @test_throws ErrorException MOI.set(algorithm, DA.OASlack(), 0.0) + @test_throws ErrorException MOI.set(algorithm, DA.OASlack(), -1.0) + @test_throws ErrorException MOI.set(algorithm, + DA.MasterReformulation(), "hull") + MOI.set(algorithm, DA.OASlack(), 5.0) + @test MOI.get(algorithm, DA.OASlack()) == 5.0 + MOI.set(algorithm, DA.MasterReformulation(), "bigm") + @test MOI.get(algorithm, DA.MasterReformulation()) == "bigm" +end + +# A genuinely nonlinear indicator row gets the curated error, not a +# raw conversion failure. +function test_nonlinear_indicator_error() + 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], z[2], x, x] in DA.DisjunctionSet([ + [MOI.GreaterThan(2.0)], [MOI.GreaterThan(5.0)]])) + @objective(model, Min, x) + err = try + optimize!(model) + nothing + catch e + e + end + @test err isa ErrorException + @test occursin("Unsupported indicator expression", err.msg) +end + +function _mock_time_limit_model(limit::Float64; kwargs...) + factory = optimizer_with_attributes( + () -> DA.Optimizer(() -> MockSolver(Ipopt.Optimizer; kwargs...), + HiGHS.Optimizer), + DA.IterationTimeLimit() => 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 + +# Strict mocks on both roles run the whole loop (master, NLP, NLPF, cuts) +# against a direct-mode solver contract: rows arrive with function +# constants and the loop must never pass one through. +function test_strict_constant_solvers() + factory = () -> DA.Optimizer( + () -> MockSolver(Ipopt.Optimizer, strict_constants = true), + () -> MockSolver(HiGHS.Optimizer, strict_constants = true)) + 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 - 7, x - 5] in DA.DisjunctionSet([ + [MOI.GreaterThan(-5.0)], [MOI.GreaterThan(0.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 2.0 atol = 1e-5 + @test value(z[1]) ≈ 1.0 atol = 1e-5 +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 second +# MIP solve of the run, the first master solve after the covering NLP, +# reports a node limit. +function test_master_abnormal_status_with_incumbent() + solves = Ref(0) + mip = () -> MockSolver(HiGHS.Optimizer, fail_from = 2, + shared_solves = solves) + factory = () -> DA.Optimizer(Ipopt.Optimizer, mip) + 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 + +# The upfront capability check fails before any subproblem work and +# names the offending inner solver and constraint type. +function test_inner_solver_support_check() + # HiGHS as the nlp_solver cannot take the quadratic disjunct row + model = Model(() -> DA.Optimizer(HiGHS.Optimizer)) + set_silent(model) + @variable(model, 0 <= x <= 4) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], x^2, x] in DA.DisjunctionSet([ + [MOI.LessThan(4.0)], [MOI.LessThan(1.0)]])) + @objective(model, Min, x) + err = try + optimize!(model) + nothing + catch e + e + end + @test err isa ErrorException + @test occursin("nlp_solver", err.msg) + @test occursin("ScalarQuadraticFunction", err.msg) + # Ipopt as the mip_solver cannot take the master's discrete parts + model = Model(() -> DA.Optimizer(Ipopt.Optimizer)) + set_silent(model) + @variable(model, 0 <= y <= 10) + @variable(model, w[1:2], Bin) + @constraint(model, [1, w[1], w[2], y, y] in DA.DisjunctionSet([ + [MOI.GreaterThan(2.0)], [MOI.GreaterThan(5.0)]])) + @objective(model, Min, y) + err = try + optimize!(model) + nothing + catch e + e + end + @test err isa ErrorException + @test occursin("mip_solver", err.msg) +end + +# A `SubproblemMethod` object routes the support check, the build, and +# every NLP solve through its own overloads; delegating back to the +# `nothing` path must reproduce the MOI-path answer exactly. +struct _RecordingSubproblems + log::Vector{Symbol} +end + +DA._check_nlp_support(m::_RecordingSubproblems, ::DA.Optimizer, + ::DA._Problem) = push!(m.log, :check) + +function DA._build_subproblem(m::_RecordingSubproblems, model::DA.Optimizer, + problem::DA._Problem) + push!(m.log, :build) + return DA._build_subproblem(nothing, model, problem) +end + +function DA._solve_nlp(m::_RecordingSubproblems, model::DA.Optimizer, + problem::DA._Problem, sub, combination, warm_start; + deadline::Float64 = Inf) + push!(m.log, :solve) + return DA._solve_nlp(nothing, model, problem, sub, combination, + warm_start; deadline) +end + +function test_subproblem_method_dispatch() + method = _RecordingSubproblems(Symbol[]) + model = Model(_loa_optimizer(DA.SubproblemMethod() => method)) + 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 + @test method.log[1] == :check + @test method.log[2] == :build + @test count(==(:solve), method.log) >= 1 + @test MOI.get(unsafe_backend(model), DA.SubproblemMethod()) === method +end + +@testset "LOA loop" begin + test_linear_disjunction() + test_row_function_constants() + test_nested_disjunction() + test_nested_disjunction_vacuous() + test_quadratic_row_shapes() + test_batched_subproblems_loa() + test_batched_matches_sequential() + test_multi_generation_cuts() + test_combination_sources() + 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_reformulation() + 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_unbounded_model() + test_nlp_failure_unresolved() + test_attribute_validation() + test_nonlinear_indicator_error() + test_strict_constant_solvers() + 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() + test_inner_solver_support_check() + test_subproblem_method_dispatch() +end + +@testset "LOA units" begin + test_set_cover_model() + test_set_cover_problem_master() + test_cut_term_directions() + test_activation_binary() + test_sense_primitives() + test_linearize_quadratic() +end diff --git a/test/mock_optimizer.jl b/test/mock_optimizer.jl new file mode 100644 index 0000000..1716a47 --- /dev/null +++ b/test/mock_optimizer.jl @@ -0,0 +1,122 @@ +################################################################################ +# 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. +# `shared_solves` counts solves across every instance handed the same +# `Ref`, so `fail_from` can address the n-th MIP solve of a run that +# instantiates the factory more than once (support probe, master, +# set-cover model). +# Deterministic triggers for the deadline and abnormal-master paths. +# `strict_constants` imitates direct wrappers (e.g. Gurobi) that reject +# scalar functions with nonzero constants. +mutable struct MockSolver <: MOI.AbstractOptimizer + inner::MOI.AbstractOptimizer + sleep_time::Float64 + fail_from::Int + fail_status::MOI.TerminationStatusCode + strict_constants::Bool + solves::Int + shared_solves::Ref{Int} + failing::Bool +end + +function MockSolver( + factory; + sleep_time::Float64 = 0.0, + fail_from::Int = typemax(Int), + fail_status::MOI.TerminationStatusCode = MOI.NODE_LIMIT, + strict_constants::Bool = false, + shared_solves::Ref{Int} = Ref(0) + ) + return MockSolver(MOI.instantiate(factory), sleep_time, fail_from, + fail_status, strict_constants, 0, shared_solves, false) +end + +function MOI.optimize!(model::MockSolver) + model.solves += 1 + model.shared_solves[] += 1 + model.sleep_time > 0 && sleep(model.sleep_time) + model.failing = model.shared_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::_WrappedIndex) = + MOI.delete(model.inner, index) +MOI.is_valid(model::MockSolver, index::_WrappedIndex) = + MOI.is_valid(model.inner, index) + +function MOI.add_constraint( + model::MockSolver, + func::MOI.AbstractFunction, + set::MOI.AbstractSet + ) + if model.strict_constants && func isa Union{ + MOI.ScalarAffineFunction{Float64}, + MOI.ScalarQuadraticFunction{Float64}} && !iszero(func.constant) + throw(MOI.ScalarFunctionConstantNotZero{ + Float64, typeof(func), typeof(set)}(func.constant)) + end + 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 diff --git a/test/moi.jl b/test/moi.jl new file mode 100644 index 0000000..68b1db2 --- /dev/null +++ b/test/moi.jl @@ -0,0 +1,34 @@ +# 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( + MOI.OptimizerWithAttributes( + () -> DA.Optimizer(Ipopt.Optimizer, HiGHS.Optimizer), + MOI.TimeLimitSec() => 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..2fd510f --- /dev/null +++ b/test/optimizer.jl @@ -0,0 +1,152 @@ +# 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(nothing) + @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(nothing) + @test MOI.get(optimizer, DA.Algorithm()) === nothing + @test MOI.supports(optimizer, DA.Algorithm()) + @test MOI.supports(optimizer, DA.NumIterationLimit()) + @test MOI.get(optimizer, DA.NumIterationLimit()) == 10 + # setting an attribute materializes the default algorithm + MOI.set(optimizer, DA.NumIterationLimit(), 5) + @test MOI.get(optimizer, DA.NumIterationLimit()) == 5 + @test MOI.get(optimizer, DA.Algorithm()) isa DA.LOA + @test MOI.get(optimizer, DA.M_Value()) == 1e9 + algorithm = DA.LOA() + MOI.set(optimizer, DA.Algorithm(), algorithm) + @test MOI.get(optimizer, DA.Algorithm()) === algorithm + @test MOI.get(optimizer, DA.NumIterationLimit()) == 10 + @test !MOI.supports(optimizer, MOI.RawOptimizerAttribute("max_iter")) + MOI.set(optimizer, MOI.Silent(), true) + @test MOI.get(optimizer, MOI.Silent()) + @test MOI.get(optimizer, MOI.TimeLimitSec()) == 3600.0 + MOI.set(optimizer, MOI.TimeLimitSec(), 100) + @test MOI.get(optimizer, MOI.TimeLimitSec()) == 100.0 + MOI.set(optimizer, MOI.TimeLimitSec(), nothing) + @test MOI.get(optimizer, MOI.TimeLimitSec()) === nothing +end + +function test_model_building() + optimizer = DA.Optimizer(nothing) + 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(nothing) + 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(nothing) + @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(nothing) + 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 + +# Only the scalar objective functions the solve demotes are supported; +# a vector objective is refused at set time. The solver version comes +# from Project.toml, not a copy. +function test_objective_function_support() + optimizer = DA.Optimizer(nothing) + @test MOI.supports(optimizer, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()) + @test MOI.supports(optimizer, + MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}()) + @test !MOI.supports(optimizer, + MOI.ObjectiveFunction{MOI.VectorAffineFunction{Float64}}()) + x = MOI.add_variable(optimizer) + func = MOI.Utilities.operate(vcat, Float64, 1.0 * x, 2.0 * x) + attr = MOI.ObjectiveFunction{typeof(func)}() + @test_throws MOI.UnsupportedAttribute MOI.set(optimizer, attr, func) + @test MOI.get(optimizer, MOI.SolverVersion()) == + string(pkgversion(DA)) +end + +@testset "Optimizer scaffold" begin + test_optimizer_scaffold() + test_optimizer_options() + test_model_building() + test_copy_to() + test_moi_forwarding() + test_unsupported_constraint_type() + test_objective_function_support() +end diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..50fb40b --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,10 @@ +using Test +import MathOptInterface as MOI +using DisjunctiveAlgorithms +const DA = DisjunctiveAlgorithms + +include("optimizer.jl") +include("mock_optimizer.jl") +include("loa.jl") +include("moi.jl") +include("integration.jl")