Updating solver __solve function for MOO

lib/OptimizationBBO/src/OptimizationBBO.jl

@@ -3,6 +3,7 @@ module OptimizationBBO
 using Reexport
 import Optimization
 import BlackBoxOptim, Optimization.SciMLBase
+import Optimization.SciMLBase: MultiObjectiveOptimizationFunction
 
 abstract type BBO end
 
@@ -15,6 +16,11 @@ for j in string.(BlackBoxOptim.SingleObjectiveMethodNames)
     eval(Meta.parse("export BBO_" * j))
 end
 
+Base.@kwdef struct BBO_borg_moea <: BBO
+    method = :borg_moea 
+end
+export BBO_borg_moea
+
 function decompose_trace(opt::BlackBoxOptim.OptRunController, progress)
     if progress
         maxiters = opt.max_steps
@@ -142,17 +148,32 @@ function SciMLBase.__solve(cache::Optimization.OptimizationCache{
     maxtime = Optimization._check_and_convert_maxtime(cache.solver_args.maxtime)
 
     _loss = function (θ)
-        if cache.callback === Optimization.DEFAULT_CALLBACK &&
-           cache.data === Optimization.DEFAULT_DATA
-            return first(cache.f(θ, cache.p))
-        elseif cache.callback === Optimization.DEFAULT_CALLBACK
-            return first(cache.f(θ, cache.p, cur...))
-        elseif cache.data !== Optimization.DEFAULT_DATA
-            x = cache.f(θ, cache.p)
-            return first(x)
+        if isa(cache.f, MultiObjectiveOptimizationFunction)
+            if cache.callback === Optimization.DEFAULT_CALLBACK &&
+               cache.data === Optimization.DEFAULT_DATA
+                return cache.f(θ, cache.p)
+            elseif cache.callback === Optimization.DEFAULT_CALLBACK
+                return cache.f(θ, cache.p, cur...)
+            elseif cache.data !== Optimization.DEFAULT_DATA
+                x = cache.f(θ, cache.p)
+                return x
+            else
+                x = cache.f(θ, cache.p, cur...)
+                return first(x)
+            end
         else
-            x = cache.f(θ, cache.p, cur...)
-            return first(x)
+            if cache.callback === Optimization.DEFAULT_CALLBACK &&
+               cache.data === Optimization.DEFAULT_DATA
+                return first(cache.f(θ, cache.p))
+            elseif cache.callback === Optimization.DEFAULT_CALLBACK
+                return first(cache.f(θ, cache.p, cur...))
+            elseif cache.data !== Optimization.DEFAULT_DATA
+                x = cache.f(θ, cache.p)
+                return first(x)
+            else
+                x = cache.f(θ, cache.p, cur...)
+                return first(x)
+            end
         end
     end
 

lib/OptimizationBBO/test/runtests.jl

@@ -1,4 +1,5 @@
-using OptimizationBBO, Optimization
+using OptimizationBBO, Optimization, BlackBoxOptim
+using Optimization.SciMLBase: MultiObjectiveOptimizationFunction
 using Test
 
 @testset "OptimizationBBO.jl" begin
@@ -46,4 +47,70 @@ using Test
             progress = true,
             maxtime = 5)
     end
+
+    # Define the initial guess and bounds
+u0 = [0.25, 0.25]
+lb = [0.0, 0.0]
+ub = [2.0, 2.0]
+
+# Define the optimizer
+opt = OptimizationBBO.BBO_borg_moea()
+
+@testset "Multi-Objective Optimization Tests" begin
+
+    # Test 1: Sphere and Rastrigin Functions
+    @testset "Sphere and Rastrigin Functions" begin
+        function multi_obj_func_1(x, p)
+            f1 = sum(x .^ 2)  # Sphere function
+            f2 = sum(x .^ 2 .- 10 .* cos.(2π .* x) .+ 10)  # Rastrigin function
+            return (f1, f2)
+        end
+
+        mof_1 = MultiObjectiveOptimizationFunction(multi_obj_func_1)
+        prob_1 = Optimization.OptimizationProblem(mof_1, u0; lb=lb, ub=ub)
+        sol_1 = solve(prob_1, opt, NumDimensions=2, FitnessScheme=ParetoFitnessScheme{2}(is_minimizing=true))
+
+        @test sol_1 ≠ nothing
+        println("Solution for Sphere and Rastrigin: ", sol_1)
+        @test sol_1.objective[1] ≈ 6.9905986e-18 atol=1e-3
+        @test sol_1.objective[2] ≈ 1.7763568e-15 atol=1e-3
+    end
+
+    # Test 2: Rosenbrock and Ackley Functions
+    @testset "Rosenbrock and Ackley Functions" begin
+        function multi_obj_func_2(x, p)
+            f1 = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2  # Rosenbrock function
+            f2 = -20.0 * exp(-0.2 * sqrt(0.5 * (x[1]^2 + x[2]^2))) - exp(0.5 * (cos(2π * x[1]) + cos(2π * x[2]))) + exp(1) + 20.0  # Ackley function
+            return (f1, f2)
+        end
+
+        mof_2 = MultiObjectiveOptimizationFunction(multi_obj_func_2)
+        prob_2 = Optimization.OptimizationProblem(mof_2, u0; lb=lb, ub=ub)
+        sol_2 = solve(prob_2, opt, NumDimensions=2, FitnessScheme=ParetoFitnessScheme{2}(is_minimizing=true))
+
+        @test sol_2 ≠ nothing
+        println("Solution for Rosenbrock and Ackley: ", sol_2)
+        @test sol_2.objective[1] ≈ 0.97438 atol=1e-3
+        @test sol_2.objective[2] ≈ 0.04088 atol=1e-3
+    end
+
+    # Test 3: ZDT1 Function
+    @testset "ZDT1 Function" begin
+        function multi_obj_func_3(x, p)
+            f1 = x[1]
+            g = 1 + 9 * sum(x[2:end]) / (length(x) - 1)
+            f2 = g * (1 - sqrt(f1 / g))
+            return (f1, f2)
+        end
+
+        mof_3 = MultiObjectiveOptimizationFunction(multi_obj_func_3)
+        prob_3 = Optimization.OptimizationProblem(mof_3, u0; lb=lb, ub=ub)
+        sol_3 = solve(prob_3, opt, NumDimensions=2, FitnessScheme=ParetoFitnessScheme{2}(is_minimizing=true))
+
+        @test sol_3 ≠ nothing
+        println("Solution for ZDT1: ", sol_3)
+        @test sol_3.objective[1] ≈ 0.273445 atol=1e-3
+        @test sol_3.objective[2] ≈ 0.477079 atol=1e-3
+    end
+end
 end

lib/OptimizationEvolutionary/src/OptimizationEvolutionary.jl

@@ -8,6 +8,7 @@ SciMLBase.allowsbounds(opt::Evolutionary.AbstractOptimizer) = true
 SciMLBase.allowsconstraints(opt::Evolutionary.AbstractOptimizer) = true
 SciMLBase.supports_opt_cache_interface(opt::Evolutionary.AbstractOptimizer) = true
 SciMLBase.requiresgradient(opt::Evolutionary.AbstractOptimizer) = false
+SciMLBase.requiresgradient(opt::Evolutionary.NSGA2) = false
 SciMLBase.requireshessian(opt::Evolutionary.AbstractOptimizer) = false
 SciMLBase.requiresconsjac(opt::Evolutionary.AbstractOptimizer) = false
 SciMLBase.requiresconshess(opt::Evolutionary.AbstractOptimizer) = false
@@ -125,8 +126,13 @@ function SciMLBase.__solve(cache::OptimizationCache{
     f = cache.f
 
     _loss = function (θ)
-        x = f(θ, cache.p, cur...)
-        return first(x)
+        if isa(f, MultiObjectiveOptimizationFunction)
+            x = f(θ, cache.p, cur...)
+            return x
+        else
+            x = f(θ, cache.p, cur...)
+            return first(x)
+        end
     end
 
     opt_args = __map_optimizer_args(cache, cache.opt; callback = _cb, cache.solver_args...,
@@ -139,9 +145,17 @@ function SciMLBase.__solve(cache::OptimizationCache{
             c = x -> (res = zeros(length(cache.lcons)); f.cons(res, x); res)
             cons = WorstFitnessConstraints(Float64[], Float64[], cache.lcons, cache.ucons,
                 c)
-            opt_res = Evolutionary.optimize(_loss, cons, cache.u0, cache.opt, opt_args)
+            if isa(f, MultiObjectiveOptimizationFunction)
+                opt_res = Evolutionary.optimize(_loss, _loss(cache.u0), cons, cache.u0, cache.opt, opt_args)
+            else
+                opt_res = Evolutionary.optimize(_loss, cons, cache.u0, cache.opt, opt_args)
+            end
         else
-            opt_res = Evolutionary.optimize(_loss, cache.u0, cache.opt, opt_args)
+            if isa(f, MultiObjectiveOptimizationFunction)
+                opt_res = Evolutionary.optimize(_loss, _loss(cache.u0), cache.u0, cache.opt, opt_args)
+            else
+                opt_res = Evolutionary.optimize(_loss, cache.u0, cache.opt, opt_args)
+            end
         end
     else
         if !isnothing(f.cons)
@@ -150,17 +164,30 @@ function SciMLBase.__solve(cache::OptimizationCache{
         else
             cons = BoxConstraints(cache.lb, cache.ub)
         end
-        opt_res = Evolutionary.optimize(_loss, cons, cache.u0, cache.opt, opt_args)
+        if isa(f, MultiObjectiveOptimizationFunction)
+                opt_res = Evolutionary.optimize(_loss, _loss(cache.u0), cons, cache.u0, cache.opt, opt_args)
+        else
+                opt_res = Evolutionary.optimize(_loss, cons, cache.u0, cache.opt, opt_args)
+        end
     end
     t1 = time()
     opt_ret = Symbol(Evolutionary.converged(opt_res))
     stats = Optimization.OptimizationStats(; iterations = opt_res.iterations,
         time = t1 - t0, fevals = opt_res.f_calls)
-    SciMLBase.build_solution(cache, cache.opt,
+    if !isa(f, MultiObjectiveOptimizationFunction)
+        SciMLBase.build_solution(cache, cache.opt,
         Evolutionary.minimizer(opt_res),
         Evolutionary.minimum(opt_res); original = opt_res,
         retcode = opt_ret,
         stats = stats)
+    else
+        ans = Evolutionary.minimizer(opt_res)
+        SciMLBase.build_solution(cache, cache.opt,
+        ans,
+        _loss(ans[1]); original = opt_res,
+        retcode = opt_ret,
+        stats = stats)
+    end
 end
 
 end

lib/OptimizationEvolutionary/test/runtests.jl

@@ -1,4 +1,5 @@
 using OptimizationEvolutionary, Optimization, Random
+using Optimization.SciMLBase: MultiObjectiveOptimizationFunction
 using Test
 
 Random.seed!(1234)
@@ -56,4 +57,111 @@ Random.seed!(1234)
     # Make sure that both the user's trace record value, as well as `curr_u` are stored in the trace.
     @test haskey(sol.original.trace[end].metadata, "TESTVAL") &&
           haskey(sol.original.trace[end].metadata, "curr_u")
+
+    # Test Suite for Different Multi-Objective Functions
+    function test_multi_objective(func, initial_guess)
+        # Define the gradient function using ForwardDiff
+        function gradient_multi_objective(x, p=nothing)
+            ForwardDiff.jacobian(func, x)
+        end
+
+        # Create an instance of MultiObjectiveOptimizationFunction
+        obj_func = MultiObjectiveOptimizationFunction(func, jac=gradient_multi_objective)
+
+        # Set up the evolutionary algorithm (e.g., NSGA2)
+        algorithm = OptimizationEvolutionary.NSGA2()
+
+        # Define the optimization problem
+        problem = OptimizationProblem(obj_func, initial_guess)
+
+        # Solve the optimization problem
+        result = solve(problem, algorithm)
+
+        return result
+    end
+
+    @testset "Multi-Objective Optimization Tests" begin
+
+        # Test 1: Sphere and Rastrigin Functions
+        @testset "Sphere and Rastrigin Functions" begin
+            function multi_objective_1(x, p=nothing)::Vector{Float64}
+                f1 = sum(x .^ 2)  # Sphere function
+                f2 = sum(x .^ 2 .- 10 .* cos.(2π .* x) .+ 10)  # Rastrigin function
+                return [f1, f2]
+            end
+            result = test_multi_objective(multi_objective_1, [0.0, 1.0])
+            @test result ≠ nothing
+            println("Solution for Sphere and Rastrigin: ", result)
+            @test result.u[1][1] ≈ 7.88866e-5 atol=1e-3
+            @test result.u[1][2] ≈ 4.96471e-5 atol=1e-3
+            @test result.objective[1] ≈ 8.6879e-9 atol=1e-3
+            @test result.objective[2] ≈ 1.48875349381683e-6 atol=1e-3
+        end
+
+        # Test 2: Rosenbrock and Ackley Functions
+        @testset "Rosenbrock and Ackley Functions" begin
+            function multi_objective_2(x, p=nothing)::Vector{Float64}
+                f1 = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2  # Rosenbrock function
+                f2 = -20.0 * exp(-0.2 * sqrt(0.5 * (x[1]^2 + x[2]^2))) - exp(0.5 * (cos(2π * x[1]) + cos(2π * x[2]))) + exp(1) + 20.0  # Ackley function
+                return [f1, f2]
+            end
+            result = test_multi_objective(multi_objective_2, [0.1, 1.0])
+            @test result ≠ nothing
+            println("Solution for Rosenbrock and Ackley: ", result)
+            @test result.u[1][1] ≈ 0.003993274873103834 atol=1e-3
+            @test result.u[1][2] ≈ 0.001433311246712721 atol=1e-3
+            @test result.objective[1] ≈ 0.9922302888530358 atol=1e-3
+            @test result.objective[2] ≈ 0.012479470703588902 atol=1e-3
+        end
+
+        # Test 3: ZDT1 Function
+        @testset "ZDT1 Function" begin
+            function multi_objective_3(x, p=nothing)::Vector{Float64}
+                f1 = x[1]
+                g = 1 + 9 * sum(x[2:end]) / (length(x) - 1)
+                sqrt_arg = f1 / g
+                f2 = g * (1 - (sqrt_arg >= 0 ? sqrt(sqrt_arg) : NaN))
+                return [f1, f2]
+            end
+            result = test_multi_objective(multi_objective_3, [0.25, 1.5])
+            @test result ≠ nothing
+            println("Solution for ZDT1: ", result)
+            @test result.u[1][1] ≈ -0.365434 atol=1e-3
+            @test result.u[1][2] ≈ 1.22128 atol=1e-3
+            @test result.objective[1] ≈ -0.365434 atol=1e-3
+            @test isnan(result.objective[2])
+        end
+
+        # Test 4: DTLZ2 Function
+        @testset "DTLZ2 Function" begin
+            function multi_objective_4(x, p=nothing)::Vector{Float64}
+                f1 = (1 + sum(x[2:end] .^ 2)) * cos(x[1] * π / 2)
+                f2 = (1 + sum(x[2:end] .^ 2)) * sin(x[1] * π / 2)
+                return [f1, f2]
+            end
+            result = test_multi_objective(multi_objective_4, [0.25, 0.75])
+            @test result ≠ nothing
+            println("Solution for DTLZ2: ", result)
+            @test result.u[1][1] ≈ 0.899183 atol=1e-3
+            @test result.u[2][1] ≈ 0.713992 atol=1e-3
+            @test result.objective[1] ≈ 0.1599915 atol=1e-3
+            @test result.objective[2] ≈ 1.001824893932647 atol=1e-3
+        end
+
+        # Test 5: Schaffer Function N.2
+        @testset "Schaffer Function N.2" begin
+            function multi_objective_5(x, p=nothing)::Vector{Float64}
+                f1 = x[1]^2
+                f2 = (x[1] - 2)^2
+                return [f1, f2]
+            end
+            result = test_multi_objective(multi_objective_5, [1.0])
+            @test result ≠ nothing
+            println("Solution for Schaffer N.2: ", result)
+            @test result.u[19][1] ≈ 0.252635 atol=1e-3
+            @test result.u[9][1] ≈ 1.0 atol=1e-3
+            @test result.objective[1] ≈ 1.0 atol=1e-3
+            @test result.objective[2] ≈ 1.0 atol=1e-3
+        end
+    end
 end