Julia and the swarm of particles

  • Tutorial


We continue to study methods of multidimensional optimization, and the next in line is the particle swarm method searching for a global minimum.


Theory



The algorithm is pretty simple:



The position of each particle at a certain moment is calculated by the formula:


$ V_ {i, t + 1} = A_cV_ {i, t} + C_pr_p (p_i-x_ {i, t}) + C_gr_g (g_i-x_ {i, t}) $


$ x_ {i, t + 1} = x_ {i, t} + V_ {i, t + 1} $


Where $ p_i $ - coordinate of the best solution of a particular particle, $ g_i $ - coordinate of the best solution for all particles for this era, $ C_p $ and $ C_g $ - weight coefficients (selected for a specific model), $ A_c $ - the coefficient of inertia, it can be made dependent on the number of the epoch, then the speed of the particles will vary smoothly.


Test functions


Since it is very pleasant to look at the work of the method, we will get more test functions:


Tucked away
parabol(x) = sum(u->u*u, x) 
# f(0,0) = 0, x_i ∈ [-10,10]
shvefel(x) = sum(u-> -u*sin(sqrt(abs(u))), x) 
# f(420.9687,420.9687) = -819?, x_i ∈ [-500,500]
rastrigin(x) = 10*length(x) + sum(u->u*u-10*cos(2*pi*u), x) 
# f(0,0) = 0, x_i ∈ [-5,5]
ekly(x) = -20exp(-0.2sqrt(0.5(x[1]*x[1]+x[2]*x[2]))) - exp(0.5(cospi(2x[1])+cospi(2x[2]))) + 20 + ℯ
# f(0,0) = 0, x_i ∈ [-5,5]
rosenbrok(x) = 100(x[2]-x[1]*x[1])^2 + (x[1]-1)^2# f(0,0) = 0, x_i ∈ [-5,5]
bill(x) = (1.5-x[1]+x[1]*x[2])^2 + (2.25-x[1]+x[1]*x[2]*x[2])^2 + (2.625-x[1]+x[1]*x[2]^3)^2# f(3,0.5) = 0, x_i ∈ [-5,5]
boot(x) = (x[1]+2x[2]-7)^2 + (2x[1]+x[2]-5)^2# f(1,3) = 0, x_i ∈ [-10,10]
bukin6(x) = 100sqrt(abs(x[2]-0.01x[1]*x[1])) + 0.01abs(x[1]+10) 
# f(-10,1) = 0, x_i ∈ [-15,-5; -3,3]
levy13(x) = sinpi(3x[1])^2 + (1+sinpi(3x[2])^2)*(x[1]-1)^2 + (1+sinpi(2x[2])^2)*(x[2]-1)^2# f(1,1) = 0, x_i ∈ [-10,10]
himmelblau(x) = (x[1]^2+x[2]-11)^2 + (x[1]+x[2]^2-7)^2# f(3,2)... = 0, x_i ∈ [-5,5]
camel3humped(x) = 2x[1]^2 - 1.05x[1]^4 + x[1]^6 /6 + x[1]*x[2] + x[2]^2# f(0,0) = 0, x_i ∈ [-5,5]
izom(x) = -cos(x[1])*cos(x[2])*exp(-( (x[1]-pi)^2 + (x[2]-pi)^2 )) 
# f(π,π) = -1, x_i ∈ [-100,100]
holdertable(x) = -abs(sin(x[1])*cos(x[2])exp(abs( 1-sqrt(x[1]^2+x[2]^2)/pi ))) 
# f(±8.05502,±9.66459) = -19.2085, x_i ∈ [-10,10]
shaffer4(x) = 0.5 + (cos(sin(abs(x[1]^2-x[2]^2)))^2-0.5) / (1+0.001(x[1]^2+x[2]^2))^2# f(0,1.25313) = 0.292579, x_i ∈ [-100,100]

And, actually, MUF itself:


function mdpso(;
       nparts = 50,
       ndimes = 2,
       ages = 50, # количество эпох
       lover = [-10-10],
       upper = [1010],
       C1 = [1.91.9], # весовые коэф-ты
       C2 = [1.81.8],
       Ac = [0.10.1],
       )
    minind = 0
    V = zeros(nparts,ndimes) # матрица нулей n на n
    X = zeros(nparts,ndimes)
    funmin = -Inf
    Fmin = Inf
    Fbest = fill(Fmin, nparts)
    funx = zeros(nparts)
    xmem = zeros(nparts,ndimes)
    xbest = zeros(ndimes) # лучшая координата# частицы разбрасываются по исследуемой областиfor i in1:nparts, j in1:ndimes
        X[i,j] = randomer(lover[j], upper[j])
    end
    for i in1:ages
        for j in1:nparts
            funx[j] = fun(X[j,:])
            if funx[j] < Fbest[j]
                Fbest[j] = funx[j]
                xmem[j,:] = X[j,:]
            end
        end
         # отрисовывает частицы на каждой эпохе
        ploter(lover, upper, X, funx, i);
        funmin = minimum(funx)
        minind = argmin(funx)
        if funmin < Fmin
            Fmin = funmin
            xbest[:] = X[minind,:]
        end
        for j in1:nparts, k in1:ndimes
            R1 = rand()
            R2 = rand()
            V[j,k] = Ac[k]*V[j,k] + C1[k]*R1*(xmem[j,k] - X[j,k]) + 
                C2[k]*R2*(xbest[k] - X[j,k])
            X[j,k] += V[j,k]
        end
        println("Age № $i\n xbest:\n $(xbest[1]) $(xbest[2])")
        println("Fmin: $Fmin\n")
    end
    f = open("$fun.txt","w") # выводим параметры модели в файл
    write(f,"C1 = $C1, C2 = $C2, Ac = $Ac, lower = $lover, upper = $upper, ages = $ages, parts = $nparts")
    close(f)
end

With drawings, though, the calculation is performed longer, but more beautiful:


using Plots
pyplot()
function ploter(l, u, xy, z, n_age )
    contour(Xs, Ys, Zs, fill = true);
    # легенду не показывать (для каждой блин частицы)# задать границы рисунка, чтоб не дергалось кода частица убежала
    scatter!(xy[:,1], xy[:,2], legend = false, xaxis=( (l[1], u[1])), yaxis=( (l[2], u[2])) );
    #savefig("$fun $n_age.png") # создаёт кадры эпох в папке с проектом
end
fun = bill
low = [-4-4]
up = [44]
    Xs = range(low[1], stop = up[1], length = 80)
    Ys = range(low[2], stop = up[2], length = 80)
    Zs = [ fun([x y]) for y in Ys, x in Xs ]
    surface(Xs, Ys, Zs)
    xaxis!( (low[1], up[1]), low[1]:(up[1]-low[1])/5:up[1] )
    yaxis!( (low[2], up[2]), low[2]:(up[2]-low[2])/5:up[2] )


mdpso(C1 = [1.21.2], C2 = [1.11.1], Ac = [0.080.08], lower = [-4-4], upper = [44], ages = 30)


fun = ekly
mdpso(C1 = [1.71.7], C2 = [1.71.7], Ac = [0.070.07], lower = [-5-5], upper = [55], ages = 15)



fun = himmelblau
mdpso(C1 = [1.11.1], C2 = [1.01.0], Ac = [0.090.09], lower = [-5-5], upper = [55], ages = 20, parts = 50)



fun = holdertable
mdpso(C1 = [1.11.1], C2 = [1.01.0], Ac = [0.090.09], lower = [-10-10], upper = [1010], ages = 20, parts = 50)



fun = levy13
mdpso(C1 = [1.11.1], C2 = [1.01.0], Ac = [0.090.09], lower = [-10-10], upper = [1010], ages = 20, parts = 50)



fun = shaffer4
mdpso(C1 = [1.11.1], C2 = [1.01.0], Ac = [0.090.09], lower = [-100-100], upper = [100100], ages = 20, parts = 50)



