用matlab实现交通分布预测方法——增长系数法

前言

对于我的第一篇文章：Matlab实现交通分布预测方法 —— 增长系数法 | 平均增长率法、底特律法、福莱特法，有不少同学私信我询问关于如何在 matlab 中调用函数、拆分代码以及需要大量调用的问题。于是我便想着把内容做一些优化，放在这篇文章里。

原来那篇文章开头有我写第一次文章时的感慨，每次打开看到感觉和下面相对理性的内容不符。。

增长系数法算法步骤

首先准备条件要有：

• 现状完整 OD 表
• 将来 OD 表中各小区发生、吸引交通量及生成总量

算法流程图如下：

matlab 代码

原来那篇文章，我是设置的交互式代码，即需要在命令行手动输入选择的方法。当需要大量使用时，有些不便。于是我把每个方法都封装成一个函数，需要哪个调用即可。

另外有些同学没有 matlab 基础，不知道如何调用函数，可以看我原来那篇文章里的使用说明部分。如果还是有困难，有需求的人多的话，我会出一期视频教程。

平均增长系数法

先上代码：

function [ res,o,d,m,T,Fo,Fd] = func_Average( od_now,U,V )
%FUNC_AVERAGE 平均增长系数法
%  res 为最终预测的od矩阵
%  o,d分别为未来各小区的总吸引量、发生量
%  m为迭代次数，T为未来出行生成总量
%  Fo,Fd分别为最终发生、吸引增长系数值
%  od_now 为现状OD矩阵（n x n） 
%  U为未来各小区发生量，V为未来各小区吸引量，均为行向量。
n=size(od_now,1);     % 小区个数
O = sum(od_now,2)';    % 现状小区发生量
D = sum(od_now,1);    % 现状小区吸引量
% 计算发生增长系数和吸引增长系数
F_o = U ./ O;   F_d = V ./ D;
alpha = 0.03;    % 误差常数
m = 0;     % 记录迭代次数
while true     % 开始循环
    m = m+1;
    % 计算系数
    f = zeros(n);
    for i = 1:n
      for j = 1:n
         f(i,j)=0.5*(F_o(i)+F_d(j)); 
      end
    end
    % 计算新od表
    od_now = f .* od_now;
    % 更新现状O和D
    O = sum(od_now,2)';    
    D = sum(od_now,1);
    % 计算新的发生和吸引系数
    F_o = U ./ O;   F_d = V ./ D;
    % 判断收敛
    if abs([F_o,F_d]-1) <= alpha
       res = od_now;
       o = sum(res,2)';
       d = sum(res,1);
       T = sum(o);
       Fo = F_o;
       Fd = F_d;
       break
    elseif m >= 100
        disp('请注意，迭代已超过100次，下面输出第100次迭代结果。');
        res = od_now;
        o = sum(res,2)';
        d = sum(res,1);
        T = sum(o);
        Fo = F_o;
        Fd = F_d;
        break
    end
end
end

我实现增长，用的是两个矩阵点乘。唯一相较困难的是增长系数放在循环结构里的位置。下面其他方法基本上如法炮制。

底特律法

function [ res,o,d,m,T,Fo,Fd] = func_Detroit( od_now,U,V )
%FUNC_Detroit 平均增长系数法
%  res 为最终预测的od矩阵
%  o,d分别为未来各小区的总吸引量、发生量
%  m为迭代次数，T为未来出行生成总量
%  Fo,Fd分别为最终发生、吸引增长系数值
%  od_now 为现状OD矩阵（n x n） 
%  U为未来各小区发生量，V为未来各小区吸引量，均为行向量。
n=size(od_now,1);     % 小区个数
O = sum(od_now,2)';    % 现状小区发生量
D = sum(od_now,1);    % 现状小区吸引量
% 计算发生增长系数和吸引增长系数
F_o = U ./ O;   F_d = V ./ D;
alpha = 0.03;    % 误差常数
m = 0;     % 记录迭代次数
while true     % 开始循环
    m = m+1;
    % 计算系数
    G = sum(O)/sum(U);  % 生成量增长率的倒数
    f = zeros(n);
    for i = 1:n
      for j = 1:n
         f(i,j)=G*F_o(i)*F_d(j); 
      end
    end
    % 计算新od表
    od_now = f .* od_now;
    % 更新现状O和D
    O = sum(od_now,2)';    
    D = sum(od_now,1);
    % 计算新的发生和吸引系数
    F_o = U ./ O;   F_d = V ./ D;
    % 判断收敛
    if abs([F_o,F_d]-1) <= alpha
       res = od_now;
       o = sum(res,2)';
       d = sum(res,1);
       T = sum(o);
       Fo = F_o;
       Fd = F_d;
       break
    elseif m >= 100
        disp('请注意，迭代已超过100次，下面输出第100次迭代结果。');
        res = od_now;
        o = sum(res,2)';
        d = sum(res,1);
        T = sum(o);
        Fo = F_o;
        Fd = F_d;
        break
    end
