生物モデルを用いた最短経路問題の解法
M2014SS012 横山 優希
指導教員：小藤 俊幸

1

はじめに

本研究で扱う粘菌には基本的な特性 [1] として, 粘菌自身の体を流動的に運動させて餌を捕食する性質を持っている. さらに, この粘菌には迷路状の地形の入り口と出口に餌を配置すると最短経路を通って迷路を進むと言う性質がある. 本研究では, その性質について述べられた論文[2] の数学モデルから, 数値解析の 1 つとして, 伝導率の更新をオイラー法を用いた近似で数値実験を行った. しかし, オイラー法ではノード数が増加すると明確な解が得られなくなる. 安定性やパスの削除の早さを求めるために γの値に注目をし, 伝導率の更新を指数オイラー法 [3] へと変えて伝導率の正値性やパスの削除の早さを比較している. その結果, 指数オイラー法による近似は安定性のあるものと分かった.

2

粘菌と最短経路問題

2.1

粘菌のもつ性質と有用性粘菌の持つ性質とは何か, それは複数箇所に餌を配置されるとすべての餌を得るために餌から餌へと体を伸縮させるということである. この性質は最短経路問題の解決法の一つとしてあげることができる. 最短経路問題を解くときにもっとも有名とされているものはダイクストラ法である. しかし, この方法では正確な最短経路を識別することは可能であるものの, 複数箇所の目的地つまりノード数が多ければ多いほど必要な計算時間が過剰になってしまう. ここで, 最短経路を見つけるときの条件を三つ述べる.

数学的モデル3

ノードの数 (餌の数) を Ni から Nj まで設置さらに, 流
れは N1 から N2 へと流れていくものとしてポアズイユ流れの法則を用いて仮定すると (1) 式が導かれる.

Qij =

 X
Qij = 0 (j ≠ 1, 2 のとき)

i
 X
Qi1 + I0 = 0iXQi2 − I0 = 0


粘菌を用いた実験方法

ここからは, 具体的に実験方法について述べていく. 粘菌により最短経路を見つける実験方法の一つとして次のものがある. 初めに迷路の中に餌を二ヶ所配置し, そしてそこに粘菌をおく. すると粘菌は様々な経路を通って餌まで辿り着こうとする. しかし, 時間が経つと最短経路を通る粘菌の管以外はなくなり, 一本の道となる. そのことで最短経路を見つけることができると言うわけだ. また,粘菌はある光や危険な化学物質を嫌うため [4] それを粘菌の通り道に配置するとその箇所を避けて最短経路を見つけようとするのである. しかし, この実験には欠点があり餌の量が均一に配置されていなければ多く餌のあるところに粘菌が行ってしまうこと, さらに迷路内での実験だと粘菌が壁を乗り越えてしまう可能性があるというものだ.

(2)iI0 は最初のノードから流れる量である. つまり, 始点と終点は一方的に流れ続け, 流れ込みつづけるため I0 を式に入れることで和を 0 としており, 途中のノードでは流れ込んでくる量, 流れていく量が同じなため総和が 0 になっている. さらに (1) 式 (2) 式を組み合わせると,(3) 式のようにまとめられる.
 −I0 (j = 1)
X Dij(pi − pj ) =I0 (j = 2)
Lij
 0(j ≠ 1, 2)i(3)また,Qij によって Dij が変動したとすると.

- 経路検索の迅速性.

2.2

(1)

この時,pi は Ni の圧力,pj は Nj の圧力を示し,Lij は Ni ,Nj間の長さ (Mij ),Dij は伝導率を示す. また Qij は流量を示している. また, 各ノードに対してキルヒホッフの法則を考慮すると, 以下のような式を得られる.

- 確実に最短経路を見つけるということ.
- 情報の更新からの最短経路の再設定に対する適応性粘菌はこの三つの条件を満たして経路を見つけることができる.

Dij(pi − pj )Lijd
Dij = f (|Qij |) − γDijdt(4)という式が得られ,f (Q) は増加関数であるため,f (0) =0 である. この式から伝導性 Dij が減少する傾向にある.
f (Q) = α|Q| として (3) 式に代入すると.

dDij = α|Qij | − γDijdt(5)

この式は Qij によって変動される Dij を求めるための式となっていて, α,γ は正の定数である.

4粘菌アルゴリズムの実装

4.1

粘菌アルゴリズム概要

最短経路問題を解く際の手順は以下の通りである.
0(a) 初期時刻の伝導率 Dij
を (適当に) 与え,Ni の圧力を基準にする. つまり pi = 0 とする.

(b) 後に示す連立 1 次方程式を解いて, 圧力 p0i を求め, 流量 Q0ij を計算する.

