Global Optimisation

REQUIRE

@@ -2,3 +2,5 @@ julia 0.5
 
 IntervalArithmetic 0.9.1
 IntervalRootFinding 0.1
+IntervalConstraintProgramming
+DataStructures

examples/optimise1.jl

@@ -0,0 +1,56 @@
+using IntervalOptimisation, IntervalArithmetic, IntervalConstraintProgramming
+
+
+# Unconstrained Optimisation
+# Example from Eldon Hansen - Global Optimisation using Interval Analysis, Chapter 11
+f(x, y) = (1.5 - x * (1 - y))^2 + (2.25 - x * (1 - y^2))^2 + (2.625 - x * (1 - y^3))^2
+# f (generic function with 2 methods)
+
+f(X) = f(X...)
+# f (generic function with 2 methods)
+
+X = (-1e6..1e6) × (-1e6..1e6)
+# [-1e+06, 1e+06] × [-1e+06, 1e+06]
+
+minimise_icp(f, X)
+# ([0, 2.39527e-07], IntervalArithmetic.IntervalBox{2,Float64}[[2.99659, 2.99748] × [0.498906, 0.499523], [2.99747, 2.99834] × [0.499198, 0.500185], [2.99833, 2.99922] × [0.499198, 0.500185], [2.99921, 3.00011] × [0.499621, 0.500415], [3.0001, 3.001] × [0.499621, 0.500415], [3.00099, 3.0017] × [0.500169, 0.500566], [3.00169, 3.00242] × [0.500169, 0.500566]])
+
+
+
+f(X) = X[1]^2 + X[2]^2
+# f (generic function with 2 methods)
+
+X = (-∞..∞) × (-∞..∞)
+# [-∞, ∞] × [-∞, ∞]
+
+minimise_icp(f, X)
+# ([0, 0], IntervalArithmetic.IntervalBox{2,Float64}[[-0, 0] × [-0, 0], [0, 0] × [-0, 0]])
+
+
+
+
+# Constrained Optimisation
+# Example 1 from http://archimedes.cheme.cmu.edu/?q=baronexamples
+f(X) = -1 * (X[1] + X[2])
+# f (generic function with 1 method)
+
+X = (-∞..∞) × (-∞..∞)
+# [-∞, ∞] × [-∞, ∞]
+
+constraints = [IntervalOptimisation.constraint(x->(x[1]), -∞..6), IntervalOptimisation.constraint(x->x[2], -∞..4), IntervalOptimisation.constraint(x->(x[1]*x[2]), -∞..4)]
+# 3-element Array{IntervalOptimisation.constraint{Float64},1}:
+#  IntervalOptimisation.constraint{Float64}(#3, [-∞, 6])
+#  IntervalOptimisation.constraint{Float64}(#4, [-∞, 4])
+#  IntervalOptimisation.constraint{Float64}(#5, [-∞, 4])
+
+minimise_icp_constrained(f, X, constraints)
+# ([-6.66676, -6.66541], IntervalArithmetic.IntervalBox{2,Float64}[[5.99918, 6] × [0.666233, 0.666758], [5.99918, 6] × [0.665717, 0.666234], [5.99887, 5.99919] × [0.666233, 0.666826], [5.99984, 6] × [0.665415, 0.665718], [5.99856, 5.99888] × [0.666233, 0.666826], [5.99969, 5.99985] × [0.665415, 0.665718]])
+
+
+
+# One-dimensional case
+minimise1d(x -> (x^2 - 2)^2, -10..11)
+# ([0, 1.33476e-08], IntervalArithmetic.Interval{Float64}[[-1.41426, -1.41358], [1.41364, 1.41429]])
+
+minimise1d_deriv(x -> (x^2 - 2)^2, -10..11)
+# ([0, 8.76812e-08], IntervalArithmetic.Interval{Float64}[[-1.41471, -1.41393], [1.41367, 1.41444]])

src/IntervalOptimisation.jl

@@ -3,12 +3,14 @@
 module IntervalOptimisation
 
 export minimise, maximise,
-       minimize, maximize
+       minimize, maximize,
+       minimise1d, minimise1d_deriv,
+       minimise_icp, minimise_icp_constrained
 
 include("SortedVectors.jl")
 using .SortedVectors
 
-using IntervalArithmetic, IntervalRootFinding
+using IntervalArithmetic, IntervalRootFinding, DataStructures, IntervalConstraintProgramming, StaticArrays, ForwardDiff
 
 if !isdefined(:sup)
     const sup = supremum
@@ -19,8 +21,9 @@ if !isdefined(:inf)
 end
 
 include("optimise.jl")
+include("optimise1.jl")
 
 const minimize = minimise
-const maximize = maximise 
+const maximize = maximise
 
 end

src/SortedVectors.jl

@@ -12,7 +12,7 @@ immutable SortedVector{T, F<:Function}
     data::Vector{T}
     by::F
 
-    function SortedVector(data::Vector{T}, by::F)
+    function SortedVector{T,F}(data::Vector{T}, by::F) where {T,F}
         new(sort(data), by)
     end
 end

src/optimise1.jl

