ドローンによる物品配送サービスの待ち行列モデル
A queueing model for goods delivery service by drones
栗田 治 Osamu Kurita
In many depopulated rural districts in Japan, it has been hard to run regional retail stores and many of them have closed down. As the result, there has been appeared many residents who feel difficulties to purchase commodities for their own. They are called “kaimono-jakusha” in Japan. To solve this social problem home delivery services by drones are proposed and several social experiments have been carried out. In this time we developed M/G/s queueing models for the drone system depending on the methods of urban operations research. Finally, we derive the number of drones to achieve the expected delivery efficiency.
Keywords : drone, home delivery service, M/G/s queueing model, distance distribution ドローン，宅配サービス，M/G/s 待ち行列モデル，距離分布 1．緒言 我が国は長期的な人口減少過程に入っている．加えて地 方から大都市圏への若年層と生産年齢人口の移住も進んでい る．そのため，多くの地方を人口減少と高齢化の双方が襲 い，利潤を確保し辛くなった小売店が次々に撤退しつつあ る．その結果，自動車を利用できない世帯は，日常の買い物 さえままならなくなり，いわゆる“ 買い物弱者 ”が増加して いる．買い物弱者が集積した地域・地区は“ food desert ” と呼ばれ，重要な都市問題と見做されている8,12,17)．さら に，地方におけるバス路線の多くが採算悪化のために廃止 されるという現象が，この問題に拍車をかけている． この都市問題を解決するための一手段として注目されて いる，ドローン (遠隔操縦あるいは自律式の無人航空機) に よる物品の宅配システムに焦点を当てる．無人で直線的に 移動できるという長所を活かし，デポ (配送拠点) に配備し た複数のドローンを用いる．ただしドローン宅配の導入に 際しては，解決すべき重要な法的・技術的課題がある： (a) 住宅地を飛行する際の安全性確保と飛行の許可 (b) ドローンの自律的な無視界飛行の技術的達成 (c) 複数ドローンの同時飛行制御システムの構築 技術的解決に向けた研究開発と実証実験は世界的に進んでい るため2,6,14,15,16)，ドローン宅配システムが近い将来に実 現する可能性は大きいものと想像する．また，ドローンの 安全性を確保しつつ導入するための法制度については，現 在国を挙げてこれをサポートする体制がとられつつあるこ とから，これも近い将来，まずは人口密集地以外の場所で 導入できる可能性を十分に秘めているものと思われる． 本論文は，上記の如き課題が達成されドローン宅配が実 現されるときの，サービス水準を見積もるためのモデルを提 案する．今回の成果は (1) ドローン基地の立地選定とドロー ————————————————————————————————————————————————– 正会員 慶應義塾大学理工学部管理工学科 (Keio University) ン配備数の決定を支援し，(2) 飛行速度の増大や (3) ドロー ン発進前の準備時間短縮等の企業努力が，サービス水準向 上に与える効果を見積もるために役立つ． 今回の数値例では円盤領域で一様な需要を与え，注文が ポアソン到着するものと想定する．そして，円盤中心から ドローン基地を偏心させて計算を行う．待ち時間の特性値 を算出する過程で，直線距離分布の理論7,9,10,11)が重要な 役割を果たす．そして，平均待ち時間の算出には木村の高 精度近似公式3,4,5)を採用した．本モデルに基づいて，一般 の都市領域の形状や現実の (離散的な) 需要分布に対応した 計算を行うことも可能である (5 章)．さらに，本研究のモ デルはオートバイや自転車等を用いて少量の物品を個別に宅 配するという従来型の配送問題にもスムーズに適用できる． なお，ドローンによる配送の先行研究は数多いが，待ち 行列現象を中心に据えたものはあまりなく，また距離分布 に基づく都市解析的な研究は皆無である．AED(自動体外式 除細動器) の基地を整数計画法で決定し，配備された複数機 のドローンによる配送を M/M/s 待ち行列系で記述した例 はあるが19)，本稿とは異なり，サービス時間が指数分布に 従うとしている点に現実からの乖離がある．また M/G/s 待ち行列系を用いた例21)は，シミュレーション分析に終 始している．しかも，これら 2 つの研究は基地の立地・ド ローン配備数・ドローンのスペック・注文到着率といった 要件が配送効率に与える影響を分析していない．その面で も本稿に一日の長がある．多くのドローン配送研究は，ト ラック輸送のみの場合に比べた優位性を示す内容や20,23)， 複数のドローン基地による単純なカバリング問題を追求す る内容であり22)，本稿とは方向性が大きく異なっている． ドローン宅配の可能性を把握するために，国内外のドロー ンによる宅配実験や現場での利用状況を概観しておく．
瀬戸内海の大三島 (愛媛県今治市) の過疎地区に住む高齢 者を支える試みとして，エネルギア・コミュニケーション ズが 2016 年 10 月 26 日にドローンで野菜を配達する実験 を行った1)．島内の起点から終点に向けてドローンを海上 に迂回させて 5 分ほど飛ばし，約 500g の野菜を送り届け た．楽天がゴルフ場内のデリバリサービス用に稼働させて いるマルチコプタ型ドローン「天空」が利用された． 福島県南相馬市の，原発からの避難指示が 2016 年 7 月 に解除された地域では，楽天とローソンが 2017 年 10 月 31 日にコンビニ商品の無人配送サービス実験を行った18)． 海外では日本よりも一足先にドローンによる宅配が実現 に近づいている．以下 3 件の情報源は参考文献1)である． Amazon は 2016 年 12 月に，英国ケンブリッジ市の Amazon 配送センター近くの客 2 人に，ドローンで注文 品を届けた．運べる荷物は重量 5 ポンド (約 2.27kg) まで であり，同社は 30 分以内の配達を目標にする，と表明した． Domino’s Pizza は，2016 年にニュージーランドのオー クランド郊外でピザを注文客の裏庭に届けた． 米国ではセブンイレブンがドローンによる配送を試みて いる．ネバダ州レノの店舗から近隣住民に注文品を配達す るものであり，2017 年 1 月時点で 77 回の配送実績をもつ． ただし FAA(連邦航空局) の規制のため視界外飛行ができな いので，複数の操縦者による有視界操縦を行っている． なお，前述では宅配サービスというシナリオで説明した が，今回の定式化で追随できる配送サービスは数多ある．例 えば大規模自然災害発生後の避難者に食料や医薬品の配送を 行うときに，避難者の分布を把握し，物品のデポにドローン を搬入し用いればよい．その際に何機のドローンを配備する 必要があるかを見積もることに意味がある．また，大規模 スポーツイベントや万博等の開催時に AED の緊急配送シ ステムを構築することにも意味がある．最新型の小型 AED の重量は 1kg 程度だからドローンによる配送が十分に可能 である．大規模イベント開催時には，広いスペースが観客 や歩行者で埋め尽くされるため，自動車・バイク・自転車 が安全かつスムーズに通行できない恐れがあり，ドローン の活躍が望まれる．さらに，2020 年開催の東京オリンピッ クでは，夏季開催ゆえの熱中症の多発が懸念される．この 場合も，熱中症治療のための水や保冷剤や薬剤を現場に配 送すべく小型のドローンを活用できる可能性が十分にある． 2．定式化 都市領域内に立地する 1 つのドローン基地に配備した s 機のドローンにより，住民に宅配サービスを提供する．た だし，客からの電話やインターネットによる注文がランダ ムに入った順に品物を用意し，1 回の飛行で 1 件の配達を こなす．例えば前述の Amazon が用いたドローンの最大積 載量 5 ポンドを基準に考えると，多くの物品を積載して複 数の配達先を巡回するのではなく，個別配送を前提とする のが現実的である．ただし将来大きな積載量のドローンが 実用に耐えうる状況になれば，本研究とは異なり，巡回配 送のモデル構築が必要になるだろう． 受注時にドローンに空きがあれば直ちに品物を積み込み 出発させる．全てのドローンが出払っていれば，最初に帰 還したドローンに直ちに積み込み出発させる． 以下を定義する： a1 = (基地での品物積込み所要時間), (1) a2 = (配達先での荷下ろし所要時間), (2) a = a1+ a2, (3) τ = (基地と配達先の往復飛行時間). (4) サービス時間を t とすると t = a1+ a2+ τ = a + τ (5) である．a1ならびに a2は正の定数で与える (勿論 a も正 定数)．客の発注は強度 λ のポアソン到着と想定する． なお，客にとっての実質的なサービス時間は a + τ /2 で あり，これに発注時刻からサービス開始時刻 (その客への納 品に向けた準備が始まる時刻) までの待ち時間 w を加えた ものがリードタイム tL(発注から納品までの時間) になる． リードタイム： tL= a + τ /2 + w. (6) ドローンのバッテリ交換と品物の積込み作業は同時並行 とし，かつバッテリ交換時間は品物積込み所要時間 a1以下 と想定する．さらに，必要な容量のバッテリが即時的に供 給できる体制が整っているものとし，バッテリ充電のため にドローンの使用開始が遅れるという事態は考えない．この とき，サービス時間はバッテリ交換時間に依存しない．こ れらを確率変数で与えて分析することも勿論可能だが，本 研究では無用の複雑化を排除すべく簡便化した． ここでドローン基地から顧客への距離を x とし f (x) = (x の確率密度関数) (xmin≤ x ≤ xmax) (7) を定義する (xminならびに xmaxは距離 x の最小値ならび に最大値)．ドローンの往復飛行距離は 2x である．サービ ス時間分布を g(t) として，これを求めておこう． ドローンの移動速度を v で与えると，飛行時間は τ = 2x/v だから，式 (5) から t = a + 2x/v (8) である． dx = v 2 dt に着目してサービス時間分布 g(t) を 求めると，次式の通りである： g(t) = v 2f (v 2(t− a) ) (9) ( a +2xmin v ≤ t ≤ a + 2xmax v ) このとき，サービス時間の期待値⟨t⟩，分散 σ2 t ならびに変 動係数 c は片道飛行距離 x の特性値によって表される： ⟨t⟩ = a + 2 v⟨x⟩, (10) σt2 = ( 2 v )2( ⟨x2⟩ − ⟨x⟩2), (11) c = σt ⟨t⟩. (12)
上で，σtは t の標準偏差である．また⟨x⟩ ならびに ⟨x2⟩ は，片道飛行距離 x の平均値と 2 乗の平均値を表す． 以上要するに，このシステムは窓口数が s，到着率が λ であって，サービス率 µ とトラフィック密度 ρ が µ = 1 ⟨t⟩, (13) ρ = λ sµ (14) の M/G/s 待ち行列系である．M/G/s 待ち行列系とは，客 の到着間隔が指数分布に従い (つまりポアソン到着)，サー ビス時間が一般分布に従う，窓口数 s の待ち行列モデルを 意味する5)．客は 1 本の行列をなし，待ちの先頭から順に 空き窓口 (今回のモデルでは空きドローン) に割り当てられ る．今回は待ち行列長に上限値がないものとする． 3．木村の平均待ち時間近似式の導入 扱い易い M/M/s 待ち行列系 (客の到着間隔とサービス 時間がそれぞれ指数分布に従う系5)) とは異なり，M/G/s 待ち行列系については平均待ち時間の厳密解 EW (M/G/s) が準備されていない．そこで今回は，この分野の最先端の研 究成果である高精度近似公式群3,4,5)から，木村の 2 モー メント近似という平均待ち時間近似公式5)を導入する： EW (M/G/s)∼ 1 + c 2 2c2+ 1− c 2 R(s, ρ) EW (M/M/s). (15) ここで c はサービス時間 t の変動係数 (12) であり，ρ はト ラフィック密度 (14) である．そして，EW (M/M/s) は問 題にしている M/G/s 待ち行列系と同じ窓口数，平均サー ビス時間ならびにトラフィック密度をもつ M/M/s 待ち行 列系の平均待ち時間であり，その解析解は周知の通り EW (M/M/s) = (sρ) s s!sµ(1− ρ)2 {s−1 ∑ k=0 (sρ)k k! + (sρ)s s!(1− ρ) }−1 (16) である．また，式 (15) の R(s, ρ) は次の様に定義される： R(s, ρ) = 1 2{1 + ϕ(s, ρ)ψ(s, ρ)}, (17) ϕ(s, ρ) = (1− ρ)(s − 1)( √ 4 + 5s− 2) 16sρ , (18) ψ(s, ρ) = 1− exp [ −(s − 1) (s + 1)ϕ(s, ρ) ] . (19) 近似公式 (15) は，M/G/1，M/M/s，M/D/s 待ち行列 に対しては厳密解を与え (M/D/s はポアソン到着でサー ビス時間が一定値であることを示す5))，s → ∞ または ρ→ 1 のとき正確な漸近的性質を持つ．加えて殆どの (s, ρ) ペアに対して，相対誤差にして 1%未満の近似を与える4)． 前述の行列をなす人数を待ち客数と呼び，これにサービ スを受けつつある人数を加えたものを系内客数と呼ぶ．待ち 客数と系内客数の平均値を各々ELqならびに EL と表す． 今回の平均待ち客数 ELq(M/G/s) と平均系内客数 EL(M/G/s) は，周知の Little 公式により次で与えられる： ELq(M/G/s) = λEW (M/G/s), (20) EL(M/G/s) = ELq(M/G/s) + sρ. (21) さらに待ち確率 Π(客が待たねばならない確率) は，木村4) が推奨するように M/M/s 型の待ち確率で近似する： Π∼ (sρ) s s!µ(1− ρ) {s−1 ∑ k=0 (sρ)k k! + (sρ)s s!(1− ρ) }−1 . (22) 4．円盤領域で一様な需要のモデル化 4.1 距離分布と特性値の準備 半径 R の円盤領域で需要点 P が一様に発生する (図 1)． そして円盤の中心点から距離 h の地点にドローン基地 F を 立地させる．ドローンの片道飛行距離は x = FP である． R R F F h h P P x x (i) h≤ R (ii) R < h O O 図 1 円盤上で一様な需要点 P とドローン基地 F の距離 x. 距離 x の確率密度関数を f (x) とすると，(i)h≤ R の場 合と (ii)R < h の場合に分けて次の様に記述される7,9,10)： (i) h≤ R の場合 f (x) = 2x R2 (0≤ x ≤ R − h) 2x πR2arccos x2+ h2− R2 2hx (R− h < x ≤ R + h) (23) (ii) R < h の場合 f (x) = 2x πR2arccos x2+ h2− R2 2hx (24) (h− R < x ≤ h + R). また，距離 x の累積分布関数を求めると次の通りである (累積分布関数は後にドローン待ち行列の確率的シミュレー ションを行うに際して乱数を発生させるために用いる)： (i) h≤ R の場合 F (x) = x2 R2 (0≤ x ≤ R − h) 1 2+ x 2πR2 { 2x arccosh 2− R2+ x2 2hx (25) −1 x ( δ + 2R2arctanh 2+ R2− x2 δ )} (R− h < x ≤ R + h)
(ii) R < h の場合 F (x) = 1 2+ x 2πR2 { 2x arccosh 2− R2+ x2 2hx −1 x ( δ + 2R2arctanh 2+ R2− x2 δ )} (26) (h− R < x ≤ h + R). ただし次を定義する： δ =√(−h + R + x)(h + R − x)(h − R + x)(h + R + x). (27) x の平均値は次の通りである7,9,10)： (i) h≤ R の場合 ⟨x⟩ = 4R 9π [{ 7 + ( h R )2} E(h R) −4 { 1− ( h R )2} K(h R) ] , (28) (ii) R < h の場合 ⟨x⟩ = 4h 9π [{ 7 + ( h R )2} E(R h) − { 2 + ( h R )2 − 3 ( R h )2} K(R h) ] . (29) ただし，K(k) は第 1 種の完全楕円積分であり，E(k) は第 2 種の完全楕円積分である13)： K(k) = ∫ π/2 0 dθ √ 1− k2sin2θ, (30) E(k) = ∫ π/2 0 √ 1− k2sin2θ dθ. (31) このように，結果が初等関数では表せないために，平均距 離の厳密な値を求めるには，数値計算が必要となる．ただ し，これらには高精度の近似公式が用意されている10,11)： (i) h < R の場合 ⟨x⟩ ∼ 2 3R + h2 2R− h4 32R3, (32) (ii) R≤ h の場合 ⟨x⟩ ∼ h +R 2 8h + R4 192h3. (33) 式 (32) の近似は正の偏差をもち，相対偏差は h = 0 のとき に 0% であって，h の増加につれて増加するものの，h = R のときでさえ約 0.32%に過ぎない．式 (33) の近似は負の 偏差をもち，相対偏差は h = R のときでさえ−0.14%に 過ぎず，h の増加につれて速やかにゼロに近づく． 一方，x の 2 乗の平均値は偏心距離 h と半径 R の大小 関係に関する場合分けを必要とせず，次の通りである： ⟨x2⟩ =R2 2 + h 2. (34) 4.2 数値例による待ち行列の計算 ここでの設定は，半径 R の円盤領域で宅配利用客が一様 に分布しており，ランダムに需要が発生するというもので ある．強度 λ のポアソン到着に s 機のドローンで対応する． 以下の数値実験では，サービス対象領域の半径を R = 3km に固定し，まず施設立地場所 h，配送品の積込み・荷 下ろし時間 a，ドローン飛行速度 v，注文の到着率 λ を表 1 の通りに様々に設定し，各場合に平均待ち時間がドロー ン配備数 s に依存して変化する様を観察する (4.2.1 節)．h はドローン基地建設時の意思決定に対応し，a と v はサー ビス提供者の運営努力やドローンとそれを取り巻くシステ ムの技術水準に依存する．そして λ はサービス対象地域に おける人口のあり様に対応する指標である． 表 1 のパラメタの設定範囲について解説する．立地点 h は 円盤中心から外部に至る 0km∼6km で与えた．積込み・荷下 ろし時間 a は基準値を 0.15h とし，常識的範囲として 0.05h ∼0.25h とした．ドローン飛行速度 v は基準値を 50km/h とし (現状の社会実験での実績は 30∼60km/h)，技術向上 により速くなる可能性を加味して，25km/h∼100km/h で 与えた．到着率 λ は基準値を 30 件 /h とし，15 件 /h∼ 60 件 /h で与えた．半径 R = 3km の領域で，例えば買い 物弱者の人口密度を 30 人 /km2とし，世帯当たり人数を 2 人とすると，(買い物弱者世帯数) = π· 32· 30/2 ≃ 425 である．仮に 2 日に 1 回の割合で発注すれば，(1 日の発注 数) = 425/2≃ 210 である．ここでドローン基地での受付 け時間を 10:00 から 16:00 までの 6 時間とし，発注分布が 一様であると仮定すると，到着率が λ = 210/6 = 35 件/h となる． ドローン基地の設計では，サービス水準確保に必要なド ローン配備数の見積もりが最重要である．そこで，平均待 ち時間をある閾値以下にするために必要なドローン配備数 に焦点を当てる (4.2.1 節)．続いて，同様の分析を平均待ち 客数ならびに平均系内客数に焦点を当てて行う (4.2.2 節)． ドローンによる宅配システムが多面的に見積もられる． 表 1 数値実験におけるパラメタ設定の一覧． パラメタ 設定値 対応する (太字は基準値) 実験 偏心距離 0, 1, 2, 3, 4, 5, 6 (a),(e) h［km］ 積込み・荷下ろし 0.05, 0.10, 0.15, (b) 時間 a［h］ 0.20, 0.25 ドローン速度 25, 50, 75, 100 (c) v［km/h］ 注文の到着率 15, 30, 45, 60 (d) λ［件 /h］ なお，複数窓口の待ち行列系が意味を持つためには，ト ラフィック密度 ρ が 1 未満でなければならない．ρ が下か ら 1 に近づくと，待ち時間や行列長の期待値が無限大に発 散してしまい，そのようなサービスに現実的な意味はない．
式 (14) の ρ を 1 未満にするための s の最小値を sminとし て求めると，次式の通りである (λ/µ 以上の最小の整数)： smin= ⌈ λ µ ⌉ . (35) ドローン配備数 s の定義域は smin以上の整数集合である． 4.2.1 平均待ち時間 4.1 節の距離の平均値⟨x⟩ と 2 乗の平均値 ⟨x2⟩ を 2 章 の式 (10) と (11) に代入すれば，平均サービス時間⟨t⟩ と 標準偏差 σtを得る．そして式 (13) と (14) でサービス率 µ とトラフィック密度 ρ が準備される．それらを式 (16)∼ (19) に代入し，さらに式 (15) の右辺に代入すれば，待ち 行列の平均待ち時間 (簡単に EW と記す) の理論値を得る． 平均距離⟨x⟩ の計算に当たり，完全楕円積分の数値計算 が即座に実行できるのなら式 (28)，(29) を用いればよい． Fortran や C といった技術計算向きの高級言語であれば完 全楕円積分のサブルーチンが容易に利用できるし，Mathe-matica，Matlab，Maple 等の高度な数式処理システムは 完全楕円積分を内生関数として持つ．本研究では，Mathe-matica Ver. 11.1.1.0 を用いた．一方，完全楕円積分のサ ブルーチンや内生関数が容易に利用できない場合は，近似 式 (32) と (33) によって遜色のない計算が可能となる． (a) 立地場所選択 h の影響 まずは，(a, v, λ) = (0.15h, 50km/h, 30 件/h) とする (基準値)．そして，ドローン基地の位置 h を表 1 の 2 行目 の様に h = 0km から h = 6km まで，1km 刻みで変化さ せる (図 2)．土地利用上の制約がなければ，勿論 h = 0km とするのが得策である (平均飛行距離が最小になる)．それ が許されないときの系のパフォーマンスを観察したい． 各 h 値毎にドローン配備数 s を動かして，平均待ち時間 を算出した (図 3)．例えば，h = 1km の場合，ドローン 8 機のとき平均待ち時間は約 5 分 30 秒である．もしも平均 待ち時間を 1 分未満 (図 3 中の 1 分の基準線より下) にし たければ，ドローンが 10 機必要である． 図 3 からは，平均待ち時間の削減に基地の立地選定がと ても重要であることがわかる．例えば h = 4km とすると， 平均待ち時間を 1 分未満にするのに 13 機のドローンが必 要となる．もしも h = 0km とすれば，10 機のドローンに よって平均待ち時間を約 30 秒にできるのである． どの場合も，配備数の増加につれて平均待ち時間が減少 し，1 機増やす減少効果は s が大きくなるにつれ逓減する． R = 3km h = 0km h = 6km 図 2 半径 3km の対象領域におけるドローン基地の配置. 7 8 9 10 11 12 13 14 0 1 2 3 4 5 6 h=0km h= 1km h=2km h=3km h=4km h=5km h=6km ドローン配備数 s [機] 平 均 待 ち 時 間 E W [分 ] 図 3 ドローン基地の立地点 h とドローン配備数 s が待ち 行列における平均待ち時間 EW に与える影響 ［パラメタ値は (a, v, λ) = (0.15h, 50km/h, 30 件/h)］． (b) 積込み・荷下ろし時間 a の影響 次に (h, v, λ) = (1km, 50km/h, 30 件/h) とする．積込 み・荷下ろし時間 a を表 1 の 3 行目のように 0.05h，0.10h， 0.15h，0.20h，0.25h と変化させる．a が小さいほど，現 場での配送品準備に関する作業効率が改善され，またドロー ン自体も積込み・荷下ろしがスムーズに運ぶように進化し ている．つまり a 値はサービス提供者の現場での改善努力 とドローンの技術革新を表す指標である． 各 a 値毎にドローン配備数 s を動かし，平均待ち時間を 求めた (図 4)．例えば a = 0.15h(つまり 9 分) の場合，ド ローン 8 機で平均待ち時間を約 5 分 30 秒にできる．平均待 ち時間を 1 分未満にするにはドローンが 10 機必要である． 積込み・荷下ろし時間を削減するための努力がとても重 要であることも読み取れる．例えば a = 0.20h(つまり 12 分) のとき，平均待ち時間を 1 分未満にするためには 12 機 のドローンが必要となる．もしも店舗内での合理化やドロー ンの技術革新によって a = 0.10h(つまり 6 分) に改善でき れば，たった 8 機のドローンで平均待ち時間を 43 秒にま で削減することが可能なのである． 5 6 7 8 9 10 11 12 0 1 2 3 4 5 6 a=0.05h a=0.10h a=0.15h a=0.20h a=0.25h ドローン配備数 s [機] 平 均 待 ち 時 間 E W [分 ] 図 4 配送時の積込み・荷下ろし時間 a とドローン配備数 s が待ち行列における平均待ち時間 EW に与える影響 ［パラメタ値は (h, v, λ) = (1km, 50km/h, 30 件/h)］．
(c) ドローン飛行速度 v の影響 今度は (h, a, λ) = (1km, 0.15h, 30 件/h) とする．そし てドローンの飛行速度 v を表 1 の 4 行目のように 25km/h， 50km/h，75km/h，100km/h と変化させる． 各 v 値毎にドローン配備数 s を動かし，平均待ち時間を 求めた (図 5)．例えば，v = 50km/h のときに，ドローン 8 機で平均待ち時間を約 5 分 30 秒にできる．平均待ち時間 を 1 分未満にするにはドローンが 10 機必要である． 飛行速度が増せばドローンを削減できる．例えば v = 50km/h のとき，平均待ち時間を 1 分未満にするためにド ローン 10 機が必要である．技術革新や法制度の改定により v = 100km/h となれば，ドローン 8 機で平均待ち時間を 54 秒にできる．ただし，こうした効果は圏域の大きさ R や 積込み・荷下ろし時間 a にも依存する．R が大きいほど，a が小さいほど，飛行速度 v の増大による効果は大きい． 6 7 8 9 10 11 12 0 1 2 3 4 5 6 v=100km/h v=75km/h v=50km/h v=25km/h ドローン配備数 s [機] 平 均 待 ち 時 間 E W [分 ] 図 5 ドローン飛行速度 v とドローン配備数 s が待ち行列 における平均待ち時間 EW に与える影響 ［パラメタ値は (h, a, λ) = (1km, 0.15h, 30 件/h)］． (d) 注文の到着率 λ の影響 外生的要因である注文の到着率 λ が平均待ち時間に与え る影響も観察せねばならない．需要量は将来の人口増減に依 存する．現状の到着率 λ を正しく見積もって適切なドロー ン数を決めたとしても，将来 λ が低下すれば，ビジネスモ デルとして成り立たなくなる恐れがある．逆に人口増加に よる λ の増大が予見されるなら，ドローンの増設計画が必 要となる．ドローンの発着には自ずと面積を必要とするか ら，ドローン基地に予めバッファ面積を確保しておかねば ならない．いわば成長・衰退管理のための分析と言える． ここでは (h, a, v) = (1km, 0.15h, 50km/h) とする．そ して発注の到着率 λ を表 1 の 5 行目のように 15 件/h， 30 件/h，45 件/h，60 件/h と変化させる． 各 λ 値毎にドローン配備数 s を動かし，平均待ち時間を 求めた (図 6)．例えば λ = 30 件/h ならば，ドローン 8 機 で平均待ち時間を約 5 分 30 秒にできる．もしも平均待ち時 間を 1 分未満にしたければ，ドローン 10 機にすればよい． もしも地域人口が減少し，将来のある時点で λ = 15 件 /h になれば，平均待ち時間を 1 分未満にするのに必要なド ローン数は 6 機になる (削減数は 4 機)．当然必要なスタッ フの数も減少する．そのときに経済性工学的にドローン基 地の運営が可能であるかどうか (すなわち損益分岐点をクリ アできるかどうか)，慎重に議論する必要がでてくるだろう． 逆に将来のある時点で λ = 60 件/h になると見積もられ たとしよう．そのときには平均待ち時間を 1 分未満にする ために 18 機のドローンが必要となる．最初からそれを見越 して，増設すべき 8 機の発着に必要な面積を，基地施設の 中に用意しておかねばならない．そして，λ が現実に増加 するのに対応して，配備数を漸増させていけばよい．ただ し，何年も経つうちにはドローンそのもののスペックも上 がり，飛行速度 v の増大や積込み・荷下ろし時間 a の減少 によって，ある程度はカバーできるかもしれない．そうし たことを多面的に考慮する必要もあるだろう． 5 10 15 20 0 1 2 3 4 5 6 λ=15 件/h λ=30 件/h λ=45 件/h λ=60 件/h ドローン配備数 s [機] 平 均 待 ち 時 間 E W [分 ] 図 6 発注の到着率 λ とドローン配備数 s が待ち行列にお ける平均待ち時間 EW に与える影響 ［パラメタ値は (h, a, v) = (1km, 0.15h, 50km/h)］． (e) 一定の平均待ち時間 EW を達成する (h, s) 対 立地点 h とドローン配備数 s が平均待ち時間 EW に与 える影響をさらに詳細に見るために，一定の平均待ち時間 をもたらす変数対 (h, s) を求めた (図 7)．図 7 では平均 待ち時間を 0.5 分，1 分，2 分，5 分の 4 通りに設定し， これらを達成するためのドローン配備数 s が基地の偏心距 離 h に依存する様子を示している．パラメタは (a, v, λ) = (0.15h, 50km/h, 30 件/h) とおいた．これを見ると，一定 0 1 2 3 4 5 6 8 10 12 14 16 EW =0.5 分 EW =1 分 EW =2 分 EW =5 分 円盤中心からの偏心距離 h [km] 必 要 な ド ロ ー ン の 数 s [機 ] 図 7 一定の平均待ち時間 EW 達成に必要な (h, s) 対 ［パラメタ値は (a, v, λ) = (0.15h, 50km/h, 30 件/h)］．
の偏心距離 (この数値例では 2km ほど) 以上の場合は，あ る平均待ち時間を達成するために必要なドローンの数 s が 偏心距離 h のほぼ線形で増加する性質が読み取れる． 以上要するに，平均待ち時間の観点から，基地の立地選 定 h，ドローンの飛行速度 v，積込み・荷下ろし時間 a，到 着率 λ が重要な役割を果たすことが示された． 4.2.2 平均待ち客数，平均系内客数と待ち確率 サービス水準を見積もるには，平均待ち客数，平均系内 8 10 12 14 0 10 20 30 40 h=0km h =1km h =2km h=3km h=4km h=5km h=6km ドローン配備数 s [機] 平 均 待 ち 客 数 E Lq [人 ] 図 8 基地立地点 h と配備数 s が平均待ち客数 ELqに与 える影響 ［(a, v, λ) = (0.15h, 50km/h, 30 件/h)］． 8 10 12 14 0 10 20 30 40 h=0km h =1km h =2km h =3km h=4km h=5km h=6km ドローン配備数 s [機] 平 均 系 内 客 数 E L [人 ] 図 9 基地立地点 h と配備数 s が平均系内客数 EL に与 える影響 ［(a, v, λ) = (0.15h, 50km/h, 30 件/h)］． 8 10 12 14 0.0 0.1 0.2 0.3 0.4 h=0km h =1km h=2kmh=3km h=4km h=5km h=6km ドローン配備数 s [機] 待 ち 確 率 Π 図 10 基地立地点 h と配備数 s が待ち確率 Π に与える影 響 ［(a, v, λ) = (0.15h, 50km/h, 30 件/h)］． 客数ならびに待ち確率も重要である．これらがドローン配 備数 s に依存する様も観察しておく (図 8，図 9，図 10)． これらを参照すれば，パラメタを適切に与えた上で，これ らの指標をある水準に保つような立地場所 h とドローン配 備数 s を決定するための基礎資料とすることができる． 前節の平均待ち時間の情報と合わせて，プロジェクトの 予算が許す範囲内で要求されるサービス水準を多面的に考 慮しつつドローン基地の設計を行うことが可能である． 4.3 確率的シミュレーションの適用 ここまで平均値に注目して数値例を与えてきた．それら はドローン基地の青写真を描くときの最初の一太刀として 利用できる．しかし，サービス実現に向けた更に詳細な分 析にも意味がある．そのためには次の方法でサービス時間 の乱数を生成して確率的シミュレーションを行えばよい． まず，[0, 1] の一様乱数 u を発生させる．そして，片道 飛行距離 x の分布関数 F (x)(式 (25) と (26)) に基づき u = F (x) (36) を x について数値的に解き，x = x∗を求める．それを式 (8) の x に代入して，サービス時間の乱数 t∗= a + 2x∗/v を発生させる．これは周知の逆関数法である． 今回は NTT データ数理システム社の汎用シミュレー ションソフトウェアである S4 Simulation System(エス クワトロ) Ver. 5 を用いた．紙面の制約から，パラメタ を (h, a, v, λ, s) = (1km, 0.15h, 50km/h, 30 件/h, 9 機) としたときの待ち時間分布を図 11 に例示するにとどめる． ドローン基地の開設時間を 10:00∼16:00(6 時間) としてシ ミュレーションを行った (1 つの) 算出結果である．シミュ レーションを丁寧に行えば，ドローンの性能・立地場所・ 到着率等に応じた現実的な見積もりが可能となる． この設定下での理論的な平均待ち時間は EW = 1 分 40 秒 である．図 11 の数値例における平均待ち時間実 測値は ¯w = 1 分 51 秒 である．理論値 EW は待ち行列の 平衡状態を想定した解であり，もとより有限時間の一出力と は異なるが，十分に意味のある参照値であることがわかる． 0 2 4 6 8 10 12 0 20 40 60 80 100 120 平均値： ¯w = 1 分 51 秒 標準偏差：σw= 1 分 46 秒 待ち時間 w [分] 観 測 数 [件 ] 図 11 シミレーションによる待ち時間の度数分布 ［(h, a, v, λ, s) = (1km, 0.15h, 50km/h, 30 件/h, 9 機)］．
5．離散的な需要分布への対応 対象地域を近似的に需要が一様と見做せる (単一あるい は複数の) 円盤で表せば，待ち行列の理論的特性値やシミュ レーションを通じてシステムの設計を支援できる (前述)． しかし，離散的な需要分布で計算したい場合もある． 例えば対象領域の n 個の町丁目に 1 つずつドローンの荷 下ろし場があるものとする．町丁目 i の一日当たり発注数を qiとし，ドローン基地から町丁目 i までの片道飛行距離を xiとする．領域全体での一日当たり発注数を Q とすると Q = q1+ q2+· · · + qn (37) である．受注の時間を T とすれば，到着率は λ = Q/T と 準備できる．片道飛行距離の平均値と 2 乗の平均値は ⟨x⟩ = 1 Q n ∑ i=1 qixi, ⟨x2⟩ = 1 Q n ∑ i=1 qix2i (38) だから，これらを式 (10) と (11) に代入して直ちに 4 章と同 様の計算ができる．確率的シミュレーションも容易である． 6．結語と提言 円盤領域で一様な需要分布を中心に据え，ドローンによ る宅配の待ち行列現象を議論した．そして実は片道飛行距 離の分布さえ与えれば同様の分析ができ，さらには確率的 シミュレーションも直ちにできることを述べた．理論的な 参照値を用いつつ，詳細はシミュレーションによればよい． ドローン基地に近い利用者は直接訪れて買い物をすると 想定したければ，距離 x の下限値 xmin(> 0) を与え，xmin でトランケイトした距離分布に基づいて計算すればよい．数 値計算の手間は増加するが，追求することが勿論可能である． 4 章の数値実験から分かるように，ドローンの望ましい 配備数はドローンのスペックと対象地域の需要量 (到着率) に依存する．今後，買い物弱者問題が顕在化した人口低密 地域に設定されたドローン特区で，宅配パイロット事業が 展開されるであろう．その際，本研究の都市解析の成果を 計画立案に役立てることが可能である．特に，ドローン配 送システムを持続可能にするためには，ドローンのスペッ ク向上と地域における需要の成長・衰退を共に視野に入れ ねばならない (4.2.1 項の (d) で述べた)．将来に向けた地 域人口の推移を予測し，そうした変化が起きても配送シス テムが経済性工学的に成り立つように，事前の準備 (需要が 減る場合は事業縮小に向けた準備，需要が増える場合はド ローン増設の施設内バッファ確保) をしておくべきである． 本研究の手法はそうした見積もりの基礎資料を与え得る． さらに付言すれば，ドローン宅配システムに都市計画上 の副産物があることを強く認識すべきである．一例を述べる と，安全に飛行できるドローンのインフラストラクチャが 整っていれば，大規模災害発生時の空撮による情報収集や 行方不明者の探索，そして被災住民への緊急物資の輸送に も，比較的容易に転用できる筈である (定常的にドローン・ システムに携わっているスタッフの存在も心強い)．本研究 が都市計画分野におけるドローン・インフラストラクチャ の計画立案に資するところがあれば誠に幸いである 謝辞 匿名の査読者から論文の質を高めるための改善案に加え て，本研究のモデル分析に基づく計画論や現行制度への示 唆も記述すべきであるという，懇切でポジティブなアドバ イスを頂戴しました．ここに深甚なる謝意を表します． 参考文献 1) ダイヤモンドオンラインの記事 (2017 年 1 月 12 日)『ドロー ン宅配，日本でも実験開始！離島暮らしの救世主になるか』 ， [http://diamond.jp/articles/-/112824]，2018 年 4 月 25 日 閲覧可． 2) 河 鐘基 (2015)：『ドローンの衝撃』，扶桑社． 3) 木村俊一 (1988a)：M/D/s 待ち行列の平均待ち時間に対する近似， 電子情報通信学会論文誌 A，Vol.71–A，No.3，pp.696–702. 4) 木村俊一 (1988b)：M/G/s 待ち行列の近似式の有効性について， オペレーションズ・リサーチ，1988 年 5 月号，pp. 20–22． 5) 木村俊一 (2016)：『待ち行列の数理モデル』，朝倉書店． 6) 小林啓倫 (2015)：『ドローン・ビジネスの衝撃―小型無人飛行機が 切り開く新たなマーケット―』，朝日新聞出版． 7) 腰塚武志 (1986)：都市平面における距離の分布，『都市計画数理』(谷 村秀彦ほか)，朝倉書店，pp. 1–55． 8) 古藤 浩 (2014)：山形県のメッシュデータによる買い物困難地域の 分析，東北芸術工科大学紀要，No.21，pp. 116–124． 9) 栗田 治 (2004)：『都市モデル読本』，共立出版． 10) 栗田 治 (2013)：『都市と地域の数理モデル―都市解析における数学 的方法―』，共立出版． 11) 栗田 治，腰塚武志 (1989)：周上で一様な点に関する平均値のある 導出法とその応用，日本オペレーションズ・リサーチ学会秋季研究 発表会アブストラクト集，1–B–2，pp. 38–39． 12) 三浦英俊，古藤 浩 (2010)：メッシュデータを用いた人口減少地域 における買い物距離の分析―山形県における食料品店を事例として―， 都市計画論文集，Vol. 45, No. 3, pp. 643–648． 13) 森口繁一他 (1956)：『数学公式 I』(岩波全書)，岩波書店． 14) 森・濱田松本法律事務所 ロボット法研究会編 (2017)：『ドローン・ ビジネスと法規制』，清文社． 15) 野波健蔵 (2017)：ドローン技術の現状と課題およびビジネス最前線， 情報管理，Vol. 59，No. 11，pp. 755–763． 16) 関口大介，岩崎覚史 (2017)：『ドローンビジネス参入ガイド』，翔 泳社． 17) 鳥海重喜 (2014)：福岡市におけるフードデザート問題の分析，都市 計画論文集，Vol. 49，No. 3，pp. 993―998． 18) Yomiuri Online の 記 事 (2017 年 11 月 01 日)： 『 コ ン ビ ニ 商 品 、ド ロ ー ン で 配 送 … 福 島 で 実 証 実 験 』 ， [http://www.yomiuri.co.jp/feature/TO000303/ 20171101-OYT1T50018.html]，2018 年 4 月 25 日閲覧可． 19) Boutilier, J.J., S.C. Brooks, A. Janmohamed, A. Byers,
J.E. Buick, C. Zhan, A.P. Schoellig, S. Cheskes, L.J. Mor-rison and T.C.Y. Chan(2017): Optimizing a drone net-work to deliver automated external defibrillators,
Circu-lation, Vol. 135, pp. 2454–2465.
20) Goodchild, A. and J. Toy(2018): Delivery by drone: An evaluation of unmanned aerial vehicle technology in reduc-ing CO2emissions in the delivery service industry,
Trans-portation Research Part D, Vol. 61, pp. 58–67.
21) Grippa, P., D.A. Behrens, F. Wall and C. Bettstet-ter(2018): Drone delivery systems: Job assignment and dimensioning, Autonomous Robots, Vol. 18 [Published online].
22) Hong, I., M. Kuby and A.T. Murray(2018): A range-restricted recharging station coverage model for drone de-livery service planning, Transportation Research Part C, Vol. 90, pp. 198–212.
23) Murray, C.C. and A.G. Chu(2015): The flying side-kick travelling salesman problem: Optimization of drone-assisted parcel delivery, Transportation Research Part C, Vol. 54, pp. 86–109.
