Merge pull request #219 from lanl-ansi/use_start_as_local

Add option to use starting value as local solution for presolve (with changes)

CHANGELOG.md

@@ -4,6 +4,8 @@
 - Re-org and clean-up in partition additions within bounding MIP model
 - Updated to cheaper `AffExpr` evaluation within linking constraints
 - Minor update in gap evaluataion after OBBT termination with fixed iterations
+- Added solver option `use_start_as_incumbent` to use warm start value as also an intial incumbent solution  in presolve without invoking a local solver. (#218)
+- Added unit tests for testing warm start point
 
 ## v0.5.0
 - New feature: Linking constraints for multilinear terms with uniform and adaptive partitions (significant speed up in run times for multilinear problems: http://www.optimization-online.org/DB_HTML/2022/07/8974.html) (@jongeunkim)
@@ -49,7 +51,7 @@
 - Bug fix in `relax_integrality` under presolve
 - Added Rosenbrock function
 - Added user-definable `presolve_bt_improv_tol` for variable range reductions in OBBT
-- Minor update on finite values for obj-bound in `create_bound_tightening_model`
+- Minor update on finite values for obj-bound in `create_obbt_model`
 - NLP solution rounding issue resolved (now compatible with default tol, 1e-6), to avoid the variable solution outside discretization (#190)
 - `circleN` instance updated to `circle_MINLPLib` to eliminate Pavito's mixed-integer cycling in unit tests (#190)
 - Clean-up in printing variable solutions (`variable_values`)

docs/src/functions.md

@@ -23,11 +23,12 @@ min_vertex_cover
 
 ## Presolve Methods
 ```@docs
-bound_tightening
+bound_tightening_wrapper
 optimization_based_bound_tightening
-create_bound_tightening_model
-solve_bound_tightening_model
+create_obbt_model
+solve_obbt_model
 resolve_var_bounds
+post_objective_bound
 ```
 
 ## Utility Methods

docs/src/parameters.md

@@ -12,35 +12,37 @@ Here are a few general solver options which control the performance of Alpine:
 
 * `time_limit (default = Inf)`: total limit on run time for Alpine in seconds.
 
-* `max_iter (default = 999)`: total number of iterations allowed in the [`global_solve`](@ref). 
+* `max_iter (default = 999)`: total number of iterations allowed in the [`global_solve`](@ref).
 
 * `rel_gap (default = 1e-4)`: relative gap considered for global convergence during [`global_solve`](@ref). Bounds are evaluated using ``\frac{UB-LB}{UB} \cdot 100 \%``.
 
-* `tol (default = 1e-6)`: numerical tolerance used during the process of [`global_solve`](@ref). 
+* `tol (default = 1e-6)`: numerical tolerance used during the process of [`global_solve`](@ref).
 
 ## Adaptive Partitioning Options
 * `apply_partitioning (default = true)`: applies Alpine's built-in MIP-based partitioning algorithm only when activated; else terminates with the presolve solution.
 
-* `partition_scaling_factor (default = 10)`: used during [`add_adaptive_partition`](@ref) for scaling the width of new partitions relative to the active partition chosen in the sequentially-solved lower-bounding MIP models. This value can substantially affect the run time for global convergence; this value can be set to different integer values (>= 4) for various classes of problems. 
+* `partition_scaling_factor (default = 10)`: used during [`add_adaptive_partition`](@ref) for scaling the width of new partitions relative to the active partition chosen in the sequentially-solved lower-bounding MIP models. This value can substantially affect the run time for global convergence; this value can be set to different integer values (>= 4) for various classes of problems.
 
 * `disc_var_pick (default = 0)`: controls Alpine's algorithm used for selecting variables for partitioning; `0` is for max-cover, `1` is for minimum-vertex-cover. This parameter allows functional inputs.
 
 * `disc_add_partition_method`: allows functional input on how new partitions could be constructed.
 
 ## Presolve Options
 
-* `presolve_track_time (default = true)`: includes/excludes presolve run time in the total Alpine's run time. 
+* `presolve_track_time (default = true)`: includes/excludes presolve run time in the total Alpine's run time.
 
-* `presolve_bt (default = false)`: performs sequential, optimization-based bound tightening (OBBT) at the presolve step. 
+* `presolve_bt (default = false)`: performs sequential, optimization-based bound tightening (OBBT) at the presolve step.
 
-* `presolve_bt_max_iter (default = 9999)`: maximum number of iterations allowed using the sequential OBBT step. 
+* `presolve_bt_max_iter (default = 9999)`: maximum number of iterations allowed using the sequential OBBT step.
 
-* `presolve_bt_width_tol (default = 1e-3)`: numerical tolerance value used in the OBBT step. Note that smaller values of this tolerance can lead to inaccurate solutions.  
+* `presolve_bt_width_tol (default = 1e-3)`: numerical tolerance value used in the OBBT step. Note that smaller values of this tolerance can lead to inaccurate solutions.
 
 * `presolve_bt_algo (default = 1)`: method chosen to perform the presolve step; choose `1` for the built-in OBBT, else `2` for user-input functions for presolve.
 
-* `presolve_bt_relax_integrality (default = false)`: relaxes the integrality of the existing integer variables in the OBBT step, if the input problem is an MINLP. 
+* `presolve_bt_relax_integrality (default = false)`: relaxes the integrality of the existing integer variables in the OBBT step, if the input problem is an MINLP.
 
-* `presolve_bt_mip_time_limit (default = Inf)`: time limit for individual MILPs solved during the sequential OBBT procedure. 
+* `presolve_bt_mip_time_limit (default = Inf)`: time limit for individual MILPs solved during the sequential OBBT procedure.
 
-Note that the above-mentioned list of solver options is not comprehensive, but can be found in [solver.jl](https://github.com/lanl-ansi/Alpine.jl/blob/master/src/solver.jl). 
+* `use_start_as_incumbent (default = false)`: if `true`, Alpine does not perform any local optimization during the presolve and uses the starting value instead (*warning*: Alpine assumes the feasibility of the starting value and does not check it).
+
+Note that the above-mentioned list of solver options is not comprehensive, but can be found in [solver.jl](https://github.com/lanl-ansi/Alpine.jl/blob/master/src/solver.jl).

examples/MINLPs/blend.jl

@@ -616,7 +616,7 @@ function blend029_gl(; solver = nothing)
         0.0,
     ]
     for i in 1:102
-        set_start_value(x[i], warmstarter[i])
+        JuMP.set_start_value(x[i], warmstarter[i])
     end
 
     @constraint(m, x[1] + x[4] + x[7] + x[10] + x[49] == 1)  #= e2: =#

src/MOI_wrapper/MOI_wrapper.jl

@@ -86,7 +86,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
     presolve_infeasible::Bool                            # Presolve infeasibility detection flag
     best_bound::Float64                                  # Best bound from MIP
     best_obj::Float64                                    # Best feasible objective value
-    initial_warmval::Vector{Float64}                     # Warmstart values set to Alpine
+    warm_start_value::Vector{Float64}                    # Warmstart values set to Alpine
     best_sol::Vector{Float64}                            # Best feasible solution
     best_bound_sol::Vector{Float64}                      # Best bound solution (arg-min)
     presolve_best_rel_gap::Float64                       # Post-OBBT relative optimality gap = |best_obj - best_bound|/|best_obj|*100
@@ -100,7 +100,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
 
     # Logging information and status
     logs::Dict{Symbol,Any}                          # Logging information
-    detected_feasible_solution::Bool
+    detected_incumbent::Bool
     detected_bound::Bool
     status::Dict{Symbol,MOI.TerminationStatusCode}   # Detailed status of every iteration in the algorithm
     alpine_status::MOI.TerminationStatusCode         # Current Alpine's status
@@ -226,7 +226,7 @@ function MOI.empty!(m::Optimizer)
         Vector{Vector{Float64}}(undef, m.options.disc_consecutive_forbid)
 
     m.best_obj = Inf
-    m.initial_warmval = Float64[]
+    m.warm_start_value = Float64[]
     m.best_sol = Float64[]
     m.best_bound = -Inf
     m.presolve_best_rel_gap = Inf
@@ -264,7 +264,7 @@ function MOI.add_variable(model::Optimizer)
     push!(model.l_var_orig, -Inf)
     push!(model.u_var_orig, Inf)
     push!(model.var_type_orig, :Cont)
-    push!(model.initial_warmval, 0.0)
+    push!(model.warm_start_value, 0.0)
     push!(model.best_sol, 0.0)
     return MOI.VariableIndex(model.num_var_orig)
 end
@@ -278,7 +278,7 @@ function MOI.set(
     vi::MOI.VariableIndex,
     value::Union{Real,Nothing},
 )
-    model.best_sol[vi.value] = model.initial_warmval[vi.value] = something(value, 0.0)
+    model.best_sol[vi.value] = model.warm_start_value[vi.value] = something(value, 0.0)
     return
 end
 

src/bounding_model.jl

@@ -9,17 +9,18 @@ will choose one specific partition as the lower bound solution. The more partiti
 better or finer bounding model relax the original MINLP while the more efforts required to solve 
 this MILP is required.
 """
-function create_bounding_mip(m::Optimizer; use_disc = nothing)  
+function create_bounding_mip(m::Optimizer; use_disc = nothing)
     if (use_disc === nothing)
-        if (m.logs[:n_iter] == 1) && (m.status[:local_solve] in STATUS_OPT || m.status[:local_solve] in STATUS_LIMIT)
+        if (m.logs[:n_iter] == 1) &&
+           (m.status[:local_solve] in union(STATUS_OPT, STATUS_LIMIT, STATUS_WARM_START))
             # Setting up an initial partition
-            Alp.add_partition(m, use_solution = m.best_sol)  
+            Alp.add_partition(m, use_solution = m.best_sol)
         elseif m.logs[:n_iter] >= 2
             # Add subsequent iteration partitions
             Alp.add_partition(m)
         end
         discretization = m.discretization
-    else 
+    else
         discretization = use_disc
     end
 
@@ -59,11 +60,11 @@ function amp_post_vars(m::Optimizer; kwargs...)
     if haskey(options, :use_disc)
         l_var = [
             options[:use_disc][i][1] for
-            i in 1:(m.num_var_orig + m.num_var_linear_mip + m.num_var_nonlinear_mip)
+            i in 1:(m.num_var_orig+m.num_var_linear_mip+m.num_var_nonlinear_mip)
         ]
         u_var = [
             options[:use_disc][i][end] for
-            i in 1:(m.num_var_orig + m.num_var_linear_mip + m.num_var_nonlinear_mip)
+            i in 1:(m.num_var_orig+m.num_var_linear_mip+m.num_var_nonlinear_mip)
         ]
     else
         l_var = m.l_var_tight
@@ -72,10 +73,10 @@ function amp_post_vars(m::Optimizer; kwargs...)
 
     JuMP.@variable(
         m.model_mip,
-        x[i = 1:(m.num_var_orig + m.num_var_linear_mip + m.num_var_nonlinear_mip)]
+        x[i = 1:(m.num_var_orig+m.num_var_linear_mip+m.num_var_nonlinear_mip)]
     )
 
-    for i in 1:(m.num_var_orig + m.num_var_linear_mip + m.num_var_nonlinear_mip)
+    for i in 1:(m.num_var_orig+m.num_var_linear_mip+m.num_var_nonlinear_mip)
         # Interestingly, not enforcing category of lifted variables is able to improve performance
         if i <= m.num_var_orig
             if m.var_type_orig[i] == :Bin

src/heuristics.jl

@@ -7,7 +7,7 @@ function update_disc_cont_var(m::Optimizer)
     length(m.candidate_disc_vars) <= 15 && return   # Algorithm Separation Point
 
     # If no feasible solution is found, do NOT update
-    if !m.detected_feasible_solution
+    if !m.detected_incumbent
         println("no feasible solution detected. No update disc var selection.")
         return
     end

src/log.jl

@@ -218,7 +218,7 @@ function create_status!(m)
 
     status[:local_solve] = MOI.OPTIMIZE_NOT_CALLED # Status of local solve
     status[:bounding_solve] = MOI.OPTIMIZE_NOT_CALLED # Status of bounding solve
-    m.detected_feasible_solution = false
+    m.detected_incumbent = false
     m.detected_bound = false
 
     return m.status = status
@@ -240,14 +240,14 @@ function summary_status(m::Optimizer)
     #               happens when lower bound problem is extremely hard to solve
     # :Unknown : termination with no exception recorded
 
-    if m.detected_bound && m.detected_feasible_solution
+    if m.detected_bound && m.detected_incumbent
         m.alpine_status =
             m.best_rel_gap > Alp.get_option(m, :rel_gap) ? MOI.OTHER_LIMIT : MOI.OPTIMAL
     elseif m.status[:bounding_solve] == MOI.INFEASIBLE
         m.alpine_status = MOI.INFEASIBLE
-    elseif m.detected_bound && !m.detected_feasible_solution
+    elseif m.detected_bound && !m.detected_incumbent
         m.alpine_status = MOI.OTHER_LIMIT
-    elseif !m.detected_bound && m.detected_feasible_solution
+    elseif !m.detected_bound && m.detected_incumbent
         m.alpine_status = MOI.LOCALLY_SOLVED
     else
         @warn "  [EXCEPTION] Indefinite Alpine status. Please report your instance (& solver configuration) as an issue (https://github.com/lanl-ansi/Alpine.jl/issues) to help us make Alpine better."

src/main_algorithm.jl

@@ -10,6 +10,9 @@ const STATUS_LIMIT = [
 ]
 const STATUS_OPT =
     [MOI.OPTIMAL, MOI.LOCALLY_SOLVED, MOI.ALMOST_OPTIMAL, MOI.ALMOST_LOCALLY_SOLVED]
+
+const STATUS_WARM_START = [MOI.OPTIMIZE_NOT_CALLED]
+
 const STATUS_INF = [MOI.INFEASIBLE, MOI.LOCALLY_INFEASIBLE]
 
 function features_available(m::Optimizer)
@@ -42,7 +45,7 @@ function load!(m::Optimizer)
     elseif m.objective_function isa Nothing
         m.obj_expr_orig = Expr(:call, :+)
     else
-        m.obj_expr_orig = _moi_function_to_expr(m.objective_function)
+        m.obj_expr_orig = Alp._moi_function_to_expr(m.objective_function)
     end
 
     # Collect original variable type and build dynamic variable type space
@@ -159,8 +162,6 @@ function MOI.optimize!(m::Optimizer)
         println(
             "====================================================================================================",
         )
-    else
-        println("  Presolve terminated with a global optimal solution")
     end
 
     Alp.summary_status(m)
@@ -203,14 +204,40 @@ function presolve(m::Optimizer)
     start_presolve = time()
     Alp.get_option(m, :log_level) > 0 && printstyled("PRESOLVE \n", color = :cyan)
     Alp.get_option(m, :log_level) > 0 && println("  Doing local search")
-    Alp.local_solve(m, presolve = true)
 
-    # Solver status
-    if m.status[:local_solve] in STATUS_OPT || m.status[:local_solve] in STATUS_LIMIT
+    if Alp.get_option(m, :use_start_as_incumbent)
+        # Check for bound feasibility of the warm start value
+        if !(m.l_var_orig <= m.warm_start_value <= m.u_var_orig)
+            error(
+                "Provide a valid, feasible, warm starting point. Else, set use_start_as_incumbent to false",
+            )
+        end
+
+        obj_warmval = if m.has_nl_objective
+            MOI.eval_objective(m.d_orig, m.warm_start_value)
+        else
+            MOI.Utilities.eval_variables(
+                vi -> m.warm_start_value[vi.value],
+                m.objective_function,
+            )
+        end
+        Alp.update_incumbent(m, obj_warmval, m.warm_start_value)
+        m.status[:local_solve] = MOI.OPTIMIZE_NOT_CALLED
         Alp.get_option(m, :log_level) > 0 && println(
-            "  Local solver returns a feasible point with value $(round(m.best_obj, digits=4))",
+            "  Using warm starting point as a local incumbent solution with value $(round(m.best_obj, digits=4))",
         )
-        Alp.bound_tightening(m, use_bound = true)    # performs bound-tightening with the local solve objective value
+    else
+        Alp.local_solve(m, presolve = true)
+    end
+
+    # Solver status
+    if m.status[:local_solve] in union(STATUS_OPT, STATUS_LIMIT, STATUS_WARM_START)
+        if !Alp.get_option(m, :use_start_as_incumbent)
+            Alp.get_option(m, :log_level) > 0 && println(
+                "  Local solver returns a feasible point with value $(round(m.best_obj, digits=4))",
+            )
+        end
+        Alp.bound_tightening_wrapper(m, use_bound = true)    # performs bound-tightening with the local solve objective value
         Alp.get_option(m, :presolve_bt) && Alp.init_disc(m)            # Re-initialize discretization dictionary on tight bounds
         Alp.get_option(m, :partition_scaling_factor_branch) && (Alp.set_option(
             m,
@@ -221,7 +248,7 @@ function presolve(m::Optimizer)
     elseif m.status[:local_solve] in STATUS_INF
         (Alp.get_option(m, :log_level) > 0) &&
             println("  Bound tightening without objective bounds (OBBT)")
-        Alp.bound_tightening(m, use_bound = false)                      # do bound tightening without objective value
+        Alp.bound_tightening_wrapper(m, use_bound = false)                      # do bound tightening without objective value
         (Alp.get_option(m, :partition_scaling_factor_branch)) && (Alp.set_option(
             m,
             :partition_scaling_factor,
@@ -230,9 +257,8 @@ function presolve(m::Optimizer)
         Alp.get_option(m, :presolve_bt) && Alp.init_disc(m)
 
     elseif m.status[:local_solve] == MOI.INVALID_MODEL
-        @warn " Warning: Presolve ends with local solver yielding $(m.status[:local_solve]). \n This may come from Ipopt's `:Not_Enough_Degrees_Of_Freedom`. \n Consider more replace equality constraints with >= and <= to resolve this."
-
-    else
+        @warn " Warning: Presolve ends with local solver yielding $(m.status[:local_solve]). \n This may come from Ipopt's `:Not_Enough_Degrees_Of_Freedom`."
+    elseif !Alp.get_option(m, :use_start_as_incumbent)
         @warn " Warning: Presolve ends with local solver yielding $(m.status[:local_solve])."
     end
 
@@ -326,8 +352,12 @@ function load_nonlinear_model(m::Optimizer, model::MOI.ModelLike, l_var, u_var)
             m.objective_function,
         )
     end
-    block = MOI.NLPBlockData(m.nl_constraint_bounds_orig, m.d_orig, m.has_nl_objective)
-    MOI.set(model, MOI.NLPBlock(), block)
+
+    if m.d_orig !== nothing
+        block =
+            MOI.NLPBlockData(m.nl_constraint_bounds_orig, m.d_orig, m.has_nl_objective)
+        MOI.set(model, MOI.NLPBlock(), block)
+    end
 
     return x
 end
@@ -393,7 +423,7 @@ function local_solve(m::Optimizer; presolve = false)
     if !presolve
         warmval = m.best_sol[1:m.num_var_orig]
     else
-        warmval = m.initial_warmval[1:m.num_var_orig]
+        warmval = m.warm_start_value[1:m.num_var_orig]
     end
     MOI.set(local_solve_model, MOI.VariablePrimalStart(), x, warmval)
 
@@ -445,7 +475,7 @@ function local_solve(m::Optimizer; presolve = false)
         end
 
         @assert length(candidate_sol) == length(sol_temp)
-        Alp.update_incumb_objective(m, candidate_obj, candidate_sol)
+        Alp.update_incumbent(m, candidate_obj, candidate_sol)
         m.status[:local_solve] = local_nlp_status
         return
 

src/multilinear.jl

@@ -583,7 +583,8 @@ function amp_post_inequalities_cont(
         return
     elseif Alp.get_option(m, :convhull_formulation) == "facet"
         for j in 1:(partition_cnt-1) # Constraint cluster of α >= f(λ)
-            sliced_indices = Alp.collect_indices(λ[ml_indices][:indices], cnt, [1:j;], dim)
+            sliced_indices =
+                Alp.collect_indices(λ[ml_indices][:indices], cnt, [1:j;], dim)
             JuMP.@constraint(
                 m.model_mip,
                 sum(α[var_ind][1:j]) >= sum(λ[ml_indices][:vars][sliced_indices])

src/nlexpr.jl

@@ -353,15 +353,15 @@ function traverse_expr_linear_to_affine(
         return lhscoeffs, lhsvars, rhs, bufferVal, bufferVar
     end
 
-    # HOT-PATCH : Special Structure Recognition
+    # PATCH : Special Structure Recognition
     start_pos = 1
     if (expr.args[1] == :*) && (length(expr.args) == 3)
         if (isa(expr.args[2], Float64) || isa(expr.args[2], Int)) &&
            (expr.args[3].head == :call)
             (coef != 0.0) ? coef = expr.args[2] * coef : coef = expr.args[2]# Patch
             # coef = expr.args[2]
             start_pos = 3
-            @warn "Special expression structure detected [*, coef, :call, ...]. Currently using a beta fix..."
+            # @warn "Special expression structure detected [*, coef, :call, ...]. Currently using a beta fix..."
         end
     end