ネットワークモデルの基礎
利用者均衡配分井料隆雅
東京大学工学系研究科（客員）東北大学情報科学研究科（本務）18 Sep 2021

Takamasa Iryo

1

講師紹介

2

• 井料隆雅（いりょう たかまさ）– 東京大学大学院工学系研究科
特定客員大講座担当教授（兼務：今日はこっち）– 東北大学大学院情報科学研究科
人間社会情報科学専攻 教授（本務：日頃はこれ）– 神戸大学大学院工学研究科客員教授（たまに）• R2/3まで：神戸大学大学院工学研究科市民工学専攻 教授

18 Sep 2021

Takamasa Iryo

交通量配分3• 所与の起終点交通量を所与のネットワーク上に割り当てる（＝配分する）• 交通需要を交通ネットワークがどのように処理するのかを計算する．
• さまざまな施策の効果を評価できる．
– 最適化問題の制約条件にすることにより（MPEC/Bi-level）最適な施策を計算できる．

18 Sep 2021

Takamasa Iryo

交通量配分の計算のポイント• ドライバー（利用者）の行動は交通状況に依存する．
• 交通状況はドライバーの行動に依存する．
• ∴ ドライバーの行動と交通状況は相互作用する．
– ドライバー同士の行動が相互作用する，ともいえる．
この「相互作用」をどうモデルかし，どう解くかが交通量配分の最大のポイントである．
※ 相互作用があるシステムでは，需要と交通行動がわかっても，正しい将来予測や施策評価はできない．

18 Sep 2021

Takamasa Iryo4

相互作用の逐次計算

5• 交通シミュレーションの例– 1台1台の車に，ネットワーク上での挙動のルールを与えて，そのとおりに動かす．
– 挙動のルールが「過去から今現在までの交通状況のみに依存する」であれば，各車両の挙動を，時間軸方向に沿って，逐次計算すればよい．
• ある時間までの車両の挙動の計算結果が，その時間における各車両の挙動を決定し，それによって，直後の時間までの交通状況を計算できる．
• 数学的に言えば常微分方程式を解いていることになる．

18 Sep 2021

Takamasa Iryo

ドライバーの将来予測？

6• ドライバーは，過去の経験に基づいて，交通状況の将来予測が可能であろう．
– 頻繁に道路を使っていれば，どの経路がより早く目的地に着くかについて，経験に基づく戦略を持つことになる，と考えられる．
– 予測を強いられることもある：到着時刻制約があれば，所要時間を予測して出発時刻を決めなくてはならない．

• 将来予測が可能である／必要とすれば，時系列に沿った逐次的な計算はできなくなる．

18 Sep 2021

Takamasa Iryo

Day-to-dayの繰り返し計算

7• 繰り返し道路を使用するドライバーの学習過程もシミュレーションしてしまえばよい？
– 複数日のシミュレーションを繰り返す．
– ドライバーは，過去の日の混雑状況を参照する．
– 今日≒昨日としてもよいし，学習モデルを考えてもよい．

• この考え方自体は正しい．
• 繰り返し計算を行った「結果」だけを先回りして計算することはできないか？

18 Sep 2021

Takamasa Iryo

均衡状態

8

• どのドライバーも，他のドライバーの行動が不変であれば，現在の選択よりもよい選択肢を持たない．
– この状態が一旦成立すれば，これ以降，どのドライバー
も，選択肢を変える動機付けを持たないので，その状態が長期的に安定して実現することが”期待”できる．

• 均衡状態– 上記の定義は「Nash均衡」に相当する．
– Wardropの第1原理も概ね同じ意味（厳密には異なる）– 均衡状態を条件とする交通量配分：利用者均衡配分

18 Sep 2021

Takamasa Iryo

均衡状態の計算

9• Day-to-dayの繰り返し計算をそのまま実行– もっとも素朴な方法– 計算時間がかかる．

• 不動点問題で定式化し，それに対する解法を適用– 比較的多くの一般的な問題に適用可能• 等価最適化問題で定式化する– もっとも利便性が高い．
– 適用可能な問題は限られる．

