元胞自动机 - ZYL-Harry/Mathematical_Modeling_Algorithms GitHub Wiki
一种时间、空间、状态都离散，空间相互作用和时间因果关系为局部的网格动力学模型，具有模拟复杂系统时空演化过程的能力。
元胞自动机的模拟仿真思路是：
- 确定相应规则
- 编程时确定规则的顺序
- 确定对象
- 对于有需要的复杂的对象用元胞来存放每个对象元素的相关数据。
- 燃烧的树将变成空地
- 如果未燃烧的树的临近周边有燃烧的树，则其也将变成燃烧的树
-
空地有概率p变为未燃烧的树
- 在临近周边没有燃烧的树的情况下，未燃烧的树有概率f被闪电击中而变成燃烧的树
figure; %利用当前属性创建窗格
axes; %建立坐标轴
rand('state',0);
set(gcf,'DoubleBuffer','on');-
rand('state',0)的作用是：使得每次产生的随机数是一致的，从而使得在包含随机数的程序每次执行时最初的样子是保持一致的 -
set(gcf,'DoubleBuffer','on')的作用是：防止在不断循环画动画的时候产生闪烁的现象
p=0.3; % 概率p，用来表示树生长的速度
f=6e-5; % 概率f,用来表示闪电击中的概率S=round(rand(300)*2); %将 S 的每个元素四舍五入为最近的整数（0、1、2）
Sk=zeros(302); %初始化Sk矩阵
Sk(2:301,2:301)=S; %加边开始的森林初值，森林的边界为一列或者一行零元素-
S=round(rand(300)*2)的作用是：随机生成一个300*300的矩阵，其中的值通过四舍五入得到——0（空地）、1（未燃烧的树）、2（燃烧的树） -
Sk(2:301,2:301)=S的作用是将生成的随机数矩阵S嵌入到全零Sk矩阵中，构成有边界的森林初始模型，边界用来处理森林原边界树的状态问题
C=zeros(302,302,3);%构造一个302*302*3的向量
R=zeros(300);%初始化R矩阵（红色通道）
G=zeros(300);%初始化G矩阵（绿色通道）
R(S==2)=1;%构成的向量的一个位置上对应的数为2，则给R赋值1
G(S==1)=1;%构成的向量的一个位置上对应的数为1，则给R赋值1
C(2:301,2:301,1)=R;%将C的第一层赋值为R矩阵
C(2:301,2:301,2)=G;%将C的第二层赋值为G矩阵
Ci=imshow(C);%在图窗中显示森林图像 I
ti=0;%初始化时间
tp=title(['T = ',num2str(ti)]);%在标题处显示时间- 构建三维向量C，这样可以使用颜色来更清楚地展示森林模型
- 需要注意处理几个条件的顺序
ti=ti+1; %时间的递增
St=Sk; %St表示下一时刻（即将）的森林情况，Sk矩阵代表着这一时刻的森林情况
St(Sk==2)=0; %正在燃烧的树变成空格位，经过一个时间间隔后-
while 1用来使得程序被不断循环执行 -
ti=ti+1表示每次循环时时间的递增 -
St=Sk：表示t时刻的森林初始情况 -
St(Sk==2)=0：燃烧的树变成空地
Sf=Sk; %Sf初始化为Sk
Sf(Sf<1.5)=0; %将森林中的空白点和有树的点去掉，只留下着火点
Sf=Sf/2; %着火点变为1，此处Sf只有着火和空格两种- 这段程序将生成一个用于判断哪些格子处于着火状态（燃烧的树），因为要进行后序（判断周边是否有着火点），所以需要Sf为302*302矩阵
Su=zeros(302); %初始化Su矩阵
Su(2:301,2:301)=Sf(1:300,1:300)+Sf(1:300,2:301)+Sf(1:300,3:302) +...
Sf(2:301,1:300)+Sf(2:301,3:302)+Sf(3:302,1:300) + ...
Sf(3:302,2:301)+Sf(3:302,3:302);
St((Su>0.5)&(Sk==1))=2; %如果绿树格位的最近邻居中有一个树在燃烧，则它变成正在燃烧的树；- 对于矩阵中的一个点来说，周围八个点只要有一个点大于0，且其本身为1那么这个点就要变成着火状态
- 所以，具体过程是：这个命令是将矩阵中的一个点的周围八个点的数相加，这个数值就是周围八个点着火点的个数，然后执行
St((Su>0.5)&(Sk==1))=2判断如果在这个位置和不为零且前一时刻这个位置有未燃烧的树，则是该未燃烧的树变成燃烧的树
Se=Sk(2:301,2:301); %将Se初始化为前一时刻的森林
Se(Se<0.5)=4; %空白地方赋值为4
Se(Se<3)=0; %有树和着火赋值为0
Se(Se>3)=1; %空白地方赋值为1
%三步操作一起将空白处变为1，不是空白处变为0
St(2:301,2:301)=St(2:301,2:301)+Se.*(rand(300)>p);- Se是用来判断哪些格子处于空地状态，因为只是判断空地，所以Se设置为300*300矩阵
-
St(2:301,2:301)=St(2:301,2:301)+Se.*(rand(300)>p)：其中Se.*(rand(300)>p)是将判断空地状态的矩阵与一个随机生成的矩阵相乘，只有当随机生成的元素大于设定的生长概率且这个位置是空地时，才能使得原有的空地变成未燃烧的树,即空地在大于p时长树，小于p时不长树
%%周边有八棵没有燃烧的树
%更新t时刻的森林St
Ss=zeros(302);%初始化Ss矩阵
Ss(Sk==1&Sk==0)=1;%%讨论绿树情况
Ss(2:301,2:301)=Ss(1:300,1:300)+Ss(1:300,2:301)+Ss(1:300,3:302) +...
Ss(2:301,1:300)+Ss(2:301,3:302)+Ss(3:302,1:300) + ...
Ss(3:302,2:301)+Ss(3:302,3:302);
%记录矩阵中的一个点的周围八个点一共有多少棵树
%这个式子求出的值就是一个点周围树的个数
Ss(Ss<7.5)=0; %周围的树小于8棵时，Ss赋值为0
Ss(Ss>7.5)=1; %周围的树等于8棵时，Ss赋值为1
d=find(Ss==1 & Sk==1);%查找非零元素的索引和值
%返回一个包含数组 X 中每个非零元素的线性索引的向量。- Ss是用来判断哪些格子周边没有燃烧的树而只有空地和未燃烧的树
for k=1:length(d)
r=rand(1);
St(d(k))=2*round((r>f));
end
%在最近的邻居中没有正在燃烧的树的情况下树在每一时步以概率f(闪电)变为正在燃烧的树-
St(d(k))=2*round((r>f))：未燃烧的树在大于f时燃烧，小于f时保持不燃烧
Sk=St; %更新t时刻的森林St
R=zeros(302); %初始化R矩阵
G=zeros(302); %初始化G矩阵
R(Sk==2)=1; %当Sk中标记为2的点，也是表示着火点则在R中标记出来
G(Sk==1)=1; %当Sk中标记为1的点，也是表示树的点在G中标记出来
C(:,:,1)=R; %把C的第一层赋值为R
C(:,:,2)=G; %把C的第二层赋值为G
set(Ci,'CData',C); %更换Ci的数据，Ci上面有定义，是显示图像，C是图像的矩阵
set(tp,'string',['T = ',num2str(ti)]); %第二行是设置图像的标题，显示T=当前时刻
pause(1);-
Sk=St：保存当前森林状态作为下次循环中的上一时刻森林状态 - 重新赋值RGB进行森林模型当前时刻的状态展示
- 对图像进行相应设置，并暂停1秒方便观察
相较于基础模型，新加入一个条件：考虑南风作用，对于燃烧的树有一定的几率会将上方距离某格子两格距离的格子的树点燃
w=0.8;- 考虑在风力影响的情况下树被点燃的概率为w
for i=2:301
for j=2:301
r0=rand(1);
if St(i,j)==2&&St(i-1,j)==1&&St(i-2,j)==1 %该树着火且前一棵有树
St(i-1,j)=round(2*(r0<=w)+(r0>w)); %风速够就着火，不够就不着火
St(i-2,j)=round(2*(r0<=w)+(r0>w)); %风速够就着火，不够就不着火
end
end
end
- 假设每辆车的车身和车距相同*
- 在绿灯亮时共有N辆车停车等待，第一辆车位于路口，其他车辆依次向后排列
- 车辆启动后，均匀的加速过程简化为每个时间单位多进一格，匀速过程为每个时间单位车辆向前移动4格【车辆启动后第1个时间单位,车向前移动一格；第2个时间单位，车向前移动2格；第3个时间单位，车向前移动3格；之后的每个时间单位，车辆移动4格】
N=20; %车辆总数
R=zeros(1,100); %初始化R矩阵
R(N+2:100)=1; %将R矩阵的N+2到100位初始化为1，“1”代表白色，“0”代表黑色，“0.5”代表灰色
R(N+1)=0.5; %将R矩阵第N+1位初始化为0.5，路口颜色
imshow(R,'InitialMagnification','fit'); %显示图像，并且图片适应屏幕的大小通过定义二维向量R来展示元胞自动机的仿真情况，其中：
- “1”（白色）代表该位置没有车
- “0”（黑色）代表该位置车
- “0.5”（灰色）代表该位置是路口停车线
C(1,:)=zeros(1,N); %表示时间变化（临时设置）
C(2,:)=(N:-1:1); %表示启动等待时间
C(3,:)=(1:N); %表示车辆的实时位置，此时表示初始位置
D=mat2cell(C,3,ones(1,N)); %生成元胞数组D，其中包含20组元胞数据，每组有3个元素，分别为实时时间、车辆的编号（按启动顺序）、车辆实时位置- 创建存储数据的元胞，包含实时时间、车辆的编号（按启动顺序）、车辆实时位置这三个数据，目前赋予的是初始状态的值，之后每次循环会对其进行更新
- 初始状态为：{0,20,1},{0,19,2},···,{0,1,20}
t=20; %总的运行时间
for k=1:t
pause(1); %暂停时间
%重新设置新一次的时间、位置并组成元胞
C(1,:)=k; %更新当前时间
temp1=C(3,:); %记录原来车的位置
D=mat2cell(C,3,ones(1,N));
temp2=cellfun(@traffic_flow_function,D); %表示实时更新车辆的位置-
t=20：设定总的循环运行次数为20次 -
pause(1)暂停1秒，方便展示 - 生成D作为当前车辆数据（当前时间、车辆编号、车辆实时位置），用来作为
traffic_flow_function函数的参数，从而确定下一时刻车辆的数据 -
temp2=cellfun(@traffic_flow_function,D)：其中cellfun的作用是将函数traffic_flow_function应用于元胞数组D的每个元胞的内容，而traffic_flow_function函数是导出下一时刻车辆数据用的。
function y=traffic_flow_function(E)
%E(1)表示当前的时间；E(2)表示车子的位置编号；E(3)是存储车辆的位置，返回的是车辆更新后的位置
k=E(1);%将当前的时间赋给k
if E(2)==k%如果E(2)=k则表明车子刚刚启动
E(3)=E(3)+1;%车辆位置加1
else
if E(2)==k-1%如果E(2)=k-1则表明车子已经启动了1秒
E(3)=E(3)+2;%车辆位置加2
else
if E(2)==k-2%如果E(2)=k-2则表明车子已经启动了2秒
E(3)=E(3)+3;%车辆位置加3
else
if E(2)<=k-3%如果E(2)<=k-1则表明车子已经启动3秒以上
E(3)=E(3)+4;%车辆位置加4
end
end
end
end
y=E(3);%输出车辆现在的位置
end- 编写的函数时根据已有的元胞自动机的条件规则：车辆启动后第1个时间单位,车向前移动一格；第2个时间单位，车向前移动2格；第3个时间单位，车向前移动3格；之后的每个时间单位，车辆移动4格。
for j=1:N
R(temp1(j))=1; %表示释放原来车的位置
R(temp2(j))=0; %表示更新车的实时位置
R(21)=0.5; %如果车到路口的位置，保持路口的颜色标志
end- 第一步：把每个位置全部置“1”（白色）【没有车】
- 第二步：把有车的位置置“0”（黑色）【有车】
- 第三步：把路口停车线位置置“0.5”（灰色）【路口停车线】
C(3,:)=temp2; %更新状态矩阵，给下一次使用
imshow(R,'InitialMagnification','fit'); %图片适应屏幕的大小-
C(3,:)=temp2：更新车辆实时位置的数据，为下一次循环时的使用做准备 -
imshow(R,'InitialMagnification','fit')：展示模型的实时情况
- 车辆从左侧路段进入现实的路段区域，为单输入的交通流模型
- 车辆可以随机选择在交叉路口处往哪个方向（向前、向左、向右）
- 当绿灯亮起时，车辆正常通行
- 当红灯亮起时，非左侧路端的车辆保持正常通行；左侧路段的车辆根据规则：如果车辆位于路口则停下；如果未到达路口但是前面有车则停下，如果未到达路口且前面没车则正常通行
%% 图框设置
figure;%利用当前属性创建窗格
axes;%建立坐标轴
rand('state',0);
set(gcf,'DoubleBuffer','on');-
rand('state',0)的作用是：使得每次产生的随机数是一致的，从而使得在包含随机数的程序每次执行时最初的样子是保持一致的 -
set(gcf,'DoubleBuffer','on')的作用是：防止在不断循环画动画的时候产生闪烁的现象
%% 交通图像初始化设置
% “1”白色——没车；“0”黑色——有车；“0.5”灰色——其他部分
road_length=50; %定义单条道路长度为100
road_width=1; %定义单条道路宽度为1
C=ones(2*road_length+1,2*road_length+1); %定义用于显示的图像矩阵
% 将除了路段的其他部分置为灰色
C(1:road_length,1:road_length)=0.5;
C(1:road_length,road_length+2:2*road_length+1)=0.5;
C(road_length+2:2*road_length+1,1:road_length)=0.5;
C(road_length+2:2*road_length+1,road_length+2:2*road_length+1)=0.5;
Ci=imshow(C); %展示模拟情况
ti=0;%初始化时间
tp=title(['T = ',num2str(ti)]);%在标题处显示时间- 图像设置将其分为两个区域：路段区域初始时置位白色，其他区域初始时置位灰色
p=0.5;
% v=2;
light_time=20;
light_gp=0.5;-
p=0.5：车辆进入该区域的概率为0.5 -
light_time=20：红绿灯的循环时长为20秒 -
light_gp=0.5：绿灯时间占总循环时间的50%
infinite=10000;
Car=cell(infinite,1); %设置一个超大的元胞数组，用于存放车辆元胞数据，提高运算速率
car_number=0; %初始化车辆的数量
car_on=0; %初始化在显示区域的道路上的车辆的数量
judge=0; %标识符：用于判断是否有车将离开该区域-
while 1：无限循环交通流模型
ti=ti+1;
%% 红绿灯判断——先绿，再红
if mod(ti,light_time)<=light_time*light_gp
flag=1; %当时间小于绿灯时间时，flag=1（绿灯）
C(road_length+3,road_length-2)=1;%绿灯时为白
else
flag=0; %当时间大于绿灯时间时，flag=0（红灯）
C(road_length+3,road_length-2)=0;%红灯时为黑
end-
ti=ti+1：每个循环时间加1 -
flag：作为红绿灯的标识符，方便为之后判断而设定，红绿灯的原则是先绿再红，“0”为红灯，“1”为绿灯 -
mod(ti,light_time)<=light_time*light_gp：判断当前时间属于绿灯时间还是红灯时间 -
C(road_length+3,road_length-2)=1：使得其他区域考路段的一个格子上通过显示白色（1）和黑色（0）来表示绿灯和红灯
根据车辆的位置上和性质（设定的方向）来控制车辆运动状态
if ti~=1
%% 根据车辆位置和性质使车辆移动
%遍历所有道路上的车辆，获得道路上车辆的编号然后根据车辆数据元胞中的坐标进行移动即可
car_ID=SearchCarData(Car); %搜寻道路上存在的车辆，获得每个车辆的编号
for i=1:size(car_ID,1) %遍历所有在显示道路上的车辆
car_data=Car{car_ID(i)};
single_car_y=car_data(4);
single_car_x=car_data(3);
single_car_ID=car_ID(i); %ID需要是正着数
judge=JudgeMarginCar(single_car_y,single_car_x,road_length);
if judge==0
[C,Car]=MoveCar(single_car_y,single_car_x,single_car_ID,Car,C,road_length,flag);
else
judge=0;
[C,Car]=DeleteCar(single_car_y,single_car_x,single_car_ID,Car,C);
car_on=car_on-1;
end
end
end-
car_ID=SearchCarData(Car)：主函数中的语句 -
SearchCarData函数：先通过遍历车辆数据的元胞数组Car的元素，找出非空且第一个元素（车辆状态：代表是否在现实路段上）为“1”的元素位置并提取出来，给到car_ID中
function car_ID=SearchCarData(Car)
%遍历所有车辆，找出子元胞的第1个元素(车辆状态)为0且第2个元素（选择方向）不为零的项
index_1=cellfun(@search_1,Car); %符合项为“1”，不符合项为“0”
car_ID=find(index_1); %车辆在所有车辆中的编号
endfunction y=search_1(Car) %找出子元胞的第1个元素为1且非空的项，并赋值“1”，其余赋值“0”
if ~isempty(Car)&&Car(1)==1
y=1;
else
y=0;
end
endfor i=1:size(car_ID,1) %遍历所有在显示道路上的车辆
car_data=Car{car_ID(i)};
single_car_y=car_data(4);
single_car_x=car_data(3);
single_car_ID=car_ID(i);- 遍历道路上的所有车辆**【注意按编号1~N的顺序，前车先走后车再动】**，在每个循环中针对每辆车进行控制，首先得到它们的位置坐标和编号
-
judge=JudgeMarginCar(single_car_y,single_car_x,road_length)：主函数中的语句 -
JudgeMarginCar函数：根据每辆车的坐标判断其是否位于显示区域的边缘，如果其在边缘处，则标识符judge置1，否则置0
function judge=JudgeMarginCar(car_y,car_x,road_length)
if car_x==(2*road_length+1)||car_y==1||car_y==(2*road_length+1)
judge=1;
else
judge=0;
end
end-
[C,Car]=MoveCar(single_car_y,single_car_x,single_car_ID,Car,C,road_length,flag)：主函数中的语句 -
MoveCar函数：
function [C,Car]=MoveCar(car_y,car_x,car_ID,Car,C,road_length,flag)
direction=Car{car_ID}(2); %拿到车辆的方向数据
C(car_y,car_x)=1;
if flag==1||(flag==0&&car_x>road_length) %如果是绿灯，或者红灯情况下在非左侧路段，则正常通行
if car_x==(road_length+1)&&car_y==(road_length+1) %判断车的位置是否在路段端口中间
if direction==1
new_y=car_y;
new_x=car_x+1;
elseif direction==2
new_y=car_y-1;
new_x=car_x;
elseif direction==3
new_y=car_y+1;
new_x=car_x;
end
elseif car_x~=(road_length+1)&&car_y==(road_length+1) %判断车的位置是否在横向的路段
new_y=car_y;
new_x=car_x+1;
elseif car_x==(road_length+1)&&car_y<(road_length+1) %判断车的位置是否在向上的路段
new_y=car_y-1;
new_x=car_x;
elseif car_x==(road_length+1)&&car_y>(road_length+1) %判断车的位置是否在向下的路段
new_y=car_y+1;
new_x=car_x;
end
else %如果是红灯，则先判断有没有到路口，再判断能不能向前通行
if car_x==road_length&&car_y==(road_length+1) %判断车的位置是否到达路口处
new_y=car_y;
new_x=car_x;
elseif car_x<road_length&&car_y==(road_length+1) %判断车的位置未到达路口而在左侧路段
if C(car_y,car_x+1)==0
new_y=car_y;
new_x=car_x;
else
new_y=car_y;
new_x=car_x+1;
end
end
end
C(new_y,new_x)=0;
Car{car_ID}(3)=new_x;
Car{car_ID}(4)=new_y;
end实现车辆移动的步骤：
- 第一步：将该车辆原来的位置置1（白色）
- 第二步：判断红绿灯表示符
flag - 第三步：若绿灯则正常通行，若红灯则需要让非左侧车道车辆外其他车辆正常通行
- 第四步：红灯时，对于左侧车道车辆，当该车位于路口时需要停下；当该车未到路口但其前面有车时需要停下；当该车未到路口且前面没有车时则可以正常通行
- 第五步：将该车移动后的位置置0（黑色）
[C,Car]=DeleteCar(single_car_y,single_car_x,single_car_ID,Car,C);
car_on=car_on-1; %显示路段上的车辆数减1-
[C,Car]=DeleteCar(single_car_y,single_car_x,single_car_ID,Car,C)：主函数中的语句 -
DeleteCar函数：将该车辆数据元胞的第一个元素（状态参数）置0——离开显示路段，将该车当前位置置1（白色）——该车离开该位置
function [C,Car]=DeleteCar(car_y,car_x,car_ID,Car,C)
Car{car_ID}(1)=0; %使车辆的状态参数变为0表示已经到显示区域边界，接下来不再在显示区域的道路上
C(car_y,car_x)=1;
end[newcar_ID,car_number,car_on,C]=CreateNewCar(ti,road_length,car_number,car_on,C); %获取车辆数据（车辆状态、直行/左转/右转、车辆坐标（横纵）、车辆出现的时间）、更新后的总车辆数、道路上的总车辆数、总的图像数据
if ~isempty(newcar_ID)
Car{car_number}=newcar_ID;
end-
[newcar_ID,car_number,car_on,C]=CreateNewCar(ti,road_length,car_number,car_on,C)：主函数中的语句 -
CreateNewCar函数：
function [car_ID,car_number,car_on,C]=CreateNewCar(ti,road_length,car_number,car_on,C)
p=0.5;
if ti==1
C(road_length+1,1)=0; %如果是第一次循环则需要直接生成新车辆
car_state=1; %车辆状态——在道路上“1”
car_number=car_number+1; %车辆总数加1
car_on=car_on+1; %显示的道路上的车辆总数加1
car_ID=[car_state,ceil(3*rand(1)),1,road_length+1,ti]; %生成车辆时，生成该车辆数据（车辆状态、1直行/2左转/3右转、车辆坐标（横纵）、车辆出现的时间）
else
% 每次循环之后根据概率来判断是否生成新车辆
in_p=rand(1); %该循环时刻车辆进入的概率
if in_p>p
C(road_length+1,1)=0; %当有车辆进入的概率达到一定值时，在进入路口处生成一辆车
car_state=1; %车辆状态——在道路上“1”
car_number=car_number+1; %车辆总数加1
car_on=car_on+1; %显示的道路上的车辆总数加1
car_ID=[car_state,ceil(3*rand(1)),1,road_length+1,ti]; %生成车辆时，生成该车辆数据（车辆状态、直行/左转/右转、车辆坐标（横纵）、车辆出现的时间）
else
car_ID=[];
end
end
end- 第一次直接生成车辆，之后每次需要根据概率进行判断是否生成车辆
- 生成车辆时，生成该车辆数据（车辆状态、1直行/2左转/3右转、车辆坐标（横纵）、车辆出现的时间）
- 注意车辆总数
car_number和道路上的车辆总数car_on都需要加1
set(Ci,'CData',C);%第一行是显示图像，C是图像的矩阵，Ci上面有定义，是显示图像。
set(tp,'string',['T = ',num2str(ti)]);%第二行是设置图像的标题，显示T=当前时刻
pause(0.1);
- 左侧、上侧、右侧、下侧路段均可以有车辆随机进入
- 红、绿、黄灯控制车辆的运动，黄灯的引入可以有效地清理交叉路口的交通
- 一般情况的车辆优先级（除了在测试过程中发现的特殊情况需要特殊处理，以防交叉路口交通非正常堵塞）：左侧车辆>上侧车辆>右侧车辆>下侧车辆
