boids algorithm 実装方法についての簡単なメモ

D3.js アルゴリズム

boids algorithm とは、群れをシミュレーションするマルチエージェントプログラム。各エージェントは、次の3つのルールで動作する。

  • 分離（Separation）エージェントが他のエージェントとぶつからないように距離をとる。
  • 整列（Alignment） エージェントが他のエージェントと概ね同じ方向に飛ぶように合わせる。
  • 結合（Cohesion) エージェントが他のエージェントが集まっている群れの中心方向へ向かうように方向を変える。

これだけのルールで、魚や鳥の群れをよくシミュレーションできるといわれていて、映画などのCGにもよく利用されているとのこと。 各エージェントはシンプルなルールながら、全体としては複雑に振る舞うことから、複雑系、自己組織化的な現象が垣間見れるといわれている。

 

いろいろな実装が考えられるようだが、特にシンプルな、各エージェントの速度が一定で、2次元空間の進行角度が変化する場合の実装について、"参考"の情報を参考にしながら、javascript & D3.jsで表現した。

事前準備

 まず、各エージェント(node)は、場所(x,y)と進行方向の角度(deg)、さらにエージェントの行動に影響を与える、距離の近い他のエージェントのリスト(neighbors)で構成する。

node =[ {x:0,y:0,deg:0,neighbors=[] },... ];

でもって、neighborsを作成するため、総当りで各エージェント間の距離を調べる。$N$個のエージェントなら、${}_n \mathrm{C}_2 = N(N-1)/2$を調べることになる。ただし、neighborsにリスト化するのは、そのエージェントが見える範囲、ということでパラメータvisionを距離の値として、それより近い距離にあるエージェントをneighborsリストに追加する。

for (i=0;i<node.length;i++) {
    for (j = i+1;j<node.length;j++) {
        d = (※node[i]とnode[j]の距離);
        if (d < node[i]のvision )
            node[i].neighbors.push({dist:d,deg:node[j].deg,x:node[j].x, y:node[j].y});
        if (d < node[j]のvision )
            node[j].neighbors.push({dist:d,deg:node[i].deg,x:node[i].x, y:node[i].y});
    }
}

node[i].neighborsのリストはあらかじめ大きい順にsortしておく。

    node[i].neighbors.sort(function(a,b){return a.dist-b.dist;});

分離の実装

さて、ここからが本題。

"分離"のルールは、「他のエージェントと近づき過ぎないように距離をとる」というもの。パラメータminimum-separationを距離の値として、エージェントのリストneighborsのうちもっとも近いエージェントがこの距離より小さかったら、離れるように動くことにする。

今回の実装では、各エージェントはスピード一定で、変更できるのは進行角度だから、この近づきすぎてしまったエージェントの進行方向と逆方向に自分の進行角度を向けば離れられる。つまり、

  • (今の進行角度 ) → (最近傍のエージェントの進行角度) + 180°
  • $\delta \theta$ = (最近傍のエージェントの進行角度) + 180° - (今の進行角度)

微少時間の動作としては、この$\delta \theta$を今の進行角度との変位として加算して推移するようにする。

ところが、この右側の角度の差が、若干厄介である。例えば、260° - 100°=160°である。では、-100°-100°=-200°か？ 260°=-100°だから同じ結果になってくれないと困るし、方向転換をするための計算なのだから、次の図のように、なるべく回りやすい方の角度(160°)を出力としたい。

レイヤ 1 a=260° b=100° a'=-100° 160° -200°

つまり、角度aと角度bの差(a-b)は、角度bを足して角度aと等しくなる、もっとも絶対値が小さいほう(つまり -180～180の間)であってほしい。(a-b)を横軸において、縦軸が角度の差の出力としたグラフは、

Layer 1 (a-b) 角度の差 360° 180° -360° -180°

 となるような実装とある。この引き算の仕方は、別に角度に限ったわけではなく、同じような周期境界条件をもつ場合(有限加群)に適用できる。360°を周期periodに変えればよいわけで、このグラフをコードにすると、次のようになる。

function diffCyclic(deg1,deg2,period) {
	var deg = (deg1 - deg2) % period;
	if ( Math.abs(deg) <= period/2 ) return deg;
	else {
		if ( deg > 0 ) return -period + deg;
		else return period + deg;
	}
}

この関数を使って、

\[ \delta \theta = \textrm{diffCyclic} ((最近傍のエージェントの進行角度) + 180° , (今の進行角度),360 ) \]

