matlab元胞元素_基于元胞自动机的二维对向行人仿真

前几天在知乎看到一个用元胞自动机做的行人仿真模型，是用python写的，自己用matlab复现了一遍，也看了那篇论文，感觉行人仿真还挺有意思的，以及仿真页面可视化的思路还挺好用且对于元胞自动机具有普适性，现在把代码分享出来。

知乎视频​www.zhihu.com

这个模型是引入决策空间，主动慢化，优势通道的异步更新模型，相比于以前的元胞自动机模型是能够在中，低密度情况下仿真出现死锁现象。具体内容可以看这篇文章介绍，和文末的那篇文章，我说一下我的理解。

该模型的可视化方法是将将平面离散化，一个元胞占据第一象限任意整数向上下左右移动0.5形成四条边组成的矩形，每个矩形中的颜色可以取{白色，共色，蓝色} ，分别代表{无人，有初始位置在左边，在右边的行人}。

定义概念决策空间：行人所在元胞向前向后各若干元胞这个范围内其他行人的数量和前进方向会对行人决策产生影响，这个范围即为该行人的决策空间。

定义概念优势通道：对一个行人而言，其决策空间内，如果一个通道上的同向行走人数 > 对向行走人数，或在其前方决策空间内没有行人，都认为这条通道为其优势通道。

更新规则如下：

在决策空间内有优势通道时，如果所在通道为优势通道时前进，否则换道至可行的最大优势通道，若无可行优势通道，原地不动

若决策空间内无优势通道时，以一定概率原地不动，若移动，优先前进，不能前进则优先左转，不行右转（左右是面向前进方向而言）

贴一下代码

clc;
clear;
close all;
dbstop if error;
% 参数设置
w = 20; %宽带
l = 40; %长度
l_wait = 10; %等待区域宽度
p_size = 60; %行人可视化尺寸
pause_time = 0.1; %图像输出间隔
p_move = 0.7; %无优势通道时移动概率
r_adv = [4,2]; %决策空间范围，宽为3行，矩阵两个元素分别为距其前后列数
num_l = 90; %左边等待区域人数
num_r = 90; %右边等待区域人数
%%  GUI
run = 0;  %值为0时暂停，值为1时开始
stop = 0;  %值为1时退出模型
hf=figure('doublebuffer','on');
uicontrol('style','pushbutton','string','run','callback','run = 1;','position',[500,50,50,25]);
uicontrol('style','pushbutton','string','pause','callback','run = 0;','position',[600,50,50,25]);
uicontrol('style','pushbutton','string','stop','callback','stop = 1;','position',[700,50,50,25]);
%% 仿真初始化
[space,ped_l,ped_r] = initialize(w,l,l_wait,num_l,num_r);
ped_plot(w,l,l_wait,ped_l,ped_r,p_size);
% 开始仿真
while ~isempty(ped_l.n) || ~isempty(ped_r.n) 
    if stop == 0
        while run == 0
            pause(0.01);  %暂停以使程序能够接受到用户输入
            if stop == 1
                close all;
                break;
            end
        end
        % 行人位置更新
        [space,ped_l,ped_r]=update(w,l,ped_l,ped_r,num_l,num_r,space,p_move,r_adv);
        % 图像更新
        cla;
        ped_plot(w,l,l_wait,ped_l,ped_r,p_size);
        pause(pause_time);
    else
        close all;
        break;
    end
end
close all;

% 行人位置图像生成函数
function ped_plot(w,l,l_wait,ped_l,ped_r,p_size)
hold on;
% 生成纵线
for id = 1 : l
    plot([id+0.5 , id+0.5] , [0.5 , w+0.5] , 'k');
end
% 生成横线
for id = 1 : w
    plot([0.5 , l+0.5] , [id+0.5 , id+0.5] , 'k');
end
% 生成等待区界线
plot([l_wait+0.5 , l_wait+0.5] , [0.5 , w+0.5] , 'c' , 'linewidth' , 3);
plot([l-l_wait+0.5 , l-l_wait+0.5] , [0.5 , w+0.5] , 'g' , 'linewidth' , 3);
% 汇出行人位置
plot(ped_l.n , w-ped_l.m+1 , '.r' , 'MarkerSize' , p_size); %矩阵储存顺序和y轴正反向相反
plot(ped_r.n , w-ped_r.m+1 , '.b' , 'MarkerSize' , p_size);
% 隐藏坐标轴
xticks([]);
yticks([]);
xlim([0.5 , l+0.5]);
ylim([0.5 , w+0.5]);
end