1(c) オイラー近似により, 時刻 t1 での伝導率 Dijを求め
nnる. 以下, 同様に時刻 tn での圧力 pi , 流量 Qij を求め,n+1オイラー近似により, 時刻 tn+1 での伝導率 Dijを求める.
(b) で記した連立 1 次方程式というのは,Qnij =nDij(pn − pnj )Lnij i(6)を以下の式に代入して,
 X nQij = 0 (j ≠ 1, 2 のとき)


i

 X n
Qi1 + I0 = 0
i

X

Qni2 − I0 = 0

図 1 最短経路問題(7)i

それで得られた式 −I0 (j = 1)nn(p−p)=I0 (j = 2)ij
Lnij 0(j ≠ 1, 2)

nX Diji(8)

これを各ノードにおいて p0i の連立 1 次方程式を立てる.
その方程式をガウス消去法で解き p0i を求める. 流量 Q0ijが求められたら, 次の時刻の伝導率を更新するために, 以下のオイラー近似を用いる.
n+1nn= Dij+ δt(α|Qnij | − γDij)
Dij(9)

ここでは α = 1.0,γ = 1.0,δ = 0.1 としている. 次に (4.4)
n+1式で得られた Dijを用いて, 同様に計算をしていく.

4.2

プログラムの実装

今回作成したプログラムは, まずあらかじめ決められたノード間の距離, 伝導率を打ち込み, 読み込ませる. a[i][j]という変数を用いるのだが, この時の a[i][j] というものは,ノード Ni にでの (4.3) における圧力 pj の係数を代入する配列として扱っている. 各ノードにおいて,p[i] の係数の値を代入していくのだが, ノード i と pj が一致する場合 (i =X Dijj). つまり, 対角成分においては a[i][i] を a[i][i] =LijjDijのように代Lij入する. その後 p[i] の初期値をノード 0 を 1.0, それ以外を 0.0 とし a[i][j] と p[i] との連立方程式をガウス消去法で解き,(4.1) 式を用いて Q を求める. その後, 伝導率 D の更新を (4.4) 式を用いて求める. この手順を 100 回ループさせたものがこのプログラムとなる.

9999
Q99

01 = 0.003378, Q04 = 0.719808, Q05 = 0.276814


9999

Q12 = 0.003140, Q14 = 0.000238


9999

Q99

23 = 0.120676, Q24 = −0.237229, Q29 = 0.119693


99

 Q39 = 0.120676
99
Q99
45 = 0.001184, Q46 = 0.481633
99

Q56 = 0.004722, Q99

57 = 0.273276

9999

Q67 = −0.000010, Q99

68 = 0.455848, Q69 = 0.030517



99

Q78 = 0.273266

 99
Q89 = 0.729114

この時の Q99
ij が相対的に大きければ大きいほど粘菌がNi ,Nj 間を最短経路とみなして移動するということである. つまりこの最短経路問題は, 0 → 4 → 6 → 8 → 9 が最短経路であると言える.

プログラムの正確性

6

オイラー法を用いたプログラムではノード数が少なければ最短経路が出る. しかし, ノード数を増やし, 乱数で伝導率や距離を与えると明確に最短経路が求められない.
そこで伝導率の更新をパラメーター γ に注目をして, 指数オイラー法を用いることで, 計算の高速化をはかる.

7プログラムの改良

7.1

指数オイラー法

プログラムの改良に当たって指数オイラー法を用いる.(5) 式から α|Qij | = g(t),Dij = D(t) と置く.

それ以外の非対角成分の時は a[i][j] = −5dD
= g(t) − γD(t)dt

このとき D(t) = C(t)e−γt と式をおいた際,D(t) を微分
すると D′ (t) = C ′ (t)e−γt − γC(t)e−γt この時,C ′ (t) =g(t)eγt ならば (10) 式を満たす. ここから実験結果

図 1 のような最短経路問題があったとする. この問題に0おいては Dijを適当に与えて, 先ほど説明したプログラム
で解くと次のような値が得られる(10)

Z tn+1g(s)eγs ds−γtn+1Z tn+1C(tn+1 ) = C(tn ) +(11)tnとなる.
−γδtD(tn+1 ) = D(tn )e+etng(s)eγs ds(12)ここで被積分関数の g(s) を定数関数 g(s) ≡ g(tn ) で近似する. また,−γtn+1e

Z tn+1eγs ds=tn1(1 − e−γδt )γ(13)であるため,(12) 式に入れて, 近似値で置き換え, 元の値に戻すとDn+1 = Dn e−γδt +α|Qij |(1 − e−γδt )γ(14)という式が得られる. プログラムへの実装はオイラー法で行った (9) 式の代わりに,(14) 式を用いることでオイラー法から指数オイラー法への変更としている.

