【元胞自动机】基于matlab元胞自动机模拟HIV传染【含Matlab源码 236期】

一、元胞自动机简介

1 元胞自动机发展历程
最初的元胞自动机是由冯 · 诺依曼在 1950 年代为模拟生物 细胞的自我复制而提出的. 但是并未受到学术界重视.
1970 年, 剑桥大学的约翰 · 何顿 · 康威设计了一个电脑游戏 “生命游戏” 后, 元胞自动机才吸引了科学家们的注意.

1983 年 S.Wolfram 发表了一系列论文. 对初等元胞机 256 种 规则所产生的模型进行了深入研究, 并用熵来描述其演化行 为, 将细胞自动机分为平稳型, 周期型, 混沌型和复杂型.

2 对元胞自动机的初步认识
元胞自动机（CA）是一种用来仿真局部规则和局部联系的方法。典型的元胞自动机是定义在网格上的，每一个点上的网格代表一个元胞与一种有限的状态。变化规则适用于每一个元胞并且同时进行。典型的变化规则，决定于元胞的状态，以及其（ 4 或 8 ）邻居的状态。

3 元胞的变化规则&元胞状态
典型的变化规则，决定于元胞的状态，以及其（ 4 或 8 ）邻居的状态。

4 元胞自动机的应用
元胞自动机已被应用于物理模拟，生物模拟等领域。

5 元胞自动机的matlab编程
结合以上，我们可以理解元胞自动机仿真需要理解三点。一是元胞，在matlab中可以理解为矩阵中的一点或多点组成的方形块，一般我们用矩阵中的一点代表一个元胞。二是变化规则，元胞的变化规则决定元胞下一刻的状态。三是元胞的状态，元胞的状态是自定义的，通常是对立的状态，比如生物的存活状态或死亡状态，红灯或绿灯，该点有障碍物或者没有障碍物等等。

6 一维元胞自动机——交通规则
定义：
6.1 元胞分布于一维线性网格上.
6.2 元胞仅具有车和空两种状态.

7 二维元胞自动机——生命游戏
定义：
7.1 元胞分布于二维方型网格上.
7.2 元胞仅具有生和死两种状态.

元胞状态由周围八邻居决定.
规则：

骷髅：死亡；笑脸：生存
周围有三个笑脸，则中间变为笑脸
少于两个笑脸或者多于三个，中间则变死亡。

8 什么是元胞自动机
离散的系统: 元胞是定义在有限的时间和空间上的, 并且元 胞的状态是有限.
动力学系统: 元胞自动机的举止行为具有动力学特征.
简单与复杂: 元胞自动机用简单规则控制相互作用的元胞 模拟复杂世界.

9 构成要素

（1）元胞 (Cell)

元胞是元胞自动机基本单元：
状态: 每一个元胞都有记忆贮存状态的功能.
离散: 简单情况下, 元胞只有两种可能状态; 较复杂情况下, 元胞具有多种状态.
更新: 元胞的状态都安照动力规则不断更新.
（2）网格 (Lattice)
不同维网格

常用二维网格

（3）邻居 (Neighborhood)

（4）边界 (Boundary)

反射型：以自己作为边界的状态
吸收型：不管边界（车开到边界就消失）

（5）规则（状态转移函数）
定义：根据元胞当前状态及其邻居状况确定下一时刻该元胞状态的动力学函数, 简单讲, 就是一个状态转移函数.
分类 ：
总和型: 某元胞下时刻的状态取决于且仅取决于它所有邻居 的当前状态以及自身的当前状态.
合法型: 总和型规则属于合法型规则. 但如果把元胞自动机 的规则限制为总和型, 会使元胞自动机具有局限性.
（6）森林火灾

绿色：树木；红色：火；黑色：空地。
三种状态循环转化：
树：周围有火或者被闪电击中就变成火。
空地：以概率p变为树木
理性分析：红为火；灰为空地；绿是树

元胞三种状态的密度和为1

