元胞自动机

元胞自动机的一些基本特性：

  1. 离散空间：元胞自动机通常定义在一个规则的网格上，如一维、二维或三维的格子。

  2. 有限状态：每个元胞在任何给定时间都处于有限个状态中的一个，状态可以是二进制的（例如0和1）或具有更多可能值。

  3. 局部规则：元胞的状态更新基于局部规则，这些规则只依赖于元胞自身和其邻居的状态。

  4. 同步更新：在每个时间步，所有元胞的状态都是同时更新的。

  5. 初始条件：元胞自动机的初始状态对系统的长期行为有很大影响。

  6. 时间演化：元胞自动机的状态随时间演化，形成一系列状态。

  7. Wolfram分类：斯蒂芬·沃尔夫勒姆（Stephen Wolfram）提出了一种著名的元胞自动机分类方法，将一维元胞自动机分为四类，从简单的规则到复杂的混沌行为。

  8. 图灵完备性：某些元胞自动机被证明具有图灵完备性，这意味着它们能够模拟任何图灵机的计算。

  9. 应用：元胞自动机被用于模拟各种自然现象，如流体动力学、细胞生长、疾病传播、交通流等。

  10. 可视化：元胞自动机的演化常常通过图形方式展示，以直观地理解其动态行为。

康威的生命游戏（Conway's Game of Life）

生命游戏是英国数学家约翰·何顿·康威在1970年发明的细胞自动机。它包括一个二维矩形世界，这个世界中的每个方格居住着一个活着的或死了的细胞。一个细胞在下一个时刻生死取决于相邻八个方格中活着的或死了的细胞的数量。通常情况，游戏的规则就是：当一个方格周围有2或3个活细胞时，方格中的活细胞在下一个时刻继续存活；即使这个时刻方格中没有活细胞，在下一个时刻也会“诞生”活细胞。 规则是：

1. 对周围的 8 个近邻的元胞状态求和

2. 如果总和为 2 的话，则下一时刻的状态不改变

3. 如果总和为 3 ，则下一时刻的状态为 1，否则状态为 0

二维元胞自动机模拟：

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
%% 设置GUI按键
plotbutton=uicontrol('style','pushbutton','string','运行', 'fontsize',12, 'position',[150,400,50,20], 'callback', 'run=1;');
erasebutton=uicontrol('style','pushbutton','string','停止','fontsize',12,'position',[250,400,50,20],'callback','freeze=1;');
quitbutton=uicontrol('style','pushbutton','string','退出','fontsize',12,'position',[350,400,50,20],'callback','stop=1;close;');
number = uicontrol('style','text','string','1','fontsize',12, 'position',[20,400,50,20]);
%% 元胞自动机设置
n=200;
%初始化各元胞状态
z = zeros(n,n);
sum = z;
cells = (rand(n,n))<.6;
% 建立图像句柄
imh = image(cat(3,cells,z,z));
set(imh, 'erasemode', 'none')
% 元胞更新的行列数设置
x = 2:n-1;
y = 2:n-1;
% 主事件循环
stop= 0; run = 0;freeze = 0; 
while stop==0
    if run==1
        % 计算邻居存活的总数
        sum(x,y) = cells(x,y-1) + cells(x,y+1) + cells(x-1, y) + cells(x+1,y)...
            + cells(x-1,y-1) + cells(x-1,y+1) + cells(x+1,y-1) + cells(x+1,y+1);
        % 按照规则更新
        cells = (sum==3) | (sum==2 & cells);
        set(imh, 'cdata', cat(3,cells,z,z) )
        stepnumber = 1 + str2double(get(number,'string'));
        set(number,'string',num2str(stepnumber))
    end
    if freeze==1
        run = 0;
        freeze = 0;
    end
    drawnow
end

元胞自动机模拟案例二

规则制定： 首先定义平面中心点的数值为1，接着在每一时间步对每一点，如果周围八个点的和为偶数，则变为0，为奇数则变为 1 

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
% 颜色控制
Map = [1 1 1; 0 0 0];
colormap(Map);    %定义所有值为1的元胞将显示为白色，值为0的元胞将显示为黑色
% 设置网格大小
S = 121;
L = zeros(S);
% 把中间一个数设置为 1 作为元胞种子
M = (S+1)/2;
L(M, M) = 1;
Temp = L;
imagesc(L);
% 计算层数
Layer = (S-1)/2 + 1;

