_is_linear test coverage

Author: dnguyen227

ext/InfiniteDisjunctiveProgramming.jl

@@ -408,25 +408,6 @@ function _transcribed_binary_refs(model::InfiniteOpt.InfiniteModel, indicator)
     return [1.0 - binary for binary in binaries]
 end
 
-# Transcribed combination for an InfiniteModel: every support of the
-# underlying transcribed binary takes the indicator's value.
-function _transcribed_combination(
-    model::InfiniteOpt.InfiniteModel,
-    combination
-    )
-    result = Dict{JuMP.VariableRef, Bool}()
-    binary_map = DP._indicator_to_binary(model)
-    for (indicator, active) in combination
-        binary_ref = binary_map[indicator]
-        value = DP._underlying_value(binary_ref, active)
-        for binary in _transcribed_refs(InfiniteOpt.transformation_variable(
-            DP._underlying_binary(binary_ref)))
-            result[binary] = value
-        end
-    end
-    return result
-end
-
 # Transcribe the inner Hull disaggregation map per support. A finite
 # variable under an infinite indicator keys its single disaggregated
 # copy by each per-support binary reference, mirroring the per-support
@@ -468,8 +449,6 @@ function DP.build_loa_problem(
                 DP._underlying_binary(binary_ref))))
     end
     unique!(binaries)
-    combinations = [_transcribed_combination(model, combination)
-        for combination in DP._set_covering_combinations(model)]
 
     disjunct_constraints = Tuple{Any, Any, Any}[]
     for (_, disjunction) in DP._disjunctions(model)
@@ -509,7 +488,7 @@ function DP.build_loa_problem(
         end
     end
 
-    return DP._LOAProblem(nlp, binaries, combinations,
+    return DP._LOAProblem(nlp, binaries,
         disjunct_constraints, global_constraints,
         _transcribed_disaggregation_map(model, disaggregation_map))
 end

src/datatypes.jl

@@ -587,6 +587,8 @@ iteration) and a master MILP accumulating OA and no-good cuts.
 - `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.
+- `set_cover_max_iter::Int`: max iterations of the set-covering seed
+  loop.
 - `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.
@@ -602,6 +604,7 @@ struct LOA{O, P, R, T} <: AbstractReformulationMethod
     mip_optimizer::P
     inner_method::R
     max_iter::Int
+    set_cover_max_iter::Int
     M_value::T
     max_slack::T
     oa_penalty::T
@@ -613,6 +616,7 @@ struct LOA{O, P, R, T} <: AbstractReformulationMethod
         nlp_optimizer::O;
         mip_optimizer::P = nlp_optimizer,
         max_iter::Int = 10,
+        set_cover_max_iter::Int = 8,
         M_value::T = 1e9,
         max_slack::T = 1e3,
         oa_penalty::T = 1e3,
@@ -626,24 +630,24 @@ struct LOA{O, P, R, T} <: AbstractReformulationMethod
             "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,
+            max_iter, set_cover_max_iter, M_value, max_slack, oa_penalty,
             convergence_tol, slack_tol,
             iteration_time_limit, time_limit)
     end
 end
 
 # The problem the LOA loop operates on: the model solved as the NLP
 # subproblem, the binary variables backing the indicators, the
-# set-covering seed combinations, the nonlinear disjunct
-# `(binary_ref, function, set)` triples (`binary_ref` is the binary or
-# its `1 - y` complement expression), the nonlinear global
-# `(function, set)` pairs, and the `(variable, binary_ref) ->
+# nonlinear disjunct `(binary_ref, function, set)` triples (`binary_ref`
+# is the binary or its `1 - y` complement expression), the nonlinear
+# global `(function, set)` pairs, and the `(variable, binary_ref) ->
 # disaggregated variable` map from an inner Hull reformulation
-# (`nothing` for Big-M / MBM). Built by `build_loa_problem`.
+# (`nothing` for Big-M / MBM). Built by `build_loa_problem`. The
+# set-covering seed generates its combinations on the fly from the
+# master, so none are stored here.
 struct _LOAProblem{M <: JuMP.AbstractModel, V <: JuMP.AbstractVariableRef}
     nlp::M
     binaries::Vector{V}
