1 2 3

生态系统建模

本文介绍如何用 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,绘制了各生物数量的变化状态以及位置变化情况。

nature system

只要保证几个核心的变量和全局变量(比如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

这个语句判断了一只狼是否是索引为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),可得到两幅图。

一幅表明生物数量的变化情况和狼在每次迭代中吃掉麋鹿的数量。

state

另一幅表明麋鹿、头狼和狼的空间位置。在最后一个时刻,位置图像如下。

state

注:在函数ecolab(30 ,100, 5, 4, 60, false, false)中,各变量分别表示:地理环境的大小,麋鹿初始数量,每个狼群中普通狼的数量,狼群数量,迭代次数,是否采用快速模式,是否保存每次迭代产生的图像为文件。