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

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)=FpiQ/Q0Faj=Pi0PijAj/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

  • 11
    点赞
  • 85
    收藏
    觉得还不错? 一键收藏
  • 13
    评论
评论 13
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值