Tag: Optimization

I posted this question in Stackoverflow, but some people suggested me to post it here.

I rewrite the SCE-UA optimization algorithm in Julia.

However, when I run my script in Julia. The optimized value is quite different from the results in MATLAB. The global best cost in MATLAB is 2.4598e-63 while the global best cost in Julia is 8.264629809290885e-8. This means something went wrong. The objective function is z=sum(x.^2) for both scripts, and the search range and parameters for the algorithm were also the same. However, I checked, again and again, to make sure the algorithms are the same, but I could not find out why the results were quite different.

Could some please check if those two scripts were different? Please give me any suggestions, thanks a lot.

The SCE-UA code in MATLAB is written by Yarpiz, which can be found in https://github.com/smkalami/ypea110-shuffled-complex-evolution/archive/master.zip

The code in Julia is as follows

using UnPack
using Plots
using StatsBase

mutable struct Pop
    position::Vector{Float64}
    cost::Float64
end

# CCE parameters
mutable struct CCE_params
    q::Int64
    alpha::Int64
    beta::Int64
    lb::Vector{Float64}
    ub::Vector{Float64}
    obj_func
end

# SCE parameters
mutable struct SCE_params
    max_iter::Int64
    n_complex::Int64
    n_complex_pop::Int64
    dim::Int64
    lb::Vector{Float64}
    ub::Vector{Float64}
    obj_func
end

import Base.isless
isless(a::Pop, b::Pop) = isless(a.cost, b.cost)

function cost_func(x)
    return sum(x.^2)
end

function uniform_rand(lb::Array{Float64, 1}, ub::Array{Float64, 1})
    dim = length(lb)
    arr = rand(dim) .* (ub .- lb) .+ lb
    return arr
end
# SCE parameters
max_iter = 500
n_complex = 5
n_complex_pop = 10
dim = 10
lb = ones(dim) * -10
ub = ones(dim) * 10
obj_func = cost_func
n_complex_pop = max(n_complex_pop, dim+1) # Nelder-Mead Standard
sce_params = SCE_params(max_iter, n_complex, n_complex_pop, dim, lb, ub, obj_func)

# CCE parameters
cce_q = max(round(Int64, 0.5*n_complex_pop), 2)
cce_alpha = 3
cce_beta = 5

cce_params = CCE_params(cce_q, cce_alpha, cce_beta, lb, ub, obj_func)

function SCE(sce_params, cce_params)
    @unpack max_iter, n_complex, n_complex_pop, dim, lb, ub, obj_func = sce_params

    n_pop = n_complex * n_complex_pop
    I = reshape(1:n_pop, n_complex, :)

    # Step 1. Generate rand_sample
    best_costs = Vector{Float64}(undef, max_iter)

    pops = ()
    for i in 1:n_pop
        pop_position = uniform_rand(lb, ub)
        pop_cost = obj_func(pop_position)
        pop = Pop(pop_position, pop_cost)
        push!(pops, pop)
    end
    complex = Array{Pop}(undef, n_complex_pop, n_complex)

    # Step 2. Rank Points
    sort!(pops)
    best_pop = pops(1)

    # Main loop
    for iter in 1:max_iter

        # Step 3. Partion into complexes
        for j in 1:n_complex
            complex(:,j) = deepcopy(pops(I(j,:)))
            # Step 4. Evolve complex, run CCE
            complex(:,j) = CCE(complex(:,j), cce_params)
            pops(I(j,:)) = deepcopy(complex(:,j))
        end
        # Step 5. Shuffle Complexes

        sort!(pops)
        best_pop = pops(1)

        best_costs(iter) = best_pop.cost

        # Show Iteration Information
        println("Iter = ", iter)
        println("The Best Cost is: ", best_costs(iter))
    end
    best_costs
end

function rand_sample(P, q)
    L = Vector{Int64}(undef, q)
    for i in 1:q
        L(i) = sample(1:length(P), Weights(P))
        # L(i) = sample(1:sizeof(P), weights(P), 1, replace=true)
    end
    return L
end

function not_in_search_space(position, lb, ub)
    return any(position .<= lb) || any(position .>= ub)
end

function CCE(complex_pops, cce_params)
    # Step 1. Initialize
    @unpack q, alpha, beta, lb, ub, obj_func = cce_params
    n_pop = length(complex_pops)

    # Step 2. Assign weights
    P = (2*(n_pop+1-i) / (n_pop*(n_pop+1)) for i in 1:n_pop)

    # Calculate Population Range (Smallest Hypercube)
    new_lb = complex_pops(1).position
    new_ub = complex_pops(1).position
    for i in 2:n_pop
        new_lb = min.(new_lb, complex_pops(i).position)
        new_ub = max.(new_ub, complex_pops(i).position)
    end

    # CCE main loop
    for it in 1:beta
        # Step 3. Select parents
        L = rand_sample(P, q)
        B = complex_pops(L)

        # Step 4. Generate Offspring
        for k in 1:alpha
            # a) Sort population
            sorted_indexs = sortperm(B)
            sort!(B)
            L(:) = L(sorted_indexs)

            # Calculate the centroid
            g = zeros(length(lb))
            for i in 1:q-1
                g .= g .+ B(i).position
            end
            g .= g ./ (q-1)

            # b) Reflection step
            reflection = deepcopy(B(end))
            reflection.position = 2 .* g .- B(end).position # newly generated point using reflection
            if not_in_search_space(reflection.position, lb, ub)
                reflection.position = uniform_rand(new_lb, new_ub)
            end
            reflection.cost = obj_func(reflection.position)

            if reflection.cost < B(end).cost
                B(end) = deepcopy(reflection)
            else # Contraction
                contraction = deepcopy(B(end))
                contraction.position = (g .+ B(end).position) ./ 2
                contraction.cost = obj_func(contraction.position)

                if contraction.cost < B(end).cost
                    B(end) = deepcopy(contraction)
                else
                    B(end).position = uniform_rand(new_lb, new_ub)
                    B(end).cost = obj_func(B(end).position)
                end
            end

        end

        complex_pops(L) = B
    end
    return complex_pops
end
best_costs = SCE(sce_params, cce_params)

plot(best_costs, yaxis=:log, label = "cost")

# savefig("Julia.png")

Both scripts could run normally with just one click.

The result of SCE-UA algorithm from MATLAB is as follows

The result from Julia is as follows

Obviously, the convergence curves are quite different. The result from Julia is not accurate enough, but I don’t know why.