四阶段法-交通分布预测方法-算法复现

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_{ij}$ 的增长与 $i$ 区产生量的增长及$j$分区吸引量的增长同时相关，而且相关的程度也相同。即 $f_{\mathrm{a}\mathrm{v}\mathrm{e}}\left(F_{pi}{F}_{aj}\right)=\frac{1}{2}\left(F_{pi}+F_{ai}\right)$

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_{ij}$的增长与i分区产生量增长率 $F_{pi}$成正比，而且还与$j$分区吸引量增长占整个区域吸引量增长的相对比率成正比。全区域在现年和规划年的吸引总量分别为${\sum }_j{A}_j/{\sum }_j{A}_j^0$。于是Detroit增长函数为 $f_D\left(F_{pi},F_{aj}\right)=F_{pi}\cdot \frac{F_{aj}}{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}$

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_{ij}$，然后平均的方法。
先从产生区考虑：
① Frator认为，$q$