元胞自动机

2019美赛D题的算法之一

一、题目分析

题目有关信息：
1、首先需要开发紧急疏散模型，以便从博物馆撤离访客，同时还允许应急人员尽快进入建筑物。重要的是找出可能限制出口移动的潜在瓶颈。
2、模型应具备很强的适应性，可用于解决广泛的考虑因素和各种类型的潜在威胁。
3、开发完成以后，验证模型并讨论卢浮宫如何实施它。
4、 根据工作成果，提出有关卢浮宫应急管理的政策和程序建议。包括认为对访客安全所必需的任何适用的人群管理和控制程序。另外，讨论如何为其他大型拥挤结构调整和实施模型
对于给出的题目信息，分析可知关键点是完成以下三个任务，首先是应急人员的进入问题，其次是如何获取管内人员的数量和位置信息，最后则是最佳路线的选定，采用什么样的疏散模型。

二、模型建立及解决

根据上面可以知道题目的重点是疏散模型的建立。那么，一般人群疏散的时候是根据什么来进行移动的呢。第一点是人群的走向移动，第二点是根据疏散通道的方向移动，第三是根据应急人员的指挥移动。因此，选用元胞自动机就可以很好的解决这一问题。把每个人看做一个单独的元胞，制定相应的元胞规则，由局部影响到整体，由此可以建立起一个简单的人员疏散模型。

三、元胞自动机

1、概念和定义

定义：
元胞自动机(cellular automata，CA) 是一种时间、空间、状态都离散，空间相互作用和时间因果关系为局部的网格动力学模型，具有模拟复杂系统时空演化过程的能力。
概念：
不同于一般的动力学模型，元胞自动机不是由严格定义的物理方程或函数确定，而是用一系列模型构造的规则构成。凡是满足这些规则的模型都可以算作是元胞自动机模型。因此，元胞自动机是一类模型的总称，或者说是一个方法框架。其特点是时间、空间、状态都离散，每个变量只取有限多个状态，且其状态改变的规则在时间和空间上都是局部的。
1)、元胞构成：
元胞：
· 元胞自动机最基本的单元
· 元胞有记忆储存状态的功能
· 所有元胞状态都按照元胞规则不断更新
格子：
· 元胞的网格空间
2)、元胞行为：
局部变化引起全局变化:
· 可简单认为元胞自动机在运动上类似与波
· 元胞的状态变化依赖于自身状态和邻居状态
元胞自动机的规则：
某元胞下时刻的状态只决定于邻居的状态以及自身的初始状态
元胞自动机的规则决定了元胞的行为特征
3)、元胞特征：离散的网格 元胞的同质 离散的状态局部的作用 离散的时间

2、算法

在和题目联系起来之前，先介绍常用的几个元胞自动机模型实例。
这里借鉴了一篇博客：https://blog.csdn.net/qq_40527086/article/details/86798384

1)生命游戏机

生命游戏是英国数学家约翰·何顿·康威在1970年发明的细胞自动机。它包括一个二维矩形世界，这个世界中的每个方格居住着一个活着的或死了的细胞。一个细胞在下一个时刻生死取决于相邻八个方格中活着的或死了的细胞的数量。通常情况，游戏的规则就是：当一个方格周围有2或3个活细胞时，方格中的活细胞在下一个时刻继续存活；即使这个时刻方格中没有活细胞，在下一个时刻也会“诞生”活细胞。
元胞规则：
· 对周围的8个近邻的元胞状态求和
· 如果总和为2的话，则下一时刻的状态不改变
· 如果总和为3，则下一时刻的状态为1，否则状态=0
动态展示元胞的变化过程，需要使用下面的代码：

h = imagesc(你的元胞矩阵)
在按照一定规则更行状态的循环中，
set(h,'cdata',你的更新后的元胞矩阵)
drawnow或者pause(时间间隔)

其余代码：

clc,clear
%% 设置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;%小于0.6计为0
% 建立图像句柄
imh = image(cat(3,cells,z,z));
set(imh, 'erasemode', 'none')%设置图形句柄为imh的图片擦写模式为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

2)扩散

元胞规则： 先把中间点置为1，每一时间步对每一点，如果周围
八个点和为偶数，则变为0，为奇数则变为1
代码：

clc,clear ,close all
%%
% 颜色控制
Map = [1 1 1; 0 0 0];
colormap(Map);
% 设置网格大小
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

3、题目相关

1)基本思想

