psorlab / eago-notebooks Goto Github PK
View Code? Open in Web Editor NEWA collection of Jupyter notebooks which use the EAGO package
A collection of Jupyter notebooks which use the EAGO package
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.
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 saysFor 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.
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 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
xvalue = Interval.(lower, upper)
F = sin(xvalue[1])*xvalue[2]^2-cos(xvalue[3])/xvalue[4]
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:
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.
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)
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?
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.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.A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.