Traffic Distribution
这里整理一下四阶段法中交通分布预测阶段中的重要预测方法。具体包括:平均系数法、底特律法、Frater法、Funess法、单约束重力模型法、双约束重力模型法,并对比一下各种方法的求解收敛速度。
import numpy as np
import matplotlib.pyplot as plt
import math
首先输入所需数据,包括现状OD表、规划年发生吸引量以及小区间阻抗矩阵
#定义初始OD矩阵
init_OD=np.array([[200, 100, 100],[150, 250, 200],[100, 150, 150]],dtype=float)
init_O=np.array([400, 600, 400],dtype=float)
init_D=np.array([450, 500, 450],dtype=float)
#定义规划年数据
later_O=np.array([1000,1000,1250])
later_D=np.array([1250,900,1100])
#阻抗矩阵
cost_mat=np.array([[14, 32, 40],[32, 16, 22],[40, 22, 16]])
stop = 0.03
再通过几种方法求解远景年OD矩阵 (各方法具体步骤及代码封装在最后)
并记录各种方法求解的收敛速度,进行比较。各个算法收敛效率比较如下图。
x_ave,y_ave=ave_para(init_OD, init_O, init_D, later_O, later_D, stop)
x_de,y_de=detroit(init_OD, init_O, init_D, later_O, later_D, stop)
x_fra,y_fra=frator(init_OD, init_O, init_D, later_O, later_D, stop)
x_fu,y_fu=funes(init_OD, init_O, init_D, later_O, later_D, stop)
x_gra_single,y_gra_single=gravity_model_single(init_OD, init_O, init_D, later_O, later_D, stop, cost_mat)
x_gra_double,y_gra_double=gravity_model_double(init_OD, init_O, init_D, later_O, later_D, stop, cost_mat)
plt.figure(figsize=(8, 6))
plt.title('收敛曲线(系数法)') # 折线图标题
plt.rcParams['font.sans-serif'] = ['SimHei'] # 显示汉字
plt.xlabel('迭代次数') # x轴标题
plt.ylabel('指标') # y轴标题
plt.xlim([0, 10])
plt.ylim([0, 2.5])
plt.plot(y_ave, label='平均系数法') # 绘制折线图,添加数据点,设置点的大小
plt.plot(y_de,label='底特律法')
plt.plot(y_fra,label='福莱特法')
plt.plot(y_fu,label='福尼斯法')
plt.plot(y_gra_single,label='单约束重力')
plt.plot(y_gra_double,label='双约束重力')
plt.axhline(0.03, color='red',linestyle='--')
plt.legend()
plt.show()
平均系数法
该方法思路为: q i j q_{ij} qij 的增长与 i i i 区产生量的增长及 j j j分区吸引量的增长同时相关,而且相关的程度也相同。即 f a v e ( F p i F a j ) = 1 2 ( F p i + F a i ) {f_{ {\rm{ave}}}}\left( { {F_{pi}}{F_{aj}}} \right) = \frac{1}{2}\left( { {F_{pi}} + {F_{ai}}} \right) fave(FpiFaj)=21(Fpi+Fai)
def ave_para(init_OD, init_O, init_D, later_O, later_D, stop):
f_O=later_O/init_O
f_D=later_D/init_D
record=[]
a=np.vstack((f_O,f_D))-1
a=np.maximum(a, -a)
record.append(np.max(a))
while np.max(a)>stop:
for i in range(init_OD.shape[0]):
init_OD[i,:]=init_OD[i,:]*(f_O[i]+f_D)/2
init_O=init_OD.sum(axis=1)
init_D=init_OD.sum(axis=0)
f_O=later_O/init_O
f_D=later_D/init_D
a=np.vstack((f_O,f_D))-1
a=np.maximum(a, -a)
record.append(np.max(a))
return init_OD,record
Detroit法
此方法思路为, q i j q_{ij} qij的增长与i分区产生量增长率 F p i F_{pi} Fpi成正比,而且还与 j j j分区吸引量增长占整个区域吸引量增长的相对比率成正比。全区域在现年和规划年的吸引总量分别为 ∑ j A j / ∑ j A j 0 \sum_{j} A_{j} / \sum_{j} A_{j}^{0} ∑jAj/∑jAj0。于是Detroit增长函数为 f D ( F p i , F a j ) = F p i ⋅ F a j Q / Q 0 = P i P i 0 ⋅ A j / A j 0 ∑ j A j / ∑ j A j 0 f_{D}\left(F_{p i}, F_{a j}\right)=F_{p i} \cdot \frac{F_{a j}}{Q / Q^{0}}=\frac{P_{i}}{P_{i}^{0}} \cdot \frac{A_{j} / A_{j}^{0}}{\sum_{j} A_{j} / \sum_{j} A_{j}^{0}} fD(Fpi,Faj)=Fpi⋅Q/Q0Faj=Pi0Pi⋅∑jAj/∑jAj0Aj/Aj0
def detroit(init_OD, init_O, init_D, later_O, later_D, stop):
f_O=later_O/init_O
f_D=later_D/init_D
record=[]
a=np.vstack((f_O,f_D))-1
a=np.maximum(a, -a)
record.append(np.max(a))
while np.max(a)>stop:
for i in range(init_OD.shape[0]):
init_OD[i,:]=init_OD[i,:]*f_O[i]*f_D/(np.sum(later_D)/np.sum(init_D))
init_O=init_OD.sum(axis=1)
init_D=init_OD.sum(axis=0)
f_O=later_O/init_O
f_D=later_D/init_D
a=np.vstack((f_O,f_D))-1
a=np.maximum(a, -a)
record.append(np.max(a))
return init_OD,record
Frator法
1954年Frator提出了分别从产生区和吸引区两个角度分析计算 q i j q_{ij} qij,然后平均的方法。
先从产生区考虑:
① Frator认为, q i j q_{ij} q