通过对题目有关数据分析，制定相应的规则，使用元胞自动机模拟多出口疏散模型。
为了模拟聚集现象，必须引入长范围作用。为了使模型简单，而尽量避免使用长范围作用，通过地磁场建立行人与系统结构(建筑物)的相互关系，反映之间的相互作用。地磁场可以看作一个元胞网格空间，用来描述行人占用的元胞位置。为了很好的模拟聚集现象，这里设计两个场，一个是动态场 ，一个是静态场 。
动态场：动态场 用来描述行人的踪迹，动态场有自己的特性，可以表示扩散和衰退现象，随着行人的运动而动态更新。动态场 用来建模行人间的相互作用。行人走过的位置会留下一些信息即踪迹，其他人倾向于按照别人走过的路线行进。踪迹信息在每一时间步长内，会按照一定的规律扩散到相临元胞，自身也会衰减。
静态场：静态场 用来定义卢浮宫内部结构，其不随着时间和行人的运动而变化。静态场可以表示展区中的一些特定区域，如出口门等。通过静态场可以引导行人向着这些区域运动。当然也可以通过具体函数定义静态场。

2)元胞规则

· 每个网格的边长为0.4m
· 地图矩阵中，0.5代表行人，0代表墙壁，1代表馆外空地，2代表馆内空地
· 四个主要通道对应矩阵的位置分别为：（96~99,56），（210,99~105），(278,99~105),(78~83,140)。 （因为大门一次性能够容纳通过的人数是不一致的，所以通道宽度也不同）
· 游客优先朝着离自己最近的大门走去
· 行人走过的位置会留下一些信息即踪迹，其他人倾向于按照别人走过的路线行进。踪迹信息在每一时间步长内，会按照一定的规律扩散到相临元胞，自身也会衰减。

3)部分代码

（只实现了静态场部分，元胞规则第四条没实现）

clc,clear,close all;
set(gcf,'DoubleBuffer','on');
n=0;
I = imread('E:\matlab2018\bin\tuxiang\firstfloor.png');
O = imread('E:\matlab2018\bin\tuxiang\firstfloor4.png');
O1 = im2bw(O)+0;imshow(O1)
I1 = im2bw(I)+0;imshow(I1)
IO = I1+O1;
N = 563;
M = 985;
S = O1;
%% 添加人数
u = rand(N,M);
for i = 1:N
    for j = 1:M
        if u(i,j)<=0.005 %改变总人数
            u(i,j) = -0.5;
        else
            u(i,j)=1;
        end
    end
end
for i = 1:N
    for j = 1:M
        if IO(i,j)==2
            S(i,j)=u(i,j)+O1(i,j);
        end
    end
end
%%
%计算每个人到出口的欧式距离
X = [916,694,324,266];%保存每个大门中心点的像素
Y = [343,337,190,464];
Dg = zeros(1,4);
Dmin = 0;
Sm = zeros(N,M);%状态矩阵
Ds = zeros(N,M);%记录所有点到大门的最短距离
%% 静态计算每个点的状态矩阵
for i = 1:N
    for j = 1:M
        %针对每个点
        if S(i,j) == 0.5 %判断是否为疏散人群
            for k = 1:4%计算到每个大门中心点的距离
                Dx = abs(i-Y(k)).^2;
                Dy = abs(j-k+X(k)).^2;
                Dg(k) = sqrt(Dx+Dy);
            end
          %保存每个点到所有大门的最小值
            Dmin = min(Dg);        
            for k =1:4
               if Dg(k)== Dmin
                   Sm(i,j)=k;%记录状态
                   Ds(i,j) = Dmin;
               end
            end
        end
    end