18 Sep 2021

Takamasa Iryo

最適化問題での定式化

10• 静的交通量配分– 交通流を連続流として扱い，時間変化を考慮しない．

• リンク旅行時間は，リンク交通量の単調増加関数として示せる．BPR関数が有名．
• 均衡状態を示す条件と等価な最適化問題が存在

18 Sep 2021

Takamasa Iryo

等価最適化問題
この最適化問題のKKT条件が，Wardropの第1原理と
等価となる
OD交通量
リンク交通量11xlMin.å ò tl ( w)dwlÎE0

subject to:

リンク旅行時間

qi = å f r for all i ÎWrÎRi

xl = åå elr f r for all l Î E
iÎW rÎRi

経路交通量

f r ³ 0 for all r Î Ri , i ÎW

Beckmann et al. 1956 イェール大のWebサイトで全文公開されている
https://cowles.yale.edu/sites/default/files/files/pub/misc/specpub-beckmann-mcguire-winsten.pdf

18 Sep 2021

Takamasa Iryo

等価最適化問題の解釈• 目的関数は「ポテンシャル関数」と解釈できる．
– 勾配（grad，∇）をとったら何が出てくるか？

• ポテンシャル関数が最小となる点が均衡状態– 交通量保存則を制約条件とした最適化問題• ゲーム理論では「ポテンシャルゲーム」と呼ぶ．
– いろいろと性質がよい．
– 均衡状態の「安定性」が保証される
Day-to-dayの繰り返し計算が必ず収束する．
（例えば，Sandholmの教科書（2010））を参照．
18 Sep 2021

Takamasa Iryo

12

最適化問題の解法

13• Frank-Wolfe法を用いる方法が一般的– 最速ではないが，実用性があり，簡単に実装できる．
– LeBlanc et al. (1975)• 最短経路探索（一般にはDijkstra法を用いる）と1次元探索（黄金分割法が多用される）を組み合わせる．
– 利用者均衡配分問題では，FW法のLP近似問題の部分を，最短経路探索問題に置き換えられる．
これによって計算が高速化されている．

18 Sep 2021

Takamasa Iryo

Frank-Wolfe法

14• xn :n回目の繰り返し計算でのリンク交通量パターン• Step 1: All-or-nothing配分– OD間の最短経路を計算し，各ODのOD交通量の全量をその最短経路に配分し，そのときのリンク交通量パターン yを計算する．
• 最短経路探索にはDijkstra法を用いることが多い．Dijkstra法は1つの起点から全ノードまでの最短経路をまとめて計算するので，
All-or-nothing配分は，起点ごとにまとめて計算するのがよい．
• 最短経路を記憶する必要がないことに注意！ 各ODのOD交通量を
各リンクに配分したら，そのODの最短経路は忘れて良い．

18 Sep 2021

Takamasa Iryo

Frank-Wolfe法

15• Step 2:1次元探索– xn+1＝(1-a) xn+ ay とし，最適化問題の目的関数
Z(xn+1)を最小化する 0以上1以下のa を計算する．
– 1次元探索の方法は何でも構わないが，黄金分割法
を用いれば効率的に計算できる．

Step1と2を繰り返す．

18 Sep 2021

Takamasa Iryo

Frank-Wolfe法

Feasible area
Direction of descend

16

1-D optimisation

Next value

Current
value

18 Sep 2021

Takamasa Iryo

Frank-Wolfe法

17

Next next value

18 Sep 2021

Takamasa Iryo

Frank-Wolfeの実装

18

Step 1

最短経路探索は
経路列挙で行なっている

18 Sep 2021

Takamasa Iryo

Frank-Wolfeの実装

19

Step 2

1次元探索は線形探索で行なっている

18 Sep 2021

Takamasa Iryo

MSA

20

• Method of Successive Average• Frank-Wolfeのステップ2の α を，目的関数の評価を行うことなく「適当に」求める．
– 適当というが，いちおう一定の方法はある．
– α = 1/n と設定することが多い 詳細はSheffiの教科書でも– α を十分小さい値に固定してもよい．
これは，Best response dynamicsと呼ばれるDay-to-dayのドライバーの行動調整モデルに相当する．