7.2

指数オイラー法のメリット

伝道率の更新を行う際に γ = 1.0 というパラメーターを用いた. この時 γ の値は不要なノード間のパスが消えるスピードとなる. このパラメーターの値を大きくすればするほど早く計算ができる. しかし, オイラー法を用いたプログラムでは, ある一定の γ の数値以上の値をとってしまうと D の値に負が出てしまい発散してしまう. 指数オイラー法を用いることで,γ の値をどれだけ大きくとってもD の値が常に正であることが (14) 式から分かる. つまり,オイラー法から指数オイラー法へとプログラムを改良するメリットは以下の通りであると考えられる.

(1) Dij の正値性が保たれる.
(2) 不要なノード間のパスを早く消すことができる.
この 2 点のメリットを実証するためにオイラー法と指数オイラー法との比較を行っていく.

8オイラー法と指数オイラー法の比較

8.1

安定性の比較

オイラー法を用いた際 γ = 16.0 より大きな値をとると,Dij の値に負の値が出てしまい, 発散してしまう. ここからは,γ = 16.0 の時のオイラー法,γ = 17.0 の時のオイラー法の,γ = 17.0 の時の指数オイラー法の 3 種類のパターンで図 1 と同じ例題を解き比較を行っていく. 最初はオイラー法を用いて γ = 16.0 の時の実験結果である.
999999
 Q01 = 0.000000, Q04 = 0.999999, Q05 = 0.000000

9999

 Q12 = 0.000000, Q14 = 0.000000

9999

 Q99
23 = 0.000001, Q24 = −0.000001, Q29 = 0.000001
99
 Q39 = 0.000001
99
Q99
45 = 0.000000, Q46 = 0.999998

99

Q99

56 = 0.000000, Q57 = 0.000000


99
9999


 Q67 = 0.000000, Q68 = 0.062500, Q69 = 0.000000
99
Q78 = 0.000001
 99
Q89 = 0.999999

このように,γ の値を大きくすることでノード間のパスが消えていることが分かる. さらに, 最短経路の部分の Q の値が大きくなりどのノードを通ると最短経路かがわかりやすくなる. また,D の値も Q と同様にパスがないところは 0 になっている. 次に γ = 17.0 の値をとったときのオイラー法での結果である.
99
9999

 Q01 = 0.752739, Q04 = 0.247261, Q05 = 0.000000

9999

 Q12 = −0.001698, Q14 = 0.754437

9999

 Q99
23 = 0.008557, Q24 = −0.027252, Q29 = 0.016997
99
 Q39 = 0.008557
99
Q99
45 = −0.025754, Q46 = 1.000200

99

Q99

56 = −0.008085, Q57 = −0.017669


9999

Q67 = 0.108650, Q68 = 0.883465, Q99

69 = 0.000000
99
=0.090981
Q

78

 99
Q89 = 0.974446

また,γ = 16.0 の時と違い D の値が一定でなくなり, 以下の通りになる.
99
9999

 D01 = −0.072053, D04 = −0.008620, D05 = 0.000000

99
99

 D12 = −0.000115, D14 = −0.072624

99
99
99

= 0.000284, D24= 0.000498, D29= 0.000563
 D23



99
 D39 = 0.008085
99
99
D45= −0.000540, D46= 0.013839

99
99

D56 = −0.000277, D57 = −0.000295

999999D
67 = −0.037545, D68 = 0.050403, D69 = 0.000000
99
D78= 0.004628

 99
D89 = 0.057850

これより γ = 16.0 を越えると実験結果と最短経路の解が合わなくなり, 安定性にかけているということが分かる. 次にオイラー法で解が出なくなったパラメーター γ = 17.0と同じ値をとり, 指数オイラー法のプログラムで解を求めると

99
99Q99

01 = 0.000000, Q04 = 0.998876, Q05 = 0.001124



99

Q99

12 = 0.000000, Q14 = 0.000000

99
99

Q23 = 0.000132, Q99

24 = −0.000263, Q29 = 0.000131


99

 Q39 = 0.000132
99
Q99
45 = 0.000000, Q46 = 0.998613
99 Q56 = 0.000000, Q99

57 = 0.001124
999999
=−0.000000,
QQ
68 = 0.998613, Q69 = 0.000000
67




Q99

78 = 0.001124
 99
Q89 = 0.999737
オイラー法との比較や正値性が保たれているかの確認のために伝導率 Dij も示す.
99
99
99
= 0.058751, D05= 0.000072= 0.000000, D04D01




9999

D12 = 0.000000, D14 = 0.000000



99
9999

D