@@ -0,0 +1,194 @@
+struct IntervalMinimum{T<:Real}
+    interval::Interval{T}
+    minimum::T
+end
+
+Base.isless(a::IntervalMinimum{T}, b::IntervalMinimum{T}) where {T<:Real} = isless(a.minimum, b.minimum)
+
+function minimise1d(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3, use_deriv=false, use_second_deriv=false) where {T<:Real}
+
+    Q = binary_minheap(IntervalMinimum{T})
+
+    global_minimum = f(interval(mid(x))).hi
+
+    arg_minima = Interval{T}[]
+
+    push!(Q, IntervalMinimum(x, global_minimum))
+
+    while !isempty(Q)
+
+        p = pop!(Q)
+
+        if isempty(p.interval)
+            continue
+        end
+
+        if p.minimum > global_minimum
+            continue
+        end
+
+        if use_deriv
+            deriv = ForwardDiff.derivative(f, p.interval)
+            if 0 ∉ deriv
+                continue
+            end
+        end
+
+        # Second derivative contractor
+        if use_second_deriv
+            doublederiv = ForwardDiff.derivative(x->ForwardDiff.derivative(f, x), p.interval)
+            if doublederiv < 0
+                continue
+            end
+        end
+
+        m = mid(p.interval)
+        current_minimum = f(interval(m)).hi
+
+        if current_minimum < global_minimum
+            global_minimum = current_minimum
+        end
+
+
+        # Contractor 1
+        if use_deriv
+            x = m .+ extended_div((interval(-∞, global_minimum) - f(m)), deriv)
+
+            x = x .∩ p.interval
+        end
+
+        if diam(p.interval) < abstol
+            push!(arg_minima, p.interval)
+        else
+
+            if use_deriv && isempty(x[2])
+                x1, x2 = bisect(x[1])
+                push!(Q, IntervalMinimum(x1, f(x1).lo), IntervalMinimum(x2, f(x2).lo))
+                continue
+            end
+
+            x1, x2 = bisect(p.interval)
+            push!(Q, IntervalMinimum(x1, f(x1).lo), IntervalMinimum(x2, f(x2).lo))
+        end
+    end
+    lb = minimum(inf.(f.(arg_minima)))
+
+    return lb..global_minimum, arg_minima
+end
+
+
+struct IntervalBoxMinimum{N, T<:Real}
+    interval::IntervalBox{N, T}
+    minimum::T
+end
+
+"""
+Datatype to provide constraints to Global Optimisation such as:
+```
+Constraint(x->(x^2 - 10), -∞..1)
+```
+"""
+struct Constraint{T<:Real}
+    f::Function
+    c::Interval{T}
+end
+
+Base.isless(a::IntervalBoxMinimum{N, T}, b::IntervalBoxMinimum{N, T}) where {N, T<:Real} = isless(a.minimum, b.minimum)
+
+function minimise_icp(f::Function, x::IntervalBox{N, T}; reltol=1e-3, abstol=1e-3) where {N, T<:Real}
+
+    Q = binary_minheap(IntervalBoxMinimum{N, T})
+
+    global_minimum = ∞
+
+    x = icp(f, x, -∞..global_minimum)
+
+    arg_minima = IntervalBox{N, T}[]
+
+    push!(Q, IntervalBoxMinimum(x, global_minimum))
+
+    while !isempty(Q)
+
+        p = pop!(Q)
+
+        if isempty(p.interval)
+            continue
+        end
+
+        if p.minimum > global_minimum
+            continue
+        end
+
+        current_minimum = f(interval.(mid(p.interval))).hi
+
+        if current_minimum < global_minimum
+            global_minimum = current_minimum
+        end
+
+        X = icp(f, p.interval, -∞..global_minimum)
+
+        if diam(p.interval) < abstol
+            push!(arg_minima, p.interval)
+
+        else
+            x1, x2 = bisect(X)
+            push!(Q, IntervalBoxMinimum(x1, f(x1).lo), IntervalBoxMinimum(x2, f(x2).lo))
+        end
+    end
+
+    lb = minimum(inf.(f.(arg_minima)))
+
+    return lb..global_minimum, arg_minima
+end
+
+function minimise_icp_constrained(f::Function, x::IntervalBox{N, T}, constraints::Vector{Constraint{T}} = Vector{Constraint{T}}(); reltol=1e-3, abstol=1e-3) where {N, T<:Real}
+
+    Q = binary_minheap(IntervalBoxMinimum{N, T})
+
+    global_minimum = ∞
+
+    for t in constraints
+        x = icp(t.f, x, t.c)
+    end
+
+    x = icp(f, x, -∞..global_minimum)
+
+    arg_minima = IntervalBox{N, T}[]
+
+    push!(Q, IntervalBoxMinimum(x, global_minimum))
+
+    while !isempty(Q)
+
+        p = pop!(Q)
+
+        if isempty(p.interval)
+            continue
+        end
+
+        if p.minimum > global_minimum
+            continue
+        end
+
+        # current_minimum = f(interval.(mid(p.interval))).hi
+        current_minimum = f(p.interval).hi
+
+        if current_minimum < global_minimum
+            global_minimum = current_minimum
+        end
+
+        X = icp(f, p.interval, -∞..global_minimum)
+
+        for t in constraints
+            X = icp(t.f, X, t.c)
+        end
+
+        if diam(p.interval) < abstol
+            push!(arg_minima, p.interval)
+        else
+            x1, x2 = bisect(X)
+            push!(Q, IntervalBoxMinimum(x1, f(x1).lo), IntervalBoxMinimum(x2, f(x2).lo))
+        end
+    end
+    lb = minimum(inf.(f.(arg_minima)))
+    return lb..global_minimum, arg_minima
+end