-    covering_combinations::Vector{Dict{V, Bool}}
     disjunct_constraints::Vector{Tuple{Any, Any, Any}}
     global_constraints::Vector{Tuple{Any, Any}}
     disaggregation_map::Any

src/loa.jl

@@ -44,11 +44,33 @@ function reformulate_model(model::JuMP.AbstractModel, method::LOA)
     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.
+    # Set-covering seed. This mimics Pyomo GDPopt's set-covering
+    # initialization: borrow the master as a covering MILP (only its
+    # objective changes), let it pick a combination that activates the
+    # nonlinear disjuncts still lacking a linearization, solve the NLP
+    # there, and seed the resulting OA and no-good cuts. Each NLP
+    # warm-starts from the last feasible primal.
     previous_result = nothing
-    for combination in problem.covering_combinations
+    cover_disjuncts = _cover_disjuncts(problem)
+    needs_cover = trues(length(cover_disjuncts))
+    num_covered = 0
+    for iteration in 1:method.set_cover_max_iter
+        (iteration == 1 || any(needs_cover)) || break
         time() < loop_deadline || break
+        # Swap in the covering objective, solve, read off the combination,
+        # then restore the OA objective (with its accumulated slack
+        # penalties) before emitting cuts against it.
+        oa_objective = JuMP.objective_function(master.model)
+        JuMP.@objective(master.model, Max,
+            _cover_objective(master, cover_disjuncts, needs_cover,
+                num_covered))
+        _cap_remaining_time(master.model, loop_deadline)
+        JuMP.optimize!(master.model)
+        solved = JuMP.is_solved_and_feasible(master.model)
+        combination = solved ? _extract_combination(problem, master) : nothing
+        JuMP.set_objective_sense(master.model, master.objective_sense)
+        JuMP.set_objective_function(master.model, oa_objective)
+        solved || break
         _set_nlp_warm_start(previous_result)
         result = _solve_nlp(problem, combination, method;
             deadline = loop_deadline)
@@ -60,6 +82,17 @@ function reformulate_model(model::JuMP.AbstractModel, method::LOA)
             best_result = result
         end
         result.feasible && (previous_result = result)
+        # A disjunct counts as covered only once it is active in a
+        # feasible NLP. An infeasible combination still contributes its
+        # no-good cut above but leaves the coverage targets untouched.
+        if result.feasible
+            for i in eachindex(cover_disjuncts)
+                needs_cover[i] || continue
+                _disjunct_active(result.combination, cover_disjuncts[i]) &&
+                    (needs_cover[i] = false)
+            end
+            num_covered = count(!, needs_cover)
+        end
     end
 
     # Master/NLP loop: `alpha_oa` is the bound, the NLP refines the
@@ -143,25 +176,50 @@ function reformulate_model(::M, ::LOA) where {M}
 end
 
 ################################################################################
-#                       SET-COVERING INITIALIZATION (simple version from pyomo)
+#                       SET-COVERING INITIALIZATION
 ################################################################################
-# `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
-    ]
+# The nonlinear disjuncts to cover: one entry per distinct indicator that
+# owns a nonlinear disjunct constraint, carrying that indicator's
+# `binary_ref`. Keyed by `(underlying binary, active value)` so the two
+# disjuncts of a single-binary disjunction (`y` and its `1 - y`
+# complement) stay distinct. Purely linear disjuncts need no cover: the
+# inner reformulation already places them in the master exactly.
+function _cover_disjuncts(problem::_LOAProblem)
+    V = eltype(problem.binaries)
+    seen = Set{Tuple{V, Bool}}()
+    disjuncts = Any[]
+    for (binary_ref, _, _) in problem.disjunct_constraints
+        key = (_underlying_binary(binary_ref),
+            _underlying_value(binary_ref, true))
+        key in seen && continue
+        push!(seen, key)
+        push!(disjuncts, binary_ref)
+    end
+    return disjuncts
+end
+
+# The set-covering objective (master space), mimicking Pyomo GDPopt:
+# maximize active disjuncts weighted `num_covered + 1` if still uncovered
+# else `1`, so one uncovered disjunct outweighs every covered one and each
+# solve must activate a new disjunct when the logic allows. Empty
+# `disjuncts` gives the zero expression: a constant objective that still
+# seeds one feasible combination (needed to bound `alpha_oa`).
+function _cover_objective(
+    master::_LOAMaster,
+    disjuncts,
+    needs_cover,
+    num_covered::Int
+    )
+    T = JuMP.value_type(typeof(master.model))
+    V = JuMP.variable_ref_type(typeof(master.model))
+    expr = JuMP.GenericAffExpr{T, V}(zero(T))
+    for i in eachindex(disjuncts)
+        weight = needs_cover[i] ? num_covered + 1 : 1
+        activation = _remap_indicator_to_binary(disjuncts[i],
+            master.variable_map)
+        JuMP.add_to_expression!(expr, T(weight), activation)
+    end
+    return expr
 end
 
 ################################################################################
