1 为什么用元胞自动机做交通流

  交通拥堵有个反直觉的地方：它并不是由某个“堵点”造成的，而是大量车辆按各自简单规则行驶时，自发涌现出来的集体现象。前车减速、后车跟着减速、减速逐级放大——这种“走走停停”的波（stop-and-go wave）在没有路口、没有事故的笔直路段上也会出现。

  要研究这种涌现，宏观连续模型（把车流当作可压缩流体）就不太够用了：它擅长描述整体的流量—密度—速度关系，却难以刻画单个驾驶员的行为如何逐级放大成拥堵。元胞自动机（Cellular Automaton）恰好补上这一环——道路被离散成一维元胞格，每格有车或没车，每辆车按固定规则更新位置与速度，规则简单到可以逐条写清楚，但迭代足够多步之后，拥堵会自己长出来。

  本文实现的是这一思路最经典的模型：NaSch 模型（Nagel 与 Schreckenberg，1992）。它的四条更新规则只有加速、减速、随机减速、前进——随机减速那一项是关键，正是它让车流从自由流相变为拥堵流。

  本文用 MATLAB 完整实现该模型，给出源代码与三组不同参数下的运行结果，用于对比不同参数对拥堵形成的影响。

2 交通流理论基础

  描述一段路上的车流，最常用的三个量是：流量 QQ（单位时间通过断面的车辆数）、速度 VV、密度 KK（单位长度上的车辆数）。三者满足一个简单关系：

Q=K×VQ=K\times V

  这个式子看着平常，但它约束了车流的基本行为：密度很小时车可以跑得快、流量随密度上升；密度增大到一定程度后车速被迫下降，流量反而减少——最终在堵死时速度与流量都归零。把 QQ 对 KK 作图得到的曲线即基本图，是判别交通运行状态的基本工具。

  围绕 QQ–KK–VV 的关系，交通流模型分两大类。宏观连续模型把车流当作可压缩流体，用密度、速度的连续方程描述整体演化，擅长给出 QQ–KK 关系的解析形式（如格林希尔茨线性模型、格林伯格对数模型），但无法刻画单个驾驶员的行为。微观模型（跟驰模型、元胞自动机）则逐个描述车辆，规则简单但能自然产生集体现象——拥堵的形成与消散正是这类模型最擅长解释的问题。

  本文采用的是元胞自动机一类中最经典的 NaSch 模型。

3 元胞自动机模型（NaSch模型）构建与优化

3.1 元胞自动机模型（NaSch模型）原理[4]

  NaSch 模型由 Nagel 与 Schreckenberg 于 1992 年提出。道路被离散成一条一维元胞链，每个元胞只有“有车”和“空”两种状态，一个元胞至多容纳一辆车；每辆车带有 0 到 vmaxv_{max} 之间的整数速度。模型每次迭代按顺序执行以下四个操作（四个步骤对全体车辆同时进行）：

  1.广泛性加速

  对于所有未达到最大速度vmaxv_{max}（道路的速度限制）的车辆，加速一个单位，即vn→min{vn+1,vmax}v_n\rightarrow min\left\{v_n+1,v_{max}\right\}。

  2.安全性减速

  在进行第a步之后，如果它的速度大于与前车间距的空元胞dnd_n，则它将速度降低到vnv_n，vn→min{dn,vn}v_n\rightarrow min\left\{d_n,v_n\right\}。其中dn=xn+1−xn−lnd_n=x_{n+1}-x_n-l_n，xnx_n代表车辆的位置，lnl_n代表车辆的长度。这里为了简化模型，假设ln=1\pmb{l_n=1}，假设各车最大减速加速度为正无穷大，即可以不耗时地减速至任意速度，避免模型中出现撞车的情形。

  3.随机减速

  对于速度大于00的汽车，速度以概率pp减少一个单位，即vn→发生的概率为pmax{vn−1,0}v_n\xrightarrow{发生的概率为p}max\left\{v_n-1,0\right\}。

  4.车辆更新

  在步骤1-3之后，基于其当前速度vnv_n，为每辆汽车分配新的位置xn→xn+vnx_n\rightarrow x_n+v_n。对于该模型在时刻t→t+1t\rightarrow t+1的演化过程如下图1所示。经过一次1-4步骤称为一次迭代。

图1 四个基本操作控制示例图

图1 四个基本操作控制示例图

3.2 参数定义

  模型共 5 个输入参数：

