Juliaで粒子群最適化 - Yoshihiro Ohta’s blog

最近，Juliaを少しずつ学んでいます．勉強がてら，粒子群最適化(Particle Swarm Optimization, PSO)のアルゴリズムを実装し，Rosenbrock関数の最適化を行ってみました．

0. 開発環境

Windows 10 20H2
Julia 1.5.3
Visual Studio Code + Julia extension

1. 目標とコンセプト

今回の目標は，単目的の粒子群最適化アルゴリズムを実装し，Rosenbrock関数を最適化することです．ただし，将来的に多目的問題・制約付き問題に拡張可能とするための考慮をします．
以前もブログに記載しましたが，私は普段Javaの多目的最適化ライブラリjMetalを用いており，多目的に拡張する箇所などはjMetalを参考にします．
ただし，JuliaはJavaのようなクラス構造を持つ言語ではなく，型の定義と多重ディスパッチによって振る舞いを変更する特徴を持つ言語なので，そこは言語に合わせた実装とするよう努力します．

2. 問題定義

最適化問題として，今回はRosenbrock関数を対象とします．

function rosenbrock(x)
    [sum([ 100(x[i + 1] - x[i]^2)^2 + (x[i] - 1)^2 for i in 1:length(x) - 1])]
end

将来的に多目的最適化でも使うため，出力はFloat64ではなくVector{Float64}としておきます．
それから，以下の工夫をするコードを追加します．

rosenbrock以外の問題にも対応できるよう，evaluateという別の関数を経由して評価関数を実行する
rosenbrockは連続値最適化問題だが，それ以外の種類の問題ではevaluate関数の動作を変更できるよう，問題を複合型とする
複合型に，変数や目的関数の上下限値，変数の数などのパラメータをもたせる
変数・目的関数値を上下限値を用いて正規化して最適化する

# 抽象型
abstract type AbstractProblem{T <: Number} end
# 具体型
mutable struct DoubleProblem{T} <: AbstractProblem{T}
    objfunc
    numobjs
    numvars
    confunc
    numcons
    limit::Matrix{Float64}
    bound::Matrix{T}
end
# コンストラクタ
function DoubleProblem(objfunc; numvars=1, numobjs=1, confunc=NaN, numcons=0, limit=fill(NaN, 1, 1), bound=repeat([0.0 1.0], numvars))
    DoubleProblem(objfunc, numobjs, numvars, confunc, numcons, limit, bound)
end

# 評価関数
norm(value, min, max) = (value - min) ./ (max - min) # 正規化
denorm(value, min, max) = value .* (max - min) + min # 非正規化
function evaluate!(problem::DoubleProblem, solutions::Vector{<:AbstractSolution})
    for s in solutions
        vars = denorm(s.variables, problem.bound[:,1], problem.bound[:,2])
        objs = problem.objfunc(vars)
        if (isnan(problem.limit[1,1]))
            s.objectives = objs
        else
            s.objectives = norm(objs, problem.limit[:,1], problem.limit[:,2])
        end
    end
end

これで最適化問題が定義できるようになりました．先程のrosenbrock()を使って，以下のように最適化問題を定義します．

numvars = 3
bound = [[-5.0 5.0];[-5.0 5.0];[-5.0 5.0]]
problem = DoubleProblem(rosenbrock, numvars=numvars, bound=bound)

3. 粒子群最適化アルゴリズム

PSOは以下の流れで最適化を行います．

初期化：複数の粒子を作り，粒子の位置=変数を初期化する
評価：粒子の位置=変数から目的関数値を計算する
ベスト粒子の更新：global best(粒子群のうち最も良い評価の位置, gbest)と，personal best(その粒子がこれまで移動した中で最も良い評価の位置, pbest)を目的関数値を用いて更新する
飛翔：粒子の速度と位置を，gbestとpbestを用いて更新する
突然変異：粒子位置に突然変異を加える
所定の繰り返し回数に到達するまで，2に戻る．

これを書き下したのが以下です．

function pso(problem::DoubleProblem; iter, nump=problem.numvars*10, c1=NaN, c2=NaN, w=NaN, vecr=false)
    # 粒子群の初期化
    swarm = [Particle(problem.numvars) for i = 1:nump]
    gbest = Particle(problem.numvars)
    evaluate!(problem, swarm)
    updategbest!(swarm, gbest)
    updatepbest!(swarm)

    # 最適化の実行
    while (iter > 0)
        updateswarm!(swarm, gbest, c1, c2, w, vecr)
        perturbation!(swarm)
        evaluate!(problem, swarm)
        updategbest!(swarm, gbest)
        updatepbest!(swarm)
        iter -= 1
    end
    return gbest
end
# 粒子の初期化
function Particle(numvars, numobjs=1, numcons=0)
    objectives = ones(numobjs) .* Inf
    variables  = rand(numvars)
    velocities = zeros(numvars)
    pbestvars  = copy(variables)
    pbestobjs = ones(numobjs) .* Inf
    constraints = zeros(numcons)
    return Particle(objectives, variables, velocities, pbestvars, pbestobjs, constraints)
