Matlab 蒙特卡洛模拟

发布于 2020-07-18  285 次阅读


Matlab 蒙特卡洛模拟

视屏 1.1

实验原理建模

  蒙特卡洛指的是使用计算机方法来模拟概率论之中的随机问题,在模拟概率论之中的问题时,只需要保持随机性即可,故一般采用假随机种子模拟随机数的产生。

  在人群疏散模型中,采用棋盘网格地图,棋子表示人员,所有棋子都需要向出口转移,统计每一轮最后一个棋子离开棋盘的时间。标准棋盘图像如下,非边缘白色为空地,边缘白色为出口,黑色为障碍物。

DefaultMap.png

建模过程

  • FlashMap函数

  这是一个构造一个新地图的函数,输入参数为地图长宽,人数,与除围墙以外的障碍物个数,输出为地图。

  使用二维数组结构模拟棋盘地图,数字0代表空地,数字1代表障碍物,负数代表出口,大于1的正整数代表人,并且移动顺序是从数字较小的人先开始移动。

  • MV函数

  这是一个移动人员的函数,输入参数为当前地图与当前在棋盘中人的列表,输出下一轮的地图与下一轮在棋盘之中人的列表。

  计算规律为:遍历当前存活列表,先检查节点邻域之中是否存在出口,若有则将当前位置至零,并且在存活列表中删除,否则计算每个节点最近的出口,在通过IF-ELSE语句判断下一步如何走。

  • 主函数

  这是负责调用每一次刷新地图与统计判断逃离棋盘时间的函数。

流程图

实验结果


注:实验数据采用多次试验的平均值,固定人数为8人,网格为10X10.
注:图像为出口—时间(平均等待时间)图像


注:出口宽度设置为3,棋盘大小设置为10X10.

代码部分

function T = Main(x,y,S,W,K)
Map = FlashMap(x,y,S,W,K);
T = 0;
live_list = 2:S+1;
while ~isempty(live_list)
    [Map,live_list] = Mv(Map,live_list);
    T  = T +1;
    %painting(Map,T,S,W);
    end
end


function Map = FlashMap(x,y,S,W,K) 
Space = zeros(x-2,y-2);
Space(randperm((x-2)*(y-2),S)) = randperm(S) + 1; 
B = zeros((x-2)*(y-2)-S,1);
B(randperm((x-2)*(y-2)-S,K)) = 1;
Space(Space==0) = B;
C = ones(x-2,1);
D = ones(1,y);
Map = [D;C Space C; D];
Map(2:W+1,1) = -1:-1:-W;
end


function [NextMap,live] = Mv(Map,live_list)
[x, ~] = size(Map);
for i = live_list
    index = find(Map == i);
    [a, b] = find(Map == i);
    Field = Map(a-1:a+1,b-1:b+1);
    if ~isempty(find(Field < 0, 1))
        Map(index) = 0;
        live_list(find(live_list == i, 1)) = [];
        continue
    end
    length = (find(Map<0) - a).^2 + (b - 1)^2;
    [~, l] = min(length);
    [ox, oy] = find(Map == -l);
    if (oy < b) && (Map(a,b-1) == 0)
        Map(index) = 0;
        Map(a,b-1) = i;
        continue
    elseif (oy > b) && (Map(a,b+1) == 0)
        Map(index) = 0;
        Map(a,b+1) = i;
        continue
    end
     if (ox < a) && (Map(a-1,b) == 0)
         if a>2 && b >2    
            if (Map(a,b-2) == 1) && (Map(a-1,b-1) == 1)
                Map(index) = 0;
                Map(a-1,b) = i;
                continue
            end
         end
                 Map(index) = 0;
                 Map(a-1,b) = i;
     elseif (ox > a) && (Map(a+1,b) == 0)
         if a>2 && b >2    
            if (Map(a,b-2) == 1) && (Map(a+1,b-1) == 1)
                Map(index) = 0;
                Map(a+1,b) = i;
                continue
            end
         end
         Map(index) = 0;
         Map(a+1,b) = i;
     else
         Map(index) = 0;
         if Map(a,b-1) == 0
             Map(a,b-1) = i;
         elseif Map(a+1,b) == 0
             Map(a+1,b) = i;
         elseif Map(a,b+1) == 0
             Map(a,b+1) = i;
         elseif Map(a-1,b) ==0 
             Map(a-1,b) = i;
         else 
             Map(index) = i;
         end         
    end
end
NextMap = Map;
live = live_list;
end



