Guaranteed global nonlinear optimization via interval methods

David P. Sanders

Department of Physics, Faculty of Sciences

Universidad Nacional Autónoma de México (UNAM)

JuMP-dev 2019, Santiago, Chile, 14 March 2019

Outline

  • Motivation and examples
  • Numerical methods for computing with sets: intervals
  • Global nonlinear optimization: spatial interval branch and bound (exhaustive search)
  • Branch and prune: Constraint propagation
In [2]:
# ] add IntervalArithmetic#fastpowers
# ] add IntervalConstraintProgramming#separators_from_contractors  # ?
# ] add ModelingToolkit#master

using IntervalArithmetic, IntervalOptimisation
using IntervalConstraintProgramming, IntervalRootFinding
using StaticArrays, ModelingToolkit, LinearAlgebra

include("constraints.jl"); include("unify_solutions.jl")
Out[2]:
unify_solutions (generic function with 1 method)

Goal

  • Find global minimum and minimizer(s) of $f: \mathbb{R}^n \to \mathbb{R}$

    i.e. $\mathbf{x}^* \colon\ f(\mathbf{x}^*) \le f(\mathbf{x}) \quad \forall \mathbf{x} \in \mathbb{R}^n $

  • ...with bounds that are guaranteed to be correct (mathematically rigorous), including bounding all floating-point errors
  • Solution: Calculate with sets or real numbers (instead of single numbers)
  • Give interval bound (lower..upper) on the range of $f$ over a set $X$
In [2]:
G(X) = 1 + sum(X.^2) / 4000 - prod( cos(X[i] / i) for i in 1:length(X) )   # Griewank
h(x) = G(x - 1.5) * G(x-1) * G(x+1) * G(x + 1.5)
Out[2]:
h (generic function with 1 method)
In [2]:
using Plots; gr()
Out[2]:
Plots.GRBackend()
In [3]:
p = plot(h, -50:0.001:50, leg=false)
UndefVarError: h not defined

Stacktrace:
 [1] top-level scope at In[3]:1
In [5]:
@format midpoint 3

using ForwardDiff

@time rts = IntervalRootFinding.roots(x->ForwardDiff.derivative(h, x), -50..50);  rts
  2.600955 seconds (8.69 M allocations: 366.117 MiB, 8.16% gc time)
Out[5]:
75-element Array{Root{Interval{Float64}},1}:
 Root(49.3 ± 7.11e-15, :unique) 
 Root(47.2 ± 7.11e-15, :unique) 
 Root(45 ± 7.11e-15, :unique)   
 Root(44 ± 7.11e-15, :unique)   
 Root(42.9 ± 7.11e-15, :unique) 
 Root(40.9 ± 7.11e-15, :unique) 
 Root(38.7 ± 7.11e-15, :unique) 
 Root(37.7 ± 7.11e-15, :unique) 
 Root(36.6 ± 7.11e-15, :unique) 
 Root(34.6 ± 7.11e-15, :unique) 
 Root(32.5 ± 7.11e-15, :unique) 
 Root(31.4 ± 1.07e-14, :unique) 
 Root(30.3 ± 3.56e-15, :unique) 
 ⋮                              
 Root(-31.4 ± 3.56e-15, :unique)
 Root(-32.5 ± 7.11e-15, :unique)
 Root(-34.6 ± 7.11e-15, :unique)
 Root(-36.6 ± 7.11e-15, :unique)
 Root(-37.7 ± 7.11e-15, :unique)
 Root(-38.7 ± 7.11e-15, :unique)
 Root(-40.9 ± 7.11e-15, :unique)
 Root(-42.9 ± 7.11e-15, :unique)
 Root(-44 ± 7.11e-15, :unique)  
 Root(-45 ± 7.11e-15, :unique)  
 Root(-47.2 ± 7.11e-15, :unique)
 Root(-49.3 ± 7.11e-15, :unique)
In [6]:
mids = mid.(interval.(rts))
second_deriv(f, x) = ForwardDiff.derivative(x->ForwardDiff.derivative(f, x), x)

scatter!(mids, h.(mids), marker_z = second_deriv.(h, mids))
Out[6]:
-40 -20 0 20 40 0.0 2.5 5.0 7.5 10.0