群集シミュレーションを用いた発生交通量推定

本問題はEC-Comp 2022 単目的部門の問題です．本問題では，実際の花火大会（関門海峡花火大会）における帰宅時のシミュレーション結果が観測データに近づくよう，帰宅者の出発時刻および経路選択を推定（最適化）します．

問題の背景

近年，将棋倒しや群集雪崩などに代表される群集災害が大勢の人が集まるイベント/施設で問題となっています．群集災害を回避するためには，標識や警備員によって適切に群集の流れを制御する必要があります．しかし，費用や安全性の観点から様々な制御方策の是非を実際の環境で評価することは困難です．

多くの人やモノの動きを群集シミュレーションによりコンピュータ上（仮想世界）で模擬することで，現実世界に比べて迅速かつ定量的に制御方策を評価することが出来ます．とくに近年，各エージェントが個々の歩行者を表現するマルチエージェントシミュレーションを用いることで，制御方策の最適化を行う研究事例が多数報告されています3,4．

適切な群集制御最適化を行うためには，シミュレータを現実に近づけることが欠かせません．現実に則したシミュレーション結果を得るためには，i) 移動経路に関する地図情報，ii) 歩行者の発生時刻や起点終点に関する情報（OD交通量），iii) 移動速度や経路選択などの歩行者の挙動を表す数理モデル，iv) 群集制御や災害による通行止めなどの外的要因といった情報を正しく与える必要があります．本コンペティションでは，なかでも計測が困難なOD交通量と歩行者モデルのパラメータに注目し，群集シミュレーション結果と観測データの誤差最小化によりこれらのパラメータを推定します．

群集シミュレーションの概要

本問題では，帰宅者の移動のシミュレーションにCrowdWalkを利用します．CrowdWalkはひとりひとりの歩行者に対応するエージェントを生成してシミュレーションを行うマルチエージェント型の群集シミュレータです．各歩行者にエージェントを対応させることで，各歩行者のミクロな動きまで想定したシミュレーションを行う事ができます．

Fig. 1: CrowdWalkで使用する一次元のネットワークマップ

多数のエージェントの移動を高速かつ低メモリで計算するために，CrowdWalk ではFig. 1のように一次元のネットワークマップを利用してエージェントの位置計算を行います．また，全てのエージェントの位置は Social Force Model に従って更新されます．一次元ネットワークマップに適した Social Force Model の詳しい計算は，山下らの論文を参照ください5．

問題定義

概要

本問題では，実際の花火大会（関門海峡花火大会）における帰宅時のシミュレーション結果が観測データに近づくよう，帰宅者の出発時刻および経路選択を推定（最適化）します．環境変数MODELをSOP-1かSOP-2に設定することで，それぞれ以下の問題を扱うことができます．

• SOP-1: 帰宅者の出発時刻のみを最適化する問題
• SOP-2: 出発時刻および経路選択を最適化する問題

定式化

Fig. 2: シミュレーション結果と観測値の誤差最小化

本問題は，CrowdWalkを用いた群集シミュレーションにおける歩行者情報推定を比較するためのベンチマーク （MAS-Bench）をもとに作成されています．MAS-Benchでは，以下の定義に従って観測データとの誤差が少なくなるよう歩行者情報が最適化されます．

