[VBA] 粘菌シミュレータで鉄道網をつくる

本内容は @STInverSpinel 様による記事を参考にさせていただきました

素敵で面白く、ためになる記事をありがとうございます
また、本記事はこの記事に書かれた概念をベースに作成しているため、前もってこの記事を読んでいただくことを強くオススメいたします

鉄道網設計問題

今年も日本人科学者の方がイグノーベル賞を受賞されましたが、その昔粘菌が構築するネットワークが迷路探索や鉄道網の設計に役立つのではないかという研究がイグノーベル賞を受賞し話題になったことを覚えていらっしゃる方も多いかと思います。

手老篤史 「粘菌の輸送ネットワークから都市構造の設計理論を構築」

特にこの「鉄道網の設計」という内容に興味を持った私は、以下の問題を解決する同様のプログラムを自作したいと考えました。

---

問

$n$個の都市それぞれを鉄道網でつなぐ。各都市の座標と人口だけが与えられているとき最適な鉄道網をデザインできるか。

表 都市データの例

---

この「最適な鉄道網」という言葉の定義に関してまだまだ考えがあいまいなのですが、今のところは「建設コストがそれなりに安く、全ての住民が人口の多い大都市ほどそれなりに短い時間でアクセスできる鉄道網」くらいに考えています。そしてこれは平面に散らばった餌を限られた体積で効率よく輸送するネットワークを作ることのできる粘菌の行動とよく似ています。

本問題に対し、今回作ったプログラムを実行すると以下のような結果が得られました。

（左：人口に応じてサイズを変えてプロットした都市、右：8400ステップ実行結果）

本プログラムからは、まずは地図の左半分にある主要４都市を結ぶ「太い（例えば列車本数の多い）」環状線を作り、残りの都市を枝分かれ路線でつないでいく鉄道網が提案されました。

これが本当に最適な鉄道網なのか、実用に堪えうるものかについては先人達によって盛んに議論されてきた¹ことなのですが、私の不勉強のため、ここではそれらについて考慮できていないことをご容赦ください。

使用する数理モデル

本記事のシミュレータモデルは基本的に @STInverSpinel 様の記事（以降「参考記事」と呼称させてください）の設定を引き継いでいるのですが、その内容についてここでも簡単に説明します。

本プログラムでは任意の都市（餌）$i,j$が$D_{ij}$の太さの線路（輸送管）でつながれています。都市$i-j$間の距離を$L_{ij}$、そこを流れる（人）流量を$Q_{ij}$、そしてそれぞれの都市における圧力を$p_{i}, p_{j}$とおくと、粘菌における物理モデルより以下のように表せるとされています（詳しくは参考記事または参考文献に示す元論文をご参照ください）。

Q_{ij}=\frac{D_{ij}}{L_{ij}}(p_{i}-p_{j})

この式から、（人）流量$Q_{ij}$は都市$i, j$における圧力差に比例しており、線路太さ（輸送管太さ）$D_{ij}$が大きいほど流れやすくなり、また都市（餌）間距離$L_{ij}$が長くなればその分輸送量は少なくなることが読み取れます。

方程式の作成

たとえば都市1を出発（流出；湧出）都市に、都市2を到着（流入；吸込）都市に設定し、そこに$I_{0}$の流れがあったとき、キルヒホッフの法則から（人）流量$Q_{ij}$に関して以下の関係が成り立ちます。