function painting(Map,T,S,W)
c = [-W S+1];
imagesc(Map,c);
print(gcf,'-dpng',[num2str(T),'.png']);
close
end

结论部分

  根据图 1-1 与 1-2,发现随着W出口宽度的增加,平均等待时间与总等待时间都是呈现下降趋势,在图 1-1 之中,做线性回归方程,发现 $R^2 = 0.9394$这表明平均等待时间与出口的线性相关程度特别大。但在图 1-2 中,发现总等待时间(取整后)线性相关程度降低,但整体上呈现负相关性。

  在图 2-1 与 2-2中,随着S人数的增加,平均等待时间呈现下降趋势,总体等待时间呈现上升趋势,在图 2-1 之中,做回归方程,得出y = 6.223x^{-0.7356} ,R^2 = 0.9978 , 这表明平均等待时间与人数的函数是一个递减幅度减小的减函数,随着人数的增加,平均等待时间递减幅度会越来越小,且最后趋近于0,当然从直观上,也可以发现,最后也会趋紧于0。

  在总体上,当人数一定时,出口越多,平均等待时间就越低,总等待时间越小;当出口宽度一定时,人数约多,平均等待时间就越低,但总等待时间会越长。

拓展部分

  笔者认为MV函数一般来说,解决普通有障碍物的情况时,可以找到最短路径,但是就某些空地上有障碍物的特殊情况,棋子会出现“摇摆不定,左右为难”的情况,无法找到正确的路径出口。

  因此笔者改进了MV函数方法,采用Search函数进行先搜索无障碍路径,再优化路径长度的方法。

  笔者在这里提供两种计算最近出口的方法,两种方法均是比较长度的平,由于同样为增函数,因此长度的平方大小比较反应了长度的大小比较。

  其中一种十分普通,采用非MATLAB独有方式计算长度,第二种代码如下:其中Father为棋子所在节点位置,譬如[3, 4],endpoint是考虑到邻域问题,将所有终点设置在其领域上,譬如出口节点为[2, 1] 则其领域节点为[1:3,2]。

    endpoint = [[1:W+1]' , ones(W+1,1)*2];
    fatherpoint = repmat(Father,W+1,1);
    length = (endpoint(:,2) - fatherpoint(:,2)).^2 + (endpoint(:,1) - fatherpoint(:,1)).^2;
    [~, l] = min(length);
    [a, b] = endpoint(l,:);

树状图

  Search函数每一次返回下一个节点,搜索方式为,30%概论随机搜索,70%向节点最近出口搜索。整体代码如下:

function next = search(Map,Father,W)
if rand()<0.35
    i = rand()
    if i<0.25
        next = Father + [0 1];
    elseif i<0.5
        next = Father + [1 0];
    elseif i<0.75
        next = Father + [0 -1];
    else
        next = Father + [-1 0];
    end 
else
    endpoint = [[1:W+1]' , ones(W+1,1)*2];
    fatherpoint = repmat(Father,W+1,1);
    length = (endpoint(:,2) - fatherpoint(:,2)).^2 + (endpoint(:,1) - fatherpoint(:,1)).^2;
    [~, l] = min(length);
    [a, b] = endpoint(l,:);
    if a < Father(1) && b< Father(2)
        if rand()<0.5
            next = Father + [-1 0];
        else
            next = Father + [0 -1];
        end
    end
    if a == Father(1) && b < Father(2)
        next = Father + [0 -1];
    end

end
end

  在主函数中改动部分代码,将live_list列表改变为 path_list,使用while循环到搜索到出口,则path_list里面节点最长的路径为总等待时间T。

  当然这样算法方式是一定能够找到出口,但不能保证路径一定是最短,除非将树状图完全拓扑,找到最短的分枝。使用完全拓扑树状图时,如果棋盘扩大,计算量是指数级增长。

  笔者曾经使用遗传算法独立解决过TSP(旅人)问题「1」,因此笔者认为是否可以引入遗传算法优化路径问题呢。

  思路如下:
  首先初始化节点路径的种群,DNA为节点路径,tRNA为计算随机中途节点和DNA末端节点离出口距离的方法,通过tRNA将DNA翻译为PRO,PRO为离出口的距离,自然选择函数Select 中,20%概率随机选择可以交配的下一代,80%概率选择离出口较近的DNA进行交配。交配函数中,8%会发生单个基因突变,2%概率发生基因段交换,返回下一代,不断迭代,可以将路径长度降低到一个可观的程度。

  笔者在TSP问题中,采用标准数据集,将总距离降低到3000左右。


你知道雪为什么是白色的吗?因为她忘记了原来的颜色