前言
对于我的第一篇文章：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 结果如下：
除了底特律法外，其他结果和书中非常接近，应该是四舍五入的问题。底特律法的结果我猜测是书中的结果有误，因为我用计算器跟着书上的例子算了一遍，最后一步迭代的时候，书上算错了。
红框部分，用计算器按一下就知道结果不对。计算器的结果和我代码跑出来的结果对得上。
我原来的文章是用的王炜老师的第二版《交通规划》中的算例。可以发现这俩算例数据和步骤一模一样，底特律法结果同样是有问题的，可能是数字太多疏忽了吧。
![[NOIP2011 普及组] 数字反转](https://img-blog.csdnimg.cn/direct/29abc83666914b5f8579b24fc746fbc9.png)