end
end

福莱特法

function [ res,o,d,m,T,Fo,Fd] = func_Fratar( od_now,U,V )
%FUNC_Detroit 平均增长系数法
%  res 为最终预测的od矩阵
%  o,d分别为未来各小区的总吸引量、发生量
%  m为迭代次数，T为未来出行生成总量
%  Fo,Fd分别为最终发生、吸引增长系数值
%  od_now 为现状OD矩阵（n x n） 
%  U为未来各小区发生量，V为未来各小区吸引量，均为行向量。
n=size(od_now,1);     % 小区个数
O = sum(od_now,2)';    % 现状小区发生量
D = sum(od_now,1);    % 现状小区吸引量
% 计算发生增长系数和吸引增长系数
F_o = U ./ O;   F_d = V ./ D;
alpha = 0.03;    % 误差常数
m = 0;     % 记录迭代次数
while true     % 开始循环
    m = m+1;
    % 计算系数
    f = zeros(n);
    for i = 1:n
      for j = 1:n
         Li = O(j)/sum((od_now(:,j)' .* F_d));
         Lj = D(i)/sum((od_now(i,:) .* F_o));
         f(i,j)=0.5*(Li+Lj)*F_o(i)*F_d(j); 
      end
    end
    % 计算新od表
    od_now = f .* od_now;
    % 更新现状O和D
    O = sum(od_now,2)';    
    D = sum(od_now,1);
    % 计算新的发生和吸引系数
    F_o = U ./ O;   F_d = V ./ D;
    % 判断收敛
    if abs([F_o,F_d]-1) <= alpha
       res = od_now;
       o = sum(res,2)';
       d = sum(res,1);
       T = sum(o);
       Fo = F_o;
       Fd = F_d;
       break
    elseif m >= 100
        disp('请注意，迭代已超过100次，下面输出第100次迭代结果。');
        res = od_now;
        o = sum(res,2)';
        d = sum(res,1);
        T = sum(o);
        Fo = F_o;
        Fd = F_d;
        break
    end
end
end

测试

用邵春福老师的第二版《交通规划原理》中的例子来测试一下。

新建一个脚本文件，放进以下代码：

clear;clc
od_now = [17 7 4;
     7 38 6;
     4 5 17;];   % 现状od
U=[38.6 91.9 36]; % 未来各小区发生量
V=[39.3 90.3 36.9]; % 未来各小区吸引量
% [ res,o,d,m,T] = func_Average(od_now, U, V);
% [ res,o,d,m,T] = func_Detroit(od_now, U, V);
[ res,o,d,m,T] = func_Fratar(od_now, U, V);

倒数三行代码是用来调用函数的，想使用哪个方法，就去掉其前面的百分号注释，像现在就是在使用福莱特法。各种方法得到的 res 结果如下：

除了底特律法外，其他结果和书中非常接近，应该是四舍五入的问题。底特律法的结果我猜测是书中的结果有误，因为我用计算器跟着书上的例子算了一遍，最后一步迭代的时候，书上算错了。

红框部分，用计算器按一下就知道结果不对。计算器的结果和我代码跑出来的结果对得上。

我原来的文章是用的王炜老师的第二版《交通规划》中的算例。可以发现这俩算例数据和步骤一模一样，底特律法结果同样是有问题的，可能是数字太多疏忽了吧。