From 2a78d4977c8be2c95119f11f5b0372471eae3036 Mon Sep 17 00:00:00 2001 From: Vaibhav Dixit Date: Mon, 12 Aug 2024 10:03:07 -0400 Subject: [PATCH] Add constraints support for NLopt --- lib/OptimizationNLopt/Project.toml | 1 + .../src/OptimizationNLopt.jl | 67 +++++++++++++++++-- lib/OptimizationNLopt/test/runtests.jl | 50 ++++++++++++-- 3 files changed, 106 insertions(+), 12 deletions(-) diff --git a/lib/OptimizationNLopt/Project.toml b/lib/OptimizationNLopt/Project.toml index c2aedb1e4..64ab699c2 100644 --- a/lib/OptimizationNLopt/Project.toml +++ b/lib/OptimizationNLopt/Project.toml @@ -6,6 +6,7 @@ version = "0.2.2" [deps] NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd" Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" [compat] diff --git a/lib/OptimizationNLopt/src/OptimizationNLopt.jl b/lib/OptimizationNLopt/src/OptimizationNLopt.jl index fe2eb9abf..660034f43 100644 --- a/lib/OptimizationNLopt/src/OptimizationNLopt.jl +++ b/lib/OptimizationNLopt/src/OptimizationNLopt.jl @@ -36,6 +36,24 @@ function SciMLBase.requiresconsjac(opt::NLopt.Algorithm) #https://github.com/Jul end end +function SciMLBase.allowsconstraints(opt::NLopt.Algorithm) + str_opt = string(opt) + if occursin("AUGLAG", str_opt) || occursin("CCSA", str_opt) || occursin("MMA", str_opt) || occursin("COBYLA", str_opt) || occursin("ISRES", str_opt) || occursin("AGS", str_opt) || occursin("ORIG_DIRECT", str_opt) || occursin("SLSQP", str_opt) + return true + else + return false + end +end + +function SciMLBase.__init(prob::SciMLBase.OptimizationProblem, opt::NLopt.Algorithm, + data = Optimization.DEFAULT_DATA; cons_tol = 1e-6, + callback = (args...) -> (false), + progress = false, kwargs...) + return OptimizationCache(prob, opt, data; cons_tol, callback, progress, + kwargs...) +end + + function __map_optimizer_args!(cache::OptimizationCache, opt::NLopt.Opt; callback = nothing, maxiters::Union{Number, Nothing} = nothing, @@ -76,7 +94,9 @@ function __map_optimizer_args!(cache::OptimizationCache, opt::NLopt.Opt; # add optimiser options from kwargs for j in kwargs - eval(Meta.parse("NLopt." * string(j.first) * "!"))(opt, j.second) + if j.first != :cons_tol + eval(Meta.parse("NLopt." * string(j.first) * "!"))(opt, j.second) + end end if cache.ub !== nothing @@ -170,14 +190,18 @@ function SciMLBase.__solve(cache::OptimizationCache{ return x[1] end - fg! = function (θ, G) - if length(G) > 0 - cache.f.grad(G, θ) + if !hasfield(typeof(cache.f), :fg) || cache.f.fg === nothing + fg! = function (θ, G) + if length(G) > 0 + cache.f.grad(G, θ) + end + return _loss(θ) end - - return _loss(θ) + else + fg! = cache.f.fg end + opt_setup = if isa(cache.opt, NLopt.Opt) if ndims(cache.opt) != length(cache.u0) error("Passed NLopt.Opt optimization dimension does not match OptimizationProblem dimension.") @@ -193,6 +217,37 @@ function SciMLBase.__solve(cache::OptimizationCache{ NLopt.min_objective!(opt_setup, fg!) end + if cache.f.cons !== nothing + eqinds = map((y) -> y[1]==y[2], zip(cache.lcons, cache.ucons)) + ineqinds = map((y) -> y[1]!=y[2], zip(cache.lcons, cache.ucons)) + if sum(ineqinds) > 0 + ineqcons = function (res, θ, J) + cons_cache = zeros(eltype(res), sum(eqinds)+sum(ineqinds)) + cache.f.cons(cons_cache, θ) + res .= @view(cons_cache[ineqinds]) + if length(J) > 0 + Jcache = zeros(eltype(J), sum(ineqinds)+sum(eqinds), length(θ)) + cache.f.cons_j(Jcache, θ) + J .= @view(Jcache[ineqinds, :])' + end + end + NLopt.inequality_constraint!(opt_setup, ineqcons, [cache.solver_args.cons_tol for i in 1:sum(ineqinds)]) + end + if sum(eqinds) > 0 + eqcons = function (res, θ, J) + cons_cache = zeros(eltype(res), sum(eqinds)+sum(ineqinds)) + cache.f.cons(cons_cache, θ) + res .= @view(cons_cache[eqinds]) + if length(J) > 0 + Jcache = zeros(eltype(res), sum(eqinds)+sum(ineqinds), length(θ)) + cache.f.cons_j(Jcache, θ) + J .= @view(Jcache[eqinds, :])' + end + end + NLopt.equality_constraint!(opt_setup, eqcons, [cache.solver_args.cons_tol for i in 1:sum(eqinds)]) + end + end + maxiters = Optimization._check_and_convert_maxiters(cache.solver_args.maxiters) maxtime = Optimization._check_and_convert_maxtime(cache.solver_args.maxtime) diff --git a/lib/OptimizationNLopt/test/runtests.jl b/lib/OptimizationNLopt/test/runtests.jl index f48a067f8..cbd79d475 100644 --- a/lib/OptimizationNLopt/test/runtests.jl +++ b/lib/OptimizationNLopt/test/runtests.jl @@ -1,5 +1,5 @@ using OptimizationNLopt, Optimization, Zygote -using Test +using Test, Random @testset "OptimizationNLopt.jl" begin rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2 @@ -16,7 +16,7 @@ using Test optprob = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) prob = OptimizationProblem(optprob, x0, _p) - sol = solve(prob, NLopt.Opt(:LN_BOBYQA, 2)) + sol = solve(prob, NLopt.Opt(:LD_LBFGS, 2)) @test sol.retcode == ReturnCode.Success @test 10 * sol.objective < l1 @@ -26,10 +26,6 @@ using Test @test sol.retcode == ReturnCode.Success @test 10 * sol.objective < l1 - sol = solve(prob, NLopt.Opt(:LD_LBFGS, 2)) - @test sol.retcode == ReturnCode.Success - @test 10 * sol.objective < l1 - sol = solve(prob, NLopt.Opt(:G_MLSL_LDS, 2), local_method = NLopt.Opt(:LD_LBFGS, 2), maxiters = 10000) @test sol.retcode == ReturnCode.MaxIters @@ -82,4 +78,46 @@ using Test #nlopt gives the last best not the one where callback stops @test sol.objective < 0.8 end + + @testset "constrained" begin + cons = (res, x, p) -> res .= [x[1]^2 + x[2]^2 - 1.0] + x0 = zeros(2) + optprob = OptimizationFunction(rosenbrock, Optimization.AutoZygote(); + cons = cons) + prob = OptimizationProblem(optprob, x0, _p, lcons = [0.0], ucons = [0.0]) + sol = solve(prob, NLopt.LN_COBYLA()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + + Random.seed!(1) + prob = OptimizationProblem(optprob, rand(2), _p, + lcons = [0.0], ucons = [0.0]) + + sol = solve(prob, NLopt.LD_SLSQP()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + + Random.seed!(1) + prob = OptimizationProblem(optprob, rand(2), _p, + lcons = [0.0], ucons = [0.0]) + sol = solve(prob, NLopt.AUGLAG(), local_method = NLopt.LD_LBFGS()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + + function con2_c(res, x, p) + res .= [x[1]^2 + x[2]^2 - 1.0, x[2] * sin(x[1]) - x[1] - 2.0] + end + + optprob = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff();cons = con2_c) + Random.seed!(1) + prob = OptimizationProblem(optprob, rand(2), _p, lcons = [0.0, -Inf], ucons = [0.0, 0.0]) + sol = solve(prob, NLopt.LD_AUGLAG(), local_method = NLopt.LD_LBFGS()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + + prob = OptimizationProblem(optprob, rand(2), _p, lcons = [-Inf, -Inf], ucons = [0.0, 0.0], lb = [-1.0, -1.0], ub = [1.0, 1.0]) + sol = solve(prob, NLopt.GN_ISRES(), maxiters = 1000) + @test sol.retcode == ReturnCode.MaxIters + @test 10 * sol.objective < l1 + end end