diff --git a/Project.toml b/Project.toml index bdec9be4..618e92a0 100644 --- a/Project.toml +++ b/Project.toml @@ -15,7 +15,10 @@ InfiniteDisjunctiveProgramming = "InfiniteOpt" [compat] Aqua = "0.8" +InfiniteOpt = "0.6.3" +Ipopt = "1.9.0" JuMP = "1.18" +Juniper = "0.9.3" Reexport = "1" julia = "1.10" Juniper = "0.9.3" @@ -25,9 +28,10 @@ InfiniteOpt = "0.6.3" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" Juniper = "2ddba703-00a4-53a7-87a5-e8b9971dde84" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["Aqua", "HiGHS", "Test", "Juniper", "Ipopt", "InfiniteOpt"] diff --git a/docs/src/index.md b/docs/src/index.md index aa1902e3..615d598e 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -61,7 +61,7 @@ A [`GDPModel`](@ref) is a `JuMP Model` with a [`GDPData`](@ref) field in the mod - `Logical Constraints`: Selector (cardinality) or proposition (Boolean) constraints describing the relationships between the logical variables. - `Disjunct Constraints`: Constraints associated with each disjunct in the model. - `Disjunctions`: Disjunction constraints. -- `Solution Method`: The reformulation technique or solution method. Currently, supported methods include Big-M, Hull, and Indicator Constraints. +- `Solution Method`: The reformulation technique or solution method. Currently, supported methods include Big-M, Multiple Big-M, Hull, P-Split, Indicator Constraints, and Logic-based Outer Approximation (LOA). - `Reformulation Variables`: List of JuMP variables created when reformulating a GDP model into a MIP model. - `Reformulation Constraints`: List of constraints created when reformulating a GDP model into a MIP model. - `Ready to Optimize`: Flag indicating if the model can be optimized. @@ -163,9 +163,15 @@ The following reformulation methods are currently supported: 1. [Big-M](https://optimization.cbe.cornell.edu/index.php?title=Disjunctive_inequalities#Big-M_Reformulation[1][2]): The [`BigM`](@ref) struct is used. -2. [Hull](https://optimization.cbe.cornell.edu/index.php?title=Disjunctive_inequalities#Convex-Hull_Reformulation[1][2]): The [`Hull`](@ref) struct is used. +2. Multiple Big-M: A tighter big-M reformulation that computes a separate big-M value for each disjunct constraint (rather than a single global value) by solving auxiliary mini-models. This is invoked with the [`MBM`](@ref) struct, which requires an optimizer. -3. [Indicator](https://jump.dev/JuMP.jl/stable/manual/constraints/#Indicator-constraints): This method reformulates each disjunct constraint into an indicator constraint with the Boolean reformulation counterpart of the Logical variable used to define the disjunct constraint. This is invoked with [`Indicator`](@ref). +3. [Hull](https://optimization.cbe.cornell.edu/index.php?title=Disjunctive_inequalities#Convex-Hull_Reformulation[1][2]): The [`Hull`](@ref) struct is used. + +4. P-Split: A partitioned reformulation that splits the variables into groups and applies a P-split formulation to each group, giving a relaxation between Big-M and Hull in tightness. This is invoked with the [`PSplit`](@ref) struct. + +5. [Indicator](https://jump.dev/JuMP.jl/stable/manual/constraints/#Indicator-constraints): This method reformulates each disjunct constraint into an indicator constraint with the Boolean reformulation counterpart of the Logical variable used to define the disjunct constraint. This is invoked with [`Indicator`](@ref). + +6. Logic-based Outer Approximation (LOA): An iterative solution method for nonlinear GDPs, rather than a single-shot reformulation. It alternates between a primal NLP (the model reformulated by an inner method — `BigM`, `MBM`, or `Hull` — with the disjunct binaries fixed at the current selection) and a master MILP that accumulates outer-approximation and no-good cuts until the bound meets the incumbent. This is invoked with the [`LOA`](@ref) struct, e.g. `optimize!(m, gdp_method = LOA(Ipopt.Optimizer))`. ## Release Notes diff --git a/ext/InfiniteDisjunctiveProgramming.jl b/ext/InfiniteDisjunctiveProgramming.jl index 5f77dce7..318e3d39 100644 --- a/ext/InfiniteDisjunctiveProgramming.jl +++ b/ext/InfiniteDisjunctiveProgramming.jl @@ -297,6 +297,767 @@ end # CUTTING PLANES FOR INFINITEMODEL ################################################################################ +# Build CP subproblem: reformulate the InfiniteModel in-place, transcribe, +# copy, and wrap in GDPSubmodel with forward variable map. +function DP.copy_and_reformulate( + model::InfiniteOpt.InfiniteModel, + decision_vars::Vector{InfiniteOpt.GeneralVariableRef}, + reform_method::DP.AbstractReformulationMethod, + method::DP.CuttingPlanes + ) + DP.reformulate_model(model, reform_method) + InfiniteOpt.build_transformation_backend!(model) + transcribed = InfiniteOpt.transformation_model(model) + transcription_fwd = Dict{InfiniteOpt.GeneralVariableRef, + Vector{JuMP.VariableRef}}() + for v in DP.collect_all_vars(model) + transcription_var = InfiniteOpt.transformation_variable(v) + var_prefs = InfiniteOpt.parameter_refs(v) + transcription_fwd[v] = isempty(var_prefs) ? + [transcription_var] : vec(transcription_var) + end + sub_copy, copy_map = JuMP.copy_model(transcribed) + fwd_map = Dict{InfiniteOpt.GeneralVariableRef, Vector{JuMP.VariableRef}}() + for v in decision_vars + haskey(transcription_fwd, v) || continue + fwd_map[v] = [copy_map[transcribed_var] for transcribed_var in transcription_fwd[v]] + end + sub = DP.GDPSubmodel(sub_copy, decision_vars, fwd_map) + JuMP.set_optimizer(sub.model, method.optimizer) + # NOTE: previously called JuMP.set_silent(sub.model) here, but + # that hides the Gurobi B&B trace even when the caller passes + # OutputFlag=1 / LogToConsole=1 via optimizer_with_attributes. + # The caller is now responsible for silencing if desired. + return sub +end + +# InfiniteOpt-internal primitive: for an InfiniteOpt variable, return +# its per-support values as a Vector. `JuMP.value` returns a scalar +# (finite vars), Vector (1 infinite param), or N-D Array (multiple +# independent param groups). `vcat` lifts a scalar to a 1-element +# Vector; `vec` flattens any N-D Array to a Vector. +_per_support_values(variable::InfiniteOpt.GeneralVariableRef) = + vec(vcat(JuMP.value(variable))) + +# Read per-support values from the transformation backend, keyed by +# InfiniteOpt vars. Skips fixed vars. The objective-side translation +# (transcribe-then-AD when the objective has aggregate refs) lives in +# the `objective_linearization_point(::InfiniteModel, ...)` override +# below — base `extract_solution` doesn't need to anticipate it. +function DP.extract_solution(model::InfiniteOpt.InfiniteModel) + return Dict(var => _per_support_values(var) + for var in DP.collect_all_vars(model) + if !JuMP.is_fixed(var)) +end + +# Add a pointwise-sum cut directly to the transformation backend and mark +# it ready so the next optimize! doesn't re-transcribe and wipe the cut. +function DP.add_cut( + model::InfiniteOpt.InfiniteModel, + decision_vars::Vector{InfiniteOpt.GeneralVariableRef}, + rBM_sol::Dict{<:JuMP.AbstractVariableRef, <:Vector{<:Number}}, + sep_sol::Dict{<:JuMP.AbstractVariableRef, <:Vector{<:Number}} + ) + transcribed = InfiniteOpt.transformation_model(model) + cut_expr = zero(JuMP.GenericAffExpr{ + JuMP.value_type(typeof(transcribed)), + JuMP.variable_ref_type(transcribed)}) + for var in decision_vars + haskey(rBM_sol, var) || continue + haskey(sep_sol, var) || continue + rbm_vals = rBM_sol[var] + sep_vals = sep_sol[var] + transcription_var = InfiniteOpt.transformation_variable(var) + transcribed_vars = transcription_var isa AbstractArray ? + vec(transcription_var) : [transcription_var] + for k in eachindex(transcribed_vars) + xi = 2 * (sep_vals[k] - rbm_vals[k]) + JuMP.add_to_expression!(cut_expr, xi, transcribed_vars[k]) + JuMP.add_to_expression!(cut_expr, -xi * sep_vals[k]) + end + end + JuMP.@constraint(transcribed, cut_expr >= 0) + InfiniteOpt.set_transformation_backend_ready(model, true) + return +end +################################################################################ +# LOA FOR INFINITEMODEL +################################################################################ +# Dispatch overrides for InfiniteModel. The base LOA in src/loa.jl is +# written for finite (scalar) models. The InfiniteOpt master and feas +# submodel are themselves InfiniteModels, with per-support handling +# via point evaluation on infinite `GeneralVariableRef`s. + +# Supports of the infinite parameter group `v` depends on. For a 1-D +# parameter, `supports(p)` is a `Vector{Float64}`; for a dependent +# group of dimension k, it is a k × N_supports `Matrix`. Returns each +# joint support point as one element — scalar for 1-D, column view for +# multi-D — so callers can iterate uniformly. Pair with `_at_support`, +# which splats vector supports into the variable's call form. +function _supports_of(v::InfiniteOpt.GeneralVariableRef) + p = only(InfiniteOpt.parameter_refs(v)) + sup = InfiniteOpt.supports(p) + ndims(sup) == 1 && return sup + # Materialize each column as a concrete `Vector{Float64}`; + # `eachcol` yields `SubArray` views that InfiniteOpt's + # `VectorTuple` constructor refuses. + return [sup[:, k] for k in axes(sup, 2)] +end + +_is_point_var(v::InfiniteOpt.GeneralVariableRef) = + InfiniteOpt.dispatch_variable_ref(v) isa InfiniteOpt.PointVariableRef + +# A `MeasureRef` collapses `f(x, t)` over an infinite parameter into +# one ref; a `ParameterFunctionRef` collapses a parameter-dependent +# function value into one ref. Either hides decision-variable +# dependence from MOI Nonlinear AD, so LOA needs to transcribe the +# enclosing expression to recover correct gradients. +_is_aggregate_ref(v::InfiniteOpt.GeneralVariableRef) = + InfiniteOpt.dispatch_variable_ref(v) isa Union{ + InfiniteOpt.MeasureRef, InfiniteOpt.ParameterFunctionRef} + +function _has_aggregate_ref(expr) + found = Ref(false) + DP._interrogate_variables(expr) do v + found[] || (found[] = _is_aggregate_ref(v)) + end + return found[] +end + +# Resolve a `PointVariableRef` (e.g., `L(0)` from a boundary +# condition) to the master's corresponding point variable: look up +# the underlying infinite var in `var_map`, then point-evaluate at +# the same support values. For all other variable refs (parameters, +# infinite/finite decision vars), fall back to direct lookup. +function DP.replace_variables_in_constraint( + v::InfiniteOpt.GeneralVariableRef, var_map::AbstractDict + ) + if _is_point_var(v) + underlying = InfiniteOpt.infinite_variable_ref(v) + return var_map[underlying](InfiniteOpt.parameter_values(v)...) + end + return var_map[v] +end + +# Bool active arises from `_set_covering_combinations`, which keys +# combinations on `LogicalVariableRef -> Bool` regardless of whether +# the indicator is infinite. Broadcast over all supports for an +# infinite indicator; for a finite (or point-variable) ref, fall +# through to the base scalar dispatch. +function DP.add_no_good_terms( + cut, binary_ref::InfiniteOpt.GeneralVariableRef, active::Bool + ) + isempty(InfiniteOpt.parameter_refs(binary_ref)) && + return invoke(DP.add_no_good_terms, + Tuple{Any, Any, Bool}, cut, binary_ref, active) + for support in _supports_of(binary_ref) + DP.add_no_good_terms(cut, _at_support(binary_ref, support), active) + end + return +end + +# AbstractVector active arises from `combination_val` after a master +# solve, which yields per-support `Vector{Bool}`. +function DP.add_no_good_terms( + cut, binary_ref::InfiniteOpt.GeneralVariableRef, + actives::AbstractVector + ) + supports = _supports_of(binary_ref) + for (k, support) in enumerate(supports) + DP.add_no_good_terms( + cut, _at_support(binary_ref, support), actives[k]) + end + return +end + +# Complement-form AffExpr (`1 - y_underlying`) with per-support active +# descriptor — same fan-out as the variable-ref form, but the +# underlying var lives inside the AffExpr's terms. Use the underlying +# var's supports to drive the loop. +function DP.add_no_good_terms( + cut, binary_ref::JuMP.GenericAffExpr, + actives::AbstractVector + ) + underlying = only(keys(binary_ref.terms)) + underlying isa InfiniteOpt.GeneralVariableRef || return invoke( + DP.add_no_good_terms, + Tuple{Any, Any, AbstractVector}, cut, binary_ref, actives) + supports = _supports_of(underlying) + for (k, support) in enumerate(supports) + DP.add_no_good_terms( + cut, _at_support(binary_ref, support), actives[k]) + end + return +end + +# Fan out across supports if either the binary or the constraint +# expression involves an infinite variable; otherwise emit one +# un-sliced site. The chosen supports come from the first infinite +# variable found in either expression. +function _infinite_cut_info( + binary_ref, + active, + constraint_func, + linearization_point, + variable_map + ) + supports = _relevant_supports(binary_ref, constraint_func) + supports === nothing && + return ((binary_ref, linearization_point, variable_map),) + actives = active isa AbstractVector ? active : + fill(active, length(supports)) + return _cut_sites(binary_ref, supports, actives, + linearization_point, variable_map) +end + +# Supports governing per-support fan-out — pick the first infinite +# variable found in the binary expression, falling back to the +# constraint expression. `nothing` means everything is finite, so the +# caller emits a single un-sliced site. +function _relevant_supports(binary_ref, constraint_func) + var = _find_infinite_var(binary_ref) + var === nothing && (var = _find_infinite_var(constraint_func)) + var === nothing && return nothing + return _supports_of(var) +end + +_find_infinite_var(v::InfiniteOpt.GeneralVariableRef) = + isempty(InfiniteOpt.parameter_refs(v)) ? nothing : v + +function _find_infinite_var(expr) + found = Ref{Any}(nothing) + DP._interrogate_variables(expr) do v + found[] === nothing || return + v isa InfiniteOpt.GeneralVariableRef || return + isempty(InfiniteOpt.parameter_refs(v)) && return + found[] = v + end + return found[] +end + +function _cut_sites( + binary_ref, + supports::AbstractVector, + actives::AbstractVector, + linearization_point::AbstractDict, + variable_map::AbstractDict + ) + sites = Any[] + for (k, support) in enumerate(supports) + actives[k] || continue + point_var_map = Dict{ + InfiniteOpt.GeneralVariableRef, + InfiniteOpt.GeneralVariableRef + }() + for (variable, master_var) in variable_map + point_var_map[variable] = _at_support(master_var, support) + end + point = Dict{InfiniteOpt.GeneralVariableRef, Float64}() + for (variable, point_value) in linearization_point + point[variable] = _at(point_value, k) + end + push!(sites, ( + _at_support(binary_ref, support), + point, + point_var_map + )) + end + return sites +end + +# Slice a per-support container at index `k`. Length-1 vectors are +# finite-shape (one value applies at every support, per +# `_per_support_values` wrapping `JuMP.value` of a finite var as +# `[scalar]`) — return the single value regardless of `k`. Scalars +# pass through. +_at(values::AbstractArray, k::Integer) = + length(values) == 1 ? values[1] : values[k] +_at(scalar, ::Integer) = scalar + +# Point-evaluate an InfiniteOpt var at `support` if it's infinite; +# return the var as-is if it's finite. `support` is one joint support +# point — a scalar for 1-D parameters, a vector for a multi-D +# dependent group — and `v(support)` matches both call forms. +_at_support(v::InfiniteOpt.GeneralVariableRef, support) = + isempty(InfiniteOpt.parameter_refs(v)) ? v : v(support) + +# Rebuild a complement-form AffExpr (`1 - y_underlying`) with its +# variables point-evaluated, so per-support fan-out can use the +# AffExpr directly in the gating term `M(1 - binary_at_support)`. +function _at_support( + expr::JuMP.GenericAffExpr{C, V}, support + ) where {C, V <: InfiniteOpt.GeneralVariableRef} + result = JuMP.GenericAffExpr{C, V}(expr.constant) + for (var, coef) in expr.terms + JuMP.add_to_expression!(result, coef, _at_support(var, support)) + end + return result +end + +# Map transcribed input JuMP var -> master point variable, used as +# `objective_ref_map`. Needed because MOI AD can't walk a `MeasureRef`, +# so a measure objective is transcribed, AD'd flat, and the gradient +# mapped back to master point vars. Per-support disjunct cuts skip this +# (they linearize natively via `_infinite_cut_info`). +function _transcribed_to_master_point( + model::InfiniteOpt.InfiniteModel, + ref_map::AbstractDict + ) + result = Dict{JuMP.VariableRef, InfiniteOpt.GeneralVariableRef}() + for v in DP.collect_all_vars(model) + # Point vars share the infinite var's transcription, so skipping + # them loses no mappings. + _is_point_var(v) && continue + master_var = ref_map[v] + transcribed = InfiniteOpt.transformation_variable(v) + prefs = InfiniteOpt.parameter_refs(v) + if isempty(prefs) + result[transcribed] = master_var + else + # `_supports_of`/`_at_support` handle 1-D (scalar) and + # multi-D (vector) supports; `vec(supports(...))` would + # flatten multi-D and break point evaluation. + for (k, support) in enumerate(_supports_of(v)) + result[vec(transcribed)[k]] = + _at_support(master_var, support) + end + end + end + return result +end + +# Build the master from `JuMP.copy_model(model)`, keeping linear +# constraints and the InfiniteOpt vars (per-support handling is done +# downstream by point evaluation). The objective branches on +# `_has_aggregate_ref`: aggregate-free objectives AD the InfiniteOpt +# expression directly; aggregate ones are transcribed flat and mapped +# back via `_transcribed_to_master_point`. +function DP.build_loa_master( + model::InfiniteOpt.InfiniteModel, + method::DP.LOA, + sink = nothing + ) + variable_type = InfiniteOpt.GeneralVariableRef + master, copy_ref_map = JuMP.copy_model( + model; + filter_constraints = function (cref) + con = JuMP.constraint_object(cref) + F = typeof(con.func) + F === variable_type && return true + DP._is_linear_F(F) || return false + return !_has_aggregate_ref(con.func) + end + ) + JuMP.set_optimize_hook(master, nothing) + JuMP.set_optimizer(master, method.mip_optimizer) + JuMP.set_silent(master) + + ref_map = Dict{InfiniteOpt.GeneralVariableRef, + InfiniteOpt.GeneralVariableRef}() + for v in DP.collect_all_vars(model) + _is_point_var(v) && continue + ref_map[v] = copy_ref_map[v] + end + for p in InfiniteOpt.all_parameters(model) + ref_map[p] = copy_ref_map[p] + end + for pfunc in InfiniteOpt.all_parameter_functions(model) + ref_map[pfunc] = copy_ref_map[pfunc] + end + + raw_objective = JuMP.objective_function(model) + objective_sense = JuMP.objective_sense(model) + if _has_aggregate_ref(raw_objective) + InfiniteOpt.build_transformation_backend!(model) + transcribed_input = InfiniteOpt.transformation_model(model) + original_objective = JuMP.objective_function(transcribed_input) + objective_ref_map = _transcribed_to_master_point( + model, ref_map) + else + original_objective = raw_objective + objective_ref_map = ref_map + end + + alpha_oa = JuMP.@variable(master, base_name = "alpha_oa") + JuMP.@objective(master, objective_sense, alpha_oa) + + binary_map = Dict{DP.LogicalVariableRef, Any}() + for (indicator, binary_ref) in DP._indicator_to_binary(model) + binary_map[indicator] = DP._remap_indicator_to_binary( + binary_ref, ref_map) + end + variable_map = Dict{InfiniteOpt.GeneralVariableRef, + InfiniteOpt.GeneralVariableRef}() + for v in DP.collect_all_vars(model) + _is_point_var(v) && continue + variable_map[v] = ref_map[v] + end + + # Aggregate-wrapped LINEAR constraints (e.g. `E(W, e) ≥ alpha`) were + # dropped at copy time and the OA-cut path skips linear `F`, so add + # each one to the master directly as a transcribed flat scalar. + _add_aggregate_linear_constraints( + master, model, ref_map, variable_type) + + return DP._LOAMaster(master, binary_map, variable_map, + objective_sense, original_objective, alpha_oa, + objective_ref_map, DP._build_disaggregator( + model, method.inner_method, variable_map, binary_map, sink)) +end + +function DP.record_disaggregation( + hull::DP._Hull, + ::InfiniteOpt.InfiniteModel, + variable, + binary, + disaggregated + ) + supports = _disaggregator_supports(variable, binary) + if supports === nothing + hull.disjunct_variables[(variable, binary)] = disaggregated + else + for support in supports + hull.disjunct_variables[(_at_support(variable, support), + _at_support(binary, support))] = + _at_support(disaggregated, support) + end + end + return +end + +# Supports the per-support fan-out gates this disjunct on (mirrors +# `_relevant_supports`): the binary's if the indicator is infinite, else +# the variable's if it is infinite, else `nothing` (both finite, one +# unsliced key). Keying off the binary lets a finite variable under an +# infinite indicator match the sliced binary the emitter looks up. +function _disaggregator_supports(master_var, master_bin) + binary_var = _find_infinite_var(master_bin) + binary_var === nothing || return _supports_of(binary_var) + isempty(InfiniteOpt.parameter_refs(master_var)) && return nothing + return _supports_of(master_var) +end + +# Transcribe each linear constraint containing a `MeasureRef` / +# `ParameterFunctionRef` and append the flat scalar form to the master. +# Already affine post-transcription, so `_linearize_at` just substitutes +# variables via `transcribed_to_master`. +function _add_aggregate_linear_constraints( + master::InfiniteOpt.InfiniteModel, + model::InfiniteOpt.InfiniteModel, + ref_map::AbstractDict, + variable_type::Type + ) + has_any = false + for (F, S) in JuMP.list_of_constraint_types(model) + F === variable_type && continue + DP._is_linear_F(F) || continue + for cref in JuMP.all_constraints(model, F, S) + con = JuMP.constraint_object(cref) + _has_aggregate_ref(con.func) || continue + has_any = true + break + end + has_any && break + end + has_any || return + + InfiniteOpt.build_transformation_backend!(model) + transcribed_to_master = _transcribed_to_master_point(model, ref_map) + # Finite parameters survive transcription as scalar JuMP vars and can + # appear in transcribed constraints, but `_transcribed_to_master_point` + # only walks decision vars. Map them to the master's parameters so the + # constraint stays parameter-relative (honors `set_value` without + # rebuild). + for p in InfiniteOpt.all_parameters(model) + InfiniteOpt.dispatch_variable_ref(p) isa + InfiniteOpt.FiniteParameterRef || continue + transcribed_p = InfiniteOpt.transformation_variable(p) + transcribed_p isa JuMP.VariableRef || continue + haskey(ref_map, p) || continue + transcribed_to_master[transcribed_p] = ref_map[p] + end + # `_linearize_at(::GenericAffExpr, ...)` ignores its + # `linearization_point` arg; pass any empty dict. + empty_point = Dict{JuMP.VariableRef, Float64}() + for (F, S) in JuMP.list_of_constraint_types(model) + F === variable_type && continue + DP._is_linear_F(F) || continue + for cref in JuMP.all_constraints(model, F, S) + con = JuMP.constraint_object(cref) + _has_aggregate_ref(con.func) || continue + con isa JuMP.ScalarConstraint || continue + transcribed_func = InfiniteOpt.transformation_expression( + con.func) + if transcribed_func isa AbstractArray + for tf in vec(transcribed_func) + master_expr = DP._linearize_at( + tf, empty_point, transcribed_to_master) + JuMP.@constraint(master, master_expr in con.set) + end + else + master_expr = DP._linearize_at( + transcribed_func, empty_point, transcribed_to_master) + JuMP.@constraint(master, master_expr in con.set) + end + end + end + return +end + +# Objective linearization point for `InfiniteModel` (the seam the shared +# base `add_oa_cuts` calls). The master objective is either the raw +# InfiniteOpt expression (non-aggregate, expects a scalar `xk[v]` keyed +# on InfiniteOpt vars) or the transcribed flat scalar expression +# (aggregate, expects a transcribed-`JuMP.VariableRef`-keyed scalar +# dict). Translate the base per-support `Vector`-valued point into +# whichever shape `original_objective` needs. +function DP.objective_linearization_point( + model::InfiniteOpt.InfiniteModel, + linearization_point::AbstractDict + ) + _has_aggregate_ref(JuMP.objective_function(model)) && + return _transcribe_linearization_point(model, linearization_point) + T = eltype(valtype(linearization_point)) + return Dict{InfiniteOpt.GeneralVariableRef, T}( + var => values[1] + for (var, values) in linearization_point + if isempty(InfiniteOpt.parameter_refs(var))) +end + +# Global OA cuts for `InfiniteModel`: route through transcription so +# per-support / aggregate-ref expressions reach `_linearize_at` as flat +# scalars over `JuMP.VariableRef`s. Aggregate-containing constraints +# (e.g. `∫f(x,t)dt ≤ 0`) transcribe to a single scalar expression. +# Constraints with infinite-parameter dependence (e.g. `x(t) ≥ 0`) +# transcribe to a per-support `AbstractArray` of scalar expressions — +# one OA cut is emitted per support. +function DP.add_global_oa_cuts( + model::InfiniteOpt.InfiniteModel, + master::DP._LOAMaster, + result::NamedTuple, + method::DP.LOA + ) + penalty_sign = DP._penalty_sign(Val(master.objective_sense)) + variable_type = InfiniteOpt.GeneralVariableRef + reform_set = DP.is_gdp_model(model) ? + Set(DP._reformulation_constraints(model)) : Set() + transcribed_xk = Ref{Any}(nothing) + transcribed_to_master = Ref{Any}(nothing) + ensure_transcribed = function () + transcribed_xk[] === nothing || return + InfiniteOpt.build_transformation_backend!(model) + transcribed_xk[] = _transcribe_linearization_point( + model, result.linearization_point) + transcribed_to_master[] = _transcribed_to_master_point( + model, master.variable_map) + return + end + for (F, S) in JuMP.list_of_constraint_types(model) + F === variable_type && continue + DP._is_linear_F(F) && continue + for cref in JuMP.all_constraints(model, F, S) + cref in reform_set && continue + con = JuMP.constraint_object(cref) + con isa JuMP.ScalarConstraint || continue + ensure_transcribed() + transcribed_func = InfiniteOpt.transformation_expression( + con.func) + if transcribed_func isa AbstractArray + for tf in vec(transcribed_func) + lin = DP._linearize_at(tf, transcribed_xk[], + transcribed_to_master[]) + DP._add_global_oa_row(master, lin, con.set, + method, penalty_sign) + end + else + lin = DP._linearize_at(transcribed_func, + transcribed_xk[], transcribed_to_master[]) + DP._add_global_oa_row(master, lin, con.set, + method, penalty_sign) + end + end + end + return +end + +# Override the disjunct-OA-cut pass for `InfiniteModel`: own the per-pass +# transcription `cache` and delegate each active disjunct constraint to +# the emitter below. Reuses the base active-disjunct walk so only the +# per-constraint emission differs from base. +function DP.add_disjunct_oa_cuts( + model::InfiniteOpt.InfiniteModel, + master::DP._LOAMaster, + result::NamedTuple, + method::DP.LOA + ) + penalty_sign = DP._penalty_sign(Val(master.objective_sense)) + cache = Ref{Any}(nothing) + DP._each_active_disjunct_constraint(model, master, + result) do binary_ref, active, constraint + _add_infinite_disjunct_constraint_oa_cuts(model, constraint, + master, binary_ref, active, result, method, penalty_sign, + cache) + end + return +end + +# Emit the OA cut(s) for one active disjunct constraint on an +# `InfiniteModel`. Non-aggregate constraints fan out per support via +# `_infinite_cut_info`; aggregate ones (`MeasureRef`) are transcribed flat +# and handed to `_add_oa_cut_for_constraint`. `cache` memoizes the +# transcription maps once per pass (lazily, on the first aggregate +# constraint). +function _add_infinite_disjunct_constraint_oa_cuts( + model, constraint, master, binary_ref, active, result, method, + penalty_sign, cache + ) + if _has_aggregate_ref(constraint.func) + if cache[] === nothing + InfiniteOpt.build_transformation_backend!(model) + cache[] = (_transcribe_linearization_point( + model, result.linearization_point), + _transcribed_to_master_point(model, master.variable_map)) + end + transcribed_xk, transcribed_to_master = cache[] + transcribed_func = InfiniteOpt.transformation_expression( + constraint.func) + transcribed_constraint = JuMP.ScalarConstraint( + transcribed_func, constraint.set) + for (point_binary, _, _) in _infinite_cut_info(binary_ref, + active, transcribed_func, result.linearization_point, + master.variable_map) + DP._add_oa_cut_for_constraint(transcribed_constraint, master, + point_binary, transcribed_xk, transcribed_to_master, + method, penalty_sign) + end + return + end + for (point_binary, point, var_map) in _infinite_cut_info(binary_ref, + active, constraint.func, result.linearization_point, + master.variable_map) + DP._add_oa_cut_for_constraint(constraint, master, point_binary, + point, var_map, method, penalty_sign) + end + return +end + +# Apply the combination's fixes and return a closure that reverses them. +# A `Bool` delegates to base `fix_indicator`; a per-support `Vector{Bool}` +# pins each support of the infinite binary with a point-equality +# constraint. The closure (used by `with_fixed_combination` to tear down +# each iteration) deletes those pins and unfixes; `commit_combination` +# ignores it to persist the committed fix. +function DP.fix_combination( + model::InfiniteOpt.InfiniteModel, combination::AbstractDict + ) + constraint_refs = InfiniteOpt.InfOptConstraintRef[] + fixed_indicators = DP.LogicalVariableRef{ + InfiniteOpt.InfiniteModel}[] + for (indicator, value) in combination + if value isa AbstractVector + binary_ref = DP._indicator_to_binary(model)[indicator] + target, target_for_active = binary_ref isa JuMP.GenericAffExpr ? + (only(keys(binary_ref.terms)), 0.0) : + (binary_ref, 1.0) + target_for_inactive = 1.0 - target_for_active + for (k, support) in enumerate(_supports_of(target)) + push!(constraint_refs, JuMP.@constraint(model, + _at_support(target, support) == + (value[k] ? target_for_active : target_for_inactive))) + end + else + DP.fix_indicator(model, indicator, value) + push!(fixed_indicators, indicator) + end + end + return function () + for ref in constraint_refs + JuMP.is_valid(model, ref) && JuMP.delete(model, ref) + end + for indicator in fixed_indicators + DP.unfix_indicator(model, indicator) + end + end +end + +# Per-support fix for an infinite logical's `Vector{Bool}` value: an +# equality at each support. `with_fixed_combination` undid the original's +# fix before the NLPF fall-through, so this is the only fix on the copy. +# Equalities (not `JuMP.fix`) avoid force-deleting the relaxed copy's +# bounds, whose refs may not survive `copy_model`. +function DP._fix_binary_on_copy(copy, binary, value::AbstractVector) + underlying = binary isa JuMP.GenericAffExpr ? + only(keys(binary.terms)) : binary + for (k, support) in enumerate(_supports_of(underlying)) + JuMP.@constraint(copy, + _at_support(binary, support) == (value[k] ? 1.0 : 0.0)) + end + return +end + +# Broadcast a per-support warm start across an infinite var's +# transcribed instances; for a finite var on an InfiniteModel, +# `transcription_variable` returns a single ref and `values` is +# length-1. +function DP.set_linearization_start( + variable::InfiniteOpt.GeneralVariableRef, + values::AbstractVector + ) + transcribed = InfiniteOpt.transformation_variable(variable) + if transcribed isa AbstractArray + for (k, ref) in enumerate(vec(transcribed)) + JuMP.set_start_value(ref, values[k]) + end + else + JuMP.set_start_value(transcribed, only(values)) + end + return +end + +# Convert an InfiniteModel-var-keyed per-support point into a +# transcribed-JuMP-var-keyed scalar point. Companion to +# `_transcribed_to_master_point`: feeds the AD walker for the +# objective OA cut. See the WHY note above +# `_transcribed_to_master_point` for the MeasureRef/AD constraint +# that motivates this whole transcribe-then-AD layer. +function _transcribe_linearization_point( + model::InfiniteOpt.InfiniteModel, + linearization_point::AbstractDict + ) + T = eltype(valtype(linearization_point)) + transcribed = Dict{JuMP.VariableRef, T}() + for (variable, values) in linearization_point + _add_to_transcribed_dict( + transcribed, + InfiniteOpt.transformation_variable(variable), + values + ) + end + return transcribed +end + +# Infinite var: per-support transcribed array, values per-support +function _add_to_transcribed_dict( + d::AbstractDict, + ts::AbstractArray, + values::AbstractVector + ) + for (k, ref) in enumerate(vec(ts)) + d[ref] = values[k] + end + return +end + +# Finite var: single transcribed ref, values 1-element +_add_to_transcribed_dict( + d::AbstractDict, + ts::JuMP.AbstractVariableRef, + values::AbstractVector + ) = + (d[ts] = values[1]; nothing) + +################################################################################ + # Build CP subproblem: reformulate the InfiniteModel in-place, transcribe, # copy, and wrap in GDPSubmodel with forward variable map. function DP.copy_and_reformulate( diff --git a/src/DisjunctiveProgramming.jl b/src/DisjunctiveProgramming.jl index 85f03455..6f1537e7 100644 --- a/src/DisjunctiveProgramming.jl +++ b/src/DisjunctiveProgramming.jl @@ -28,6 +28,7 @@ include("print.jl") include("extension_api.jl") include("utilities.jl") include("psplit.jl") +include("loa.jl") # Define additional stuff that should not be exported const _EXCLUDE_SYMBOLS = [Symbol(@__MODULE__), :eval, :include] diff --git a/src/cuttingplanes.jl b/src/cuttingplanes.jl index dbc818e5..63692258 100644 --- a/src/cuttingplanes.jl +++ b/src/cuttingplanes.jl @@ -8,25 +8,17 @@ function collect_cutting_planes_vars(model::JuMP.AbstractModel) return collect_all_vars(model) end -# Extract solution from a solved model (in-place). Extensions -# override for models where values live on a backend. +# Read primal values from a solved model. Returns +# `Dict{var, Vector{value}}` — per-support shape uniformly: finite +# models trivially have one "support" (length-1 Vector), the +# InfiniteOpt extension overrides this dispatch to populate +# multi-support Vectors. Skips fixed vars. function extract_solution(model::JuMP.AbstractModel) dvars = collect_cutting_planes_vars(model) V = eltype(dvars) T = JuMP.value_type(typeof(model)) return Dict{V, Vector{T}}( - v => [JuMP.value(v)] for v in dvars) -end - -# Extract solution from a GDPSubmodel (SEP path). -function extract_solution(sub::GDPSubmodel) - V = eltype(sub.decision_vars) - T = JuMP.value_type(typeof(sub.model)) - sol = Dict{V, Vector{T}}() - for var in sub.decision_vars - sol[var] = JuMP.value.(sub.fwd_map[var]) - end - return sol + v => [JuMP.value(v)] for v in dvars if !JuMP.is_fixed(v)) end # Set quadratic separation objective: min Σ (x_k - rBM_k)². diff --git a/src/datatypes.jl b/src/datatypes.jl index bdcd4c2e..4d8c9208 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -417,8 +417,11 @@ constraints. """ struct Hull{T} <: AbstractReformulationMethod value::T - function Hull(ϵ::T = 1e-6) where {T} - new{T}(ϵ) + # Internal: LOA installs a Dict here to collect the disaggregation + # map during reformulation; `nothing` for every other use. + sink::Any + function Hull(ϵ::T = 1e-6; sink = nothing) where {T} + new{T}(ϵ, sink) end end @@ -427,11 +430,13 @@ mutable struct _Hull{V <: JuMP.AbstractVariableRef, T} <: AbstractReformulationM value::T disjunction_variables::Dict{V, Vector{V}} disjunct_variables::Dict{Tuple{V, Union{V, JuMP.GenericAffExpr{T, V}}}, V} + sink::Any function _Hull(method::Hull{T}, vrefs::Set{V}) where {T, V <: JuMP.AbstractVariableRef} new{V, T}( method.value, - Dict{V, Vector{V}}(vref => V[] for vref in vrefs), - Dict{Tuple{V, Union{V, JuMP.GenericAffExpr{T, V}}}, V}() + Dict{V, Vector{V}}(vref => V[] for vref in vrefs), + Dict{Tuple{V, Union{V, JuMP.GenericAffExpr{T, V}}}, V}(), + method.sink ) end end @@ -473,25 +478,24 @@ end """ GDPSubmodel{M, V, W} -A unified submodel wrapper used by MBM and cutting plane -reformulations. It encapsulates a flat JuMP optimization -submodel built from a single disjunct's feasible region, -along with mappings back to the original model's variables. +A JuMP submodel paired with the parent model's decision variables +and a map from each parent variable to its submodel representation. ## Fields -- `model::M`: The JuMP submodel representing a disjunct's - feasible region (constraints and variable bounds). -- `decision_vars::Vector{V}`: Ordered decision variables in - the submodel, matching the original model's ordering. -- `fwd_map::Dict{V, Vector{W}}`: Forward map from original - model variables to their submodel counterparts. +- `model::M`: The JuMP submodel. +- `decision_vars::Vector{V}`: Parent-model decision variables in + the order used by the submodel. +- `fwd_map::Dict{V, W}`: Map from each parent-model variable to + its submodel counterpart. `W` can be a single variable + reference (1:1 copy) or a collection of references (e.g. one + per transcription support). """ struct GDPSubmodel{M <: JuMP.AbstractModel, V <: JuMP.AbstractVariableRef, - W <: JuMP.AbstractVariableRef} + W} model::M decision_vars::Vector{V} - fwd_map::Dict{V, Vector{W}} + fwd_map::Dict{V, W} end """ @@ -567,6 +571,82 @@ A type for using indicator constraint approach for linear disjunctive constraint """ struct Indicator <: AbstractReformulationMethod end +################################################################################ +# LOA +################################################################################ +""" + LOA{O, P, R, T} <: AbstractReformulationMethod + +Logic-based Outer Approximation solver for GDP models. Iterates a primary +NLP (original model reformulated by `inner_method`, binaries fixed per +iteration) and a master MILP accumulating OA and no-good cuts. +`inner_method` is `BigM` (default), `MBM`, or `Hull`. + +## Fields +- `nlp_optimizer::O`: solver for the primary NLP. +- `mip_optimizer::P`: solver for the master MILP (default `nlp_optimizer`). +- `inner_method::R`: NLP reformulation — `BigM`, `MBM`, or `Hull`. +- `max_iter::Int`: max iterations after set-covering seeding. +- `M_value::T`: big-M for the disjunct OA cut gating term. +- `max_slack::T`: upper bound per slack variable. +- `oa_penalty::T`: penalty on slacks in the master objective. +- `convergence_tol::Float64`: relative gap tolerance for the early stop. +- `slack_tol::Float64`: max total slack for which the bound still counts + as converged (positive slack = nonconvex crossing, keep iterating). +- `iteration_time_limit::Float64`: budget (s) for the iteration loop. +- `time_limit::Float64`: overall budget (s) incl. the final solve + (default 3600; `Inf` disables). +""" +struct LOA{O, P, R, T} <: AbstractReformulationMethod + nlp_optimizer::O + mip_optimizer::P + inner_method::R + max_iter::Int + M_value::T + max_slack::T + oa_penalty::T + convergence_tol::Float64 + slack_tol::Float64 + iteration_time_limit::Float64 + time_limit::Float64 + function LOA( + nlp_optimizer::O; + mip_optimizer::P = nlp_optimizer, + max_iter::Int = 10, + M_value::T = 1e9, + max_slack::T = 1e3, + oa_penalty::T = 1e3, + inner_method::R = BigM(M_value), + convergence_tol::Float64 = 1e-6, + slack_tol::Float64 = 1e-4, + iteration_time_limit::Float64 = Inf, + time_limit::Float64 = 3600.0 + ) where {O, P, R <: AbstractReformulationMethod, T} + R <: Union{BigM, MBM, Hull} || error( + "LOA inner_method must be BigM, MBM, or Hull (got $R). " * + "PSplit is not yet supported.") + new{O, P, R, T}(nlp_optimizer, mip_optimizer, inner_method, + max_iter, M_value, max_slack, oa_penalty, + convergence_tol, slack_tol, + iteration_time_limit, time_limit) + end +end + +# The LOA master MILP plus maps from original- to master-model refs. +# `objective_ref_map` splits from `variable_map` so an extension can map +# objective vars separately (identical in base); `disaggregator` is the +# master `_Hull` for Hull cuts, `nothing` for Big-M / MBM. +mutable struct _LOAMaster{M <: JuMP.AbstractModel, OF, AO, BM, VM, RM, DG} + model::M + binary_map::BM + variable_map::VM + objective_sense::_MOI.OptimizationSense + original_objective::OF + alpha_oa::AO + objective_ref_map::RM + disaggregator::DG +end + ################################################################################ # GDP Data ################################################################################ diff --git a/src/hull.jl b/src/hull.jl index 779369e9..9075799f 100644 --- a/src/hull.jl +++ b/src/hull.jl @@ -42,6 +42,8 @@ function _disaggregate_variable( #temp storage push!(method.disjunction_variables[vref], dvref) method.disjunct_variables[vref, bvref] = dvref + #record into the LOA sink when one is installed (LOA-Hull only) + method.sink === nothing || (method.sink[(vref, lvref)] = dvref) #create bounding constraints dvname = JuMP.name(dvref) lbname = isempty(dvname) ? "" : "$(dvname)_lower_bound" @@ -254,7 +256,8 @@ function reformulate_disjunction(model::JuMP.AbstractModel, disj::Disjunction, m return ref_cons end function reformulate_disjunction(model::JuMP.AbstractModel, disj::Disjunction, method::_Hull) - return reformulate_disjunction(model, disj, Hull(method.value)) + return reformulate_disjunction(model, disj, + Hull(method.value; sink = method.sink)) end function reformulate_disjunct_constraint( diff --git a/src/loa.jl b/src/loa.jl new file mode 100644 index 00000000..14d6e629 --- /dev/null +++ b/src/loa.jl @@ -0,0 +1,760 @@ +################################################################################ +# LOGIC-BASED OUTER APPROXIMATION (LOA) +################################################################################ + +################################################################################ +# SENSE HANDLERS +################################################################################ +_penalty_sign(::Val{_MOI.MIN_SENSE}) = 1 +_penalty_sign(::Val{_MOI.MAX_SENSE}) = -1 +_worst_objective(::Val{_MOI.MIN_SENSE}) = Inf +_worst_objective(::Val{_MOI.MAX_SENSE}) = -Inf +_is_better(::Val{_MOI.MIN_SENSE}, new, best) = new < best +_is_better(::Val{_MOI.MAX_SENSE}, new, best) = new > best +_gap(::Val{_MOI.MIN_SENSE}, best, bound) = best - bound +_gap(::Val{_MOI.MAX_SENSE}, best, bound) = bound - best + +################################################################################ +# MAIN LOOP +################################################################################ +# LOA entry: set-covering seeds, the master/NLP loop, commit the best combo +function reformulate_model(model::JuMP.AbstractModel, method::LOA) + _clear_reformulations(model) + combinations = _set_covering_combinations(model) + # Hull needs the disaggregation map if its used as inner function + inner, sink = _loa_inner_method(model, method.inner_method) + reformulate_model(model, inner) + + master = build_loa_master(model, method, sink) + # The original is only ever the NLP (build_loa_master keeps the + # master's binaries), so relax it once; fixing overwrites in place. + relax_logical_vars(model) + JuMP.set_optimizer(model, method.nlp_optimizer) + JuMP.set_silent(model) + t_start = time() + overall_deadline = t_start + method.time_limit + loop_deadline = + min(t_start + method.iteration_time_limit, overall_deadline) + original_time_limit = JuMP.time_limit_sec(model) + sense_token = Val(JuMP.objective_sense(model)) + best_objective = _worst_objective(sense_token) + best_result = nothing + master_bound = nothing + + # Seed the master with an OA cut per set-covering combination, each + # NLP warm-started from the last feasible primal. + previous_result = nothing + for combination in combinations + time() < loop_deadline || break + _set_nlp_warm_start(previous_result) + result = _solve_nlp(model, combination, method; + deadline = loop_deadline) + avoid_combination(master.model, combination, master.binary_map) + add_oa_cuts(model, master, result, method) + if result.feasible && + _is_better(sense_token, result.objective, best_objective) + best_objective = result.objective + best_result = result + end + result.feasible && (previous_result = result) + end + + # Master/NLP loop: `alpha_oa` is the bound, the NLP refines the + # incumbent. Exit on convergence, infeasible master, max_iter, or time. + converged = false + for _ in 1:method.max_iter + time() < loop_deadline || break + _cap_remaining_time(master.model, loop_deadline) + JuMP.optimize!(master.model) + JuMP.is_solved_and_feasible(master.model) || break + master_bound = JuMP.value(master.alpha_oa) + if best_result !== nothing + gap = _gap(sense_token, best_objective, master_bound) + total_slack = abs(JuMP.objective_value(master.model) - + master_bound) / method.oa_penalty + tol = method.convergence_tol * max(abs(best_objective), 1.0) + if gap <= tol && total_slack <= method.slack_tol + converged = true + break + end + end + combination = _extract_combination(model, master) + _set_nlp_warm_start(previous_result) + result = _solve_nlp(model, combination, method; + deadline = loop_deadline) + avoid_combination(master.model, combination, master.binary_map) + add_oa_cuts(model, master, result, method) + if result.feasible && + _is_better(sense_token, result.objective, best_objective) + best_objective = result.objective + best_result = result + end + result.feasible && (previous_result = result) + end + + if best_result !== nothing + commit_combination(model, best_result.combination, + best_result.linearization_point) + end + if isfinite(method.time_limit) + # Overall cap: the final committed solve gets the budget left. + JuMP.set_time_limit_sec(model, max(0.0, overall_deadline - time())) + elseif isfinite(method.iteration_time_limit) + # Loop-only budget: restore so the final solve isn't crippled. + _restore_time_limit(model, original_time_limit) + end + _report_loa_gap(best_objective, best_result, master_bound, sense_token, + converged) + _set_solution_method(model, method) + _set_ready_to_optimize(model, true) + return +end + +# Report the final gap. The bound is rigorous only for a convex inner +# problem; on a nonconvex one it can cross the incumbent (negative gap). +function _report_loa_gap( + best_objective::Real, + best_result, + master_bound, + sense_token::Val, + converged::Bool + ) + if best_result === nothing + @warn "LOA finished: no feasible incumbent found." + return + end + if master_bound === nothing + @info "LOA finished [no bound]: incumbent $best_objective " * + "(master produced no bound; best seed kept)." + return + end + gap = _gap(sense_token, best_objective, master_bound) + relative = abs(best_objective) > 1e-10 ? + gap / abs(best_objective) : gap + label = converged ? "converged" : "limit hit" + @info "LOA finished [$label]: incumbent $best_objective, master " * + "bound $master_bound, gap $gap (relative $relative)." + return +end + +# Error fallback for unsupported model types. +function reformulate_model(::M, ::LOA) where {M} + error("reformulate_model not implemented for model type `$(M)` with LOA.") +end + +################################################################################ +# SET-COVERING INITIALIZATION (simple version from pyomo) +################################################################################ +# `K = max disjunction size` combinations that activate every indicator +# at least once: combination `k` activates the `k`-th indicator of each +# disjunction, cycling via `mod1`. Inconsistent nested combinations are +# caught by the no-good cut from the infeasible NLP. +function _set_covering_combinations(model::JuMP.AbstractModel) + LogicalRef = LogicalVariableRef{typeof(model)} + indicator_lists = [collect(d.constraint.indicators) + for (_, d) in _disjunctions(model)] + isempty(indicator_lists) && return Dict{LogicalRef, Bool}[] + K = maximum(length, indicator_lists) + return [ + Dict{LogicalRef, Bool}( + indicator => (indicator == indicators[mod1(k, length(indicators))]) + for indicators in indicator_lists + for indicator in indicators) + for k in 1:K + ] +end + +################################################################################ +# MASTER CONSTRUCTION +################################################################################ +# True for linear constraint function types (variable refs / affine). +# Nonlinear functions enter the master as OA cuts, not at copy time. +_is_linear_F(::Type{<:JuMP.AbstractVariableRef}) = true +_is_linear_F(::Type{<:JuMP.GenericAffExpr}) = true +_is_linear_F(::Type{<:AbstractVector{<:JuMP.AbstractVariableRef}}) = true +_is_linear_F(::Type{<:AbstractVector{<:JuMP.GenericAffExpr}}) = true +_is_linear_F(::Type) = false + +""" + build_loa_master(model, method::LOA, sink = nothing)::_LOAMaster + +Build the LOA master MILP: copy `model`'s variables and linear +constraints, install the `alpha_oa` objective auxiliary, and record the +original-to-master reference maps. Nonlinear objective and disjunct +constraints are not copied; they enter later as OA cuts per NLP solve. +`sink` is the `(variable, indicator) -> disaggregated variable` map from +an inner Hull reformulation (`nothing` for Big-M / MBM). + +## Returns +- `_LOAMaster`: the master model with its reference maps and + disaggregator. +""" +function build_loa_master( + model::JuMP.AbstractModel, + method::LOA, + sink = nothing + ) + original_objective = JuMP.objective_function(model) + objective_sense = JuMP.objective_sense(model) + variable_type = JuMP.variable_ref_type(typeof(model)) + + master = _copy_model(model) + JuMP.set_optimizer(master, method.mip_optimizer) + JuMP.set_silent(master) + + variable_map = copy_variables_onto_model(master, model) + + for (F, S) in JuMP.list_of_constraint_types(model) + F === variable_type && continue + _is_linear_F(F) || continue + for constraint_ref in JuMP.all_constraints(model, F, S) + constraint = JuMP.constraint_object(constraint_ref) + new_func = replace_variables_in_constraint( + constraint.func, variable_map) + JuMP.@constraint(master, new_func in constraint.set) + end + end + + binary_map = Dict{LogicalVariableRef, Any}() + for (indicator, binary_ref) in _indicator_to_binary(model) + binary_map[indicator] = _remap_indicator_to_binary( + binary_ref, variable_map) + end + + alpha_oa = JuMP.@variable(master, base_name = "alpha_oa") + JuMP.@objective(master, objective_sense, alpha_oa) + + disaggregator = _build_disaggregator(model, method.inner_method, + variable_map, binary_map, sink) + return _LOAMaster(master, binary_map, variable_map, objective_sense, + original_objective, alpha_oa, variable_map, disaggregator) +end + +# The inner reformulation method LOA runs, plus the disaggregation sink to +# collect (Big-M / MBM need none). +#For Hull, return a sink-carrying copy so the reformulation records its `(variable, indicator) +# -> disaggregated variable` map into a fresh LOA-owned `Dict`. +_loa_inner_method(::JuMP.AbstractModel, inner::Union{BigM, MBM}) = + (inner, nothing) +function _loa_inner_method(model::JuMP.AbstractModel, inner::Hull) + V = JuMP.variable_ref_type(typeof(model)) + sink = Dict{Tuple{V, LogicalVariableRef{typeof(model)}}, V}() + return (Hull(inner.value; sink = sink), sink) +end + +# Master-space Hull disaggregator for the convex-hull cut emitter +# (`nothing` for Big-M / MBM). Remaps the collected `(variable, indicator) +# -> disaggregated variable` sink into master refs. +_build_disaggregator(::JuMP.AbstractModel, ::Union{BigM, MBM}, + variable_map, binary_map, sink) = nothing +function _build_disaggregator( + model::JuMP.AbstractModel, + inner::Hull, + variable_map::AbstractDict, + binary_map::AbstractDict, + sink + ) + hull = _Hull(inner, Set{valtype(variable_map)}()) + for ((variable, indicator), disaggregated) in sink + record_disaggregation(hull, model, variable_map[variable], + binary_map[indicator], variable_map[disaggregated]) + end + return hull +end + +""" + record_disaggregation(hull::_Hull, model, variable, binary, + disaggregated) + +Record one master-space disaggregation in the Hull `hull` used to gate +Hull OA cuts, keyed by `(variable, binary)`. +""" +function record_disaggregation( + hull::_Hull, + ::JuMP.AbstractModel, + variable, + binary, + disaggregated + ) + hull.disjunct_variables[(variable, binary)] = disaggregated + return +end + +################################################################################ +# NLP SUBPROBLEM +################################################################################ +function _cap_remaining_time(target::JuMP.AbstractModel, deadline::Float64) + isfinite(deadline) || return + JuMP.set_time_limit_sec(target, max(0.0, deadline - time())) + return +end +_restore_time_limit(model::JuMP.AbstractModel, ::Nothing) = + JuMP.unset_time_limit_sec(model) +_restore_time_limit(model::JuMP.AbstractModel, seconds::Real) = + JuMP.set_time_limit_sec(model, seconds) + +# Solve the primary NLP at a fixed combination. 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::M, + combination, + method::LOA; + deadline::Float64 = Inf + ) where {M <: JuMP.AbstractModel} + _cap_remaining_time(model, deadline) + primary = with_fixed_combination(model, combination) do + JuMP.optimize!(model, ignore_optimize_hook = true) + if JuMP.is_solved_and_feasible(model) + return (combination = combination, + linearization_point = extract_solution(model), + objective = JuMP.objective_value(model), feasible = true) + end + return nothing + end + primary === nothing || return primary + + # Primary NLP infeasible — try NLPF. + nlpf = _solve_nlpf(model, combination, method; deadline = deadline) + nlpf === nothing || return nlpf + + return (combination = combination, + linearization_point = Dict{JuMP.AbstractVariableRef, Any}(), + objective = Inf, feasible = false) +end + +################################################################################ +# NLPF (FEASIBILITY SUBPROBLEM) +################################################################################ +# When the primary NLP is infeasible, copy the model, slack every scalar +# inequality with one nonnegative `u`, minimize `u`, and return the +# primal as the linearization site (equalities/bounds exact, inequalities +# violated by at most `u`). +function _solve_nlpf( + model::M, + combination, + method::LOA; + deadline::Float64 = Inf + ) where {M <: JuMP.AbstractModel} + # The original's binaries are permanently relaxed, so this copy is + # continuous — any NLP solver (e.g. Ipopt) can solve NLPF. + copy, ref_map = JuMP.copy_model(model) + JuMP.set_optimizer(copy, method.nlp_optimizer) + JuMP.set_silent(copy) + _cap_remaining_time(copy, deadline) + + var_type = JuMP.variable_ref_type(typeof(copy)) + u = JuMP.@variable(copy, lower_bound = 0.0, base_name = "_nlpf_u") + + for (F, S) in JuMP.list_of_constraint_types(copy) + F === var_type && continue + _nlpf_should_slack(S) || continue + for cref in JuMP.all_constraints(copy, F, S) + con = JuMP.constraint_object(cref) + JuMP.delete(copy, cref) + new_func = _nlpf_slacked_func(con.func, u, con.set) + JuMP.@constraint(copy, new_func in con.set) + end + end + + JuMP.@objective(copy, Min, u) + + # Fix the combination on the copy so NLPF solves at the chosen + # binaries, not by leaning on the copy inheriting the original's fix. + fix_combination_on_copy(copy, model, combination, ref_map) + + JuMP.optimize!(copy, ignore_optimize_hook = true) + # Use the primal only at a genuine feasible point; a solver can report + # `has_values` with a nonfeasible/NaN primal that would poison the cut. + JuMP.is_solved_and_feasible(copy) || return nothing + + linearization_point = _nlpf_extract_primal(model, ref_map) + return (combination = combination, + linearization_point = linearization_point, + objective = Inf, feasible = false) +end + +# Fix each indicator's binary on the NLPF copy at its combination value, +# mapped into copy space. Explicit so NLPF is self-contained rather than +# relying on the copy inheriting the original's fix. +function fix_combination_on_copy(copy, model, combination, ref_map) + for (indicator, value) in combination + binary = _binary_on_copy( + _indicator_to_binary(model)[indicator], ref_map) + _fix_binary_on_copy(copy, binary, value) + end + return +end + +# Map an original binary (or its `1 - y` complement) into NLPF-copy space. +_binary_on_copy(binary::JuMP.AbstractVariableRef, ref_map) = ref_map[binary] +_binary_on_copy(binary::JuMP.GenericAffExpr, ref_map) = + 1.0 - ref_map[only(keys(binary.terms))] + +# OVERRIDABLE. Fix one binary on the copy at a scalar value via an +# equality (the copy is discarded, so the fix needs no teardown). The +# InfiniteOpt extension adds a per-support method for a `Vector{Bool}`. +_fix_binary_on_copy(copy, binary, value::Bool) = + JuMP.@constraint(copy, binary == (value ? 1.0 : 0.0)) + +_nlpf_should_slack(::Type{<:_MOI.LessThan}) = true +_nlpf_should_slack(::Type{<:_MOI.GreaterThan}) = true +_nlpf_should_slack(::Type) = false + +_nlpf_slacked_func(func, u, ::_MOI.LessThan) = func - u +_nlpf_slacked_func(func, u, ::_MOI.GreaterThan) = func + u + +# Read primal values from `copy` keyed by the original model's +# variables, in the same vector shape `extract_solution` would have +# produced on the original. +function _nlpf_extract_primal(model::JuMP.AbstractModel, ref_map) + V = JuMP.variable_ref_type(typeof(model)) + T = JuMP.value_type(typeof(model)) + result = Dict{V, Vector{T}}() + for v in collect_all_vars(model) + JuMP.is_fixed(v) && continue + target = ref_map[v] + val = JuMP.value(target) + result[v] = val isa AbstractArray ? vec(val) : [val] + end + return result +end + +# Fix `combination`, run `f()`, then undo. The finite undo is a no-op +# (fixes overwrite in place); the infinite `fix_combination` clears its +# per-support fixed values. + +# TODO: Refactor to avoid this. +function with_fixed_combination( + f, + model::JuMP.AbstractModel, + combination::AbstractDict + ) + undo = fix_combination(model, combination) + try + return f() + finally + undo() + end +end + +# Iter-to-iter NLP warm start: seed the next NLP from the last FEASIBLE +# primal (no-op on the first seed and after an NLPF fall-through, whose +# primal is slack-distorted). Routes through `set_linearization_start`. +function _set_nlp_warm_start(previous) + previous === nothing && return + for (variable, values) in previous.linearization_point + set_linearization_start(variable, values) + end + return +end + +# Finalize the model with the LOA-optimal combination: fix indicators +# at the committed values (in place; the model is already permanently +# relaxed by reformulate_model) and warm-start from the linearization +# point. After this, the model is no longer in a state suitable for +# re-running with a different `gdp_method`. +function commit_combination( + model::JuMP.AbstractModel, + combination::AbstractDict, + linearization_point::AbstractDict + ) + fix_combination(model, combination) + for (variable, values) in linearization_point + set_linearization_start(variable, values) + end + return +end + +""" + fix_combination(model, combination::AbstractDict)::Function + +Fix each indicator in `combination` to its value on `model`, in place, +and return a zero-argument closure that undoes the fixes. + +## Returns +- `Function`: an undo closure called by `with_fixed_combination`. +""" +function fix_combination(model::JuMP.AbstractModel, combination::AbstractDict) + for (indicator, value) in combination + fix_indicator(model, indicator, value) + end + return () -> nothing +end + +""" + set_linearization_start(variable, values::AbstractVector) + +Warm-start `variable` from the linearization point, unwrapping the single +element of `values` as its start value. +""" +set_linearization_start(variable, values::AbstractVector) = + JuMP.set_start_value(variable, only(values)) + +################################################################################ +# COMBO EXTRACTION (master -> NLP) +################################################################################ +# Read each indicator's binary value from the master MILP solution. +# `combination_val` returns a `Bool` (finite indicator) or per-support +# `Vector{Bool}` (infinite); downstream dispatches on that. +function _extract_combination( + model::M, + master::_LOAMaster + ) where {M <: JuMP.AbstractModel} + combination = Dict{LogicalVariableRef{M}, Any}() + for (_, disjunction_data) in _disjunctions(model) + for indicator in disjunction_data.constraint.indicators + combination[indicator] = + combination_val(master.binary_map[indicator]) + end + end + return combination +end + +# Read the master binary and round to `Bool` (a MILP solver returns +# values within its integer tolerance, never exactly 0/1). Branches on +# the value: scalar (finite) or per-support array (infinite). +function combination_val(binary_ref) + val = JuMP.value(binary_ref) + return val isa AbstractArray ? vec(round.(Bool, val)) : round(Bool, val) +end + +################################################################################ +# OA CUT EMISSION +################################################################################ +function add_oa_cuts( + model::JuMP.AbstractModel, + master::_LOAMaster, + result::NamedTuple, + method::LOA + ) + isempty(result.linearization_point) && return + linearization = _linearize_at(master.original_objective, + objective_linearization_point(model, result.linearization_point), + master.objective_ref_map) + _add_objective_cut( + Val(master.objective_sense), master, linearization, method) + add_global_oa_cuts(model, master, result, method) + add_disjunct_oa_cuts(model, master, result, method) + return +end + +""" + objective_linearization_point(model, linearization_point) + +Return the point at which the master objective is linearized. Returns +`linearization_point` unchanged. +""" +objective_linearization_point(::JuMP.AbstractModel, linearization_point) = + linearization_point + +# Slacked objective cut. MIN: `lin <= alpha_oa + sigma`; MAX symmetric. +# The slack (like the disjunct/global cuts) keeps a nonconvex objective +# linearization from making the master infeasible. +function _add_objective_cut( + sense_token::Val, + master::_LOAMaster, + lin, + method::LOA + ) + penalty_sign = _penalty_sign(sense_token) + slack = _penalized_slack(master, method, penalty_sign) + _add_objective_cut_body(sense_token, master, lin, slack) +end +_add_objective_cut_body(::Val{_MOI.MIN_SENSE}, master, lin, slack) = + JuMP.@constraint(master.model, lin <= master.alpha_oa + slack) +_add_objective_cut_body(::Val{_MOI.MAX_SENSE}, master, lin, slack) = + JuMP.@constraint(master.model, lin >= master.alpha_oa - slack) + +""" + add_global_oa_cuts(model, master::_LOAMaster, result::NamedTuple, + method::LOA) + +Add a slacked OA cut to `master` for every nonlinear global (non-disjunct) +constraint of `model`, linearized at `result.linearization_point`. +Variable bounds, linear functions, and reformulation constraints are +skipped. +""" +function add_global_oa_cuts( + model::JuMP.AbstractModel, + master::_LOAMaster, + result::NamedTuple, + method::LOA + ) + penalty_sign = _penalty_sign(Val(master.objective_sense)) + variable_type = JuMP.variable_ref_type(typeof(model)) + reform_set = is_gdp_model(model) ? + Set(_reformulation_constraints(model)) : Set() + for (F, S) in JuMP.list_of_constraint_types(model) + F === variable_type && continue + _is_linear_F(F) && continue + for cref in JuMP.all_constraints(model, F, S) + cref in reform_set && continue + con = JuMP.constraint_object(cref) + con isa JuMP.ScalarConstraint || continue + linearization = _linearize_at(con.func, + result.linearization_point, master.variable_map) + _add_global_oa_row(master, linearization, con.set, + method, penalty_sign) + end + end + return +end + +# Iterate each active disjunct's disjunct constraints, invoking +# `f(binary_ref, active, constraint)` per constraint. Shared by the base +# and InfiniteOpt `add_disjunct_oa_cuts` drivers so neither reimplements +# the active-disjunct walk. +function _each_active_disjunct_constraint(f, model, master, result) + for (indicator, active) in result.combination + is_active(active) || continue + haskey(_indicator_to_constraints(model), indicator) || continue + for cref in _indicator_to_constraints(model)[indicator] + cref isa DisjunctConstraintRef || continue + constraint = _disjunct_constraints(model)[ + JuMP.index(cref)].constraint + f(master.binary_map[indicator], active, constraint) + end + end + return +end + +""" + add_disjunct_oa_cuts(model, master::_LOAMaster, result::NamedTuple, + method::LOA) + +Add a slacked OA cut to `master` for each active disjunct's nonlinear +constraints, linearized at `result.linearization_point` and gated to +match the inner reformulation. Emits one cut per constraint. +""" +function add_disjunct_oa_cuts( + model::JuMP.AbstractModel, + master::_LOAMaster, + result::NamedTuple, + method::LOA + ) + penalty_sign = _penalty_sign(Val(master.objective_sense)) + _each_active_disjunct_constraint(model, master, + result) do binary_ref, _active, constraint + _add_oa_cut_for_constraint(constraint, master, binary_ref, + result.linearization_point, master.variable_map, method, + penalty_sign) + end + return +end + +# Fresh penalized slack added to the master objective. Shared by the +# disjunct, global, and objective cuts so a nonconvex (invalid) +# linearization can't make the master infeasible. +function _penalized_slack(master::_LOAMaster, method::LOA, penalty_sign::Int) + slack = JuMP.@variable(master.model, + lower_bound = 0.0, upper_bound = method.max_slack) + JuMP.set_objective_function(master.model, + JuMP.objective_function(master.model) + + penalty_sign * method.oa_penalty * slack) + return slack +end + +# Linearize the constraint at `linearization_point`, then emit slacked OA +# cut(s) gated to match the inner reformulation (Big-M / MBM or Hull). +# Linear constraints are exact via the reformulation and skipped. +function _add_oa_cut_for_constraint( + constraint::JuMP.AbstractConstraint, + master::_LOAMaster, + binary_ref, + linearization_point::AbstractDict, + var_map::AbstractDict, + method::LOA, + penalty_sign::Int + ) + _is_linear_F(typeof(constraint.func)) && return + linearization = _linearize_at( + constraint.func, linearization_point, var_map) + _emit_disjunct_oa_cut(method.inner_method, + constraint.set, master, binary_ref, linearization, method, + penalty_sign) + return +end + +# Extract RHS from an MOI set. +_set_rhs(s::Union{_MOI.LessThan, _MOI.GreaterThan, _MOI.EqualTo}) = + _MOI.constant(s) +_set_rhs(::Any) = 0.0 + +# The `<= 0` directions of an OA cut for `set`: `lin - rhs` for LessThan, +# `rhs - lin` for GreaterThan, both for EqualTo / Interval. The caller +# adds the gating, slack, and (for Hull) disaggregation per direction. +_oa_cut_terms(set::_MOI.GreaterThan, lin) = (_set_rhs(set) - lin,) +_oa_cut_terms(set::_MOI.EqualTo, lin) = + (lin - _MOI.constant(set), _MOI.constant(set) - lin) +_oa_cut_terms(set::_MOI.Interval, lin) = (lin - set.upper, set.lower - lin) +_oa_cut_terms(set, lin) = (lin - _set_rhs(set),) + +# Big-M / MBM disjunct cut: each direction slacked and gated by +# `M(1 - binary)`. +function _emit_disjunct_oa_cut( + ::Union{BigM, MBM}, + set, + master::_LOAMaster, + binary_ref, + linearization, + method::LOA, + penalty_sign::Int + ) + slack = _penalized_slack(master, method, penalty_sign) + for term in _oa_cut_terms(set, linearization) + JuMP.@constraint(master.model, + term - slack <= method.M_value * (1 - binary_ref)) + end + return +end + +# Convex-hull disjunct cut: `disaggregate_expression` rewrites each +# direction into disaggregated space (variables -> per-disjunct copies, +# constant scaled by the binary), so the cut switches off at `y = 0` with +# no big-M. It is linear in its argument, so the two-sided sets need no +# special casing. +function _emit_disjunct_oa_cut( + ::Hull, + set, + master::_LOAMaster, + binary_ref, + linearization, + method::LOA, + penalty_sign::Int + ) + slack = _penalized_slack(master, method, penalty_sign) + for term in _oa_cut_terms(set, linearization) + body = disaggregate_expression( + master.model, term, binary_ref, master.disaggregator) + JuMP.@constraint(master.model, body - slack <= 0) + end + return +end + +# Slacked global OA row(s): each direction carries a penalized slack so a +# nonconvex linearization can't make the master infeasible. Unknown sets +# fall back to a hard cut. +function _add_global_oa_row( + master::_LOAMaster, + lin, + set::Union{_MOI.LessThan, _MOI.GreaterThan, _MOI.EqualTo, _MOI.Interval}, + method::LOA, + penalty_sign::Int + ) + slack = _penalized_slack(master, method, penalty_sign) + for term in _oa_cut_terms(set, lin) + JuMP.@constraint(master.model, term <= slack) + end + return +end +function _add_global_oa_row(master::_LOAMaster, lin, set, ::LOA, ::Int) + JuMP.@constraint(master.model, lin in set) + return +end + +# Is this indicator active anywhere? Finite: a `Bool`; infinite: a +# per-support `Vector{Bool}`. Skips disjuncts that are off everywhere. +is_active(active::Bool) = active +is_active(active::AbstractVector{Bool}) = any(active) + diff --git a/src/mbm.jl b/src/mbm.jl index 23c5096e..cb726fc1 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -508,7 +508,8 @@ function replace_variables_in_constraint( W = _var_ref_type(T, var_map) new_aff = zero(JuMP.GenericAffExpr{C, W}) for (var, coef) in fun.terms - JuMP.add_to_expression!(new_aff, coef, var_map[var]) + JuMP.add_to_expression!(new_aff, coef, + replace_variables_in_constraint(var, var_map)) end new_aff.constant = new_aff.constant + fun.constant return new_aff @@ -522,7 +523,9 @@ function replace_variables_in_constraint( new_quad = zero(JuMP.GenericQuadExpr{C, W}) for (vars, coef) in fun.terms JuMP.add_to_expression!(new_quad, - coef * var_map[vars.a] * var_map[vars.b]) + coef * + replace_variables_in_constraint(vars.a, var_map) * + replace_variables_in_constraint(vars.b, var_map)) end new_aff = replace_variables_in_constraint(fun.aff, var_map) JuMP.add_to_expression!(new_quad, new_aff) diff --git a/src/utilities.jl b/src/utilities.jl index bca1b748..7e61fdec 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -8,6 +8,23 @@ function _copy_model( return M() end +# OVERRIDABLE. Copy `source`'s decision variables into `target` (with +# bounds and integrality) and return a `source → target` ref map. +# Constraints and objective are NOT copied — the caller adds whichever +# subset it wants. InfiniteOpt overrides to additionally copy +# parameters, derivatives, and parameter functions. +function copy_variables_onto_model( + target::JuMP.AbstractModel, + source::JuMP.AbstractModel + ) + V = JuMP.variable_ref_type(typeof(source)) + ref_map = Dict{V, V}() + for variable in JuMP.all_variables(source) + ref_map[variable] = variable_copy(target, variable) + end + return ref_map +end + """ copy_and_reformulate(model, decision_vars, reform_method, method) @@ -62,6 +79,120 @@ function reformulate_and_relax( return sub, undo_relax end +""" + extract_solution(sub::GDPSubmodel) + +Read the primal solution of `sub.model` after a solve, keyed by the +parent-model decision variables via `sub.fwd_map`. Shape follows +`fwd_map` values: `Vector`-valued fwd_maps (CP/MBM) yield per-support +`Vector`s; scalar fwd_maps (LOA feas) yield scalars. +""" +function extract_solution(sub::GDPSubmodel) + return Dict( + var => JuMP.value.(sub.fwd_map[var]) for var in sub.decision_vars) +end + +################################################################################ +# INDICATOR FIXING +################################################################################ +""" + fix_indicator(model, indicator::LogicalVariableRef, value::Bool) + fix_indicator(binary_ref, value::Bool) + +Fix a logical indicator's binary backing variable to `value` (true → +1.0, false → 0.0). The 3-arg form is the user-facing API: pass the +model and the `LogicalVariableRef`. The 2-arg form takes the binary +backing reference directly (or its complement-form expression +`1 - other_binary` when the indicator was declared as a complement). + +Mirrors [`relax_logical_vars`](@ref) for selectively fixing rather +than relaxing. +""" +fix_indicator(model::JuMP.AbstractModel, + indicator::LogicalVariableRef, value::Bool) = + fix_indicator(_indicator_to_binary(model)[indicator], value) +fix_indicator(binary_ref::JuMP.AbstractVariableRef, value::Bool) = + JuMP.fix(binary_ref, value ? 1.0 : 0.0; force = true) +function fix_indicator( + binary_ref::JuMP.GenericAffExpr, value::Bool + ) + underlying, _ = only(binary_ref.terms) + JuMP.fix(underlying, value ? 0.0 : 1.0; force = true) + return +end + +""" + unfix_indicator(model, indicator::LogicalVariableRef) + unfix_indicator(binary_ref) + +Undo [`fix_indicator`](@ref). No-op if not currently fixed. +""" +unfix_indicator(model::JuMP.AbstractModel, + indicator::LogicalVariableRef) = + unfix_indicator(_indicator_to_binary(model)[indicator]) +function unfix_indicator(binary_ref) + JuMP.is_fixed(binary_ref) && JuMP.unfix(binary_ref) + return +end +function unfix_indicator(binary_ref::JuMP.GenericAffExpr) + underlying = only(keys(binary_ref.terms)) + JuMP.is_fixed(underlying) && JuMP.unfix(underlying) + return +end + +################################################################################ +# NO-GOOD CUTS +################################################################################ +""" + avoid_combination(model, combination) + avoid_combination(model, combination, binary_map) + +Add a no-good cut to `model` excluding `combination` from any future +solution. The added constraint + Σ_{j active} (1 - y_j) + Σ_{j inactive} y_j ≥ 1 +forces at least one indicator to differ from `combination`. + +`combination` maps `LogicalVariableRef` → `Bool` (whether each +indicator was active). The 3-arg form lets you supply an explicit +`binary_map` (defaulting to `_indicator_to_binary(model)`) when the +binaries used in the cut belong to a copy of the original model +(e.g. an LOA master). + +Returns the constraint reference of the added cut. +""" +avoid_combination(model::JuMP.AbstractModel, combination) = + avoid_combination( + model, combination, _indicator_to_binary(model)) +function avoid_combination(model::JuMP.AbstractModel, + combination, binary_map + ) + V = JuMP.variable_ref_type(typeof(model)) + T = JuMP.value_type(typeof(model)) + cut = JuMP.GenericAffExpr{T, V}(zero(T)) + for (indicator, active) in combination + haskey(binary_map, indicator) || continue + add_no_good_terms(cut, binary_map[indicator], active) + end + return JuMP.@constraint(model, cut >= one(T)) +end + +""" + add_no_good_terms(cut, binary_ref, active::Bool) + +Append one indicator's contribution to a no-good cut expression: +adds `1 - y_j` if the indicator was active in the excluded +combination, or `y_j` otherwise. Used by [`avoid_combination`](@ref). +""" +function add_no_good_terms(cut, binary_ref, active::Bool) + if active + JuMP.add_to_expression!(cut, -1.0, binary_ref) + JuMP.add_to_expression!(cut, 1.0) + else + JuMP.add_to_expression!(cut, 1.0, binary_ref) + end + return +end + ################################################################################ # LOGICAL VARIABLE RELAXATION ################################################################################ @@ -399,5 +530,106 @@ function _remap_constraint_to_indicator( ::Dict{DisjunctConstraintRef{M}, DisjunctConstraintRef{M}}, disj_map::Dict{DisjunctionRef{M}, DisjunctionRef{M}} ) where {M <: JuMP.AbstractModel} - return disj_map[con_ref] + return disj_map[con_ref] +end + +################################################################################ +# LINEARIZATION & EXPRESSION CONVERSION +################################################################################ +# First-order Taylor for a single var or affine expression mapped into +# master space. +function _linearize_at( + variable::JuMP.AbstractVariableRef, + ::AbstractDict, + ref_map::AbstractDict + ) + target = ref_map[variable] + return JuMP.GenericAffExpr{Float64, typeof(target)}(0.0, target => 1.0) end +function _linearize_at( + func::JuMP.GenericAffExpr, + ::AbstractDict, + ref_map::AbstractDict + ) + V = valtype(ref_map) + T = JuMP.value_type(V) <: Number ? JuMP.value_type(V) : Float64 + result = JuMP.GenericAffExpr{T, V}(T(func.constant)) + for (variable, coefficient) in func.terms + JuMP.add_to_expression!(result, coefficient, ref_map[variable]) + end + return result +end + +# Convert JuMP expression trees to Julia Expr with +# MOI.VariableIndex leaves for MOI.Nonlinear evaluation. +function _to_nlp_expr(expr::JuMP.GenericNonlinearExpr, idx::Dict) + args = Any[_to_nlp_expr(a, idx) for a in expr.args] + return Expr(:call, expr.head, args...) +end +function _to_nlp_expr(expr::JuMP.GenericAffExpr, idx::Dict) + parts = Any[expr.constant] + for (var, coef) in expr.terms + push!(parts, Expr(:call, :*, coef, _MOI.VariableIndex(idx[var]))) + end + length(parts) == 1 && return parts[1] + return Expr(:call, :+, parts...) +end +function _to_nlp_expr(expr::JuMP.GenericQuadExpr, idx::Dict) + parts = Any[_to_nlp_expr(expr.aff, idx)] + for (pair, coef) in expr.terms + push!(parts, Expr(:call, :*, coef, + _MOI.VariableIndex(idx[pair.a]), + _MOI.VariableIndex(idx[pair.b]))) + end + length(parts) == 1 && return parts[1] + return Expr(:call, :+, parts...) +end +function _to_nlp_expr(var::JuMP.AbstractVariableRef, idx::Dict) + return _MOI.VariableIndex(idx[var]) +end +_to_nlp_expr(x::Number, ::Dict) = x + +# First-order Taylor linearization of a quadratic or nonlinear +# expression at point xk via MOI.Nonlinear reverse-mode AD. +function _linearize_at( + func::Union{JuMP.GenericQuadExpr, JuMP.GenericNonlinearExpr}, + xk::Dict, + ref_map + ) + vars = JuMP.AbstractVariableRef[] + _interrogate_variables(v -> push!(vars, v), func) + unique!(vars) + isempty(vars) && return JuMP.AffExpr(JuMP.value(v -> 0.0, func)) + + n = length(vars) + T = JuMP.value_type(typeof(JuMP.owner_model(vars[1]))) + idx = Dict(vars[i] => i for i in 1:n) + nlp = _MOI.Nonlinear.Model() + _MOI.Nonlinear.set_objective(nlp, _to_nlp_expr(func, idx)) + ord = [_MOI.VariableIndex(i) for i in 1:n] + evaluator = _MOI.Nonlinear.Evaluator( + nlp, _MOI.Nonlinear.SparseReverseMode(), ord) + _MOI.initialize(evaluator, [:Grad]) + + xk_vec = [_unwrap_scalar(get(xk, v, zero(T))) for v in vars] + f_xk = _MOI.eval_objective(evaluator, xk_vec) + grad = zeros(T, n) + _MOI.eval_objective_gradient(evaluator, grad, xk_vec) + + constant = T(f_xk) + for i in 1:n + constant -= grad[i] * xk_vec[i] + end + V = typeof(ref_map[vars[1]]) + result = JuMP.GenericAffExpr{T, V}(constant) + for i in 1:n + iszero(grad[i]) && continue + JuMP.add_to_expression!(result, grad[i], ref_map[vars[i]]) + end + return result +end + +# Unwrap a 1-element `Vector` to its scalar; scalars pass through +# (`extract_solution` returns length-1 vectors for scalar variables). +_unwrap_scalar(v::Real) = v +_unwrap_scalar(v::AbstractVector) = only(v) diff --git a/test/constraints/loa.jl b/test/constraints/loa.jl new file mode 100644 index 00000000..c599102a --- /dev/null +++ b/test/constraints/loa.jl @@ -0,0 +1,503 @@ +using HiGHS, Ipopt, Juniper + +function test_loa_datatype() + method = LOA(HiGHS.Optimizer) + @test method.nlp_optimizer == HiGHS.Optimizer + @test method.mip_optimizer == HiGHS.Optimizer + @test method.max_iter == 10 + @test method.M_value == 1e9 + @test method.max_slack == 1000.0 + @test method.oa_penalty == 1000.0 + @test method.convergence_tol == 1e-6 + @test method.slack_tol == 1e-4 + @test method.inner_method isa BigM + @test method.inner_method.value == 1e9 + + method = LOA(HiGHS.Optimizer; max_iter = 50, M_value = 1e6, + max_slack = 500.0, oa_penalty = 200.0, + convergence_tol = 1e-4, slack_tol = 1e-3) + @test method.max_iter == 50 + @test method.M_value == 1e6 + @test method.max_slack == 500.0 + @test method.oa_penalty == 200.0 + @test method.convergence_tol == 1e-4 + @test method.slack_tol == 1e-3 + @test method.inner_method isa BigM + @test method.inner_method.value == 1e6 + + method = LOA(HiGHS.Optimizer; inner_method = MBM(HiGHS.Optimizer)) + @test method.inner_method isa MBM +end + +function test_set_covering_combos() + model = GDPModel() + @variable(model, x) + @variable(model, Y[1:2], Logical) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x >= 5, Disjunct(Y[2])) + @disjunction(model, Y) + + combos = DP._set_covering_combinations(model) + + # Should cover both Y[1] and Y[2] + all_active = Set() + for combo in combos + for (ind, active) in combo + active && push!(all_active, ind) + end + end + @test length(all_active) == 2 +end + +function test_no_good_cut() + model = GDPModel() + @variable(model, x) + @variable(model, Y[1:2], Logical) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x >= 5, Disjunct(Y[2])) + @disjunction(model, Y) + + DP.reformulate_model(model, BigM(1e9)) + master = DP.build_loa_master( + model, LOA(HiGHS.Optimizer)) + master_model = master.model + + combo = Dict(Y[1] => true, Y[2] => false) + + num_cons_before = length(JuMP.all_constraints( + master_model; + include_variable_in_set_constraints = false)) + DP.avoid_combination(master.model, combo, master.binary_map) + num_cons_after = length(JuMP.all_constraints( + master_model; + include_variable_in_set_constraints = false)) + + @test num_cons_after == num_cons_before + 1 +end + +function test_loa_reformulate_simple() + model = GDPModel(HiGHS.Optimizer) + @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) + + method = LOA(HiGHS.Optimizer) + DP.reformulate_model(model, method) + + @test DP._ready_to_optimize(model) +end + +function test_loa_solve_simple() + # Simple GDP: max x s.t. (x <= 3) OR (x <= 7), 0 <= x <= 10 + # Optimal: x = 7 (select second disjunct) + model = GDPModel(HiGHS.Optimizer) + 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 = LOA(HiGHS.Optimizer)) + @test termination_status(model) == MOI.OPTIMAL + @test objective_value(model) ≈ 7.0 atol=1e-4 +end + +function test_loa_solve_simple_with_mbm() + # Same GDP as test_loa_solve_simple, but inner_method = MBM so the + # NLP model is built with per-constraint tight Ms instead of BigM. + # Optimum unchanged: x = 7 via the second disjunct. + model = GDPModel(HiGHS.Optimizer) + 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 = LOA(HiGHS.Optimizer; + inner_method = MBM(HiGHS.Optimizer))) + @test termination_status(model) == MOI.OPTIMAL + @test objective_value(model) ≈ 7.0 atol=1e-4 +end + +function test_loa_solve_two_disjunctions() + # Two disjunctions: max x + z + # D1: (x <= 3) OR (x <= 7) + # D2: (z <= 2) OR (z <= 5) + # Optimal: x = 7, z = 5 + model = GDPModel(HiGHS.Optimizer) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, 0 <= z <= 10) + @variable(model, Y[1:2], Logical) + @variable(model, W[1:2], Logical) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x <= 7, Disjunct(Y[2])) + @disjunction(model, Y) + @constraint(model, z <= 2, Disjunct(W[1])) + @constraint(model, z <= 5, Disjunct(W[2])) + @disjunction(model, W) + @objective(model, Max, x + z) + + optimize!(model, gdp_method = LOA(HiGHS.Optimizer)) + @test termination_status(model) == MOI.OPTIMAL + @test objective_value(model) ≈ 12.0 atol=1e-4 +end + +function test_loa_error_fallback() + method = LOA(HiGHS.Optimizer) + @test_throws ErrorException DP.reformulate_model(42, method) +end + +function test_loa_nonlinear_global() + # max x s.t. x^2 <= 25 (global), (x <= 3) ∨ (x <= 8), 0 <= x <= 10. + # Disjunct Y[2] permits x up to 8 but the global x^2 <= 25 bounds + # x to 5. Verifies that `add_global_oa_cuts` emits the + # linearization of the global into the master without breaking + # the loop. Result must hit the global-binding optimum. + ipopt = optimizer_with_attributes(Ipopt.Optimizer, + "print_level" => 0, "sb" => "yes") + juniper = optimizer_with_attributes(Juniper.Optimizer, + "nl_solver" => ipopt, "log_levels" => []) + model = GDPModel(juniper) + set_silent(model) + @variable(model, 0 <= x <= 10) + @constraint(model, x^2 <= 25) + @variable(model, Y[1:2], Logical) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x <= 8, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, x) + optimize!(model, + gdp_method = LOA(juniper; mip_optimizer = HiGHS.Optimizer)) + @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + @test objective_value(model) ≈ 5.0 atol = 1e-3 +end + +function test_loa_complement_indicator_nonlinear_disjunct() + # Regression: complement-form indicators store `1 - y_base` (an + # AffExpr) in `binary_map`. When the complement disjunct has a + # nonlinear constraint, `add_disjunct_oa_cuts` feeds that AffExpr + # straight into the OA cut gating term `M(1 - binary)`, which must + # accept an AffExpr binary (not just a plain variable ref). + # + # Setup: 0 <= x <= 10, Y2 ≡ ¬Y1. + # Disjunct(Y1): x <= 3 (linear) + # Disjunct(Y2): x^2 <= 64 (nonlinear, optimum: x = 8 with Y2) + ipopt = optimizer_with_attributes(Ipopt.Optimizer, + "print_level" => 0, "sb" => "yes") + juniper = optimizer_with_attributes(Juniper.Optimizer, + "nl_solver" => ipopt, "log_levels" => []) + model = GDPModel(juniper) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, Y1, Logical) + @variable(model, Y2, Logical, logical_complement = Y1) + @constraint(model, x <= 3, Disjunct(Y1)) + @constraint(model, x^2 <= 64, Disjunct(Y2)) + @disjunction(model, [Y1, Y2]) + @objective(model, Max, x) + optimize!(model, + gdp_method = LOA(juniper; mip_optimizer = HiGHS.Optimizer)) + @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + @test objective_value(model) ≈ 8.0 atol = 1e-3 +end + +function test_loa_nlpf_infeasible_disjunct() + # Y1 disjunct constraint x^2 >= 200 is NLP-infeasible against the + # variable bound x in [0, 10] (max x^2 = 100). The primary NLP at + # the Y1=true seed therefore fails and NLPF (V&G 1990 eq. 8) must + # kick in: slack `u` on the GreaterThan, minimize u, return the + # x = 10 / u = 100 primal as the linearization site for the OA cut + # on x^2 >= 200. With the master's per-cut penalized slack + # absorbing the resulting (invalid in isolation) linearization and + # the no-good cut forbidding Y1, LOA should still converge to the + # Y2=true optimum x = 5. + ipopt = optimizer_with_attributes(Ipopt.Optimizer, + "print_level" => 0, "sb" => "yes") + juniper = optimizer_with_attributes(Juniper.Optimizer, + "nl_solver" => ipopt, "log_levels" => []) + model = GDPModel(juniper) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, Y[1:2], Logical) + @constraint(model, x^2 >= 200, Disjunct(Y[1])) + @constraint(model, x <= 5, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, x) + optimize!(model, + gdp_method = LOA(juniper; mip_optimizer = HiGHS.Optimizer)) + @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + @test objective_value(model) ≈ 5.0 atol = 1e-3 + @test value(Y[2]) ≈ 1.0 atol = 1e-6 + @test value(Y[1]) ≈ 0.0 atol = 1e-6 +end + +function test_loa_sense_primitives() + # Sense-dispatched primitives the main loop reads through. The + # finite/infinite solve tests cover every branch except the + # minimize-sense gap (no minimizing model reaches a master bound in + # those tests), so pin all six here for both senses directly. + minv = Val(MOI.MIN_SENSE) + maxv = Val(MOI.MAX_SENSE) + @test DP._penalty_sign(minv) == 1 + @test DP._penalty_sign(maxv) == -1 + @test DP._worst_objective(minv) == Inf + @test DP._worst_objective(maxv) == -Inf + @test DP._is_better(minv, 3.0, 5.0) + @test !DP._is_better(minv, 5.0, 3.0) + @test DP._is_better(maxv, 5.0, 3.0) + @test !DP._is_better(maxv, 3.0, 5.0) + @test DP._gap(minv, 5.0, 3.0) == 2.0 + @test DP._gap(maxv, 3.0, 5.0) == 2.0 +end + +function test_loa_iteration_loop() + # Force the master/NLP refinement loop to run its body rather than + # exhaust the master on the seeds or converge on the first master + # solve. Two disjunctions over disjoint x/z half-lines put the + # optimum on an OFF-diagonal combination the set-covering seeds (the + # diagonal (Y1,W1) and (Y2,W2)) never try. After seeding, the master + # stays feasible and its bound beats the seed incumbent, so LOA + # enters the loop body, extracts the off-diagonal combination, solves + # its NLP, and updates the incumbent before converging. `Min` sense + # also drives the minimize-sense gap path. + # min x - z + # D1: (x <= 4) [Y1] v (x >= 6) [Y2] + # D2: (z <= 4) [W1] v (z >= 6) [W2] + # Optimum: Y1 & W2 -> x = 0, z = 10, obj = -10 (off-diagonal). + model = GDPModel(HiGHS.Optimizer) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, 0 <= z <= 10) + @variable(model, Y[1:2], Logical) + @variable(model, W[1:2], Logical) + @constraint(model, x <= 4, Disjunct(Y[1])) + @constraint(model, x >= 6, Disjunct(Y[2])) + @disjunction(model, Y) + @constraint(model, z <= 4, Disjunct(W[1])) + @constraint(model, z >= 6, Disjunct(W[2])) + @disjunction(model, W) + @objective(model, Min, x - z) + optimize!(model, gdp_method = LOA(HiGHS.Optimizer)) + @test termination_status(model) == MOI.OPTIMAL + @test objective_value(model) ≈ -10.0 atol = 1e-4 + @test value(Y[1]) ≈ 1.0 atol = 1e-6 + @test value(W[2]) ≈ 1.0 atol = 1e-6 +end + +function test_loa_time_limits() + # Exercise the wall-clock budget paths. A finite `time_limit` + # (overall cap) drives `_cap_remaining_time` during the subproblem + # solves and the final-solve budget reset; a finite + # `iteration_time_limit` (loop-only budget) drives the post-loop + # time-limit restore from no prior limit (the `::Nothing` path). + # Both limits are generous: they govern the path taken, not the + # result. + function build_model() + model = GDPModel(HiGHS.Optimizer) + 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) + return model + end + + model = build_model() + optimize!(model, gdp_method = LOA(HiGHS.Optimizer; time_limit = 60.0)) + @test termination_status(model) == MOI.OPTIMAL + @test objective_value(model) ≈ 7.0 atol = 1e-4 + + model = build_model() + optimize!(model, gdp_method = LOA(HiGHS.Optimizer; + iteration_time_limit = 60.0, time_limit = Inf)) + @test termination_status(model) == MOI.OPTIMAL + @test objective_value(model) ≈ 7.0 atol = 1e-4 +end + +function test_loa_no_feasible_incumbent() + # Every disjunct combination is NLP-infeasible: the global x == 20 + # contradicts the x <= 10 bound under any selection. NLPF does not + # slack equalities, so NLPF also fails (no values) and `_solve_nlp` + # returns its infeasible fallback. No seed or loop NLP yields a + # feasible incumbent, so `best_result` stays `nothing`, the master + # is infeasible (no bound), and LOA exits on the "no feasible + # incumbent" warning from `_report_loa_gap`. + model = GDPModel(HiGHS.Optimizer) + set_silent(model) + @variable(model, 0 <= x <= 10) + @constraint(model, x == 20) + @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) + method = LOA(HiGHS.Optimizer) + @test_logs (:warn, r"no feasible incumbent") match_mode = :any begin + DP.reformulate_model(model, method) + end + @test DP._ready_to_optimize(model) +end + +function test_loa_hull_linear() + # Same GDP as test_loa_solve_simple, but inner_method = Hull so the + # NLP and master carry disaggregated variables and the disjunct + # cuts are convex-hull cuts rather than big-M. Optimum unchanged: + # x = 7 via the second disjunct. + model = GDPModel(HiGHS.Optimizer) + 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 = LOA(HiGHS.Optimizer; + inner_method = Hull())) + @test termination_status(model) == MOI.OPTIMAL + @test objective_value(model) ≈ 7.0 atol = 1e-4 +end + +function test_loa_hull_two_disjunctions() + # Two independent disjunctions under Hull; optimum x = 7, z = 5. + model = GDPModel(HiGHS.Optimizer) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, 0 <= z <= 10) + @variable(model, Y[1:2], Logical) + @variable(model, W[1:2], Logical) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x <= 7, Disjunct(Y[2])) + @disjunction(model, Y) + @constraint(model, z <= 2, Disjunct(W[1])) + @constraint(model, z <= 5, Disjunct(W[2])) + @disjunction(model, W) + @objective(model, Max, x + z) + optimize!(model, gdp_method = LOA(HiGHS.Optimizer; + inner_method = Hull())) + @test termination_status(model) == MOI.OPTIMAL + @test objective_value(model) ≈ 12.0 atol = 1e-4 +end + +function test_loa_hull_nonlinear_global() + # max x s.t. x^2 <= 25 (global), (x <= 3) ∨ (x <= 8) under Hull. + # Linear disjuncts enter the master as exact Hull perspectives; the + # global nonlinear enters via global OA cuts. Optimum: x = 5. + ipopt = optimizer_with_attributes(Ipopt.Optimizer, + "print_level" => 0, "sb" => "yes") + juniper = optimizer_with_attributes(Juniper.Optimizer, + "nl_solver" => ipopt, "log_levels" => []) + model = GDPModel(juniper) + set_silent(model) + @variable(model, 0 <= x <= 10) + @constraint(model, x^2 <= 25) + @variable(model, Y[1:2], Logical) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x <= 8, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, x) + optimize!(model, gdp_method = LOA(juniper; + mip_optimizer = HiGHS.Optimizer, inner_method = Hull())) + @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + @test objective_value(model) ≈ 5.0 atol = 1e-3 +end + +function test_loa_hull_nonlinear_disjunct() + # Nonlinear disjunct constraint under Hull: this is the case big-M + # and Hull genuinely differ on. 0 <= x <= 10, Y1: x <= 3, + # Y2: x^2 <= 64. The convex-hull OA cut disaggregates the + # linearization of x^2. Optimum: x = 8 via Y2. + ipopt = optimizer_with_attributes(Ipopt.Optimizer, + "print_level" => 0, "sb" => "yes") + juniper = optimizer_with_attributes(Juniper.Optimizer, + "nl_solver" => ipopt, "log_levels" => []) + model = GDPModel(juniper) + 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 = LOA(juniper; + mip_optimizer = HiGHS.Optimizer, inner_method = Hull())) + @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + @test objective_value(model) ≈ 8.0 atol = 1e-3 +end + +function test_loa_hull_complement_nonlinear() + # Hull with a complement indicator Y2 ≡ ¬Y1 and a nonlinear + # disjunct under Y2. The disaggregator is keyed by the complement + # binary `1 - y1` (an AffExpr); the cut must disaggregate against + # it. Optimum: x = 8 via Y2. + ipopt = optimizer_with_attributes(Ipopt.Optimizer, + "print_level" => 0, "sb" => "yes") + juniper = optimizer_with_attributes(Juniper.Optimizer, + "nl_solver" => ipopt, "log_levels" => []) + model = GDPModel(juniper) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, Y1, Logical) + @variable(model, Y2, Logical, logical_complement = Y1) + @constraint(model, x <= 3, Disjunct(Y1)) + @constraint(model, x^2 <= 64, Disjunct(Y2)) + @disjunction(model, [Y1, Y2]) + @objective(model, Max, x) + optimize!(model, gdp_method = LOA(juniper; + mip_optimizer = HiGHS.Optimizer, inner_method = Hull())) + @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + @test objective_value(model) ≈ 8.0 atol = 1e-3 +end + +function test_loa_hull_nested_sink() + # The inner disjunction (over y) is a disjunct of the outer (gated by + # z[1]). Hull's nested redispatch must propagate the LOA sink so the + # inner disaggregations are recorded, else LOA-Hull emits the nested + # cut on the aggregated variable. + model = GDPModel() + @variable(model, 0 <= x <= 10) + @variable(model, y[1:2], Logical) + @variable(model, z[1:2], Logical) + @constraint(model, x <= 3, Disjunct(y[1])) + @constraint(model, x <= 8, Disjunct(y[2])) + @disjunction(model, y, Disjunct(z[1])) + @constraint(model, x <= 1, Disjunct(z[2])) + @disjunction(model, z) + sink = Dict{Any, Any}() + DP.reformulate_model(model, Hull(; sink = sink)) + @test haskey(sink, (x, y[1])) + @test haskey(sink, (x, y[2])) +end + +@testset "LOA" begin + test_loa_datatype() + test_set_covering_combos() + test_no_good_cut() + test_loa_reformulate_simple() + test_loa_solve_simple() + test_loa_solve_simple_with_mbm() + test_loa_solve_two_disjunctions() + test_loa_error_fallback() + test_loa_nonlinear_global() + test_loa_complement_indicator_nonlinear_disjunct() + test_loa_nlpf_infeasible_disjunct() + test_loa_sense_primitives() + test_loa_iteration_loop() + test_loa_time_limits() + test_loa_no_feasible_incumbent() + test_loa_hull_linear() + test_loa_hull_two_disjunctions() + test_loa_hull_nonlinear_global() + test_loa_hull_nonlinear_disjunct() + test_loa_hull_complement_nonlinear() + test_loa_hull_nested_sink() +end diff --git a/test/constraints/mbm.jl b/test/constraints/mbm.jl index 11f164da..ab9e8287 100644 --- a/test/constraints/mbm.jl +++ b/test/constraints/mbm.jl @@ -52,7 +52,7 @@ function test__replace_variables_quad_numeric_map() @test result3.aff.terms[y] ≈ 3.0 end -function test_replace_variables_in_constraint() +function testreplace_variables_in_constraint() model = Model() sub_model = Model() @variable(model, x[1:3]) @@ -809,7 +809,7 @@ end test_mbm() test__var_ref_type_numeric_map() test__replace_variables_quad_numeric_map() - test_replace_variables_in_constraint() + testreplace_variables_in_constraint() test_prepare_max_M_objective() test_raw_M() test_maximize_M() diff --git a/test/extensions/InfiniteDisjunctiveProgramming.jl b/test/extensions/InfiniteDisjunctiveProgramming.jl index 47bb884b..c34b27e6 100644 --- a/test/extensions/InfiniteDisjunctiveProgramming.jl +++ b/test/extensions/InfiniteDisjunctiveProgramming.jl @@ -513,6 +513,76 @@ function test_mbm_finite_and_integer_var() [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] end +function test_mbm_infinite_simple() + model = InfiniteGDPModel(HiGHS.Optimizer) + set_silent(model) + model = InfiniteGDPModel(HiGHS.Optimizer) + set_silent(model) + K = 4 + @infinite_parameter(model, t ∈ [0, 1], num_supports = K) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y[1:2], InfiniteLogical(t)) + @constraint(model, x >= 5, Disjunct(Y[1])) + @constraint(model, x <= 3, Disjunct(Y[2])) + @disjunction(model, Y) + JuMP.fix(Y[2], true) # force disj 2 active + @objective(model, Min, ∫(x, t)) + DP.reformulate_model(model, BigM(10.0)) + set_optimizer(model, HiGHS.Optimizer) + set_silent(model) + optimize!(model, ignore_optimize_hook = true) + sol = DP.extract_solution(model) + @test haskey(sol, x) + @test length(sol[x]) == K + @test all(v -> isapprox(v, 0.0; atol=1e-6), sol[x]) +end + +# add_cut adds one pointwise-sum cut to the transformation backend and +# marks the backend ready so the next optimize! does NOT re-transcribe. +function test_add_cut_infinite() + model = InfiniteGDPModel(HiGHS.Optimizer) + set_silent(model) + K = 3 + @infinite_parameter(model, t ∈ [0, 1], num_supports = K) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y[1:2], InfiniteLogical(t)) + @constraint(model, x >= 5, Disjunct(Y[1])) + @constraint(model, x <= 3, Disjunct(Y[2])) + @disjunction(model, Y) + DP.reformulate_model(model, BigM(10.0)) + InfiniteOpt.build_transformation_backend!(model) + transcribed = InfiniteOpt.transformation_model(model) + n_before = JuMP.num_constraints(transcribed; + count_variable_in_set_constraints = false) + rBM_sol = Dict(x => [1.0, 2.0, 3.0]) + sep_sol = Dict(x => [0.5, 1.5, 2.5]) + DP.add_cut(model, [x], rBM_sol, sep_sol) + n_after = JuMP.num_constraints(transcribed; + count_variable_in_set_constraints = false) + @test n_after == n_before + 1 + # set_transformation_backend_ready(true) — next optimize! should + # reuse without re-transcribing (otherwise our cut would be lost) + @test InfiniteOpt.transformation_backend_ready(model) +end + +# MBM with finite + integer variables in an InfiniteModel. +function test_mbm_finite_and_integer_var() + model = InfiniteGDPModel(HiGHS.Optimizer) + set_silent(model) + @infinite_parameter(model, t ∈ [0, 1], num_supports = 10) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, 0 <= w <= 5, Int) + @variable(model, Y[1:2], InfiniteLogical(t)) + @constraint(model, x + w >= 5, Disjunct(Y[1])) + @constraint(model, x + w <= 3, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Min, ∫(x, t) + w) + @test optimize!(model, + gdp_method = MBM(HiGHS.Optimizer)) isa Nothing + @test termination_status(model) in + [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] +end + function test_mbm_infinite_simple() model = InfiniteGDPModel(HiGHS.Optimizer) set_silent(model) @@ -791,6 +861,324 @@ function test_CuttingPlanes_multiparameter() [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] end +function test_loa_infinite_nonlinear_global() + # max ∫x dt s.t. x(t)^2 <= 25 (global, per-support nonlinear), + # (x <= 3) ∨ (x <= 8), 0 <= x <= 10 over t ∈ [0, 1]. + # Disjunct Y[2] permits x up to 8 but the global x^2 <= 25 caps + # x at 5. The per-support global transcribes to an `AbstractArray` + # of scalar constraints, so this exercises the array branch of + # `add_global_oa_cuts`. Without the global cut the + # master would allow x = 8 and report 8.0; the binding optimum + # is ∫5 dt = 5. + ipopt = optimizer_with_attributes(Ipopt.Optimizer, + "print_level" => 0, "sb" => "yes") + model = InfiniteGDPModel(ipopt) + set_silent(model) + @infinite_parameter(model, t ∈ [0, 1], num_supports = 10) + @variable(model, 0 <= x <= 10, Infinite(t)) + @constraint(model, x^2 <= 25) + @variable(model, Y[1:2], InfiniteLogical(t)) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x <= 8, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, ∫(x, t)) + optimize!(model, + gdp_method = LOA(ipopt; mip_optimizer = HiGHS.Optimizer)) + @test termination_status(model) in + (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + @test objective_value(model) ≈ 5.0 atol = 1e-2 +end + +function test_loa_infinite_complement_indicator() + # Regression: `_indicator_to_binary(model)[Y2]` returns the AffExpr + # `1 - binary(Y1)` for a logical-complement indicator. Earlier the + # extension's `fix_combination` called `JuMP.fix` directly on that + # AffExpr and crashed; now it delegates to base `fix_indicator` + # which unwraps and inverts. + # + # Linear disjuncts only — a nonlinear disjunct constraint on an + # infinite variable would also hit a separate gap where + # `_infinite_cut_info` for a finite indicator does not slice the + # per-support linearization point. Out of scope for this regression. + # + # max ∫ x dt with 0 ≤ x ≤ 10 over t ∈ [0,1]: + # Y1: x ≤ 3 Y2 ≡ ¬Y1: x ≤ 8 + # Y2 is the maximizer (x = 8 across t), giving ∫8 dt = 8. + ipopt = optimizer_with_attributes(Ipopt.Optimizer, + "print_level" => 0, "sb" => "yes") + model = InfiniteGDPModel(ipopt) + set_silent(model) + @infinite_parameter(model, t ∈ [0, 1], num_supports = 5) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y1, Logical) + @variable(model, Y2, Logical, logical_complement = Y1) + @constraint(model, x <= 3, Disjunct(Y1)) + @constraint(model, x <= 8, Disjunct(Y2)) + @disjunction(model, [Y1, Y2]) + @objective(model, Max, ∫(x, t)) + optimize!(model, + gdp_method = LOA(ipopt; mip_optimizer = HiGHS.Optimizer)) + @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + @test objective_value(model) ≈ 8.0 atol = 1e-2 +end + +function test_loa_infinite_complement_nonlinear_disjunct() + # Regression: a finite (or complement-form AffExpr) indicator + # gating a NONLINEAR constraint on an infinite variable. Earlier + # `_infinite_cut_info` decided fan-out from the indicator only — finite + # indicator → single un-sliced site → `_linearize_at` received + # per-support `Vector{Float64}` values and tripped on + # `_unwrap_scalar`'s `only(v)`. Fixed by inspecting the constraint + # expression too: any infinite var found there governs the fan-out + # supports, even when the indicator is scalar. + # + # max ∫ x dt with 0 ≤ x ≤ 10 over t ∈ [0,1]: + # Y1: x ≤ 3 Y2 ≡ ¬Y1: x² ≤ 64 + # Y2 (complement) is the maximizer (x = 8 across t), giving 8. + ipopt = optimizer_with_attributes(Ipopt.Optimizer, + "print_level" => 0, "sb" => "yes") + juniper = optimizer_with_attributes(Juniper.Optimizer, + "nl_solver" => ipopt, "log_levels" => []) + model = InfiniteGDPModel(juniper) + set_silent(model) + @infinite_parameter(model, t ∈ [0, 1], num_supports = 5) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y1, Logical) + @variable(model, Y2, Logical, logical_complement = Y1) + @constraint(model, x <= 3, Disjunct(Y1)) + @constraint(model, x^2 <= 64, Disjunct(Y2)) + @disjunction(model, [Y1, Y2]) + @objective(model, Max, ∫(x, t)) + optimize!(model, + gdp_method = LOA(juniper; mip_optimizer = HiGHS.Optimizer)) + @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + @test objective_value(model) ≈ 8.0 atol = 1e-2 +end + +function test_loa_infinite_multidim_parameter() + # Regression: multi-D dependent parameter group. `supports(ξ)` is + # a 2 × N matrix, not a vector. Earlier the extension called + # `vec(supports(ξ))`, flattened to 2·N scalars, then called + # `binary_ref(scalar)` on a 2-D infinite var — dim mismatch. + # Fixed by returning `eachcol(supports(ξ))` and splatting vector + # supports through `_at_support`. + # + # max w with 0 ≤ x ≤ 10 over ξ[1:2] ∈ [0,1]², w ≤ x at every joint + # support: + # Y[1]: x ≤ 3 Y[2]: x ≤ 8 + # Y[2] is the maximizer (x = 8 across ξ allows w = 8). + ipopt = optimizer_with_attributes(Ipopt.Optimizer, + "print_level" => 0, "sb" => "yes") + model = InfiniteGDPModel(ipopt) + set_silent(model) + @infinite_parameter(model, ξ[1:2] ∈ [0, 1], num_supports = 4) + @variable(model, 0 <= x <= 10, Infinite(ξ)) + @variable(model, 0 <= w <= 10) + @constraint(model, w <= x) + @variable(model, Y[1:2], InfiniteLogical(ξ)) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x <= 8, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, w) + optimize!(model, + gdp_method = LOA(ipopt; mip_optimizer = HiGHS.Optimizer)) + @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + @test objective_value(model) ≈ 8.0 atol = 1e-2 +end + +function test_loa_infinite_nlpf_infeasible_disjunct() + # InfiniteOpt LOA: Y[1]'s per-support constraint x >= 200 is + # NLP-infeasible against the bound x in [0, 10]. With infinite + # indicators `Y[1:2], InfiniteLogical(t)`, the infeasible + # Y[1]-everywhere seed makes the primary NLP fail, so NLPF runs: it + # copies the model, re-pins the binaries on the copy via + # `nlpf_fix_on_copy`, slacks the inequality, and returns a + # linearization point. The disjunct is deliberately LINEAR: + # the model is then convex, so LOA's master bound is rigorous and + # it converges to Y[2] active everywhere (x = 5), giving ∫x dt = 5. + # (A nonconvex infeasible disjunct such as x^2 >= 200 would make LOA + # a heuristic and let the local NLP solver report a constraint- + # violating point as solved — out of scope for the NLPF path here.) + ipopt = optimizer_with_attributes(Ipopt.Optimizer, + "print_level" => 0, "sb" => "yes") + juniper = optimizer_with_attributes(Juniper.Optimizer, + "nl_solver" => ipopt, "log_levels" => []) + model = InfiniteGDPModel(juniper) + set_silent(model) + @infinite_parameter(model, t ∈ [0, 1], num_supports = 3) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y[1:2], InfiniteLogical(t)) + @constraint(model, x >= 200, Disjunct(Y[1])) + @constraint(model, x <= 5, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, ∫(x, t)) + optimize!(model, + gdp_method = LOA(juniper; mip_optimizer = HiGHS.Optimizer)) + @test termination_status(model) in + (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + @test objective_value(model) ≈ 5.0 atol = 1e-2 +end + +function test_loa_infinite_aggregate_global() + # min ∫x dt s.t. ∫(x^2, t) <= 4 (aggregate global), x >= y, + # (y >= 1) ∨ (y >= 3), 0 <= x, y <= 10 over t ∈ [0, 1]. + # The aggregate global transcribes to a single scalar (the + # measure is flattened), exercising the non-array branch of + # `add_global_oa_cuts`. Y[1] (y >= 1) is the cheaper + # disjunct: x = 1 satisfies x >= y and ∫x^2 = 1 <= 4, giving + # objective ∫1 dt = 1. + ipopt = optimizer_with_attributes(Ipopt.Optimizer, + "print_level" => 0, "sb" => "yes") + model = InfiniteGDPModel(ipopt) + set_silent(model) + @infinite_parameter(model, t ∈ [0, 1], num_supports = 10) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, 0 <= y <= 10) + @constraint(model, ∫(x^2, t) <= 4) + @constraint(model, x >= y) + @variable(model, Y[1:2], Logical) + @constraint(model, y >= 1, Disjunct(Y[1])) + @constraint(model, y >= 3, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Min, ∫(x, t)) + optimize!(model, + gdp_method = LOA(ipopt; mip_optimizer = HiGHS.Optimizer)) + @test termination_status(model) in + (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + @test objective_value(model) ≈ 1.0 atol = 1e-2 +end + +function test_loa_infinite_iteration_loop() + # Force the InfiniteModel LOA main loop to fix per-support + # Vector{Bool} combinations, exercising the in-place fix-constraint + # stash: create on the first main-loop fix, update via + # set_normalized_rhs on the committed re-fix. Two disjunctions over + # disjoint x/z half-lines put the optimum on an off-diagonal + # combination the scalar set-covering seeds never try, so the master + # stays feasible and the loop body runs with per-support values. + # min ∫(x - z) dt + # D1: (x <= 4) [Y1] ∨ (x >= 6) [Y2] + # D2: (z <= 4) [W1] ∨ (z >= 6) [W2] + # Optimum: Y1 & W2 everywhere → x = 0, z = 10, ∫(x - z) = -10. + model = InfiniteGDPModel(HiGHS.Optimizer) + set_silent(model) + @infinite_parameter(model, t ∈ [0, 1], num_supports = 3) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, 0 <= z <= 10, Infinite(t)) + @variable(model, Y[1:2], InfiniteLogical(t)) + @variable(model, W[1:2], InfiniteLogical(t)) + @constraint(model, x <= 4, Disjunct(Y[1])) + @constraint(model, x >= 6, Disjunct(Y[2])) + @disjunction(model, Y) + @constraint(model, z <= 4, Disjunct(W[1])) + @constraint(model, z >= 6, Disjunct(W[2])) + @disjunction(model, W) + @objective(model, Min, ∫(x - z, t)) + optimize!(model, gdp_method = LOA(HiGHS.Optimizer)) + @test termination_status(model) in + (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + @test objective_value(model) ≈ -10.0 atol = 1e-3 +end + +function test_loa_infinite_hull_linear() + # InfiniteOpt LOA with inner_method = Hull. Linear disjuncts enter + # the master as exact per-support Hull perspectives over + # disaggregated infinite variables. max ∫x dt, Y1: x <= 3, + # Y2: x <= 8. Optimum: x = 8 across t → ∫8 dt = 8. + model = InfiniteGDPModel(HiGHS.Optimizer) + set_silent(model) + @infinite_parameter(model, t ∈ [0, 1], num_supports = 5) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y[1:2], InfiniteLogical(t)) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x <= 8, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, ∫(x, t)) + optimize!(model, gdp_method = LOA(HiGHS.Optimizer; + inner_method = Hull())) + @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + @test objective_value(model) ≈ 8.0 atol = 1e-2 +end + +function test_loa_infinite_hull_nonlinear_disjunct() + # The infinite Hull case big-M differs from: a NONLINEAR constraint + # on an infinite variable inside a disjunct, under a plain infinite + # indicator. The convex-hull OA cut is built per support over the + # disaggregated infinite variable (point-evaluated). max ∫x dt, + # Y1: x <= 3, Y2: x^2 <= 64. Optimum: x = 8 via Y2. + ipopt = optimizer_with_attributes(Ipopt.Optimizer, + "print_level" => 0, "sb" => "yes") + juniper = optimizer_with_attributes(Juniper.Optimizer, + "nl_solver" => ipopt, "log_levels" => []) + model = InfiniteGDPModel(juniper) + set_silent(model) + @infinite_parameter(model, t ∈ [0, 1], num_supports = 5) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y[1:2], InfiniteLogical(t)) + @constraint(model, x <= 3, Disjunct(Y[1])) + @constraint(model, x^2 <= 64, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, ∫(x, t)) + optimize!(model, gdp_method = LOA(juniper; + mip_optimizer = HiGHS.Optimizer, inner_method = Hull())) + @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + @test objective_value(model) ≈ 8.0 atol = 1e-2 +end + +function test_loa_infinite_hull_complement_nonlinear() + # Hull with a complement indicator Y2 ≡ ¬Y1 gating a nonlinear + # constraint on an infinite variable. The per-support disaggregator + # is keyed by the point-evaluated complement binary `1 - y1(t_k)` + # (an AffExpr). max ∫x dt, Y1: x <= 3, Y2: x^2 <= 64. Optimum 8. + ipopt = optimizer_with_attributes(Ipopt.Optimizer, + "print_level" => 0, "sb" => "yes") + juniper = optimizer_with_attributes(Juniper.Optimizer, + "nl_solver" => ipopt, "log_levels" => []) + model = InfiniteGDPModel(juniper) + set_silent(model) + @infinite_parameter(model, t ∈ [0, 1], num_supports = 5) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y1, Logical) + @variable(model, Y2, Logical, logical_complement = Y1) + @constraint(model, x <= 3, Disjunct(Y1)) + @constraint(model, x^2 <= 64, Disjunct(Y2)) + @disjunction(model, [Y1, Y2]) + @objective(model, Max, ∫(x, t)) + optimize!(model, gdp_method = LOA(juniper; + mip_optimizer = HiGHS.Optimizer, inner_method = Hull())) + @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + @test objective_value(model) ≈ 8.0 atol = 1e-2 +end + +function test_loa_infinite_hull_finite_var() + # A FINITE variable disaggregated inside a nonlinear disjunct + # constraint gated by an INFINITE indicator. The fan-out slices the + # binary per support, so the disaggregator must key the finite + # variable's copy by the per-support binary `y(t_k)` (driven off the + # binary, not the variable). max ∫x dt + w over t ∈ [0,1]: + # Y1: x^2 + w^2 <= 25 Y2: x <= 0 + # Y1 is active; x(t) = w = sqrt(12.5), so the optimum is + # ∫sqrt(12.5) dt + sqrt(12.5) = 2 sqrt(12.5) ≈ 7.0711. + ipopt = optimizer_with_attributes(Ipopt.Optimizer, + "print_level" => 0, "sb" => "yes") + juniper = optimizer_with_attributes(Juniper.Optimizer, + "nl_solver" => ipopt, "log_levels" => []) + model = InfiniteGDPModel(juniper) + set_silent(model) + @infinite_parameter(model, t ∈ [0, 1], num_supports = 4) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, 0 <= w <= 10) + @variable(model, Y[1:2], InfiniteLogical(t)) + @constraint(model, x^2 + w^2 <= 25, Disjunct(Y[1])) + @constraint(model, x <= 0, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, ∫(x, t) + w) + optimize!(model, gdp_method = LOA(juniper; + mip_optimizer = HiGHS.Optimizer, inner_method = Hull())) + @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + @test objective_value(model) ≈ 7.0711 atol = 1e-2 +end + @testset "InfiniteDisjunctiveProgramming" begin @testset "Model" begin @@ -855,4 +1243,29 @@ end test_CuttingPlanes_multiparameter() end + # LOA on InfiniteOpt needs `JuMP.copy_model(::InfiniteModel)` + # (build_loa_master). That lives in the InfiniteOpt fork and is + # not in the registered release, so guard the suite: it runs + # (and validates the global-cut paths) when the fork is dev'd + # into the env, and skips cleanly otherwise. + @testset "LOA" begin + if hasmethod(JuMP.copy_model, Tuple{InfiniteModel}) + test_loa_infinite_nonlinear_global() + test_loa_infinite_aggregate_global() + test_loa_infinite_complement_indicator() + test_loa_infinite_complement_nonlinear_disjunct() + test_loa_infinite_multidim_parameter() + test_loa_infinite_nlpf_infeasible_disjunct() + test_loa_infinite_iteration_loop() + test_loa_infinite_hull_linear() + test_loa_infinite_hull_nonlinear_disjunct() + test_loa_infinite_hull_complement_nonlinear() + test_loa_infinite_hull_finite_var() + else + @info "Skipping InfiniteOpt LOA tests: " * + "JuMP.copy_model(::InfiniteModel) unavailable " * + "(registered InfiniteOpt; dev the fork to run)." + end + end + end diff --git a/test/runtests.jl b/test/runtests.jl index 06e8813d..ca34f816 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,6 +18,7 @@ include("constraints/mbm.jl") include("constraints/bigm.jl") include("constraints/psplit.jl") include("constraints/cuttingplanes.jl") +include("constraints/loa.jl") include("constraints/hull.jl") include("constraints/fallback.jl") include("constraints/disjunction.jl") diff --git a/test/solve.jl b/test/solve.jl index c3ca9441..c53cf04a 100644 --- a/test/solve.jl +++ b/test/solve.jl @@ -62,7 +62,15 @@ function test_linear_gdp_example(m, use_complements = false) @test value(Y[2]) @test !value(W[1]) @test !value(W[2]) - + + @test optimize!(m, gdp_method = LOA(HiGHS.Optimizer)) isa Nothing + @test termination_status(m) == MOI.OPTIMAL + @test objective_value(m) ≈ 11 atol=1e-3 + @test value.(x) ≈ [9,2] atol=1e-3 + @test !value(Y[1]) + @test value(Y[2]) + @test !value(W[1]) + @test !value(W[2]) m_copy, ref_map = JuMP.copy_model(m) lv_map = DP.copy_gdp_data(m, m_copy, ref_map) @@ -130,11 +138,20 @@ function test_quadratic_gdp_example(use_complements = false) #psplit does not wo @test optimize!(m, gdp_method = PSplit(2,m)) isa Nothing @test termination_status(m) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] - @test objective_value(m) ≈ 6.1237 atol=1e-3 - @test value.(x) ≈ [4.0825, 2.0412] atol=1e-3 - @test !value(Y[1]) + @test objective_value(m) ≈ 6.1237 atol=1e-3 + @test value.(x) ≈ [4.0825, 2.0412] atol=1e-3 + @test !value(Y[1]) @test value(Y[2]) - @test !value(W[1]) + @test !value(W[1]) + @test !value(W[2]) + + @test optimize!(m, gdp_method = LOA(optimizer)) isa Nothing + @test termination_status(m) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] + @test objective_value(m) ≈ 6.1237 atol=1e-3 + @test value.(x) ≈ [4.0825, 2.0412] atol=1e-3 + @test !value(Y[1]) + @test value(Y[2]) + @test !value(W[1]) @test !value(W[2]) end diff --git a/test/utilities.jl b/test/utilities.jl index 2b87e6e1..2a98366d 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -238,6 +238,77 @@ function test_get_constant_variable() @test DP.get_constant(x) == 0.0 end +function test_linearize_nonlinear_exp() + # exp(x) + y at (1, 2): + # f = e + 2, ∇f = [e, 1] + # linear: e*(x-1) + 1*(y-2) + (e+2) = e*x + y + model = GDPModel() + @variable(model, x) + @variable(model, y) + func = @expression(model, exp(x) + y) + xk = Dict{JuMP.AbstractVariableRef, Float64}( + x => 1.0, y => 2.0) + id_map = Dict(x => x, y => y) + lin = DP._linearize_at(func, xk, id_map) + @test JuMP.constant(lin) ≈ 0.0 atol = 1e-8 + @test JuMP.coefficient(lin, x) ≈ exp(1.0) atol = 1e-8 + @test JuMP.coefficient(lin, y) ≈ 1.0 atol = 1e-8 +end + +function test_linearize_nonlinear_sin() + # sin(x) at x = π/6: + # f = 0.5, f' = cos(π/6) = √3/2 + # linear: 0.5 + (√3/2)(x - π/6) + model = GDPModel() + @variable(model, x) + func = @expression(model, sin(x)) + xk = Dict{JuMP.AbstractVariableRef, Float64}( + x => π / 6) + id_map = Dict(x => x) + lin = DP._linearize_at(func, xk, id_map) + expected_const = 0.5 - (√3 / 2) * (π / 6) + @test JuMP.constant(lin) ≈ expected_const atol = 1e-8 + @test JuMP.coefficient(lin, x) ≈ √3 / 2 atol = 1e-8 +end + +function test_linearize_nonlinear_multivar() + # exp(x) * sin(y) at (1, π/2): + # f = e*1 = e, ∂f/∂x = e*sin(π/2) = e, ∂f/∂y = e*cos(π/2) = 0 + # linear: e + e*(x-1) + 0*(y-π/2) = e*x + model = GDPModel() + @variable(model, x) + @variable(model, y) + func = @expression(model, exp(x) * sin(y)) + xk = Dict{JuMP.AbstractVariableRef, Float64}( + x => 1.0, y => π / 2) + id_map = Dict(x => x, y => y) + lin = DP._linearize_at(func, xk, id_map) + @test JuMP.constant(lin) ≈ 0.0 atol = 1e-8 + @test JuMP.coefficient(lin, x) ≈ exp(1.0) atol = 1e-8 + @test JuMP.coefficient(lin, y) ≈ 0.0 atol = 1e-8 +end + +function test_to_nlp_expr() + model = GDPModel() + @variable(model, x) + @variable(model, y) + idx = Dict(x => 1, y => 2) + + # NonlinearExpr + nl = @expression(model, exp(x)) + e = DP._to_nlp_expr(nl, idx) + @test e == Expr(:call, :exp, MOI.VariableIndex(1)) + + # AffExpr + aff = @expression(model, 2x + 3y + 1) + e = DP._to_nlp_expr(aff, idx) + @test e isa Expr + @test e.head == :call && e.args[1] == :+ + + # Number + @test DP._to_nlp_expr(42, idx) == 42 +end + @testset "Utility Functions" begin test_all_variables() test_collect_all_vars() @@ -245,4 +316,8 @@ end test_get_constant_quadratic() test_get_constant_number() test_get_constant_variable() + test_linearize_nonlinear_exp() + test_linearize_nonlinear_sin() + test_linearize_nonlinear_multivar() + test_to_nlp_expr() end