元胞自动机——应用于森林火灾和传染病场景
最近接触了元胞自动机模型,做了一些资料搜查,并进行学习,推荐这篇文章澳洲变燠洲,考拉成烤拉!澳大利亚山火为什么难以控制?
以下对所学进行记录。
森林火灾元胞自动机原理
在元胞自动机模型中,空间被离散成网格,每一个网格被称为元胞。森林火灾元胞有三种状态:树,火(正在燃烧的树)和空(空地)状态。元胞下一时刻状态的更新规则如下:
树变火:一棵树,其上下左右若有一个状态为火,下一刻就会变成火。或者一棵树遇上闪电,下一刻就会变成火。由于遇上闪电着火的概率Plight很小。
火变空:火在下一时刻会变成空。
空变树:空地下一时刻会以一个很小的概率Pgrowth长出新树。
改进模型会考虑树的对角位置有没有着火。或者会考虑风向(比如吹西风(火从东吹向西),火的西边着火的机率会变大(顺风),火的东边着火的几率变小(逆风)),这里盗个图:
图a是基础元胞自动机,图b是考虑对角的元胞自动机,图c是吹西风的元胞自动机。
本人这里对基础元胞自动机,考虑对角情况的元胞自动机,考虑吹西风的元胞自动机三种模型进行仿真,仿真结果如下:
代码如下:
m = 300;
n = 300; % 表示森林的矩阵行列 m x n
Plight = 5e-6; % 闪电概率
Pgrowth = 1e-2; % 生长概率
% 邻居方位 d 和点燃概率 p
d1 = {[1,0], [0,1], [-1,0], [0,-1]};%最下一行往最上,最右一列往最左,最前一列往最下,最左一列往最右,circshift(S,d{j})
p1 = [ 1, 1, 1, 1];
% % 改进元胞自动机
d2 = {[1,0], [0,1], [-1,0], [0,-1], [1,1], [-1,1], [-1,-1], [1,-1]};%考虑对角
p2 = [ones(1,4), ones(1,4)*(sqrt(1/2)-1/2)];
% % 考虑风的情况,西风,因此西边方向概率为0.3,东方向概率为1,增加对角(0.3),东二格子(0.8)
d3 = {[1,0], [0,1], [-1,0], [0,-1], [1,1], [-1,1], [-1,-1], [1,-1], [0,-2]};
p3 = [ 0.80, 0.30, 0.80, 1.00, 0.12, 0.12, 0.30, 0.30, 0.8];
% 空=0, 火=1, 树=2
E = 0;
F = 1;
T = 2;
S = T*(ones(m,n));%全树
S(m/2,n/2)=1;
imh = image(cat(3, S==F, S==T, zeros(m,n)));
axis image;
%plot box fits tightly around the data
for t = 1:3000
% 计算邻居中能传播着火的个数
sum = zeros(size(S));
note=' 吹西风';
p=p3;
d=d3;
for j = 1:length(d)
sum = sum + p(j) * (circshift(S,d{j})==F);
end
isE = (S==E);
isF = (S==F);
isT = (S==T); % 找出三种不同的状态
ignite = rand(m,n)<sum | (rand(m,n)<Plight); % 着火条件
% 规则 1: 着火F
Rule1 = F*(isT & ignite);
% 规则 2: 烧尽E
Rule2 = F*isF - F*isF;%这里的结果是0矩阵,除去新着火,新生及原来没烧的树木,剩下的就是烧尽的空(0)
% 规则 3: 新生及原来没烧的树木T
Rule3 = T*(isE & rand(m,n)<Pgrowth)+T*(isT & ~ignite) ;
S = Rule1 + Rule2 + Rule3;
Fnum=length(find(S==F));
Tnum=length(find(S==T));
Enum=length(find(S==0));%返回true的索引,一列一列来
figure(1);
set(imh, 'cdata', cat(3, isF, isT, zeros(m,n)) )%红绿蓝
str1=['树:',num2str(Tnum),' 火:',num2str(Fnum),' 空:',num2str(Enum)];
str2=['时间:',num2str(t),'次',note];
title({str1,str2});%title{[],[]}换行显示
drawnow%Use this command if you modify graphics objects and want to see the updates on the screen immediately.
pause(0.1);
end
传染病元胞自动机原理
传染病元胞自动机中包含四种元胞状态,分别为易感、潜伏、感染和移除(死亡或治愈)。元胞自动机的规则如下:
如果易感元胞的东、南、西、北四个邻居中有潜伏元胞,则易感元胞下一时刻以概率P0转变为潜伏元胞 。
潜伏元胞经过T0个时间步后转变为感染元胞 。
感染元胞经过C个时间步后转变为移除元胞 。
估计T0=7;C=5;
gamma=1/C=1/5,gamma是治愈率,C是治愈周期。
beta=kb,k是每个病人每天接触的人,b是传染概率,beta是每个病人每天有效接触的人,P0=beta/(gamma4*T0),参数估计内容更详细请见武汉加油——传染病模型拟合以及大揭秘!武汉肺炎传播四大数学模型(附代码)
仿真结果如下:
代码如下:
%%%%%%%%%%%%%%%%参数估计###########33
d=0.03;%死亡率
global C;%治愈时间
C=5;
global N gamma
gamma=1/C;
N=13000000;
t0=0;
t1=35;
t2=41;
I0=0;
I1=log(1182);
I2=log(2758);
global k;
k=1;
global k1;
k1=1.5;
t=[t0,t1,t2];
I=[I0,I1,I2];
p=polyfit(t,I,1);
b=(p(1)+gamma)/k;
global beta beta2;
beta=b*k;
beta2=b*k1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m = 300;
n = 300; % 表示森林的矩阵行列 m x n
%Plight = 5e-6; % 闪电概率
%Pgrowth = 1e-2; % 生长概率
% 邻居方位 d 和点燃概率 p0
global T0;
T0=7;%潜伏天数,R0=beta2/gamma,周围四个有患病危险
p0=(beta2/gamma)/(4*T0);
d1 = {[1,0], [0,1], [-1,0], [0,-1]};%最下一行往最上,最右一列往最左,最前一列往最下,最左一列往最右,circshift(S,d{j})
p1 = [p0,p0,p0,p0];
% % 改进元胞自动机
d2 = {[1,0], [0,1], [-1,0], [0,-1], [1,1], [-1,1], [-1,-1], [1,-1]};%考虑对角
p2 = [ones(1,4), ones(1,4)*(sqrt(1/2)-1/2)];
% 易感=0, 潜伏=1, 患病=2,移除=3
S = 0;
E = 1;
I = 2;
R = 3;
people = zeros(m,n);
people(ceil(m/2),ceil(n/2))=1;
imh = image(cat(3, people==I, people==R, people==E));
axis image;
%plot box fits tightly around the data
E_time=(people==E);
I_time=(people==I);
for t = 1:3000
% 计算邻居中能传播着火的个数
sum = zeros(size(people));
%note=' 传染';
p=p1;
d=d1;
for j = 1:length(d)
sum = sum + p(j) * (circshift(people,d{j})==E);
end
%s上一步
isS = (people==S);
isE = (people==E);
isI = (people==I); % 找出三种不同的状态
isR = (people==R);
if(E*isE+S*isS+I*isI+R*isR~=people)
disp('error');
end
%本步
incubate = rand(m,n)<sum;
E_time=E_time+isE+incubate;
I_time=I_time+isI+(E_time>=T0);
E_time=E_time-T0*(E_time>=T0);
remove=(I_time>=C);%移除条件;
I_time=I_time-C*(I_time>=C) ;
% 变潜伏
Rule1 = E*(E_time~=0);
% 规则 2:变感染
Rule2 =I*(I_time~=0);
% 规则 3: 变移除
Rule3 =R*(remove)+R*(isR); %isR是上一步的
people = Rule1 + Rule2 + Rule3;
Snum=length(find(people==S));
Enum=length(find(people==E));
Inum=length(find(people==I));%返回true的索引,一列一列来
Rnum=length(find(people==R));
figure(1);
set(imh, 'cdata', cat(3, isI, isR, isE) )%红绿蓝
str1=['健康者(黑):',num2str(Snum),' 潜伏者(蓝):',num2str(Enum),' 感染者(红):',num2str(Inum),' 移除者(绿):',num2str(Rnum)];
str2=['时间:',num2str(t),'次'];
title({str1,str2});%title{[],[]}换行显示
drawnow%Use this command if you modify graphics objects and want to see the updates on the screen immediately.
pause(0.1);
end