$${\theta }^{\ast }={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{\theta \in \mathrm{\Theta }}{ϵ}^{\mathrm{o}\mathrm{b}\mathrm{s}}(O,S(\theta ))\backslash theta^*\; =\; \{\backslash rm\; argmin\}\_\{\backslash theta\; \backslash in\; \backslash Theta\}\; \backslash epsilon^\{\backslash rm\; obs\}(O,\; S(\backslash theta))$$

ここで， $\theta \backslash theta$ は出発時刻や経路選択に関する歩行者パラメータ（決定変数）， $OO$ はカメラやGPSから得た観測データ， $S(\theta )S(\backslash theta)$ はパラメータ $\theta \backslash theta$ を用いたシミュレーションの結果， ${ϵ}^{\mathrm{o}\mathrm{b}\mathrm{s}}(O,S(\theta ))\backslash epsilon^\{\backslash rm\; obs\}(O,\; S(\backslash theta))$ は観測データとシミュレーション結果の誤差評価値（目的関数）です．

最適化の流れ

Fig. 3: MAS-Benchのシステム構成

MAS-Benchのシステム構成はFig. 3のようになっています．ユーザ側が操作する内容（最適化計算）は黒波線で囲まれており，MAS-Bench側が操作する内容（群集シミュレーション計算）は赤波線で囲まれています．
MAS-Bench側は決定変数 $\theta \backslash theta$ をユーザ側から受け取り，目的関数 ${ϵ}^{\mathrm{o}\mathrm{b}\mathrm{s}}(O,S(\theta ))\backslash epsilon^\{\backslash rm\; obs\}(O,\; S(\backslash theta))$ をユーザ側に返します．
つまり，群集シミュレーションの実行とシミュレーション結果に対する評価計算はMAS-Bench側で行われており，ユーザ側は歩行者情報の推定をブラックボックス最適化として行うことが可能です．

シミュレーションの構成要素

地図情報

Fig. 6: 実際の花火大会における交通制御の様子

本問題では，関門海峡花火大会終了時におけるメイン会場から最寄り駅までの群集移動を対象としました．関門海峡花火大会の会場周辺マップ（Google航空写真）をFig. 4に示します．会場周辺マップにおいて，Oは花火イベント会場, Dは門司港駅を表します．
関門海峡花火大会では，交通規制により観客の帰宅経路は大きく3種類に分けられます．また，それぞれの経路上には警備員が配置されており，警備員はルートを変更する分岐制御と歩行者の動きを止める停止制御を行います．シミュレーションにおいてエージェントが移動するルートの情報を簡略化したネットワークマップがFig. 5になります．各ルートの詳細は以下の通りです．

• Route1
  最も短い経路（Route1: 320[m]）．配置されている警備員の数は分岐制御1人，停止制御3人です．
• Route2
  次に短い経路（Route2: 600[m]）．配置されている警備員の数は分岐制御2人，停止制御2人です．
• Route3
  最も長い経路（Route3: 770[m]）．配置されている警備員の数は分岐制御2人，停止制御2人です．

多くの歩行者の帰宅経路は，分岐制御における警備員の指示によって決定します．必ずRoute1を通る（警備員がRoute1を指示するまで立ち止まる）といった指示に従わない歩行者も一部存在します．詳しくは歩行者モデルをご覧ください．また，シミュレータ上でも実際の制御と同様にエージェントの移動制御を行っています．

計測機器

Fig. 7: 2種類の計測機器

本問題では，関門海峡花火大会で実際に得られている２種類の観測データを用いてシミュレーション結果と実際の帰宅者情報の誤差を求めます．

• カメラ
  駅(歩行者の目的地)に到着した人数を計測し，時間あたりの到着人数を記録します．
• GPS
  GPSトラッカーを持ったアルバイトに一定の時間間隔（10分ごと）でメイン会場から駅まで移動してもらい，それぞれ時間あたりの積算移動距離を記録します．群集の流れを計測するために，アルバイトは周りの流れに合わせて移動しています．

歩行者モデル

定点カメラを用いた群集計測から花火イベントでは，誘導に従う参加者と従わない参加者がいることが分かっています．本問題では参加者を，誘導に従うエージェント（Guided），駅への最短経路を選択するエージェント（Busy），お祭り会場となっている長い経路を選択するエージェント（Slow）の３種類に分類します．また，これら全てのエージェントは停止信号に従うものとします．

• Guided
  分岐信号にしたがって歩行するエージェントです．このエージェントは全ての経路を歩く可能性があります．
• Busy
  Guidedより少し速く歩行するエージェントです．このエージェントは必ず最短経路Route1を歩きます．
• Slow
  Guidedより少し遅く歩行するエージェントです．このエージェントはRoute1への誘導を無視します．その後は誘導に従い，Route2,3のどちらかを歩きます．

決定変数

Fig. 8: 出発人数ヒストグラムと混合ガウス分布による近似

本問題では，全ての歩行者が花火会場を出発し駅まで移動すると仮定します．出発地・目的地が固定となるため，決定変数となる歩行者パラメータは時刻ごとの各エージェントタイプの出発人数として表現できます．人数は1分ごとに集計されるものとし，19時から24時までの各時刻の出発人数を整数値で指定します．すなわち，決定変数の数は300分✖エージェントタイプ数となります．歩行者の合計（総数）は駅に設置したカメラから分かるため，決定変数の合計値と歩行者の総数は一致する必要があります．
30分ごとに人数を集計した場合の出発人数ヒストグラムとそれに対応した決定変数を図に示します．各時刻におけるエージェントタイプごとの出発人数を推定することで，「いつ」に加えて「どのような人が」に関する歩行者情報も推定しています．

---

（Tips）混合ガウス分布で実数最適化問題として解く方法
上記の決定変数を直接最適化する場合，300分×3種類=900個の変数を最適化する必要があります．出題者の研究グループでは出発人数ヒストグラムを混合ガウスモデル (GMM) を用いて近似することで実数最適化問題として解きましたので，実数値を用いる最適化手法を使う参加者の参考として説明いたします．

GMM法による出発人数推定では，エージェントタイプごとの出発人数をひとつのガウス分布に対応させて推定します．その後，ガウス分布を度数分布に変換し，合計が総人数と一致するよう補正することで時間あたりの出発人数を計算します．

エージェントタイプ $kk$ に対応するガウス分布を $f_k(t;{\mu }_k,{\sigma }_k,{\pi }_k)f\_\{k\}(t;\; \backslash mu\_k,\; \backslash sigma\_k,\; \backslash pi\_k)$ とします． ${\mu }_k\backslash mu\_k$ ， ${\sigma }_k\backslash sigma\_k$ はガウス分布の平均および分散， ${\pi }_k\backslash pi\_k$ は総人数に対するエージェントタイプ $kk$ の比率です．全てのエージェントの出発人数分布 $f(t;\mu ,\sigma ,\pi )f(t;\; \backslash mu,\; \backslash sigma,\; \backslash pi)$ は次のように計算できます．

$f_k(t;{\mu }_k,{\sigma }_k,{\pi }_k)={\displaystyle \frac{C_k{\pi }_k}{\sqrt{2\pi {\sigma }_k^2}}\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{(t-{\mu }_k{)}^2}{2{\sigma }_k^2})}f\_\{k\}(t;\; \backslash mu\_k,\; \backslash sigma\_k,\; \backslash pi\_k)\; =\; \{\backslash displaystyle\; \backslash frac\{C\_\{k\}\backslash pi\_\{k\}\}\{\backslash sqrt\{2\backslash pi\backslash sigma\_\{k\}^2\}\}\{\backslash rm\; exp\}\backslash biggl(-\backslash frac\{(t\; -\; \backslash mu\_\{k\})^2\}\{2\backslash sigma\_\{k\}^2\}\; \backslash biggl)\}$
$f(t;\mu ,\sigma ,\pi )={\sum }_{k=1}^K{f}_k(t;{\mu }_k,{\sigma }_k,{\pi }_k)f(t;\; \backslash mu,\; \backslash sigma,\; \backslash pi)\; =\; \backslash sum\_\{k=1\}^\{K\}\; f\_\{k\}(t;\; \backslash mu\_k,\; \backslash sigma\_k,\; \backslash pi\_k)$