\sum_{j} Q_{ij} = \left\{
\begin{array}{ll}
I_{0} & (i=1) \\
-I_{0} & (i=2) \\
0 & (i \neq 1,2) \\
\end{array}
\right.

この式を理解するために簡単な例を考えます。

先ほどの設定に加え、二つの都市1,2が都市3ともつながっており、都市1から都市2に$\frac{3}{4} I_{0}$、都市3を経由するルートにも$\frac{1}{4} I_{0}$の流れがあったとします。このときの流量を整理すると以下のようになります。

\displaylines{
Q_{12}+Q_{13} = \frac{3}{4} I_{0} + \frac{1}{4} I_{0} = I_{0} \\
Q_{21}+Q_{23} = -\frac{3}{4} I_{0} - \frac{1}{4} I_{0} = -I_{0} \\
Q_{31}+Q_{32} = -\frac{1}{4} I_{0} + \frac{1}{4} I_{0} = 0
}

このことから、出発都市からは$I_{0}$の流出、到着都市には$I_{0}$の流入、そしてそれ以外の都市には流れがあったとしても総和は$0$になるというのがキルヒホッフの式であり、先ほどの式はそれを一般化したものであることがわかります。

ここで出発都市を1、到着都市を2、（人）流量を$I_{0}$としたまま、都市が全部で$n$個ある場合を想定すると、以下のような式が成り立ちます。ここで、同じ都市間の流量はゼロ、すなわち$Q_{ii}=0$であることに注意してください。

\displaylines{
0+\frac{D_{12}}{L_{12}} (p_{1}-p_{2}) + \frac{D_{13}}{L_{13}} (p_{1}-p_{3}) + \cdots + \frac{D_{1n}}{L_{1n}} (p_{1}-p_{n}) = I_{0} \\
\frac{D_{21}}{L_{21}} (p_{2}-p_{1})+ 0 + \frac{D_{23}}{L_{23}} (p_{2}-p_{3}) + \cdots + \frac{D_{1n}}{L_{1n}} (p_{2}-p_{n}) = -I_{0} \\
\frac{D_{31}}{L_{31}} (p_{3}-p_{1}) + \frac{D_{32}}{L_{32}} (p_{3}-p_{2}) + \cdots + \frac{D_{3n}}{L_{3n}} (p_{3}-p_{n}) = 0 \\
\frac{D_{41}}{L_{41}} (p_{4}-p_{1}) + \frac{D_{42}}{L_{42}} (p_{4}-p_{2}) + \cdots + \frac{D_{3n}}{L_{3n}} (p_{4}-p_{n}) = 0 \\
\vdots
}

この関係は行列を用いて以下のように表示されます。

\begin{pmatrix}
0+\frac{D_{12}}{L_{12}}+\frac{D_{13}}{L_{13}}+ \cdots +\frac{D_{1n}}{L_{1n}} & -\frac{D_{12}}{L_{12}} & \cdots & -\frac{D_{1n}}{L_{1n}} \\
-\frac{D_{21}}{L_{21}} & \frac{D_{21}}{L_{21}}+0+\frac{D_{23}}{L_{23}}+ \cdots +\frac{D_{2n}}{L_{2n}} & \cdots & -\frac{D_{2n}}{L_{2n}} \\
\vdots & \vdots & \ddots & \vdots \\
-\frac{D_{n1}}{L_{n1}} & -\frac{D_{n2}}{L_{n2}} & \cdots & \frac{D_{n1}}{L_{n1}}+\frac{D_{n2}}{L_{n2}}+ \cdots +\frac{D_{n(n-1)}}{L_{n(n-1)}}+0 \\
\end{pmatrix}
\begin{pmatrix}
p_{1}\\
p_{2}\\
\vdots\\
p_{n}\\
\end{pmatrix}
=
\begin{pmatrix}
I_{0}\\
-I_{0}\\
\vdots\\
0\\
\end{pmatrix}

以降は上の式を以下の行列で簡潔に表すことにします。

\boldsymbol{A}\boldsymbol{p}=\boldsymbol{b}

また、都市$i, j$をつなぐ線路（輸送管）は流れの向きによらず同じであることから

\displaylines{
D_{ij}=D_{ji}\\
L_{ij}=L_{ji}
}

であり、$D_{ij}, L_{ij}$を要素に持つ行列$\boldsymbol{D},\boldsymbol{L}$および$\boldsymbol{A}$はすべて対称行列になっています（そのことを私はあまり活用できていないのですが、もしかしたら計算コストを減らせるかもしれません）。

プログラムの設計

本プログラムの大まかな設計として、１つのステップごとに粘菌を模した流れ現象が発生し、その結果に基づいて線路太さ（輸送管太さ）$D_{ij}$がわずかに変わるようにして、十分なステップ数を繰り返せば目的の形に近づくようなプログラムを考えます。

人口に応じた出発、到着都市の選択

本プログラムでは、１ステップにおける出発都市、到着都市を人口に応じた確率で選択するようにしています。下に示すコードは一見ややこしそうですが、やっていることは袋の中に都市の名前が書かれた球をその都市の人口分ずつ入れていった後、その中から無作為に選択した球を当ステップにおける出発or到着都市と決めていくような作業です。これにより、都市の人口が多ければ出発or到着都市として選ばれる確率が高くなるような選択が行われます²。

Do While (RouFr >= RouTo)
    'TotalPopu: 各都市の人口を合計した値
    'an, bn: TotalPopu以下で無作為に選ばれる自然数
    an = Int(TotalPopu * Rnd + 1)
    bn = Int(TotalPopu * Rnd + 1)

    'popl(k): k番目の都市における人口 

    '出発都市の選出
    RouFr = 1
    Do While (an > 0)
        an = an - popl(RouFr)
        RouFr = RouFr + 1
    Loop

    '到着都市の選出
    RouTo = 1
    Do While (bn > 0)
        bn = bn - popl(RouTo)
        RouTo = RouTo + 1
    Loop

    '計算上の操作
    RouFr = RouFr - 1
    RouTo = RouTo - 1
Loop

こうして選ばれた都市をもとに一次方程式$\boldsymbol{A}\boldsymbol{p}=\boldsymbol{b}$の右辺$\boldsymbol{b}$について、選出された都市番号に対応する行の流入量、流出量が設定されたベクトルを作成することができます。

\boldsymbol{b} = \begin{pmatrix}
\vdots\\
I_{0}\\
\vdots\\
0\\
0\\
\vdots\\
-I_{0}\\
\vdots\\
\end{pmatrix}

連立一次方程式の計算

それぞれの都市をつなぐ線路の太さ（輸送管太さ）については、初期状態で$D_{ij}=1$などと適当な値で設定し、それ以降のステップでは直前のステップで得られた$D_{ij}$を引き継ぐようにします。こうして$D_{ij}$が決まったことで行列$\boldsymbol{A}$の中身が決まります。直前の節で当ステップにおける$\boldsymbol{b}$も決めたので、いよいよ連立一次方程式$\boldsymbol{A}\boldsymbol{p}=\boldsymbol{b}$を解く作業に入ることができます。ただしこの問題は理工系大学生ならプログラミングの練習課題として解かされるほど有名な問題なので、ここではその説明を割愛させていただきます³。
この作業により、当ステップにおける各都市の圧力ベクトル$\boldsymbol{p}$が求まります。

線路太さ（輸送管太さ）の変化

都市ごとの圧力$p_{i}$が求まりましたので、これをもとにそれぞれのエッジ（都市間をつなぐ管）における（人）流量$Q_{ij}$を算出することができます。

Q_{ij}=\frac{D_{ij}}{L_{ij}}(p_{i}-p_{j})

ここで、算出した（人）流量$Q_{ij}$をもとに線路太さ（輸送管太さ）$D_{ij}$を変化させていきます。これは流量強化関数$f(|Q|)$によって行われ、こちらも参考記事の手法をそのままお借りしています。

\displaylines{
f(|Q_{ij}|)=\frac{|Q_{ij}|^\gamma}{1+|Q_{ij}|^\gamma}\\
D_{ij}^{(n+1)}=D_{ij}^{(n)}+ \lbrack f(|Q_{ij}|)-D_{ij}^{(n)} \rbrack dt\\
}

ここで$\gamma$は1.8などの定数、$dt$は0.1などの微小量、$D_{ij}^{(n)}$は$n$ステップ目における$D_{ij}$を表します。

$f(|Q|)$は$Q$の絶対値に応じて1に漸近する関数です。

したがって、次ステップにおける線路太さ（輸送管太さ）$D_{ij}$は流量が大きければ太くなり、流量が小さければより細くなることがわかります。しかしながら、$f(|Q|)$は常に1以下の値をとるため、本プログラムのように初期値を$D_{ij}=1$と設定した場合、次のステップの$f(|Q|)-D_{ij}$は必ず負の値となってしまい、同様の計算を繰り返すといずれはすべての$D_{ij}$が0に落ち着いてしまいます。これは初期値をこのような値にしたために起こる問題なのですが、この問題を解決するため、本プログラムでは以下のような「正規化」を行いました。

D_{ij}^{(n+1) '}=\frac{\sum_{i,j}D_{ij}^{(n)}}{\sum_{i,j}D_{ij}^{(n+1)}} D_{ij}^{(n+1)}

これは全ステップを通じて$D_{ij}^{(n)}$の総量が初期値から変わらないように値を調整するものです（体積一定の粘菌の物理モデルに即したものかはわかりませんが、本プログラムでは十分機能したのでこの式を採用しています）。この操作により、流量の少ない輸送管の太さ$D_{ij}$が減る分、流量の多い一定数のルートはその分$D_{ij}$の値が増え、太くなっていくことになります。
この値でもって$D_{ij}$を更新することで次のステップへと入り、再び人口に応じた確率で出発都市と到着都市を決定し、同様の計算を行っていきます。

実行例

以下の８都市が与えられたとします。

X	Y	人口
18	80	4
91	87	4
85	93	1
64	79	1
50	64	8
87	15	1
15	36	1
64	25	1

これに対し、本プログラムを実行すると以下のような結果が得られました。

この問題では、最も人口の多い都市から放射状に線路を伸ばす鉄道網が提案されました。

ソースコード

最後に、作成したエクセルシートおよびVBAコードを以下に示します。

Sub InitializeTheGraph()

'[参考文献]
'https://qiita.com/STInverSpinel/items/8ced06ea7881613a3e2c

Range("H4:L2400").Select
Selection.ClearContents

End Sub

Sub GraphMaking()

'都市の座標および人口データからグラフを作成

Dim X(100), Y(100) As Double
Dim popl(100) As Double
Dim L(100, 100) As Double
Dim D(100, 100) As Double

'都市数の取得
NoN_X = Cells(4, 2).End(xlDown).Row - 3
NoN_Y = Cells(4, 3).End(xlDown).Row - 3
If (NoN_X = NoN_Y) Then
    NoN = NoN_X
End If

'都市座標の取得
For i = 1 To NoN
    X(i) = Cells(3 + i, 2)
    Y(i) = Cells(3 + i, 3)
Next

'人口の取得
For i = 1 To NoN
    popl(i) = Cells(3 + i, 4)
Next

'グラフ描画
k = 1
For i = 1 To NoN
    For j = i + 1 To NoN
        '番号
        Cells(3 + k, 8) = i
        Cells(3 + k + 1, 8) = j
        
        '座標
        Cells(3 + k, 9) = X(i)
        Cells(3 + k, 10) = Y(i)
        Cells(3 + k + 1, 9) = X(j)
        Cells(3 + k + 1, 10) = Y(j)
        
        '距離L
        L(i, j) = Sqr((X(i) - X(j)) ^ 2 + (Y(i) - Y(j)) ^ 2)
        Cells(3 + k, 11) = L(i, j)
        
        '太さD
        D(i, j) = 1
        Cells(3 + k, 12) = D(i, j)
        
        k = k + 3      
    Next
Next

Call GraphUpdate

'最初は都市間ノードを消しておく
ActiveSheet.ChartObjects("グラフ 1").Activate
ActiveChart.FullSeriesCollection(1).Select
Selection.Format.Line.Visible = msoFalse

End Sub

Sub Step()

'1ステップ進める

Dim popl(100) As Double
Dim F(100) As Double
Dim p(100) As Double
Dim L(100, 100) As Double
Dim D(100, 100) As Double
Dim a(100, 100) As Double
Dim Q(100, 100) As Double

'都市数の取得
NoN_X = Cells(4, 2).End(xlDown).Row - 3
NoN_Y = Cells(4, 3).End(xlDown).Row - 3
If (NoN_X = NoN_Y) Then
    NoN = NoN_X
End If

'現在のLとDの取得（上三角部分）
For i = 1 To NoN
    For j = i + 1 To NoN

        L(i, j) = Cells(3 + k, 11)
        D(i, j) = Cells(3 + k, 12)    
        k = k + 3

    Next
Next

'下三角部分への反映
For i = 1 To NoN
    For j = i + 1 To NoN
        L(j, i) = L(i, j)
        D(j, i) = D(i, j)
    Next
Next

'現在のDの総和
TD = 0
k = 1
For i = 1 To NoN
    For j = i + 1 To NoN
        TD = TD + D(i, j)
    Next
Next

'出発（流出；湧出）都市と到着（流入；吸込）都市

'連立一次方程式"Ap=F"の"F"
'初期化
For i = 1 To NoN
    F(i) = 0
Next
'人口の取得
'popl(k): k番目の都市における人口
'TotalPop:地図全体の総人口
    TotalPopu = 0
For i = 1 To NoN
    popl(i) = Cells(3 + i, 4)
    TotalPopu = TotalPopu + popl(i)
Next

'人口に応じた出発都市と到着都市の選出
Do While (RouFr >= RouTo)
    'TotalPopu以下の自然数をランダムに取得
    an = Int(TotalPopu * Rnd + 1)
    bn = Int(TotalPopu * Rnd + 1)

    '出発都市番号
    RouFr = 1
    Do While (an > 0)
        an = an - popl(RouFr)
        RouFr = RouFr + 1
    Loop

    '到着都市番号
    RouTo = 1
    Do While (bn > 0)
        bn = bn - popl(RouTo)
        RouTo = RouTo + 1
    Loop

    '計算上の操作
    RouFr = RouFr - 1
    RouTo = RouTo - 1
Loop

'出発都市から流量I_0の流出
'到着都市に流量I_0の流入
'それ以外は流入流出の総和が0であるようにする
F(RouFr) = 1.2
F(RouTo) = -1.2

'行列Aの作成
'初期化
For i = 1 To NoN
    For j = 1 To NoN
        a(i, j) = 0
    Next
Next
'行列の作成
For i = 1 To NoN
    For j = 1 To NoN
        If (i = j) Then
            For k = 1 To NoN
                If (i <> k) Then
                    a(i, j) = a(i, j) + D(i, k) / L(i, k)
                End If
            Next
        Else
            a(i, j) = -D(i, j) / L(i, j)
        End If
    Next
Next

'ガウス=ザイデル法で連立一次方程式Ap=Fを解く
'参考文献
'http://nkl.cc.u-tokyo.ac.jp/14n/SolverIterative.pdf
BNRM = 0#
For i = 1 To NoN
    BNRM = BNRM + F(i) ^ 2
Next
BNRM = Sqr(BNRM)

For iter = 1 To NoN * 2400
    For i = 1 To NoN
        RESID = F(i)
        For j = 1 To NoN
           If (j <> i) Then
              RESID = RESID - a(i, j) * p(j)
           End If
        Next
        p(i) = RESID / a(i, i)
    Next
    
    '収束判定
    VALU = 0#
    For i = 1 To NoN
        RESID = F(i)
        For j = 1 To NoN
            RESID = RESID - a(i, j) * p(j)
        Next
        VALU = VALU + (RESID) ^ 2
    Next
    
    VALU = Sqr(VALU) / BNRM
    
    If (VALU < 0.0000001) Then
        Exit For
    End If
Next

'次世代のDを計算
For i = 1 To NoN
    For j = i + 1 To NoN
        Q(i, j) = (D(i, j) / L(i, j)) * (p(i) - p(j))
        
        'f(|Q|)の計算
        fQ = (Abs(Q(i, j)) ^ 1.8) / (1 + Abs(Q(i, j)) ^ 1.8)
        
        'Dの更新
        D(i, j) = D(i, j) + (fQ - D(i, j)) * 0.1
    Next
Next

'次世代Dの総和
TDNG = 0
For i = 1 To NoN
    For j = i + 1 To NoN  
        TDNG = TDNG + D(i, j)
    Next
Next

'「正規化」したD値を更新
k = 1
For i = 1 To NoN
    For j = i + 1 To NoN
        D(i, j) = D(i, j) * (TD / TDNG)
        Cells(3 + k, 12) = D(i, j)
        k = k + 3
    Next
Next

End Sub

Sub SeveralStep()

'上記ステップを反復して呼び出す
For iter = 1 To 4200
    Call Step
Next

'グラフへの反映
Call GraphUpdate

End Sub

Sub GraphUpdate()

'グラフ更新作業

Dim popl(100) As Double

'都市数の取得
NoN_X = Cells(4, 2).End(xlDown).Row - 3
NoN_Y = Cells(4, 3).End(xlDown).Row - 3
If (NoN_X = NoN_Y) Then
    NoN = NoN_X
End If

'人口の取得
For i = 1 To NoN
    popl(i) = Cells(3 + i, 4)
Next

'グラフのアップデート
ActiveSheet.ChartObjects("グラフ 1").Activate
ActiveChart.FullSeriesCollection(1).Select
Selection.Format.Line.Visible = msoTrue

For m = 1 To ((NoN) ^ 2) * 2

    numb = Cells(3 + m, 8)
    D_ij = Cells(3 + m - 1, 12)
    
    'データ点サイズを人口に応じて変化させる
    If (numb <> "") Then
        ActiveChart.FullSeriesCollection(1).Points(m).Select
        With Selection
            .MarkerStyle = 1
            .MarkerSize = popl(numb) + 4
        End With
    End If
    
    'エッジ太さおよび透明度をD_ij値によって変化させる
    If (numb <> "" And Cells(3 + m + 1, 8) = "") Then
        ActiveChart.FullSeriesCollection(1).Points(m).Select
        With Selection.Format.Line
            .Weight = Sqr(D_ij)
            .Transparency = (1 - ((D_ij ^ 1.8) / (1 + D_ij ^ 1.8))) * 0.1
        End With
    End If
    
Next

End Sub

参考文献

• Atsushi Tero, Ryo Kobayashi, and Toshiyuki Nakagaki, "Physarum solver: A biologically inspired method of road-network navigation", Physica A: Statistical Mechanics and its Applications Volume 363, Issue 1, 15 April 2006, Pages 115-119
• 手老篤史「粘菌の輸送ネットワークから都市構造の設計理論を構築」―都市間を結ぶ最適な道路・鉄道網の法則確立に期待― 科学技術振興機構報, 第708号, 2010年1月22日
• Atsushi Tero, Seiji Takagi, Tetsu Saigusa, Kentaro Ito, Dan P. Bebber, Mark D. Fricker, Kenji Yumiki, Ryo Kobayashi, Toshiyuki Nakagaki, "Rules for Biologically Inspired Adaptive Network Design", Science 327, 22 January 2010, Pages 439-442

脚注

1. 例えば手老氏らはネットワーク性能の評価として「総距離（構築・維持コスト）vs輸送効率（ある２都市間をネットワーク上を通って移動する効率の平均）」、「事故や災害による断線への耐性」などの観点から評価を行っています。 ↩

2. そしてこの操作は、例えば地方都市から大都市へと買い物に出かけたり、大都市から地方都市へと働きに出かける人の流れを１ステップにつき１人ずつ観察していくプロセスを模しているとも捉えることができ、地方から地方への移動よりも人口の多い都市における人の出入りにスポットが当たりやすいことにも対応しています。 ↩

3. 私は反復法の一つであるGauss-Seidel法を用いました。行列の対称性を利用したうまいやり方があるのかもしれませんが、今のところは見つかっていません。 ↩