火转化为空地的密度等于空地转换为树的密度（新长出来的树等于烧没的树）

f是闪电的概率:远远小于树生成的概率；T s m a x T_{smax}T smax
​是一大群树被火烧的时间尺度
程序实现
周期性边界条件
购进啊

其中的数字为编号
构建邻居矩阵

上面矩阵中的数字编号，对应原矩阵相同位置编号的上邻居编号，一 一对应
同样道理：

（7）交通概念
车距和密度

流量方程

守恒方程

时空轨迹（横轴是空间纵轴为时间）

红线横线与蓝色交点表示每个时间车的位置。
如果是竖线则表示车子在该位置对应的时间

宏观连续模型：

最常用的规则：

红色条表示速度是满的。

1 加速规则：不能超过v m a x （ 2 格 / s ） v_{max}（2格/s）v
max（2格/s）
2 防止碰撞：不能超过车距

理论分析：

结果分析: 密度与流量

第一个图：横坐标是归一化后的密度，纵坐标是车流量。第二个图：理论值与CA的结果

结果分析: 时空轨迹

中间的深色区域是交通堵塞的区域。

二、部分源代码

%% NOTE
 
% Total runtime for this code was about 7.5 hours on the computer used by
% the original researcher (for 1000 run average of 3 different adherence 
% values).
% This runtime is expected to vary from device to device, and is being used
% here to obtain a general comparison for runtime for the different
% models. It is not for a universal runtime value
 
% When run, this code does not display anything at the beginning. However,
% it will produce 2 graphs for each P_adh value every 2.5 hours approx.
 
%% Documentation
 
% n             size of each dimension of our square cellular automata grid
% P_HIV         fraction (probability) of cells initially infected by virus
% P_i           probability of a healthy cell becoming infected if its
%               neighborhood contains 1 I1 cell or X I2 cells
% P_v           probability of a healthy cell becoming infected by coming
%               in contact with a virus randomly (not from its
%               neighborhood)
% P_RH          Probability of a dead cell becoming replaced by a healthy
%               cell
% P_RI          Probability of a dead cell becoming replaced by an infected
%               cell
% P_T1          Probability of a healthy cell receiving therapy 1
% P_iT1         Probability of a healthy cell receiving therapy 1 of
%               becoming infected
% P_T2          Probability of a healthy cell receiving therapy 2
% P_iT2         Probability of a healthy cell receiving therapy 2 of
%               becoming infected
% X             Number of I2 cells in the neighborhood of an H cell that
%               can cause it to become infected
% tau1          tau1 is the number of timesteps it takes for an acute
%               infected cell to become latent.
% tau2          tau2 is the number of timesteps it takes for a latent
%               infected cell to become dead.
% tau3          tau3 is the number of timesteps after which a healthy
%               cell receiving dual therapy becomes a healthy cell again
% totalsteps    totalsteps is the total number of steps of the CA (the
%               total number of weeks of simulations)
% grid          our cellular automata (CA) grid
% tempgrid      tempgrid is a temporary grid full of random numbers that is
%               used to randomly add different states to our CA grid.
% taugrid       taugrid is a grid the same size as our CA grid that stores
%               the number of timesteps that a cell has been in state I_1.
%               If the number reaches tau1, then the state changes to I_2.
% state         state is a [9 x totalsteps] size matrix that stores
%               the total number of cells in each state at each timestep
%               and the last 2 rows store total healthy and total infected
%               cells
% timestep      each simulation step of the cellular automata
%               1 timestep = 1 week of time in the real world
% nextgrid      nextgrid is a temporary grid. It is a copy of the CA grid
%               from the previous simulation. It stores all the CA rule
%               updates of the current timestep and stores it all back to
%               the grid to display.
 
%% Clean-up
 
clc;            % clears command window
clear all;      % clears workspace and deletes all variables
% close all;      % closes all open figures
 
%% Parameters
 