18 Sep 2021

Takamasa Iryo

MSA

21

• 目的関数の計算が難しい場合に利用できる．
– SUEの目的関数は経路交通量が入る• 経路なしの定式化も可能ではある（赤松他1991）– 最適化問題に変換できない問題でもいちおう使える．
• この場合は収束性が保証されないことに注意．

18 Sep 2021

Takamasa Iryo

均衡状態の一意性と安定性

22

ここは正しい（学習過程のモデルが正しければ）

ここは，本当にそれでいいのか？

18 Sep 2021

Takamasa Iryo

均衡状態の一意性と安定性

23• 均衡状態が複数ある場合は，どこに収束するかがわからない．計算の初期状態に依存する．
• 安定な均衡状態がなければ，そもそも均衡状態を将来実現する状態とすることはできない．

Beckmann et al. 1956 https://cowles.yale.edu/sites/default/files/files/pub/misc/specpub-beckmann-mcguire-winsten.pdf
18 Sep 2021

Takamasa Iryo

不安定な均衡状態

24t indicates date (continuous number)Within-day time is discretised into 100 time slots.

18 Sep 2021

Takamasa Iryo

不安定な均衡状態

18 Sep 2021

25

Takamasa Iryo

不安定な均衡状態

26t indicates date (continuous number)# travellers choosing
each time slot

18 Sep 2021

Takamasa Iryo

不安定な均衡状態

18 Sep 2021

27

Takamasa Iryo

複雑な相互作用の均衡28

It is well known that uniqueness andstability are guaranteed properties oftraffic equilibria in static user-equilibrium
traffic assignment problems, if the link
travel utilities are assumed to be strictly
monotonically decreasing with respect tothe link traffic volumes. However, these
preferable properties may not necessarily
hold in a wide range of transport
problems with complex interactions, e.g.
asymmetric interactions (including
dynamic traffic assignment), social
interactions, or with economies of scale.

18 Sep 2021

Takamasa Iryo

複雑な相互作用の動学

29• … the existence of positive interactions. In the
traditional context of transport studies, interactions
between users are considered to be mostly originatedby congestion phenomena, which creates negativeinteractions (snip). In contrast, positive interaction ischaracterised by a situation in which the increased
patronage of a transport facility increases its
attractiveness for other users.

18 Sep 2021

Takamasa Iryo

複雑な相互作用の動学

two alternatives and two user groups
two-dimensional problem

18 Sep 2021

Takamasa Iryo30

複雑な相互作用の動学

18 Sep 2021

Takamasa Iryo31

複雑な相互作用の動学

可能性のある全ケースを9種類に場合わけして性質を特定
18 Sep 2021

一意な安定均衡解があるケースは1個しかない！
Takamasa Iryo32

複雑な相互作用の動学

33

There exist two groups oftravellers, called leaders andfollowers in a city. Followers
always observe leadersbehaviour because they want
to imitate it and leaders also
have slight social interactions
between them.

Leader-follower model
端点解か，内点解か，不安定解かはDay-to-dayの動学の特性に依存する．
18 Sep 2021

Takamasa Iryo

まとめ

34

• 交通システムの将来予測や施策立案のためには，
システム内での利用者の相互作用の影響を考えなくてはならない．
• 均衡状態：相互作用の結果，最終的に安定して実現する（しそう）な状態．
– それを条件とした交通量配分：利用者均衡配分• 等価最適化問題による解法：Frank-Wolfe法• 均衡解の一意性と安定性の問題– 複雑な相互作用を考えると，これらが必ずしも成立しなくなる．
– Day-to-dayの選択行動調整過程をどうモデリングするか？

18 Sep 2021

Takamasa Iryo

付録：1次元探索

18 Sep 2021

Takamasa Iryo

35

How to solve 1D optimisation problem• The golden-section search is a well-known
method for a convex (or concave) function.
• Consider:min 𝑓𝑓 𝑥𝑥 sub. to: 0 ≤ 𝑥𝑥 ≤ 1.