% 初始化函数
function [space,ped_l,ped_r] = initialize(w,l,l_wait,num_l,num_r)
% 随机生成行人初始位置
space = zeros(w , l); %空间平面矩阵
mid = 1 + fix(rand(1)*(l_wait*w - 1));
for id = 1 : num_l+num_r
    if id <= num_l
        while ismember(mid , find(space == 1))
            mid = 1 + fix(rand(1)*(l_wait*w - 1)); %左边等待区随机生成行人
        end
    else
        while ismember(mid , find(space == 1))
            mid = 1+(l-l_wait)*w + fix(rand(1)*(l*w - (l-l_wait)*w)); %右边等待区随机生成行人
        end
    end
    space(mid) = 1;
end
% 生成行人结构体
[y1 , x1]=find(space(1:w , 1:l_wait)); %y = m,表示在m行
[y2 , x2]=find(space(1:w , l-l_wait+1:l)); %x = n,表示在n列
x2 = x2 + l-l_wait; %切片使用find函数要补上被切去的列数
ped_l = struct('m' , y1' , 'n' , x1' , 'p' , y1*w+x1); %初始在左边的行人信息结构体
ped_r = struct('m' , y2' , 'n' , x2' , 'p' , y2*w+x2); %行人所在行，列，空间位置
end

% 行人位置更新函数
function [space,ped_l,ped_r]=update(w,l,ped_l,ped_r,num_l,num_r,space,p_move,r_adv)
id_in_l = zeros(num_l); %储存初始在左边未完成疏散的行人id
id_count = 1; %计数器
% 左边行人更新
for id = 1:length(ped_l.n)
    space(ped_l.m(id) , ped_l.n(id))= 0;
    % 计算下一时刻位置
    if ped_l.n(id) ~= l %下一时刻完成疏散行人将被移除空间不参与更新
        % 优势通道情况获取
        [adv] = advantage(r_adv,w,l,ped_l,id,space,1);
        % 下一时刻位置更新
        if ~isempty(find(adv(1,:) .* adv(3,:))) %有可行的优势通道
            if adv(1,2) * adv(3,2) %优先前进，行人偏好原通道
                ped_l.n(id) = ped_l.n(id) + 1;
            else
                [~,m_mid] = max(adv(1,:) .* adv(3,:));
                ped_l.m(id) = ped_l.m(id) + m_mid - 2; %向更多同行行人通道换道
            end
        else
            if rand(1) < p_move %无可行优势通道概率移动
                % 若移动，优先前进，不能前进则优先左转，不行右转（左右是面向前进方向而言）
                if adv(3,2)
                    ped_l.n(id) = ped_l.n(id) + 1;
                elseif adv(3,1)
                    ped_l.m(id) = ped_l.m(id) - 1;
                elseif adv(3,3)
                    ped_l.m(id) = ped_l.m(id) + 1;
                end
            end
        end
        % 位置矩阵更新
        space(ped_l.m(id) , ped_l.n(id))= 1;
        % 未完成疏散行人标记
        id_in_l(id_count) = id;
        id_count = id_count + 1;
    end
    % 空间位置更新
    ped_l.p = ped_l.n*w + ped_l.m;