for t=2:Layer
    for x=M-t+1:M+t-1
       if x==M-t+1 || x==M+t-1

          for y=M-t+1:M+t-1
            SUM = 0;
            for m=-1:1
               for n=-1:1
                  if x+m>0 && x+m<=S && y+n>0 && y+n<=S
                     SUM = SUM + L(x+m, y+n); 
                  end
               end
            end
            SUM = SUM - L(x, y);
            Temp(x, y) = mod(SUM, 2);
          end
          
       else
            y = M-t+1;
            SUM = 0;
            for m=-1:1
               for n=-1:1
                  if x+m>0 && x+m<=S && y+n>0 && y+n<=S
                     SUM = SUM + L(x+m, y+n); 
                  end
               end
            end
            SUM = SUM - L(x, y);
            Temp(x, y) = mod(SUM, 2);
            
            y = M+t-1;
            SUM = 0;
            for m=-1:1
               for n=-1:1
                  if x+m>0 && x+m<=S && y+n>0 && y+n<=S
                     SUM = SUM + L(x+m, y+n); 
                  end
               end
            end
            SUM = SUM - L(x, y);
            Temp(x, y) = mod(SUM, 2);
       end
    end
    L = Temp;
    imagesc(L);
    % 速度控制
    pause(0.2);
end

元胞自动机的应用实例---Nagel Schreckenberg 模型模拟

模型介绍

Nagel-Schreckenberg模型是一种用于模拟交通流的简单而有效的元胞自动机模型，由德国物理学家Reinhard Nagel和Michael Schreckenberg在1992年提出。该模型主要用于描述单向道路上车辆流动的行为，能够展示出交通堵塞的形成、消散以及密度波的传播等现象。它通过规则定义车辆的加速、减速、保持速度以及换道行为，每个时间步长内，模型中的每一辆车根据以下规则更新其速度和位置：

  1. 加速规则：如果一辆车的速度v小于最大速度v_max，并且前方距离为d空白格（无车占据），则该车的速度增加1，直到达到v_max或者前方空位不足。

  2. 减速规则：考虑安全距离，如果前方只有一辆车且两车距离小于安全距离d_0（通常是车速v的一倍，即d_0=v），则必须减速直到两车距离满足安全要求，但最多减到0（即停车）。

  3. 随机刹车：即使前车距离足够，也有一定概率p降低速度，模拟驾驶员的不可预测行为或道路状况的变化，通常取p较小，如0.15。

  4. 速度限制：确保更新后的速度不会超过最大速度v_max，也不会低于0（即不会出现负速度）。

模型中的每辆车占据一格空间，每一步更新时，车辆按照其当前速度前移相应的格数，然后根据上述规则调整速度，以此循环模拟交通流随时间的动态变化。Nagel-Schreckenberg模型因其简明性和能够捕捉交通流宏观特性而广泛应用于交通工程、复杂系统科学和物理学的研究中。

代码实现

python简单模拟 

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
# 参考链接：http://t.csdnimg.cn/2QgNz
import numpy as np
from matplotlib import pyplot as plt
import time
 
start = time.time()
path = 5000  # 道路长度
n = 100  # 车辆数目
v0 = 90  # 初始速度
p = 0.3  # 车辆减速概率
Times = 6000  # 模拟的时刻数目,时刻越长所耗时间越长
 
 
np.random.seed(0)
x = np.random.rand(n) * path  # 保存每辆车在道路上的位置,随机进行初始化
x.sort()
v = np.tile([v0], n).astype(np.float64)  # 保存每辆车的速度,并且初速度相同
 
plt.figure(figsize=(8, 6), facecolor='w')
# 模拟每个时刻
for t in range(Times):
    plt.scatter(x, [t] * n, s=1, c='k', alpha=0.05)
    # 模拟每辆车
    for i in range(n):
        # 计算当前车与前车的距离（是环形车道）
        if x[(i + 1) % n] > x[i]:
            d = x[(i + 1) % n] - x[i]
        else:
            d = path - x[i] + x[(i + 1) % n]
        # 根据距离计算下一秒的速度
        if v[i] < d:
            if np.random.rand() > p:
                v[i] += 1
            else:
                v[i] -= 1
        else:
            v[i] = d - 1
    # 对速度进行限制(最高速度不得超过120)
    v = v.clip(0, 120)
    # 一秒后,车辆的位置发生变化
    x += v
    x %= path
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.xlim(0, path)
plt.ylim(0, Times)
plt.xlabel('车辆位置')
plt.ylabel('模拟时间')
plt.title(f'Nagel-Schreckenberg模型模拟(车道长度:{path},车辆数:{n},初速度:{v0},减速概率:{p})')
end = time.time()
plt.text(1800, 6400, s=f'Running time: {(end - start)} seconds', fontsize=15, color='red')
 
plt.show()