f(x)f(x) is a convex functionx
18 Sep 2021

Takamasa Iryo

36

How to solve 1D optimisation problem
f(x)

0

x1a

x1b

1

x

Set two intermediate points, denoted by x1a and x1b
and calculate f(x1a) and f(x1b) .

18 Sep 2021

Takamasa Iryo

37

How to solve 1D optimisation problem

38

f(x)

0

x1a

x1b

1

x

If f(x1a) > f(x1b) , the optimal point is between x1a and 1
at all times because the function is convex.

18 Sep 2021

Takamasa IryoHow to solve 1D optimisation problem

39

f(x)

0

x1a

x2a

x2b

1

x

Do the same thing between x1 and 1 instead of 0 and 1.
Repeating this process, the section where the optimal
point must exist is getting shorter.
18 Sep 2021

Takamasa Iryo

Golden section search

40

f(x)

x

x1b x2b 1=x2a
Now find x1a and x1b so that x1b = x2a. If we can do so,
we can reduce the number of calculations of f(x).
0

18 Sep 2021

x1a

Takamasa Iryo

Golden section search!!" "!!#•= 𝑥𝑥#$#"!!#

3− 5𝑥𝑥#$ =
2• 𝑥𝑥#% = 1 − 𝑥𝑥#$

0

18 Sep 2021

41

x1a

x1b x2b

Takamasa Iryo

1xGolden section search

42

3− 5
3− 5
1+ 5
𝑎𝑎: 𝑏𝑏 =
:1 −
= 1:
2
2
2

a

0

18 Sep 2021

Golden ratio

b

x1a

1

Takamasa IryoxGolden section search: Algorithm

18 Sep 2021

Takamasa Iryo

43

44

18 Sep 2021

Takamasa Iryo

45

18 Sep 2021

Takamasa Iryo

46

18 Sep 2021

Takamasa Iryo

import mathdef func(x):return (x-0.7)**2# FW（1O Dペア専用，最短経路は経路列挙による簡便な手法）# Frank-Wolfe algorithm, with shortest path search enumerating all routes.

def goldenSearch(func,xa,xb,tolerance,isExact,isPrint):# 分割比（区間全体を1とした値）Golden ratio (entire section is set to 1)
gr = (3-math.sqrt(5))/2

import numpy as np
import pandas as pd
from scipy.optimize import minimize
import matplotlib.pyplot as plt# 区間の幅 With of the entire section# 経路交通量からO D交通量を計算する Calculating O D traffic volume from route traffic volume
def calcO DTraf(route_traf):
od_traf = 0
for route in routes:
od_traf += route_traf[route]return od_traf

wd = xb - xa
# 左側分割位置の初期値 Initial valuel of the left-hand side dividing point
xa_int = xa + wd*gr
fa = func(xa_int)# 経路交通量からリンク交通量を計算する Calculating link traffic volume from route traffic volume
def calcLinkTraf(route_traf):
link_traf = {}for link in links:
link_traf[link] = 0
for route in routes:
for link in routes[route]:
link_traf[link] += route_traf[route]return link_traf# 右側分割位置の初期値 Initial valuel of the right-hand side dividing point
xb_int = xb - wd*gr
fb = func(xb_int)
while True:# リンク交通量から利用者均衡配分用の目的関数を計算する# Calculating the objective function from link traffic volume
def calcO bjUE(link_traf):
obj = 0for link in links:
obj += 0.5*links[link][1]*link_traf[link]**2 ¥+links[link][0]*link_traf[link]
return obj# 計算結果の出力 Printing the calculation process
if isPrint:
print ("(%6.4f,%6.4f,%6.4f,%6.4f) (%8.6f,%8.6f)" % (xa,xa_int,xb_int,xb,fa,fb))# リンク交通量から各経路の旅行時間を計算する Calculating route travel time from link traffic volume
def calcRouteTT(link_traf):route_TT={}for route in routes:
route_TT[route] = 0
for route in routes:
for link in routes[route]:
route_TT[route] = route_TT[route] ¥+ links[link][1]*link_traf[link]+links[link][0]return route_TT# 2つの分割間で右下がりであれば，左側分割より左側を切り落とす
if (fa > fb):
xa = xa_int
wd = xb - xa
# 新しい左側分割の計算 Caculating new left-hand side dividing point# ここからメインルーチン Main routine# 黄金分割の特性を用いずに，あえて計算しなおしている（チェック用）# すべてのリンクとそのリンク旅行時間パラメータ（[自由流旅行時間，交通量との比例係数]）# All links and their link travel time parameters ([free-flow travel time, coeff of traffic volume])
links = {'A':[10,10],'B':[60,1],'C':[60,1],'D':[10,10],'E':[10,1]}# The following calculation does not utilse the advantahe of the golden section method,# for checking purpose only.

