2011-12-20

駆け足で読む『Swarm Intelligence; From Natural to Artificial Systems』Chapter 2 アリの食糧探索行動、組み合わせ最適化、意思伝達ネットワークによる経路決定

２経路の選択
- ２つの要素
- ２択の場合
- $P_A=\frac{Ka}{Ka+Kb}=1-P_B;Ka=(S+a)^n,Kb=(S+b)^n$ 、ここで $a,b$ は選択肢A,Bのそれぞれに、「実績として起きた事象の数(アリで言えば、アリが落としたフェロモン量)」
- $n$ は非線形の程度を決める。 $S$ はフェロモンがないときの「得点」→ランダム性のための下駄
複数の選択肢でどれかが選ばれていく様子

library(MCMCpack)
Niter<-1000

n<-4
k<-20

A<-4

M<-matrix(0,Niter,A)

for(i in 2:Niter){
	M[i,]<-M[i-1,]
	prob<-(M[i-1,]+k)^n
	prob<-prob/sum(prob)
	s<-sample(1:A,1,prob=prob)
	M[i,s]<-M[i,s]+1
}

matplot(M/((1:Niter)-1),type="l")

StarLogoというソフト(こちら)
平面の探索
- 平面を正方格子で考えると、それは、正方格子というグラフを空間とした行動であって、「グラフ」を空間にすることで一般化できるだろう
- 以下のソースは次のように考えている
- 有限平面を正方格子として、それを表すグラフについて、隣接ノードを列挙する
- アリは、このグラフ上を確率的に移動する
- アリは、複数の状態を持つ(エサを巣に持ち帰るときには、エサを持っている状態と持っていない状態の２状態)
- アリは、自身の状態によって、進行エッジの選択に関して異なる確率ベクトルを持つ
- アリは、複数の進行エッジの先のノードのフェロモン量に応じて、エッジを選ぶが、かなりの乱雑さを持たせることが、集団としての探索がうまく行くコツのようである(この乱雑さは、フェロモンの総量によらずに下駄をはかせる方がよさそうである(どんなにフェロモンが多くても、一定の確率で、フェロモンのない格子点へも進めるようにする)
- アリは自身の状態により、異なるフェロモンを現在地に置く
- アリの置くフェロモンは、特定の状態のアリのエッジ選択情報となる
- フェロモンは時間とともに消失する
- 巣とエサの位置は指定する(複数の格子点を一塊として指定するのがよい)
- エサはアリに持ち帰らせる(結果としていつか消滅する)ようにしてもよいし無尽蔵にしもよい(以下のソースは無尽蔵)
- 巣とエサとは、生物学的な意味であるが、これは、アリの状態を変更する場所であると定義すればよい
解く課題は
- 経路の最適化・分布の最適化
- セールスマン問題
- ２次割り付け問題(Quadratic assignment problem(Wikiはこちら))は(こちらの用量反応曲線に関して用いる、Isotonic regression(こちら)の２次計画問題と同様に最適化の困難問題)

# 正方格子
# 格子の１辺の点の数
N<-100
# 全格子点の２次元座標
XY<-expand.grid(1:N,1:N)
# 格子点の数
Nloc<-length(XY[,1])
# すべての格子点について隣接点を全列挙する
Neighbors<-list()

difX<-outer(XY[,1],XY[,1],"-")
difY<-outer(XY[,2],XY[,2],"-")
for(i in 1:Nloc){
	Neighbors[[i]]<-which(difX[i,]^2+difY[i,]^2<=1)
	Neighbors[[i]]<-Neighbors[[i]][which(Neighbors[[i]]!=i)]
	
	#Neighbors[[i]]<-Neighbors[[i]][sample(1:length(Neighbors[[i]]),2)]
	
}
# アリの数
Nant<-300
# アリの状態数
Nstate<-2
# 初期状態で、すべてのアリは１か所(巣の中)にいることにする
half<-floor(N^2/2+N/2)
Ax<-sample(half:(half+1),Nant,replace=TRUE)
# アリの初期状態も１つ(エサを探しに行く状態)
Sx<-sample(1:1,Nant,replace=TRUE)
# 状態別のアリ分布ベクトル
V<-matrix(0,Nloc,Nstate)

for(i in 1:Nant){
	V[Ax[i],Sx[i]]<-V[Ax[i],Sx[i]]+1
}

# 状態変化指定行列
# ある場所は、ある状態をある状態に変える。それを表す行列
C<-matrix(rep(1:Nstate,Nloc),byrow=TRUE,ncol=Nstate)
# 巣の中心とその半径
Center.1to2<-c(N*2/9,N*2/9)
Radius.1to2<-N/20
# エサのある場所の中心とその半径
Center.2to1<-c(N/2,N/2)
Radius.2to1<-N/20
for(i in 1:Nloc){
	if(sum(XY[i,]-Center.1to2)^2<=Radius.1to2){
		C[i,1]<-2
	}
	if(sum(XY[i,]-Center.2to1)^2<=Radius.2to1){
		C[i,2]<-1
	}
}

# 状態別フェロモン分布ベクトル
P<-matrix(0,Nloc,Nstate)
# 状態別フェロモン下駄値
Pgeta<-rep(0.5,Nstate)
# 状態別フェロモン指数
Pexp<-rep(5,Nstate)
# 状態別フェロモン投下量
Pv<-c(1,1)
# 状態別フェロモン投下対象状態
Pt<-c(2,1)
# 単位時間当たりの減衰率
Pr<-c(0.95,0.95)
# 試行時間
Niter<-10000
par(mfcol=c(2,3))
for(i in 1:Niter){
	# 移動
	tmpV<-matrix(0,Nloc,Nstate)
	for(j in 1:Nstate){
		for(k in 1:Nloc){
			destination<-Neighbors[[k]]
			#prob<-(P[destination,j]+Pgeta[j])^Pexp[j]
			prob<-(P[destination,j])^Pexp[j]
			if(sum(prob)==0){
				prob<-rep(1,length(prob))
			}
			prob<-prob/sum(prob)
			prob<-prob+Pgeta[j]
			prob<-prob/sum(prob)
			s<-sample(destination,V[k,j],replace=TRUE,prob=prob)
			for(l in s){
				tmpV[l,j]<-tmpV[l,j]+1
			}
		}
	}
	# 状態変化
	tmpV2<-matrix(0,Nloc,Nstate)
	for(j in 1:Nstate){
		for(k in 1:Nstate){
			tmp<-which(C[,j]==k)
			tmpV2[tmp,k]<-tmpV2[tmp,k]+tmpV[tmp,j]
		}
		
	}
	V<-tmpV2
	#print(apply(V,2,sum))
	# フェロモン
	# フェロモン減衰
	for(j in 1:Nstate){
		P[,j]<-P[,j]*Pr[j]
	}

	for(j in 1:Nstate){
		P[,Pt[j]]<-P[,Pt[j]]+V[,j]*Pv[j]
	}
	
	#par(ask=FALSE)
	image(matrix(log(P[,1]+1),N,N),main="pheromone1")
	image(matrix(log(P[,2]+1),N,N),main="pheromone2")
	image(matrix(apply(log(P+1),1,sum),N,N),main="pheromone1+2")
	image(matrix(sign(V[,1]),N,N),main="ant1")
	image(matrix(sign(V[,2]),N,N),main="ant2")
	image(matrix(apply(sign(V),1,sum),N,N),main="ant1+2")
	#Sys.sleep(0.1)
	#par(ask=TRUE)
}