元胞自动机基于matlab元胞自动机单边教室疏散含Matlab源码 1207期
Posted 紫极神光
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了元胞自动机基于matlab元胞自动机单边教室疏散含Matlab源码 1207期相关的知识，希望对你有一定的参考价值。
一、元胞自动机简介
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(~)
%画出动图
A=zeros(113,26);
M=ceil(sqrt(113^2+26^2));
A1=zeros(29,22);
A1(:,[1 22])=M;
A1([1 29],:)=M;
A1([4 5],[10 11 12 13])=M;%给出讲台的障碍量
A1(8:2:24,[4 5 6 10 11 12 13 17 18 19])=M;%给出桌子的障碍量
A1(9:2:25,[5 11 12 18])=A(9:2:25,[5 11 12 18]);%给出位于邻接座位中间位子的距门口的距离,这里去除了+1的情况
A1(5,1)=0;%给出门的位置
A([1,29,57,85],5:26)=M;
A(113,:)=M;
A(:,[1,5,26])=M;
A(1:29,5:26)=A1;
A(29:57,5:26)=A1;
A(57:85,5:26)=A1;
A(85:113,5:26)=A1;%至此给出了教室的形状，障碍物的分布，门口的位置
C=zeros(113,26);
C1=zeros(29,22);
C1(9:2:25,[4 5 6 10 11 12 13 17 18 19])=M;
C(85:113,5:26)=C1;%至此给出了整个楼层的人员分布情况表如矩阵C所示，此时给出的是初始情况，每个人都坐在座位上
for g=1:340;%200次肯定可以结束，可以尝试用while，这里懒得了……
imshow(max(A,C)==0,'InitialMagnification','fit')%取最大的作为人员在教室的分布，以及给出教室轮廓和障碍物图像
%其后两个参数是为了调整图像的大小
pause(0.25);
C=run(C);
end
function [D]=choic(D)
A=zeros(29,22);%生成一个29行22列的零矩阵，并且处理好后的A矩阵将表示每个元胞对门口的渴望程度
B=ones(25,22);%生成25行22列的每个位置是1的矩阵，这个是为了处理
A(5:29,:)=sqrt((cumsum(B,1)-1).^2+(cumsum(B,2)-1).^2);%运用计算技巧计算好每个元胞距门口的距离
C=ones(4,22);%生成4行22列的每个位置是1的矩阵
A(1:4,:)=sqrt((cumsum(C,1)-5).^2+(cumsum(C,2)-1).^2);%运用计算技巧计算好每个元胞距门口的距离
M=ceil(sqrt(113^2+26^2));%取整个楼层的对角线的向上取整作为最远的距离，这里也是编写完走廊之后为了将两个程序拼写在一起才做的微调，否则取房间的对角线向上取整即可
A(:,[1 22])=M;
A([1 29],:)=M;%给出墙壁的障碍量
A([4 5],[10 11 12 13])=M;%给出讲台的障碍量
A(8:2:24,[4 5 6 10 11 12 13 17 18 19])=M;%给出桌子的障碍量
A(9:2:25,[5 11 12 18])=A(9:2:25,[5 11 12 18]);%给出位于邻接座位中间位子的距门口的距离,这里去除了+1的情况
A(5,1)=0;%给出门口的位置
%D=zeros(29,22);%另给出一张表，表示在同样大小的教室
%D(9:2:25,[4 5 6 10 11 12 13 17 18 19])=M;%给出教室中人员的位置
%在这里给出矩阵D是为了验证程序写的正确性时便于调用
X=zeros(29,22); %给出目标矩阵，即人作出的选择，下一刻，即将要走向的位置
Y=zeros(29,22); %给出想要运动到X中的同一目标（x,y）的人所在的位置,其上值为x*sqrt(2)+y*sqrt(3)表示目标要到（x,y）去
Z=rand(29,22);%生成一张随机数表，用于最后的决策判断，注意这个表所应该放置的位置
H=zeros(29,22);%生成矩阵H式为了作为最后成功抢到位置的人所在的位置
%A(5,1)=0;%给出初始值
%D(5,1)=0;%给出D的门口的位置
%X(5,1)=0;
%Y(5,1)=0;
%Z(5,1)=0;
%这些量是在单独考虑一个房间时使房间门口的人走出去，如果考虑整个楼层，则使教室的人走出去的工作放在run函数中了
for x=2:28;%给出x的循环范围
for y=2:21;%给出y的循环范围
if D(x,y)==M %九宫格的中心如果是人的话才进行如下操作
E=max(A,D);%取A和D中对应值的最大的表示此刻人也在教室，同时除本人所在的位置之外的每个元胞的值表示此人距门口的距离
F=E((x-1):1:(x+1),(y-1):(y+1));%取出以（x,y）为中心的九宫格
F(2,2)=A(x,y);%九宫格的中心对别人有影响，对自己没有，所以改回原来的值
G=sort(F);%对列排序
b=G(1,:);%取行向量
d=sort(b);%对行向量排序
if length(find(F==d(1)))==1 %先考虑找到的最小值为1个的情况
[r,c]=find(F==d(1)); %找到最小值的位置
if r==2&&c==2 %最小值恰好在中心
X(x,y)=D(x,y);%目标矩阵是其本身，保持不变
Y(x,y)=x*sqrt(2)+y*sqrt(3);%给出到达目标矩阵（x,y）的位置上的值
else p=x-2+r; %最小值不在中心的时候找到的最小值位置即为目标位置
q=y-2+c;
if p~=1
X(p,q)=M;%给出目标矩阵
Y(x,y)=p*sqrt(2)+q*sqrt(3);%给出到达目标矩阵（x,y）的位置上的值
else X(x,y)=D(x,y);
Y(x,y)=x*sqrt(2)+y*sqrt(3);
end
end
else length(find(F==d(1)))==2; %考虑同时出现两个最小值得情况
[r,c]=find(F==d(1));
s=rand(1);
if s>0.5
p=x-2+r(1);
q=y-2+c(1);
if p~=1;
X(p,q)=M;
Y(x,y)=p*sqrt(2)+q*sqrt(3);
else X(x,y)=D(x,y);
Y(x,y)=x*sqrt(2)+y*sqrt(3);
end
else
p=x-2+r(2);
q=y-2+c(2);
if p~=1;
X(p,q)=M;
Y(x,y)=p*sqrt(2)+q*sqrt(3);
else X(x,y)=D(x,y);
Y(x,y)=x*sqrt(2)+y*sqrt(3);
end
end
end
end
end
end
for x=2:28;%给出x的循环范围
for y=2:22;%给出y的循环范围
if X(x,y)>0;
Y1=Y((x-1):1:(x+1),(y-1):(y+1));%对目标（x,y）给出九宫格
Y2=Z((x-1):1:(x+1),(y-1):(y+1));%给出九宫格的随机数作为之后的判断
w=x*sqrt(2)+y*sqrt(3);
%(Y1==w)*Y2;%从九宫格中找到以中心为目标的人所对应的随机数
t=max(max((Y1==w).*Y2));%找到以中心为目标的人所对应位置的最大值
[r1,c1]=find((Y1==w).*Y2==t);%找到最大值的位置
end
end
%对以（5,1）为目标的应该如何讨论
for x=2:28;%给出x的循环范围
for y=1;%给出y的循环范围
if X(x,y)>0;
元胞自动机基于matlab元胞自动机甲板火灾下人群疏散含Matlab源码 1291期