と書けば進行角度の変位$\delta \theta$が求められる。

※今回の場合はx,y空間も周期的な境界(つまりトーラス空間)を持つことにした。なので、事前準備で、各エージェントの距離を算出するときもこの差の関数diffCyclicを利用する。

整列の実装

次に、整列の実装である。"整列"は、「周辺のエージェントの角度と同じ方向にあわせる」わけなので、リストneighborsの各エージェントの進行角度の平均を取得して、その差分を変位とすればよい。

ところがこれまた厄介で、角度の平均というのは一般に難しい。例えば、10°と20°の平均は(10+20)/2 = 15°でよいが、345°と5°の平均は、(345+5)/2=180°となってしまう。

そこで、いったん角度を半径1の円上の点の(x,y)に展開して、xとyについて合計を取り、最後に(sum x,sum y)に対して$\tan^{-1}$をとるようにして、求めたい、各エージェントの進行角度の平均を取得する。取得した平均進行角度にあわせるように変位$\delta \theta$を決める。これを三角法と言うらしい。

var ax = 0,ay = 0;
for (j=0;j<node[i].neighbors.length;j++) {
    ax += Math.cos(node[i].neighbors[j].deg/degrees); // degrees = 180 / Math.PIで割戻し
    ay += Math.sin(node[i].neighbors[j].deg/degrees);
}
da = diffCyclic(Math.atan2(ay,ax)*degrees,node[i].deg,360); //例の差の関数を使って角度の差分を求める。

結合の実装

最後に結合の実装。"結合"は「群れの中心方向へ向かうように方向を変える」ので、neighborsの平均位置の方へ向くように進行角度を変更すればよい。

x,y空間は周期境界条件をつけていた。

境界をはさんだエージェントの位置の平均をそのまま計算してしまうと、えらい間違いがおこってしまう。ここでは、先ほどの差の関数diffCyclicを使って、自分自身のエージェントとの相対位置の平均を取ることにすればよい。

var cx = 0,cy = 0;
for (j=0;j<node[i].neighbors.length;j++) {
	cx += diffCyclic(node[i].neighbors[j].x,node[i].x,width); // 自分エージェントとの相対位置を差の関数で得る
	cy += diffCyclic(node[i].neighbors[j].y,node[i].y,height);
}
dc = diffCyclic(Math.atan2(cy,cx)*degrees,node[i].deg,360);  // cohere rule

相対位置の合計(cx,cy)から、$\tan^{-1}$で角度を求めている。

パラメータ

さて、最後に、それぞれのルールに対する強度のパラメータを定義することで、パラメータの差異によるエージェントの振る舞いの違いを検討できるようにする。

  • max-separate-turn... "分離(Separation)"が発生したときの変更される進行角度の最大絶対値。これが大きいと"分離"しやすくなる。
  • max-align-turn... "整列（Alignment)"ルールで変更される進行角度の最大絶対値。これが大きいと"整列"しやすくなる。
  • max-cohere-turn... "結合（Cohesion)"ルールで変更される進行角度の最大絶対値。これが大きいと"結合"しやすくなる。

それぞれのルールで算出した$\delta \theta$の絶対値を上記のパラメータで比較して、変位が上記のパラメータより大きくならないように制限する。

...
deltad = compDegree(max-separate-turn,deltad);
...
function compDegree(maxd,deg) {
	if (Math.abs(deg) < maxd) return deg;
	else {
		if (deg > 0 ) return maxd;
		else return -maxd;
	}
}

以上の実装をD3.jsを加味したり省略整理などしたものが以下。

こんな感じ

そのまんまだと面白みが無いので、ここではエージェントにグループを複数設定できて、グループごとにエージェントの数と各種パラメータを設定できるようにした。デフォルトの設定だとグループ0のminimum-separationが小さくて、だんだんグループ0が固まりとなってクラスタ化していき、それに対してグループ1,2は比較的ばらける様子が伺える。(ちゃんと定量的な評価が必要なんでしょうね)

参考

複雑系科学で有名なサンタフェ研究所 (Santa Fe Institude)が主催しているComplexity Explorer というHPがある。ここでは、複雑系科学に関する free online course を受講できる。

Complexity Explorer

実装は講義 "Introduction To Complexity"で紹介されていたboids algorithmを参考にした。NetLogoでのboidsライブラリも紹介されていて、パラメータの設定はその実装を参考にした。(アルゴリズムも実はほぼ同じ)

NetLogo Home Page