ここで， $C_k=1\mathrm{/}{\int }_0^T{f}_k(t;{\mu }_k,{\sigma }_k)dtC\_\{k\}=1/\backslash int\_\{0\}^\{T\}f\_\{k\}(t;\; \backslash mu\_k,\; \backslash sigma\_k)dt$ は正規化定数， $t\in Tt\; \backslash in\; T$ はシミュレーション時間を示しています．各ガウス分布の比率 ${\pi }_k\backslash pi\_\{k\}$ は，必ず ${\sum }_{k=1}^K{\pi }_k=1\backslash sum\_\{k=1\}^\{K\}\; \backslash pi\_\{k\}\; =\; 1$ を満たします．

次に，時間あたりの出発人数を計算するためにガウス分布を度数分布に変換します．

$N_k^{\mathrm{o}\mathrm{g}\mathrm{n}}(t;{\mu }_k,{\sigma }_k,{\pi }_k)=\lfloor \alpha {\int }_t^{t+\mathrm{\Delta }t}{f}_k(t;{\mu }_k,{\sigma }_k,{\pi }_k)dt\rfloor N^\{\backslash rm\; ogn\}\_\{k\}(t;\backslash mu\_k,\; \backslash sigma\_k,\; \backslash pi\_k)\; =\; \backslash biggl\backslash lfloor\; \backslash alpha\; \backslash int\_\{t\}^\{t+\backslash Delta\; t\}\; f\_\{k\}(t;\; \backslash mu\_k,\; \backslash sigma\_k,\; \backslash pi\_k)\; dt\backslash biggl\backslash rfloor$
$N^{\mathrm{o}\mathrm{g}\mathrm{n}}(t;\mu ,\sigma ,\pi )={\sum }_{k=1}^K{N}_k^{\mathrm{o}\mathrm{g}\mathrm{n}}(t;{\mu }_k,{\sigma }_k,{\pi }_k)N^\{\backslash rm\; ogn\}(t;\; \backslash mu,\; \backslash sigma,\; \backslash pi)\; =\; \backslash sum\_\{k=1\}^\{K\}\; N^\{\backslash rm\; ogn\}\_\{k\}(t;\backslash mu\_k,\; \backslash sigma\_k,\; \backslash pi\_k)$

