粒子群最適化を使って関数最小化問題を解いてみる

粒子群最適化

粒子群最適化(PSO:Particle Swarm Optimization)とは, 群知能の一種であり, 解探索手法として組合せ最適化問題に使用されます.
粒子の速度と位置という2つの情報の更新を繰り返して, 探索を進めていく流れになります.
下の図は粒子群最適化のイメージです.

##アルゴリズム
粒子群最適化のアルゴリズムは以下のようになります.

##更新式
以下の更新式より粒子の速度と位置を更新します.
簡単に説明すると, 速度は粒子を進化させる方向を表し, 位置は粒子自身のパラメータを表します.

実験

実際に最適化問題を解いていきましょう.
今回は

x^2 + y^2

の最小化問題を解いていきましょう.
ですので, (x,y)=(0,0)が最適解となります.
使用したコードは以下のようになります.

# -*- coding: utf-8 -*-
import numpy as np
import random

# 評価関数
def evaluate(particle):
    z = 0
    for i in range(len(particle)):
        z += particle[i] ** 2
    return z

# 位置更新
def update_position(particle, velocity):
    new_particle = particle + velocity
    return new_particle

# 速度更新
def update_velocity(particle, velocity, pbest, gbest, w=0.5, max=0.15):
    new_velocity = np.array([0.0 for i in range(len(particle))])
    #new_velocity = [0.0 for i in range(len(particle))]
    r1 = random.uniform(0, max)
    r2 = random.uniform(0, max)
    for i in range(len(particle)):
        new_velocity[i] = (w * float(velocity[i]) + r1 * (float(pbest[i]) - float(particle[i])) + r2 * (float(gbest[0]) - float(particle[i])))

    return new_velocity

def main():
    N = 100 # 粒子の数
    length = 2  # 次元数
    para_max = 100  #パラメータの最大値
    # 粒子位置の初期化
    ps = [[random.uniform(-para_max, para_max) for j in range(length)] for i in range(N)]
    vs = [[0.0 for j in range(length)] for i in range(N)]
    # パーソナルベスト
    personal_best_position = ps
    # パーソナルベストの評価
    personal_best_scores = [evaluate(p) for p in ps]
    # 評価値が最も小さい粒子のインデックス
    best_particle = np.argmin(personal_best_scores)
    # グローバルベスト
    global_best_position = personal_best_position[best_particle]

    generation = 30 # 最大世代数
    # 世代数分ループ
    for t in range(generation):
        file = open("data/pso/pso" + str(t+1) + ".txt", "w")
        # 粒子数分ループ
        for n in range(N):
            # ファイル書き込み
            file.write(str(ps[n][0]) + " " + str(ps[n][1]) + "\n")
            # 粒子の速度の更新
            vs[n] = update_velocity(ps[n], vs[n], personal_best_position[n], global_best_position)
            # 粒子の位置の更新
            ps[n] = update_position(ps[n], vs[n])

            # 評価値計算をしてパーソナルベストを求める
            score = evaluate(ps[n])
            if score < personal_best_scores[n]:
                personal_best_scores[n] = score
                personal_best_position[n] = ps[n]
        # グローバルベストの更新をする
        best_particle = np.argmin(personal_best_scores)
        global_best_position = personal_best_position[best_particle]
        file.close()

    print(global_best_position)
    print(min(personal_best_scores))

if __name__ == '__main__':
    main()

##実験結果
1,10,20,30世代目の個体を色分けして図にプロットしています.
世代が進むに連れて, 粒子が(0,0)に収束しているのが確認できますね.