@@ -176,20 +234,6 @@ _underlying_binary(binary_ref::JuMP.GenericAffExpr) =
 _underlying_value(::JuMP.AbstractVariableRef, active::Bool) = active
 _underlying_value(::JuMP.GenericAffExpr, active::Bool) = !active
 
-# Translate an indicator-level combination onto the underlying binary
-# variables, inverting the value for complement-form indicators.
-function _binary_combination(model::JuMP.AbstractModel, combination)
-    V = JuMP.variable_ref_type(typeof(model))
-    binary_map = _indicator_to_binary(model)
-    result = Dict{V, Bool}()
-    for (indicator, active) in combination
-        binary_ref = binary_map[indicator]
-        result[_underlying_binary(binary_ref)] =
-            _underlying_value(binary_ref, active)
-    end
-    return result
-end
-
 """
     build_loa_problem(
         model::JuMP.AbstractModel,
@@ -199,11 +243,11 @@ end
 
 Build the problem the LOA loop operates on from `model` (already
 reformulated by the LOA inner method): the NLP subproblem, the binary
-variables backing the indicators, the set-covering seed combinations,
-the nonlinear disjunct `(binary_ref, function, set)` triples, the
-nonlinear global `(function, set)` pairs, and the Hull disaggregation
-map. The NLP is `model` itself; the InfiniteOpt extension overloads
-this to build the problem from the transcribed backend instead.
+variables backing the indicators, the nonlinear disjunct
+`(binary_ref, function, set)` triples, the nonlinear global
+`(function, set)` pairs, and the Hull disaggregation map. The NLP is
+`model` itself; the InfiniteOpt extension overloads this to build the
+problem from the transcribed backend instead.
 
 ## Returns
 - `_LOAProblem`: the problem.
@@ -219,8 +263,6 @@ function build_loa_problem(
     binaries = V[_underlying_binary(binary_ref)
         for (_, binary_ref) in binary_map]
     unique!(binaries)
-    combinations = [_binary_combination(model, combination)
-        for combination in _set_covering_combinations(model)]
 
     disjunct_constraints = Tuple{Any, Any, Any}[]
     for (_, disjunction) in _disjunctions(model)
@@ -255,7 +297,7 @@ function build_loa_problem(
     binary_disaggregations = disaggregation_map === nothing ? nothing :
         Dict((variable, binary_map[indicator]) => disaggregated
             for ((variable, indicator), disaggregated) in disaggregation_map)
-    return _LOAProblem(model, binaries, combinations, disjunct_constraints,
+    return _LOAProblem(model, binaries, disjunct_constraints,
         global_constraints, binary_disaggregations)
 end
 

test/constraints/loa.jl

@@ -5,6 +5,7 @@ function test_loa_datatype()
     @test method.nlp_optimizer == HiGHS.Optimizer
     @test method.mip_optimizer == HiGHS.Optimizer
     @test method.max_iter == 10
+    @test method.set_cover_max_iter == 8
     @test method.M_value == 1e9
     @test method.max_slack == 1000.0
     @test method.oa_penalty == 1000.0
@@ -29,22 +30,31 @@ function test_loa_datatype()
     @test method.inner_method isa MBM
 end
 
-function test_set_covering_combos()
+function test_cover_disjuncts()
+    # Only disjuncts owning a nonlinear constraint need covering: the two
+    # `x^2` disjuncts do, the two linear `x` disjuncts do not.
     model = GDPModel()
-    @variable(model, x)
+    @variable(model, -5 <= x <= 5)
     @variable(model, Y[1:2], Logical)
-    @constraint(model, x <= 3, Disjunct(Y[1]))
-    @constraint(model, x >= 5, Disjunct(Y[2]))
+    @constraint(model, x^2 <= 3, Disjunct(Y[1]))
+    @constraint(model, x^2 >= 5, Disjunct(Y[2]))
     @disjunction(model, Y)
+    @variable(model, W[1:2], Logical)
+    @constraint(model, x <= 1, Disjunct(W[1]))
+    @constraint(model, x >= -1, Disjunct(W[2]))
+    @disjunction(model, W)
 
-    combos = DP._set_covering_combinations(model)
+    DP.reformulate_model(model, BigM(1e9))
+    method = LOA(HiGHS.Optimizer)
+    problem = DP.build_loa_problem(model, method)
+    disjuncts = DP._cover_disjuncts(problem)
 
-    # K = 2 combinations, each activating exactly one indicator, and
-    # together covering both.
-    @test length(combos) == 2
-    @test all(count(values(combo)) == 1 for combo in combos)
-    @test combos[1][Y[1]] && !combos[1][Y[2]]
-    @test combos[2][Y[2]] && !combos[2][Y[1]]
+    # Two nonlinear disjuncts to cover; the linear W disjuncts omitted.
+    @test length(disjuncts) == 2
+    binary_map = DP._indicator_to_binary(model)
+    underlying = Set(DP._underlying_binary(dref) for dref in disjuncts)
+    @test underlying == Set([DP._underlying_binary(binary_map[Y[1]]),
+        DP._underlying_binary(binary_map[Y[2]])])
 end
 
 function test_oa_cut_terms()
@@ -57,6 +67,17 @@ function test_oa_cut_terms()
     @test DP._oa_cut_terms(MOI.Interval(1.0, 4.0), 5.0) == (1.0, -4.0)
 end
 
+function test_is_linear_F()
+    # Scalar and vector variable-ref / affine function types are linear
+    # (copied into the master at build time); anything else, e.g. a
+    # quadratic, is nonlinear and enters only as an OA cut.
+    @test DP._is_linear_F(JuMP.VariableRef)
+    @test DP._is_linear_F(JuMP.AffExpr)
+    @test DP._is_linear_F(Vector{JuMP.VariableRef})
+    @test DP._is_linear_F(Vector{JuMP.AffExpr})
+    @test !DP._is_linear_F(JuMP.QuadExpr)
+end
+
 function test_no_good_cut()
     model = GDPModel()
     @variable(model, x)
@@ -71,7 +92,12 @@ function test_no_good_cut()
     master = DP._build_loa_master(problem, method)
     master_model = master.model
 
-    combo = DP._binary_combination(model, Dict(Y[1] => true, Y[2] => false))
+    binary_map = DP._indicator_to_binary(model)
+    combo = Dict(
+        DP._underlying_binary(binary_map[Y[1]]) =>
+            DP._underlying_value(binary_map[Y[1]], true),
+        DP._underlying_binary(binary_map[Y[2]]) =>
+            DP._underlying_value(binary_map[Y[2]], false))
 
     num_cons_before = length(JuMP.all_constraints(
         master_model;
@@ -84,7 +110,6 @@ function test_no_good_cut()
     @test num_cons_after == num_cons_before + 1
     # The cut is (1 - y1) + y2 >= 1, i.e. normalized -y1 + y2 >= 0:
     # excludes exactly the (Y1 active, Y2 inactive) combination.
-    binary_map = DP._indicator_to_binary(model)
     y1 = master.variable_map[binary_map[Y[1]]]
     y2 = master.variable_map[binary_map[Y[2]]]
     @test JuMP.normalized_coefficient(cref, y1) == -1.0
@@ -645,8 +670,9 @@ end
 
 @testset "LOA" begin
     test_loa_datatype()
-    test_set_covering_combos()
+    test_cover_disjuncts()
     test_oa_cut_terms()
+    test_is_linear_F()
     test_no_good_cut()
     test_loa_reformulate_simple()
     test_loa_solve_simple()