ここで $\mathrm{\Delta }t\backslash Delta\; t$ はシミュレーションのタイムステップサイズ， $N_k^{\mathrm{o}\mathrm{g}\mathrm{n}}(t)N^\{\backslash rm\; ogn\}\_\{k\}(t)$ はエージェントタイプ $kk$ の時刻 $tt$ における出発人数です．なお，歩行者の総人数は分かっているため， ${\sum }_{m=1}^M{N}^{\mathrm{o}\mathrm{g}\mathrm{n}}(t_m)=N\backslash sum\_\{m\; =\; 1\}^\{M\}\; N^\{\backslash rm\; ogn\}(t\_m)\; =\; N$ となるよう調整パラメータ $\alpha \backslash alpha$ によって人数を補正しています．

以上の変換を行う事で，出発人数ヒストグラムにおける整数値の最適化を実数パラメータ ${\mu }_k\backslash mu\_k$ ， ${\sigma }_k\backslash sigma\_k$ ， ${\pi }_k\backslash pi\_k$ の最適化に近似することができます．

目的関数

本問題では，実観測データとして駅に設置したカメラで計測した到着人数（ヒストグラム）とアルバイトの持つGPSによる移動軌跡を用います．シミュレーション結果と実観測データの誤差を評価するために，到着人数とGPS移動軌跡はそれぞれ二乗平均平方根誤差(RMSE)を用いて比較され，最終的な評価は二つの誤差値の線形結合によって表されます．計算式は次の通りです．