What is really depressing is the fussing with the parameters and the element of chance: if not a single particle flew past the global minimum, then the whole thing could fall into a local one:


fun = rastrigin
mdpso(mdpso(C1 = [0.10.1], C2 = [11], Ac = [0.080.08], lover = low, upper = up, ages = 30))



And old Not really Good Rosenbrock still does not allow himself to be comprehended:


fun = rosenbrok
mdpso(C1 = [1.71.7], C2 = [1.51.5], Ac = [0.150.15], lover = low, upper = up, ages = 20, nparts = 50)
...
Age № 20
 xbest:
 0.377964213418868660.12799160066705667
Fmin: 0.409026370833564



But as I said earlier, it is possible to use the MF to find a good approximation to the global minimum, and then specify, for example, Nelder-Mead:


Simplex method
vecl(x) = sqrt( sum(u -> u*u, x) )
function sortcoord(Mx)
    N = size(Mx,2)
    f = [fun(Mx[:,i]) for i in1:N] # значение функции в вершинах
    Mx[:, sortperm(f)]
end
function normx(Mx)
    m = size(Mx,2)
    D = zeros(m-1,m)
    for i = 1:m, j = i+1:m
        D[i,j] = vecl(Mx[:,i] - Mx[:,j]) # считает длину разности столбцов
    end
    D
    sqrt(maximum(D))
end
function ofNelderMid(; ndimes = 2, ε = 1e-4, fit = [.1, .1], low = [-1-1], up = [11])
    k = 0
    N = ndimes
    dz = zeros(N, N+1)
    Xx = zeros(N, N+1)
    for i = 1:N+1
        Xx[:,i] = fit
    end
    for i = 1:N
        dz[i,i] = 0.5*vecl(fit)
    end
    Xx += dz
    p = normx(Xx)
    while p > ε
        k += 1
        Xx = sortcoord(Xx)
        Xo = [ sum(Xx[i,1:N])/N for i = 1:N ] # среднее эл-тов i-й строки
        Ro = 2Xo - Xx[:,N+1]
        FR = fun(Ro)
        if FR > fun(Xx[:,N+1])
            for i = 2:N+1
                Xx[:,i] = 0.5(Xx[:,1] + Xx[:,i])
            end     
        elseif FR < fun(Xx[:,1])
                Eo = Xo + 2(Xo - Xx[:,N+1])
                if FR > fun(Eo)
                    Xx[:,N+1] = Eo
                else
                    Xx[:,N+1] = Ro
                end
            elseif FR <= fun(Xx[:,N])
                    Xx[:,N+1] = Ro
                else
                    Co = Xo + 0.5(Xo - Xx[:,N+1])
                    if FR > fun(Co)
                        Xx[:,N+1] = Co
                    else
                        Xx[:,N+1] = Ro
                    end
                end
            end
        end
        println(k, " ", p, " ", Xx[:,1])
        p = normx(Xx)
    end #while 
    fit = Xx[:,1]
end

ofNelderMid(fit = [0.377960.127992])
...
920.00022610400555036366 [1.0, 1.0]
930.00015987967588703512 [1.0, 1.0]
940.00011305200343052599 [1.0, 1.0]
2-element Array{Float64,1}:
 0.99999999966459730.9999999995466575

For a classic SMF, not everything is clear with the stopping criterion: the best point can hold a position for several epochs, the distances between some particles can also not change for a long time. Therefore, a limit on the number of epochs is used. In order to increase the chances of finding a global minimum, it is necessary to increase the number of particles and epochs, which is very expensive in terms of memory and, especially, time.


Morality


  • If you do not want or can not use complex and modern methods, you can use the composition of methods easier
  • Very often there is a narrow method for the task (for example, for gully functions)
  • It is advisable to have on hand several different MOs for consistent comparison.
  • It is not always possible to stick your task into the method and immediately get the right answer - you should spend time researching, varying the parameters, if you can’t see the relief, you should at least print out intermediate calculations to track convergence.

That's all for today, thank you for your attention and wish you all a good optimization!


Also popular now: