元胞自动机

元胞自动机是复杂系统研究的一个典型方法，特别适合用于空间复杂系统的时空动态模拟研究。不同于一般的动力学的模型，元胞自动机不是由严格定义的物理方程或函数组成，而是用一系列模型构造的规则构成，或者说是方法框架

五个最基本条件

• 包含元胞$cell$的规则和网格$Lattice$ $Grid$。空间被离散成网格，元胞分布在网格中的各个格子里.
• 离散的时间步长。元胞状态的改变是定义在两个时间步之间的，同一个时间步内，元胞的状态的固定的。
• 每个元胞的状态数是有限的。也就是说元胞只能在固定几种状态中切换。
• 所有元胞状态的更新都遵循同样的规则。不存在享受特殊待遇的元胞
• 元胞自动机的规则定义在局部。也就是说元胞状态只受其周围元胞状态的影响

满足上述五个条件的都算是元胞自动机

要素

• 元胞：作为规则网格中的一个格子，元胞可以有很多形状。元胞根据它周围元胞的状态来改变自身的状态，状态可以是 { 0 , 1 } \{0,1\} {0,1}的二进制形式，或是 { s 0 , s 1 , s 2 , . . . s n } \{s_0,s_1,s_2,...s_n\} {s0​,s1​,s2​,...sn​}整数形式的离散集，严格意义上，元胞自动机的元胞只能有一个状态变量，但在实际应用中，往往将其进行了拓展，每个元胞可以有多个状态

• 网格：所有元胞都是放在网格中的，从纬度上分，网格主要有一维，二维，三维及更高维。最简单的情况就是一维网格，一维网格就是一系列的元胞排成一条线，在模拟单条车道的交通流时就是用的一维网格。下面是常见的二维网格：
  

• 邻居：以上的元胞以及元胞空间都只是系统的静态成分，为将”动态“引入系统，必须加入演化规则。在元胞自动机中，这些规则是定义在空间局部范围内的, 即一个元胞下一时刻的状态决定于本身状态和它的邻居元胞的状态。因而，在指定规划前，必须先定义邻居规则，即明确哪些元胞属于该元胞的邻居，在一维网格中，通常以半径大小确定邻居，例如距离一个半径，即前后。二维元胞自动机的邻居定义比较复杂，在二维的方形网格中，常用的邻居有如下三种：
  

• 边界条件：

  💡 在理论上, 元胞空间通常是在各维向上是无限延展的, 这有利于在理论上的推理和研究。但是在实际应用过程中，我们无法在计算机上实现这一理想条件，因此实际模拟中就会有边界。

  在定义元胞的邻居时，会出现一个问题：边界上的元胞邻居不足。元胞自动机对每个元胞施加同样的规则，因此需要设置边界条件，使得边界上的元胞和其他元胞具有相同的邻居数。下图表示了4种常见的边界条件：

• 规则：根据元胞当前状态及其邻居状况确定下一时刻该元胞状态的动力学函数，简单来讲就是一个状态转移函数。我们称该函数为元胞自动机的局部映射或局部规则。规则是支配整个元胞自动机的动力学行为的。一般规则都是定义在局部。通过局部间元胞的相互作用而引起全局变化。

比如说流行病的传播, 只有相互靠的近的人才能进行直接的传播, 这就是局部作用, 但当流行病传播一定时间时就会引起大范围的变化, 这就是全局变化. 而这个全局变化是由一个个局部作用引起的. 元胞自动机的规则一般为总和型. 就是说某个元胞下时刻的状态只决定于它所有邻居当前状态及自身的当前状态. 一个简单的例子是足球场观众台上的人浪, 当你最靠近的人即邻居站起来时, 你才能站, 当他们坐下后你才能坐下, 这就是规则。

CA模拟森林火灾Code1

clear,clc;
n = 200;    % 森林的尺寸：n x n
pltg = 5e-6;    % 闪电的概率
pgrw = 1e-2;    % 生长的概率
NW = [n 1:n-1]; % 表示北面邻居和西面邻居，第一行的北面是第n行，第一列的西面的第n列
SE = [2:n,1];   % 表示南面邻居和东面邻居，第一行的南面是第2行，第一列的东面是第2列
veg = zeros(n); % veg={0表示空，1表示火，2表示树}
imh = image(cat(3,(veg==1),(veg==2),zeros(n)));

for i=1:3000
    % 周围四个邻居中状态为火的数量 生成的num是n x n的logical数组
    num =               (veg(NW,:)==1) + ...
            (veg(:,NW)==1)      +       (veg(:,SE)==1) + ...
                        (veg(SE,:)==1);
    % 状态更新方程（规则）
    veg = 2*( (veg==2) | (veg==0 & rand(n)<pgrw)) - ... % 表示下一状态是树的两种状态-veg=2本身就是树和veg=0是空地且符合生长条件
            ( (veg==2) & (num >0 | rand(n)<pltg));      % 表示下一状态是火的状态-veg=2要有树，且要么周围有火num>0要么有闪电

    set(imh,'cdata',cat(3,(veg==1),(veg==2),zeros(n)))
    drawnow
end

Code2

% 元胞自动机模拟森林火灾
% 规则1：正在燃烧的树下一刻会变成空格
% 规则2：如果绿树格位的最近邻居中有一个树在燃烧，则下一刻也会变成燃烧的树
% 规则3：在空格处，有概率p的可能长出一棵树
% 规则4：不失一般性，如果考虑南风对森林火灾的影响，对于燃烧的树有w的概率会将上方距离有两格的树引燃

%% 准备工作
close all;
clc,clear;
figure;
p=0.01; % 概率p，用来表示树生成的速度
f=6e-5; % 概率f，用来表示闪电击中的概率
w=0.8;
axes;
% 刷新rand的state，避免第二次运行是接着上一次运行的结果
rand('state',0);
% 设置的目的是为了防止在不断循环画动画的时候会产生闪烁的现象
set(gcf,'DoubleBuffer','on');

%% 初始化
% 300x300的矩阵森林中，1表示树木，2表示火灾，0表示空
S=ones(300);    % 初始化森林，全为树木
S(151,151)=2;   % 假设森林的中心出现火灾
% 用302x302的Sk矩阵储存森林 给周围设定边界
Sk=zeros(302);  
Sk(2:301,2:301)=S;
% 画布设置用RGB色彩空间表示
C=zeros(302,302,3); % 构造一个三维张量，用来表示RGB色彩空间
R=zeros(300);       % 初始化R矩阵
G=zeros(300);       % 初始化G矩阵
R(S==2)=1;          % S矩阵上某个点的位置若为2，则给R矩阵对应位置赋值1
G(S==1)=1;          % S矩阵上某个点的位置若为1，则給G矩阵对应位置赋值1
C(2:301,2:301,1)=R; % 将C的第一层设置为R
C(2:301,2:301,2)=G; % 将C的第二层设置为G
Ci=imshow(C);
% 另外一种表示方法
% imh = imshow(cat(3,Sk==2,Sk==1,zeros(302));

%% CA模拟过程
ti=0; % 初始化时间
tp=title(['T = ',num2str(ti)]);
while 1
    ti=ti+1;
    St=Sk;  % St表示t时刻的森林情况，Sk表示随机出来的森林情况
    % 规则1
    St(Sk==2)=0;    % ti-1时刻正在燃烧的树在ti时刻烧成空
    Su=zeros(302);

    Sf = Sk;
    Sf(Sf<=1)=0;% 将森林中的空白点和树都去掉
    Sf=Sf/2;    % 着火点变成1
    % 统计周围八个邻居的状态,Su>0表示周围有树木燃烧
    Su(2:301,2:301)=Sf(1:300,1:300)+Sf(1:300,2:301)+Sf(1:300,3:302)+...
        Sf(2:301,1:300)+Sf(2:301,3:302)+...
        Sf(3:302,1:300)+Sf(3:302,2:301)+Sf(3:302,3:302);
    % 规则2
    St(Su>0 & Sk==1) = 2; % 经过上述统计，附近有数目燃烧 & 自身是树木的状态修改为燃烧

    % 规则4
    for i=2:301
        for j=2:301
            r0=rand(1);
            if St(i,j)==2 && St(i-1,j)==1
                St(i-1,j)=round(2*(r0<=w)+(r0>w));
                St(i-2,j)=round(2*(r0<=w)+(r0>w));
            end
        end
    end
    
    % 规则3
    Se=Sk(2:301,2:301); % Se初始化为最初的森林
    Se(Se==0)=4;    % 空白的地方赋值4
    Se(Se<=2)=0;    % 有树1，有火灾2的地方赋值0
    Se(Se==4)=1;    % 空白的地方再改回1
    % 经过上述操作，最初森林中，空白的地方值为1，不空白的地方值为0
    St(2:301,2:301)=St(2:301,2:301) + Se.*(rand(300)<p);
    
    % 把当前时刻的森林传递传递给Sk
    Sk = St;

    %% 更新Image
    R = zeros(302);
    G = zeros(302);
    R(Sk==2)=1;
    G(Sk==1)=1;
    C(:,:,1)=R;
    C(:,:,2)=G;
    set(Ci,'CData',C);
    set(tp,'string',['T = ',num2str(ti)]);
    pause(0.2);
end