粒子群最適化(PSO)のパラメータについて

はじめに

これは「群知能と最適化2019 Advent Calender」5日目の記事です。
Qiitaに記事を投稿するのは初めてなのでつたない部分もあるかもですが、読んでいただけると嬉しいです。
また、この記事はPSOのパラメータに関する記事なので、そもそもPSOについて知らないという方は昨日のganariyさんの記事「粒子群最適化法(PSO)を救いたい」を先に読んでいただくとより理解できるかと思います。

粒子群最適化(PSO)とパラメータについて

PSOの詳細については昨日の記事でganariyさんが述べてくださった通り、鳥や魚などの群れで行動する習性を基にしたメタヒューリスティックアルゴリズムの一つです。
1995年に基となる論文が発表されて以降、様々なモデルが提案されてきました。
PSOは各粒子が現在の位置$x$と速度$v$を持ち、それらを更新していくことで最適解を探索します。移動にかかる各粒子の速度と位置は以下の式で更新されます。
$$
v_{id}(t+1) = w_{id}(t) + c_p × r_p × (pBest_{id} - x_{id}(t)) + c_g × r_g × (gBest_{id} - x_{id}(t))
$$
$$
x_{id}(t+1) = x_{id}(t) + v_{id}(t+1)
$$
ここで、それぞれが見つけた最良の位置をパーソナルベスト(pBest)、群全体の最良の位置をグローバルベスト(gBest)と呼びます。
また、$v_{id}$はd次元のi番目の粒子の速度、$x_{id}$はd次元のi番目の粒子の位置、$w$は重み、$c_p$、$c_g$は加速係数と呼ばれる定数、$r_p$、$r_g$は0~1の乱数です。
この中で人が決めることのできるパラメータとして、重み$w$、加速係数$c_p$$c_g$がありますが、特に$w$は粒子の収束に大きな影響を与えるパラメータとして多くの研究が行われています。

重みwに関して行われている研究

重みwに関するレビュー論文として、2011年に「A novel particle swarm optimization algorithm with adaptive inertia weight」という論文が報告されています。
この論文の中で、重みwは主に3つのタイプに分類されています。

1. 重みが定数orランダム
   定数を与える論文では、様々な定数での重みが試され、$w>1.2$では浅い探索となってしまう一方、$w<0.8$では局所解に収束する傾向があることを示しました。
   また、ランダム値を与える論文では、最適な重みを動的な環境で決めるのは難しいため、以下の式のようなランダムな値を与えることを提案しています。
   $$
   w = 0.5 + rand()/2
   $$

2. 時変する重み
   PSOを用いた多くの手法では、重みの値を反復回数に基づいて決定する手法を用いています。
   なかでも特に多くの論文で用いられている手法として、線形減少重みがあります[参考: 1, 2, 3]。これは、最大反復回数$iter_{max}$と反復回数$iter$、初期重み$w_{max}$と最終値$w_{min}$を用いて以下の式で表されます。
   $$
   w(iter) = ((iter_{max} - iter)/iter_{max})(w_{max} - w_{min}) + w_{min}
   $$
   時変する重みに関しては、これ以外にもカオスモデルを使用した手法や非線形減少重みなど様々な重みが提案されています。

3. 適応重み
   これは適宜状況を監視し、フィードバックパラメータに基づいて重みを決定していくものです。非常に様々な手法が提案されているので、興味があればレビュー論文を読んでみてください。

LDWPSOの実装

上記した中でも、最も有名で多く使われる線形減少重みPSO (LDWPSO) を実装してみたいと思います。
評価関数にはSphere functionを用います。
プログラムは以下のようになります。
実行すると4ループ目まで可視化して表示されるようにしています。

import numpy as np
import random
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import time

SWARM_SIZE = 100 # 粒子数
ITERATION = 30 # PSOループ回数
C1 = C2 = 2.0 # 加速係数
MIN_X, MIN_Y = -5.0, -5.0 # 探索開始時の範囲最小値
MAX_X, MAX_Y = 5.0, 5.0 # 探索開始時の範囲最大値
W_MAX, W_MIN = 0.9, 0.4

# 評価関数(z = x^2+y^2)
def evaluation(x, y):
    return x**2 + y**2 # Sphere function

# PSO
def pso(position, velocity, personal_best_scores, personal_best_positions, best_position, search_space_x, search_space_y):
    # ループ回数分Particle移動
    for i in range(ITERATION):
        for s in range(SWARM_SIZE):
            # 変更前の情報の代入
            x, y = position[s]["x"], position[s]["y"]
            vx, vy = velocity[s]["x"], velocity[s]["y"]
            p = personal_best_positions[s]
            
            # 粒子の位置更新
            new_x, new_y = update_position(x, y, vx, vy, search_space_x, search_space_y)
            position[s] = {"x": new_x, "y": new_y}
            # 粒子の速度更新
            new_vx, new_vy = update_velocity(new_x, new_y, vx, vy, p, best_position, i)
            velocity[s] = {"x": new_vx, "y": new_vy}

            # 評価値を求める
            score = evaluation(new_x, new_y)
            # update personal best
            if score < personal_best_scores[s]:
                personal_best_scores[s] = score
                personal_best_positions[s] = {"x": new_x, "y": new_y}
        # update global best
        best_particle = np.argmin(personal_best_scores)
        best_position = personal_best_positions[best_particle]

        # PSO Visualization
        
        if i == 0 or i == 1 or i == 2 or i == 3:
            print("ITERATION = " + str(i+1))
            visualization(personal_best_positions)
        