23 = 0.000009, D24 = 0.000017, D29 = 0.000009
99

 D39 = 0.000009
99
99= 0.058734
D45 = 0.000000, D46
9999

D=0.000000,
D=0.000072
5657
999999= 0.000000
 D67 = 0.000000, D68 = 0.058734, D69
99
 D78 = 0.000072

 99
D89 = 0.058806

γ = 17.0 の時のオイラー法と指数オイラー法の実験結果の値を比較すると D の値 Q の値共に指数オイラー法の方がより良い値をとっていることが分かる. 指数オイラー法の時の Dij の値に負の数が出ずに正値性が保たれていることがわかり, 安定性があるといえる.

10.8

死滅のスピードの比較

0.6
D8.2

γ = 16.0
γ = 100.0

次に, 死滅のスピードの比較を行う. 図 2 はオイラー法
で計算可能な値 γ = 1.0 と γ = 16.0 を比較した際, 伝導率 D56 に注目をしその減少のスピードをグラフ化したものである. x 軸は繰り返しの回数 t の値をとり,y 軸にはその時の D56 の値をとっている. オイラー法を用いた計0.40.2000.5

11.522.533.54t4

γ = 1.0
γ = 16.0

図 4 γ = 16.0 とγ = 100.0

3.539D

2.521.510.5002468101214161820t

図 2 γ = 1.0 とγ = 16.0
算だけでも, ノード間のパスの死滅の早さは目に見えて分かる. 図 3 ではオイラー法で計算可能な γ = 16.0 と計算不可能な γ = 17.0 のパラメーターをとった時のグラフを示している. 図 2 でも D56 を選択したが, 理由としては
γ = 16.0,γ = 17.0 どちらの値をとっても手順を 100 回繰り返せば 0 へとなる点からである. 図 3 より, どちらも約0.5

γ = 16.0
γ = 17.00.4

D0.30.20.100246810t

図 3 γ = 16.0 とγ = 17.0

8 回目の繰り返しでパスの削除がされている. 回数はどちらも大差が無いように見えるが,γ の値を 1 変えるだけで,死滅のスピードが微量であるが早くなっていることが分かる. さらにオイラー法で計算可能な範囲 γ = 16.0 と,指数オイラー法で計算可能な値をわかりやすく比較したものが, 図 4 のグラフである.γ = 100.0 まで数値を上げる
ことで, プログラムの手順を 1 度繰り返しただけで, 伝導率が 0 へと急激に近づいている. 図 3 ではオイラー法と指数オイラー法の違いが微量であったが, 指数オイラー法はここまで, さらにはこれ以上のスピードを出すことができる.

おわりに我々が目にするであろうカーナビ, 最短経路のルートを見つけられるだけではなく,「事故による情報の更新といった条件」も粘菌を使った最短経路問題の解法はクリアできる. さらに, 光や, 化学物質を避けるといった性質を利用す
れば鉄道の線路を効率よく敷くための最短経路を導き出せる. 日本地図状の山や川など線路を敷けない地形に光や化学物質を配置し, 駅をおきたい場所に餌を並べておけば粘菌が効率の良い線路の敷き方を示してくれるといった具合である. ただし, 実地での調査を行わなければ山道など急激なカーブや劣悪な地面の状況などが分からないため, 大きな事故を引き起こす可能性があるという欠点もある. また, この数学モデルからオイラー法を学びプログラムを組み実験を行うことで, プログラムを用いて最短経路問題を解くことができるということが分かり粘菌アルゴリズムに対して理解が深まった. それだけではなく, 実験結果を出す早さや安定性を追求するために, 指数オイラー法を学び先に組んだプログラムを改良することで, より安定性のあるプログラムにすることができた. 今後の課題としてはダイクストラ法での解き方をプログラムに組み込み, その結果を一致するか自動で判別するプログラムを組むことがあげられる. また, 最短経路が 2 パターン得られる問題と直面した際, どのような基準で経路が決定されるか調べることが課題である.

参考文献[1] 中垣俊之.『粘菌 その驚くべき知性』, 株式会社 PHP研究所,2010.
[2] A. Tero,R. Kobayashi,T. Nakagaki, ”Physarum
solver: A biologically inspired method of road- network navigationb” Phisiva A 363 (2006), 115-119.
[3] M. Hochbruck, A. Ostermann, Exponential integrators, Acta Numerica 19 (2010), 209-286[4] 築 地 克 弥, 松 田 千 夏, 高 橋 和 成.『 忌 避 化 学 物 質を 利 用 し た 粘 菌 変 形 体 の ネット ワ ー ク 形 成 』.
http://www.ous.ac.jp/garden/kenkyuhoukoku/
16/Naturalistae16-39-45.pdf