参数 含义 作用
NN 车辆数 与道路长度共同决定初始密度 KK
RR 初始状态最右端车辆位置 即道路长度（元胞总数）
vmaxv_{max} 最大速度 道路限速；决定自由流速度
pp 随机减速概率 模型的关键参数，刻画驾驶行为的随机性
TT 模拟迭代次数 演化总步数

  其中 pp 是 NaSch 模型的精髓所在。若取 p=0p=0（完全确定性），车流最终会演化成一个稳定的周期状态，永远不会自发形成拥堵；只有当 p>0p>0 时，随机减速造成的微小扰动才会逐级放大，最终触发自由流到拥堵流的相变。

3.3 建模思路

3.3.1 整体思路

  建立一个长度NN为的行向量xltxlt

xlt=[x1,x2,...,xN]xlt=[x_1,x_2,...,x_N]

来表示这辆车的位置。其中，的元素在满足

1≤xi≤R,xi∈N∗1\le x_i\le R,x_i\in N^*

的条件下，初始随机化。

  再建立一个长度为NN的行向量vltvlt

vlt=[v1,v2,...,vN]vlt=[v_1,v_2,...,v_N]

来表示这辆车的速度。其中，vltvlt的元素viv_i在满足

1≤vi≤vmax,vi∈N∗1\le v_i\le v_{max},v_i\in N^*

  在xltxlt和vltvlt的基础上进行一次迭代，更新xltxlt和vltvlt内各元素的值。

  重复迭代TT次，输出最终的xltxlt和vltvlt。

3.3.2 对于一次迭代

  对于一次迭代，可视为一个函数。其输入值是迭代前的（初始的或上一次迭代后的）和，其输出值是迭代后的xltxlt和vltvlt。

  对于一次迭代，有四次操作。每次操作，遍历nn从11到NN：

vn→min{vn+1,vmax}v_n\rightarrow min\left\{v_n+1,v_{max}\right\}

vn→min{xn+1−xn−ln,vn}v_n\rightarrow min\left\{x_{n+1}-x_n-l_n,v_n\right\}

vn→发生的概率为pmax{vn−1,0}v_n\xrightarrow{发生的概率为p}max\left\{v_n-1,0\right\}

xn→xn+vnx_n\rightarrow x_n+v_n

  注意到，前三次操作，都只牵扯到vltvlt的修改与更新，并没有修改xltxlt。并且vnv_n的值与vn+1v_{n+1}无关。故可将前三次操作分别进行三次遍历简化为一次遍历下的三次操作。即简化为：

  遍历nn从11到NN：

vn→min{vn+1,vmax}v_n\rightarrow min\left\{v_n+1,v_{max}\right\}

vn→min{xn+1−xn−ln,vn}v_n\rightarrow min\left\{x_{n+1}-x_n-l_n,v_n\right\}

vn→发生的概率为pmax{vn−1,0}v_n\xrightarrow{发生的概率为p}max\left\{v_n-1,0\right\}

  遍历nn从11到NN：

xn→xn+vnx_n\rightarrow x_n+v_n

  此次迭代结束，输出迭代后的xltxlt和vltvlt。

3.3.3 MATLAB源代码

  本文中使用的MATLAB版本是R2024b。

  根据3.3 建模思路，可以写出以下的MATLAB源代码。源代码已上传Github，可以点击链接，下载源代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
clc;
clear;
%输入参数
N = input('车辆数N:');
R = input('初始状态最右端车辆位置R：');
vmax = input('最大速度vmax:');
p = input('随机减速概率p:');
T = input('模拟迭代次数T:');

disp('初始状态：')
%随机生成1xN大小的递增的行向量xlt,要求所有元素大于等于1，小于等于R
xlt = sort(randperm(R, N))
%随机生成1xN大小的行向量vlt
vlt = randi([0, vmax], 1, N)

%进行T次迭代
for i = 1:T
    [xlt, vlt] = lteration(xlt, vlt, N, vmax, p);
end

%输出结果
disp('最终状态：')
xlt
vlt

%一次迭代函数定义
function [x, v] = lteration(xlt0, vlt0, n, vmax, p)
    x = xlt0;
    v = vlt0;
    x(n + 1) = inf;
    for i = 1:n
        v(i) = min(v(i) + 1, vmax);
        v(i) = min(x(i + 1) - x(i) - 1, v(i));
        % 以概率p执行
        if rand < p
            v(i) = max(v(i) - 1, 0);
        end
    end
    for i = 1:n
        x(i) = x(i) + v(i);
    end
    x(n + 1) = [];
end

3.3.4 样例输出