# 粒子の位置更新関数
def update_position(x, y, vx, vy, search_space_x, search_space_y):
    new_x = x + vx
    new_y = y + vy
    # 探索範囲内か確認
    if new_x < search_space_x["min"] or new_y < search_space_y["min"]:
        new_x, new_y = search_space_x["min"], search_space_y["min"]
    elif new_x > search_space_x["max"] or new_y > search_space_y["max"]:
        new_x, new_y = search_space_x["max"], search_space_y["max"]
    return new_x, new_y

# 粒子の速度更新関数
def update_velocity(x, y, vx, vy, p, g, i):
    #random parameter (0~1)
    r1 = random.random()
    r2 = random.random()

    # 線形減少重み
    L = (ITERATION - i) / ITERATION
    D = W_MAX - W_MIN
    w = L * D + W_MIN
    print(w)

    # 速度更新
    new_vx = w * vx + C1 * (g["x"] - x) * r1 + C2 * (p["x"] - x) * r2
    new_vy = w * vy + C1 * (g["y"] - y) * r1 + C2 * (p["y"] - y) * r2
    return new_vx, new_vy

# 可視化関数
def visualization(positions):
    fig = plt.figure()
    ax = Axes3D(fig)
    # Mesh
    #mesh_x = np.arange(-0.5e-3, 0.5e-3, 0.1e-4)
    #mesh_y = np.arange(-0.5e-3, 0.5e-3, 0.1e-4)
    mesh_x = np.arange(-5.0, 5.0, 0.5)
    mesh_y = np.arange(-5.0, 5.0, 0.5)
    mesh_X, mesh_Y = np.meshgrid(mesh_x, mesh_y)
    mesh_Z = evaluation(mesh_X, mesh_Y)
    ax.plot_wireframe(mesh_X, mesh_Y, mesh_Z)
    # Particle
    for j in range(SWARM_SIZE):
        z = evaluation(positions[j]["x"], positions[j]["y"])
        ax.scatter(positions[j]["x"], positions[j]["y"], z)
    ax.set_xlim([-5.0, 5.0])
    ax.set_ylim([-5.0, 5.0])
    ax.set_zlim([-1.0, 30.0])
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("f(x, y)")
    plt.show()

def run(position, velocity, personal_best_scores, personal_best_positions, best_position, search_space_x, search_space_y):
    pso(position, velocity, personal_best_scores, personal_best_positions, best_position, search_space_x, search_space_y)
    return best_position, personal_best_scores

def main():
    # 時間計測開始
    start = time.time()

    # 各粒子の初期位置, 速度, personal best, global best 及びsearch space設定
    position = []
    velocity = []
    personal_best_scores = []
    # 初期位置, 初期速度
    for s in range(SWARM_SIZE):
        position.append({"x": random.uniform(MIN_X, MAX_X), "y": random.uniform(MIN_Y, MAX_Y)})
        velocity.append({"x": random.uniform(0, 1), "y": random.uniform(0, 1)})
    # personal best
    personal_best_positions = list(position)
    for p in position:
        personal_best_scores.append(evaluation(p["x"], p["y"]))
    # global best
    best_particle = np.argmin(personal_best_scores)
    best_position = personal_best_positions[best_particle]
    # search space
    search_space_x, search_space_y = {'min': MIN_X, "max": MAX_X}, {"min": MIN_Y, "max": MAX_Y}

    # run
    best_position, personal_best_scores = run(position, velocity, personal_best_scores, personal_best_positions, best_position, search_space_x, search_space_y)

    # 時間計測終了
    process_time = time.time() - start

    # Optimal solution
    print("Best Position:", best_position)
    print("Score:", min(personal_best_scores))
    print("time:", process_time)

if __name__ == '__main__':
    main()

本当はシンプルなPSOと比較するところまでできればいいのですが、気力が尽きたのでここまでで、、、$w_{min}$の値を0.9に変更するだけでシンプルなPSOになるので興味があれば試してみてください。

まとめ

かなりアカデミックな記事になってしまいましたが、まとめると、
・PSOの重み$w$は最適解への収束に大きな影響を与える
・重みは主に3つのタイプがある
・最もよく使われる重みとして、時変する重みの一つである線形減少重みがある
となります。
このように、パラメータ一つでも様々な研究が行われているということを知っていただければ嬉しいです。