【元胞自动机】基于matlab保守策略元胞自动机三车道（不开放辅路，软件园不影响）交通流模型【含Matlab源码 1293期】

一、获取代码方式

获取代码方式1：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。

获取代码方式2：
通过紫极神光博客主页开通CSDN会员，凭支付凭证，私信博主，可获得此代码。

获取代码方式3：
完整代码已上传我的资源：【元胞自动机】基于matlab保守策略元胞自动机三车道（不开放辅路，软件园不影响）交通流模型【含Matlab源码 1293期】

备注：开通CSDN会员，仅只能免费获得1份代码（有效期为开通日起，三天内有效）；
订阅紫极神光博客付费专栏，可免费获得2份代码（有效期为订阅日起，三天内有效）；

二、元胞自动机简介

1 元胞自动机发展历程
最初的元胞自动机是由冯 · 诺依曼在 1950 年代为模拟生物细胞的自我复制而提出的. 但是并未受到学术界重视.
1970 年, 剑桥大学的约翰 · 何顿 · 康威设计了一个电脑游戏 “生命游戏” 后, 元胞自动机才吸引了科学家们的注意.

1983 年 S.Wolfram 发表了一系列论文. 对初等元胞机 256 种规则所产生的模型进行了深入研究, 并用熵来描述其演化行为, 将细胞自动机分为平稳型, 周期型, 混沌型和复杂型.

2 对元胞自动机的初步认识
元胞自动机（CA）是一种用来仿真局部规则和局部联系的方法。典型的元胞自动机是定义在网格上的，每一个点上的网格代表一个元胞与一种有限的状态。变化规则适用于每一个元胞并且同时进行。典型的变化规则，决定于元胞的状态，以及其（ 4 或 8 ）邻居的状态。

3 元胞的变化规则&元胞状态
典型的变化规则，决定于元胞的状态，以及其（ 4 或 8 ）邻居的状态。

4 元胞自动机的应用
元胞自动机已被应用于物理模拟，生物模拟等领域。

5 元胞自动机的matlab编程
结合以上，我们可以理解元胞自动机仿真需要理解三点。一是元胞，在matlab中可以理解为矩阵中的一点或多点组成的方形块，一般我们用矩阵中的一点代表一个元胞。二是变化规则，元胞的变化规则决定元胞下一刻的状态。三是元胞的状态，元胞的状态是自定义的，通常是对立的状态，比如生物的存活状态或死亡状态，红灯或绿灯，该点有障碍物或者没有障碍物等等。

6 一维元胞自动机——交通规则
定义：
6.1 元胞分布于一维线性网格上.
6.2 元胞仅具有车和空两种状态.

7 二维元胞自动机——生命游戏
定义：
7.1 元胞分布于二维方型网格上.
7.2 元胞仅具有生和死两种状态.

元胞状态由周围八邻居决定.
规则：

骷髅：死亡；笑脸：生存
周围有三个笑脸，则中间变为笑脸
少于两个笑脸或者多于三个，中间则变死亡。

8 什么是元胞自动机
离散的系统: 元胞是定义在有限的时间和空间上的, 并且元胞的状态是有限.
动力学系统: 元胞自动机的举止行为具有动力学特征.
简单与复杂: 元胞自动机用简单规则控制相互作用的元胞模拟复杂世界.

9 构成要素

（1）元胞 (Cell)

元胞是元胞自动机基本单元：
状态: 每一个元胞都有记忆贮存状态的功能.
离散: 简单情况下, 元胞只有两种可能状态; 较复杂情况下, 元胞具有多种状态.
更新: 元胞的状态都安照动力规则不断更新.
（2）网格 (Lattice)
不同维网格

常用二维网格

（3）邻居 (Neighborhood)

（4）边界 (Boundary)

反射型：以自己作为边界的状态
吸收型：不管边界（车开到边界就消失）

（5）规则（状态转移函数）
定义：根据元胞当前状态及其邻居状况确定下一时刻该元胞状态的动力学函数, 简单讲, 就是一个状态转移函数.
分类：
总和型: 某元胞下时刻的状态取决于且仅取决于它所有邻居的当前状态以及自身的当前状态.
合法型: 总和型规则属于合法型规则. 但如果把元胞自动机的规则限制为总和型, 会使元胞自动机具有局限性.
（6）森林火灾

绿色：树木；红色：火；黑色：空地。
三种状态循环转化：
树：周围有火或者被闪电击中就变成火。
空地：以概率p变为树木
理性分析：红为火；灰为空地；绿是树

元胞三种状态的密度和为1

火转化为空地的密度等于空地转换为树的密度（新长出来的树等于烧没的树）

