【元胞自动机】基于matlab激进策略元胞自动机三车道(不开放辅路,软件园不影响)交通流模型【含Matlab源码 1296期】

688 阅读8分钟

一、元胞自动机简介

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矩阵
%              123 边界 
%   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;


alltime=0;
%mp(25,comy)=2;
cy2=130;
mp(25,cy2)=2;
for i = 1:T
    sumtime=sumtime+dt;
%   刷车
    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_road2(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)
 function [mp, v, time] = chos_road(mp,v,time)
%
% 
%  
% 
% 
% 
%

[N, M] = size(mp);
for y=M-1:-1:2
    for x=26:28
        if(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;
                else
                    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;
            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] = chos_road(mp,v,time)
%
% 
%  
% 
% 
% 
%

[N, M] = size(mp);
dis=zeros(N,M)-1;
for y=M:-1:2
    for x=26:28 
        if (mp(x,y)==2|mp(x,y)==3) 
            dis(x,y) = 0;
        end;
    end;
end;
for y=M-1:-1:2
    for x=26:28
        if (dis(x,y)==-1 && dis(x,y+1)~=-1 ) 
            dis(x,y) = dis(x,y+1)+1;
        end;
    end;
end;
for y=2:M-1
    for x=26:28
        if(mp(x,y)==2)
            if(dis(x,y+1)==-1)dis(x,y)=123;
            else dis(x,y)=dis(x,y+1)+1;
            end;
        end;
    end;
end;

for y=2:M-1
    for x=26:28
        if(mp(x,y)==2 & (dis(x,y)-1<v(x,y)||mp(x,y+1)==2))
            if(mp(x-1,y)==1 & dis(x-1,y)>dis(x,y)  & mp(x+1,y)==1 & dis(x+1,y)>=dis(x,y))
                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;
                else
                    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;
            elseif(mp(x+1,y)==1 & dis(x+1,y)>=dis(x,y))
                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 & dis(x-1,y)>dis(x,y))
                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;
     

三、运行结果

在这里插入图片描述

四、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)