# すべての経路とそれに属するリンク# All routes and links included by each route.
routes={'a':['A','B'], 'b':['C','D'], 'c':['A','E','D']}

if (isExact):
xa_int = xa + wd*gr# 初期リンク交通量パターン．とりあえず経路aにO D交通量6を全部置く# 本来のFWであれば，自由流旅行時間に基づいて最短経路探索して経路を決めるべき# Initial link traffic volume pattern. Assign all OD traffic to route 'a' here.
# It must be replaced by the shortest path search with free-flow link travel times given.
O D_traf = 6
link_traf_x = calcLinkTraf({'a':OD_traf,'b':0,'c':0})

fa = func(xa_int)# 黄金分割による計算 Calculation based on the golden section method.
else:

print ("Initial link traffic volume: ",link_traf_x)
print ()

xa_int = xb_int
fa = fb# FWをとりあえず5回繰り返す FW is going to repeat for five times here.
for k in range(1,5):# 経路旅行時間を計算（本来のFWであれば，リンク旅行時間を計算してから最短経路探索をする）# Calculate route travel time, note that it must be replaced by the shortest path search# with link travel times calculated.
route_TT = calcRouteTT(link_traf_x)
minTT = 10000000
minRoute = ''for route in route_TT:
if minTT > route_TT[route]:
minTT = route_TT[route]
minRoute = route
print ("Shortest route: ",minRoute,minTT)# 新しい右側分割の計算 Caculating new right-hand side dividing point
xb_int = xb - wd*gr
fb = func(xb_int)else:
xb = xb_int
wd = xb - xa
if (isExact):
xb_int = xb - wd*gr# All-or-nothing 配分
link_traf_y = calcLinkTraf({'a':O D_traf*('a'==minRoute),'b':O D_traf*('b'==minRoute),'c':O D_traf*('c'==minRoute)})
print ("All-or-nothing assigmnent: ",link_traf_y)

fb = func(xb_int)
else:# 一次元探索（単純な方法による）# O ne-dimension search (Naive method)
minO bj = 100000000
minAl = 0
for alpha in range(0,30001):
al = alpha / 30000
link_mix = {}for link in links:
link_mix[link] = al * link_traf_y[link] + (1-al) * link_traf_x[link]obj = calcO bjUE(link_mix)
if (minO bj > obj):
minO bj = obj
minAl = alxb_int = xa_int
fb = faxa_int = xa + wd*gr
fa = func(xa_int)# 収束判定 Convergence checkif (xb_int - xa_int) < tolerance:
break

print ("1D search alpha = ",minAl," O bj = ",minO bj)# 一次元探索の最適解を用いて，リンク交通量を更新# Update link traffic volume pattern using the optimal solution of 1-D search.
for link in links:
link_traf_x[link] = minAl * link_traf_y[link] + (1-minAl) * link_traf_x[link]

return [0.5*(xb_int + xa_int),0.5*(fa+fb)]

print ("Current link traffic volume: ",link_traf_x)
print ()

print ("Answer:", goldenSearch(func,0,1,0.01,True,True))
print()
print ("Answer:", goldenSearch(func,0,1,0.01,False,True))

18 Sep 2021

Takamasa Iryo

47