end
% 右边行人更新
id_in_r = zeros(num_r); %储存初始在右边未完成疏散的行人id
id_count = 1; %计数器
for id = 1:length(ped_r.n)
    space(ped_r.m(id) , ped_r.n(id))= 0;
    % 计算下一时刻位置
    if ped_r.n(id) ~= 1 %下一时刻完成疏散行人将被移除空间
        % 优势通道情况获取
        [adv] = advantage(r_adv,w,l,ped_r,id,space,0);
        if ~isempty(find(adv(1,:) .* adv(3,:))) %有可行的优势通道
            if adv(1,2) * adv(3,2) %优先留在原通道
                ped_r.n(id) = ped_r.n(id) - 1;
            else
                [~,m_mid] = max(adv(1,:) .* adv(3,:));
                ped_r.m(id) = ped_r.m(id) + m_mid - 2; %向更多同行行人通道换道
            end
        else
            if rand(1) < p_move %无可行优势通道概率移动
                % 若移动，优先前进，不能前进则优先左转，不行右转（左右是面向前进方向而言）
                if adv(3,2)
                    ped_r.n(id) = ped_r.n(id) - 1;
                elseif adv(3,3)
                    ped_r.m(id) = ped_r.m(id) + 1;
                elseif adv(3,1)
                    ped_r.m(id) = ped_r.m(id) - 1;
                end
            end
        end
        % 位置矩阵更新
        space(ped_r.m(id) , ped_r.n(id))= 1;
        % 未完成疏散行人标记
        id_in_r(id_count) = id;
        id_count = id_count + 1;
    end
    % 空间位置更新
    ped_r.p = ped_r.n*w + ped_r.m;
end
% 剔除完成疏散行人
id_in_l = id_in_l(id_in_l ~= 0);
id_in_r = id_in_r(id_in_r ~= 0);
ped_l.m =ped_l.m(id_in_l);
ped_l.n =ped_l.n(id_in_l);
ped_r.m =ped_r.m(id_in_r);
ped_r.n =ped_r.n(id_in_r);
end

% 优势通道判断及三个移动方向可行性判断函数
function [adv] = advantage(r_adv,w,l,ped,id,space,flag)
% flag用于区分左边行人和右边行人，1 左 0 右
adv = zeros(3,3); %储存决策空间内三条通道是否为优势通道的逻辑值数组，该通道同向行人数，该方向可行性
adv_count = 0; %同向车辆计数器
for h = max(1 , ped.m(id)-1) : min(w , ped.m(id)+1)
    % 依次判断决策范围内每条通道是否为优势通道
    for x = max(1 , ped.n(id)-r_adv(2)) : min(l , ped.n(id)+r_adv(1))
        % 数出该行决策空间内同向行人数（包括自己,因为不与其他两通道比较无所谓）
        if ismember(x*w+h , ped.p)
            adv_count = adv_count + 1;
        end
    end
    % 计算对向行人数
    disadv = length(find(space(h,max(1,ped.n(id)-r_adv(2)):min(l,ped.n(id)+r_adv(1))) == 1)) - adv_count; %方向行人数
    if adv_count > disadv || isempty(space(h , min(l , ped.n(id)+r_adv(1)))) %决策空间内同向行人占优或前方决策空间无人
        if h == ped.m(id)
            adv(1 , 2) = 1;
            adv(2 , 2) = adv_count;
        elseif h < ped.m(id)
            adv(1 , 1) = 1;
            adv(2 , 1) = adv_count;
        else
            adv(1 , 3) = 1;
            adv(2 , 3) = adv_count;
        end
    end
    % 三个移动方向可行性判断
    if h == ped.m(id)
        if flag
            if space(h , ped.n(id) + 1) == 0 %初始在左边行人可以前进
                adv(3 , 2) = 1;
            end
        else
            if space(h , ped.n(id) - 1) == 0 %初始在右边行人可以前进
                adv(3 , 2) = 1;
            end
        end
    elseif h < ped.m(id)
        if space(h, ped.n(id)) == 0 %可以向左（对于左边行人）
            adv(3 , 1) = 1;
        end
    else
        if space(h, ped.n(id)) == 0 %可以向右（对于左边行人）
            adv(3 , 3) = 1;
        end
    end
end
end

运行模型我发现两个现象

知乎视频​www.zhihu.com

知乎视频​www.zhihu.com

1. 如果模型中允许行人可以后退，那么对向行人在对峙时就没有礼让，会出现死锁。

2. 会出现边界整个平面只有一个可通行行，这行走完之后，可通行行向两边扩展，使得缺口从此向两边扩展的现象。即，两个方向行人优势通道交叉排列现象，这时候只要创造条件使得其中任意一行行人走完毕就可以使得所有人顺序通过。

只要将行人换为车辆，通过不同颜色来区分不同类型的车辆就能很直观地区分开来了。