皆さんは、郡知能（Swarm Intelligence）という言葉をご存知でしょうか？
自然界の中で、群れ（Swarm）は個体からは想像できない行動をとることがあります。たとえば、ガチョウが群れで飛ぶ時は、なぜか三角形の形状になります。
こうした自然界の群れ現象を解明し、アルゴリズムに応用したものが群知能です。
こうした群れが構成されるのは、指揮命令個体が指示しているわけはなく、また群れ全体が何らかの方法でコミュニケーションしているわけではありません。個体が持つちょっとしたルールによって群れ現象が構成されます。コンピュータ・グラフィックス（CG）の世界では、このちょっとしたルールを再現することで群れのCGを作り出しいるそうです。（専門外なので詳しくは知りません）
アルゴリズムの世界では、この群れの行動を主として最適化アルゴリズムに利用しています。元々は、人工生命（Artificial Life）分野で研究されていたそうですが、最近は群知能という分野名になっているようです。その中で、最も有名なのは人工生命で利用されていた蟻の巣最適化アルゴリズム（Ant Corony Optimization）かと思われます。他にも生物を模したアルゴリズムが多数あり、現在でも新しいアルゴリズムが発見されています。
Iztok Fister Jr. らによれば、自然界に触発された最適化アルゴリズム（Nature-inspired optimization algorithm）は、次のように分類されています。
- swarm intelligence-based
- bio-inspired (not swarm)
- physics/chemistry-based
- others
つまり、自然界の生き物の中で群れに着目したアルゴリズムという位置付けです。
さて、このような群知能分野では、新しいアルゴリズムがどんどん生み出されているそうなので、それらをまとめた本でもないかなぁと思い、だいぶ前に以下の書籍を購入しました。
参考文献
Swarm Intelligence: Principles, Advances, and Applications
この本では、最初にこうもりを模したアルゴリズムの概略が記載されていたので、これを実装してみました。
こうもりの捕獲メカニズム
有名な話だと思いますが、こうもりは人間の耳には聞こえないパルス音波を発し、その周囲から反射された音を感知して、壁や餌の位置を把握（エコロケーション）しています。これを理想化し、３つのシンプルなルールにします。
- 全てのこうもりは、エコロケーションを使って距離を測り、餌や壁の違いも把握している
- こうもりは、位置x(i)で速度v(i)でランダムに飛んでおり、周波数fmin（あるいは波長λ）、振幅A0の音波で餌を探している。この時、放出されるパルスの波長と放出量は捕獲対象に合わせて自動的に調整（アジャスト）される
- 振幅は様々な方法で変化するが、ここでは正の大きなA0から、最小値Aminまで変化すると仮定する
こうもりの位置と速度の更新
こうもりは、ランダムに飛んでいるので、ブラウン運動のように乱数を使って位置と速度を更新します。ここでは、こうもりの番号をi、時間ステップをtで表すことにします。
１．周波数の計算
f(i) = fmin + (fmax-fmin)*beta,
betaは[0,1]の一様乱数
２．速度の計算
v(i,t) = v(i,t-1) + (x(i,t)-xbest)*f(i),
xbestは現在の餌の最良位置
３．位置の計算
x(i,t) = x(i,t-1) + v(i,t)
betaのランダム性がランダムに飛んでいることを表しています。位置の計算では、時間ステップdt=1を利用しているため、次元も問題ありません。（本来なら、x+vdtのはず）
こうもりの獲物探索
獲物探索は、ランダムウォークを使って局所探索（ローカルサーチ）をするとは書いてあるものの、具体的な方法は書いていないため、独自の方法を作りました。ただし、探索範囲については振幅（音の大きさ）によって変わるので、本の通りに次の公式を用いています。
X_new = X_old + eps * A_t
ここで、epsは[-1,1]の正規分布による乱数、A_tは全こうもりの振幅の平均値です。
実装上では、目的関数g(x(i))に対して、次のようにしました。
- gbest=g(x(i))を計算する
- 乱数ベクトルepsを生成する
- gnew=g(x(i)+eps*A_t)を計算する
- gnew < gbestなら、gbest=gnewとする
- ２〜４を100回繰り返す
こうすると、100回の反復の中でg(x)が最小の位置x(i,new)が求まります。
パルスの放出量と振幅の更新
パルスの放出量rと振幅Aは、書籍の通りに下記の公式で実装しました。
r(i,t+1) = r(i,0)*{ 1 - exp(-gamma*t) }
A(i,t+1) = alpha * A(i,t)
ここで、r(i,0)は初期放出量、alpha(<1)とgammaは定数です。
テスト関数
今回は、目的関数を３次元の二重井戸関数にし、その最小値が半径1の球になるような問題設定にしました。
g(x(i)) = x(i)^4 - 2x(i)^2 + 1
G(X) = g(x(0)) + g(x(1)) + ... + g(x(n))
１次元の関数g(x)は、以下のような図になります。見ての通り、x=-1,1でg(x)は最小になります。
バットアルゴリズム
バットアルゴリズムは、こうもりをエージェントとした最適化アルゴリズムです。アルゴリズムの概要は次のようになります。
こうもり数nとし、初期位置xをランダムに生成する。
最初の最良位置（g(x)が最小の位置）をxbestとする。
for tが最大反復maxiteより小さい do
foreach こうもりIがこうもり数nより小さい do
位置と速度を更新する
endfor
if Aの平均値 < 0.85 then
反復終了
endif
乱数randを生成する
if rand > r(i) then
g(x(i,t))を使って局所探索し、x(i,new)を決める
endif
if rand < A(i) then
if 新しいx(i,new)を使って、G(Xnew) < G(Xold)となる then
こうもりiの位置を更新する：x(i,t)=x(i,new)
パルス放出量r(i,t)を更新する
パルス振幅A(i,t)を更新する
endif
endif
最良な位置xbestを更新する
endfor
実行結果
下記のアニメーションは、上記の二重井戸関数に対して、こうもり５０匹で解を探索したものになります。こうもりを青ドット、最良位置を赤ドットで表しています。時間が経つと、こうもりが球形に集まっていくのが分かると思います。
アニメーション出力は、３D散布図をGIF回転アニメーションにするを参考させてもらいました。
ソースコード
最後にソースコードを記載しておきます。プログラミングには、Python 3.8を使用しました。アルゴリズム自体をmain関数にしていますが、これは意味がないので関数にしなくても構いません。
import numpy as np
import copy as cp
import math
# To create an animation gif
import matplotlib.pyplot as plt
import matplotlib.animation as ani
from mpl_toolkits.mplot3d import Axes3D
from io import BytesIO
from PIL import Image
# objective function ( Double-well Potential Energy V(x) )
def local_func( x ):
# V(x) = a*x^4 + b*x^2 +c
a = 1
b = -2
c = 1
x2 = np.linalg.norm(x,ord=2)
f = a*x2*x2 + b*x2 + c
return f
def find_best(x):
n = len(x)
ybest = local_func(x[0])
ibest = 0
for i in range(1,n):
ytmp = local_func(x[i])
if ytmp < ybest :
ybest = ytmp
ibest = i
return ibest
def func( x ):
n = len(x)
r = 0
for i in range(0,n):
r += local_func(x[i])
return r
def main():
# loop controllers
maxite = 1000
maxite_local = 100
threshold = 0.85 # under limit of roudness
# scalars
n = 50 # number of bats
dim = 3 # size of dimension
fmin = 0.00 # minimum frequency
fmax = 1.00 # maximum frequency
Amax = 10 # maximum roudness
# list of vectors
x = -1 + ( 1 + 1 ) * np.random.rand(n,dim) # initial bat's positions
v = np.zeros((n,dim)) # initial bat's velocities
r = np.random.rand(n) # initial pulse emission rates
A = Amax*np.random.rand(n) # initial roudness
r0= cp.deepcopy(r)
# constant parameters
alpha = 0.8 # roundness decreasing rate (constant)
gamma = 0.5 # pulse rate control parameter
maxeps = 0.9 # sound amplitude control parameter
# figure
fig = plt.figure(figsize=(8,6))
imgs= []
# Find the best solution
ibest = find_best(x)
xbest = cp.deepcopy(x[ibest])
# Bat Algorithm
xnew = cp.deepcopy(x)
vnew = cp.deepcopy(v)
for t in range(0,maxite):
# Generate new solution by adjusting frequency(xnew)
beta = np.random.rand(n)
f = fmin + ( fmax - fmin ) * beta
for i in range(0,n):
v[i] = v[i] + ( x[i] - xbest ) * f[i]
xnew[i] = x[i] + v[i]
# show logs
Aavg = np.mean(A)
ravg = np.mean(r)
ybest= local_func(xbest)
ytotal = func(x)
rbest = np.linalg.norm(xbest,ord=2)
print(t, ybest, ytotal, rbest, Aavg, ravg)
# create images
fig.clear()
ax = fig.add_subplot(111,projection='3d')
ax.scatter(x[:,0],x[:,1],x[:,2],color='blue')
ax.scatter(xbest[0],xbest[1],xbest[2],color='red')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_xlim(-1.5,1.5)
ax.set_ylim(-1.5,1.5)
ax.set_zlim(-1.5,1.5)
fig.text(0,0,'t={:6}'.format(t))
buf = BytesIO()
fig.savefig(buf)
imgs.append(Image.open(buf))
# stoping criteria
if Aavg < threshold : break
rand = np.random.randn();
vnew = v
for i in range(0,n):
# Perform local search around global best best(xnew)
if rand > r[i]:
# local search
yold = local_func(x[i])
xibest= cp.deepcopy(x[i])
for k in range(0,maxite_local):
epsv = np.random.randn(dim)
leps = np.linalg.norm(epsv,ord=2)
if leps > maxeps : epsv = (maxeps/leps)*epsv
xtmp = cp.deepcopy(x[i]) + epsv*A[i]
ynew = local_func(xtmp)
if ynew < yold :
yold = ynew
xibest = xtmp
xnew[i] = xibest
#print(i,rand,r[i],A[i],yold,xnew)
for i in range(0,n):
if rand < A[i]:
yold = func(x)
ynew = func(xnew)
if ynew < yold:
x[i] = xnew[i]
r[i] = r0[i]*(1-math.exp(-gamma*t))
A[i] = alpha*A[i]
else:
xnew[i] = x[i]
# Update the best solution
ibest = find_best(x)
ybest_old = local_func(xbest)
ybest_new = local_func(x[ibest])
if ybest_new < ybest_old :
#xbest = cp.deepcopy(x[ibest])
xbest = x[ibest]
print(xbest)
imgs[0].save('output.gif',save_all=True,append_images=imgs[1:],duration=100,loop=0)
main()
ちょっとハマった点としては、最適化の途中で急に配列の値が変わってしまうというバグがありました。これは、Pythonの既定コピーが参照コピー（ウィークコピー）だということを自分が知らなかったためでした。いくつかのコピー処理を、値をコピー（ディープコピー）するcopy.deepcopy()関数に変更すれば、問題なく最適化が進むようになりました。
まとめ
群知能バットアルゴリズムを実装しました。
最適解が球形になるテスト関数では、こうもりが球形に集まり、ちゃんと最適化できていることがわかりました。
群知能は、今後何かに使えるような気がしています。