基于MATLAB编写交通流模拟程序，通过计算机仿真来探索和分析交通流量与道路密度之间的关系，这通常被称为“流量-密度曲线”（Fundamental Diagram of Traffic Flow）。

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
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
% 初始化参数
clear; 
clc;   
vmax = 10; % 最大速度
p = 0.6;       % 随机制动的概率（原始数据为0.8）
road_length = 1000; % 道路长度
simulation_steps = 1000; % 模拟步数
render_on = 0; % 是否渲染图形界面
pause_on = 0; % 是否在模拟中暂停
delay_on = 0; % 是否在模拟中延迟
delay_length = 0.02;    % 延迟时间，用于控制帧率（10 FPS）

% 初始化道路和速度状态
road = zeros(1,road_length); % 道路占用状态数组
road_next = road; % 下一状态的数组
velocities = zeros(1,road_length); % 车辆速度状态数组
velocities_next = velocities; % 下一速度状态的数组

% 采样设置
num_samples = 2000; % 采样数量
samples = zeros(2,num_samples); % 存储密度和流量的数组
density_step = 1/num_samples; % 密度步长

% 初始化历史记录数组
history = zeros(simulation_steps, road_length); % 记录道路状态历史
velocity_history = zeros(simulation_steps, road_length); % 记录速度状态历史

% 打开图形界面
figure;

% 采样循环
for g=1:num_samples
    % 生成不同密度的交通流
    road = zeros(1,road_length); % 重置道路状态
    road_next = road; % 重置下一状态
    density = g/num_samples; % 当前密度

    % 根据密度在道路上随机放置车辆
    for i=1:road_length
        if rand < density
            road(i) = 1;
        end
    end

    % 如果开启渲染，显示当前道路状态
    if render_on
        imshow(road);
        drawnow;
    end

    % 模拟循环
    for i=1:simulation_steps
        % 记录当前状态
        history(i, :) = road;
        velocity_history(i,:) = velocities;

        % 速度更新
        for j=1:road_length
            if road(j) == 1
                % 寻找前方 vmax 个单元格内没有车辆的最远距离
                distance = 0;
                bf = 0; % 标志变量，表示是否找到车辆
                for k=1:vmax
                    distance = k;
                    if j+k <= road_length
                        index = j+k;
                    else
                        index = j+k-road_length; % 处理循环道路的边界情况
                    end
                    if road(index) == 1
                        bf = 1;
                    end
                    if bf == 1, break, end
                end

                % 根据规则更新车辆速度
                if velocities(j) < vmax
                    velocities(j) = velocities(j) + 1; % 加速
                end
                if (velocities(j) > distance - 1) && bf == 1
                    velocities(j) = distance - 1; % 避免碰撞
                end
                if rand < p && velocities(j) > 0
                    velocities(j) = velocities(j) - 1; % 随机减速
                end
            end
        end

        % 移动更新
        for j=1:road_length
            if road(j) == 1
                if j+velocities(j) <= road_length
                    index = j+velocities(j);
                else
                    index = j+velocities(j) - road_length; % 处理循环道路的边界情况
                end
                % 检测碰撞
                if road_next(index) == 1
                    disp('Collision detected');
                end
                road_next(index) = 1;
                velocities_next(index) = velocities(j);
            end
        end

        % 更新道路和速度状态
        velocities = velocities_next;
        road = road_next;
        road_next = zeros(1,road_length);

        % 如果开启渲染，显示当前道路状态
        if render_on
            imshow(road);
            drawnow;
        end

        % 如果开启暂停，等待用户操作
        if pause_on
            pause;
        end

        % 如果开启延迟，暂停一段时间
        if delay_on
            pause(delay_length);
        end
    end

    % 记录当前密度下的流量和密度
    velocity_history = velocity_history.*history; % 计算流量
    samples(:,g) = [mean2(history) (sum(velocity_history(:))/sum(history(:)))*mean2(history)]; % 计算平均密度和流量
    disp('Sample step:')
    g % 显示当前采样步骤
end

% 绘制流量-密度曲线
scatter(samples(1,:), samples(2,:));
axis([0 1 0 1]); % 设置坐标轴范围
xlabel('Density'); % X轴标签
ylabel('Flow (normalized)'); % Y轴标签
title('Flow-density Curve'); % 图形标题

