はじめに

最適化アルゴリズムの実装シリーズです。
まずは概要を見てください。

コードはgithubにあります。

粒子群最適化

概要

粒子群最適化(Particle Swarm Optimization：PSO)とは、スズメやイワシといった小さい個体が大きな群れを作って効率よく餌を探す行動に着目して作られたアルゴリズムです。
群れに属する個体は以下のような行動モデルに基づいているといわれています。

近くにいる個体に影響されて行動する
他の個体の近くにとどまろうとするが、一定以上は近づかない
他の個体の速度に合わせて移動する

参考
・粒子群最適化法(PSO)を救いたい
・進化計算アルゴリズム入門生物の行動科学から導く最適解(amazonリンク)

アルゴリズム

各個体を粒子と呼びます。
全粒子の中で今までに一番良かった位置をグローバルベスト、自分自身がこれまでに見つけた位置の中で一番良かった位置をパーソナルベストとして覚えておきます。
これら2つの情報を元に次の自分の位置を更新していく探索が粒子群最適化となります。

アルゴリズムのフロー

問題との対応

問題	粒子群最適化
入力値の配列	粒子
入力値	粒子の座標
評価値	粒子の評価値

コード内で使う変数の意味

変数名	数式	意味	所感
problem		任意の問題のクラス（問題編を参照）
particle_max		粒子の数
inertia	$I$	慣性係数	加速度の減少率
global_acceleration	$A_g$	加速係数(グローバルベスト)	グローバルベストに近づく割合
personal_acceleration	$A_p$	加速係数(パーソナルベスト)	パーソナルベストに近づく割合
particles		粒子の配列
particles[i]["personal"]		i番目の粒子のパーソナルベスト
particles[i]["v"]		i番目の粒子の速度

粒子の位置と加速度の更新

加速度と座標の初期値は空間内の範囲でランダムに初期化します。

位置を求めるにあたりまずは加速度が更新されます。
時刻(t)から時刻(t+1)への速度はグローバルベストとパーソナルベストを用いて以下の式となります。

$$
\vec{v_i}(t+1) = I \vec{v_i}(t)

A_g ( \vec{g}(t) - \vec{x_i}(t) ) \times rand[0,1]
A_p ( \vec{p}(t) - \vec{x_i}(t) ) \times rand[0,1]
$$

$I$ は慣性係数、$A_g$と$A_p$は加速係数で3つとも1より小さい値を指定します。
また、$\vec{g}$はグローバルベストの位置を表し、$\vec{p}$はその粒子が持っているパーソナルベストの位置を表します。
$rand[0,1]$ は0～1の乱数です。

コード全体

コード全体です。

import math
import random

import numpy as np

class PSO():
    def __init__(self,
        particle_max,
        inertia=0.9,
        global_acceleration=0.9,
        personal_acceleration=0.9,
    ):
        self.particle_max = particle_max
        self.inertia = inertia
        self.global_acceleration = global_acceleration
        self.personal_acceleration = personal_acceleration


    def init(self, problem):
        self.problem = problem

        # 初期粒子群を生成
        self.global_best = None
        self.particles = []
        for _ in range(self.particle_max):
            o = problem.create()  # ランダムな位置に粒子を作成

            # 初期加速度
            v = [(problem.MAX_VAL - problem.MIN_VAL) * random.uniform(-1, 1) for _ in range(problem.size)]

            # パーソナルベストと速度の情報を付与する
            d = {
                "particle": o,
                "personal": None,
                "v": np.asarray(v),
            }
            self.particles.append(d)

            # パーソナルベストとグローバルベストの更新
            self._updateBest(d)
    

    def step(self):

        # 各粒子に対して
        for particle in self.particles:
            # 各座標を出力(numpy化してベクトル計算をしやすくしています)
            pos = np.asarray(particle["particle"].getArray())
            g_pos = np.asarray(self.global_best.getArray())
            p_pos = np.asarray(particle["personal"].getArray())

            # 加速度を計算
            v = particle["v"]
            v = self.inertia * v
            v += self.global_acceleration * (g_pos - pos) * random.random()
            v += self.personal_acceleration * (p_pos - pos) * random.random()
            particle["v"] = v

            # 座標を更新
            particle["particle"].setArray(pos + v)

            # パーソナルベストとグローバルベストの更新
            self._updateBest(particle)


    def _updateBest(self, particle):

        # パーソナルベストの更新
        if particle["personal"] is None or particle["personal"]["particle"].getScore() < particle["particle"].getScore():
            # パーソナルベストとしてコピーして保存する
            particle["personal"] = particle["particle"].copy()

        # グローバルベストの更新
        if self.global_best is None or self.global_best.getScore() < particle["particle"].getScore():
            # グローバルベストとしてコピーして保存する
            self.global_best = particle["particle"].copy()

ハイパーパラメータ例

各問題に対して optuna でハイパーパラメータを最適化した結果です。
最適化の1回の試行は、探索時間を2秒間として結果を出しています。
これを100回実行し、最適なハイパーパラメータを optuna に探してもらいました。

問題	global_acceleration	inertia	particle_max	personal_acceleration
EightQueen	0.7839287773369192	0.8235075101639827	28	0.9467161852490191
function_Ackley	0.11778673996028691	0.062427889313103786	47	0.6739352477292235
function_Griewank	0.1939265329859151	0.9886894890970368	21	0.9870275206300417
function_Michalewicz	0.5004800811479669	0.7828732926170527	23	0.5089894615381799
function_Rastrigin	0.7067489271622391	0.8580154745241855	31	0.2988574441297245
function_Schwefel	0.6251927739536115	0.9030347446658787	42	0.2540225969044799
function_StyblinskiTang	0.970834173658411	0.9365843106326938	31	0.9298455482443944
LifeGame	0.539526858214075	0.08845815509172461	14	0.4249882950339423
OneMax	0.5056801912085691	0.6273453999411102	42	0.29656995228071314
TSP	0.7018381808726923	0.888427895281042	48	0.6464768696714887

実際の動きの可視化

1次元は6個体、2次元は20個体で100step実行した結果です。
赤い丸がそのstepでの最高スコアを持っている個体です。

パラメータは以下で実行しました。

PSO(N, inertia=0.2, global_acceleration=0.2, personal_acceleration=0.2)

OneMax

1次元

2次元

function_Ackley

1次元

2次元

function_Griewank

1次元

2次元

function_Michalewicz

1次元

2次元

function_Rastrigin

1次元
2次元

function_Schwefel

1次元

2次元

function_StyblinskiTang

1次元

2次元

function_XinSheYang

1次元

2次元

あとがき

現状の最適解にみんなが集まってくる様子が粒子群っぽくていいですね。
アルゴリズム的には、一度局所解に入ってしまうとみんなそこに集まってしまうので抜け出せなさそうなイメージです。
また現状の最適解周辺でうろうろしている個体はパーソナルベストとグローバルベストの両方に引っ張られて釣り合っている個体ですね。

今回実装した内容は一番シンプルなPSOなので、改良されたPSOならこれらの問題点が解決されているかもしれません。