Particle Swarm Optimization

PSO(particle swarm optimization) を C++ で実装しました。

1. PSOについて

日本語だと、粒子群最適化と訳されたりする、メタヒューリスティクス手法です。

多点で、勾配情報を用いない探索をすることが特徴に挙げられます。ちなみに、鳥の群れが餌を探すときの行動を模擬したモデルだそうな。

このアルゴリズムのメリット/デメリットは、以下。

●メリット

多点探索を行い、大域的な探索が可能
勾配情報を用いない為、微分不可能な関数でも探索が可能

●デメリット

複数パラメータとして乱数を用いるため、確率的な収束になる
解析的な手法でないため、最適解に収束する保証はない

2. 更新式

その点の今までの過去情報の中で、Bestな(最も関数値が小さい)ものp_bestと、グループ全体の中での、Bestな(最も関数値が小さい)ものg_bsetを保存しておいて、イテレーションの移動量に影響を与える。

p は点の数、λ に関しては、5.b.で説明。

3. 諸々の条件と結果

変数の次元	2
初期点数	5
ベンチマーク関数	2N-minima関数

結果の図をどど～んと載せる。

綺麗に、左下の(-3,-3)付近の大域的最適解が求められました。

4. コード

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
#include <iostream>
#include <fstream>
#include <math.h>

#define N 2  // 次元
#define I 5  // 初期点数
#define K 30 // イテレーション数
#define v_max 3 // 最大の移動量

using namespace std;

double __2n_minima (double *x);

int main()
{
	// file出力の為の定義
	ofstream fout;
	fout.open("pso/pso-example.txt");
	// 見出し
	fout << "k\tx1_1\tx1_2\tx2_1\tx2_2\tx3_1\tx3_2\tx4_1\tx4_2\tx5_1\tx5_2" << endl;

	// 初期点郡
	double x[I][N] = {
		{-4, 0},
		{-1, 3},
		{0, -2},
		{-4, -3},
		{3, -1},
	};
	// 各点での関数値とBest値
	double f_temp, f_pbest[I], f_gbest;
	// 各点の実際の値
	double x_pbest[I][N], x_gbest[N];
	// 移動量
	double v[I][N];

	// イテレーション0回目の設定
	fout << 0 << "\t";
	for (int i = 0; i < I; i++) {
		f_temp = __2n_minima(x[i]);
		f_pbest[i] = f_temp;
		for (int n = 0; n < N; n++) {
			x_pbest[i][n] = x[i][n];
			v[i][n] = 0;
			fout << x[i][n] << "\t";
		}
		if (i == 0 || f_gbest > f_temp) {
			f_gbest = f_temp;
			for (int n = 0; n < N; n++) {
				x_gbest[n] = x[i][n];
			}
		}
	}
	fout << endl;

	// 各パラメータ
	double c = 0.1;
	double λ, r1, r2;
	double λ_max = 1.0, λ_min = 0.6;

	for(int k = 1; k <= K; k++) {
		// Inertia Weight Approach
		λ = λ_max - (λ_max - λ_min) * k / K;
		r1 = ((double)rand() / ((double)RAND_MAX + 1)) * 1.0;
		r2 = ((double)rand() / ((double)RAND_MAX + 1)) * 1.0;

		// グループ内のBestを計算
		for (int i = 0; i < I; i++) {
			f_temp = __2n_minima(x[i]);
			if (f_pbest[i] > f_temp) {
				for (int n = 0; n < N; n++) {
					x_pbest[i][n] = x[i][n];
				}
				f_pbest[i] = f_temp;
				if (f_gbest > f_temp) {
					f_gbest = f_temp;
					for (int n = 0; n < N; n++) {
						x_gbest[n] = x[i][n];
					}
				}
			}
		}

		fout << k << "\t";
		// イテレーション毎のgbest を用いるため、for 文は分ける
		for (int i = 0; i < I; i++) {
			for (int n = 0; n < N; n++) {
				v[i][n] = λ * v[i][n] + r1 * (x_pbest[i][n] - x[i][n]) + r2 * (x_gbest[n] - x[i][n]);
				// 発散を防ぐ
				if (v[i][n] < -1 * v_max) {
					v[i][n] = -1 * v_max;
				} else if (v[i][n] > v_max) {
					v[i][n] = v_max;
				}

				x[i][n] = x[i][n] + c * v[i][n];
				fout << x[i][n] << "\t";
			}
		}
		fout << endl;
	}
}

/*
 * 目的関数 2N-minima 関数
 */
double __2n_minima (double x[N])
{
	double f;
	f = pow(x[0], 4) - 16 * pow(x[0], 2) + 5 * x[0] + pow(x[1], 4) - 16 * pow(x[1], 2) + 5 * x[1];
	return f;
}

5. 補足

a. イテレーションの移動量の制限を設ける

1
#define v_max 3 // 最大の移動量

あまり必要ないですが、一度のイテレーションでの移動幅のMaxを設定してます。

b. 慣性係数 λ について

1イテレーション前の移動量項の係数 λ の決定を乱数で行うことも出来るのですが、より精度を上げるために、Inertia Weight Approach(IWA)という手法を用います。

これはイテレーション経過とともに、λを小さくしていき、 探索序盤では、大域的探索が期待できる不安定な探索を実行し、探索終盤は、局所的な安定な探索を行えるようにする手法です。