LoginSignup
17
12

More than 3 years have passed since last update.

最適化アルゴリズムを実装していくぞ(粒子群最適化)

Posted at

はじめに

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

コードはgithubにあります。

粒子群最適化

概要

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

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

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

アルゴリズム

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

  • アルゴリズムのフロー

draw2-PSO.png

  • 問題との対応
問題 粒子群最適化
入力値の配列 粒子
入力値 粒子の座標
評価値 粒子の評価値
  • コード内で使う変数の意味
変数名 数式 意味 所感
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次元

OneMax_PSO_2.gif

  • 2次元

OneMax_PSO_3.gif

function_Ackley

  • 1次元
  • function_Ackley_PSO_2.gif

  • 2次元

function_Ackley_PSO_3.gif

function_Griewank

  • 1次元

function_Griewank_PSO_2.gif

  • 2次元

function_Griewank_PSO_3.gif

function_Michalewicz

  • 1次元

function_Michalewicz_PSO_2.gif

  • 2次元

function_Michalewicz_PSO_3.gif

function_Rastrigin

  • 1次元
    function_Rastrigin_PSO_2.gif

  • 2次元

function_Rastrigin_PSO_3.gif

function_Schwefel

  • 1次元

function_Schwefel_PSO_2.gif

  • 2次元

function_Schwefel_PSO_3.gif

function_StyblinskiTang

  • 1次元

function_StyblinskiTang_PSO_2.gif

  • 2次元

function_StyblinskiTang_PSO_3.gif

function_XinSheYang

  • 1次元

function_XinSheYang_PSO_2.gif

  • 2次元

function_XinSheYang_PSO_3.gif

あとがき

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

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

17
12
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
17
12