psorlab / eago-notebooks Goto Github PK

Hi, @mewilhel, this issue follows the previous one #3, but I think it is more appropriate to open two issues.

I am new to McCormick relaxation and learn the basics from the paper [1], especially the affine relaxation with the subgradient. I also read your paper on EAGO [2]. I played with the notebook demonstrating customized routines in EAGO, since I have a similar parameter identification problem, i.e., a nonlinear least-squares problem with only box constraints. I have however some unexpected findings and would like to consult you for more details. The following discussion is constrained to box constrained problems only. The modified notebook can be found here.

1. Affine relaxation in lower_problem!
In lower_problem! of your original notebook, the lower bounding is underestimated by the natural interval extension. However, it is often stated in the literature (e.g., [1, 2]) that the natural interval extension is not as tight as McCormick relaxation. I thus examined the simplest single-point affine relaxation in lower_problem! (thanks for the excellent McCormick.jl package). Additionally, I also check which one is a tighter underestimation. The code is pasted below (hope nothing wrong).

import EAGO.lower_problem!
function lower_problem!(t::IntervalExt, x::EAGO.Optimizer)
    
    # retrieve bounds at current node
    n = x._current_node
    lower = n.lower_variable_bounds
    upper = n.upper_variable_bounds
    
    # McCormick relaxation
    box = Interval.(lower, upper)
    xmid = mid.(box)
    xMC = [MC{4, NS}(xmid[i], box[i], i) for i =1:4]
    fMC = obj_func(xMC)
    
    # affine relaxation at `xmid`
    s = fMC.cv_grad
    xopt = similar(lower)  # store minimum solution to the affine relaxation
    @inbounds for (i, ele) in enumerate(s)
        xopt[i] = ele > 0 ? lower[i] : upper[i]
    end
    fopt = fMC.cv + s ⋅ (xopt .- xmid)
    
    # report lower bounding results
    x._lower_objective_value = fopt
    x._lower_solution = xopt
    x._lower_feasibility = true
    x._cut_add_flag = false
    
    # can we get a tigher lower bound than natural interval extension?
    itv_nat = fMC.Intv
    if fopt > lo(itv_nat)
        @show (lo(itv_nat), fopt, n.depth, n.id)
    end
    return
end

The above obj_func is just obj_func(x) = sin(x[1])*x[2]^2 - cos(x[3])/x[4] for simplicity.

• It is surprising that fopt > lo(itv_nat) never becomes true, which means the natural interval extension is always tighter for lower bounding. However, in Section 3.6 of [2], it says

For nonlinear expressions, an affine relaxation is generated via an affine approximation of the expression at the midpoint
of the domain using subgradient information.

Since I am new to this field, does it often happen that the affine approximate is even worse than natural extension? If so, we may combine the two and fetch the smaller one (I am not sure whether EAGO has already done it).

In addition, is it worth attempting multi-point affine relaxation in contrast to the single midpoint as done above? For constrained problems, multi-point affine relaxation leads to a larger LP problem; for non-constrained problems as we discuss here, the minimum value of the relaxation is min max g_i(x) where g_i is an affine function. There exists faster algorithm than LP to solve this minimization.

• Though there happens to be always fopt <= lo(itv_nat), it actually took less iterations than the original method using natural interval extension (557 vs. 734), which seems strange to me. A tighter lower bounding is expected to reduce rather than increase the number of iterations in branch and bound, right? Is it due to the branching mechanism in EAGO? (since the domain reduction techniques have been disabled and the upper bounding method is the same)

2. Inspection of the obtained global optimal solutions

In #3, I noticed that there is discrimination between the reported result and the actual one. This issue still exists here:

The optimal value is:    -1.498837828202948, the solution found is [-1.573774954304099, 0.9998533502221107, -0.006667515262961388, 2.003519594669342].
Verify obj value @ xsol: -1.4988128398246658

• The digits differ after 1e-5 as #3.