${ϵ}^{\mathrm{o}\mathrm{b}\mathrm{s}}=\beta {ϵ}^{\mathrm{d}\mathrm{e}\mathrm{s}\mathrm{t}}+(1-\beta ){\overline{ϵ}}^{\mathrm{p}\mathrm{a}\mathrm{c}}\backslash epsilon^\{\backslash rm\; obs\}\; =\; \backslash beta\; \backslash epsilon^\{\backslash rm\; dest\}\; +\; (1\; -\; \backslash beta)\; \backslash bar\{\backslash epsilon\}^\{\backslash rm\; pac\}$

ここで， $\beta \backslash beta$ は誤差評価にカメラとGPSが与える重要度を決める[0, 1]の重みであり，本問題では0.5とします．

到着人数に関する誤差評価は，単位時間 $\mathrm{\Delta }t\backslash Delta\; t$ あたりの到着人数に対するRMSEによって計算されます．本問題ではシミュレーション開始から終了までを$M$回に分けて取り扱います．花火大会では19:00から24:00までシミュレーションするため， $\mathrm{\Delta }t=300min\mathrm{.}\mathrm{/}M\backslash Delta\; t=300min./M$ となります．到着人数に関する誤差評価の計算式は次の通りです．

${ϵ}^{\mathrm{d}\mathrm{e}\mathrm{s}\mathrm{t}}=\sqrt{\frac{1}{M}{\sum }_{m=1}^M(N^{\mathrm{d}\mathrm{e}\mathrm{s}\mathrm{t}}(t_m)-{\widehat{N}}^{\mathrm{d}\mathrm{e}\mathrm{s}\mathrm{t}}(t_m){)}^2}where{t}_m=m\mathrm{\Delta }t\backslash epsilon^\{\backslash rm\; dest\}\; =\; \backslash sqrt\{\backslash frac\{1\}\{M\}\; \backslash sum\_\{m\; =\; 1\}^\{M\}\; (N^\{\backslash rm\; dest\}(t\_m)\; -\; \backslash hat\{N\}^\{\backslash rm\; dest\}(t\_m))^2\}\; \backslash qquad\; where\backslash qquad\; t\_\{m\}\; =\; m\backslash Delta\; t$

ここで， $N^{\mathrm{d}\mathrm{e}\mathrm{s}\mathrm{t}}(t_m)N^\{\backslash rm\; dest\}(t\_m)$ は時刻 $t_{m}t\_\{m\}$ において実際に観測された到着人数，${\widehat{N}}^{\mathrm{d}\mathrm{e}\mathrm{s}\mathrm{t}}(t_m)\backslash hat\; N^\{\backslash rm\; dest\}(t\_m)$ はシミュレーションで得られた到着人数です．

現実世界の測定における不確実性を考慮するため， ${\widehat{N}}^{\mathrm{d}\mathrm{e}\mathrm{s}\mathrm{t}}(t_m)\backslash hat\; N^\{\backslash rm\; dest\}(t\_m)$ を計算する際には駅に到着したエージェントの数に対してガウシアンフィルタによる平滑化を行っています．ガウシアンフィルタによる平滑化は以下の式により行います．

$G(t)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp (-\frac{t^2}{2{\sigma }^2})G(t)\; =\; \backslash frac\{1\}\{\backslash sqrt\{2\; \backslash pi\; \backslash sigma^2\}\}\; \backslash exp\{\backslash biggl(-\; \backslash frac\{t^2\}\{2\; \backslash sigma^2\}\backslash biggl)\}$

${\widehat{N}}^{\mathrm{d}\mathrm{e}\mathrm{s}\mathrm{t}}(t_m)={\sum }_{q=1}^M{\widehat{N}}^{\mathrm{a}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{t}}(t_q)\cdot G(t_q-t_m)\backslash hat\{N\}^\{\backslash rm\; dest\}(t\_m)\; =\; \backslash sum\_\{q\; =\; 1\}^\{M\}\backslash hat\{N\}^\{\backslash rm\; agent\}(t\_q)\; \backslash cdot\; G(t\_q\; -\; t\_m)$