end
Ii = imshow(S);
%% 动态更新
while(1)
for i=1:N
    for j = 1:M 
        Skm = Sm+S;
        if Sm(i,j)==1 %第一个出口
            %距离
            D1 = sqrt(abs(i+1-X(1)).^2+abs(j-Y(1)).^2);%右方向
            D2 = sqrt(abs(i-1-X(1)).^2+abs(j-Y(1)).^2);%左方
            D3 = sqrt(abs(i-X(1)).^2+abs(j+1-Y(1)).^2);%上方
            D4 = sqrt(abs(i-X(1)).^2+abs(j-1-Y(1)).^2);%下方
            D5 = sqrt(abs(i+1-X(1)).^2+abs(j+1-Y(1)).^2);%右上
            D6 = sqrt(abs(i+1-X(1)).^2+abs(j-1-Y(1)).^2);%右下
            D7 = sqrt(abs(i-1-X(1)).^2+abs(j+1-Y(1)).^2);%左上
            D8 = sqrt(abs(i-1-X(1)).^2+abs(j-1-Y(1)).^2); %左下       
            %判断周围情况
            if Skm(i+1,j)~=2
                D1 = inf;  
            end
            if Skm(i-1,j)~=2
                D2 = inf;   
            end
            if Skm(i,j+1)~=2
                D3 = inf;  
            end
            if Skm(i,j-1)~=2
                D4 = inf;  
            end
            if Skm(i+1,j+1)~=2
                D5 = inf;   
            end
            if Skm(i+1,j-1)~=2
                D6 = inf;     
            end
            if Skm(i-1,j+1)~=2
                D7 = inf;   
            end
            if Skm(i-1,j-1)~=2
                D8 = inf;
            end      
            
            Dt = [D1,D2,D3,D4,D5,D6,D7,D8];
            Dmin1 = min(Dt);
            for k = 1:8
                if Dt(k)==Dmin1
                    break
                end
                if Dt(k)==inf
                    k = 9;
                end
            end           
            switch k
                case k==1
                    tran = Sm(i,j);%状态矩阵改变
                    Sm(i,j) = Sm(i+1,j);
                    Sm(i+1,j) = tran;  
                    
                    tran2 = Ds(i,j);%距离矩阵改变
                    Ds(i,j) = Ds(i+1,j);
                    Ds(i+1,j) = tran2; 
                    
                    tran3 = S(i,j);%状态矩阵改变
                    S(i,j) = S(i+1,j);
                    S(i+1,j) = tran3; 
               case k==2
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i-1,j);
                    Sm(i-1,j) = tran; 
                    
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i-1,j);
                    Ds(i-1,j) = tran2;
                    tran3 = S(i,j);
                    S(i,j) = S(i-1,j);
                    S(i-1,j) = tran3;                 
               case k==3
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i,j+1);
                    Sm(i,j+1) = tran;
                    
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i,j+1);
                    Ds(i,j+1) = tran2;
                    
                    tran3 = S(i,j);
                    S(i,j) = S(i,j+1);
                    S(i,j+1) = tran3;
                    
               case k==4
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i,j-1);
                    Sm(i,j-1) = tran; 
                    
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i,j-1);
                    Ds(i,j-1) = tran2;  
                    tran3 = S(i,j);
                    S(i,j) = S(i,j-1);
                    S(i,j-1) = tran3; 
                    
               case k==5
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i+1,j+1);
                    Sm(i+1,j+1) = tran;    
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i+1,j+1);
                    Ds(i+1,j+1) = tran2; 
                    tran3 = S(i,j);
                    S(i,j) = S(i+1,j+1);
                    S(i+1,j+1) = tran3;                   
               case k==6
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i+1,j-1);
                    Sm(i+1,j-1) = tran;  
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i+1,j-1);
                    Ds(i+1,j-1) = tran2;   
                    tran3 = S(i,j);
                    S(i,j) = S(i+1,j-1);
                    S(i+1,j-1) = tran3;                     
               case k==7
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i-1,j+1);
                    Sm(i-1,j+1) = tran; 
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i-1,j+1);
                    Ds(i-1,j+1) = tran2;    
                    tran3 = S(i,j);
                    S(i,j) = S(i-1,j+1);
                    S(i-1,j+1) = tran3;                     
               case k==8
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i-1,j-1);
                    Sm(i-1,j-1) = tran; 
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i-1,j-1);
                    Ds(i-1,j-1) = tran2;  
                    tran3 = S(i,j);
                    S(i,j) = S(i-1,j-1);
                    S(i-1,j-1) = tran3; 
                otherwise
                    Sm(i,j) = Sm(i,j);
                    S(i,j) = S(i,j);
                    Ds(i,j) = Ds(i,j);              
            end           
        end
        if Sm(i,j)==3 %第三个出口
            %距离
            D1 = sqrt(abs(i+1-X(3)).^2+abs(j-Y(3)).^2);%右方向
            D2 = sqrt(abs(i-1-X(3)).^2+abs(j-Y(3)).^2);%左方
            D3 = sqrt(abs(i-X(3)).^2+abs(j+1-Y(3)).^2);%上方
            D4 = sqrt(abs(i-X(3)).^2+abs(j-1-Y(3)).^2);%下方
            D5 = sqrt(abs(i+1-X(3)).^2+abs(j+1-Y(3)).^2);%右上
            D6 = sqrt(abs(i+1-X(3)).^2+abs(j-1-Y(3)).^2);%右下
            D7 = sqrt(abs(i-1-X(3)).^2+abs(j+1-Y(3)).^2);%左上
            D8 = sqrt(abs(i-1-X(3)).^2+abs(j-1-Y(3)).^2); %左下       
            %判断周围情况
            if Skm(i+1,j)~=2
                D1 = inf;  
            end
            if Skm(i-1,j)~=2
                D2 = inf;   
            end
            if Skm(i,j+1)~=2
                D3 = inf;  
            end
            if Skm(i,j-1)~=2
                D4 = inf;  
            end
            if Skm(i+1,j+1)~=2
                D5 = inf;   
            end
            if Skm(i+1,j-1)~=2
                D6 = inf;     
            end
            if Skm(i-1,j+1)~=2
                D7 = inf;   
            end
            if Skm(i-1,j-1)~=2
                D8 = inf;
            end      
            
            Dt = [D1,D2,D3,D4,D5,D6,D7,D8];
            Dmin1 = min(Dt);
            for k = 1:8
                if Dt(k)==Dmin1
                    break
                end
                if Dt(k)==inf
                    k = 9;
                end
            end           
            switch k
                case k==1
                    tran = Sm(i,j);%状态矩阵改变
                    Sm(i,j) = Sm(i+1,j);
                    Sm(i+1,j) = tran;  
                    
                    tran2 = Ds(i,j);%距离矩阵改变
                    Ds(i,j) = Ds(i+1,j);
                    Ds(i+1,j) = tran2; 
                    
                    tran3 = S(i,j);%状态矩阵改变
                    S(i,j) = S(i+1,j);
                    S(i+1,j) = tran3; 
               case k==2
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i-1,j);
                    Sm(i-1,j) = tran; 
                    
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i-1,j);
                    Ds(i-1,j) = tran2;
                    tran3 = S(i,j);
                    S(i,j) = S(i-1,j);
                    S(i-1,j) = tran3;                 
               case k==3
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i,j+1);
                    Sm(i,j+1) = tran;
                    
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i,j+1);
                    Ds(i,j+1) = tran2;
                    
                    tran3 = S(i,j);
                    S(i,j) = S(i,j+1);
                    S(i,j+1) = tran3;
                    
               case k==4
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i,j-1);
                    Sm(i,j-1) = tran; 
                    
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i,j-1);
                    Ds(i,j-1) = tran2;  
                    tran3 = S(i,j);
                    S(i,j) = S(i,j-1);
                    S(i,j-1) = tran3; 
                    
               case k==5
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i+1,j+1);
                    Sm(i+1,j+1) = tran;    
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i+1,j+1);
                    Ds(i+1,j+1) = tran2; 
                    tran3 = S(i,j);
                    S(i,j) = S(i+1,j+1);
                    S(i+1,j+1) = tran3;                   
               case k==6
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i+1,j-1);
                    Sm(i+1,j-1) = tran;  
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i+1,j-1);
                    Ds(i+1,j-1) = tran2;   
                    tran3 = S(i,j);
                    S(i,j) = S(i+1,j-1);
                    S(i+1,j-1) = tran3;                     
               case k==7
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i-1,j+1);
                    Sm(i-1,j+1) = tran; 
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i-1,j+1);
                    Ds(i-1,j+1) = tran2;    
                    tran3 = S(i,j);
                    S(i,j) = S(i-1,j+1);
                    S(i-1,j+1) = tran3;                     
               case k==8
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i-1,j-1);
                    Sm(i-1,j-1) = tran; 
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i-1,j-1);
                    Ds(i-1,j-1) = tran2;  
                    tran3 = S(i,j);
                    S(i,j) = S(i-1,j-1);
                    S(i-1,j-1) = tran3; 
                otherwise
                    Sm(i,j) = Sm(i,j);
                    S(i,j) = S(i,j);
                    Ds(i,j) = Ds(i,j);              
            end           
        end
        if Sm(i,j)==2 %第二个出口
            %距离
            D1 = sqrt(abs(i+1-X(2)).^2+abs(j-Y(2)).^2);%右方向
            D2 = sqrt(abs(i-1-X(2)).^2+abs(j-Y(2)).^2);%左方
            D3 = sqrt(abs(i-X(2)).^2+abs(j+1-Y(2)).^2);%上方
            D4 = sqrt(abs(i-X(2)).^2+abs(j-1-Y(2)).^2);%下方
            D5 = sqrt(abs(i+1-X(2)).^2+abs(j+1-Y(2)).^2);%右上
            D6 = sqrt(abs(i+1-X(2)).^2+abs(j-1-Y(2)).^2);%右下
            D7 = sqrt(abs(i-1-X(2)).^2+abs(j+1-Y(2)).^2);%左上
            D8 = sqrt(abs(i-1-X(2)).^2+abs(j-1-Y(2)).^2); %左下       
            %判断周围情况
            if Skm(i+1,j)~=2
                D1 = inf;  
            end
            if Skm(i-1,j)~=2
                D2 = inf;   
            end
            if Skm(i,j+1)~=2
                D3 = inf;  
            end
            if Skm(i,j-1)~=2
                D4 = inf;  
            end
            if Skm(i+1,j+1)~=2
                D5 = inf;   
            end
            if Skm(i+1,j-1)~=2
                D6 = inf;     
            end
            if Skm(i-1,j+1)~=2
                D7 = inf;   
            end
            if Skm(i-1,j-1)~=2
                D8 = inf;
            end      
            
            Dt = [D1,D2,D3,D4,D5,D6,D7,D8];
            Dmin1 = min(Dt);
            for k = 1:8
                if Dt(k)==Dmin1
                    break
                end
                if Dt(k)==inf
                    k = 9;
                end
            end           
            switch k
                case k==1
                    tran = Sm(i,j);%状态矩阵改变
                    Sm(i,j) = Sm(i+1,j);
                    Sm(i+1,j) = tran;  
                    
                    tran2 = Ds(i,j);%距离矩阵改变
                    Ds(i,j) = Ds(i+1,j);
                    Ds(i+1,j) = tran2; 
                    
                    tran3 = S(i,j);%状态矩阵改变
                    S(i,j) = S(i+1,j);
                    S(i+1,j) = tran3; 
               case k==2
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i-1,j);
                    Sm(i-1,j) = tran; 
                    
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i-1,j);
                    Ds(i-1,j) = tran2;
                    tran3 = S(i,j);
                    S(i,j) = S(i-1,j);
                    S(i-1,j) = tran3;                 
               case k==3
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i,j+1);
                    Sm(i,j+1) = tran;
                    
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i,j+1);
                    Ds(i,j+1) = tran2;
                    
                    tran3 = S(i,j);
                    S(i,j) = S(i,j+1);
                    S(i,j+1) = tran3;
                    
               case k==4
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i,j-1);
                    Sm(i,j-1) = tran; 
                    
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i,j-1);
                    Ds(i,j-1) = tran2;  
                    tran3 = S(i,j);
                    S(i,j) = S(i,j-1);
                    S(i,j-1) = tran3; 
                    
               case k==5
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i+1,j+1);
                    Sm(i+1,j+1) = tran;    
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i+1,j+1);
                    Ds(i+1,j+1) = tran2; 
                    tran3 = S(i,j);
                    S(i,j) = S(i+1,j+1);
                    S(i+1,j+1) = tran3;                   
               case k==6
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i+1,j-1);
                    Sm(i+1,j-1) = tran;  
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i+1,j-1);
                    Ds(i+1,j-1) = tran2;   
                    tran3 = S(i,j);
                    S(i,j) = S(i+1,j-1);
                    S(i+1,j-1) = tran3;                     
               case k==7
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i-1,j+1);
                    Sm(i-1,j+1) = tran; 
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i-1,j+1);
                    Ds(i-1,j+1) = tran2;    
                    tran3 = S(i,j);
                    S(i,j) = S(i-1,j+1);
                    S(i-1,j+1) = tran3;                     
               case k==8
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i-1,j-1);
                    Sm(i-1,j-1) = tran; 
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i-1,j-1);
                    Ds(i-1,j-1) = tran2;  
                    tran3 = S(i,j);
                    S(i,j) = S(i-1,j-1);
                    S(i-1,j-1) = tran3; 
                otherwise
                    Sm(i,j) = Sm(i,j);
                    S(i,j) = S(i,j);
                    Ds(i,j) = Ds(i,j);              
            end           
        end
        if Sm(i,j)==4 %第四个出口
            %距离
            D1 = sqrt(abs(i+1-X(4)).^2+abs(j-Y(4)).^2);%右方向
            D2 = sqrt(abs(i-1-X(4)).^2+abs(j-Y(4)).^2);%左方
            D3 = sqrt(abs(i-X(4)).^2+abs(j+1-Y(4)).^2);%上方
            D4 = sqrt(abs(i-X(4)).^2+abs(j-1-Y(4)).^2);%下方
            D5 = sqrt(abs(i+1-X(4)).^2+abs(j+1-Y(4)).^2);%右上
            D6 = sqrt(abs(i+1-X(4)).^2+abs(j-1-Y(4)).^2);%右下
            D7 = sqrt(abs(i-1-X(4)).^2+abs(j+1-Y(4)).^2);%左上
            D8 = sqrt(abs(i-1-X(4)).^2+abs(j-1-Y(4)).^2); %左下       
            %判断周围情况
            if Skm(i+1,j)~=2
                D1 = inf;  
            end
            if Skm(i-1,j)~=2
                D2 = inf;   
            end
            if Skm(i,j+1)~=2
                D3 = inf;  
            end
            if Skm(i,j-1)~=2
                D4 = inf;  
            end
            if Skm(i+1,j+1)~=2
                D5 = inf;   
            end
            if Skm(i+1,j-1)~=2
                D6 = inf;     
            end
            if Skm(i-1,j+1)~=2
                D7 = inf;   
            end
            if Skm(i-1,j-1)~=2
                D8 = inf;
            end      
            
            Dt = [D1,D2,D3,D4,D5,D6,D7,D8];
            Dmin1 = min(Dt);
            for k = 1:8
                if Dt(k)==Dmin1
                    break
                end
                if Dt(k)==inf
                    k = 9;
                end
            end           
            switch k
                case k==1
                    tran = Sm(i,j);%状态矩阵改变
                    Sm(i,j) = Sm(i+1,j);
                    Sm(i+1,j) = tran;  
                    
                    tran2 = Ds(i,j);%距离矩阵改变
                    Ds(i,j) = Ds(i+1,j);
                    Ds(i+1,j) = tran2; 
                    
                    tran3 = S(i,j);%状态矩阵改变
                    S(i,j) = S(i+1,j);
                    S(i+1,j) = tran3; 
               case k==2
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i-1,j);
                    Sm(i-1,j) = tran; 
                    
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i-1,j);
                    Ds(i-1,j) = tran2;
                    tran3 = S(i,j);
                    S(i,j) = S(i-1,j);
                    S(i-1,j) = tran3;                 
               case k==3
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i,j+1);
                    Sm(i,j+1) = tran;
                    
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i,j+1);
                    Ds(i,j+1) = tran2;
                    
                    tran3 = S(i,j);
                    S(i,j) = S(i,j+1);
                    S(i,j+1) = tran3;
                    
               case k==4
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i,j-1);
                    Sm(i,j-1) = tran; 
                    
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i,j-1);
                    Ds(i,j-1) = tran2;  
                    tran3 = S(i,j);
                    S(i,j) = S(i,j-1);
                    S(i,j-1) = tran3; 
                    
               case k==5
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i+1,j+1);
                    Sm(i+1,j+1) = tran;    
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i+1,j+1);
                    Ds(i+1,j+1) = tran2; 
                    tran3 = S(i,j);
                    S(i,j) = S(i+1,j+1);
                    S(i+1,j+1) = tran3;                   
               case k==6
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i+1,j-1);
                    Sm(i+1,j-1) = tran;  
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i+1,j-1);
                    Ds(i+1,j-1) = tran2;   
                    tran3 = S(i,j);
                    S(i,j) = S(i+1,j-1);
                    S(i+1,j-1) = tran3;                     
               case k==7
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i-1,j+1);
                    Sm(i-1,j+1) = tran; 
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i-1,j+1);
                    Ds(i-1,j+1) = tran2;    
                    tran3 = S(i,j);
                    S(i,j) = S(i-1,j+1);
                    S(i-1,j+1) = tran3;                     
               case k==8
                    tran = Sm(i,j);
                    Sm(i,j) = Sm(i-1,j-1);
                    Sm(i-1,j-1) = tran; 
                    tran2 = Ds(i,j);
                    Ds(i,j) = Ds(i-1,j-1);
                    Ds(i-1,j-1) = tran2;  
                    tran3 = S(i,j);
                    S(i,j) = S(i-1,j-1);
                    S(i-1,j-1) = tran3; 
                otherwise
                    Sm(i,j) = Sm(i,j);
                    S(i,j) = S(i,j);
                    Ds(i,j) = Ds(i,j);              
            end           
        end
    n = n+1    
    set(Ii,'CData',S);%更新显示的数据
    pause(0.2);%暂停一下显示出动画效果 
    end
end
if sum(sum(Ds))<=10000
    break;
end
end