本文介绍如何用 Agent-based model 的方法，对一个包含狼、麋鹿、草三种生物的生态系统建模。
GitHub 项目地址：nature-system
Agent指个体，Agent-base model即基于个体的模型。在Agent-based model中，虚拟实体(virtual entity)和现实个体一一对应。一般来说，虚拟实体的属性是对现实个体的抽象，虚拟实体的行为则是依照已经写好的规则迭代。这种建模方式往往涉及个体详细的行为模式，因此适合以研究个体行为规则为目的的建模。
需要注意的是，利用Agent-based model方法建模，并不是变量越多越好，过多的变量会造成调试上的困难。选取的变量足够支撑起我们研究的问题就可以了。
Note: 除了Agent-base model之外，Equation-based model也是一种比较常见的建模方式。Agent-based model适合表达每个个体身上发生了什么，比如个体的速度、位置的变化等。Equation-based model则而更适合表达系统平均发生了什么，比如人口模型中人口数量变化就更适合用Equation-based model建模。
模型假设
（一）生物的属性
模型内有三种生物：狼(wolf)、麋鹿(moose)、草(grass)。此外，狼还分为头狼和普通狼。下面这张表展示了头狼、普通狼、麋鹿的属性：
| # | age | food | pos | speed | last_breed | population | packNo | range | eaten |
|---|---|---|---|---|---|---|---|---|---|
| 头狼 | √ | √ | √ | √ | √ | √ | √ | √ | |
| 普通狼 | √ | √ | √ | √ | √ | ||||
| 麋鹿 | √ | √ | √ | √ | √ |
从上表可以看出，年龄、食物、位置、速度是头狼、普通狼和麋鹿的共有属性。last_breed指和上一次生育的时间间隔，在模型内用一次迭代表示一次时间间隔，因为迭代必须是整数，因此last_breed也是整型变量。头狼特有population属性，用来表示该头狼所在狼群中狼的数量；头狼特有range属性，表示该狼群活动的范围；头狼特有eaten属性，表示狼群在一次迭代中，全体成员吃到的麋鹿的总数。普通狼特有packNo属性，用来标记其所属狼群的头狼在列表中的索引号。普通狼没有last_breed属性，是因为其生育是以整个狼群为单位的。由于狼群是由头狼定义的，狼群相关的属性和行为也绑定在头狼身上，因此在代码实现上，狼群的生育是通过头狼实现的。
草没有上述生物那么复杂的模型，本模型中的草均匀平铺在环境中，且不会生长，吃完就没了。因此麋鹿为了获取新的食物来源，需要经常迁徙。
（二）生物的行为
头狼、普通狼和麋鹿都具有四种行为：死亡、迁徙、生育、捕猎。这四种行为在不同的生物种类上有不同的定义。比如麋鹿的迁徙规则是等到自己的位置没草吃了，才开始迁徙。而头狼迁徙规则则比较复杂，头狼会根据狼群内的种群数量，判断狼群的散布范围，以及调整迁徙的积极性。
1.死亡(die)
狼有两种死法：饿死或老死。麋鹿则有三种死法：饿死、老死或者被吃。
每一种生物都有固定的年龄上限值，超过该值将在本次迭代中死亡。本模型中，狼的最大年龄为25。如果一只狼如果没有饿死，它将在第25次迭代时老死。而饿死在生物food属性为0时发生。
对于麋鹿，每只麋鹿都有一个food属性。对于狼，整个狼群共享一个food属性，当狼群的food属性值为0时，会有一只普通狼死亡。
如果头狼死去时，狼群内普通狼的数量大于0，则会有一只普通狼晋升为头狼。
2.生育(breed)
不同的物种有不同的生育规则。
每只麋鹿一次只生一个。生育时母亲将把自己food属性值的一半分给孩子。生育有时间间隔，不能持续生育。在本模型中，麋鹿经历一次生育后，起码要经历10次迭代，才能进行第二次生育。此外，生育还有食物限制，当food属性高于10，麋鹿才会生育。如果同时满足时间间隔和食物条件，那么麋鹿将在本轮迭代进行生育。
本模型中，一个狼群为一个生育单位。和麋鹿的生育相似，狼群也仅在食物有富余时进行生育。生育的数量与狼群群内个体数量成正比。在不同群内生物数量范围内，比例系数略有不同。生育后，母体的食物也会按比例减少，并将此部分食物赋予孩子。同样，狼群的生育也有时间间隔限制。一次生育完成后，起码要经历5次迭代的时间，才能进行下一次生育。
3.迁徙(migrate)
狼群的迁徙取决于头狼，头狼相当于狼群的神经中枢。它接收上一轮狼群的捕猎信息，来判断下一步狼群的迁徙方向。由于头狼无法获知整张地图的麋鹿分布信息，它仅仅只能从已知信息推断，所以这可以视为局部优化问题。
具体来说，狼群的迁徙规则分两类：如果在上一次迭代中，狼群没有捕到一只麋鹿，那么狼群将随机移动以寻找麋鹿。如果起码吃到了一只，那么就用一个for循环，找出是狼群中的哪只狼吃到了麋鹿。在下一步迭代中，头狼将往这只狼的方向移动。
头狼移动之后，狼群内的普通狼将在头狼的range属性值代表的活动范围内随机分布。range的数值取决于狼群群内个体数量，群内个体数量大的狼群，其range值也较大。
麋鹿当且仅当在所在位置的草吃完以后菜进行移动。麋鹿吃完草以后，如果在搜索范围内有草，则向该方向移动，否则随机移动以进行搜索。
PS: 本模型有地理范围限制，所有生物均不可以移动到地图之外的地方。因此在代码实现时，所有生物在迁徙之前，都必须验证移动位置的有效性。
4.捕猎(eat)
无论是狼还是麋鹿，食物水平都会随着时间流逝而下降。
对于麋鹿，每次迭代中如果没有吃到食物，它的food属性值将会减1。吃到食物时，food属性值则会加2，但减去随时间流逝下降的1点食物值，在数值变动上实际只增加了1。
对于狼，它捕食麋鹿为生。在本模型中，狼的速度值与搜索范围值相等。如果一只麋鹿出现在狼的搜索范围内，则麋鹿和狼的距离越近，越容易被判定为被狼捕获。代码实现如下：
pk = 1 - (d/spd)
if pk > rand
{麋鹿死，狼食物加1，狼移动到麋鹿原先的位置}
end
其中，pk是麋鹿被捕杀的概率。d是麋鹿和狼之间的距离，spd是狼的速度，rand是一个介于0到1之间的随机数。距离(d)越短，狼的速度(spd)越大，pk值就越大，pk > rand的可能性就越大，因此麋鹿被捕获的概率就越高。
吃到麋鹿时，狼的food属性值加2。同时因为本轮迭代下降的1点food值，实际只增长1点food值。
整体架构
本模型包含有34个m文件，总计1846行代码，难以详述，具体代码实现请看GitHub仓库，本文只能说明大致框架。
下面用流程图表示了各函数之间的调用关系。其中ecolab.m是主函数。左边一系列以create打头的函数搭建了模型的基本参数和数据结构。agent_solve是一个核心函数，它通过调用各生物的eat, breed, migrate, die函数，实现了对各生物行为和状态的建模。此外，它还通过调用update_messages，在每轮迭代后更新各生物在全局变量MESSAGES中的对应值。最后，它通过调用plot_reasult，绘制了各生物数量的变化状态以及位置变化情况。
只要保证几个核心的变量和全局变量（比如agent和MESSAGES）在每轮迭代中正常更新，那么即使其他部分有错，依旧能保证模型的正常运行。
代码实现
1.下面这段代码用于维护各生物的属性数据，每轮迭代前都要执行。
function [agent] = initial_iteration(agent,nf)
% This a new function to initial the property at the beginning of each iteration
global MESSAGES PARAM
%% initial alphaWolf.eaten at the beginning of each iteration
%eaten is the number of moose be eaten at each iteration, this
%number will be useful to decision making process of migrate.
for p = 1:length(agent)
if isa(agent{p},'alphaWolf')
agent{p}.eaten = 0;
end
end
%% initial wolf.migration at the beginning of each iteration
% wolf.migration
for q = 1:length(agent)
if isa(agent{q},'wolf')
agent{q}.migration = 0;
end
end
%% update the alphaWolf.population of wolf pack for each iteration
for m = 1:length(agent)
currAgent = agent{m};
if isa(currAgent,'alphaWolf')
count = 0;
for n = 1:length(agent)
if isa(agent{n},'wolf') & agent{n}.packNo == m & MESSAGES.atype(n) ~= 0
count = count + 1;
end
end
currAgent.population = count;
agent{m} = currAgent;
end
end
%% updatet alphaWolf.range at the beginning of each iteration
for m = 1:length(agent)
currAgent = agent{m};
if isa(currAgent,'alphaWolf')
if currAgent.population <= PARAM.F_SIZE
nrange = ceil(PARAM.F_SPD/2 * (1 + nf/5));
else
nrange = PARAM.F_SPD;
end
currAgent.range = nrange;
agent{m} = currAgent;
end
end
从这段代码中可以看到一些整个模型都在反复使用的语句。比如：
if isa(agent{p},'alphaWolf')
对于$isa()$，MATLAB帮助文档是这么解释的：tf=isa(obj, ClassName) returns true if obj is an instance of the class specified by ClassName, and false otherwise. isa also returns true if obj is an instance of a class that is derived from ClassName.
也就是说，$isa()$用于判断一个变量是否是一个类的实例，或是否是一个类的子类的实例。
还有一个常用的判断是：
if agent{n}.packNo == m
2.下面这段代码出自从agent_solve中。
for cn=n:-1:1
curr=agent{cn};
if isa(curr,'moose')|isa(curr,'wolf')|isa(curr,'alphaWolf')
%% eat rules
% moose
if isa(curr,'moose')
[curr,moose_eaten]=eat(curr,cn);
end
% wolf
if isa(curr,'wolf') %alphaWolf can not chase prey
[curr,agent]=eat(curr,cn,agent); %eating rules (mooses eat food, wolves eat mooses)
end
if isa(curr,'alphaWolf')
eaten=curr.eaten;
end
%% migrate rule
% if population is large, wolf pack will be active
flag = 0;
if isa(curr,'alphaWolf') & curr.population > PARAM.F_SIZE & eaten < ceil(curr.population * 1/2)
[curr,agent_migrate]=migrate(curr,cn,agent,eaten); % if no food was eaten, then migrate in search of some
agent = agent_migrate; %up date cell array with modified agent data structure
flag =1;
% if population is small, wolf pack will not move so frequently
elseif isa(curr,'alphaWolf') & curr.population <= PARAM.F_SIZE & eaten < ceil(curr.population * 1/4)
[curr,agent_migrate]=migrate(curr,cn,agent,eaten); % if no food was eaten, then migrate in search of some
agent = agent_migrate; %up date cell array with modified agent data structure
flag =1;
elseif isa(curr,'moose') & moose_eaten==0
curr=migrate(curr,cn);
end
% if the alphawolf did not migrate,flag is 0
% in this case,wolf should move randomly to find the food
if isa(curr,'alphaWolf') & flag == 0
for p = 1:length(agent)
if isa(agent{p},'wolf') & agent{p}.packNo == cn & MESSAGES.atype(p) ~= 0 & agent{p}.migration == 0
[agent{p}] = migrate(agent{p},p,agent);
end
end
end
%% death rule
% wolf (from old age)
% moose (from starvation or old age)
if isa(curr,'moose')|isa(curr,'wolf')
[curr,klld]=die(curr,cn);
end
% alphaWolf (from old age) and wolf pack (from starvation)
% PS: alphaWolf represent for the whole wolf pack here
if isa(curr,'alphaWolf')
[curr,klld,agent]=die(curr,cn,agent);
end
%% breed rule
% moose
if isa(curr,'moose')&klld==0
new=[];
[curr,new]=breed(curr,cn); %breeding rule
if ~isempty(new) %if current agent has bred during this iteration
n_new=n_new+1; %increase new agent number
agent{n+n_new}=new; %add new to end of agent list
end
end
% alphaWolf
% only alphaWolf can breed, wolf can not
if isa(curr,'alphaWolf')&klld==0
new=[];
[curr,new]=breed(curr,cn); %breeding rule
if ~isempty(new) %if current agent has bred during this iteration
for ite = 1:length(new)
n_new=n_new+1; %increase new agent number
agent{n+n_new}=new{ite}; %add new to end of agent list
end
end
end
agent{cn}=curr; %up date cell array with modified agent data structure
end
end
由于一些函数的输出是另一些函数的输入，因此各个函数的调用顺序需要谨慎调整。本模型按照eat, migrate, die, breed的顺序依次调用上述函数。
基于同样的原因，上面这段代码中第一句for循环，写成了倒序形式：
for cn=n:-1:1
运行结果
在命令窗口中输入ecolab(30 ,100, 5, 4, 60, false, false)，可得到两幅图。
一幅表明生物数量的变化情况和狼在每次迭代中吃掉麋鹿的数量。
另一幅表明麋鹿、头狼和狼的空间位置。在最后一个时刻，位置图像如下。
注：在函数ecolab(30 ,100, 5, 4, 60, false, false)中，各变量分别表示：地理环境的大小，麋鹿初始数量，每个狼群中普通狼的数量，狼群数量，迭代次数，是否采用快速模式，是否保存每次迭代产生的图像为文件。