f是闪电的概率:远远小于树生成的概率；T s m a x T_{smax}T smax
是一大群树被火烧的时间尺度
程序实现
周期性边界条件
购进啊

其中的数字为编号
构建邻居矩阵

上面矩阵中的数字编号，对应原矩阵相同位置编号的上邻居编号，一一对应
同样道理：

（7）交通概念
车距和密度

流量方程

守恒方程

时空轨迹（横轴是空间纵轴为时间）

红线横线与蓝色交点表示每个时间车的位置。
如果是竖线则表示车子在该位置对应的时间

宏观连续模型：

最常用的规则：

红色条表示速度是满的。

1 加速规则：不能超过v m a x （ 2 格 / s ） v_{max}（2格/s）v
max（2格/s）
2 防止碰撞：不能超过车距

理论分析：

结果分析: 密度与流量

第一个图：横坐标是归一化后的密度，纵坐标是车流量。第二个图：理论值与CA的结果

结果分析: 时空轨迹

中间的深色区域是交通堵塞的区域。

三、部分源代码

function main
%   元胞自动机模拟三车道干路中间有大企业两边有一一车道辅路高峰时段
%   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Sidex1        Sidex2
%         +---------------+     Sidey = 2
%         |               |
%         |               |
%   ======+===============+==== Mainy
%   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   T               总模拟实际时间
%   Arrival         总车流量
%   Mainlength      干道长
%   Mainy           干道纵坐标
%   Mainwide        干道宽
%   Sidewide        辅路宽
%   Sidex1 Sidex2   辅路横坐标
%   Sidey           辅路纵坐标
%   Sidelengthh     辅路南北长
%   Sidelengthw     辅路东西长
%   MP              map矩阵
%              1 空 2 车 3 边界
%   lightmp         红绿灯
%   lighttime       红绿灯间隔
%   v               车速度矩阵
%   vmainmax        干道最高速
%   vsidemax        辅路最高速
%   time            车花费时间矩阵
%   dt              单位时间
%   tt              图像刷新时间
%   h               handle
%   sumtime         车总时间
%   alltime         车总花费时间
%   changepro       右侧车道进入铺路概率
clear;clc
T = 3000;
Mainlength = 200; Mainwide = 3; vmainmax = 6;
%Sidelengthh = 25; Sidelengthw = 100; Sidewide = 1;
changepro = 0;
vsidemax = 4;
Sidex1 = 30; Sidex2 = 150; Sidey = 2; Mainy = 26;
[mp, v, time,lightmp] = init_mp(Sidex1,Sidex2,Sidey,Mainy,Mainlength,Mainwide);h = print_mp(mp, NaN, 0.1);
dt = 1.2; tt=0.1;
sumtime = 0; sumpassed =0;
lighttime = 20;%   刷车if(rem(i,2)&(mp(27,8)~=2||v(27,8)~=0))[mp, v] = new_cars(dt,mp,v,vmainmax); end;
%   软件园if(rem(i,2))mp(24,cy2)=1;else mp(24,cy2)=2;end;mp(26,cy2)=2;v(26,cy2)=0;%   变道 chos_road2 激进 chos_road 保守[mp, v, time] = chos_road(mp,v,time);
%   变速[mp, v, time] = change_speed(mp,v,time,vmainmax);flag = (mp(Mainy,Sidex2)==0)
%   辅道    if(flag) [mp, v, time] = fulu(mp,v,time,vsidemax,Sidex1,Sidex2,Sidey,changepro);end;
%   位移&消车[alltime, mp, v, time , sumpassed ] = move(alltime,sumpassed,mp,v,time,dt);
%   辅道 if(~flag)[mp, v, time] = fulu(mp,v,time,vsidemax,Sidex1,Sidex2,Sidey,changepro);end;
%   输出h = print_mp(mp,h,tt);xlabel({strcat(strcat('time = ',num2str(sumtime,'%.1f')),'s'),strcat('total number of vehicles =',num2str(sumpassed)),strcat('average passing time =',num2str(alltime/sumpassed,'%.4f')),'不开放辅路，软件园影响，保守策略'});
%     w= getframe;
%     imind  =frame2im(w);
%     [imind ,cm]=rgb2ind(imind,256);
%     if(i~=1)
%         imwrite(imind, cm, 'example2.gif','gif','WriteMode','append','DelayTime',0.01);
%     else
%         imwrite(imind, cm, 'example2.gif','gif','Loopcount',inf,'DelayTime',0.01);
%     end;
endfunction [mp, v, time] = chos_road(mp,v,time)
%
%
%
%
%
%
%[N, M] = size(mp);
for y=M-1:-1:2for x=26:28if(mp(x,y)==2 & mp(x,y+1)~=1)if(mp(x-1,y)==1 & mp(x+1,y)==1 )if(rand<0.5)mp(x+1,y)=2;mp(x,y)=1;v(x+1,y)=v(x,y);v(x,y)=0;time(x+1,y)=time(x,y);time(x,y)=0;elsemp(x-1,y)=2;mp(x,y)=1;v(x-1,y)=v(x,y);v(x,y)=0;time(x-1,y)=time(x,y);time(x,y)=0;end;elseif(mp(x+1,y)==1)mp(x+1,y)=2;mp(x,y)=1;v(x+1,y)=v(x,y);v(x,y)=0;time(x+1,y)=time(x,y);time(x,y)=0;elseif(mp(x-1,y)==1)mp(x-1,y)=2;mp(x,y)=1;v(x-1,y)=v(x,y);v(x,y)=0;time(x-1,y)=time(x,y);time(x,y)=0;end;end;end;
end;function [mp, v, time] = fulu(mp,v,time,vsidemax,Sidey1,Sidey2,Sidex,changepro)
%
%
%
speedpro =0.9;
[N, M] = size(mp);
step=0; x=25; y=Sidey2;
while(1)step=step+1;if(mp(x,y)==2) break;end;if(x==Sidex) break;end;x=x-1;
end;
if(mp(x,y)==2&&v(x,y)>=step)v(25,y)=1; mp(25,y)=2; time(25,y)=time(x,y)+1.2;v(x,y)=0; mp(x,y)=1; time(x,y)=0;
end;
dis=zeros(N,M)-1;
x=25;y=Sidey2;
dis(26,Sidey2)=0;
while(1)if(x==26) break;end;if(mp(x,y)==2) dis(x,y)=0;end;if(y==Sidey2&&x~=Sidex)x=x-1;elseif(x==Sidex&&y~=Sidey1)y=y-1;elseif(y==Sidey1&&x~=26)x=x+1;end;
end;
lastx=26; lasty=Sidey2; x=25; y=Sidey2;
while(1)if(x==26) break;end;if(dis(x,y)==-1) dis(x,y)=dis(lastx,lasty)+1;end;lastx=x; lasty=y;if(y==Sidey2&&x~=Sidex)x=x-1;  elseif(x==Sidex&&y~=Sidey1)y=y-1;elseif(y==Sidey1&&x~=26)x=x+1;end;
end;
x=25; y=Sidey1; nxtx=0; nxty=0;
while(1)if(x==26) break;end;if(y==Sidey1&&x~=Sidex)nxtx=x-1; nxty=y;elseif(x==Sidex&&y~=Sidey2)nxty=y+1; nxtx=x;elseif(y==Sidey2&&x~=26)nxtx=x+1; nxty=y;end;if(dis(x,y)==0)dis(x,y)=dis(nxtx,nxty)+1;end;x=nxtx; y=nxty;
end;
x=25; y=Sidey2;
tot=0;
X=1:500;Y=1:500;
while(1)if(x==26) break;end; tot=tot+1;if(mp(x,y)==2)if(rand<speedpro) rlv = min(vsidemax,min(v(x,y)+1,dis(x,y)-1));else rlv= min(vsidemax,min(max(0,v(x,y)-1),dis(x,y)-1));end;if(rlv~=0) tox = X(tot-rlv); toy=Y(tot-rlv);mp(tox,toy)=2; v(tox,toy)=rlv;time(tox,toy)=time(x,y)+1.2; time(x,y)=0;mp(x,y)=1; v(x,y)=0;end;end;X(tot)=x; Y(tot)=y;if(y==Sidey2&&x~=Sidex)x=x-1;  elseif(x==Sidex&&y~=Sidey1)y=y-1;elseif(y==Sidey1&&x~=26)x=x+1;end;
end;

四、运行结果

五、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 蔡利梅.MATLAB图像处理——理论、算法与实例分析[M].清华大学出版社，2020.
[2]杨丹,赵海滨,龙哲.MATLAB图像处理实例详解[M].清华大学出版社，2013.
[3]周品.MATLAB图像处理与图形用户界面设计[M].清华大学出版社，2013.
[4]刘成龙.精通MATLAB图像处理[M].清华大学出版社，2015.
[5]孟逸凡,柳益君.基于PCA-SVM的人脸识别方法研究[J].科技视界. 2021,(07)
[6]张娜,刘坤,韩美林,陈晨.一种基于PCA和LDA融合的人脸识别算法研究[J].电子测量技术. 2020,43(13)
[7]陈艳.基于BP神经网络的人脸识别方法分析[J].信息与电脑(理论版). 2020,32(23)
[8]戴骊融,陈万米,郭盛.基于肤色模型和SURF算法的人脸识别研究[J].工业控制计算机. 2014,27(02)