n = 100;            % meaning that our grid will have the dimensions n x n
P_HIV = 0.05;       % initial grid will have P_hiv acute infected cells
P_i = 0.997;        % probability of infection by neighbors
P_v = 0.00001;      % probability of infection by random viral contact
P_RH = 0.99;        % probability of dead cell being replaced by healthy
P_RI = 0.00001;     % probability of dead cell being replaced by infected
P_T1 = 0.70;        % probability of cell receiving therapy 1
P_iT1 = 0.07;       % probability of infection of healthy with therapy 1
P_T2 = 0.50;        % probability of cell receiving therapy 2
P_iT2 = 0.05;       % probability of infection of healthy with therapy 2
X = 4;              % there must be at least X I_2 neighbors to infect cell
tau1 = 4;           % time delay for I_1 cell to become I_2 cell
tau2 = 1;           % time delay for I_2 cell to become D cell
tau3 = 1;           % time delay for H_Tb cell to become H cell
totalsteps = 600;   % total number of weeks of simulation to be performed
T_start = 20;       % The medication therapy will start on week T_start
totalruns = 1000;   % total number of times to run the simulation to get an
% average
 
%% States
 
% State 1: H:       Healthy                     (Color- Green)
% State 2: H_T1:    Healthy with therapy 1      (Color- Red)
% State 3: H_T2:    Healthy with therapy 2      (Color- Red)
% State 4: H_Tb:    Healthy with dual therapy   (Color- Red)
% State 5: I_1:     Active Infected             (Color- Cyan)
% State 6: I_2:     Latent Infected             (Color- Blue)
% State 7: D:       Dead                        (Color- Black)
 
%% Simulation
 
for P_adh = 0.5:0.2:0.9
    
    state1 = zeros(totalruns,totalsteps);
    state2 = zeros(totalruns,totalsteps);
    state3 = zeros(totalruns,totalsteps);
    state4 = zeros(totalruns,totalsteps);
    state5 = zeros(totalruns,totalsteps);
    state6 = zeros(totalruns,totalsteps);
    state7 = zeros(totalruns,totalsteps);
    
    for run = 1:totalruns
    
    
    state1dev = std(state1);
    state2dev = std(state2);
    state3dev = std(state3);
    state4dev = std(state4);
    state5dev = std(state5);
    state6dev = std(state6);
    state7dev = std(state7);
    state8dev = state1dev + state2dev + state3dev + state4dev;
    state9dev = state5dev + state6dev;
    
    state1mean = mean(state1);
    state2mean = mean(state2);
    state3mean = mean(state3);
    state4mean = mean(state4);
    state5mean = mean(state5);
    state6mean = mean(state6);
    state7mean = mean(state7);
    state8mean = state1mean + state2mean + state3mean + state4mean;
    state9mean = state5mean + state6mean;
    
    % The following lines of code are to display a graph of each state of
    % cells during simulation
    set(figure, 'OuterPosition', [200 100 700 500]) % sets figure window size
    plot( 1:totalsteps , state1mean, 'g', ...
        1:totalsteps , state2mean, 'r', ...
        1:totalsteps , state3mean, 'm', ...
        1:totalsteps , state4mean, 'y', ...
        1:totalsteps , state5mean, 'c', ...
        1:totalsteps , state6mean, 'b', ...
        1:totalsteps , state7mean, 'k' , 'linewidth', 2 );
    xlim([0 600]);
    ylim([0 10000]);
    hold on;
    gridxy(20,'Color',[0.8 0.5 0.0],'linewidth',5) ;
    errorbar( 1:15:totalsteps , state1mean(1:15:totalsteps), state1dev(1:15:totalsteps), 'g');
    errorbar( 1:15:totalsteps , state2mean(1:15:totalsteps), state2dev(1:15:totalsteps), 'r');
    errorbar( 1:15:totalsteps , state3mean(1:15:totalsteps), state3dev(1:15:totalsteps), 'm');
    errorbar( 1:15:totalsteps , state4mean(1:15:totalsteps), state4dev(1:15:totalsteps), 'y');
    errorbar( 1:15:totalsteps , state5mean(1:15:totalsteps), state5dev(1:15:totalsteps), 'c');
    errorbar( 1:15:totalsteps , state6mean(1:15:totalsteps), state6dev(1:15:totalsteps), 'b');
    errorbar( 1:15:totalsteps , state7mean(1:15:totalsteps), state7dev(1:15:totalsteps), 'k');
    legend( 'Therapy start week', 'Healthy', 'Healthy with Therapy1', 'Healthy with Therapy2', ...
        'Healthy with dual Therapy', 'Acute Infected', 'Latent Infected', ...
        'Dead', 'Location' ,'NorthEast' );
    saveas(gcf,strcat('Model3withAdherencePadh',num2str(P_adh*100),'Graph1.pdf'));
    
    
    % The following lines of code are to display a graph of each state of
    % cells during simulation
    set(figure, 'OuterPosition', [200 100 700 500]) % sets figure window size
    plot( 1:totalsteps , state8mean, 'g', ...
        1:totalsteps , state9mean, 'b', ...
        1:totalsteps , state7mean, 'k' , 'linewidth', 2 );
    xlim([0 600]);
    ylim([0 10000]);
    hold on;
    gridxy(20,'Color',[0.8 0.5 0.0],'linewidth',5) ;
    errorbar( 1:15:totalsteps , state8mean(1:15:totalsteps), state8dev(1:15:totalsteps), 'g');
    errorbar( 1:15:totalsteps , state9mean(1:15:totalsteps), state9dev(1:15:totalsteps), 'b');
    errorbar( 1:15:totalsteps , state7mean(1:15:totalsteps), state7dev(1:15:totalsteps), 'k');
    legend( 'Therapy start week', 'Healthy', 'Infected', 'Dead', ...
        'Location' ,'NorthEast' );
    saveas(gcf,strcat('Model3withAdherencePadh',num2str(P_adh*100),'Graph2.pdf'));
    
