Handle constraints by evaluating jac and hess once

SciML · Oct 20, 2023 · 14812ba · 14812ba
1 parent 4357b5d
commit 14812ba
Show file tree

Hide file tree

Showing 3 changed files with 39 additions and 32 deletions.
diff --git a/lib/OptimizationPRIMA/Project.toml b/lib/OptimizationPRIMA/Project.toml
@@ -4,7 +4,6 @@ authors = ["Vaibhav Dixit <[email protected]> and contributors"]
 version = "1.0.0-DEV"
 
 [deps]
-ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 PRIMA = "0a7d04aa-8ac2-47b3-b7a7-9dbd6ad661ed"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

diff --git a/lib/OptimizationPRIMA/src/OptimizationPRIMA.jl b/lib/OptimizationPRIMA/src/OptimizationPRIMA.jl
@@ -16,17 +16,6 @@ SciMLBase.allowsconstraints(::Union{LINCOA, COBYLA}) = true
 SciMLBase.allowsbounds(opt::Union{BOBYQA, LINCOA, COBYLA}) = true
 SciMLBase.requiresconstraints(opt::COBYLA) = true
 
-function SciMLBase.__init(prob::SciMLBase.OptimizationProblem, opt::PRIMASolvers,
-    data = Optimization.DEFAULT_DATA;
-    callback = (args...) -> (false),
-    progress = false, kwargs...)
-    if opt isa COBYLA
-        sys = ModelingToolkit.modelingtoolkitize(prob)
-        prob = OptimizationProblem(sys, prob.u0, prob.p; lb = prob.lb, ub = prob.ub, cons_sparse = true)
-    end
-    return OptimizationCache(prob, opt, data; callback, progress,
-        kwargs...)
-end
 
 function get_solve_func(opt::PRIMASolvers)
     if opt isa UOBYQA
@@ -132,33 +121,34 @@ function SciMLBase.__solve(cache::OptimizationCache{
 
     t0 = time()
     if cache.opt isa COBYLA
-        lininds = Int[]
+        lineqsinds = Int[]
+        linineqsinds = Int[]
         nonlininds = Int[]
-        symboliclinexpr = []
-        for i in eachindex(cache.f.cons_expr)
-            println(cache.f.cons_expr[i])
-            println(isempty(cache.f.cons_hess_prototype[i]))
-            if isempty(cache.f.cons_hess_prototype[i])
-                push!(lininds, i)
-                symbolicexpr = Symbolics.parse_expr_to_symbolic(cache.f.cons_expr[i], (@__MODULE__,))
-                println(symbolicexpr)
-                push!(symboliclinexpr, symbolicexpr)
+        H = [zeros(length(cache.u0), length(cache.u0)) for i in 1:length(cache.lcons)]
+        J = zeros(length(cache.lcons), length(cache.u0))
+
+        cache.f.cons_h(H, ones(length(cache.u0)))
+        cache.f.cons_j(J, ones(length(cache.u0)))
+        for i in eachindex(cache.lcons)
+            if iszero(H[i]) && cache.lcons[i] == cache.ucons[i]
+                push!(lineqsinds, i)
+            elseif iszero(H[i]) && cache.lcons[i] != cache.ucons[i]
+                push!(linineqsinds, i)
             else
                 push!(nonlininds, i)
             end
         end
-        res1 = zeros(length(cache.f.cons_expr))
+        res1 = zeros(length(cache.lcons))
         nonlincons = (res, θ) -> (cache.f.cons(res1, θ); res .= res1[nonlininds])
-        println(symboliclinexpr)
-        A = hcat(res1[lininds]...)
-        println(A)
-        b = cache.lcons[lininds]
-        println(b)
+        A₁ = J[lineqsinds, :]
+        b₁ = cache.lcons[lineqsinds]
+        A₂ = J[linineqsinds, :]
+        b₂ = cache.ucons[linineqsinds]
         function fwcons(θ, res)
             nonlincons(res, θ)
             return _loss(θ)
         end
-        (minx, minf, nf, rc, cstrv) = optfunc(fwcons, cache.u0; linear_eq = (A, b), nonlinear_ineq = length(nonlininds), kws...)
+        (minx, minf, nf, rc, cstrv) = optfunc(fwcons, cache.u0; linear_eq = (A₁, b₁), linear_ineq = (A₂, b₂), nonlinear_ineq = length(nonlininds), kws...)
     elseif cache.opt isa LINCOA
         (minx, minf, nf, rc, cstrv) = optfunc(_loss, cache.u0; kws...)
     else

diff --git a/lib/OptimizationPRIMA/test/runtests.jl b/lib/OptimizationPRIMA/test/runtests.jl
@@ -1,4 +1,4 @@
-using OptimizationPRIMA, Optimization
+using OptimizationPRIMA, Optimization, ForwardDiff, ModelingToolkit
 using Test
 
 @testset "OptimizationPRIMA.jl" begin
@@ -21,6 +21,24 @@ using Test
     function con2_c(res, x, p)
         res .= [x[1] + x[2], x[2] * sin(x[1]) - x[1]]
     end
-    optprob = OptimizationFunction(rosenbrock, cons = con2_c)
-    prob = OptimizationProblem(, x0, _p, lcons = [-Inf, -Inf], ucons = [Inf, Inf])
+    optprob = OptimizationFunction(rosenbrock, AutoForwardDiff(), cons = con2_c)
+    prob = OptimizationProblem(optprob, x0, _p, lcons = [1, -100], ucons = [1, 100])
+    sol = Optimization.solve(prob, COBYLA(), maxiters = 1000)
+    @test sol.objective < l1
+
+    function con2_c(res, x, p)
+        res .= [x[1] + x[2]]
+    end
+    optprob = OptimizationFunction(rosenbrock, AutoForwardDiff(), cons = con2_c)
+    prob = OptimizationProblem(optprob, x0, _p, lcons = [1], ucons = [1])
+    sol = Optimization.solve(prob, COBYLA(), maxiters = 1000)
+    @test sol.objective < l1
+
+    function con2_c(res, x, p)
+        res .= [x[2] * sin(x[1]) - x[1]]
+    end
+    optprob = OptimizationFunction(rosenbrock, AutoModelingToolkit(), cons = con2_c)
+    prob = OptimizationProblem(optprob, x0, _p, lcons = [10], ucons = [50])
+    sol = Optimization.solve(prob, COBYLA(), maxiters = 1000)
+    @test 10*sol.objective < l1
 end