% 可以选择显示整个模拟过程中的道路状态或流量状态
%imshow(history);
%ts(simulation_steps,history);
  1. 初始化变量:
    • vmax = 10;：设置车辆的最大速度。
    • p = 0.6;：随机制动的概率。
    • road_length = 1000;：道路长度，单位为车道单元。
    • simulation_steps = 1000;：仿真步数。
    • render_on, pause_on, delay_on: 控制是否渲染图形、暂停和延迟的布尔变量。
    • delay_length = 0.02;：控制每帧渲染间隔时间，这里对应大约10帧每秒的刷新率。
  2. 定义数组:
    • 初始化道路状态数组road和road_next，以及速度状态数组velocities和velocities_next，用于存储车道上车辆的分布和速度。
  3. 采样循环 (for g=1:num_samples):
    • 这个循环用于在不同的车辆密度下进行多次模拟。
    • 计算当前的密度并随机生成初始交通分布。
    • 如果render_on为真，则显示当前的道路状态。
  4. 仿真循环 (for i=1:simulation_steps):
    • 速度更新: 根据前方空闲空间调整车辆速度，考虑加速、碰撞避免和随机制动。
    • 移动更新: 根据当前速度移动车辆，并处理道路两端的环绕逻辑。
    • 检查碰撞，更新road_next和velocities_next，然后将它们应用于当前状态。
    • 再次检查是否需要渲染、暂停或延迟。
  5. 后处理:
    • 在每次采样结束后，计算并记录平均密度和流量（流率），存储于samples矩阵中。
    • 最后，使用散点图展示所有样本点，横坐标为密度，纵坐标为归一化后的流量，得到流量-密度曲线。
  6. 可视化:
    • 使用MATLAB的图形功能绘制流量-密度曲线，帮助直观理解不同密度下的平均流量变化。

注意点：

这里的"密度"（density）是指在模拟的道路长度单位内单位时间内存在的车辆数，或者说是单位长度车道上的车辆数。具体到代码实现中，它是通过变量density = g/num_samples;来定义的，其中g是当前的采样序号，num_samples是总的采样数量，因此density实质上代表了一个从0到1范围内的比例值，用来模拟从完全没有车辆（密度为0）到道路完全被车辆占据（密度接近1）的各种情况。这个密度值决定了初始时在模拟道路上放置车辆的随机概率，即if rand < density这一条件控制下的车辆分布。

ts 的函数，用于绘制交通模拟的时间-空间图。时间-空间图是一种展示随时间变化的空间分布情况的图表，常用于交通流模拟中，以观察车辆在道路上随时间的移动情况。

1
2
3
4
5
6
7
8
9
10
function ts(simulation_steps,history)
figTS=figure;
axes1 = axes('Parent',figTS,'FontSize',12,'FontName','Arial');
for i=1:1:simulation_steps
    [~,col]=find(history(i,:)==1);
    plot(i*ones(1,length(col)),col,'k.');
    hold on
end
xlabel('Time','FontSize',12,'FontName','Arial')
ylabel('Space','FontSize',12,'FontName','Arial')

  1. function ts(simulation_steps,history)：定义了一个名为 ts 的函数，它接受两个参数：simulation_steps（模拟步数）和 history（记录每一步道路状态的矩阵）。

  2. figTS=figure;：创建一个新的图形窗口。

  3. axes1 = axes('Parent',figTS,'FontSize',12,'FontName','Arial');：在图形窗口中创建一个坐标轴，设置字体大小为12，字体类型为 Arial。

  4. for i=1:1:simulation_steps：开始一个循环，从1到 simulation_steps（模拟步数）。

  5. [~,col]=find(history(i,:)==1);：在每次模拟步 i 的历史记录中，找到所有值为1（表示有车辆）的位置，并将其存储在 col 变量中。find 函数返回的首个输出（在这里被忽略，用 ~ 表示）是行索引，第二个输出是列索引。

  6. plot(i*ones(1,length(col)),col,'k.');：在坐标轴上绘制点。i*ones(1,length(col)) 创建一个与 col 长度相同的数组，所有元素都是当前的模拟步数 i，代表时间。col 代表空间位置。'k.' 指定了点的颜色和样式，这里是黑色圆点。

  7. hold on：保持当前的坐标轴和图形，以便在同一个图上绘制更多的数据。

  8. xlabel('Time','FontSize',12,'FontName','Arial')：设置 X 轴的标签为 "Time"，并设置字体大小和类型。

  9. ylabel('Space','FontSize',12,'FontName','Arial')：设置 Y 轴的标签为 "Space"，并设置字体大小和类型。

上述函数的功能是绘制出每个时间步长上车辆的位置，从而形成时间-空间图。图中的每个点代表在特定时间步长下道路上的一个车辆位置。通过观察这个图，可以直观地看到车辆随时间在道路上的分布和移动情况。

数模
元胞自动机
http://jrhu0048.github.io/2024/05/01/yuan-bao-zi-dong-ji/
作者
JR.HU
发布于
2024年5月1日
更新于
2024年10月15日
许可协议