end

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
  
8
  
9
  
10
  
11
  
12
  
13
  
14
  
15
  
16
  
17
  
18
  
19
  
20
  
21
  
22
  
23
  
24
  
25
  
26
  
27
  
28
  
29
  
30
  
31
  
32
  
33
  
34
  
35
  
36
  
37
  
38
  
39
  
40
  
41
  
42
  
43
  
44
  
45
  
46
  
47
  
48
  
49
  
50
  
51
  
52
  
53
  
54
  
55
  
56
  
57
  
58
  
59
  
60
  
61
  
62
  
63
  
64
  
65
  
66
  
67
  
68
  
69
  
70
  
71
  
72
  
73
  
74
  
75
  
76
  
77
  
78
  
79
  
80
  
81
  
82
  
83
  
84
  
85
  
86
  
87
  
88
  
89
  
90
  
91
  
92
  
93
  
94
  
95
  
96
  
97
  
98
  
99
  
100
  
101
  
102
  
103
  
104
  
105
  
106
  
107
  
108
  
109
  
110
  
111
  
112
  
113
  
114
  
115
  
116
  
117
  
118
  
119
  
120
  
121
  
122
  
123
  
124
  
125
  
126
  
127
  
128
  
129
  
130
  
131
  
132
  
133
  
134
  
135
  
136
  
137
  
138
  
139
  
140
  
141
  
142
  
143
  
144
  
145
  
146
  
147
  
148
  
149
  
150
  
151
  
152
  
153
  
154
  
155
  
156
  
157
  
158
  
159
  
160
  
161
  
162
  
163
  
164
  
165
  
166
  
167
  
168
  
169
  
170
  
171
  
172
  
173
  
174
  
175
 

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例（第2版）[M].电子工业出版社，2016.
[2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社，2017.
[3]【数学建模】元胞自动机.博主：二进制 人工智能

文章来源: qq912100926.blog.csdn.net，作者：海神之光，版权归原作者所有，如需转载，请联系作者。

原文链接：qq912100926.blog.csdn.net/article/details/113724445