人工蜂コロニーアルゴリズムで関数の最小値探索を行う
概要
大学の授業で群知能についていくつか扱ったのでその復習に人工蜂コロニーアルゴリズムを実装します。
この記事の対象
本編
群知能とは
自然界には多くのアイデアが眠っています。蓮の葉から水をはじく構造のヒントを得たとか、オナモミからマジックテープのヒントを得たとか、有名な話だと思います。
群知能もそんな風に自然界から着想された一つです。群知能とは「鳥の群れのように個体間の局所的な関係を合わせて、集団として高度な動きを見せる現象を模した計算手法」です。粒子群最適化（PSO）や蟻コロニー最適化（ACO）などいろいろありますが、今回は人工蜂コロニーアルゴリズム（artificial bee colony algorhythm , 以下ABCと表記）について実装します。
アルゴリズム
収穫蜂（employed bee）、追従蜂（onlooker bee）、偵察蜂（scout bee）の３つの種類の蜂の動きを模した操作で探索点を更新していきます。
準備
収穫蜂N匹をアルゴリズムを適応する空間（D次元）にランダムに初期配置する。全ての収穫蜂に対してカウンタTCというものをそれぞれ定義し、初期値0にしておく。Step1 ー収穫蜂ー
(1) 収穫蜂iについて順番に次の更新式で探索点を更新していく。
(2) このとき新探索点候補が現在の探索点より適合度が高ければ（最小値問題なら関数値が小さければ）、探索点を更新しカウンタを0にする。反対に現在の探索点の方が適合度が高ければ更新はせず代わりにカウンタを+1する。 数式で書くと以下のようになる。1
この作業を全ての収穫蜂について繰り返す。
Step2 ー追従蜂ー
追従蜂N匹を使ってさらに探索を進める。
追従蜂は、以下のように求められる確率を使って収穫蜂１匹をルーレット選択し、Step1の(1)(2)の作業を行う。
この作業を全ての追従蜂について繰り返す。Step3 ー偵察蜂ー
カウンタTCがあらかじめ設定した閾値を超えている場合、その収穫蜂の探索点を初期化する。Step1~3を規定回数行ったら終了とする。（あるいはあらかじめ決めておいた終了条件を満たした、などでも良い）
探索の中で見つかった最良の探索点が近似解である。
Pythonで実装
import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import math # 初期設定 N = 100 # 個体数 d = 2 # 次元 TC = np.zeros(N) #更新カウント lim = 30 xmax = 5 xmin = -5 G = 100 # 繰り返す周期 x = np.zeros((N,d)) for i in range(N): x[i] = (xmax-xmin)*np.random.rand(d) + xmin x_best = x[0] #x_bestの初期化 best = 100 # ルーレット選択用関数 def roulette_choice(w): tot = [] acc = 0 for e in w: acc += e tot.append(acc) r = np.random.random() * acc for i, e in enumerate(tot): if r <= e: return i #目的関数（例としてsphere関数） def func(x): return np.sum(x**2) # アルゴリズム best_value = [] for g in range(G): best_value.append(func(x_best)) # employee bee step for i in range(N): v = x.copy() k = i while k == i: k = np.random.randint(N) for j in range(d): r = np.random.rand()*2-1 #-1から1までの一様乱数 v[i,j] = x[i,j] + r * (x[i,j] - x[k,j]) if func(x[i]) > func(v[i]): x[i] = v[i] TC[i] = 0 else: TC[i] += 1 # onlooker bee step for i in range(N): w = [] for j in range(N): w.append(np.exp(-func(x[j]))) i = roulette_choice(w) for j in range(d): r = np.random.rand()*2-1 #-1から1までの一様乱数 v[i,j] = x[i,j] + r * (x[i,j] - x[k,j]) if func(x[i]) > func(v[i]): x[i] = v[i] TC[i] = 0 else: TC[i] += 1 # scout bee step for i in range(N): if TC[i] > lim: for j in range(d): x[i,j] = np.random.rand()*(xmax-xmin) + xmin TC[i] = 0 # 最良個体の発見 for i in range(N): if best > func(x[i]): x_best = x[i] best = func(x_best) print(x_best,func(x_best)) plt.plot(range(G),best_value) plt.yscale('log') plt.title("ABC") plt.show()
今回は２次元のsphere関数（z = x2 + y2,-5<x,y<5)の最小値を求めてみます。アルゴリズム自体は上で説明した通りです。工夫というか疑問点として、適合度関数を決めるときに、最小値探索なので値が低いほど適合度関数は高くなるように設定するべきだと思うのですが、確率計算の都合上全て正の値をとる関数にしたくて、結局exp(-x)を採用しました。（これでよかったのか…？という疑問が残りますが）
追記：結局exp(-x)ではなくシンプルに「関数値が小さい＝適合度が高い」としてます。
実行と軽い考察
Sphere関数
実行してみると、以下のように最小値が更新されていく様子が見えます（縦軸がログスケールであることに注意）。
得られた最良点は[7.25e-15 9.89e-15]で、そのときの目的関数の値は1.50e-28と求められました。
z = x2 + y2の最小値はx=y=0のとき0なので、ほぼ正しい解が求められていることがわかります。
どのようにミツバチが動いているのかもう少し確認してみます。以下は20周期ごとの収穫蜂の探索点の三次元プロットです。
縦軸の尺度の変化に注意すれば、最小値付近にすごいスピードで集まっていることがわかると思います。
次にこのsphere関数の最小値をアルゴリズムと同時に動かしてみます。
こんな具合です。上の黒い点が動くのを追って周りの収穫蜂たちがわらわらと群がる感じのアニメーションです2。（現実に蜂にこんなに群がられたらショック死する自信あります…）
Griewank関数
他の関数も試してみましょう。最小値問題のアルゴリズム評価によく使われるベンチマーク関数の一つ、Griewank関数というものを試してみます。
どういう関数かここには貼りませんので検索してみてください。たわしみたいな関数です（語彙力）。
これは多峰性関数といって局所的最適解が非常に多く存在しています。これにABCを用いると以下のようになりました。
得られた最適解は[-253 -0.1]で最小値は1.21e-06でした。Griewank関数は[0,0]で最小値0を取るので、今回のABCだと局所解にはまってしまったようです。学習はうまく進んでいるんですけどね。
そもそもこれは相手が強すぎるというか大域的最適解探索をとことん妨害するために作られた関数だと思うのでしょうがないかなと思います。
まとめ
今回はABCを実装しました。他にも群知能アルゴリズムはいろいろあるので、次の機会にこんな感じでまた実装したいと思います。
ABCについても、他にも「偵察蜂のSTEP3を消してしまう」や「カウンタの閾値を変化させる」、「他の様々なベンチマーク関数に試す」、「次元を増やす」などいろいろ遊んでみると発見があって面白いと思います。
おまけ
ミツバチに追われる黒い点のgifをぜひ皆さんにも実行して欲しいと思ったのでソースコードを置いておきます。整理してないので汚いですが貼ってpythonで実行するだけなので是非やってみてください。 蜂に追われる恐怖が少しわかります。
import numpy as np import matplotlib.pyplot as plt import math from matplotlib import animation # 関数の設定 # xはnp.array #sphere def func(x): return np.sum(x**2) # 初期設定 N = 100 # 個体数 d = 2 # 次元 TC = np.zeros(N) #更新カウント lim = 30 xmax = 5 xmin = -5 G = 300 # 繰り返す回数 x = np.zeros((N,d)) for i in range(N): x[i] = (xmax-xmin)*np.random.rand(d) + xmin # ルーレット選択用関数 def roulette_choice(w): tot = [] acc = 0 for e in w: acc += e tot.append(acc) r = np.random.random() * acc for i, e in enumerate(tot): if r <= e: return i x_best = x[0] #x_bestの初期化 best = 100 # 繰り返し fig = plt.figure() best_value = [] ims = [] for g in range(G): centroid = 5*math.sin(0.07*g) def func(x): return np.sum((x-centroid)**2) best_value.append(best) z = np.zeros(N) for i in range(N): z[i] = func(x[i]) img2 = plt.scatter(centroid,centroid,color='black') img1 = plt.scatter(x[:,0],x[:,1],marker='h',color='gold') ims.append([img1]+[img2]) # employee bee step for i in range(N): v = x.copy() k = i while k == i: k = np.random.randint(N) for j in range(d): r = np.random.rand()*2-1 #-1から1までの一様乱数 v[i,j] = x[i,j] + r * (x[i,j] - x[k,j]) if func(x[i]) > func(v[i]): x[i] = v[i] TC[i] = 0 else: TC[i] += 1 # onlooker bee step for i in range(N): w = [] for j in range(N): w.append(np.exp(-func(x[j]))) i = roulette_choice(w) for j in range(d): r = np.random.rand()*2-1 #-1から1までの一様乱数 v[i,j] = x[i,j] + r * (x[i,j] - x[k,j]) if func(x[i]) > func(v[i]): x[i] = v[i] TC[i] = 0 else: TC[i] += 1 # scout bee step for i in range(N): if TC[i] > lim: for j in range(d): x[i,j] = np.random.rand()*(xmax-xmin) + xmin TC[i] = 0 # 最良個体の発見 for i in range(N): if best > func(x[i]): x_best = x[i] best = func(x_best) # 結果 print(x_best,func(x_best)) ani = animation.ArtistAnimation(fig, ims, interval=70) plt.show()