The above observations may, unfortunately, imply that EAGO has numerical issues somewhere, which asks for improvement.

3. Finally, as stated in [1], the McCormick based relaxation is particularly suitable for a parameter identification problem (like the example in Section 5.1 of [1]), where the number of variables is small but the number of factors can be large due to multiple data points. For such an unconstrained NLP, do you have suggestions on the lower_problem! based on your experience? I have not tested the multi-point affine relaxation yet.

[1] Mitsos, Alexander, Benoît Chachuat, and Paul I. Barton. "McCormick-based relaxations of algorithms." SIAM Journal on Optimization 20.2 (2009): 573-601.

[2] Wilhelm, M. E., and M. D. Stuber. "EAGO. jl: easy advanced global optimization in Julia." Optimization Methods and Software (2020): 1-26.

I ran into two issues when reproducing nlpopt_interval_bnb example

1. In "lower_problem!", x is defined twice (an EAGO.Optimizer and an array of intervals). I changed the latter one as

xvalue = Interval.(lower, upper)
F = sin(xvalue[1])*xvalue[2]^2-cos(xvalue[3])/xvalue[4]

2. After resolving the previous issue, I got

ERROR: LoadError: UndefVarError: intv_lo not defined
Stacktrace:
 [1] cut_condition(::IntervalExt, ::Optimizer{GLPK.Optimizer,Ipopt.Optimizer}) at ......\.julia\packages\EAGO\Fr4Op\src\eago_optimizer\subroutines.jl:718

Versions:

• Julia: 1.1
• EAGO: v0.3.0

Received
ERROR: LoadError: UndefVarError: NS not defined
while running the example code in
https://github.com/PSORLab/EAGO-notebooks/blob/master/notebooks/relaxation_plot_1d_explicit.ipynb

OS: Windows 10
Julia Version: 1.2.0
Package versions:
EAGO v0.3.0
IntervalArithmetic v0.16.1

Hi, @mewilhel, thanks for your excellent work. I played with this notebook, as indicated in the title, and found some problems that need your help.

1. x._lower_solution = IntervalArithmetic.mid.(x_value) in lower_problem!
   Though there is no documentation for _lower_solution, I guess it is the optimal solution to the lower bound problem. In the line above it x._lower_objective_value = F.lo, the minimum objective value is set the lower bound of the natural interval extension. However, how do we know that this value is attained by mid.(x_value)? Of course, from my limited knowledge of branch and bound (BB), the actual solution of the lower bound problem may not matter (but it may be used in branching though I am not sure how EAGO does)

2. After optimization, the actual objective value differs from the reported one.
   I changed the last cell to (added the last line)

fval = JuMP.objective_value(m)
xsol = JuMP.value.(x)
status_term = JuMP.termination_status(m)
status_prim = JuMP.primal_status(m)

println("EAGO terminated with a status of $status_term and a result code of $status_prim")
println("The optimal value is: $fval, the solution found is $xsol.")

println("Verify objective value: ", sin(xsol[1])xsol[2]^2-cos(xsol[3])/xsol[4])

The output is

EAGO terminated with a status of OPTIMAL and a result code of FEASIBLE_POINT
The optimal value is: -1.498986093897677, the solution found is [-1.56982421875, -0.999755859375, -0.00244140625, 2.002197265625].
Verify objective value: -1.4989611040367379

Note that JuMP.objective_value(m) leads to -1.49898 but the one computed using the optimal solution xsol is -1.49896. What is the reason for this perceptible difference that cannot be ignored?

3. About logging
   In EAGO.EAGOParameters, there is a field output_iterations which controls "Display summary of iteration to console every output_iterations (default = 10)". It seems the correct default value is actually 1000, and that's why no display in this notebook.
   After setting opt_dict[:output_iterations] = 10, I saw the display every 10 iterations. My question is how we can access this logging data instead of just printing it. For example, I want to plot the convergence curve of the lower and upper bounds during optimization after the optimization finishes.