3.3.4.1 样例1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
车辆数N:4
初始状态最右端车辆位置R：7
最大速度vmax:2
随机减速概率p:0.3
模拟迭代次数T:50
初始状态：
xlt =
3 5 6 7
vlt =
1 2 2 0
最终状态：
xlt =
68 71 75 82
vlt =
2 1 2 2
3.3.4.2 样例2
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
车辆数N:4
初始状态最右端车辆位置R：7
最大速度vmax:2
随机减速概率p:0.3
模拟迭代次数T:50
初始状态：
xlt =
2 3 4 6
vlt =
1 1 0 0
最终状态：
xlt =
73 79 83 89
vlt =
2 2 2 1
3.3.4.3 样例3
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
车辆数N:100
初始状态最右端车辆位置R：1000
最大速度vmax:5
随机减速概率p:0.1
模拟迭代次数T:30
初始状态：
xlt =
列 1 至 10
1 9 18 24 47 95 108 125 127 160
列 11 至 20
166 175 176 179 181 208 220 221 223 228
列 21 至 30
239 247 261 264 265 270 284 301 313 323
列 31 至 40
346 356 388 389 391 411 454 461 467 495
列 41 至 50
503 504 508 524 532 538 547 552 556 582
列 51 至 60
583 622 645 646 654 657 679 717 720 732
列 61 至 70
737 743 759 760 761 762 777 789 790 797
列 71 至 80
801 802 809 815 819 832 833 843 853 856
列 81 至 90
858 859 881 890 916 917 918 921 943 948
列 91 至 100
955 956 959 964 967 968 982 989 995 998
vlt =
列 1 至 10
0 2 1 1 4 1 1 0 1 4
列 11 至 20
2 3 5 2 3 5 2 3 1 5
列 21 至 30
1 1 5 2 1 0 1 1 5 4
列 31 至 40
4 5 2 4 2 4 2 5 0 5
列 41 至 50
1 0 0 2 5 1 4 1 4 1
列 51 至 60
0 5 4 5 5 3 5 2 5 0
列 61 至 70
0 4 2 2 3 0 2 0 4 3
列 71 至 80
0 2 2 3 1 0 2 4 3 5
列 81 至 90
1 1 1 4 4 5 4 3 2 3
列 91 至 100
5 5 2 0 2 3 0 5 5 2
最终状态：
xlt =
列 1 至 5
140 149 157 163 192
列 6 至 10
237 250 261 267 282
列 11 至 15
297 303 311 317 325
列 16 至 20
333 341 350 358 368
列 21 至 25
374 380 389 396 402
列 26 至 30
408 426 440 461 471
列 31 至 35
494 500 517 523 529
列 36 至 40
556 582 592 599 617
列 41 至 45
626 634 641 658 665
列 46 至 50
675 682 691 704 712
列 51 至 55
721 770 778 784 792
列 56 至 60
799 819 827 833 841
列 61 至 65
847 856 866 872 882
列 66 至 70
889 896 905 912 917
列 71 至 75
923 928 934 940 945
列 76 至 80
951 959 965 973 981
列 81 至 85
991 997 1011 1017 1024
列 86 至 90
1037 1043 1048 1057 1063
列 91 至 95
1070 1077 1085 1092 1097
列 96 至 100
1111 1117 1127 1133 1141
vlt =
列 1 至 10
5 4 5 5 5 5 5 5 5 5
列 11 至 20
5 5 5 4 5 5 5 5 5 5
列 21 至 30
5 4 5 5 5 5 5 5 5 5
列 31 至 40
5 5 5 5 5 5 5 5 5 5
列 41 至 50
5 4 5 5 5 5 5 5 5 5
列 51 至 60
5 5 5 5 5 5 5 5 5 5
列 61 至 70
5 5 5 5 5 5 5 4 5 4
列 71 至 80
5 4 5 5 4 5 5 5 5 5
列 81 至 90
5 5 5 5 5 5 5 4 5 5
列 91 至 100
5 5 5 5 4 5 5 5 5 5

致谢

  感谢Matlab老师这一个学期的耐心教导！

  感谢父母与朋友在生活中无微不至的照顾！

  在我这么多朋友中，我尤其感谢远在河北某本科上学的异父异母的亲兄弟——Hana，感谢他在编程方面对我的帮助。感谢他凌晨一点半还热心地教我怎样在Github上上传文件，我的博客网站得以建成，他功不可没。

力学Matlab
许可协议

本文由Primerusse原创，采用CC BY-NC-SA 4.0 许可协议，转载请注明出处。

分享文章