From cbaf2c50f45333695a9e2c85b0b8f1fc6a15b1ee Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Tue, 11 Aug 2026 12:59:09 -0400 Subject: [PATCH 01/18] Initial package --- .github/dependabot.yml | 24 + .github/workflows/CI.yml | 72 +++ .github/workflows/CompatHelper.yml | 18 + .github/workflows/TagBot.yml | 19 + Project.toml | 23 + README.md | 40 +- docs/Project.toml | 3 + docs/make.jl | 22 + docs/src/index.md | 14 + src/DisjunctiveAlgorithms.jl | 15 + src/cuts.jl | 202 ++++++++ src/loa.jl | 257 ++++++++++ src/master.jl | 181 +++++++ src/nlp.jl | 168 ++++++ src/optimizer.jl | 380 ++++++++++++++ src/problem.jl | 159 ++++++ test/integration.jl | 121 +++++ test/loa.jl | 791 +++++++++++++++++++++++++++++ test/moi.jl | 35 ++ test/optimizer.jl | 127 +++++ test/runtests.jl | 9 + 21 files changed, 2679 insertions(+), 1 deletion(-) create mode 100644 .github/dependabot.yml create mode 100644 .github/workflows/CI.yml create mode 100644 .github/workflows/CompatHelper.yml create mode 100644 .github/workflows/TagBot.yml create mode 100644 Project.toml create mode 100644 docs/Project.toml create mode 100644 docs/make.jl create mode 100644 docs/src/index.md create mode 100644 src/DisjunctiveAlgorithms.jl create mode 100644 src/cuts.jl create mode 100644 src/loa.jl create mode 100644 src/master.jl create mode 100644 src/nlp.jl create mode 100644 src/optimizer.jl create mode 100644 src/problem.jl create mode 100644 test/integration.jl create mode 100644 test/loa.jl create mode 100644 test/moi.jl create mode 100644 test/optimizer.jl create mode 100644 test/runtests.jl diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000..5f06744 --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,24 @@ +# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates +version: 2 +updates: + - package-ecosystem: "github-actions" + directory: "/" + schedule: + interval: "weekly" + # To group all GitHub Actions updates into a single PR, uncomment the following: + # groups: + # github-actions: + # patterns: + # - "*" + - package-ecosystem: "julia" + directories: + - "/" + - "/docs" + - "/test" + schedule: + interval: "weekly" + # To group all Julia dependency updates into a single PR, uncomment the following: + # groups: + # julia-dependencies: + # patterns: + # - "*" diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml new file mode 100644 index 0000000..2e521a4 --- /dev/null +++ b/.github/workflows/CI.yml @@ -0,0 +1,72 @@ +name: CI +on: + push: + branches: + - main + tags: ['*'] + pull_request: + workflow_dispatch: +concurrency: + # Skip intermediate builds: always. + # Cancel intermediate builds: only if it is a pull request build. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} +jobs: + test: + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} + runs-on: ${{ matrix.os }} + timeout-minutes: 60 + permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created + actions: write + contents: read + strategy: + fail-fast: false + matrix: + version: + - '1.10' + - '1.11' + - 'pre' + os: + - ubuntu-latest + arch: + - x64 + steps: + - uses: actions/checkout@v6 + - uses: julia-actions/setup-julia@v3 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v2 + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-runtest@v1 + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v6 + with: + files: lcov.info + token: ${{ secrets.CODECOV_TOKEN }} + fail_ci_if_error: false + docs: + name: Documentation + runs-on: ubuntu-latest + permissions: + actions: write # needed to allow julia-actions/cache to proactively delete old caches that it has created + contents: write + statuses: write + steps: + - uses: actions/checkout@v6 + - uses: julia-actions/setup-julia@v3 + with: + version: '1' + - uses: julia-actions/cache@v2 + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-docdeploy@v1 + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} + - name: Run doctests + shell: julia --project=docs --color=yes {0} + run: | + using Documenter: DocMeta, doctest + using DisjunctiveAlgorithms + DocMeta.setdocmeta!(DisjunctiveAlgorithms, :DocTestSetup, :(using DisjunctiveAlgorithms); recursive=true) + doctest(DisjunctiveAlgorithms) diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml new file mode 100644 index 0000000..ff25ebc --- /dev/null +++ b/.github/workflows/CompatHelper.yml @@ -0,0 +1,18 @@ +name: CompatHelper +on: + schedule: + - cron: 0 0 * * * + workflow_dispatch: +jobs: + CompatHelper: + runs-on: ubuntu-latest + steps: + - name: Install CompatHelper + run: using Pkg; Pkg.add("CompatHelper") + shell: julia --color=yes {0} + - name: Run CompatHelper + run: using CompatHelper; CompatHelper.main() + shell: julia --color=yes {0} + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }} diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml new file mode 100644 index 0000000..3f42eb1 --- /dev/null +++ b/.github/workflows/TagBot.yml @@ -0,0 +1,19 @@ +name: TagBot +on: + issue_comment: + types: + - created + workflow_dispatch: +jobs: + TagBot: + if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot' + runs-on: ubuntu-latest + steps: + - uses: JuliaRegistries/TagBot@v1 + with: + token: ${{ secrets.GITHUB_TOKEN }} + # For commits that modify workflow files: SSH key enables tagging, but + # releases require manual creation. For full automation of such commits, + # use a PAT with `workflow` scope instead of GITHUB_TOKEN. + # See: https://github.com/JuliaRegistries/TagBot#commits-that-modify-workflow-files + ssh: ${{ secrets.DOCUMENTER_KEY }} diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..ff640ed --- /dev/null +++ b/Project.toml @@ -0,0 +1,23 @@ +name = "DisjunctiveAlgorithms" +uuid = "af3b8e3a-4098-461f-b6c0-2c9706adb26e" +authors = ["Daniel Nguyen"] +version = "0.1.0" + +[deps] +DisjunctiveProgramming = "0d27d021-0159-4c7d-b4a7-9ccb5d9366cf" +MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" + +[compat] +DisjunctiveProgramming = "0.6" +MathOptInterface = "1" +julia = "1.10" + +[extras] +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" +InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["HiGHS", "InfiniteOpt", "Ipopt", "JuMP", "Test"] diff --git a/README.md b/README.md index d6ac54f..1d82799 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,40 @@ # DisjunctiveAlgorithms.jl -An optimizer suite for generalized disjunctive programming. + +[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://infiniteopt.github.io/DisjunctiveAlgorithms.jl/stable/) +[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://infiniteopt.github.io/DisjunctiveAlgorithms.jl/dev/) +[![Build Status](https://github.com/infiniteopt/DisjunctiveAlgorithms.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/infiniteopt/DisjunctiveAlgorithms.jl/actions/workflows/CI.yml?query=branch%3Amain) +[![Coverage](https://codecov.io/gh/infiniteopt/DisjunctiveAlgorithms.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/infiniteopt/DisjunctiveAlgorithms.jl) + +An optimizer suite for generalized disjunctive programming (GDP). + +DisjunctiveAlgorithms.jl is an MOI-layer solver for models that contain +disjunctions encoded as vector constraints in +`DisjunctiveProgramming.DisjunctionSet`. Disjunction-aware algorithms +(currently logic-based outer approximation) solve the model by +dispatching subproblems to user-provided MIP and NLP solvers. + +The design follows +[MultiObjectiveAlgorithms.jl](https://github.com/jump-dev/MultiObjectiveAlgorithms.jl): +one `Optimizer` that wraps inner solvers, with the algorithm and its +options selected through optimizer attributes. + +## Usage with DisjunctiveProgramming.jl + +```julia +using DisjunctiveProgramming, DisjunctiveAlgorithms, HiGHS, Ipopt +import DisjunctiveAlgorithms as DA + +model = GDPModel(() -> DA.Optimizer(nlp_solver = Ipopt.Optimizer, + mip_solver = HiGHS.Optimizer)) +@variable(model, 0 <= x <= 10) +@variable(model, Y[1:2], Logical) +@constraint(model, x <= 3, Disjunct(Y[1])) +@constraint(model, x^2 == 64, Disjunct(Y[2])) +@disjunction(model, Y) +@objective(model, Max, x) +optimize!(model, gdp_method = MOIDisjunction()) +``` + +`MOIDisjunction()` lowers each disjunction to a single +`DisjunctionSet` constraint that this package consumes directly; no +Big-M or Hull reformulation is performed on the modeling side. diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..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..6831d27 --- /dev/null +++ b/src/DisjunctiveAlgorithms.jl @@ -0,0 +1,15 @@ +module DisjunctiveAlgorithms + +import MathOptInterface as MOI +import DisjunctiveProgramming: DisjunctionSet, num_disjuncts, + activation_index, indicator_indices, row_indices, + _SupportedInnerSet + +include("optimizer.jl") +include("problem.jl") +include("master.jl") +include("nlp.jl") +include("cuts.jl") +include("loa.jl") + +end diff --git a/src/cuts.jl b/src/cuts.jl new file mode 100644 index 0000000..079cc2e --- /dev/null +++ b/src/cuts.jl @@ -0,0 +1,202 @@ +################################################################################ +# LINEARIZATION +################################################################################ +# First-order Taylor expansions of the nonlinear rows, sharing one +# `MOI.Nonlinear` reverse-mode evaluator per row across iterations +# (only the evaluation point changes). +struct _Linearizer + evaluators::Dict{UInt64, + Tuple{MOI.Nonlinear.Evaluator, Vector{MOI.VariableIndex}}} +end + +_Linearizer() = _Linearizer(Dict{UInt64, + Tuple{MOI.Nonlinear.Evaluator, Vector{MOI.VariableIndex}}}()) + +function _append_variables( + variables::Vector{MOI.VariableIndex}, + func::MOI.VariableIndex + ) + return push!(variables, func) +end +function _append_variables( + variables::Vector{MOI.VariableIndex}, + func::MOI.ScalarAffineFunction{Float64} + ) + return append!(variables, term.variable for term in func.terms) +end +function _append_variables( + variables::Vector{MOI.VariableIndex}, + func::MOI.ScalarQuadraticFunction{Float64} + ) + append!(variables, term.variable for term in func.affine_terms) + for term in func.quadratic_terms + push!(variables, term.variable_1, term.variable_2) + end + return variables +end +function _append_variables( + variables::Vector{MOI.VariableIndex}, + func::MOI.ScalarNonlinearFunction + ) + for arg in func.args + arg isa Real || _append_variables(variables, arg) + end + return variables +end + +function _evaluator(linearizer::_Linearizer, func) + return get!(linearizer.evaluators, objectid(func)) do + variables = _append_variables(MOI.VariableIndex[], func) + unique!(variables) + nonlinear = MOI.Nonlinear.Model() + MOI.Nonlinear.set_objective(nonlinear, func) + evaluator = MOI.Nonlinear.Evaluator(nonlinear, + MOI.Nonlinear.SparseReverseMode(), variables) + MOI.initialize(evaluator, [:Grad]) + return (evaluator, variables) + end +end + +# Exact for affine functions, first-order Taylor at `point` otherwise. +# The result stays in cache space; callers remap when adding cuts. +_linearize(::_Linearizer, func::MOI.ScalarAffineFunction{Float64}, point) = func +_linearize(::_Linearizer, func::MOI.VariableIndex, point) = _to_affine(func) +function _linearize( + linearizer::_Linearizer, + func::Union{MOI.ScalarQuadraticFunction{Float64}, + MOI.ScalarNonlinearFunction}, + point::AbstractDict + ) + evaluator, variables = _evaluator(linearizer, func) + x = [point[vi] for vi in variables] + value = MOI.eval_objective(evaluator, x) + gradient = zeros(length(x)) + MOI.eval_objective_gradient(evaluator, gradient, x) + constant = value - sum(gradient[i] * x[i] for i in eachindex(x); + init = 0.0) + terms = [MOI.ScalarAffineTerm(gradient[i], variables[i]) + for i in eachindex(variables) if gradient[i] != 0.0] + return MOI.ScalarAffineFunction(terms, constant) +end + +################################################################################ +# OA CUT EMISSION +################################################################################ +_penalty_sign(sense::MOI.OptimizationSense) = sense == MOI.MAX_SENSE ? -1 : 1 + +# The `<= 0` directions of an OA cut for `set`: `lin - rhs` for +# LessThan, `rhs - lin` for GreaterThan, both for EqualTo / Interval. +_oa_cut_terms(set::MOI.LessThan{Float64}, lin) = + (MOI.Utilities.operate(-, Float64, lin, set.upper),) +_oa_cut_terms(set::MOI.GreaterThan{Float64}, lin) = + (MOI.Utilities.operate(-, Float64, set.lower, lin),) +_oa_cut_terms(set::MOI.EqualTo{Float64}, lin) = + (MOI.Utilities.operate(-, Float64, lin, set.value), + MOI.Utilities.operate(-, Float64, set.value, lin)) +_oa_cut_terms(set::MOI.Interval{Float64}, lin) = + (MOI.Utilities.operate(-, Float64, lin, set.upper), + MOI.Utilities.operate(-, Float64, set.lower, lin)) + +# Emit all OA cuts for one NLP result: the objective cut, a slacked row +# per nonlinear global, and a gated cut per active nonlinear disjunct +# row. Every slack keeps a nonconvex linearization from making the +# master infeasible. +function _add_oa_cuts( + model::Optimizer, + problem::_Problem, + master::_Master, + linearizer::_Linearizer, + result::NamedTuple + ) + result.point === nothing && return + sign = _penalty_sign(master.sense) + _add_objective_cut(model, master, linearizer, result.point, sign) + for (func, set) in problem.nonlinear_rows + _is_linear(func) && continue + lin = _master_linearization(master, linearizer, func, result.point) + _add_global_oa_row(model, master, lin, set, sign) + end + for disjunction in problem.disjunctions, disjunct in disjunction.disjuncts + _disjunct_active(result.combination, disjunct) || continue + for (func, set) in zip(disjunct.functions, disjunct.sets) + _is_linear(func) && continue + lin = _master_linearization(master, linearizer, func, result.point) + _add_disjunct_oa_cut(model, master, disjunct, lin, set, sign) + end + end + return +end + +function _master_linearization( + master::_Master, + linearizer::_Linearizer, + func, + point + ) + return _map_to(master.variable_map, _linearize(linearizer, func, point)) +end + +# Slacked objective cut. MIN: `lin <= alpha_oa + slack`; MAX symmetric. +function _add_objective_cut( + model::Optimizer, + master::_Master, + linearizer::_Linearizer, + point, + sign::Int + ) + lin = _master_linearization(master, linearizer, master.objective, point) + slack = _add_penalized_slack(master, model.options, sign) + alpha = _to_affine(master.alpha_oa) + if master.sense == MOI.MAX_SENSE + body = MOI.Utilities.operate(-, Float64, + MOI.Utilities.operate(-, Float64, alpha, lin), _to_affine(slack)) + else + body = MOI.Utilities.operate(-, Float64, + MOI.Utilities.operate(-, Float64, lin, alpha), _to_affine(slack)) + end + MOI.Utilities.normalize_and_add_constraint(master.model, body, + MOI.LessThan(0.0)) + return +end + +# Slacked global OA row(s): each direction `term - slack <= 0`. +function _add_global_oa_row( + model::Optimizer, + master::_Master, + lin::MOI.ScalarAffineFunction{Float64}, + set::MOI.AbstractScalarSet, + sign::Int + ) + slack = _add_penalized_slack(master, model.options, sign) + for term in _oa_cut_terms(set, lin) + body = MOI.Utilities.operate(-, Float64, term, _to_affine(slack)) + MOI.Utilities.normalize_and_add_constraint(master.model, body, + MOI.LessThan(0.0)) + end + return +end + +# Gated disjunct cut: each direction `term - slack <= M * (1 - z)` for +# an activation `z` (its complement gates on `M * z`). +function _add_disjunct_oa_cut( + model::Optimizer, + master::_Master, + disjunct::_Disjunct, + lin::MOI.ScalarAffineFunction{Float64}, + set::MOI.AbstractScalarSet, + sign::Int + ) + M = Float64(_option(model, "M_value")) + activation = _map_to(master.variable_map, disjunct.activation) + # M * (1 - activation) moved left: term - slack + M * activation - M + gate = MOI.Utilities.operate(-, Float64, + MOI.Utilities.operate(*, Float64, M, activation), M) + slack = _add_penalized_slack(master, model.options, sign) + for term in _oa_cut_terms(set, lin) + body = MOI.Utilities.operate(+, Float64, + MOI.Utilities.operate(-, Float64, term, _to_affine(slack)), gate) + MOI.Utilities.normalize_and_add_constraint(master.model, body, + MOI.LessThan(0.0)) + end + return +end diff --git a/src/loa.jl b/src/loa.jl new file mode 100644 index 0000000..d68fdd6 --- /dev/null +++ b/src/loa.jl @@ -0,0 +1,257 @@ +################################################################################ +# LOGIC-BASED OUTER APPROXIMATION +################################################################################ +_worst_objective(sense::MOI.OptimizationSense) = + sense == MOI.MAX_SENSE ? -Inf : Inf +_is_better(sense::MOI.OptimizationSense, new, best) = + sense == MOI.MAX_SENSE ? new > best : new < best +_gap(sense::MOI.OptimizationSense, best, bound) = + sense == MOI.MAX_SENSE ? bound - best : best - bound + +# one record per NLP solve, for convergence traces +function _log_progress( + model::Optimizer, + t_start::Float64, + best_objective, + master_bound + ) + model.silent && return + bound = master_bound === nothing ? NaN : master_bound + @info "LOA progress: elapsed=$(time() - t_start) " * + "incumbent=$best_objective bound=$bound" + return +end + +# one target per (binary, value); linear disjuncts need no cover +# since the master already carries them exactly +function _cover_disjuncts(problem::_Problem) + seen = Set{Tuple{MOI.VariableIndex, Bool}}() + cover = _Disjunct[] + for disjunction in problem.disjunctions, disjunct in disjunction.disjuncts + _is_nonlinear_disjunct(disjunct) || continue + key = (disjunct.binary, disjunct.active_value) + key in seen && continue + push!(seen, key) + push!(cover, disjunct) + end + return cover +end + +# GDPopt's covering weights: an uncovered disjunct outweighs all +# covered ones +function _cover_objective( + master::_Master, + cover::Vector{_Disjunct}, + needs_cover, + num_covered::Int + ) + objective = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0) + for i in eachindex(cover) + weight = Float64(needs_cover[i] ? num_covered + 1 : 1) + activation = _map_to(master.variable_map, cover[i].activation) + objective = MOI.Utilities.operate(+, Float64, objective, + MOI.Utilities.operate(*, Float64, weight, activation)) + end + return objective +end + +# variable starts warm start the first NLP; later iterations use +# the last feasible primal +function _user_start_values(model::Optimizer, problem::_Problem) + point = Dict{MOI.VariableIndex, Float64}() + for vi in problem.variables + start = MOI.get(model.cache, MOI.VariablePrimalStart(), vi) + start === nothing || (point[vi] = start) + end + return isempty(point) ? nothing : (point = point,) +end + +function _set_master_objective(master::_Master, sense, objective) + MOI.set(master.model, MOI.ObjectiveSense(), sense) + MOI.set(master.model, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + objective) + return +end + +################################################################################ +# MAIN LOOP +################################################################################ +function MOI.optimize!(model::Optimizer) + t_start = time() + _reset_results(model) + problem = _build_problem(model) + master = _build_master(model, problem) + subproblem = _build_subproblem(model, problem) + linearizer = _Linearizer() + sense = problem.sense + overall_deadline = t_start + Float64(_option(model, "time_limit")) + loop_deadline = min(overall_deadline, + t_start + Float64(_option(model, "iteration_time_limit"))) + + best_objective = _worst_objective(sense) + best_result = nothing + previous_result = _user_start_values(model, problem) + master_bound = nothing + master_status = nothing + converged = false + + # Shared iteration tail: no-good cut, OA cuts, incumbent update. + process_result = result -> begin + _avoid_combination(master, result.combination) + _add_oa_cuts(model, problem, master, linearizer, result) + if result.feasible && + _is_better(sense, result.objective, best_objective) + best_objective = result.objective + best_result = result + end + result.feasible && (previous_result = result) + _log_progress(model, t_start, best_objective, master_bound) + return + end + warm_start = () -> + previous_result === nothing ? nothing : previous_result.point + + # set covering: reuse the master with a coverage objective so + # every nonlinear disjunct gets visited once + cover = _cover_disjuncts(problem) + needs_cover = trues(length(cover)) + num_covered = 0 + for iteration in 1:_option(model, "set_cover_max_iter") + (iteration == 1 || any(needs_cover)) || break + time() < loop_deadline || break + _set_master_objective(master, MOI.MAX_SENSE, + _cover_objective(master, cover, needs_cover, num_covered)) + _cap_remaining_time(master.model, loop_deadline) + MOI.optimize!(master.model) + solved = _solved_and_feasible(master.model) + # capture the status before the objective restore invalidates it + status = MOI.get(master.model, MOI.TerminationStatus()) + combination = solved ? _extract_combination(problem, master) : nothing + _set_master_objective(master, master.sense, master.oa_objective) + if !solved + master_status = status + break + end + result = _solve_nlp(model, problem, subproblem, combination, + warm_start(); deadline = loop_deadline) + process_result(result) + # covered only once active in a feasible NLP; infeasible + # combinations just leave their no-good cut + if result.feasible + for i in eachindex(cover) + needs_cover[i] || continue + _disjunct_active(result.combination, cover[i]) && + (needs_cover[i] = false) + end + num_covered = count(!, needs_cover) + end + end + + # main loop: alpha_oa gives the bound, the NLP the incumbent + if master_status === nothing + for _ in 1:_option(model, "max_iter") + time() < loop_deadline || break + _cap_remaining_time(master.model, loop_deadline) + MOI.optimize!(master.model) + if !_solved_and_feasible(master.model) + master_status = MOI.get(master.model, MOI.TerminationStatus()) + break + end + master_bound = MOI.get(master.model, MOI.VariablePrimal(), + master.alpha_oa) + if best_result !== nothing + gap = _gap(sense, best_objective, master_bound) + total_slack = abs(MOI.get(master.model, + MOI.ObjectiveValue()) - master_bound) / + Float64(_option(model, "oa_penalty")) + tol = Float64(_option(model, "convergence_tol")) * + max(abs(best_objective), 1.0) + if gap <= tol && + total_slack <= Float64(_option(model, "slack_tol")) + converged = true + break + end + end + combination = _extract_combination(problem, master) + result = _solve_nlp(model, problem, subproblem, combination, + warm_start(); deadline = loop_deadline) + process_result(result) + end + end + + _store_results(model, sense, best_objective, best_result, master_bound, + master_status, converged, loop_deadline) + model.solve_time = time() - t_start + return +end + +################################################################################ +# RESULT SYNTHESIS +################################################################################ +# the OA bound is valid only for convex problems, so report +# LOCALLY_SOLVED, never OPTIMAL +function _store_results( + model::Optimizer, + sense::MOI.OptimizationSense, + best_objective::Float64, + best_result, + master_bound, + master_status, + converged::Bool, + loop_deadline::Float64 + ) + timed_out = time() >= loop_deadline + if best_result === nothing + model.primal_status = MOI.NO_SOLUTION + model.objective_value = NaN + model.raw_status = "No feasible incumbent found." + if master_status == MOI.INFEASIBLE + model.termination_status = MOI.INFEASIBLE + model.raw_status = "No feasible incumbent: the master " * + "problem is infeasible." + elseif timed_out + model.termination_status = MOI.TIME_LIMIT + elseif master_status === nothing + model.termination_status = MOI.ITERATION_LIMIT + else + model.termination_status = MOI.OTHER_LIMIT + model.raw_status = "No feasible incumbent: the master " * + "solve finished with status $master_status." + end + return + end + model.primal_status = MOI.FEASIBLE_POINT + model.incumbent = best_result.point + model.objective_value = best_objective + model.objective_bound = master_bound === nothing ? nothing : + Float64(master_bound) + if converged + model.termination_status = MOI.LOCALLY_SOLVED + elseif timed_out + model.termination_status = MOI.TIME_LIMIT + elseif master_status == MOI.INFEASIBLE + # all combinations visited; the incumbent is best over all + # of them, but the OA bound is gone + model.termination_status = MOI.LOCALLY_SOLVED + elseif master_status !== nothing + model.termination_status = MOI.OTHER_LIMIT + else + model.termination_status = MOI.ITERATION_LIMIT + end + label = converged ? "converged" : (master_status == MOI.INFEASIBLE ? + "combinations exhausted" : "limit hit") + if master_bound === nothing + model.raw_status = "LOA finished [$label]: incumbent " * + "$best_objective (master produced no bound)." + else + gap = _gap(sense, best_objective, master_bound) + relative = abs(best_objective) > 1e-10 ? + gap / abs(best_objective) : gap + model.relative_gap = relative + model.raw_status = "LOA finished [$label]: incumbent " * + "$best_objective, master bound $master_bound, gap $gap " * + "(relative $relative)." + end + return +end diff --git a/src/master.jl b/src/master.jl new file mode 100644 index 0000000..7e35d78 --- /dev/null +++ b/src/master.jl @@ -0,0 +1,181 @@ +################################################################################ +# MASTER CONSTRUCTION +################################################################################ +# `alpha_oa` carries the objective; `oa_objective` also tracks the +# slack penalties so set covering can swap objectives out and back +mutable struct _Master + model::MOI.ModelLike + variable_map::Dict{MOI.VariableIndex, MOI.VariableIndex} + sense::MOI.OptimizationSense + objective::MOI.AbstractScalarFunction + alpha_oa::MOI.VariableIndex + oa_objective::MOI.ScalarAffineFunction{Float64} +end + +function _instantiate(factory, name::String) + factory === nothing && + error("DisjunctiveAlgorithms requires the `$name` optimizer factory.") + solver = MOI.instantiate(factory; + with_cache_type = Float64, with_bridge_type = Float64) + MOI.set(solver, MOI.Silent(), true) + return solver +end + +_map_to(variable_map::AbstractDict, func) = + MOI.Utilities.map_indices(vi -> variable_map[vi], func) + +function _solved_and_feasible(solver::MOI.ModelLike) + status = MOI.get(solver, MOI.TerminationStatus()) + return status in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) && + MOI.get(solver, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT +end + +function _cap_remaining_time(solver::MOI.ModelLike, deadline::Float64) + isfinite(deadline) || return + MOI.set(solver, MOI.TimeLimitSec(), max(0.0, deadline - time())) + return +end + +function _build_master(model::Optimizer, problem::_Problem) + mip = _instantiate(model.mip_solver, "mip_solver") + variable_map = Dict{MOI.VariableIndex, MOI.VariableIndex}( + vi => MOI.add_variable(mip) for vi in problem.variables) + for ci in problem.variable_cis + vi = MOI.get(model.cache, MOI.ConstraintFunction(), ci) + MOI.add_constraint(mip, variable_map[vi], + MOI.get(model.cache, MOI.ConstraintSet(), ci)) + end + for ci in problem.linear_cis + func = MOI.get(model.cache, MOI.ConstraintFunction(), ci) + MOI.add_constraint(mip, _map_to(variable_map, func), + MOI.get(model.cache, MOI.ConstraintSet(), ci)) + end + # nonlinear-typed rows that demoted to affine stay in the master + for (func, set) in problem.nonlinear_rows + _is_linear(func) || continue + MOI.add_constraint(mip, _map_to(variable_map, _to_affine(func)), set) + end + for disjunction in problem.disjunctions + _add_exactly_one(mip, variable_map, disjunction) + for disjunct in disjunction.disjuncts + _add_gated_rows(mip, variable_map, disjunct, model.options) + end + end + alpha_oa = MOI.add_variable(mip) + sense = problem.sense + oa_objective = MOI.ScalarAffineFunction( + [MOI.ScalarAffineTerm(1.0, alpha_oa)], 0.0) + MOI.set(mip, MOI.ObjectiveSense(), sense) + MOI.set(mip, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + oa_objective) + return _Master(mip, variable_map, sense, problem.objective, alpha_oa, + oa_objective) +end + +# Indicators sum to the activation (1 top-level, parent indicator +# nested); a complement pair normalizes to a trivial `1 == 1` row +function _add_exactly_one( + mip::MOI.ModelLike, + variable_map::AbstractDict, + disjunction::_Disjunction + ) + total = MOI.Utilities.operate(+, Float64, + (_map_to(variable_map, disjunct.activation) + for disjunct in disjunction.disjuncts)...) + body = MOI.Utilities.operate(-, Float64, total, + _map_to(variable_map, disjunction.activation)) + MOI.Utilities.normalize_and_add_constraint(mip, body, MOI.EqualTo(0.0)) + return +end + +# split Interval rows; the indicator bridge takes one-sided sets only +_indicator_sets(set::MOI.Interval{Float64}) = + (MOI.LessThan(set.upper), MOI.GreaterThan(set.lower)) +_indicator_sets(set::MOI.AbstractScalarSet) = (set,) + +# Gate each linear row with an indicator constraint or big-M per +# `master_gating`; nonlinear rows enter the master only as OA cuts +function _add_gated_rows( + mip::MOI.ModelLike, + variable_map::AbstractDict, + disjunct::_Disjunct, + options::Dict{String, Any} + ) + gating = options["master_gating"] + gating in ("indicator", "bigm") || + error("Unknown `master_gating` value `$gating`.") + activate = disjunct.active_value ? MOI.ACTIVATE_ON_ONE : + MOI.ACTIVATE_ON_ZERO + binary = variable_map[disjunct.binary] + activation = _map_to(variable_map, disjunct.activation) + M = Float64(options["M_value"]) + # M * (1 - activation) moved left, as in the disjunct OA cuts + gate = MOI.Utilities.operate(-, Float64, + MOI.Utilities.operate(*, Float64, M, activation), M) + for (func, set) in zip(disjunct.functions, disjunct.sets) + _is_linear(func) || continue + row = _map_to(variable_map, _to_affine(func)) + if gating == "bigm" + for term in _oa_cut_terms(set, row) + body = MOI.Utilities.operate(+, Float64, term, gate) + MOI.Utilities.normalize_and_add_constraint(mip, body, + MOI.LessThan(0.0)) + end + else + for inner in _indicator_sets(set) + gated = MOI.Utilities.operate(vcat, Float64, binary, row) + MOI.add_constraint(mip, gated, + MOI.Indicator{activate}(inner)) + end + end + end + return +end + +# bounded penalized slack so an invalid cut cannot blow up the master +function _add_penalized_slack( + master::_Master, + options::Dict{String, Any}, + penalty_sign::Int + ) + slack = MOI.add_variable(master.model) + MOI.add_constraint(master.model, slack, MOI.GreaterThan(0.0)) + MOI.add_constraint(master.model, slack, + MOI.LessThan(Float64(options["max_slack"]))) + penalty = penalty_sign * Float64(options["oa_penalty"]) + master.oa_objective = MOI.Utilities.operate(+, Float64, + master.oa_objective, MOI.ScalarAffineFunction( + [MOI.ScalarAffineTerm(penalty, slack)], 0.0)) + MOI.set(master.model, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + master.oa_objective) + return slack +end + +# round to Bool; MILP values are only within integer tolerance +function _extract_combination(problem::_Problem, master::_Master) + return Dict{MOI.VariableIndex, Bool}( + binary => round(Bool, MOI.get(master.model, MOI.VariablePrimal(), + master.variable_map[binary])) + for binary in problem.binaries) +end + +# No-good cut: active `1 - z` plus inactive `z` terms must reach 1 +function _avoid_combination(master::_Master, combination::AbstractDict) + terms = MOI.ScalarAffineTerm{Float64}[] + constant = 0.0 + for (binary, value) in combination + mapped = master.variable_map[binary] + if value + push!(terms, MOI.ScalarAffineTerm(-1.0, mapped)) + constant += 1.0 + else + push!(terms, MOI.ScalarAffineTerm(1.0, mapped)) + end + end + cut = MOI.ScalarAffineFunction(terms, constant) + MOI.Utilities.normalize_and_add_constraint(master.model, cut, + MOI.GreaterThan(1.0)) + return +end diff --git a/src/nlp.jl b/src/nlp.jl new file mode 100644 index 0000000..eaff857 --- /dev/null +++ b/src/nlp.jl @@ -0,0 +1,168 @@ +################################################################################ +# NLP SUBPROBLEM +################################################################################ +# Built once; each iteration overwrites the binary fixes in place +# and swaps the active disjuncts' rows. No big-M anywhere. +struct _Subproblem + model::MOI.ModelLike + variable_map::Dict{MOI.VariableIndex, MOI.VariableIndex} + fixes::Dict{MOI.VariableIndex, + MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{Float64}}} + rows::Vector{MOI.ConstraintIndex} +end + +function _build_subproblem(model::Optimizer, problem::_Problem) + nlp = _instantiate(model.nlp_solver, "nlp_solver") + variable_map = Dict{MOI.VariableIndex, MOI.VariableIndex}( + vi => MOI.add_variable(nlp) for vi in problem.variables) + indicators = Set(problem.binaries) + for ci in problem.variable_cis + vi = MOI.get(model.cache, MOI.ConstraintFunction(), ci) + vi in indicators && continue + MOI.add_constraint(nlp, variable_map[vi], + MOI.get(model.cache, MOI.ConstraintSet(), ci)) + end + fixes = Dict(binary => MOI.add_constraint(nlp, variable_map[binary], + MOI.EqualTo(0.0)) for binary in problem.binaries) + for ci in problem.linear_cis + func = MOI.get(model.cache, MOI.ConstraintFunction(), ci) + MOI.add_constraint(nlp, _map_to(variable_map, func), + MOI.get(model.cache, MOI.ConstraintSet(), ci)) + end + for (func, set) in problem.nonlinear_rows + MOI.add_constraint(nlp, _map_to(variable_map, func), set) + end + MOI.set(nlp, MOI.ObjectiveSense(), problem.sense) + objective = _map_to(variable_map, problem.objective) + MOI.set(nlp, MOI.ObjectiveFunction{typeof(objective)}(), objective) + return _Subproblem(nlp, variable_map, fixes, MOI.ConstraintIndex[]) +end + +function _set_warm_start(nlp::MOI.ModelLike, variable_map::AbstractDict, point) + point === nothing && return + for (vi, value) in point + MOI.set(nlp, MOI.VariablePrimalStart(), variable_map[vi], value) + end + return +end + +function _extract_point( + nlp::MOI.ModelLike, + problem::_Problem, + variable_map::AbstractDict + ) + return Dict{MOI.VariableIndex, Float64}( + vi => MOI.get(nlp, MOI.VariablePrimal(), variable_map[vi]) + for vi in problem.variables) +end + +# Solve the NLP at a fixed combination: overwrite the binary fixes, +# swap the active disjuncts' rows, and optimize. If infeasible, fall +# through to NLPF (a slacked version that always solves) so the master +# still gets a linearization site, not just a no-good cut. +function _solve_nlp( + model::Optimizer, + problem::_Problem, + sub::_Subproblem, + combination::AbstractDict, + warm_start; + deadline::Float64 = Inf + ) + for (binary, value) in combination + MOI.set(sub.model, MOI.ConstraintSet(), sub.fixes[binary], + MOI.EqualTo(value ? 1.0 : 0.0)) + end + for ci in sub.rows + MOI.delete(sub.model, ci) + end + empty!(sub.rows) + for disjunction in problem.disjunctions, disjunct in disjunction.disjuncts + _disjunct_active(combination, disjunct) || continue + for (func, set) in zip(disjunct.functions, disjunct.sets) + push!(sub.rows, MOI.add_constraint(sub.model, + _map_to(sub.variable_map, func), set)) + end + end + _set_warm_start(sub.model, sub.variable_map, warm_start) + _cap_remaining_time(sub.model, deadline) + MOI.optimize!(sub.model) + if _solved_and_feasible(sub.model) + return (combination = combination, + point = _extract_point(sub.model, problem, sub.variable_map), + objective = MOI.get(sub.model, MOI.ObjectiveValue()), + feasible = true) + end + if Bool(_option(model, "use_nlpf")) + result = _solve_nlpf(model, problem, combination, warm_start; + deadline = deadline) + result === nothing || return result + end + return (combination = combination, + point = nothing, objective = Inf, feasible = false) +end + +################################################################################ +# NLPF (FEASIBILITY SUBPROBLEM) +################################################################################ +_nlpf_slacked(func, u, ::MOI.LessThan{Float64}) = + MOI.Utilities.operate(-, Float64, func, u) +_nlpf_slacked(func, u, ::MOI.GreaterThan{Float64}) = + MOI.Utilities.operate(+, Float64, func, u) +_nlpf_slacked(func, u, ::MOI.AbstractScalarSet) = nothing + +# The slacked feasibility NLP: one nonnegative `u` relaxes every scalar +# inequality row (bounds and equalities stay exact) and is minimized. +# Its solution is a linearization site for an infeasible combination. +function _solve_nlpf( + model::Optimizer, + problem::_Problem, + combination::AbstractDict, + warm_start; + deadline::Float64 = Inf + ) + nlp = _instantiate(model.nlp_solver, "nlp_solver") + variable_map = Dict{MOI.VariableIndex, MOI.VariableIndex}( + vi => MOI.add_variable(nlp) for vi in problem.variables) + u = MOI.add_variable(nlp) + MOI.add_constraint(nlp, u, MOI.GreaterThan(0.0)) + for ci in problem.variable_cis + vi = MOI.get(model.cache, MOI.ConstraintFunction(), ci) + haskey(combination, vi) && continue + MOI.add_constraint(nlp, variable_map[vi], + MOI.get(model.cache, MOI.ConstraintSet(), ci)) + end + for (binary, value) in combination + MOI.add_constraint(nlp, variable_map[binary], + MOI.EqualTo(value ? 1.0 : 0.0)) + end + rows = Tuple{MOI.AbstractScalarFunction, MOI.AbstractScalarSet}[] + for ci in problem.linear_cis + push!(rows, (MOI.get(model.cache, MOI.ConstraintFunction(), ci), + MOI.get(model.cache, MOI.ConstraintSet(), ci))) + end + append!(rows, problem.nonlinear_rows) + for disjunction in problem.disjunctions, disjunct in disjunction.disjuncts + _disjunct_active(combination, disjunct) || continue + append!(rows, zip(disjunct.functions, disjunct.sets)) + end + for (func, set) in rows + mapped = _map_to(variable_map, func) + slacked = _nlpf_slacked(mapped, u, set) + if slacked === nothing + MOI.add_constraint(nlp, mapped, set) + else + MOI.add_constraint(nlp, slacked, set) + end + end + MOI.set(nlp, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.set(nlp, MOI.ObjectiveFunction{MOI.VariableIndex}(), u) + _set_warm_start(nlp, variable_map, warm_start) + _cap_remaining_time(nlp, deadline) + MOI.optimize!(nlp) + # Use the primal only at a genuine feasible point; a solver can + # report values at a nonfeasible/NaN primal that poisons the cut. + _solved_and_feasible(nlp) || return nothing + return (combination = combination, + point = _extract_point(nlp, problem, variable_map), + objective = Inf, feasible = false) +end diff --git a/src/optimizer.jl b/src/optimizer.jl new file mode 100644 index 0000000..b55dd60 --- /dev/null +++ b/src/optimizer.jl @@ -0,0 +1,380 @@ +################################################################################ +# OPTIMIZER +################################################################################ +const _Cache = MOI.Utilities.UniversalFallback{MOI.Utilities.Model{Float64}} + +# Raw options, mirroring DisjunctiveProgramming.jl's LOA defaults. +const _DEFAULT_OPTIONS = Dict{String, Any}( + "max_iter" => 10, + "set_cover_max_iter" => 8, + "M_value" => 1e9, + "master_gating" => "indicator", + "max_slack" => 1e3, + "oa_penalty" => 1e3, + "use_nlpf" => true, + "convergence_tol" => 1e-6, + "slack_tol" => 1e-4, + "iteration_time_limit" => Inf, + "time_limit" => 3600.0, +) + +""" + Optimizer(; nlp_solver, mip_solver = nlp_solver, kwargs...) + +Logic-based outer approximation solver for models containing +[`DisjunctionSet`](@ref) constraints. `nlp_solver` and `mip_solver` +are optimizer factories as accepted by `MOI.instantiate`. The +remaining keyword arguments set raw options (also reachable through +`MOI.RawOptimizerAttribute`): + +- `max_iter = 10`: master/NLP iterations after the set-covering seed. +- `set_cover_max_iter = 8`: set-covering initialization iterations. +- `M_value = 1e9`: big-M gating the disjunct OA cuts in the master. +- `master_gating = "indicator"`: how linear disjunct rows enter the + master, `"indicator"` (constraint gated by the binary) or `"bigm"` + (rows relaxed by `M_value * (1 - z)`, tighter for solvers that + cannot strengthen indicators, e.g. with presolve disabled). +- `max_slack = 1e3`: upper bound of each OA cut slack. +- `oa_penalty = 1e3`: objective penalty per unit of cut slack. +- `use_nlpf = true`: solve a slacked feasibility NLP when the primary + NLP is infeasible, so its point still seeds OA cuts. +- `convergence_tol = 1e-6`: relative incumbent/bound gap tolerance. +- `slack_tol = 1e-4`: total cut slack tolerance for convergence. +- `iteration_time_limit = Inf`: seconds allotted to the LOA loop. +- `time_limit = 3600.0`: overall seconds budget. +""" +mutable struct Optimizer <: MOI.AbstractOptimizer + nlp_solver::Any + mip_solver::Any + cache::_Cache + options::Dict{String, Any} + silent::Bool + # results (filled by MOI.optimize!) + termination_status::MOI.TerminationStatusCode + primal_status::MOI.ResultStatusCode + incumbent::Dict{MOI.VariableIndex, Float64} + objective_value::Float64 + objective_bound::Union{Nothing, Float64} + relative_gap::Float64 + raw_status::String + solve_time::Float64 +end + +function Optimizer(; + nlp_solver = nothing, + mip_solver = nlp_solver, + kwargs... + ) + options = copy(_DEFAULT_OPTIONS) + for (key, value) in kwargs + haskey(options, string(key)) || throw(ArgumentError( + "Unknown option `$key`.")) + options[string(key)] = value + end + return Optimizer(nlp_solver, mip_solver, + MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), + options, false, MOI.OPTIMIZE_NOT_CALLED, MOI.NO_SOLUTION, + Dict{MOI.VariableIndex, Float64}(), NaN, nothing, NaN, "", NaN) +end + +_option(model::Optimizer, name::String) = model.options[name] + +MOI.get(::Optimizer, ::MOI.SolverName) = "DisjunctiveAlgorithms" +MOI.get(::Optimizer, ::MOI.SolverVersion) = "0.1.0" + +MOI.is_empty(model::Optimizer) = MOI.is_empty(model.cache) + +function MOI.empty!(model::Optimizer) + MOI.empty!(model.cache) + _reset_results(model) + return +end + +# The non-cache part of MOI.empty!, also run at the top of every solve +# so a re-optimize cannot leak the previous solve's bound, gap, or +# point. +function _reset_results(model::Optimizer) + model.termination_status = MOI.OPTIMIZE_NOT_CALLED + model.primal_status = MOI.NO_SOLUTION + empty!(model.incumbent) + model.objective_value = NaN + model.objective_bound = nothing + model.relative_gap = NaN + model.raw_status = "" + model.solve_time = NaN + return +end + +MOI.supports_incremental_interface(::Optimizer) = true + +function MOI.copy_to(model::Optimizer, src::MOI.ModelLike) + return MOI.Utilities.default_copy_to(model, src) +end + +################################################################################ +# MODEL-BUILDING FORWARDING +################################################################################ +# Everything model-building lands in the cache; the solve partitions it +# in MOI.optimize!, so incremental additions need no invalidation. +const _Index = Union{MOI.VariableIndex, MOI.ConstraintIndex} + +MOI.add_variable(model::Optimizer) = MOI.add_variable(model.cache) + +MOI.add_variables(model::Optimizer, n::Int) = MOI.add_variables(model.cache, n) + +function MOI.add_constrained_variable( + model::Optimizer, + set::MOI.AbstractScalarSet + ) + return MOI.add_constrained_variable(model.cache, set) +end + +function MOI.add_constrained_variables( + model::Optimizer, + set::MOI.AbstractVectorSet + ) + return MOI.add_constrained_variables(model.cache, set) +end + +function MOI.add_constraint( + model::Optimizer, + func::MOI.AbstractFunction, + set::MOI.AbstractSet + ) + return MOI.add_constraint(model.cache, func, set) +end + +# Honest constraint support (the cache would claim everything, which +# disables the bridges that rewrite e.g. vector cones into supported +# scalar rows): only what `_build_problem` actually partitions. +const _ScalarFunction = Union{MOI.ScalarAffineFunction{Float64}, + MOI.ScalarQuadraticFunction{Float64}, MOI.ScalarNonlinearFunction} +const _VectorFunction = Union{MOI.VectorOfVariables, + MOI.VectorAffineFunction{Float64}, + MOI.VectorQuadraticFunction{Float64}, MOI.VectorNonlinearFunction} + +function MOI.supports_constraint( + ::Optimizer, + ::Type{<:MOI.AbstractFunction}, + ::Type{<:MOI.AbstractSet} + ) + return false +end + +# Non-indicator discrete variables pass through to the subproblems +# with their integrality intact, so the nlp_solver must handle them. +function MOI.supports_constraint( + ::Optimizer, + ::Type{MOI.VariableIndex}, + ::Type{<:Union{_SupportedInnerSet, MOI.ZeroOne, MOI.Integer}} + ) + return true +end + +function MOI.supports_constraint( + ::Optimizer, + ::Type{<:_ScalarFunction}, + ::Type{<:_SupportedInnerSet} + ) + return true +end + +function MOI.supports_constraint( + ::Optimizer, + ::Type{<:_VectorFunction}, + ::Type{DisjunctionSet} + ) + return true +end + +# A `DisjunctionSet` never constrains variables on creation: without +# this, `copy_to` turns a `VectorOfVariables` disjunction (which may +# repeat a variable across rows) into `add_constrained_variables`. +function MOI.supports_add_constrained_variables( + ::Optimizer, + ::Type{DisjunctionSet} + ) + return false +end + +MOI.is_valid(model::Optimizer, index::_Index) = + MOI.is_valid(model.cache, index) + +MOI.delete(model::Optimizer, index::_Index) = MOI.delete(model.cache, index) + +function MOI.get(model::Optimizer, T::Type{<:_Index}, name::String) + return MOI.get(model.cache, T, name) +end + +# Non-result attributes forward to the cache; the result attributes +# below override these generic methods, and any solve-set attribute +# without an override (duals, ConstraintPrimal) is refused rather than +# forwarded to the cache, which never holds results. +# Forward to the wrapped inner model (honest), not the UniversalFallback, +# which claims support for every attribute and so never rejects an +# unsupported one on copy_to. +MOI.supports(model::Optimizer, attr::MOI.AbstractModelAttribute) = + MOI.supports(model.cache.model, attr) + +function MOI.set(model::Optimizer, attr::MOI.AbstractModelAttribute, value) + MOI.supports(model, attr) || throw(MOI.UnsupportedAttribute(attr)) + return MOI.set(model.cache, attr, value) +end + +function MOI.get(model::Optimizer, attr::MOI.AbstractModelAttribute) + MOI.is_set_by_optimize(attr) && throw(MOI.GetAttributeNotAllowed(attr)) + return MOI.get(model.cache, attr) +end + +function MOI.supports( + model::Optimizer, + attr::MOI.AbstractVariableAttribute, + ::Type{MOI.VariableIndex} + ) + return MOI.supports(model.cache.model, attr, MOI.VariableIndex) +end + +# The LOA loop consumes warm starts, so accept them even though the +# honest inner model does not store them (the fallback does). +function MOI.supports( + ::Optimizer, + ::MOI.VariablePrimalStart, + ::Type{MOI.VariableIndex} + ) + return true +end + +function MOI.set( + model::Optimizer, + attr::MOI.AbstractVariableAttribute, + vi::MOI.VariableIndex, + value + ) + MOI.supports(model, attr, MOI.VariableIndex) || + throw(MOI.UnsupportedAttribute(attr)) + return MOI.set(model.cache, attr, vi, value) +end + +function MOI.get( + model::Optimizer, + attr::MOI.AbstractVariableAttribute, + vi::MOI.VariableIndex + ) + MOI.is_set_by_optimize(attr) && throw(MOI.GetAttributeNotAllowed(attr)) + return MOI.get(model.cache, attr, vi) +end + +function MOI.supports( + model::Optimizer, + attr::MOI.AbstractConstraintAttribute, + C::Type{<:MOI.ConstraintIndex} + ) + return MOI.supports(model.cache.model, attr, C) +end + +function MOI.set( + model::Optimizer, + attr::MOI.AbstractConstraintAttribute, + ci::MOI.ConstraintIndex, + value + ) + MOI.supports(model, attr, typeof(ci)) || + throw(MOI.UnsupportedAttribute(attr)) + return MOI.set(model.cache, attr, ci, value) +end + +function MOI.get( + model::Optimizer, + attr::MOI.AbstractConstraintAttribute, + ci::MOI.ConstraintIndex + ) + MOI.is_set_by_optimize(attr) && throw(MOI.GetAttributeNotAllowed(attr)) + return MOI.get(model.cache, attr, ci) +end + +################################################################################ +# OPTIMIZER ATTRIBUTES +################################################################################ +MOI.supports(::Optimizer, ::MOI.Silent) = true + +function MOI.set(model::Optimizer, ::MOI.Silent, value::Bool) + model.silent = value + return +end + +MOI.get(model::Optimizer, ::MOI.Silent) = model.silent + +MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true + +function MOI.set( + model::Optimizer, + ::MOI.TimeLimitSec, + value::Union{Nothing, Real} + ) + model.options["time_limit"] = value === nothing ? Inf : Float64(value) + return +end + +function MOI.get(model::Optimizer, ::MOI.TimeLimitSec) + limit = model.options["time_limit"] + return isfinite(limit) ? limit : nothing +end + +function MOI.supports(model::Optimizer, attr::MOI.RawOptimizerAttribute) + return haskey(model.options, attr.name) +end + +function MOI.set(model::Optimizer, attr::MOI.RawOptimizerAttribute, value) + MOI.supports(model, attr) || throw(MOI.UnsupportedAttribute(attr)) + model.options[attr.name] = value + return +end + +function MOI.get(model::Optimizer, attr::MOI.RawOptimizerAttribute) + MOI.supports(model, attr) || throw(MOI.UnsupportedAttribute(attr)) + return model.options[attr.name] +end + +################################################################################ +# RESULT ATTRIBUTES +################################################################################ +MOI.get(model::Optimizer, ::MOI.TerminationStatus) = model.termination_status + +MOI.get(model::Optimizer, ::MOI.RawStatusString) = model.raw_status + +function MOI.get(model::Optimizer, ::MOI.ResultCount) + return model.primal_status == MOI.NO_SOLUTION ? 0 : 1 +end + +function MOI.get(model::Optimizer, attr::MOI.PrimalStatus) + return attr.result_index == 1 ? model.primal_status : MOI.NO_SOLUTION +end + +MOI.get(model::Optimizer, ::MOI.DualStatus) = MOI.NO_SOLUTION + +function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue) + MOI.check_result_index_bounds(model, attr) + return model.objective_value +end + +function MOI.get(model::Optimizer, ::MOI.ObjectiveBound) + bound = model.objective_bound + if bound === nothing + sense = MOI.get(model.cache, MOI.ObjectiveSense()) + return sense == MOI.MAX_SENSE ? Inf : -Inf + end + return bound +end + +function MOI.get( + model::Optimizer, + attr::MOI.VariablePrimal, + vi::MOI.VariableIndex + ) + MOI.check_result_index_bounds(model, attr) + return model.incumbent[vi] +end + +MOI.get(model::Optimizer, ::MOI.RelativeGap) = model.relative_gap + +MOI.get(model::Optimizer, ::MOI.SolveTimeSec) = model.solve_time diff --git a/src/problem.jl b/src/problem.jl new file mode 100644 index 0000000..ebce607 --- /dev/null +++ b/src/problem.jl @@ -0,0 +1,159 @@ +################################################################################ +# PROBLEM PARTITION +################################################################################ +# activation is `z` or `1 - z`; rows stay in cache space +struct _Disjunct + activation::MOI.ScalarAffineFunction{Float64} + binary::MOI.VariableIndex + active_value::Bool + functions::Vector{MOI.AbstractScalarFunction} + sets::Vector{MOI.AbstractScalarSet} +end + +# activation: 1 at top level, the parent indicator when nested +struct _Disjunction + activation::MOI.ScalarAffineFunction{Float64} + disjuncts::Vector{_Disjunct} +end + +# Constraint indices stay in cache space; the builders remap them +struct _Problem + variables::Vector{MOI.VariableIndex} + binaries::Vector{MOI.VariableIndex} + disjunctions::Vector{_Disjunction} + variable_cis::Vector{MOI.ConstraintIndex} + linear_cis::Vector{MOI.ConstraintIndex} + # stable function objects so the linearizer can cache evaluators + nonlinear_rows::Vector{Tuple{MOI.AbstractScalarFunction, + MOI.AbstractScalarSet}} + sense::MOI.OptimizationSense + objective::MOI.AbstractScalarFunction +end + +_is_linear(::Union{MOI.VariableIndex, MOI.ScalarAffineFunction}) = true +_is_linear(::MOI.AbstractScalarFunction) = false + +_to_affine(func::MOI.ScalarAffineFunction{Float64}) = func +function _to_affine(func::MOI.AbstractScalarFunction) + return convert(MOI.ScalarAffineFunction{Float64}, func) +end + +# a raw `VariableIndex` row becomes a bound and collides with the +# variable's own bounds in the subproblem +_as_row(func::MOI.VariableIndex) = _to_affine(func) +_as_row(func::MOI.AbstractScalarFunction) = func + +# demote rows the enclosing vector function promoted +_demote(func::MOI.AbstractScalarFunction) = func +function _demote(func::MOI.ScalarQuadraticFunction{Float64}) + return _try_convert(MOI.ScalarAffineFunction{Float64}, func) +end +function _demote(func::MOI.ScalarNonlinearFunction) + return _try_convert(MOI.ScalarAffineFunction{Float64}, + _try_convert(MOI.ScalarQuadraticFunction{Float64}, func)) +end + +function _try_convert(T::Type, func) + return try + convert(T, func) + catch + func + end +end + +_scalarize(func::MOI.AbstractVectorFunction) = + collect(MOI.Utilities.eachscalar(func)) + +# `z` -> (z, true), `1 - z` -> (z, false) +function _activation_binary(activation::MOI.ScalarAffineFunction{Float64}) + canonical = MOI.Utilities.canonical(activation) + if length(canonical.terms) == 1 + term = only(canonical.terms) + if term.coefficient == 1.0 && canonical.constant == 0.0 + return term.variable, true + elseif term.coefficient == -1.0 && canonical.constant == 1.0 + return term.variable, false + end + end + return error("Unsupported indicator expression `$activation`: each " * + "`DisjunctionSet` indicator must be a binary variable `z` or " * + "its complement `1 - z`.") +end + +function _parse_disjunction( + cache::_Cache, + ci::MOI.ConstraintIndex{F, DisjunctionSet} + ) where {F} + set = MOI.get(cache, MOI.ConstraintSet(), ci) + rows = _scalarize(MOI.get(cache, MOI.ConstraintFunction(), ci)) + disjunction_activation = _to_affine( + _demote(rows[activation_index(set)])) + disjuncts = _Disjunct[] + for (i, j) in enumerate(indicator_indices(set)) + activation = _to_affine(_demote(rows[j])) + binary, active_value = _activation_binary(activation) + zero_one = MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}( + binary.value) + MOI.is_valid(cache, zero_one) || error("The indicator variable " * + "of a `DisjunctionSet` disjunct must be `MOI.ZeroOne`.") + functions = MOI.AbstractScalarFunction[ + _as_row(_demote(rows[k])) for k in row_indices(set, i)] + push!(disjuncts, _Disjunct(activation, binary, active_value, + functions, set.inner_sets[i])) + end + return _Disjunction(disjunction_activation, disjuncts) +end + +# partition the cache for the LOA loop +function _build_problem(model::Optimizer) + cache = model.cache + sense = MOI.get(cache, MOI.ObjectiveSense()) + if sense == MOI.FEASIBILITY_SENSE + # feasibility model: minimize a constant zero objective + sense = MOI.MIN_SENSE + objective = MOI.ScalarAffineFunction( + MOI.ScalarAffineTerm{Float64}[], 0.0) + else + F = MOI.get(cache, MOI.ObjectiveFunctionType()) + objective = _demote(MOI.get(cache, MOI.ObjectiveFunction{F}())) + end + variables = MOI.get(cache, MOI.ListOfVariableIndices()) + disjunctions = _Disjunction[] + variable_cis = MOI.ConstraintIndex[] + linear_cis = MOI.ConstraintIndex[] + nonlinear_rows = Tuple{MOI.AbstractScalarFunction, + MOI.AbstractScalarSet}[] + for (FC, S) in MOI.get(cache, MOI.ListOfConstraintTypesPresent()) + cis = MOI.get(cache, MOI.ListOfConstraintIndices{FC, S}()) + if S === DisjunctionSet + append!(disjunctions, + _parse_disjunction(cache, ci) for ci in cis) + elseif FC === MOI.VariableIndex && S <: MOI.AbstractScalarSet + append!(variable_cis, cis) + elseif FC === MOI.ScalarAffineFunction{Float64} && + S <: MOI.AbstractScalarSet + append!(linear_cis, cis) + elseif FC <: Union{MOI.ScalarQuadraticFunction{Float64}, + MOI.ScalarNonlinearFunction} && S <: MOI.AbstractScalarSet + append!(nonlinear_rows, + (_demote(MOI.get(cache, MOI.ConstraintFunction(), ci)), + MOI.get(cache, MOI.ConstraintSet(), ci)) + for ci in cis) + else + error("DisjunctiveAlgorithms does not support `$FC`-in-`$S` constraints.") + end + end + binaries = unique!([disjunct.binary for disjunction in disjunctions + for disjunct in disjunction.disjuncts]) + # non-indicator discrete variables keep their integrality; the + # nlp_solver must handle whatever remains + return _Problem(variables, binaries, disjunctions, variable_cis, + linear_cis, nonlinear_rows, sense, objective) +end + +function _disjunct_active(combination::AbstractDict, disjunct::_Disjunct) + return combination[disjunct.binary] == disjunct.active_value +end + +_is_nonlinear_disjunct(disjunct::_Disjunct) = + any(!_is_linear(func) for func in disjunct.functions) diff --git a/test/integration.jl b/test/integration.jl new file mode 100644 index 0000000..41067e0 --- /dev/null +++ b/test/integration.jl @@ -0,0 +1,121 @@ +using DisjunctiveProgramming, HiGHS, Ipopt, InfiniteOpt + +function _optimizer_factory() + return () -> DA.Optimizer(nlp_solver = Ipopt.Optimizer, + mip_solver = HiGHS.Optimizer) +end + +# The same GDP solved through a BigM reformulation and through the +# lowering into DisjunctiveAlgorithms must agree: max x with +# (x <= 3) or (x <= 7). +function test_lowering_solve_linear() + model = GDPModel(_optimizer_factory()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, Y[1:2], Logical) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x <= 7, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, x) + optimize!(model, gdp_method = MOIDisjunction()) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 7.0 atol = 1e-4 + @test value(x) ≈ 7.0 atol = 1e-4 + @test value(Y[2]) + + reference = GDPModel(HiGHS.Optimizer) + set_silent(reference) + @variable(reference, 0 <= x2 <= 10) + @variable(reference, Y2[1:2], Logical) + @constraint(reference, x2 <= 3, Disjunct(Y2[1])) + @constraint(reference, x2 <= 7, Disjunct(Y2[2])) + @disjunction(reference, Y2) + @objective(reference, Max, x2) + optimize!(reference, gdp_method = BigM()) + @test objective_value(model) ≈ objective_value(reference) atol = 1e-6 +end + +# Nonlinear disjunct row: max x with (x <= 3) or (x^2 == 64); the +# unique optimum 8 is checked directly +function test_lowering_solve_nonlinear() + model = GDPModel(_optimizer_factory()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, Y[1:2], Logical) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x^2 == 64, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, x) + optimize!(model, gdp_method = MOIDisjunction()) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 8.0 atol = 1e-3 + @test value(Y[2]) +end + +# Nested solve: max x with x <= 2, or a nested mode choice between +# x^2 <= 25 and x <= 8. +function test_lowering_solve_nested() + model = GDPModel(_optimizer_factory()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, Y[1:2], Logical) + @variable(model, W[1:2], Logical) + @constraint(model, x <= 2, Disjunct(Y[1])) + @constraint(model, x^2 <= 25, Disjunct(W[1])) + @constraint(model, x <= 8, Disjunct(W[2])) + @disjunction(model, W, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, x) + optimize!(model, gdp_method = MOIDisjunction()) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 8.0 atol = 1e-3 + @test value(Y[2]) + @test value(W[2]) +end + +# An InfiniteGDPModel lowers through the same method: the disjunction +# transcribes to one DisjunctionSet constraint per support. +function test_lowering_infinite() + model = InfiniteGDPModel(_optimizer_factory()) + set_silent(model) + @infinite_parameter(model, t in [0, 1], num_supports = 3) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y[1:2], InfiniteLogical(t)) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x >= 5, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, integral(x, t)) + optimize!(model, gdp_method = MOIDisjunction()) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 10.0 atol = 1e-4 + @test all(value(x) .>= 5.0 .- 1e-4) +end + +# A nested infinite disjunction transcribes with the parent indicator +# as its activation, so the rebuilt constraint function is a uniform +# vector of scalar expressions (no constant to promote). +function test_lowering_infinite_nested() + model = InfiniteGDPModel(_optimizer_factory()) + set_silent(model) + @infinite_parameter(model, t in [0, 1], num_supports = 3) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y[1:2], InfiniteLogical(t)) + @variable(model, W[1:2], InfiniteLogical(t)) + @constraint(model, x <= 2, Disjunct(Y[1])) + @constraint(model, x >= 5, Disjunct(W[1])) + @constraint(model, x >= 3, Disjunct(W[2])) + @disjunction(model, W, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, integral(x, t)) + optimize!(model, gdp_method = MOIDisjunction()) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 10.0 atol = 1e-4 +end + +@testset "DisjunctiveProgramming integration" begin + test_lowering_solve_linear() + test_lowering_solve_nonlinear() + test_lowering_solve_nested() + test_lowering_infinite() + test_lowering_infinite_nested() +end diff --git a/test/loa.jl b/test/loa.jl new file mode 100644 index 0000000..5102e92 --- /dev/null +++ b/test/loa.jl @@ -0,0 +1,791 @@ +using JuMP +import HiGHS, Ipopt + +function _loa_optimizer(; kwargs...) + return () -> DA.Optimizer(; nlp_solver = Ipopt.Optimizer, + mip_solver = HiGHS.Optimizer, kwargs...) +end + +# min x with x >= 2 (disjunct 1) or x >= 5 (disjunct 2). The loop +# enumerates both combinations and keeps the incumbent at 2. +function test_linear_disjunction() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.GreaterThan(2.0)], [MOI.GreaterThan(5.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test primal_status(model) == MOI.FEASIBLE_POINT + @test objective_value(model) ≈ 2.0 atol = 1e-5 + @test value(x) ≈ 2.0 atol = 1e-5 + @test value(z[1]) ≈ 1.0 atol = 1e-5 + @test value(z[2]) ≈ 0.0 atol = 1e-5 + @test occursin("LOA finished", raw_status(model)) + @test solve_time(model) > 0.0 + @test dual_status(model) == MOI.NO_SOLUTION +end + +# Convex quadratic objective over linear disjuncts: y >= x or +# y >= 2 - x. The optimum sits at x = 3, y = 0 in disjunct 2. +function test_quadratic_objective() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 4) + @variable(model, 0 <= y <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], y - x, y - (2 - x)] in + DA.DisjunctionSet([ + [MOI.GreaterThan(0.0)], [MOI.GreaterThan(0.0)]])) + @objective(model, Min, (x - 3)^2 + y) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 0.0 atol = 1e-5 + @test value(z[2]) ≈ 1.0 atol = 1e-5 + @test value(x) ≈ 3.0 atol = 1e-4 + @test value(y) ≈ 0.0 atol = 1e-5 + # exhaustion path: the bound is the last master's, still valid + @test objective_bound(model) <= objective_value(model) + 1e-6 + @test !isnan(relative_gap(model)) +end + +# max x with (x <= 3) or (x <= 7): port of DP.jl test_loa_solve_simple. +function test_max_sense_linear() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(7.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 7.0 atol = 1e-4 + @test value(z[2]) ≈ 1.0 atol = 1e-5 +end + +# Two disjunctions: port of DP.jl test_loa_solve_two_disjunctions. +function test_two_disjunctions() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, 0 <= w <= 10) + @variable(model, zx[1:2], Bin) + @variable(model, zw[1:2], Bin) + @constraint(model, [1, zx[1], zx[2], x, x] in DA.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(7.0)]])) + @constraint(model, [1, zw[1], zw[2], w, w] in DA.DisjunctionSet([ + [MOI.LessThan(2.0)], [MOI.LessThan(5.0)]])) + @objective(model, Max, x + w) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 12.0 atol = 1e-4 +end + +# Nonlinear global x^2 <= 25 binds before the chosen disjunct: port of +# DP.jl test_loa_nonlinear_global (no Juniper needed - the layer's NLP +# has its binaries fixed). +function test_nonlinear_global() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, x^2 <= 25) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(8.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 5.0 atol = 1e-3 + @test value(z[2]) ≈ 1.0 atol = 1e-5 + @test isfinite(relative_gap(model)) +end + +# Nonlinear global equality: one seed combination is NLP-infeasible, so +# the NLPF path supplies the linearization site. Port of DP.jl +# test_loa_nonlinear_equality_global. +function test_nonlinear_equality_global() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, x^2 == 25) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(8.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 5.0 atol = 1e-3 + @test value(x) ≈ 5.0 atol = 1e-3 + @test value(z[2]) ≈ 1.0 atol = 1e-5 +end + +# Nonlinear equality inside a disjunct: set covering must activate it +# and the cut emits both gated directions. Port of DP.jl +# test_loa_nonlinear_equality_disjunct. +function test_nonlinear_equality_disjunct() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, 1.0z[1], 1.0z[2], 1.0x, x^2] in + DA.DisjunctionSet([[MOI.LessThan(3.0)], [MOI.EqualTo(64.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 8.0 atol = 1e-3 + @test value(z[2]) ≈ 1.0 atol = 1e-5 + # deterministic gap convergence: set covering visits the nonlinear + # disjunct first, then the master proves the other is worse + @test occursin("converged", raw_status(model)) + @test relative_gap(model) <= 1e-5 +end + +# Nonlinear Interval row inside a disjunct: port of DP.jl +# test_loa_nonlinear_interval_disjunct. +function test_nonlinear_interval_disjunct() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, 1.0z[1], 1.0z[2], 1.0x, x^2] in + DA.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.Interval(36.0, 64.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 8.0 atol = 1e-3 + @test value(z[2]) ≈ 1.0 atol = 1e-5 +end + +# A single binary drives both disjuncts through its complement: port of +# DP.jl test_loa_complement_indicator_nonlinear_disjunct. +function test_complement_indicator() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z, Bin) + @constraint(model, [1, 1.0z, 1 - z, 1.0x, x^2] in + DA.DisjunctionSet([[MOI.LessThan(3.0)], [MOI.EqualTo(64.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 8.0 atol = 1e-3 + @test value(z) ≈ 0.0 atol = 1e-5 +end + +# `use_nlpf = false` still solves by enumeration when the seeds are +# feasible. +function test_nlpf_disabled() + model = Model(_loa_optimizer(use_nlpf = false)) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, 1.0z[1], 1.0z[2], 1.0x, x^2] in + DA.DisjunctionSet([[MOI.LessThan(3.0)], [MOI.EqualTo(64.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 8.0 atol = 1e-3 +end + +# Big-M master gating solves the same instances as indicator gating +# (interval split included). +function test_bigm_master_gating() + model = Model(_loa_optimizer(master_gating = "bigm", M_value = 100.0)) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, 1.0z[1], 1.0z[2], 1.0x, 1.0x] in + DA.DisjunctionSet([ + [MOI.LessThan(2.0)], [MOI.Interval(3.0, 7.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 7.0 atol = 1e-4 +end + +# Integer-typed options convert at their use sites, including the +# Bool read of use_nlpf on the infeasible-seed path. +function test_integer_options() + model = Model(_loa_optimizer(use_nlpf = 0, M_value = 10^9, + time_limit = 3600)) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, 1.0z[1], 1.0z[2], 1.0x, x^2] in + DA.DisjunctionSet([[MOI.LessThan(3.0)], [MOI.EqualTo(64.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 8.0 atol = 1e-3 +end + +# A binary that is not a disjunction indicator keeps its integrality +# in the subproblem, which the nlp_solver must then handle (HiGHS +# both roles here since the model is linear). +function test_non_indicator_binary() + factory = () -> DA.Optimizer(nlp_solver = HiGHS.Optimizer, + mip_solver = HiGHS.Optimizer) + model = Model(factory) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @variable(model, w, Bin) + @constraint(model, 1.0x + w <= 8) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(7.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 7.0 atol = 1e-5 +end + +# Globally infeasible model: the master is infeasible before any +# incumbent exists. +function test_infeasible_no_incumbent() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, 1.0x >= 20) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(7.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.INFEASIBLE + @test primal_status(model) == MOI.NO_SOLUTION + @test result_count(model) == 0 +end + +# A zero time limit exits before any solve. +function test_time_limit() + model = Model(_loa_optimizer(time_limit = 0.0)) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(7.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.TIME_LIMIT + @test result_count(model) == 0 +end + +# Nonconvex objective: the OA bound is no certificate, so the layer +# must never report OPTIMAL and the raw status carries the gap record. +function test_nonconvex_never_optimal() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 2) + @variable(model, z[1:2], Bin) + @constraint(model, [1, 1.0z[1], 1.0z[2], 1.0x, 1.0x] in + DA.DisjunctionSet([[MOI.LessThan(1.0)], [MOI.LessThan(2.0)]])) + @objective(model, Min, -x^2) + optimize!(model) + @test termination_status(model) != MOI.OPTIMAL + @test termination_status(model) in + (MOI.LOCALLY_SOLVED, MOI.ITERATION_LIMIT, MOI.TIME_LIMIT) + @test primal_status(model) == MOI.FEASIBLE_POINT + @test occursin("LOA finished", raw_status(model)) +end + +# A feasibility model (no objective) minimizes a constant zero. +function test_feasibility_sense() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + set_start_value(x, 6.0) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.GreaterThan(5.0)], [MOI.LessThan(1.0)]])) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test primal_status(model) == MOI.FEASIBLE_POINT + @test value(x) >= 5.0 - 1e-6 || value(x) <= 1.0 + 1e-6 +end + +# A linear Interval disjunct row reaches the master as two one-sided +# indicator constraints (Indicator{A}(Interval) has no MILP bridge). +function test_linear_interval_disjunct() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, 1.0z[1], 1.0z[2], 1.0x, 1.0x] in + DA.DisjunctionSet([ + [MOI.LessThan(2.0)], [MOI.Interval(3.0, 7.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 7.0 atol = 1e-4 + @test value(z[2]) ≈ 1.0 atol = 1e-5 +end + +# A second solve on the same optimizer must not leak the first solve's +# bound, gap, or incumbent. +function test_reoptimize_resets_results() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.GreaterThan(2.0)], [MOI.GreaterThan(5.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 2.0 atol = 1e-5 + @constraint(model, 1.0x >= 20) + optimize!(model) + @test termination_status(model) == MOI.INFEASIBLE + @test result_count(model) == 0 + @test objective_bound(model) == -Inf + @test isnan(relative_gap(model)) +end + +# An unsupported constraint type is rewritten by the JuMP bridge layer +# into supported rows instead of erroring inside the solve. +function test_bridged_vector_constraint() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1.0x - 5.0] in MOI.Nonnegatives(1)) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(7.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 5.0 atol = 1e-4 + @test value(z[2]) ≈ 1.0 atol = 1e-5 +end + +# Genuinely nonlinear (non-quadratic) global rows stay +# `ScalarNonlinearFunction`s through demotion and the linearizer walks +# their argument trees (variable and affine arguments). +function test_nonlinear_exp_global() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, exp(x) <= exp(5.0)) + @constraint(model, exp(x + 1.0) <= exp(6.5)) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(8.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 5.0 atol = 1e-3 + @test value(z[2]) ≈ 1.0 atol = 1e-5 +end + +# `use_nlpf = false` with an infeasible combination: the solve keeps +# only the no-good cut and still finishes from the other disjunct. +function test_nlpf_disabled_infeasible_combination() + model = Model(_loa_optimizer(use_nlpf = false)) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, 1.0x >= 5) + @constraint(model, [1, 1.0z[1], 1.0z[2], x^2, x^2] in + DA.DisjunctionSet([[MOI.LessThan(9.0)], [MOI.LessThan(64.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 5.0 atol = 1e-3 + @test value(z[2]) ≈ 1.0 atol = 1e-5 +end + +# The NLPF slacks GreaterThan rows and copies the global linear rows +# when restoring an infeasible combination. +function test_nlpf_slacked_rows() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, 1.0x >= 4) + @constraint(model, x^2 >= 25) + @constraint(model, [1, 1.0z[1], 1.0z[2], x^2, x^2] in + DA.DisjunctionSet([[MOI.LessThan(9.0)], [MOI.LessThan(64.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 5.0 atol = 1e-3 + @test value(z[2]) ≈ 1.0 atol = 1e-5 +end + +# A quadratic-typed global row with no quadratic terms demotes to +# affine and lands in the master directly. +function test_demoted_affine_global() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, zero(QuadExpr) + 1.0x <= 6) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(8.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 6.0 atol = 1e-4 +end + +# An all-variable disjunction arrives as `VectorOfVariables`; its raw +# variable rows are promoted to affine rows rather than bounds. +function test_vector_of_variables_disjunction() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, zout[1:2], Bin) + @variable(model, zin[1:2], Bin) + @constraint(model, [1, zout[1], zout[2], x] in DA.DisjunctionSet([ + [MOI.LessThan(2.0)], MOI.AbstractScalarSet[]])) + @constraint(model, [zout[2], zin[1], zin[2], x, x] in + DA.DisjunctionSet([[MOI.LessThan(5.0)], [MOI.LessThan(8.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 8.0 atol = 1e-4 + @test value(zout[2]) ≈ 1.0 atol = 1e-5 + @test value(zin[2]) ≈ 1.0 atol = 1e-5 +end + +# Zero iteration budgets exit before any solve, without an incumbent. +function test_iteration_limit_no_incumbent() + model = Model(_loa_optimizer(set_cover_max_iter = 0, max_iter = 0)) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(7.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.ITERATION_LIMIT + @test result_count(model) == 0 +end + +# `max_iter = 0` keeps the set-covering incumbent but produces no +# master bound. +function test_iteration_limit_with_incumbent() + model = Model(_loa_optimizer(max_iter = 0)) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.GreaterThan(2.0)], [MOI.GreaterThan(5.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.ITERATION_LIMIT + @test primal_status(model) == MOI.FEASIBLE_POINT + @test occursin("limit hit", raw_status(model)) + @test occursin("no bound", raw_status(model)) + @test objective_bound(model) == -Inf +end + +# An unbounded master (no OA cuts yet bound alpha_oa) surfaces its +# status instead of looping. +function test_master_abnormal_status() + model = Model(_loa_optimizer(set_cover_max_iter = 0)) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.LessThan(3.0)], [MOI.LessThan(7.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.OTHER_LIMIT + @test result_count(model) == 0 + @test occursin("master solve finished", raw_status(model)) +end + +################################################################################ +# MOCK SOLVER +################################################################################ +# Delegates every MOI call to a wrapped optimizer, but sleeps +# `sleep_time` seconds in each solve and, from solve `fail_from` on, +# skips the inner solve and reports `fail_status` with no solution. +# Deterministic triggers for the deadline and abnormal-master paths. +mutable struct MockSolver <: MOI.AbstractOptimizer + inner::MOI.AbstractOptimizer + sleep_time::Float64 + fail_from::Int + fail_status::MOI.TerminationStatusCode + solves::Int + failing::Bool +end + +function MockSolver( + factory; + sleep_time::Float64 = 0.0, + fail_from::Int = typemax(Int), + fail_status::MOI.TerminationStatusCode = MOI.NODE_LIMIT + ) + return MockSolver(MOI.instantiate(factory), sleep_time, fail_from, + fail_status, 0, false) +end + +function MOI.optimize!(model::MockSolver) + model.solves += 1 + model.sleep_time > 0 && sleep(model.sleep_time) + model.failing = model.solves >= model.fail_from + model.failing || MOI.optimize!(model.inner) + return +end + +const _WrappedAttr = Union{MOI.AbstractModelAttribute, + MOI.AbstractOptimizerAttribute} +const _WrappedIndexAttr = Union{MOI.AbstractVariableAttribute, + MOI.AbstractConstraintAttribute} +const _WrappedIndex = Union{MOI.VariableIndex, MOI.ConstraintIndex} + +function MOI.get(model::MockSolver, attr::_WrappedAttr) + if model.failing + attr isa MOI.TerminationStatus && return model.fail_status + attr isa MOI.PrimalStatus && return MOI.NO_SOLUTION + end + return MOI.get(model.inner, attr) +end + +MOI.is_empty(model::MockSolver) = MOI.is_empty(model.inner) +MOI.empty!(model::MockSolver) = MOI.empty!(model.inner) +MOI.supports_incremental_interface(::MockSolver) = true +MOI.copy_to(model::MockSolver, src::MOI.ModelLike) = + MOI.copy_to(model.inner, src) +MOI.add_variable(model::MockSolver) = MOI.add_variable(model.inner) +MOI.delete(model::MockSolver, index) = MOI.delete(model.inner, index) +MOI.is_valid(model::MockSolver, index) = MOI.is_valid(model.inner, index) + +function MOI.add_constraint( + model::MockSolver, + func::MOI.AbstractFunction, + set::MOI.AbstractSet + ) + return MOI.add_constraint(model.inner, func, set) +end + +function MOI.supports_constraint( + model::MockSolver, + F::Type{<:MOI.AbstractFunction}, + S::Type{<:MOI.AbstractSet} + ) + return MOI.supports_constraint(model.inner, F, S) +end + +MOI.supports(model::MockSolver, attr::_WrappedAttr) = + MOI.supports(model.inner, attr) + +MOI.set(model::MockSolver, attr::_WrappedAttr, value) = + MOI.set(model.inner, attr, value) + +function MOI.supports( + model::MockSolver, + attr::_WrappedIndexAttr, + I::Type{<:_WrappedIndex} + ) + return MOI.supports(model.inner, attr, I) +end + +function MOI.get( + model::MockSolver, + attr::_WrappedIndexAttr, + index::_WrappedIndex + ) + return MOI.get(model.inner, attr, index) +end + +function MOI.set( + model::MockSolver, + attr::_WrappedIndexAttr, + index::_WrappedIndex, + value + ) + return MOI.set(model.inner, attr, index, value) +end + +function _mock_time_limit_model(limit::Float64; kwargs...) + factory = () -> DA.Optimizer( + nlp_solver = () -> MockSolver(Ipopt.Optimizer; kwargs...), + mip_solver = HiGHS.Optimizer, + iteration_time_limit = limit) + model = Model(factory) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.GreaterThan(2.0)], [MOI.GreaterThan(5.0)]])) + @objective(model, Min, x) + return model +end + +# The loop deadline passes right after the covering incumbent: the +# mock NLP sleeps past `iteration_time_limit`, so the main loop never +# starts and the incumbent is reported against the time limit. The +# first solve warms up the mock stack so compilation latency cannot +# eat the deadline before the covering pass runs. +function test_time_limit_with_incumbent() + warmup = _mock_time_limit_model(Inf) + optimize!(warmup) + @test objective_value(warmup) ≈ 2.0 atol = 1e-5 + model = _mock_time_limit_model(0.5, sleep_time = 1.5) + optimize!(model) + @test termination_status(model) == MOI.TIME_LIMIT + @test primal_status(model) == MOI.FEASIBLE_POINT + @test occursin("limit hit", raw_status(model)) +end + +# The master finishes abnormally after an incumbent exists: the mock +# master reports a node limit on its second solve. +function test_master_abnormal_status_with_incumbent() + factory = () -> DA.Optimizer( + nlp_solver = Ipopt.Optimizer, + mip_solver = () -> MockSolver(HiGHS.Optimizer, fail_from = 2)) + model = Model(factory) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, z[1:2], Bin) + @constraint(model, [1, z[1], z[2], x, x] in DA.DisjunctionSet([ + [MOI.GreaterThan(2.0)], [MOI.GreaterThan(5.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.OTHER_LIMIT + @test primal_status(model) == MOI.FEASIBLE_POINT + @test occursin("limit hit", raw_status(model)) +end + +################################################################################ +# UNIT TESTS +################################################################################ +function test_cut_term_directions() + x = MOI.VariableIndex(1) + lin = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(2.0, x)], 1.0) + @test length(DA._oa_cut_terms(MOI.LessThan(5.0), lin)) == 1 + @test length(DA._oa_cut_terms(MOI.GreaterThan(5.0), lin)) == 1 + @test length(DA._oa_cut_terms(MOI.EqualTo(5.0), lin)) == 2 + @test length(DA._oa_cut_terms(MOI.Interval(0.0, 5.0), lin)) == 2 + less = only(DA._oa_cut_terms(MOI.LessThan(5.0), lin)) + @test less.constant == -4.0 + greater = only(DA._oa_cut_terms(MOI.GreaterThan(5.0), lin)) + @test greater.constant == 4.0 + @test only(greater.terms).coefficient == -2.0 +end + +function test_activation_binary() + z = MOI.VariableIndex(1) + saf(a, c) = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(a, z)], c) + @test DA._activation_binary(saf(1.0, 0.0)) == (z, true) + @test DA._activation_binary(saf(-1.0, 1.0)) == (z, false) + @test_throws ErrorException DA._activation_binary(saf(2.0, 0.0)) + @test_throws ErrorException DA._activation_binary(saf(1.0, 3.0)) +end + +function test_sense_primitives() + @test DA._penalty_sign(MOI.MIN_SENSE) == 1 + @test DA._penalty_sign(MOI.MAX_SENSE) == -1 + @test DA._worst_objective(MOI.MIN_SENSE) == Inf + @test DA._worst_objective(MOI.MAX_SENSE) == -Inf + @test DA._is_better(MOI.MIN_SENSE, 1.0, 2.0) + @test DA._is_better(MOI.MAX_SENSE, 2.0, 1.0) + @test DA._gap(MOI.MIN_SENSE, 5.0, 3.0) == 2.0 + @test DA._gap(MOI.MAX_SENSE, 3.0, 5.0) == 2.0 +end + +function test_linearize_quadratic() + x = MOI.VariableIndex(1) + func = MOI.ScalarQuadraticFunction( + [MOI.ScalarQuadraticTerm(2.0, x, x)], + MOI.ScalarAffineTerm{Float64}[], 0.0) # x^2 + linearizer = DA._Linearizer() + lin = DA._linearize(linearizer, func, Dict(x => 3.0)) + # x^2 at x = 3: 9 + 6 (x - 3) = 6 x - 9 + @test lin.constant ≈ -9.0 + @test only(lin.terms).coefficient ≈ 6.0 + @test length(linearizer.evaluators) == 1 + DA._linearize(linearizer, func, Dict(x => 4.0)) + @test length(linearizer.evaluators) == 1 +end + +# Nested disjunction: the inner disjunction's activation component is +# the outer indicator, so it selects a mode only while the outer +# disjunct is active and is vacuous otherwise. +function test_nested_disjunction() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, zout[1:2], Bin) + @variable(model, zin[1:2], Bin) + @constraint(model, [1, zout[1], zout[2], x] in DA.DisjunctionSet([ + [MOI.LessThan(2.0)], MOI.AbstractScalarSet[]])) + @constraint(model, [1.0zout[2], zin[1], zin[2], x^2, x] in + DA.DisjunctionSet([[MOI.LessThan(25.0)], [MOI.LessThan(8.0)]])) + @objective(model, Max, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 8.0 atol = 1e-3 + @test value(zout[2]) ≈ 1.0 atol = 1e-5 + @test value(zin[2]) ≈ 1.0 atol = 1e-5 +end + +# When the parent disjunct is off, the inner indicators sum to zero +# and the inner rows impose nothing. +function test_nested_disjunction_vacuous() + model = Model(_loa_optimizer()) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, zout[1:2], Bin) + @variable(model, zin[1:2], Bin) + @constraint(model, [1, zout[1], zout[2], x] in DA.DisjunctionSet([ + [MOI.LessThan(2.0)], MOI.AbstractScalarSet[]])) + @constraint(model, [1.0zout[2], zin[1], zin[2], x, x^2] in + DA.DisjunctionSet([ + [MOI.GreaterThan(5.0)], [MOI.GreaterThan(36.0)]])) + @objective(model, Min, x) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 0.0 atol = 1e-4 + @test value(zout[1]) ≈ 1.0 atol = 1e-5 + @test value(zin[1]) + value(zin[2]) ≈ 0.0 atol = 1e-5 +end + +@testset "LOA loop" begin + test_linear_disjunction() + test_nested_disjunction() + test_nested_disjunction_vacuous() + test_quadratic_objective() + test_max_sense_linear() + test_two_disjunctions() + test_nonlinear_global() + test_nonlinear_exp_global() + test_nonlinear_equality_global() + test_nonlinear_equality_disjunct() + test_nonlinear_interval_disjunct() + test_complement_indicator() + test_nlpf_disabled() + test_nlpf_disabled_infeasible_combination() + test_nlpf_slacked_rows() + test_demoted_affine_global() + test_vector_of_variables_disjunction() + test_bigm_master_gating() + test_integer_options() + test_non_indicator_binary() + test_infeasible_no_incumbent() + test_time_limit() + test_iteration_limit_no_incumbent() + test_iteration_limit_with_incumbent() + test_master_abnormal_status() + test_time_limit_with_incumbent() + test_master_abnormal_status_with_incumbent() + test_nonconvex_never_optimal() + test_feasibility_sense() + test_linear_interval_disjunct() + test_reoptimize_resets_results() + test_bridged_vector_constraint() +end + +@testset "LOA units" begin + test_cut_term_directions() + test_activation_binary() + test_sense_primitives() + test_linearize_quadratic() +end diff --git a/test/moi.jl b/test/moi.jl new file mode 100644 index 0000000..25fd8b1 --- /dev/null +++ b/test/moi.jl @@ -0,0 +1,35 @@ +# MOI contract conformance. MOI.Test cannot build a `DisjunctionSet`, so it +# only exercises the API-contract surface (model API + optimizer attributes), +# not the LOA solving logic -- that is covered by the other test files. The +# solving families and the dual/basis attributes are out of scope by design. +import HiGHS +import Ipopt + +@testset "MOI contract (model API + attributes)" begin + optimizer = MOI.instantiate( + () -> DA.Optimizer( + nlp_solver = Ipopt.Optimizer, + mip_solver = HiGHS.Optimizer, + time_limit = 20.0, + ); + with_cache_type = Float64, + with_bridge_type = Float64, + ) + config = MOI.Test.Config( + atol = 1e-3, + rtol = 1e-3, + optimal_status = MOI.LOCALLY_SOLVED, # never reports OPTIMAL + exclude = Any[ + MOI.ConstraintDual, + MOI.DualObjectiveValue, + MOI.ConstraintBasisStatus, + MOI.VariableBasisStatus, + ], + ) + MOI.Test.runtests( + optimizer, + config; + include = Regex[r"test_model_", r"test_attribute_"], + warn_unsupported = false, + ) +end diff --git a/test/optimizer.jl b/test/optimizer.jl new file mode 100644 index 0000000..e5f004c --- /dev/null +++ b/test/optimizer.jl @@ -0,0 +1,127 @@ +# A two-disjunct toy stored straight into a model: binaries z1, z2 and +# continuous x with x <= 0 (disjunct 1) or x >= 1 (disjunct 2). +function _build_toy_disjunction(model) + x = MOI.add_variable(model) + z = MOI.add_variables(model, 2) + for zi in z + MOI.add_constraint(model, zi, MOI.ZeroOne()) + end + set = DA.DisjunctionSet([[MOI.LessThan(0.0)], [MOI.GreaterThan(1.0)]]) + func = MOI.VectorAffineFunction( + [MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, z[1])), + MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(1.0, z[2])), + MOI.VectorAffineTerm(4, MOI.ScalarAffineTerm(1.0, x)), + MOI.VectorAffineTerm(5, MOI.ScalarAffineTerm(1.0, x))], + [1.0, 0.0, 0.0, 0.0, 0.0]) # constant activation in component 1 + ci = MOI.add_constraint(model, func, set) + return x, z, set, ci +end + +function test_optimizer_scaffold() + optimizer = DA.Optimizer() + @test MOI.is_empty(optimizer) + @test MOI.get(optimizer, MOI.SolverName()) == "DisjunctiveAlgorithms" + @test MOI.get(optimizer, MOI.TerminationStatus()) == + MOI.OPTIMIZE_NOT_CALLED + @test MOI.get(optimizer, MOI.ResultCount()) == 0 + @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.NO_SOLUTION + @test MOI.get(optimizer, MOI.DualStatus()) == MOI.NO_SOLUTION +end + +function test_optimizer_options() + optimizer = DA.Optimizer(max_iter = 25) + @test MOI.get(optimizer, MOI.RawOptimizerAttribute("max_iter")) == 25 + MOI.set(optimizer, MOI.RawOptimizerAttribute("max_iter"), 5) + @test MOI.get(optimizer, MOI.RawOptimizerAttribute("max_iter")) == 5 + @test MOI.get(optimizer, MOI.RawOptimizerAttribute("M_value")) == 1e9 + bogus = MOI.RawOptimizerAttribute("bogus") + @test !MOI.supports(optimizer, bogus) + @test_throws MOI.UnsupportedAttribute MOI.get(optimizer, bogus) + @test_throws MOI.UnsupportedAttribute MOI.set(optimizer, bogus, 1) + @test_throws ArgumentError DA.Optimizer(bogus = 1) + MOI.set(optimizer, MOI.Silent(), true) + @test MOI.get(optimizer, MOI.Silent()) + MOI.set(optimizer, MOI.TimeLimitSec(), 100) + @test MOI.get(optimizer, MOI.TimeLimitSec()) == 100.0 + @test MOI.get(optimizer, MOI.RawOptimizerAttribute("time_limit")) == 100.0 + MOI.set(optimizer, MOI.TimeLimitSec(), nothing) + @test MOI.get(optimizer, MOI.TimeLimitSec()) === nothing +end + +function test_model_building() + optimizer = DA.Optimizer() + x, z, set, ci = _build_toy_disjunction(optimizer) + @test MOI.is_valid(optimizer, ci) + @test MOI.get(optimizer, MOI.NumberOfVariables()) == 3 + @test (MOI.VectorAffineFunction{Float64}, DA.DisjunctionSet) in + MOI.get(optimizer, MOI.ListOfConstraintTypesPresent()) + @test MOI.get(optimizer, MOI.ConstraintSet(), ci) == set + objective = MOI.ScalarAffineFunction( + [MOI.ScalarAffineTerm(1.0, x)], 0.0) + MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.set(optimizer, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + objective) + @test MOI.get(optimizer, MOI.ObjectiveSense()) == MOI.MIN_SENSE + @test !MOI.is_empty(optimizer) + MOI.empty!(optimizer) + @test MOI.is_empty(optimizer) +end + +function test_copy_to() + src = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) + x, z, set, ci = _build_toy_disjunction(src) + MOI.set(src, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.set(src, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0)) + dest = DA.Optimizer() + index_map = MOI.copy_to(dest, src) + @test MOI.get(dest, MOI.NumberOfVariables()) == 3 + @test (MOI.VectorAffineFunction{Float64}, DA.DisjunctionSet) in + MOI.get(dest, MOI.ListOfConstraintTypesPresent()) + @test MOI.get(dest, MOI.ConstraintSet(), index_map[ci]) == set + dest_func = MOI.get(dest, MOI.ConstraintFunction(), index_map[ci]) + @test dest_func.terms[1].scalar_term.variable == index_map[z[1]] +end + +# The incremental MOI surface forwards to the cache: names, starts, +# constrained variables, and deletion. +function test_moi_forwarding() + optimizer = DA.Optimizer() + @test MOI.supports_incremental_interface(optimizer) + x = MOI.add_variable(optimizer) + MOI.set(optimizer, MOI.VariableName(), x, "x") + @test MOI.get(optimizer, MOI.VariableIndex, "x") == x + MOI.set(optimizer, MOI.VariablePrimalStart(), x, 2.5) + @test MOI.get(optimizer, MOI.VariablePrimalStart(), x) == 2.5 + y, ci_y = MOI.add_constrained_variables(optimizer, MOI.Nonnegatives(2)) + @test length(y) == 2 + @test MOI.is_valid(optimizer, ci_y) + func = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0) + ci = MOI.add_constraint(optimizer, func, MOI.LessThan(1.0)) + MOI.set(optimizer, MOI.ConstraintName(), ci, "c") + @test MOI.get(optimizer, typeof(ci), "c") == ci + MOI.delete(optimizer, ci) + @test !MOI.is_valid(optimizer, ci) + @test !MOI.supports_add_constrained_variables(optimizer, + DA.DisjunctionSet) +end + +# A constraint type the partition cannot place errors at solve time. +function test_unsupported_constraint_type() + optimizer = DA.Optimizer() + x = MOI.add_variable(optimizer) + func = MOI.VectorAffineFunction( + [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x))], [0.0]) + MOI.add_constraint(optimizer, func, MOI.Nonnegatives(1)) + @test_throws ErrorException MOI.optimize!(optimizer) +end + +@testset "Optimizer scaffold" begin + test_optimizer_scaffold() + test_optimizer_options() + test_model_building() + test_copy_to() + test_moi_forwarding() + test_unsupported_constraint_type() +end diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..f659abc --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,9 @@ +using Test +import MathOptInterface as MOI +using DisjunctiveAlgorithms +const DA = DisjunctiveAlgorithms + +include("optimizer.jl") +include("loa.jl") +include("moi.jl") +include("integration.jl") From 07be35e031f01202c8c7eeb1ef60c58d7259622a Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Sat, 15 Aug 2026 06:17:56 -0400 Subject: [PATCH 02/18] Dependency update --- .github/workflows/CI.yml | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 2e521a4..993afe0 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -37,6 +37,11 @@ jobs: 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 @@ -58,6 +63,13 @@ jobs: 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: From 517985fa7c8237e1ef27034956e513012b99bb2b Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Sat, 15 Aug 2026 07:41:22 -0400 Subject: [PATCH 03/18] Docstring fix --- src/optimizer.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/optimizer.jl b/src/optimizer.jl index b55dd60..29e44af 100644 --- a/src/optimizer.jl +++ b/src/optimizer.jl @@ -22,10 +22,10 @@ const _DEFAULT_OPTIONS = Dict{String, Any}( Optimizer(; nlp_solver, mip_solver = nlp_solver, kwargs...) Logic-based outer approximation solver for models containing -[`DisjunctionSet`](@ref) constraints. `nlp_solver` and `mip_solver` -are optimizer factories as accepted by `MOI.instantiate`. The -remaining keyword arguments set raw options (also reachable through -`MOI.RawOptimizerAttribute`): +`DisjunctiveProgramming.DisjunctionSet` constraints. `nlp_solver` +and `mip_solver` are optimizer factories as accepted by +`MOI.instantiate`. The remaining keyword arguments set raw options +(also reachable through `MOI.RawOptimizerAttribute`): - `max_iter = 10`: master/NLP iterations after the set-covering seed. - `set_cover_max_iter = 8`: set-covering initialization iterations. From ebdf742368f8ea9146495dcd987740e209844813 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Sat, 15 Aug 2026 11:58:47 -0400 Subject: [PATCH 04/18] Shift row function constants --- src/problem.jl | 10 +++++++--- test/loa.jl | 20 ++++++++++++++++++++ 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/src/problem.jl b/src/problem.jl index ebce607..3f7b959 100644 --- a/src/problem.jl +++ b/src/problem.jl @@ -96,10 +96,14 @@ function _parse_disjunction( binary.value) MOI.is_valid(cache, zero_one) || error("The indicator variable " * "of a `DisjunctionSet` disjunct must be `MOI.ZeroOne`.") - functions = MOI.AbstractScalarFunction[ - _as_row(_demote(rows[k])) for k in row_indices(set, i)] + # 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, - functions, set.inner_sets[i])) + MOI.AbstractScalarFunction[f for (f, _) in normalized], + MOI.AbstractScalarSet[s for (_, s) in normalized])) end return _Disjunction(disjunction_activation, disjuncts) end diff --git a/test/loa.jl b/test/loa.jl index 5102e92..df9a662 100644 --- a/test/loa.jl +++ b/test/loa.jl @@ -28,6 +28,25 @@ function test_linear_disjunction() @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. function test_quadratic_objective() @@ -750,6 +769,7 @@ end @testset "LOA loop" begin test_linear_disjunction() + test_row_function_constants() test_nested_disjunction() test_nested_disjunction_vacuous() test_quadratic_objective() From 7f04fc4766504730d71aef2c7274e85a06f19389 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Sun, 16 Aug 2026 15:53:42 -0400 Subject: [PATCH 05/18] More style changes --- src/DisjunctiveAlgorithms.jl | 2 +- src/{loa.jl => algorithms/LOA.jl} | 85 ++++++++++-- src/cuts.jl | 8 +- src/master.jl | 19 +-- src/nlp.jl | 2 +- src/optimizer.jl | 221 +++++++++++++++++++++--------- test/integration.jl | 3 +- test/loa.jl | 193 ++++++++------------------ test/mock_optimizer.jl | 113 +++++++++++++++ test/moi.jl | 7 +- test/optimizer.jl | 37 ++--- test/runtests.jl | 1 + 12 files changed, 446 insertions(+), 245 deletions(-) rename src/{loa.jl => algorithms/LOA.jl} (77%) create mode 100644 test/mock_optimizer.jl diff --git a/src/DisjunctiveAlgorithms.jl b/src/DisjunctiveAlgorithms.jl index 6831d27..0985c40 100644 --- a/src/DisjunctiveAlgorithms.jl +++ b/src/DisjunctiveAlgorithms.jl @@ -10,6 +10,6 @@ include("problem.jl") include("master.jl") include("nlp.jl") include("cuts.jl") -include("loa.jl") +include("algorithms/LOA.jl") end diff --git a/src/loa.jl b/src/algorithms/LOA.jl similarity index 77% rename from src/loa.jl rename to src/algorithms/LOA.jl index d68fdd6..05089e9 100644 --- a/src/loa.jl +++ b/src/algorithms/LOA.jl @@ -1,6 +1,73 @@ ################################################################################ # 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, seeded by a +set-covering pass over the nonlinear disjuncts. The OA bound is only +valid on convex problems, so convergence reports `LOCALLY_SOLVED`, +never `OPTIMAL`. + +## Supported optimizer attributes + +- [`NumIterationLimit`](@ref) +- [`SetCoverIterationLimit`](@ref) +- [`M_Value`](@ref) +- [`MasterReformulation`](@ref) +- [`MaxSlack`](@ref) +- [`OASlack`](@ref) +- [`UseNLPF`](@ref) +- [`ConvergenceTolerance`](@ref) +- [`SlackTolerance`](@ref) +- [`IterationTimeLimit`](@ref) +""" +mutable struct LOA <: AbstractAlgorithm + num_iteration_limit::Union{Nothing, Int} + set_cover_iteration_limit::Union{Nothing, Int} + 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} + + LOA() = new(nothing, nothing, nothing, nothing, nothing, nothing, + nothing, nothing, nothing, nothing) +end + +_default(::Algorithm) = LOA() + +# An unset field (`nothing`) reads back as the attribute default. +for (attr, field) in ( + (NumIterationLimit, :num_iteration_limit), + (SetCoverIterationLimit, :set_cover_iteration_limit), + (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), + ) + @eval begin + MOI.supports(::LOA, ::$attr) = true + function MOI.set(algorithm::LOA, ::$attr, value) + algorithm.$field = value + return + end + function MOI.get(algorithm::LOA, attr::$attr) + return something(algorithm.$field, _default(algorithm, attr)) + end + end +end + _worst_objective(sense::MOI.OptimizationSense) = sense == MOI.MAX_SENSE ? -Inf : Inf _is_better(sense::MOI.OptimizationSense, new, best) = @@ -77,7 +144,7 @@ end ################################################################################ # MAIN LOOP ################################################################################ -function MOI.optimize!(model::Optimizer) +function _optimize!(algorithm::LOA, model::Optimizer) t_start = time() _reset_results(model) problem = _build_problem(model) @@ -85,9 +152,9 @@ function MOI.optimize!(model::Optimizer) subproblem = _build_subproblem(model, problem) linearizer = _Linearizer() sense = problem.sense - overall_deadline = t_start + Float64(_option(model, "time_limit")) + overall_deadline = t_start + something(model.time_limit_sec, Inf) loop_deadline = min(overall_deadline, - t_start + Float64(_option(model, "iteration_time_limit"))) + t_start + Float64(MOI.get(algorithm, IterationTimeLimit()))) best_objective = _worst_objective(sense) best_result = nothing @@ -117,7 +184,7 @@ function MOI.optimize!(model::Optimizer) cover = _cover_disjuncts(problem) needs_cover = trues(length(cover)) num_covered = 0 - for iteration in 1:_option(model, "set_cover_max_iter") + for iteration in 1:MOI.get(algorithm, SetCoverIterationLimit()) (iteration == 1 || any(needs_cover)) || break time() < loop_deadline || break _set_master_objective(master, MOI.MAX_SENSE, @@ -150,7 +217,7 @@ function MOI.optimize!(model::Optimizer) # main loop: alpha_oa gives the bound, the NLP the incumbent if master_status === nothing - for _ in 1:_option(model, "max_iter") + for _ in 1:MOI.get(algorithm, NumIterationLimit()) time() < loop_deadline || break _cap_remaining_time(master.model, loop_deadline) MOI.optimize!(master.model) @@ -164,11 +231,11 @@ function MOI.optimize!(model::Optimizer) gap = _gap(sense, best_objective, master_bound) total_slack = abs(MOI.get(master.model, MOI.ObjectiveValue()) - master_bound) / - Float64(_option(model, "oa_penalty")) - tol = Float64(_option(model, "convergence_tol")) * + Float64(MOI.get(algorithm, OASlack())) + tol = Float64(MOI.get(algorithm, ConvergenceTolerance())) * max(abs(best_objective), 1.0) - if gap <= tol && - total_slack <= Float64(_option(model, "slack_tol")) + if gap <= tol && total_slack <= + Float64(MOI.get(algorithm, SlackTolerance())) converged = true break end diff --git a/src/cuts.jl b/src/cuts.jl index 079cc2e..0c79bf8 100644 --- a/src/cuts.jl +++ b/src/cuts.jl @@ -145,7 +145,7 @@ function _add_objective_cut( sign::Int ) lin = _master_linearization(master, linearizer, master.objective, point) - slack = _add_penalized_slack(master, model.options, sign) + slack = _add_penalized_slack(master, model, sign) alpha = _to_affine(master.alpha_oa) if master.sense == MOI.MAX_SENSE body = MOI.Utilities.operate(-, Float64, @@ -167,7 +167,7 @@ function _add_global_oa_row( set::MOI.AbstractScalarSet, sign::Int ) - slack = _add_penalized_slack(master, model.options, sign) + 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, @@ -186,12 +186,12 @@ function _add_disjunct_oa_cut( set::MOI.AbstractScalarSet, sign::Int ) - M = Float64(_option(model, "M_value")) + 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.options, sign) + 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) diff --git a/src/master.jl b/src/master.jl index 7e35d78..a909839 100644 --- a/src/master.jl +++ b/src/master.jl @@ -58,7 +58,7 @@ function _build_master(model::Optimizer, problem::_Problem) for disjunction in problem.disjunctions _add_exactly_one(mip, variable_map, disjunction) for disjunct in disjunction.disjuncts - _add_gated_rows(mip, variable_map, disjunct, model.options) + _add_gated_rows(mip, variable_map, disjunct, model) end end alpha_oa = MOI.add_variable(mip) @@ -95,21 +95,22 @@ _indicator_sets(set::MOI.Interval{Float64}) = _indicator_sets(set::MOI.AbstractScalarSet) = (set,) # Gate each linear row with an indicator constraint or big-M per -# `master_gating`; nonlinear rows enter the master only as OA cuts +# `MasterReformulation`; nonlinear rows enter the master only as OA +# cuts function _add_gated_rows( mip::MOI.ModelLike, variable_map::AbstractDict, disjunct::_Disjunct, - options::Dict{String, Any} + model::Optimizer ) - gating = options["master_gating"] + gating = MOI.get(model, MasterReformulation()) gating in ("indicator", "bigm") || - error("Unknown `master_gating` value `$gating`.") + error("Unknown `MasterReformulation` value `$gating`.") activate = disjunct.active_value ? MOI.ACTIVATE_ON_ONE : MOI.ACTIVATE_ON_ZERO binary = variable_map[disjunct.binary] activation = _map_to(variable_map, disjunct.activation) - M = Float64(options["M_value"]) + M = 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) @@ -136,14 +137,14 @@ end # bounded penalized slack so an invalid cut cannot blow up the master function _add_penalized_slack( master::_Master, - options::Dict{String, Any}, + 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(options["max_slack"]))) - penalty = penalty_sign * Float64(options["oa_penalty"]) + 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)) diff --git a/src/nlp.jl b/src/nlp.jl index eaff857..0aded0b 100644 --- a/src/nlp.jl +++ b/src/nlp.jl @@ -92,7 +92,7 @@ function _solve_nlp( objective = MOI.get(sub.model, MOI.ObjectiveValue()), feasible = true) end - if Bool(_option(model, "use_nlpf")) + if Bool(MOI.get(model, UseNLPF())) result = _solve_nlpf(model, problem, combination, warm_start; deadline = deadline) result === nothing || return result diff --git a/src/optimizer.jl b/src/optimizer.jl index 29e44af..9d321b7 100644 --- a/src/optimizer.jl +++ b/src/optimizer.jl @@ -1,53 +1,154 @@ +################################################################################ +# 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. Defaults to `8`. +""" +struct SetCoverIterationLimit <: AbstractAlgorithmAttribute end +_default(::SetCoverIterationLimit) = 8 + +""" + 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 + ################################################################################ # OPTIMIZER ################################################################################ const _Cache = MOI.Utilities.UniversalFallback{MOI.Utilities.Model{Float64}} -# Raw options, mirroring DisjunctiveProgramming.jl's LOA defaults. -const _DEFAULT_OPTIONS = Dict{String, Any}( - "max_iter" => 10, - "set_cover_max_iter" => 8, - "M_value" => 1e9, - "master_gating" => "indicator", - "max_slack" => 1e3, - "oa_penalty" => 1e3, - "use_nlpf" => true, - "convergence_tol" => 1e-6, - "slack_tol" => 1e-4, - "iteration_time_limit" => Inf, - "time_limit" => 3600.0, -) - """ - Optimizer(; nlp_solver, mip_solver = nlp_solver, kwargs...) + Optimizer(nlp_solver, mip_solver = nlp_solver) -Logic-based outer approximation solver for models containing +Solver for models containing `DisjunctiveProgramming.DisjunctionSet` constraints. `nlp_solver` and `mip_solver` are optimizer factories as accepted by -`MOI.instantiate`. The remaining keyword arguments set raw options -(also reachable through `MOI.RawOptimizerAttribute`): - -- `max_iter = 10`: master/NLP iterations after the set-covering seed. -- `set_cover_max_iter = 8`: set-covering initialization iterations. -- `M_value = 1e9`: big-M gating the disjunct OA cuts in the master. -- `master_gating = "indicator"`: how linear disjunct rows enter the - master, `"indicator"` (constraint gated by the binary) or `"bigm"` - (rows relaxed by `M_value * (1 - z)`, tighter for solvers that - cannot strengthen indicators, e.g. with presolve disabled). -- `max_slack = 1e3`: upper bound of each OA cut slack. -- `oa_penalty = 1e3`: objective penalty per unit of cut slack. -- `use_nlpf = true`: solve a slacked feasibility NLP when the primary - NLP is infeasible, so its point still seeds OA cuts. -- `convergence_tol = 1e-6`: relative incumbent/bound gap tolerance. -- `slack_tol = 1e-4`: total cut slack tolerance for convergence. -- `iteration_time_limit = Inf`: seconds allotted to the LOA loop. -- `time_limit = 3600.0`: overall seconds budget. +`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 = MOIDisjunction()) +``` """ mutable struct Optimizer <: MOI.AbstractOptimizer nlp_solver::Any mip_solver::Any cache::_Cache - options::Dict{String, Any} + algorithm::Union{Nothing, AbstractAlgorithm} + time_limit_sec::Union{Nothing, Float64} silent::Bool # results (filled by MOI.optimize!) termination_status::MOI.TerminationStatusCode @@ -60,24 +161,15 @@ mutable struct Optimizer <: MOI.AbstractOptimizer solve_time::Float64 end -function Optimizer(; - nlp_solver = nothing, - mip_solver = nlp_solver, - kwargs... - ) - options = copy(_DEFAULT_OPTIONS) - for (key, value) in kwargs - haskey(options, string(key)) || throw(ArgumentError( - "Unknown option `$key`.")) - options[string(key)] = value - end +function Optimizer(nlp_solver, mip_solver = nlp_solver) return Optimizer(nlp_solver, mip_solver, MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), - options, false, MOI.OPTIMIZE_NOT_CALLED, MOI.NO_SOLUTION, + nothing, 3600.0, false, MOI.OPTIMIZE_NOT_CALLED, MOI.NO_SOLUTION, Dict{MOI.VariableIndex, Float64}(), NaN, nothing, NaN, "", NaN) end -_option(model::Optimizer, name::String) = model.options[name] +_algorithm(model::Optimizer) = + something(model.algorithm, _default(Algorithm())) MOI.get(::Optimizer, ::MOI.SolverName) = "DisjunctiveAlgorithms" MOI.get(::Optimizer, ::MOI.SolverVersion) = "0.1.0" @@ -111,6 +203,8 @@ 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 ################################################################################ @@ -311,28 +405,33 @@ function MOI.set( ::MOI.TimeLimitSec, value::Union{Nothing, Real} ) - model.options["time_limit"] = value === nothing ? Inf : Float64(value) + model.time_limit_sec = value === nothing ? nothing : Float64(value) return end -function MOI.get(model::Optimizer, ::MOI.TimeLimitSec) - limit = model.options["time_limit"] - return isfinite(limit) ? limit : nothing +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::MOI.RawOptimizerAttribute) - return haskey(model.options, attr.name) +function MOI.supports(model::Optimizer, attr::AbstractAlgorithmAttribute) + return MOI.supports(_algorithm(model), attr) end -function MOI.set(model::Optimizer, attr::MOI.RawOptimizerAttribute, value) - MOI.supports(model, attr) || throw(MOI.UnsupportedAttribute(attr)) - model.options[attr.name] = value +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::MOI.RawOptimizerAttribute) - MOI.supports(model, attr) || throw(MOI.UnsupportedAttribute(attr)) - return model.options[attr.name] +function MOI.get(model::Optimizer, attr::AbstractAlgorithmAttribute) + return MOI.get(_algorithm(model), attr) end ################################################################################ diff --git a/test/integration.jl b/test/integration.jl index 41067e0..f7732b8 100644 --- a/test/integration.jl +++ b/test/integration.jl @@ -1,8 +1,7 @@ using DisjunctiveProgramming, HiGHS, Ipopt, InfiniteOpt function _optimizer_factory() - return () -> DA.Optimizer(nlp_solver = Ipopt.Optimizer, - mip_solver = HiGHS.Optimizer) + return () -> DA.Optimizer(Ipopt.Optimizer, HiGHS.Optimizer) end # The same GDP solved through a BigM reformulation and through the diff --git a/test/loa.jl b/test/loa.jl index df9a662..2ceed65 100644 --- a/test/loa.jl +++ b/test/loa.jl @@ -1,9 +1,9 @@ using JuMP import HiGHS, Ipopt -function _loa_optimizer(; kwargs...) - return () -> DA.Optimizer(; nlp_solver = Ipopt.Optimizer, - mip_solver = HiGHS.Optimizer, kwargs...) +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 @@ -37,9 +37,8 @@ function test_row_function_constants() 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)]])) + @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 @@ -195,10 +194,10 @@ function test_complement_indicator() @test value(z) ≈ 0.0 atol = 1e-5 end -# `use_nlpf = false` still solves by enumeration when the seeds are -# feasible. +# `UseNLPF() => false` still solves by enumeration when the seeds +# are feasible. function test_nlpf_disabled() - model = Model(_loa_optimizer(use_nlpf = false)) + model = Model(_loa_optimizer(DA.UseNLPF() => false)) set_silent(model) @variable(model, 0 <= x <= 10) @variable(model, z[1:2], Bin) @@ -212,8 +211,9 @@ end # Big-M master gating solves the same instances as indicator gating # (interval split included). -function test_bigm_master_gating() - model = Model(_loa_optimizer(master_gating = "bigm", M_value = 100.0)) +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) @@ -226,11 +226,11 @@ function test_bigm_master_gating() @test objective_value(model) ≈ 7.0 atol = 1e-4 end -# Integer-typed options convert at their use sites, including the -# Bool read of use_nlpf on the infeasible-seed path. +# Integer-typed attribute values convert on set, including the Bool +# `UseNLPF` on the infeasible-seed path. function test_integer_options() - model = Model(_loa_optimizer(use_nlpf = 0, M_value = 10^9, - time_limit = 3600)) + 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) @@ -246,8 +246,7 @@ end # in the subproblem, which the nlp_solver must then handle (HiGHS # both roles here since the model is linear). function test_non_indicator_binary() - factory = () -> DA.Optimizer(nlp_solver = HiGHS.Optimizer, - mip_solver = HiGHS.Optimizer) + factory = () -> DA.Optimizer(HiGHS.Optimizer) model = Model(factory) set_silent(model) @variable(model, 0 <= x <= 10) @@ -281,7 +280,7 @@ end # A zero time limit exits before any solve. function test_time_limit() - model = Model(_loa_optimizer(time_limit = 0.0)) + model = Model(_loa_optimizer(MOI.TimeLimitSec() => 0.0)) set_silent(model) @variable(model, 0 <= x <= 10) @variable(model, z[1:2], Bin) @@ -400,10 +399,11 @@ function test_nonlinear_exp_global() @test value(z[2]) ≈ 1.0 atol = 1e-5 end -# `use_nlpf = false` with an infeasible combination: the solve keeps -# only the no-good cut and still finishes from the other disjunct. +# `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(use_nlpf = false)) + model = Model(_loa_optimizer(DA.UseNLPF() => false)) set_silent(model) @variable(model, 0 <= x <= 10) @variable(model, z[1:2], Bin) @@ -473,7 +473,8 @@ end # Zero iteration budgets exit before any solve, without an incumbent. function test_iteration_limit_no_incumbent() - model = Model(_loa_optimizer(set_cover_max_iter = 0, max_iter = 0)) + model = Model(_loa_optimizer(DA.SetCoverIterationLimit() => 0, + DA.NumIterationLimit() => 0)) set_silent(model) @variable(model, 0 <= x <= 10) @variable(model, z[1:2], Bin) @@ -485,10 +486,10 @@ function test_iteration_limit_no_incumbent() @test result_count(model) == 0 end -# `max_iter = 0` keeps the set-covering incumbent but produces no -# master bound. +# `NumIterationLimit() => 0` keeps the set-covering incumbent but +# produces no master bound. function test_iteration_limit_with_incumbent() - model = Model(_loa_optimizer(max_iter = 0)) + model = Model(_loa_optimizer(DA.NumIterationLimit() => 0)) set_silent(model) @variable(model, 0 <= x <= 10) @variable(model, z[1:2], Bin) @@ -506,7 +507,7 @@ end # An unbounded master (no OA cuts yet bound alpha_oa) surfaces its # status instead of looping. function test_master_abnormal_status() - model = Model(_loa_optimizer(set_cover_max_iter = 0)) + model = Model(_loa_optimizer(DA.SetCoverIterationLimit() => 0)) set_silent(model) @variable(model, 0 <= x <= 10) @variable(model, z[1:2], Bin) @@ -519,115 +520,11 @@ function test_master_abnormal_status() @test occursin("master solve finished", raw_status(model)) end -################################################################################ -# MOCK SOLVER -################################################################################ -# Delegates every MOI call to a wrapped optimizer, but sleeps -# `sleep_time` seconds in each solve and, from solve `fail_from` on, -# skips the inner solve and reports `fail_status` with no solution. -# Deterministic triggers for the deadline and abnormal-master paths. -mutable struct MockSolver <: MOI.AbstractOptimizer - inner::MOI.AbstractOptimizer - sleep_time::Float64 - fail_from::Int - fail_status::MOI.TerminationStatusCode - solves::Int - failing::Bool -end - -function MockSolver( - factory; - sleep_time::Float64 = 0.0, - fail_from::Int = typemax(Int), - fail_status::MOI.TerminationStatusCode = MOI.NODE_LIMIT - ) - return MockSolver(MOI.instantiate(factory), sleep_time, fail_from, - fail_status, 0, false) -end - -function MOI.optimize!(model::MockSolver) - model.solves += 1 - model.sleep_time > 0 && sleep(model.sleep_time) - model.failing = model.solves >= model.fail_from - model.failing || MOI.optimize!(model.inner) - return -end - -const _WrappedAttr = Union{MOI.AbstractModelAttribute, - MOI.AbstractOptimizerAttribute} -const _WrappedIndexAttr = Union{MOI.AbstractVariableAttribute, - MOI.AbstractConstraintAttribute} -const _WrappedIndex = Union{MOI.VariableIndex, MOI.ConstraintIndex} - -function MOI.get(model::MockSolver, attr::_WrappedAttr) - if model.failing - attr isa MOI.TerminationStatus && return model.fail_status - attr isa MOI.PrimalStatus && return MOI.NO_SOLUTION - end - return MOI.get(model.inner, attr) -end - -MOI.is_empty(model::MockSolver) = MOI.is_empty(model.inner) -MOI.empty!(model::MockSolver) = MOI.empty!(model.inner) -MOI.supports_incremental_interface(::MockSolver) = true -MOI.copy_to(model::MockSolver, src::MOI.ModelLike) = - MOI.copy_to(model.inner, src) -MOI.add_variable(model::MockSolver) = MOI.add_variable(model.inner) -MOI.delete(model::MockSolver, index) = MOI.delete(model.inner, index) -MOI.is_valid(model::MockSolver, index) = MOI.is_valid(model.inner, index) - -function MOI.add_constraint( - model::MockSolver, - func::MOI.AbstractFunction, - set::MOI.AbstractSet - ) - return MOI.add_constraint(model.inner, func, set) -end - -function MOI.supports_constraint( - model::MockSolver, - F::Type{<:MOI.AbstractFunction}, - S::Type{<:MOI.AbstractSet} - ) - return MOI.supports_constraint(model.inner, F, S) -end - -MOI.supports(model::MockSolver, attr::_WrappedAttr) = - MOI.supports(model.inner, attr) - -MOI.set(model::MockSolver, attr::_WrappedAttr, value) = - MOI.set(model.inner, attr, value) - -function MOI.supports( - model::MockSolver, - attr::_WrappedIndexAttr, - I::Type{<:_WrappedIndex} - ) - return MOI.supports(model.inner, attr, I) -end - -function MOI.get( - model::MockSolver, - attr::_WrappedIndexAttr, - index::_WrappedIndex - ) - return MOI.get(model.inner, attr, index) -end - -function MOI.set( - model::MockSolver, - attr::_WrappedIndexAttr, - index::_WrappedIndex, - value - ) - return MOI.set(model.inner, attr, index, value) -end - function _mock_time_limit_model(limit::Float64; kwargs...) - factory = () -> DA.Optimizer( - nlp_solver = () -> MockSolver(Ipopt.Optimizer; kwargs...), - mip_solver = HiGHS.Optimizer, - iteration_time_limit = limit) + 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) @@ -638,6 +535,26 @@ function _mock_time_limit_model(limit::Float64; kwargs...) 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 @@ -657,9 +574,8 @@ end # The master finishes abnormally after an incumbent exists: the mock # master reports a node limit on its second solve. function test_master_abnormal_status_with_incumbent() - factory = () -> DA.Optimizer( - nlp_solver = Ipopt.Optimizer, - mip_solver = () -> MockSolver(HiGHS.Optimizer, fail_from = 2)) + factory = () -> DA.Optimizer(Ipopt.Optimizer, + () -> MockSolver(HiGHS.Optimizer, fail_from = 2)) model = Model(factory) set_silent(model) @variable(model, 0 <= x <= 10) @@ -786,7 +702,7 @@ end test_nlpf_slacked_rows() test_demoted_affine_global() test_vector_of_variables_disjunction() - test_bigm_master_gating() + test_bigm_master_reformulation() test_integer_options() test_non_indicator_binary() test_infeasible_no_incumbent() @@ -794,6 +710,7 @@ end test_iteration_limit_no_incumbent() test_iteration_limit_with_incumbent() test_master_abnormal_status() + test_strict_constant_solvers() test_time_limit_with_incumbent() test_master_abnormal_status_with_incumbent() test_nonconvex_never_optimal() diff --git a/test/mock_optimizer.jl b/test/mock_optimizer.jl new file mode 100644 index 0000000..3897989 --- /dev/null +++ b/test/mock_optimizer.jl @@ -0,0 +1,113 @@ +################################################################################ +# MOCK SOLVER +################################################################################ +# Delegates every MOI call to a wrapped optimizer, but sleeps +# `sleep_time` seconds in each solve and, from solve `fail_from` on, +# skips the inner solve and reports `fail_status` with no solution. +# Deterministic triggers for the deadline and abnormal-master paths. +# `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 + 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 + ) + return MockSolver(MOI.instantiate(factory), sleep_time, fail_from, + fail_status, strict_constants, 0, false) +end + +function MOI.optimize!(model::MockSolver) + model.solves += 1 + model.sleep_time > 0 && sleep(model.sleep_time) + model.failing = model.solves >= model.fail_from + model.failing || MOI.optimize!(model.inner) + return +end + +const _WrappedAttr = Union{MOI.AbstractModelAttribute, + MOI.AbstractOptimizerAttribute} +const _WrappedIndexAttr = Union{MOI.AbstractVariableAttribute, + MOI.AbstractConstraintAttribute} +const _WrappedIndex = Union{MOI.VariableIndex, MOI.ConstraintIndex} + +function MOI.get(model::MockSolver, attr::_WrappedAttr) + if model.failing + attr isa MOI.TerminationStatus && return model.fail_status + attr isa MOI.PrimalStatus && return MOI.NO_SOLUTION + end + return MOI.get(model.inner, attr) +end + +MOI.is_empty(model::MockSolver) = MOI.is_empty(model.inner) +MOI.empty!(model::MockSolver) = MOI.empty!(model.inner) +MOI.supports_incremental_interface(::MockSolver) = true +MOI.copy_to(model::MockSolver, src::MOI.ModelLike) = + MOI.copy_to(model.inner, src) +MOI.add_variable(model::MockSolver) = MOI.add_variable(model.inner) +MOI.delete(model::MockSolver, index) = MOI.delete(model.inner, index) +MOI.is_valid(model::MockSolver, index) = MOI.is_valid(model.inner, index) + +function MOI.add_constraint( + model::MockSolver, + func::MOI.AbstractFunction, + set::MOI.AbstractSet + ) + 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 index 25fd8b1..68b1db2 100644 --- a/test/moi.jl +++ b/test/moi.jl @@ -7,10 +7,9 @@ import Ipopt @testset "MOI contract (model API + attributes)" begin optimizer = MOI.instantiate( - () -> DA.Optimizer( - nlp_solver = Ipopt.Optimizer, - mip_solver = HiGHS.Optimizer, - time_limit = 20.0, + MOI.OptimizerWithAttributes( + () -> DA.Optimizer(Ipopt.Optimizer, HiGHS.Optimizer), + MOI.TimeLimitSec() => 20.0, ); with_cache_type = Float64, with_bridge_type = Float64, diff --git a/test/optimizer.jl b/test/optimizer.jl index e5f004c..5764c61 100644 --- a/test/optimizer.jl +++ b/test/optimizer.jl @@ -18,7 +18,7 @@ function _build_toy_disjunction(model) end function test_optimizer_scaffold() - optimizer = DA.Optimizer() + optimizer = DA.Optimizer(nothing) @test MOI.is_empty(optimizer) @test MOI.get(optimizer, MOI.SolverName()) == "DisjunctiveAlgorithms" @test MOI.get(optimizer, MOI.TerminationStatus()) == @@ -29,27 +29,32 @@ function test_optimizer_scaffold() end function test_optimizer_options() - optimizer = DA.Optimizer(max_iter = 25) - @test MOI.get(optimizer, MOI.RawOptimizerAttribute("max_iter")) == 25 - MOI.set(optimizer, MOI.RawOptimizerAttribute("max_iter"), 5) - @test MOI.get(optimizer, MOI.RawOptimizerAttribute("max_iter")) == 5 - @test MOI.get(optimizer, MOI.RawOptimizerAttribute("M_value")) == 1e9 - bogus = MOI.RawOptimizerAttribute("bogus") - @test !MOI.supports(optimizer, bogus) - @test_throws MOI.UnsupportedAttribute MOI.get(optimizer, bogus) - @test_throws MOI.UnsupportedAttribute MOI.set(optimizer, bogus, 1) - @test_throws ArgumentError DA.Optimizer(bogus = 1) + 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 - @test MOI.get(optimizer, MOI.RawOptimizerAttribute("time_limit")) == 100.0 MOI.set(optimizer, MOI.TimeLimitSec(), nothing) @test MOI.get(optimizer, MOI.TimeLimitSec()) === nothing end function test_model_building() - optimizer = DA.Optimizer() + 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 @@ -74,7 +79,7 @@ function test_copy_to() 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() + 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 @@ -87,7 +92,7 @@ end # The incremental MOI surface forwards to the cache: names, starts, # constrained variables, and deletion. function test_moi_forwarding() - optimizer = DA.Optimizer() + optimizer = DA.Optimizer(nothing) @test MOI.supports_incremental_interface(optimizer) x = MOI.add_variable(optimizer) MOI.set(optimizer, MOI.VariableName(), x, "x") @@ -109,7 +114,7 @@ end # A constraint type the partition cannot place errors at solve time. function test_unsupported_constraint_type() - optimizer = DA.Optimizer() + optimizer = DA.Optimizer(nothing) x = MOI.add_variable(optimizer) func = MOI.VectorAffineFunction( [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x))], [0.0]) diff --git a/test/runtests.jl b/test/runtests.jl index f659abc..50fb40b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,6 +4,7 @@ using DisjunctiveAlgorithms const DA = DisjunctiveAlgorithms include("optimizer.jl") +include("mock_optimizer.jl") include("loa.jl") include("moi.jl") include("integration.jl") From 24db6803b9bdda03780942e69bbe45e4a80076eb Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Mon, 17 Aug 2026 12:50:36 -0400 Subject: [PATCH 06/18] Checking inner solvers for support --- README.md | 7 ++- src/algorithms/LOA.jl | 1 + src/nlp.jl | 2 +- src/optimizer.jl | 4 +- src/problem.jl | 99 +++++++++++++++++++++++++++++++++++++++++++ test/integration.jl | 10 ++--- test/loa.jl | 39 +++++++++++++++++ 7 files changed, 150 insertions(+), 12 deletions(-) diff --git a/README.md b/README.md index 1d82799..c2b0b16 100644 --- a/README.md +++ b/README.md @@ -24,17 +24,16 @@ options selected through optimizer attributes. using DisjunctiveProgramming, DisjunctiveAlgorithms, HiGHS, Ipopt import DisjunctiveAlgorithms as DA -model = GDPModel(() -> DA.Optimizer(nlp_solver = Ipopt.Optimizer, - mip_solver = HiGHS.Optimizer)) +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 = MOIDisjunction()) +optimize!(model, gdp_method = Direct()) ``` -`MOIDisjunction()` lowers each disjunction to a single +`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/src/algorithms/LOA.jl b/src/algorithms/LOA.jl index 05089e9..585454b 100644 --- a/src/algorithms/LOA.jl +++ b/src/algorithms/LOA.jl @@ -148,6 +148,7 @@ function _optimize!(algorithm::LOA, model::Optimizer) t_start = time() _reset_results(model) problem = _build_problem(model) + _check_inner_support(model, problem) master = _build_master(model, problem) subproblem = _build_subproblem(model, problem) linearizer = _Linearizer() diff --git a/src/nlp.jl b/src/nlp.jl index 0aded0b..c96f8bf 100644 --- a/src/nlp.jl +++ b/src/nlp.jl @@ -160,7 +160,7 @@ function _solve_nlpf( _cap_remaining_time(nlp, deadline) MOI.optimize!(nlp) # Use the primal only at a genuine feasible point; a solver can - # report values at a nonfeasible/NaN primal that poisons the cut. + # 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), diff --git a/src/optimizer.jl b/src/optimizer.jl index 9d321b7..e5a92cc 100644 --- a/src/optimizer.jl +++ b/src/optimizer.jl @@ -140,7 +140,7 @@ attributes it documents. A new optimizer starts with a 3600 second julia> model = GDPModel(() -> DisjunctiveAlgorithms.Optimizer( Ipopt.Optimizer, HiGHS.Optimizer)); -julia> optimize!(model, gdp_method = MOIDisjunction()) +julia> optimize!(model, gdp_method = Direct()) ``` """ mutable struct Optimizer <: MOI.AbstractOptimizer @@ -329,7 +329,7 @@ function MOI.supports( end # The LOA loop consumes warm starts, so accept them even though the -# honest inner model does not store them (the fallback does). +# inner model does not store them (the fallback does). function MOI.supports( ::Optimizer, ::MOI.VariablePrimalStart, diff --git a/src/problem.jl b/src/problem.jl index 3f7b959..1f96938 100644 --- a/src/problem.jl +++ b/src/problem.jl @@ -161,3 +161,102 @@ 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. +function _check_inner_support(model::Optimizer, problem::_Problem) + mip = _instantiate(model.mip_solver, "mip_solver") + _check_support(mip, "mip_solver", "master problem", + _master_constraint_types(model, problem)) + nlp = _instantiate(model.nlp_solver, "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 index f7732b8..c628981 100644 --- a/test/integration.jl +++ b/test/integration.jl @@ -16,7 +16,7 @@ function test_lowering_solve_linear() @constraint(model, x <= 7, Disjunct(Y[2])) @disjunction(model, Y) @objective(model, Max, x) - optimize!(model, gdp_method = MOIDisjunction()) + 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 @@ -45,7 +45,7 @@ function test_lowering_solve_nonlinear() @constraint(model, x^2 == 64, Disjunct(Y[2])) @disjunction(model, Y) @objective(model, Max, x) - optimize!(model, gdp_method = MOIDisjunction()) + 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]) @@ -65,7 +65,7 @@ function test_lowering_solve_nested() @disjunction(model, W, Disjunct(Y[2])) @disjunction(model, Y) @objective(model, Max, x) - optimize!(model, gdp_method = MOIDisjunction()) + 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]) @@ -84,7 +84,7 @@ function test_lowering_infinite() @constraint(model, x >= 5, Disjunct(Y[2])) @disjunction(model, Y) @objective(model, Max, integral(x, t)) - optimize!(model, gdp_method = MOIDisjunction()) + 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) @@ -106,7 +106,7 @@ function test_lowering_infinite_nested() @disjunction(model, W, Disjunct(Y[2])) @disjunction(model, Y) @objective(model, Max, integral(x, t)) - optimize!(model, gdp_method = MOIDisjunction()) + optimize!(model, gdp_method = Direct()) @test termination_status(model) == MOI.LOCALLY_SOLVED @test objective_value(model) ≈ 10.0 atol = 1e-4 end diff --git a/test/loa.jl b/test/loa.jl index 2ceed65..42c4d0f 100644 --- a/test/loa.jl +++ b/test/loa.jl @@ -683,6 +683,44 @@ function test_nested_disjunction_vacuous() @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 + @testset "LOA loop" begin test_linear_disjunction() test_row_function_constants() @@ -718,6 +756,7 @@ end test_linear_interval_disjunct() test_reoptimize_resets_results() test_bridged_vector_constraint() + test_inner_solver_support_check() end @testset "LOA units" begin From 890b547bb2b12f6773c47a71dbf626d5e3bad070 Mon Sep 17 00:00:00 2001 From: dnguyen227 Date: Sun, 23 Aug 2026 14:41:42 -0400 Subject: [PATCH 07/18] solve statuses and set-time attribute validation --- Project.toml | 2 +- README.md | 14 ++++++ src/DisjunctiveAlgorithms.jl | 5 +-- src/algorithms/LOA.jl | 82 +++++++++++++++++++++++++++++++++--- src/master.jl | 2 - src/nlp.jl | 7 +-- src/optimizer.jl | 16 +++++-- src/problem.jl | 25 ++++++++--- test/loa.jl | 80 +++++++++++++++++++++++++++++++++++ test/mock_optimizer.jl | 6 ++- test/optimizer.jl | 20 +++++++++ 11 files changed, 232 insertions(+), 27 deletions(-) diff --git a/Project.toml b/Project.toml index ff640ed..7bd9b24 100644 --- a/Project.toml +++ b/Project.toml @@ -8,7 +8,7 @@ DisjunctiveProgramming = "0d27d021-0159-4c7d-b4a7-9ccb5d9366cf" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" [compat] -DisjunctiveProgramming = "0.6" +DisjunctiveProgramming = "0.6.2" MathOptInterface = "1" julia = "1.10" diff --git a/README.md b/README.md index c2b0b16..5b2300c 100644 --- a/README.md +++ b/README.md @@ -18,6 +18,20 @@ The design follows one `Optimizer` that wraps inner solvers, with the algorithm and its options selected through optimizer attributes. +## Installation + +DisjunctiveAlgorithms requires DisjunctiveProgramming v0.6.2, which is +not registered yet: it lives on the `gdp_optimizer` branch (it adds +`DisjunctionSet` and the `Direct()` reformulation). Until that release +ships, 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 diff --git a/src/DisjunctiveAlgorithms.jl b/src/DisjunctiveAlgorithms.jl index 0985c40..b25fdc3 100644 --- a/src/DisjunctiveAlgorithms.jl +++ b/src/DisjunctiveAlgorithms.jl @@ -1,9 +1,8 @@ module DisjunctiveAlgorithms import MathOptInterface as MOI -import DisjunctiveProgramming: DisjunctionSet, num_disjuncts, - activation_index, indicator_indices, row_indices, - _SupportedInnerSet +import DisjunctiveProgramming: DisjunctionSet, activation_index, + indicator_indices, row_indices, SupportedInnerSet include("optimizer.jl") include("problem.jl") diff --git a/src/algorithms/LOA.jl b/src/algorithms/LOA.jl index 585454b..8575814 100644 --- a/src/algorithms/LOA.jl +++ b/src/algorithms/LOA.jl @@ -43,6 +43,23 @@ 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(::OASlack, value) + value > 0 || + error("`OASlack` must be positive (got `$value`).") + return nothing +end + # An unset field (`nothing`) reads back as the attribute default. for (attr, field) in ( (NumIterationLimit, :num_iteration_limit), @@ -58,7 +75,8 @@ for (attr, field) in ( ) @eval begin MOI.supports(::LOA, ::$attr) = true - function MOI.set(algorithm::LOA, ::$attr, value) + function MOI.set(algorithm::LOA, attr::$attr, value) + _validate(attr, value) algorithm.$field = value return end @@ -75,6 +93,14 @@ _is_better(sense::MOI.OptimizationSense, 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, @@ -163,9 +189,20 @@ function _optimize!(algorithm::LOA, model::Optimizer) 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 -> begin + if !result.feasible && _nlp_unbounded(result.status) + unbounded = true + return + end + result.feasible || _nlp_infeasible(result.status) || + (num_unresolved += 1) _avoid_combination(master, result.combination) _add_oa_cuts(model, problem, master, linearizer, result) if result.feasible && @@ -204,6 +241,7 @@ function _optimize!(algorithm::LOA, model::Optimizer) result = _solve_nlp(model, problem, subproblem, combination, warm_start(); deadline = loop_deadline) process_result(result) + unbounded && break # covered only once active in a feasible NLP; infeasible # combinations just leave their no-good cut if result.feasible @@ -217,7 +255,7 @@ function _optimize!(algorithm::LOA, model::Optimizer) end # main loop: alpha_oa gives the bound, the NLP the incumbent - if master_status === nothing + if master_status === nothing && !unbounded for _ in 1:MOI.get(algorithm, NumIterationLimit()) time() < loop_deadline || break _cap_remaining_time(master.model, loop_deadline) @@ -245,11 +283,12 @@ function _optimize!(algorithm::LOA, model::Optimizer) result = _solve_nlp(model, problem, subproblem, combination, warm_start(); deadline = loop_deadline) process_result(result) + unbounded && break end end _store_results(model, sense, best_objective, best_result, master_bound, - master_status, converged, loop_deadline) + master_status, converged, unbounded, num_unresolved, loop_deadline) model.solve_time = time() - t_start return end @@ -267,21 +306,43 @@ function _store_results( 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 + 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 " * @@ -294,21 +355,28 @@ function _store_results( 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 master_status == MOI.INFEASIBLE + 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" : (master_status == MOI.INFEASIBLE ? - "combinations exhausted" : "limit hit") + 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)." diff --git a/src/master.jl b/src/master.jl index a909839..5bf8dbb 100644 --- a/src/master.jl +++ b/src/master.jl @@ -104,8 +104,6 @@ function _add_gated_rows( model::Optimizer ) gating = MOI.get(model, MasterReformulation()) - gating in ("indicator", "bigm") || - error("Unknown `MasterReformulation` value `$gating`.") activate = disjunct.active_value ? MOI.ACTIVATE_ON_ONE : MOI.ACTIVATE_ON_ZERO binary = variable_map[disjunct.binary] diff --git a/src/nlp.jl b/src/nlp.jl index c96f8bf..44a52bf 100644 --- a/src/nlp.jl +++ b/src/nlp.jl @@ -86,19 +86,20 @@ function _solve_nlp( _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) + 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 + result === nothing || return (; result..., status = status) end return (combination = combination, - point = nothing, objective = Inf, feasible = false) + point = nothing, objective = Inf, feasible = false, status = status) end ################################################################################ diff --git a/src/optimizer.jl b/src/optimizer.jl index e5a92cc..aa463d3 100644 --- a/src/optimizer.jl +++ b/src/optimizer.jl @@ -172,7 +172,7 @@ _algorithm(model::Optimizer) = something(model.algorithm, _default(Algorithm())) MOI.get(::Optimizer, ::MOI.SolverName) = "DisjunctiveAlgorithms" -MOI.get(::Optimizer, ::MOI.SolverVersion) = "0.1.0" +MOI.get(::Optimizer, ::MOI.SolverVersion) = string(pkgversion(@__MODULE__)) MOI.is_empty(model::Optimizer) = MOI.is_empty(model.cache) @@ -260,7 +260,7 @@ end function MOI.supports_constraint( ::Optimizer, ::Type{MOI.VariableIndex}, - ::Type{<:Union{_SupportedInnerSet, MOI.ZeroOne, MOI.Integer}} + ::Type{<:Union{SupportedInnerSet, MOI.ZeroOne, MOI.Integer}} ) return true end @@ -268,7 +268,7 @@ end function MOI.supports_constraint( ::Optimizer, ::Type{<:_ScalarFunction}, - ::Type{<:_SupportedInnerSet} + ::Type{<:SupportedInnerSet} ) return true end @@ -310,6 +310,16 @@ end 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) diff --git a/src/problem.jl b/src/problem.jl index 1f96938..f7e9358 100644 --- a/src/problem.jl +++ b/src/problem.jl @@ -64,6 +64,21 @@ 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) @@ -75,9 +90,7 @@ function _activation_binary(activation::MOI.ScalarAffineFunction{Float64}) return term.variable, false end end - return error("Unsupported indicator expression `$activation`: each " * - "`DisjunctionSet` indicator must be a binary variable `z` or " * - "its complement `1 - z`.") + return _indicator_error(activation) end function _parse_disjunction( @@ -86,11 +99,11 @@ function _parse_disjunction( ) where {F} set = MOI.get(cache, MOI.ConstraintSet(), ci) rows = _scalarize(MOI.get(cache, MOI.ConstraintFunction(), ci)) - disjunction_activation = _to_affine( - _demote(rows[activation_index(set)])) + disjunction_activation = _activation_affine( + rows[activation_index(set)]) disjuncts = _Disjunct[] for (i, j) in enumerate(indicator_indices(set)) - activation = _to_affine(_demote(rows[j])) + activation = _activation_affine(rows[j]) binary, active_value = _activation_binary(activation) zero_one = MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}( binary.value) diff --git a/test/loa.jl b/test/loa.jl index 42c4d0f..3c7f423 100644 --- a/test/loa.jl +++ b/test/loa.jl @@ -520,6 +520,82 @@ function test_master_abnormal_status() @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...), @@ -748,6 +824,10 @@ end 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() diff --git a/test/mock_optimizer.jl b/test/mock_optimizer.jl index 3897989..b21c826 100644 --- a/test/mock_optimizer.jl +++ b/test/mock_optimizer.jl @@ -56,8 +56,10 @@ MOI.supports_incremental_interface(::MockSolver) = true MOI.copy_to(model::MockSolver, src::MOI.ModelLike) = MOI.copy_to(model.inner, src) MOI.add_variable(model::MockSolver) = MOI.add_variable(model.inner) -MOI.delete(model::MockSolver, index) = MOI.delete(model.inner, index) -MOI.is_valid(model::MockSolver, index) = MOI.is_valid(model.inner, index) +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, diff --git a/test/optimizer.jl b/test/optimizer.jl index 5764c61..2fd510f 100644 --- a/test/optimizer.jl +++ b/test/optimizer.jl @@ -122,6 +122,25 @@ function test_unsupported_constraint_type() @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() @@ -129,4 +148,5 @@ end test_copy_to() test_moi_forwarding() test_unsupported_constraint_type() + test_objective_function_support() end From a579c04a7a241df24d78c2c53c5ffef5bcb632e1 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Tue, 8 Sep 2026 14:36:42 -0400 Subject: [PATCH 08/18] Remove try/catch demotion and dead factory guard --- src/master.jl | 6 ++---- src/nlp.jl | 4 ++-- src/problem.jl | 37 ++++++++++++++++++++++++------------- 3 files changed, 28 insertions(+), 19 deletions(-) diff --git a/src/master.jl b/src/master.jl index 5bf8dbb..e2146a0 100644 --- a/src/master.jl +++ b/src/master.jl @@ -12,9 +12,7 @@ mutable struct _Master oa_objective::MOI.ScalarAffineFunction{Float64} end -function _instantiate(factory, name::String) - factory === nothing && - error("DisjunctiveAlgorithms requires the `$name` optimizer factory.") +function _instantiate(factory) solver = MOI.instantiate(factory; with_cache_type = Float64, with_bridge_type = Float64) MOI.set(solver, MOI.Silent(), true) @@ -37,7 +35,7 @@ function _cap_remaining_time(solver::MOI.ModelLike, deadline::Float64) end function _build_master(model::Optimizer, problem::_Problem) - mip = _instantiate(model.mip_solver, "mip_solver") + 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 diff --git a/src/nlp.jl b/src/nlp.jl index 44a52bf..aa09244 100644 --- a/src/nlp.jl +++ b/src/nlp.jl @@ -12,7 +12,7 @@ struct _Subproblem end function _build_subproblem(model::Optimizer, problem::_Problem) - nlp = _instantiate(model.nlp_solver, "nlp_solver") + nlp = _instantiate(model.nlp_solver) variable_map = Dict{MOI.VariableIndex, MOI.VariableIndex}( vi => MOI.add_variable(nlp) for vi in problem.variables) indicators = Set(problem.binaries) @@ -121,7 +121,7 @@ function _solve_nlpf( warm_start; deadline::Float64 = Inf ) - nlp = _instantiate(model.nlp_solver, "nlp_solver") + nlp = _instantiate(model.nlp_solver) variable_map = Dict{MOI.VariableIndex, MOI.VariableIndex}( vi => MOI.add_variable(nlp) for vi in problem.variables) u = MOI.add_variable(nlp) diff --git a/src/problem.jl b/src/problem.jl index f7e9358..401c496 100644 --- a/src/problem.jl +++ b/src/problem.jl @@ -43,22 +43,33 @@ end _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 + # demote rows the enclosing vector function promoted _demote(func::MOI.AbstractScalarFunction) = func function _demote(func::MOI.ScalarQuadraticFunction{Float64}) - return _try_convert(MOI.ScalarAffineFunction{Float64}, func) + return isempty(func.quadratic_terms) ? _to_affine(func) : func end function _demote(func::MOI.ScalarNonlinearFunction) - return _try_convert(MOI.ScalarAffineFunction{Float64}, - _try_convert(MOI.ScalarQuadraticFunction{Float64}, func)) -end - -function _try_convert(T::Type, func) - return try - convert(T, func) - catch - func - end + degree = _degree(func) + degree in (0, 1) && return _to_affine(func) + degree == 2 && return convert(MOI.ScalarQuadraticFunction{Float64}, func) + return func end _scalarize(func::MOI.AbstractVectorFunction) = @@ -261,10 +272,10 @@ end # Fail before any subproblem work when an inner solver cannot take the # constraint types routed to it, naming the solver and the type. function _check_inner_support(model::Optimizer, problem::_Problem) - mip = _instantiate(model.mip_solver, "mip_solver") + mip = _instantiate(model.mip_solver) _check_support(mip, "mip_solver", "master problem", _master_constraint_types(model, problem)) - nlp = _instantiate(model.nlp_solver, "nlp_solver") + nlp = _instantiate(model.nlp_solver) _check_support(nlp, "nlp_solver", "NLP subproblems", _nlp_constraint_types(model, problem)) F = typeof(problem.objective) From 1db01ef3a46cbc7f2397a95db9a69714e1f1d87c Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Thu, 27 Aug 2026 11:51:56 -0400 Subject: [PATCH 09/18] Add SubproblemMethod extension point and solver attribute guards --- src/algorithms/LOA.jl | 26 +++++++++++++++++++------- src/master.jl | 6 ++++-- src/nlp.jl | 6 +++++- src/optimizer.jl | 15 +++++++++++++++ src/problem.jl | 13 +++++++++++-- test/loa.jl | 43 +++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 97 insertions(+), 12 deletions(-) diff --git a/src/algorithms/LOA.jl b/src/algorithms/LOA.jl index 8575814..445b8b9 100644 --- a/src/algorithms/LOA.jl +++ b/src/algorithms/LOA.jl @@ -24,6 +24,7 @@ never `OPTIMAL`. - [`ConvergenceTolerance`](@ref) - [`SlackTolerance`](@ref) - [`IterationTimeLimit`](@ref) +- [`SubproblemMethod`](@ref) """ mutable struct LOA <: AbstractAlgorithm num_iteration_limit::Union{Nothing, Int} @@ -36,9 +37,10 @@ mutable struct LOA <: AbstractAlgorithm convergence_tolerance::Union{Nothing, Float64} slack_tolerance::Union{Nothing, Float64} iteration_time_limit::Union{Nothing, Float64} + subproblem_method::Any LOA() = new(nothing, nothing, nothing, nothing, nothing, nothing, - nothing, nothing, nothing, nothing) + nothing, nothing, nothing, nothing, nothing) end _default(::Algorithm) = LOA() @@ -86,6 +88,15 @@ for (attr, field) in ( end end +# `nothing` is this attribute's default, so it skips 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 + _worst_objective(sense::MOI.OptimizationSense) = sense == MOI.MAX_SENSE ? -Inf : Inf _is_better(sense::MOI.OptimizationSense, new, best) = @@ -174,9 +185,10 @@ function _optimize!(algorithm::LOA, model::Optimizer) t_start = time() _reset_results(model) problem = _build_problem(model) - _check_inner_support(model, problem) + method = MOI.get(algorithm, SubproblemMethod()) + _check_inner_support(method, model, problem) master = _build_master(model, problem) - subproblem = _build_subproblem(model, problem) + subproblem = _build_subproblem(method, model, problem) linearizer = _Linearizer() sense = problem.sense overall_deadline = t_start + something(model.time_limit_sec, Inf) @@ -238,8 +250,8 @@ function _optimize!(algorithm::LOA, model::Optimizer) master_status = status break end - result = _solve_nlp(model, problem, subproblem, combination, - warm_start(); deadline = loop_deadline) + result = _solve_nlp(method, model, problem, subproblem, + combination, warm_start(); deadline = loop_deadline) process_result(result) unbounded && break # covered only once active in a feasible NLP; infeasible @@ -280,8 +292,8 @@ function _optimize!(algorithm::LOA, model::Optimizer) end end combination = _extract_combination(problem, master) - result = _solve_nlp(model, problem, subproblem, combination, - warm_start(); deadline = loop_deadline) + result = _solve_nlp(method, model, problem, subproblem, + combination, warm_start(); deadline = loop_deadline) process_result(result) unbounded && break end diff --git a/src/master.jl b/src/master.jl index e2146a0..dbda8d8 100644 --- a/src/master.jl +++ b/src/master.jl @@ -15,7 +15,8 @@ end function _instantiate(factory) solver = MOI.instantiate(factory; with_cache_type = Float64, with_bridge_type = Float64) - MOI.set(solver, MOI.Silent(), true) + MOI.supports(solver, MOI.Silent()) && + MOI.set(solver, MOI.Silent(), true) return solver end @@ -30,7 +31,8 @@ end function _cap_remaining_time(solver::MOI.ModelLike, deadline::Float64) isfinite(deadline) || return - MOI.set(solver, MOI.TimeLimitSec(), max(0.0, deadline - time())) + MOI.supports(solver, MOI.TimeLimitSec()) && + MOI.set(solver, MOI.TimeLimitSec(), max(0.0, deadline - time())) return end diff --git a/src/nlp.jl b/src/nlp.jl index aa09244..1e068bc 100644 --- a/src/nlp.jl +++ b/src/nlp.jl @@ -1,6 +1,9 @@ ################################################################################ # 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 @@ -11,7 +14,7 @@ struct _Subproblem rows::Vector{MOI.ConstraintIndex} end -function _build_subproblem(model::Optimizer, problem::_Problem) +function _build_subproblem(::Nothing, model::Optimizer, problem::_Problem) nlp = _instantiate(model.nlp_solver) variable_map = Dict{MOI.VariableIndex, MOI.VariableIndex}( vi => MOI.add_variable(nlp) for vi in problem.variables) @@ -61,6 +64,7 @@ end # 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, diff --git a/src/optimizer.jl b/src/optimizer.jl index aa463d3..01003d2 100644 --- a/src/optimizer.jl +++ b/src/optimizer.jl @@ -119,6 +119,21 @@ budget is the standard `MOI.TimeLimitSec` attribute. struct IterationTimeLimit <: AbstractAlgorithmAttribute end _default(::IterationTimeLimit) = Inf +""" + 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 ################################################################################ diff --git a/src/problem.jl b/src/problem.jl index 401c496..d9ed868 100644 --- a/src/problem.jl +++ b/src/problem.jl @@ -270,11 +270,20 @@ function _check_support(solver, name::String, destination::String, required) end # Fail before any subproblem work when an inner solver cannot take the -# constraint types routed to it, naming the solver and the type. -function _check_inner_support(model::Optimizer, problem::_Problem) +# 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)) diff --git a/test/loa.jl b/test/loa.jl index 3c7f423..69ec8cc 100644 --- a/test/loa.jl +++ b/test/loa.jl @@ -797,6 +797,48 @@ function test_inner_solver_support_check() @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() @@ -837,6 +879,7 @@ end test_reoptimize_resets_results() test_bridged_vector_constraint() test_inner_solver_support_check() + test_subproblem_method_dispatch() end @testset "LOA units" begin From dafe768fa3a7c9909c1ff290edba4f99ec150c91 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Thu, 10 Sep 2026 09:54:38 -0400 Subject: [PATCH 10/18] Evaluate degree <= 2 operator trees into typed functions MOI's own `convert` is 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` threw there. `_polynomial` walks the tree instead, and `_to_affine` gains a constant method. --- src/problem.jl | 34 ++++++++++++++++++++++++++++++++-- test/loa.jl | 23 +++++++++++++++++++++++ 2 files changed, 55 insertions(+), 2 deletions(-) diff --git a/src/problem.jl b/src/problem.jl index d9ed868..a745ae5 100644 --- a/src/problem.jl +++ b/src/problem.jl @@ -37,6 +37,10 @@ _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 @@ -60,6 +64,32 @@ function _degree(func::MOI.ScalarNonlinearFunction) 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}) @@ -67,8 +97,8 @@ function _demote(func::MOI.ScalarQuadraticFunction{Float64}) end function _demote(func::MOI.ScalarNonlinearFunction) degree = _degree(func) - degree in (0, 1) && return _to_affine(func) - degree == 2 && return convert(MOI.ScalarQuadraticFunction{Float64}, func) + degree in (0, 1) && return _to_affine(_polynomial(func)) + degree == 2 && return _demote(_polynomial(func)) return func end diff --git a/test/loa.jl b/test/loa.jl index 69ec8cc..42e8bab 100644 --- a/test/loa.jl +++ b/test/loa.jl @@ -48,6 +48,28 @@ 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 + function test_quadratic_objective() model = Model(_loa_optimizer()) set_silent(model) @@ -844,6 +866,7 @@ end test_row_function_constants() test_nested_disjunction() test_nested_disjunction_vacuous() + test_quadratic_row_shapes() test_quadratic_objective() test_max_sense_linear() test_two_disjunctions() From 401dba844eec43120543b2c31980b34ff742b2e0 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Thu, 10 Sep 2026 13:49:28 -0400 Subject: [PATCH 11/18] Add multi-generation cuts and batched NLP subproblems `MultiGenerationSize` evaluates several indicator combinations per master solve. They come from the `mip_solver`'s solution pool where it exposes one, and otherwise from re-solving the master behind a no-good cut per combination already taken. `BatchedSubproblems` stacks them into one NLP, with the sequential subproblem as the fallback for a failed batch. Adds `MasterSolveCount` and `NLPSolveCount`. --- src/algorithms/LOA.jl | 76 ++++++++++--- src/master.jl | 6 +- src/nlp.jl | 244 +++++++++++++++++++++++++++++++++++++----- src/optimizer.jl | 41 ++++++- test/loa.jl | 108 +++++++++++++++++++ 5 files changed, 434 insertions(+), 41 deletions(-) diff --git a/src/algorithms/LOA.jl b/src/algorithms/LOA.jl index 445b8b9..2b3b582 100644 --- a/src/algorithms/LOA.jl +++ b/src/algorithms/LOA.jl @@ -24,6 +24,7 @@ never `OPTIMAL`. - [`ConvergenceTolerance`](@ref) - [`SlackTolerance`](@ref) - [`IterationTimeLimit`](@ref) +- [`MultiGenerationSize`](@ref) - [`SubproblemMethod`](@ref) """ mutable struct LOA <: AbstractAlgorithm @@ -37,10 +38,11 @@ mutable struct LOA <: AbstractAlgorithm convergence_tolerance::Union{Nothing, Float64} slack_tolerance::Union{Nothing, Float64} iteration_time_limit::Union{Nothing, Float64} + multi_generation_size::Union{Nothing, Int} subproblem_method::Any LOA() = new(nothing, nothing, nothing, nothing, nothing, nothing, - nothing, nothing, nothing, nothing, nothing) + nothing, nothing, nothing, nothing, nothing, nothing) end _default(::Algorithm) = LOA() @@ -62,6 +64,12 @@ function _validate(::OASlack, 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), @@ -74,6 +82,7 @@ for (attr, field) in ( (ConvergenceTolerance, :convergence_tolerance), (SlackTolerance, :slack_tolerance), (IterationTimeLimit, :iteration_time_limit), + (MultiGenerationSize, :multi_generation_size), ) @eval begin MOI.supports(::LOA, ::$attr) = true @@ -170,6 +179,44 @@ function _user_start_values(model::Optimizer, problem::_Problem) return isempty(point) ? nothing : (point = point,) end +function _solve_master(model::Optimizer, master::_Master, deadline::Float64) + _cap_remaining_time(master.model, deadline) + MOI.optimize!(master.model) + model.num_master_solves += 1 + return _solved_and_feasible(master.model) +end + +# The master's solution pool first (result indices past 1), then +# re-solves behind a no-good cut for whatever the pool did not supply. +# Those cuts are the ones `process_result` would add later, so the +# second return value tells the caller to skip them. +function _extract_combinations( + model::Optimizer, + problem::_Problem, + master::_Master, + count::Int, + deadline::Float64 + ) + combinations = [_extract_combination(problem, master)] + for index in 2:min(count, MOI.get(master.model, MOI.ResultCount())) + combination = _extract_combination(problem, master, index) + combination in combinations || push!(combinations, combination) + end + excluded = false + while length(combinations) < count && time() < deadline + if !excluded + foreach(c -> _avoid_combination(master, c), combinations) + excluded = true + end + _solve_master(model, master, deadline) || break + combination = _extract_combination(problem, master) + combination in combinations && break + push!(combinations, combination) + _avoid_combination(master, combination) + end + return combinations, excluded +end + function _set_master_objective(master::_Master, sense, objective) MOI.set(master.model, MOI.ObjectiveSense(), sense) MOI.set(master.model, @@ -208,14 +255,14 @@ function _optimize!(algorithm::LOA, model::Optimizer) # 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 -> begin + 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) - _avoid_combination(master, result.combination) + nogood && _avoid_combination(master, result.combination) _add_oa_cuts(model, problem, master, linearizer, result) if result.feasible && _is_better(sense, result.objective, best_objective) @@ -239,9 +286,7 @@ function _optimize!(algorithm::LOA, model::Optimizer) time() < loop_deadline || break _set_master_objective(master, MOI.MAX_SENSE, _cover_objective(master, cover, needs_cover, num_covered)) - _cap_remaining_time(master.model, loop_deadline) - MOI.optimize!(master.model) - solved = _solved_and_feasible(master.model) + solved = _solve_master(model, master, loop_deadline) # capture the status before the objective restore invalidates it status = MOI.get(master.model, MOI.TerminationStatus()) combination = solved ? _extract_combination(problem, master) : nothing @@ -252,6 +297,7 @@ function _optimize!(algorithm::LOA, model::Optimizer) end result = _solve_nlp(method, model, problem, subproblem, combination, warm_start(); deadline = loop_deadline) + model.num_nlp_solves += 1 process_result(result) unbounded && break # covered only once active in a feasible NLP; infeasible @@ -270,9 +316,7 @@ function _optimize!(algorithm::LOA, model::Optimizer) if master_status === nothing && !unbounded for _ in 1:MOI.get(algorithm, NumIterationLimit()) time() < loop_deadline || break - _cap_remaining_time(master.model, loop_deadline) - MOI.optimize!(master.model) - if !_solved_and_feasible(master.model) + if !_solve_master(model, master, loop_deadline) master_status = MOI.get(master.model, MOI.TerminationStatus()) break end @@ -291,10 +335,16 @@ function _optimize!(algorithm::LOA, model::Optimizer) break end end - combination = _extract_combination(problem, master) - result = _solve_nlp(method, model, problem, subproblem, - combination, warm_start(); deadline = loop_deadline) - process_result(result) + combinations, excluded = _extract_combinations(model, problem, + master, MOI.get(algorithm, MultiGenerationSize()), + loop_deadline) + results = _solve_nlps(method, model, problem, subproblem, + combinations, warm_start(); deadline = loop_deadline) + model.num_nlp_solves += length(results) + for result in results + process_result(result; nogood = !excluded) + unbounded && break + end unbounded && break end end diff --git a/src/master.jl b/src/master.jl index dbda8d8..f6f13bf 100644 --- a/src/master.jl +++ b/src/master.jl @@ -153,10 +153,10 @@ function _add_penalized_slack( end # round to Bool; MILP values are only within integer tolerance -function _extract_combination(problem::_Problem, master::_Master) +function _extract_combination(problem::_Problem, master::_Master, index = 1) return Dict{MOI.VariableIndex, Bool}( - binary => round(Bool, MOI.get(master.model, MOI.VariablePrimal(), - master.variable_map[binary])) + binary => round(Bool, MOI.get(master.model, + MOI.VariablePrimal(index), master.variable_map[binary])) for binary in problem.binaries) end diff --git a/src/nlp.jl b/src/nlp.jl index 1e068bc..0916505 100644 --- a/src/nlp.jl +++ b/src/nlp.jl @@ -14,8 +14,10 @@ struct _Subproblem rows::Vector{MOI.ConstraintIndex} end -function _build_subproblem(::Nothing, model::Optimizer, problem::_Problem) - nlp = _instantiate(model.nlp_solver) +# 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) @@ -35,6 +37,36 @@ function _build_subproblem(::Nothing, model::Optimizer, problem::_Problem) 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) @@ -72,21 +104,11 @@ function _solve_nlp( warm_start; deadline::Float64 = Inf ) - for (binary, value) in combination - MOI.set(sub.model, MOI.ConstraintSet(), sub.fixes[binary], - MOI.EqualTo(value ? 1.0 : 0.0)) - end for ci in sub.rows MOI.delete(sub.model, ci) end - empty!(sub.rows) - for disjunction in problem.disjunctions, disjunct in disjunction.disjuncts - _disjunct_active(combination, disjunct) || continue - for (func, set) in zip(disjunct.functions, disjunct.sets) - push!(sub.rows, MOI.add_constraint(sub.model, - _map_to(sub.variable_map, func), set)) - end - end + 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) @@ -106,6 +128,172 @@ function _solve_nlp( 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) ################################################################################ @@ -118,14 +306,14 @@ _nlpf_slacked(func, u, ::MOI.AbstractScalarSet) = nothing # The slacked feasibility NLP: one nonnegative `u` relaxes every scalar # inequality row (bounds and equalities stay exact) and is minimized. # Its solution is a linearization site for an infeasible combination. -function _solve_nlpf( +# 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, - warm_start; - deadline::Float64 = Inf + combination::AbstractDict ) - nlp = _instantiate(model.nlp_solver) variable_map = Dict{MOI.VariableIndex, MOI.VariableIndex}( vi => MOI.add_variable(nlp) for vi in problem.variables) u = MOI.add_variable(nlp) @@ -153,12 +341,20 @@ function _solve_nlpf( for (func, set) in rows mapped = _map_to(variable_map, func) slacked = _nlpf_slacked(mapped, u, set) - if slacked === nothing - MOI.add_constraint(nlp, mapped, set) - else - MOI.add_constraint(nlp, slacked, set) - end + 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) diff --git a/src/optimizer.jl b/src/optimizer.jl index 01003d2..f4aa7f3 100644 --- a/src/optimizer.jl +++ b/src/optimizer.jl @@ -119,6 +119,21 @@ 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`), 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 + """ SubproblemMethod() <: AbstractAlgorithmAttribute -> Any @@ -174,13 +189,15 @@ mutable struct Optimizer <: MOI.AbstractOptimizer relative_gap::Float64 raw_status::String solve_time::Float64 + num_master_solves::Int + num_nlp_solves::Int 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) + Dict{MOI.VariableIndex, Float64}(), NaN, nothing, NaN, "", NaN, 0, 0) end _algorithm(model::Optimizer) = @@ -209,6 +226,8 @@ function _reset_results(model::Optimizer) model.relative_gap = NaN model.raw_status = "" model.solve_time = NaN + model.num_master_solves = 0 + model.num_nlp_solves = 0 return end @@ -466,6 +485,26 @@ 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 + +MOI.is_set_by_optimize(::Union{MasterSolveCount, NLPSolveCount}) = true +MOI.get(model::Optimizer, ::MasterSolveCount) = model.num_master_solves +MOI.get(model::Optimizer, ::NLPSolveCount) = model.num_nlp_solves + function MOI.get(model::Optimizer, ::MOI.ResultCount) return model.primal_status == MOI.NO_SOLUTION ? 0 : 1 end diff --git a/test/loa.jl b/test/loa.jl index 42e8bab..8b3422c 100644 --- a/test/loa.jl +++ b/test/loa.jl @@ -70,6 +70,111 @@ function test_quadratic_row_shapes() @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 + +# 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()) + 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 + function test_quadratic_objective() model = Model(_loa_optimizer()) set_silent(model) @@ -867,6 +972,9 @@ end test_nested_disjunction() test_nested_disjunction_vacuous() test_quadratic_row_shapes() + test_batched_subproblems_loa() + test_batched_matches_sequential() + test_multi_generation_cuts() test_quadratic_objective() test_max_sense_linear() test_two_disjunctions() From 5dcfd78062885b67acc3bbce4de7a8d2bbb75ea5 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Thu, 10 Sep 2026 19:17:33 -0400 Subject: [PATCH 12/18] Add a CombinationSource attribute for multi-generation `combination_sources.jl` holds the default pool-then-re-solve source alongside the experimental `Neighborhood`, `LPRounding`, `RandomCombinations` and `CutoffResolve`, none of which need anything from the solver. A source extends `_candidate_combinations`. --- Project.toml | 1 + src/DisjunctiveAlgorithms.jl | 2 + src/algorithms/LOA.jl | 54 +++--- src/combination_sources.jl | 323 +++++++++++++++++++++++++++++++++++ src/optimizer.jl | 14 ++ test/loa.jl | 34 ++++ 6 files changed, 401 insertions(+), 27 deletions(-) create mode 100644 src/combination_sources.jl diff --git a/Project.toml b/Project.toml index 7bd9b24..4c65d74 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ 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.2" diff --git a/src/DisjunctiveAlgorithms.jl b/src/DisjunctiveAlgorithms.jl index b25fdc3..c04d840 100644 --- a/src/DisjunctiveAlgorithms.jl +++ b/src/DisjunctiveAlgorithms.jl @@ -1,12 +1,14 @@ 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("master.jl") +include("combination_sources.jl") include("nlp.jl") include("cuts.jl") include("algorithms/LOA.jl") diff --git a/src/algorithms/LOA.jl b/src/algorithms/LOA.jl index 2b3b582..075137c 100644 --- a/src/algorithms/LOA.jl +++ b/src/algorithms/LOA.jl @@ -25,6 +25,7 @@ never `OPTIMAL`. - [`SlackTolerance`](@ref) - [`IterationTimeLimit`](@ref) - [`MultiGenerationSize`](@ref) +- [`CombinationSource`](@ref) - [`SubproblemMethod`](@ref) """ mutable struct LOA <: AbstractAlgorithm @@ -39,10 +40,11 @@ mutable struct LOA <: AbstractAlgorithm 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, nothing, nothing, nothing, nothing, nothing) end _default(::Algorithm) = LOA() @@ -97,7 +99,7 @@ for (attr, field) in ( end end -# `nothing` is this attribute's default, so it skips the +# `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) @@ -106,6 +108,13 @@ function MOI.set(algorithm::LOA, ::SubproblemMethod, value) 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) = @@ -186,35 +195,24 @@ function _solve_master(model::Optimizer, master::_Master, deadline::Float64) return _solved_and_feasible(master.model) end -# The master's solution pool first (result indices past 1), then -# re-solves behind a no-good cut for whatever the pool did not supply. -# Those cuts are the ones `process_result` would add later, so the -# second return value tells the caller to skip them. +# 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 ) - combinations = [_extract_combination(problem, master)] - for index in 2:min(count, MOI.get(master.model, MOI.ResultCount())) - combination = _extract_combination(problem, master, index) - combination in combinations || push!(combinations, combination) - end - excluded = false - while length(combinations) < count && time() < deadline - if !excluded - foreach(c -> _avoid_combination(master, c), combinations) - excluded = true - end - _solve_master(model, master, deadline) || break - combination = _extract_combination(problem, master) - combination in combinations && break - push!(combinations, combination) - _avoid_combination(master, combination) - end - return combinations, excluded + 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::_Master, sense, objective) @@ -335,9 +333,11 @@ function _optimize!(algorithm::LOA, model::Optimizer) break end end - combinations, excluded = _extract_combinations(model, problem, - master, MOI.get(algorithm, MultiGenerationSize()), - loop_deadline) + combinations, excluded = _extract_combinations( + MOI.get(algorithm, CombinationSource()), model, problem, + master, best_result === nothing ? nothing : + best_result.combination, + MOI.get(algorithm, MultiGenerationSize()), loop_deadline) results = _solve_nlps(method, model, problem, subproblem, combinations, warm_start(); deadline = loop_deadline) model.num_nlp_solves += length(results) diff --git a/src/combination_sources.jl b/src/combination_sources.jl new file mode 100644 index 0000000..4e4aae2 --- /dev/null +++ b/src/combination_sources.jl @@ -0,0 +1,323 @@ +################################################################################ +# 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, THEN RE-SOLVES +################################################################################ +# Result indices past 1 are the solver's pool (Gurobi's PoolSolutions); +# whatever it does not supply comes from 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 + length(candidates) >= count && 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) + for (F, S) in MOI.get(lp, MOI.ListOfConstraintTypesPresent()) + 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 index f4aa7f3..ba6b08c 100644 --- a/src/optimizer.jl +++ b/src/optimizer.jl @@ -134,6 +134,20 @@ 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 re-solves the master behind no-good +cuts for the rest. 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 diff --git a/test/loa.jl b/test/loa.jl index 8b3422c..007b4bb 100644 --- a/test/loa.jl +++ b/test/loa.jl @@ -175,6 +175,39 @@ function test_multi_generation_cuts() 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) @@ -975,6 +1008,7 @@ end 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() From c56e09e371fd9b26823630466d3fd48b1ef1f012 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Fri, 11 Sep 2026 18:42:08 -0400 Subject: [PATCH 13/18] Loosen the DisjunctiveProgramming compat bound to 0.6 The README install note pinned a patch version, when what the branch requirement is really about is `DisjunctionSet` and `Direct()` not being in a registered release yet. Say that instead. --- Project.toml | 2 +- README.md | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 4c65d74..8b4b711 100644 --- a/Project.toml +++ b/Project.toml @@ -9,7 +9,7 @@ MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [compat] -DisjunctiveProgramming = "0.6.2" +DisjunctiveProgramming = "0.6" MathOptInterface = "1" julia = "1.10" diff --git a/README.md b/README.md index 5b2300c..143013b 100644 --- a/README.md +++ b/README.md @@ -20,10 +20,10 @@ options selected through optimizer attributes. ## Installation -DisjunctiveAlgorithms requires DisjunctiveProgramming v0.6.2, which is -not registered yet: it lives on the `gdp_optimizer` branch (it adds -`DisjunctionSet` and the `Direct()` reformulation). Until that release -ships, install the branch first: +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 From 2f408ddb843022f65c3afc45f7fdadd01f51f246 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Fri, 11 Sep 2026 19:41:14 -0400 Subject: [PATCH 14/18] Add master and NLP solve-time attributes `MasterSolveTime` and `NLPSolveTime` report where a run spent its time. The default combination source now takes as many pool entries as the solver holds, and re-solves the master only when there is no pool at all. --- src/algorithms/LOA.jl | 6 +++++- src/combination_sources.jl | 8 +++---- src/optimizer.jl | 44 ++++++++++++++++++++++++++++++-------- test/loa.jl | 4 ++++ 4 files changed, 48 insertions(+), 14 deletions(-) diff --git a/src/algorithms/LOA.jl b/src/algorithms/LOA.jl index 075137c..a4cc586 100644 --- a/src/algorithms/LOA.jl +++ b/src/algorithms/LOA.jl @@ -190,7 +190,7 @@ end function _solve_master(model::Optimizer, master::_Master, deadline::Float64) _cap_remaining_time(master.model, deadline) - MOI.optimize!(master.model) + model.master_time += @elapsed MOI.optimize!(master.model) model.num_master_solves += 1 return _solved_and_feasible(master.model) end @@ -293,8 +293,10 @@ function _optimize!(algorithm::LOA, model::Optimizer) master_status = status break end + 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 @@ -338,8 +340,10 @@ function _optimize!(algorithm::LOA, model::Optimizer) 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) diff --git a/src/combination_sources.jl b/src/combination_sources.jl index 4e4aae2..d00d0ed 100644 --- a/src/combination_sources.jl +++ b/src/combination_sources.jl @@ -36,11 +36,11 @@ _push_candidate!(candidates, combination, proposal) = push!(candidates, combination) ################################################################################ -# DEFAULT: POOL, THEN RE-SOLVES +# DEFAULT: POOL, OR RE-SOLVES ################################################################################ # Result indices past 1 are the solver's pool (Gurobi's PoolSolutions); -# whatever it does not supply comes from re-solving behind a no-good -# cut per combination taken. +# 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, @@ -56,7 +56,7 @@ function _candidate_combinations( _push_candidate!(candidates, _extract_combination(problem, master, index), proposal) end - length(candidates) >= count && return candidates, false + MOI.get(master.model, MOI.ResultCount()) > 1 && return candidates, false return _resolve_candidates(model, problem, master, proposal, candidates, count, deadline) end diff --git a/src/optimizer.jl b/src/optimizer.jl index ba6b08c..6d7dd88 100644 --- a/src/optimizer.jl +++ b/src/optimizer.jl @@ -126,8 +126,9 @@ 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`), and otherwise from re-solving the master behind a -no-good cut per combination already taken. All of them go through +`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. """ @@ -139,11 +140,12 @@ _default(::MultiGenerationSize) = 1 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 re-solves the master behind no-good -cuts for the rest. 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`. +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 @@ -205,13 +207,16 @@ mutable struct Optimizer <: MOI.AbstractOptimizer solve_time::Float64 num_master_solves::Int num_nlp_solves::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) + Dict{MOI.VariableIndex, Float64}(), NaN, nothing, NaN, "", NaN, 0, 0, + 0.0, 0.0) end _algorithm(model::Optimizer) = @@ -242,6 +247,8 @@ function _reset_results(model::Optimizer) model.solve_time = NaN model.num_master_solves = 0 model.num_nlp_solves = 0 + model.master_time = 0.0 + model.nlp_time = 0.0 return end @@ -515,9 +522,28 @@ Indicator combinations evaluated by NLP subproblems in the last """ struct NLPSolveCount <: MOI.AbstractModelAttribute end -MOI.is_set_by_optimize(::Union{MasterSolveCount, NLPSolveCount}) = true +""" + 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, + 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, ::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 diff --git a/test/loa.jl b/test/loa.jl index 007b4bb..d472a3c 100644 --- a/test/loa.jl +++ b/test/loa.jl @@ -166,6 +166,10 @@ function test_multi_generation_cuts() @test MOI.get(optimizer, DA.MasterSolveCount()) >= 1 @test MOI.get(model, DA.NLPSolveCount()) == MOI.get(optimizer, 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) From 93c65de15262b895be6a696546e5c2f417caad2c Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Mon, 14 Sep 2026 14:08:40 -0400 Subject: [PATCH 15/18] Seed LOA from a logic-only set cover The covering pass was solving the full OA master, cuts and continuous variables included, so initialization was already MILP/NLP instead of the NLP-only pass Turkay and Grossmann describe. `_SetCoverModel` is the disjunct binaries under the exactly-one and propositional rows only, with its own no-good cuts, and stops once a solve covers nothing new. The master MILP now first solves in the main loop. Also adds `NLPInfeasibleCount`, and drops indicator rows from the LP copy in `_relaxed_activations`, where a gated row has no relaxation without its binary. --- src/algorithms/LOA.jl | 48 +++++++++++++++++++++++--------------- src/combination_sources.jl | 7 ++++++ src/master.jl | 43 ++++++++++++++++++++++++++++++++-- src/optimizer.jl | 20 +++++++++++++--- test/loa.jl | 45 +++++++++++++++++++++++++++++++---- test/mock_optimizer.jl | 13 ++++++++--- 6 files changed, 145 insertions(+), 31 deletions(-) diff --git a/src/algorithms/LOA.jl b/src/algorithms/LOA.jl index a4cc586..5eb710c 100644 --- a/src/algorithms/LOA.jl +++ b/src/algorithms/LOA.jl @@ -7,8 +7,9 @@ `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, seeded by a -set-covering pass over the nonlinear disjuncts. The OA bound is only +NLP subproblems solved at fixed indicator combinations, seeded by +the set-covering pass: a logic-only MIP over the indicators that +activates every nonlinear disjunct once. The OA bound is only valid on convex problems, so convergence reports `LOCALLY_SOLVED`, never `OPTIMAL`. @@ -162,7 +163,7 @@ end # GDPopt's covering weights: an uncovered disjunct outweighs all # covered ones function _cover_objective( - master::_Master, + set_cover::_SetCoverModel, cover::Vector{_Disjunct}, needs_cover, num_covered::Int @@ -170,7 +171,7 @@ function _cover_objective( objective = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0) for i in eachindex(cover) weight = Float64(needs_cover[i] ? num_covered + 1 : 1) - activation = _map_to(master.variable_map, cover[i].activation) + activation = _map_to(set_cover.variable_map, cover[i].activation) objective = MOI.Utilities.operate(+, Float64, objective, MOI.Utilities.operate(*, Float64, weight, activation)) end @@ -188,7 +189,11 @@ function _user_start_values(model::Optimizer, problem::_Problem) return isempty(point) ? nothing : (point = point,) end -function _solve_master(model::Optimizer, master::_Master, deadline::Float64) +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 @@ -215,7 +220,11 @@ function _extract_combinations( return [proposal; extra], excluded end -function _set_master_objective(master::_Master, sense, objective) +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}}(), @@ -260,6 +269,7 @@ function _optimize!(algorithm::LOA, model::Optimizer) 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 && @@ -274,25 +284,25 @@ function _optimize!(algorithm::LOA, model::Optimizer) warm_start = () -> previous_result === nothing ? nothing : previous_result.point - # set covering: reuse the master with a coverage objective so - # every nonlinear disjunct gets visited once + # set covering: the logic-only problem picks combinations that + # activate every nonlinear disjunct once; the master only receives + # their NLPs' cuts cover = _cover_disjuncts(problem) needs_cover = trues(length(cover)) num_covered = 0 + set_cover = _build_set_cover(model, problem) for iteration in 1:MOI.get(algorithm, SetCoverIterationLimit()) (iteration == 1 || any(needs_cover)) || break time() < loop_deadline || break - _set_master_objective(master, MOI.MAX_SENSE, - _cover_objective(master, cover, needs_cover, num_covered)) - solved = _solve_master(model, master, loop_deadline) - # capture the status before the objective restore invalidates it - status = MOI.get(master.model, MOI.TerminationStatus()) - combination = solved ? _extract_combination(problem, master) : nothing - _set_master_objective(master, master.sense, master.oa_objective) - if !solved - master_status = status - break - end + _set_master_objective(set_cover, MOI.MAX_SENSE, + _cover_objective(set_cover, cover, needs_cover, num_covered)) + _solve_master(model, set_cover, loop_deadline) || break + combination = _extract_combination(problem, set_cover) + # 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 + _avoid_combination(set_cover, combination) t_nlp = time() result = _solve_nlp(method, model, problem, subproblem, combination, warm_start(); deadline = loop_deadline) diff --git a/src/combination_sources.jl b/src/combination_sources.jl index d00d0ed..fb12ef2 100644 --- a/src/combination_sources.jl +++ b/src/combination_sources.jl @@ -242,7 +242,14 @@ function _relaxed_activations( ) 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) diff --git a/src/master.jl b/src/master.jl index f6f13bf..ba2beb4 100644 --- a/src/master.jl +++ b/src/master.jl @@ -12,6 +12,14 @@ mutable struct _Master 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 + function _instantiate(factory) solver = MOI.instantiate(factory; with_cache_type = Float64, with_bridge_type = Float64) @@ -73,6 +81,30 @@ function _build_master(model::Optimizer, problem::_Problem) 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( @@ -153,7 +185,11 @@ function _add_penalized_slack( end # round to Bool; MILP values are only within integer tolerance -function _extract_combination(problem::_Problem, master::_Master, index = 1) +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])) @@ -161,7 +197,10 @@ function _extract_combination(problem::_Problem, master::_Master, index = 1) end # No-good cut: active `1 - z` plus inactive `z` terms must reach 1 -function _avoid_combination(master::_Master, combination::AbstractDict) +function _avoid_combination( + master::Union{_Master, _SetCoverModel}, + combination::AbstractDict + ) terms = MOI.ScalarAffineTerm{Float64}[] constant = 0.0 for (binary, value) in combination diff --git a/src/optimizer.jl b/src/optimizer.jl index 6d7dd88..5304e5a 100644 --- a/src/optimizer.jl +++ b/src/optimizer.jl @@ -45,7 +45,9 @@ _default(::NumIterationLimit) = 10 """ SetCoverIterationLimit() <: AbstractAlgorithmAttribute -> Int -Set-covering initialization iterations. Defaults to `8`. +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 @@ -207,6 +209,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer solve_time::Float64 num_master_solves::Int num_nlp_solves::Int + num_nlp_infeasible::Int master_time::Float64 nlp_time::Float64 end @@ -216,7 +219,7 @@ function Optimizer(nlp_solver, mip_solver = nlp_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, 0.0, 0.0) end _algorithm(model::Optimizer) = @@ -247,6 +250,7 @@ function _reset_results(model::Optimizer) 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 @@ -522,6 +526,15 @@ Indicator combinations evaluated by NLP subproblems in the last """ 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 @@ -538,10 +551,11 @@ feasibility restoration and stacked batches included. struct NLPSolveTime <: MOI.AbstractModelAttribute end const _SolveStatistic = Union{MasterSolveCount, NLPSolveCount, - MasterSolveTime, NLPSolveTime} + 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 diff --git a/test/loa.jl b/test/loa.jl index d472a3c..af065ad 100644 --- a/test/loa.jl +++ b/test/loa.jl @@ -88,6 +88,37 @@ function _batched_test_model(attrs::Pair...) 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 @@ -166,6 +197,8 @@ function test_multi_generation_cuts() @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 @@ -811,11 +844,14 @@ function test_time_limit_with_incumbent() @test occursin("limit hit", raw_status(model)) end -# The master finishes abnormally after an incumbent exists: the mock -# master reports a node limit on its second solve. +# 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() - factory = () -> DA.Optimizer(Ipopt.Optimizer, - () -> MockSolver(HiGHS.Optimizer, fail_from = 2)) + 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) @@ -1052,6 +1088,7 @@ end end @testset "LOA units" begin + test_set_cover_model() test_cut_term_directions() test_activation_binary() test_sense_primitives() diff --git a/test/mock_optimizer.jl b/test/mock_optimizer.jl index b21c826..1716a47 100644 --- a/test/mock_optimizer.jl +++ b/test/mock_optimizer.jl @@ -4,6 +4,10 @@ # 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. @@ -14,6 +18,7 @@ mutable struct MockSolver <: MOI.AbstractOptimizer fail_status::MOI.TerminationStatusCode strict_constants::Bool solves::Int + shared_solves::Ref{Int} failing::Bool end @@ -22,16 +27,18 @@ function MockSolver( sleep_time::Float64 = 0.0, fail_from::Int = typemax(Int), fail_status::MOI.TerminationStatusCode = MOI.NODE_LIMIT, - strict_constants::Bool = false + strict_constants::Bool = false, + shared_solves::Ref{Int} = Ref(0) ) return MockSolver(MOI.instantiate(factory), sleep_time, fail_from, - fail_status, strict_constants, 0, false) + 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.solves >= model.fail_from + model.failing = model.shared_solves[] >= model.fail_from model.failing || MOI.optimize!(model.inner) return end From 2f1a41582e31d75c768e06df241e5639d8d6e1db Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Tue, 15 Sep 2026 12:37:01 -0400 Subject: [PATCH 16/18] Move the LOA files under src/algorithms/LOA master.jl, nlp.jl and cuts.jl are specific to logic-based outer approximation rather than general, so they belong beside LOA.jl; `_instantiate` moves to optimizer.jl. --- src/DisjunctiveAlgorithms.jl | 6 +++--- src/{ => algorithms/LOA}/cuts.jl | 0 src/{ => algorithms/LOA}/master.jl | 8 -------- src/{ => algorithms/LOA}/nlp.jl | 0 src/optimizer.jl | 9 +++++++++ 5 files changed, 12 insertions(+), 11 deletions(-) rename src/{ => algorithms/LOA}/cuts.jl (100%) rename src/{ => algorithms/LOA}/master.jl (97%) rename src/{ => algorithms/LOA}/nlp.jl (100%) diff --git a/src/DisjunctiveAlgorithms.jl b/src/DisjunctiveAlgorithms.jl index c04d840..60155c3 100644 --- a/src/DisjunctiveAlgorithms.jl +++ b/src/DisjunctiveAlgorithms.jl @@ -7,10 +7,10 @@ import DisjunctiveProgramming: DisjunctionSet, activation_index, include("optimizer.jl") include("problem.jl") -include("master.jl") +include("algorithms/LOA/master.jl") include("combination_sources.jl") -include("nlp.jl") -include("cuts.jl") +include("algorithms/LOA/nlp.jl") +include("algorithms/LOA/cuts.jl") include("algorithms/LOA.jl") end diff --git a/src/cuts.jl b/src/algorithms/LOA/cuts.jl similarity index 100% rename from src/cuts.jl rename to src/algorithms/LOA/cuts.jl diff --git a/src/master.jl b/src/algorithms/LOA/master.jl similarity index 97% rename from src/master.jl rename to src/algorithms/LOA/master.jl index ba2beb4..42a735e 100644 --- a/src/master.jl +++ b/src/algorithms/LOA/master.jl @@ -20,14 +20,6 @@ struct _SetCoverModel variable_map::Dict{MOI.VariableIndex, MOI.VariableIndex} end -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 - _map_to(variable_map::AbstractDict, func) = MOI.Utilities.map_indices(vi -> variable_map[vi], func) diff --git a/src/nlp.jl b/src/algorithms/LOA/nlp.jl similarity index 100% rename from src/nlp.jl rename to src/algorithms/LOA/nlp.jl diff --git a/src/optimizer.jl b/src/optimizer.jl index 5304e5a..2d4d184 100644 --- a/src/optimizer.jl +++ b/src/optimizer.jl @@ -225,6 +225,15 @@ 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__)) From e70fdfa2712858f4e6a92fdef284d21c20c1b440 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Thu, 24 Sep 2026 13:16:15 -0400 Subject: [PATCH 17/18] Initialization options --- src/algorithms/LOA.jl | 38 +++++++++++++++++++++++++++--------- src/algorithms/LOA/master.jl | 9 +++++++-- src/optimizer.jl | 13 ++++++++++++ test/loa.jl | 28 ++++++++++++++++++++++++++ 4 files changed, 77 insertions(+), 11 deletions(-) diff --git a/src/algorithms/LOA.jl b/src/algorithms/LOA.jl index 5eb710c..dd1911c 100644 --- a/src/algorithms/LOA.jl +++ b/src/algorithms/LOA.jl @@ -17,6 +17,7 @@ never `OPTIMAL`. - [`NumIterationLimit`](@ref) - [`SetCoverIterationLimit`](@ref) +- [`SetCoverProblem`](@ref) - [`M_Value`](@ref) - [`MasterReformulation`](@ref) - [`MaxSlack`](@ref) @@ -32,6 +33,7 @@ never `OPTIMAL`. mutable struct LOA <: AbstractAlgorithm num_iteration_limit::Union{Nothing, Int} set_cover_iteration_limit::Union{Nothing, Int} + set_cover_problem::Union{Nothing, String} m_value::Union{Nothing, Float64} master_reformulation::Union{Nothing, String} max_slack::Union{Nothing, Float64} @@ -45,7 +47,8 @@ mutable struct LOA <: AbstractAlgorithm subproblem_method::Any LOA() = new(nothing, nothing, nothing, nothing, nothing, nothing, - nothing, nothing, nothing, nothing, nothing, nothing, nothing) + nothing, nothing, nothing, nothing, nothing, nothing, nothing, + nothing) end _default(::Algorithm) = LOA() @@ -61,6 +64,12 @@ function _validate(::MasterReformulation, value) return nothing end +function _validate(::SetCoverProblem, value) + value in ("logic", "master") || error( + "`SetCoverProblem` must be \"logic\" or \"master\" " * + "(got `$value`).") + return nothing +end function _validate(::OASlack, value) value > 0 || error("`OASlack` must be positive (got `$value`).") @@ -77,6 +86,7 @@ end 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), @@ -163,7 +173,7 @@ end # GDPopt's covering weights: an uncovered disjunct outweighs all # covered ones function _cover_objective( - set_cover::_SetCoverModel, + set_cover::Union{_Master, _SetCoverModel}, cover::Vector{_Disjunct}, needs_cover, num_covered::Int @@ -284,25 +294,35 @@ function _optimize!(algorithm::LOA, model::Optimizer) warm_start = () -> previous_result === nothing ? nothing : previous_result.point - # set covering: the logic-only problem picks combinations that - # activate every nonlinear disjunct once; the master only receives - # their NLPs' cuts + # set covering: the covering problem picks combinations that + # activate every nonlinear disjunct once. Under `SetCoverProblem() + # = "logic"` it is the logic-only problem and the master only + # receives the NLPs' cuts; under `"master"` the master itself is + # solved with the covering objective swapped in and out, so the + # seeds satisfy the linear rows (GDPopt) and the no-good cuts come + # from `process_result`. cover = _cover_disjuncts(problem) needs_cover = trues(length(cover)) num_covered = 0 - set_cover = _build_set_cover(model, problem) + cover_master = MOI.get(algorithm, SetCoverProblem()) == "master" + set_cover = cover_master ? master : _build_set_cover(model, problem) 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)) - _solve_master(model, set_cover, loop_deadline) || break - combination = _extract_combination(problem, set_cover) + solved = _solve_master(model, set_cover, loop_deadline) + # the master loses its solution when its objective is restored + combination = solved ? _extract_combination(problem, set_cover) : + nothing + cover_master && _set_master_objective(master, sense, + master.oa_objective) + 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 - _avoid_combination(set_cover, combination) + cover_master || _avoid_combination(set_cover, combination) t_nlp = time() result = _solve_nlp(method, model, problem, subproblem, combination, warm_start(); deadline = loop_deadline) diff --git a/src/algorithms/LOA/master.jl b/src/algorithms/LOA/master.jl index 42a735e..a67462c 100644 --- a/src/algorithms/LOA/master.jl +++ b/src/algorithms/LOA/master.jl @@ -29,10 +29,15 @@ function _solved_and_feasible(solver::MOI.ModelLike) 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()) && - MOI.set(solver, MOI.TimeLimitSec(), max(0.0, deadline - time())) + 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 diff --git a/src/optimizer.jl b/src/optimizer.jl index 2d4d184..b5dcc53 100644 --- a/src/optimizer.jl +++ b/src/optimizer.jl @@ -52,6 +52,19 @@ disjunct has been active in a feasible NLP. Defaults to `8`. struct SetCoverIterationLimit <: AbstractAlgorithmAttribute end _default(::SetCoverIterationLimit) = 8 +""" + SetCoverProblem() <: AbstractAlgorithmAttribute -> String + +Which problem the set-covering pass solves for its seed combinations: +`"logic"` (the disjunct binaries under the exactly-one and +propositional rows only, Turkay and Grossmann's covering problem) or +`"master"` (the full MILP master with the covering objective swapped +in, as Pyomo GDPopt does, so the seeds also satisfy the linear rows +and their NLPs are far less often infeasible). Defaults to `"logic"`. +""" +struct SetCoverProblem <: AbstractAlgorithmAttribute end +_default(::SetCoverProblem) = "logic" + """ M_Value() <: AbstractAlgorithmAttribute -> Float64 diff --git a/test/loa.jl b/test/loa.jl index af065ad..df07eb7 100644 --- a/test/loa.jl +++ b/test/loa.jl @@ -668,6 +668,33 @@ function test_vector_of_variables_disjunction() @test value(zin[2]) ≈ 1.0 atol = 1e-5 end +# `SetCoverProblem() => "master"` seeds 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 seed draws a combination whose +# NLP is infeasible, while the master seed lands on a feasible one and +# the run ends with the same incumbent. +function test_set_cover_problem_master() + for cover in ("logic", "master") + 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 == "master" && @test infeasible == 0 + end + @test_throws ErrorException Model(_loa_optimizer( + DA.SetCoverProblem() => "cuts")) +end + # Zero iteration budgets exit before any solve, without an incumbent. function test_iteration_limit_no_incumbent() model = Model(_loa_optimizer(DA.SetCoverIterationLimit() => 0, @@ -1089,6 +1116,7 @@ end @testset "LOA units" begin test_set_cover_model() + test_set_cover_problem_master() test_cut_term_directions() test_activation_binary() test_sense_primitives() From 2d2c1328c42fcf34fc3384b2082cb3d58d436047 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Thu, 24 Sep 2026 14:28:26 -0400 Subject: [PATCH 18/18] . --- Project.toml | 2 +- src/algorithms/LOA.jl | 57 ++++++++++++++++++++++++++----------------- src/optimizer.jl | 39 +++++++++++++++++++++++------ test/loa.jl | 14 +++++------ 4 files changed, 74 insertions(+), 38 deletions(-) diff --git a/Project.toml b/Project.toml index 8b4b711..e922c08 100644 --- a/Project.toml +++ b/Project.toml @@ -9,7 +9,7 @@ MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [compat] -DisjunctiveProgramming = "0.6" +DisjunctiveProgramming = "0.6, 0.7" MathOptInterface = "1" julia = "1.10" diff --git a/src/algorithms/LOA.jl b/src/algorithms/LOA.jl index dd1911c..93b4867 100644 --- a/src/algorithms/LOA.jl +++ b/src/algorithms/LOA.jl @@ -7,9 +7,12 @@ `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, seeded by -the set-covering pass: a logic-only MIP over the indicators that -activates every nonlinear disjunct once. The OA bound is only +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`. @@ -33,7 +36,7 @@ never `OPTIMAL`. mutable struct LOA <: AbstractAlgorithm num_iteration_limit::Union{Nothing, Int} set_cover_iteration_limit::Union{Nothing, Int} - set_cover_problem::Union{Nothing, String} + set_cover_problem::Union{Nothing, AbstractSetCover} m_value::Union{Nothing, Float64} master_reformulation::Union{Nothing, String} max_slack::Union{Nothing, Float64} @@ -65,9 +68,8 @@ function _validate(::MasterReformulation, value) end function _validate(::SetCoverProblem, value) - value in ("logic", "master") || error( - "`SetCoverProblem` must be \"logic\" or \"master\" " * - "(got `$value`).") + value isa AbstractSetCover || error("`SetCoverProblem` must be an " * + "`AbstractSetCover`, e.g. `LogicOnlyCover()` (got `$value`).") return nothing end function _validate(::OASlack, value) @@ -170,8 +172,25 @@ function _cover_disjuncts(problem::_Problem) return cover end -# GDPopt's covering weights: an uncovered disjunct outweighs all -# covered ones +# 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}, @@ -294,35 +313,29 @@ function _optimize!(algorithm::LOA, model::Optimizer) warm_start = () -> previous_result === nothing ? nothing : previous_result.point - # set covering: the covering problem picks combinations that - # activate every nonlinear disjunct once. Under `SetCoverProblem() - # = "logic"` it is the logic-only problem and the master only - # receives the NLPs' cuts; under `"master"` the master itself is - # solved with the covering objective swapped in and out, so the - # seeds satisfy the linear rows (GDPopt) and the no-good cuts come - # from `process_result`. + # 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 - cover_master = MOI.get(algorithm, SetCoverProblem()) == "master" - set_cover = cover_master ? master : _build_set_cover(model, problem) + 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) - # the master loses its solution when its objective is restored combination = solved ? _extract_combination(problem, set_cover) : nothing - cover_master && _set_master_objective(master, sense, - master.oa_objective) + _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 - cover_master || _avoid_combination(set_cover, combination) + _exclude_combination(covering_algorithm, set_cover, combination) t_nlp = time() result = _solve_nlp(method, model, problem, subproblem, combination, warm_start(); deadline = loop_deadline) diff --git a/src/optimizer.jl b/src/optimizer.jl index b5dcc53..e460a2c 100644 --- a/src/optimizer.jl +++ b/src/optimizer.jl @@ -53,17 +53,40 @@ struct SetCoverIterationLimit <: AbstractAlgorithmAttribute end _default(::SetCoverIterationLimit) = 8 """ - SetCoverProblem() <: AbstractAlgorithmAttribute -> String + AbstractSetCover -Which problem the set-covering pass solves for its seed combinations: -`"logic"` (the disjunct binaries under the exactly-one and -propositional rows only, Turkay and Grossmann's covering problem) or -`"master"` (the full MILP master with the covering objective swapped -in, as Pyomo GDPopt does, so the seeds also satisfy the linear rows -and their NLPs are far less often infeasible). Defaults to `"logic"`. +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) = "logic" +_default(::SetCoverProblem) = LogicOnlyCover() """ M_Value() <: AbstractAlgorithmAttribute -> Float64 diff --git a/test/loa.jl b/test/loa.jl index df07eb7..a4f4264 100644 --- a/test/loa.jl +++ b/test/loa.jl @@ -668,13 +668,13 @@ function test_vector_of_variables_disjunction() @test value(zin[2]) ≈ 1.0 atol = 1e-5 end -# `SetCoverProblem() => "master"` seeds 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 seed draws a combination whose -# NLP is infeasible, while the master seed lands on a feasible one and +# `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 ("logic", "master") + for cover in (DA.LogicOnlyCover(), DA.MasterCover()) model = Model(_loa_optimizer(DA.SetCoverProblem() => cover)) set_silent(model) @variable(model, 0 <= x <= 10) @@ -689,10 +689,10 @@ function test_set_cover_problem_master() @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 == "master" && @test infeasible == 0 + cover isa DA.MasterCover && @test infeasible == 0 end @test_throws ErrorException Model(_loa_optimizer( - DA.SetCoverProblem() => "cuts")) + DA.SetCoverProblem() => "master")) end # Zero iteration budgets exit before any solve, without an incumbent.