${\widehat{N}}^{\mathrm{a}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{t}}(t_q)\backslash hat\{N\}^\{\backslash rm\; agent\}(t\_q)$ は時刻 $t_{q}t\_q$ で駅に到着したエージェントの数です．

移動軌跡に関する誤差評価は，合計移動距離を単位時間ごとに比較して求められます．GPSによる移動軌跡に対応するデータを得るため，シミュレーションでは周りのエージェントの平均速度で移動するエージェント（Tracked Agent）を実際のアルバイトと同じ時間間隔で投入し，それぞれ時間あたりの積算移動距離を記録しました．
RMSEの計算はTracked Agentごとに行うため，全Tracked AgentのRMSEの平均値を誤差評価とします．計算式は次の通りです．

${\overline{ϵ}}^{\mathrm{p}\mathrm{a}\mathrm{c}}=\frac{1}{A}{\sum }_{a=1}^A\sqrt{\frac{1}{M_a}{\sum }_{m=1}^{M_a}(P(a,t_m)-\widehat{P}(a,t_m){)}^2}\backslash bar\{\backslash epsilon\}^\{\backslash rm\; pac\}\; =\; \backslash frac\{1\}\{A\}\backslash sum\_\{a=1\}^\{A\}\; \backslash sqrt\{\backslash frac\{1\}\{M\_a\}\; \backslash sum\_\{m\; =\; 1\}^\{M\_a\}\; (P(a,\; t\_m)\; -\; \backslash hat\{P\}(a,\; t\_m))^2\}$

ここで，AはGPS所有者の数，総移動時間 $M_{a}M\_a$ はGPS所有者 $aa$ の出発時刻 $t_a^{\mathrm{o}\mathrm{g}\mathrm{n}}t\_\{a\}^\{\backslash rm\; ogn\}$ と到着時刻 $t_a^{\mathrm{d}\mathrm{e}\mathrm{s}\mathrm{t}}t\_\{a\}^\{\backslash rm\; dest\}$ により $M_a=(t_a^{\mathrm{d}\mathrm{e}\mathrm{s}\mathrm{t}}-t_a^{\mathrm{o}\mathrm{g}\mathrm{n}})\mathrm{/}\mathrm{\Delta }tM\_a\; =\; (t\_\{a\}^\{\backslash rm\; dest\}\; -\; t\_\{a\}^\{\backslash rm\; ogn\})\; /\; \backslash Delta\; t$ である. $P(a,t_m)P(a,\; t\_m)$ は時刻 $t_{m}t\_m$ におけるGPS所有者 $aa$ の合計移動距離， $\widehat{P}(a,t_m)\backslash hat\{P\}(a,\; t\_m)$ は時刻 $t_{m}t\_m$ における $aa$ に対応したTracked Agentの合計移動距離です．

---

参考文献

1. C. Antoniou et al. (2016), Towards a generic benchmarking platform for origin–destination flows estimation/updating algorithms: Design, demonstration and validation, Transportation Research Part C: Emerging Technologies, vol. 66, pp. 79-98.
   
2. D. Wolinski et al. (2014), Parameter estimation and comparative evaluation of crowd simulations, Computer Graphics Forum, vol. 33, no. 2, pp. 303-312.
   
3. R. K. Dubey et al. (2020), AUTOSIGN: A multi-criteria optimization approach to computer aided design of signage layouts in complex buildings, Computers & Graphics, vol. 88, pp. 13-23.
   
4. A. Abdelghany et al. (2014), Modeling framework for optimal evacuation of large-scale crowded pedestrian facilities, European Journal of Operational Research, vol. 237, no. 3, pp. 1105-1118.
   
5. 山下倫央 et al. (2012), 一次元歩行者モデルを用いた高速避難シミュレータの開発とその応用, 情報処理学会論文誌, vol. 53, no. 7, pp. 1732-1744.