Global Optimization via NLopt

The build_loss_objective function builds an objective function compatible with MathOptInterface-associated solvers. This includes packages like IPOPT, NLopt, MOSEK, etc. Building off of the previous example, we can build a cost function for the single parameter optimization problem like:

using DifferentialEquations, DiffEqParamEstim, Optimization, OptimizationMOI,
      OptimizationNLopt, NLopt, RecursiveArrayTools
using Random

function f(du, u, p, t)
    du[1] = p[1] * u[1] - u[1] * u[2]
    du[2] = -3 * u[2] + u[1] * u[2]
end

u0 = [1.0; 1.0]
tspan = (0.0, 10.0)
p = [1.5]
prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob, Tsit5())

t = collect(range(0, stop = 10, length = 200))
Random.seed!(1234)
randomized = VectorOfArray([(sol(t[i]) + 0.01randn(2)) for i in 1:length(t)])
data = convert(Array, randomized)

obj = build_loss_objective(prob, Tsit5(), L2Loss(t, data), Optimization.AutoForwardDiff())

SciMLBase.OptimizationFunction{true, ADTypes.AutoForwardDiff{nothing, Nothing}, DiffEqParamEstim.var"#37#38"{Nothing, typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR), Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}, SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEqTsit5.Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), FastBroadcast.Serial}, L2Loss{Vector{Float64}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing}, Nothing, Tuple{}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}(DiffEqParamEstim.var"#37#38"{Nothing, typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR), Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}, SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEqTsit5.Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), FastBroadcast.Serial}, L2Loss{Vector{Float64}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing}, Nothing, Tuple{}}(nothing, DiffEqParamEstim.STANDARD_PROB_GENERATOR, Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}(), SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}, SciMLBase.StandardODEProblem}(SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}(Main.f, LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, SciMLBase.DEFAULT_OBSERVED, nothing, nothing, nothing, nothing), [1.0, 1.0], (0.0, 10.0), [1.5], Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}(), SciMLBase.StandardODEProblem()), OrdinaryDiffEqTsit5.Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), FastBroadcast.Serial}(OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, FastBroadcast.Serial()), L2Loss{Vector{Float64}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing}([0.0, 0.05025125628140704, 0.10050251256281408, 0.1507537688442211, 0.20100502512562815, 0.25125628140703515, 0.3015075376884422, 0.35175879396984927, 0.4020100502512563, 0.45226130653266333  …  9.547738693467336, 9.597989949748744, 9.64824120603015, 9.698492462311558, 9.748743718592964, 9.798994974874372, 9.849246231155778, 9.899497487437186, 9.949748743718592, 10.0], [0.9964027109317651 1.0237452817740549 … 0.9884636708095201 1.033500062415095; 1.0108720849242858 0.9121858132876625 … 1.0081722075883954 0.9044870156093264], nothing, nothing, nothing, nothing, nothing), nothing, ()), ADTypes.AutoForwardDiff(), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, SciMLBase.DEFAULT_OBSERVED_NO_TIME, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing)

You can either use the NLopt package directly or through either the OptimizationNLopt or OptimizationMOI which provides an interface for all MathOptInterface compatible non-linear solvers.

We can now use this obj as the objective function with MathProgBase solvers. For our example, we will use NLopt. To use the local derivative-free Constrained Optimization BY Linear Approximations algorithm, we can simply do:

opt = Opt(:LN_COBYLA, 1)
maxeval!(opt, 10)
optprob = Optimization.OptimizationProblem(obj, [1.3])
res = solve(optprob, opt)
@assert isapprox(res.u, p; atol = 0.05)

For a modified evolutionary algorithm, we can use:

opt = Opt(:GN_ESCH, 1)
lower_bounds!(opt, [0.0])
upper_bounds!(opt, [5.0])
xtol_rel!(opt, 1e-3)
maxeval!(opt, 100_000)
NLopt.srand(1234)
res = solve(optprob, opt)
@assert isapprox(res.u, p; atol = 0.05)

We can even use things like the Improved Stochastic Ranking Evolution Strategy (and add constraints if needed). Let's use this through OptimizationMOI:

optprob = Optimization.OptimizationProblem(obj, [0.2], lb = [-1.0], ub = [5.0])
res = solve(optprob,
    OptimizationMOI.MOI.OptimizerWithAttributes(NLopt.Optimizer,
        "algorithm" => :GN_ISRES,
        "xtol_rel" => 1e-3,
        "seed" => 1234,
        "maxeval" => 10_000))
@assert isapprox(res.u, p; atol = 0.05)

which is very robust to the initial condition. We can also directly use the NLopt interface as below. The fastest result comes from the following algorithm choice:

opt = Opt(:LN_BOBYQA, 1)
lower_bounds!(opt, [0.0])
upper_bounds!(opt, [5.0])
xtol_rel!(opt, 1e-6)
maxeval!(opt, 100)
min_objective!(opt, (x, _) -> obj(x))
minf, minx, ret = NLopt.optimize(opt, [1.3])
@assert isapprox(minx, p; atol = 0.01)
(minf, minx, ret)

(0.03895090457147877, [1.500007012746243], :XTOL_REACHED)

For more information, see the NLopt documentation for more details. And give IPOPT or MOSEK a try!