在Julia中并行进行粒子群优化

如何解决在Julia中并行进行粒子群优化

我用Julia编写了粒子群优化算法，并尝试使用Threads加速每个总体的计算。 这是代码

function uniform(dim::Int,lb::Array{Float64,1},ub::Array{Float64,1})
    arr = rand(Float64,dim)
    @inbounds for i in 1:dim; arr[i] = arr[i] * (ub[i] - lb[i]) + lb[i] end
    return arr
end

function initialize_particles(problem,population)
    dim = problem.dim
    lb = problem.lb
    ub = problem.ub
    cost_func = problem.cost_func

    gbest_position = uniform(dim,lb,ub)
    gbest = Gbest(gbest_position,cost_func(gbest_position))

    particles = []
    for i in 1:population
        position = uniform(dim,ub)
        velocity = zeros(dim)
        cost = cost_func(position)
        best_position = copy(position)
        best_cost = copy(cost)
        push!(particles,Particle(position,velocity,cost,best_position,best_cost))

        if best_cost < gbest.cost
            gbest.position = copy(best_position)
            gbest.cost = copy(best_cost)
        end
    end
    return gbest,particles
end

function PSO(problem; max_iter=100,population=100,c1=1.4962,c2=1.4962,w=0.7298,wdamp=1.0)
    dim = problem.dim
    lb = problem.lb
    ub = problem.ub
    cost_func = problem.cost_func

    gbest,particles = initialize_particles(problem,population)

    # main loop
    for iter in 1:max_iter
        @threads for i in 1:population
            particles[i].velocity .= w .* particles[i].velocity .+
                c1 .* rand(dim) .* (particles[i].best_position .- particles[i].position) .+
                c2 .* rand(dim) .* (gbest.position .- particles[i].position)

            particles[i].position .= particles[i].position .+ particles[i].velocity
            particles[i].position .= max.(particles[i].position,lb)
            particles[i].position .= min.(particles[i].position,ub)

            particles[i].cost = cost_func(particles[i].position)

            if particles[i].cost < particles[i].best_cost
                particles[i].best_position = copy(particles[i].position)
                particles[i].best_cost = copy(particles[i].cost)
        
                if particles[i].best_cost < gbest.cost
                    gbest.position = copy(particles[i].best_position)
                    gbest.cost = copy(particles[i].best_cost)
                end
            end
        end
        w = w * wdamp
        if iter % 50 == 1
            println("Iteration " * string(iter) * ": Best Cost = " * string(gbest.cost))
            println("Best Position = " * string(gbest.position))
            println()
        end
    end
    gbest,particles
end

很明显，它在PSO的主循环中不是线程安全的。这是因为计算时需要修改gbest的位置。 我尝试使用Atomic，但是每个粒子的位置都是一个数组。所以这种方法不合适。 锁需要初始化，因此也将无法使用。 我也尝试过使用Floop，但是出了错。

那么有没有一种方法可以确保在进行并行计算时数据被锁定？

谢谢您的时间！

王

解决方法

无需锁定变量。当您进行并行计算时，您总是可以从Split-Combine策略的角度来思考（抱歉，与DataFrames计算类似，只是创造了这个术语，您也可以将其称为Map-Reduce方法）。这样做的想法是，您将计算拆分成不同的线程，然后独立进行，最后阶段，将结果合并到一个线程中。因此您的代码看起来像这样（可能存在语法错误，因为如果没有Particle和其他定义，我将无法运行代码，但是我希望它足以给出想法）

function minby(itr; by=identity,init=nothing)
     init = isnothing(init) ? pop!(itr) : init
     mapreduce(x->(by(x)=>x),(x,y)->(first(x)>first(y)) ? y : x,itr; init=by(init)=>init).second
end

function PSO(problem; max_iter=100,population=100,c1=1.4962,c2=1.4962,w=0.7298,wdamp=1.0)
    dim = problem.dim
    lb = problem.lb
    ub = problem.ub
    cost_func = problem.cost_func

    gbest,particles = initialize_particles(problem,population)

    # main loop
    for iter in 1:max_iter
        gbests = fill((gbest.cost,0),Threads.nthreads())
        @threads for i in 1:population
            particles[i].velocity .= w .* particles[i].velocity .+
                c1 .* rand(dim) .* (particles[i].best_position .- particles[i].position) .+
                c2 .* rand(dim) .* (gbest.position .- particles[i].position)

            particles[i].position .= particles[i].position .+ particles[i].velocity
            particles[i].position .= max.(particles[i].position,lb)
            particles[i].position .= min.(particles[i].position,ub)

            particles[i].cost = cost_func(particles[i].position)

            if particles[i].cost < particles[i].best_cost
                particles[i].best_position = copy(particles[i].position)
                particles[i].best_cost = copy(particles[i].cost)
        
                if particles[i].best_cost < gbests[Threads.threadid()][1]
                    gbests[Threads.threadid()] = (particles[i].best_cost,i)
                end
            end
        end
        gbest1 = minby(gbests,by = x -> x[1])
        if gbest1[2] != 0
            idx = gbest1[2]
            gbest.position = copy(particles[idx].best_position)
            gbest.cost = copy(particles[idx].best_cost)
        end
        w = w * wdamp
        if iter % 50 == 1
            println("Iteration " * string(iter) * ": Best Cost = " * string(gbest.cost))
            println("Best Position = " * string(gbest.position))
            println()
        end
    end
    gbest,particles
end

顺便说一句，您可能会发现软件包UnPack.jl非常方便。除了手动分配，您还可以

using UnPack
function PSO(problem; max_iter=100,wdamp=1.0)
    @unpack dim,lb,ub,cost_func = problem
    ....