end
# 粒子速度・位置の更新
function updateswarm!(swarm::Vector{Particle{Float64}}, gbest::Particle, c1::Float64=NaN, c2::Float64=NaN, w::Float64=NaN, vecr::Bool=false)
    for p in swarm
        if(vecr)
            r1 = rand(length(p.variables))
            r2 = rand(length(p.variables))
        else
            r1 = rand()
            r2 = rand()
        end
        if(isnan(c1)) c1 = rand() * 0.5 + 1.5 end
        if(isnan(c2)) c2 = rand() * 0.5 + 1.5 end
        if(isnan(w))  w = rand() * 0.4 + 0.1 end
        
        # 位置と速度の更新
        p.velocities = w * p.velocities + c1 * r1 * (p.pbestvars - p.variables) + c2 * r2 * (gbest.variables - p.variables)
        p.variables += p.velocities

        # 境界チェック
        for i = 1:length(p.variables)
            if (p.variables[i] > 1.0)
                p.variables[i] = 1.0
                p.velocities[i] *= -1.0
            elseif (p.variables[i] < 0.0)
                p.variables[i] = 0.0
                p.velocities[i] *= -1.0
            end
        end
    end
end

# 突然変異(未実装)
function perturbation!(swarm::Vector{Particle{Float64}})

end

# pbestの更新
function updatepbest!(swarm::Vector{Particle{Float64}})
    for p in swarm
        if (p.objectives < p.pbestobjs)
            p.pbestvars = copy(p.variables)
            p.pbestobjs = copy(p.objectives)
        end
    end
end

# gbestの更新
function updategbest!(swarm::Vector{Particle{Float64}}, gbest::Particle{Float64})
    for p in swarm
        if (p.objectives < gbest.objectives)
            gbest.variables = copy(p.variables)
            gbest.objectives = copy(p.objectives)
        end
    end
end

ここで，粒子は，他の問題と共通の解の抽象型(AbstractSolution)の具体型としています．粒子は，上述のように位置・速度・評価値と，pbestの位置・評価値を持つ必要があるので，これらを型のメンバとしています．

# 抽象型
abstract type AbstractSolution{T <: Number} end
# 具体型
mutable struct Particle{T} <: AbstractSolution{T}
    objectives::Vector{Float64} # 評価値
    variables::Vector{T} # 位置
    velocities::Vector{T} # 速度
    pbestvars::Vector{T} # pbestの位置
    pbestobjs::Vector{Float64} # pbestの評価値
    constraints::Vector{T} # 制約違反量，未使用．
end

PSOの実行の際には，以下のパラメータを必要とします．
* 粒子の数: nump * 繰り返し回数: iter

また，粒子の位置・速度の更新式は次式の様になっているのですが，ここでもパラメータを指定する必要があります．

$v_{t+1} = wv_t +c_1 r_1(x_{pbest}-x_t)+c_2 r_2 (x_{gbest}-x_t)$
$x_{t+1} = x_t + v_{t+1}$

重みw
重み(加速係数)c1, c2

さらに，乱数r1, r2を，1つの解で1つの値とするか，変数ごとに異なる値とする(ベクトル化する)かでも探索の傾向が異なってくるので，これもパラメータとします．

r1, r2をベクトル化するフラグ: vecr

これらのパラメータは，指定する必要がなければデフォルトのパラメータで動くようにします．粒子の数numpは，変数の数×10としています．また，重みw, c1, c2は所定の範囲で乱数としています．今回は，OMOPSOで用いられる乱数範囲と同一としました．

4. 最適化の結果

以下のpso()関数を実行して，最適化を行います．

@time res = pso(problem, iter=1000, nump=100)

何度か実行したところ，以下のような結果となりました．何回かに1回，近似解で止まってしまいますが，ほとんどの試行で最適解(評価値0)に収束しています．
また，Juliaは初回実行時はプリコンパイルがあるので実行時間が長いという話を聞きますが，初回は0.3秒，その後は0.1秒ほどと，その傾向がそのまま出ている気がします．

objectives: [0.0]
variables: [1.0, 1.0, 1.0]
  0.260429 seconds (1.75 M allocations: 165.164 MiB, 14.49% gc time)

objectives: [0.0]
variables: [1.0, 1.0, 1.0]
  0.164344 seconds (1.65 M allocations: 159.511 MiB, 35.80% gc time)

objectives: [2.1172833042747135]
variables: [0.04515555334956112, 0.0031784243458048422, -0.04600862315194654]
  0.088791 seconds (1.83 M allocations: 177.379 MiB, 22.27% gc time)

5. 課題

PSOの実装とRosenbrock関数の最適化ができることを確認しました．今回，アルゴリズムとしては以下ができていないので，今後，少しずつ取り組みたいと思っています．

多目的および制約付き最適化への拡張
冒頭に述べたように，拡張できるような配慮をしましたので，多目的・制約付き最適化でも適用できるようにしたいです．
突然変異の実装
オリジナルのPSOには突然変異はありませんが，何らかの手法を使えるようにしたいです．
w, c1, c2のパラメータ
w, c1, c2の与え方には様々な検討がなされています．別途指定できるようにしたいです．

また，Juliaの学習という意味では以下も行いたいと思います．

Juliaのデザインパターンの学習
今回のような設計が良かったかどうかの感覚があまりつかめていません．他の人のコードを見たりデザインパターンを調べて，実装もしながら慣れていきたいと思います．
探索経過を確認できるような可視化・グラフ描画
グラフ描画の練習も兼ねて，探索経過を可視化する機能を作りたいです．
モジュール分割
プログラム自体がやや長くなってしまったので，問題とアルゴリズムのような形でモジュール化し分割したいと思います．