再探猎食者猎物模型

在本节中,我们将重温第一章的猎食者猎物系统。不过,这一次我们将利用单独的组件创建系统模型。在重复第一章得到的行为后,我们将扩大考虑的范围。我们会重新配置这些组件模型,以创建系统模型去演示不同的动态等。

经典猎食者猎物

我们首先观察经典的猎食者猎物系统。为了使用组件模型来建立这样一个系统,我们将需要不同模型。这包括兔子和狐狸的人口模型,以及代表繁殖、饥饿和捕食的模型。

连接器

但是,我们在对连接器的讨论中了解到了一点。在开始建立组件模型前,我们必须先定义相互作用的部件间所交换的信息。要做到这点,我们要定义连接器。在本节我们会使用的connectorSpecies连接器。它的定义如下:

within ModelicaByExample.Components.LotkaVolterra.Interfaces;
connector Species "Used to represent the population of a specific species"
  Real population "Animal population";
  flow Real rate "Flows that affect animal population";
end Species;

此连接器的定义很有趣。因为这些定义并非来自工程领域。相反,这些定义其实来源于生态学。在这种情况下,我们的横跨变量是population。这代表特定物种动物的实际数量。flow限定词所标示的通过变量是rate。此变量表示,连接器所在部件的新动物“输入”率。

区域种群大小

要追踪一个地区内一个特定物种的数目,我们将使用RegionalPopulation模型。这个模型有几个值得注意的方面。所以我们会分步介绍模型。首先:

within ModelicaByExample.Components.LotkaVolterra.Components;
model RegionalPopulation "Population of animals in a specific region"
  encapsulated type InitializationOptions = enumeration(
      Free "No initial conditions",
      FixedPopulation "Specify initial population",
      SteadyState "Population initially in steady state");

The first two lines are as expected.前两行没有意外。但此后,我们看到模型定义了一个名为InitializationOptions的类型。类型定义前添加了encapsulated关键字。这是必要的。因为该类型是在模型而不是包内定义的。Modelica有这样的一条规则。如果我们想从model定义外引用此类型,类型定义前必须添加encapsulated。从定义可以看到,enumeration定义了三个不同的值:FreeFixedPopulation以及SteadyState。我们将很快看到如何使用这些值。

RegionalPopulation模型的第一个声明为:

  parameter InitializationOptions init=InitializationOptions.Free
    annotation(choicesAllMatching=true);

需要注意的是init,也就是第一个parameter,利用了InitializationOptions枚举值。不但以此指定了其类型(为枚举值),而且也以此指定了起初始值Free。还要注意这里有choicesAllMatching标注。我们会在本章的后面回顾概念时以及随后的章节里更多地讨论annotation

紧接的声明是:

  parameter Real initial_population
    annotation(Dialog(group="Initialization",
    enable=init==InitializationOptions.FixedPopulation));

initial_population参数是用来表示在该区域在仿真开始时种群大小的初始值。然而,我们在不久后就会在等式区域里看到,只有在init值设为FixedPopulation时,模型才会使用上述值。出于这个原因,parameterenable标注设定为init==InitializationOptions.FixedPopulation。这个标注把上述关系告知了Modelica工具。在建立与该模型相关联的图形用户界面(例如参数的对话框)时,工具就可以使用该信息了。

另外值得注意的一点是initial_population定义内的Dialog标注。该标注允许模型开发者把参数组织成不同的类别。在这种情况下参数在“初始化”类别里。Modelica工具通常会使用这些信息去帮助组织参数对话框。

模型的最后一次公有声明为connector实例。该连接器可以让组件与外界其他组件进行相互作用:

  Interfaces.Species species
    annotation (Placement(transformation(extent={{-10,90},{10,110}}),
        iconTransformation(extent={{-10,90},{10,110}})));

这里我们再次看到了组件位置标注。不过,我们还是暂时不对其进行讨论。这样就剩下最后一条声明。这条声明恰好为protected

protected
  Real population(start=10) = species.population "Population in this region";

这个变量表示在这个地区的动物总数。变量的start值非零。这样做是为了避免前面对猎食者猎物系统稳定状态初始化讨论里所看到的平凡解。从这个声明我们也可以看到,该声明令局部变量populationspecies连接器内横跨变量species.population的值相等。结果,population变量其实是species.population表达式的别名。

以上就是所有声明。现在,让我们来看看RegionalPopulation模型的相关公式:

initial equation
  if init==InitializationOptions.FixedPopulation then
    population = initial_population;
  elseif init==InitializationOptions.SteadyState then
    der(population) = 0;
  else
  end if;
equation
  der(population) = species.rate;
  assert(population>=0, "Population must be greater than zero");
end RegionalPopulation;

initial equation展示了init参数值对组件行为的作用。在一些情况下,init值等于在InitializationOptions枚举内的FixedPopulation值。此时模型会引入一个方程,以指定在仿真开始时的population变量的值等于initial_population参数。在另一些情况下,init值等于SteadyState枚举值。这时模型引入等式,以指定仿真开始时的种群大小变化速率必须为零。很显然,如果init等于Free(最后剩下的可能选择),那模型就没有初始方程。

equation区域,我们看到population变化率等于species连接器内flow变量species.rate的值。再一次,我们要记得变量的正负号约定。flow变量正值意味着流入组件内。上述公式与正负号约定相一致(即若species.rate为正,则会让区域内population增加)。

值得一提的最后一点是模型是在equation区域内的assert调用。这是用来定义模型的边界(即令模型中的方程无效的点)。该语句用于实行区域内人口不能小于零这一约束。若某个解被发现违反约束,模型就会产生“人口必须大于零”的错误信息。

这个模型在定义里也有相关的Icon标注。像往常一样,Icon标注不包含在显示。但是,当这个组件模型是在系统模型内渲染后,其图标看起来就像这样:

繁殖

我们将研究的第一个真正的效应是繁殖。从前面的讨论我们知道,由繁殖带来的种群数量增长正比于该区域中动物的数量。因此,我们可以非常简洁地将繁殖描述为:

within ModelicaByExample.Components.LotkaVolterra.Components;
model Reproduction "Model of reproduction"
  extends Interfaces.SinkOrSource;
  parameter Real alpha "Birth rate proportionality constant";
equation
  growth = alpha*species.population "Growth is proporational to population";
end Reproduction;

其中alpha为比例常数。不过,这个模型之所以简单明了,主要要归功于其对SinkOrSource模型的继承。此前的“DRY”电气部件也得益于对TwoPin模型的继承。因此在这点上,两个继承的作用别无二致。

SinkOrSource模型任何增加或减少动物数的模型的基点。其定义如下:

within ModelicaByExample.Components.LotkaVolterra.Interfaces;
partial model SinkOrSource "Used to describe single species effects"

  Species species
    annotation (Placement(transformation(extent={{-10,90},{10,110}})));
protected
  Real growth "Growth in the population (if positive)";
  Real decline "Decline in the population (if positive)";
equation
  decline = -growth;
  species.rate = decline;
end SinkOrSource;

要理解这些方程,首先要理解一点。任何以SinkOrSource为基础extends的模型,通常都将连接到RegionalPopulation实例(但本身不会一个RegionalPopulation模型)。这意味着,如果实例的flow变量species.rate为正时,模型将动物从RegionalPopulation模型内抽出来。通过观察SinkOrSource模型,我们可以看到变量decline只是species.rate的一个别名。换句话说,若decline为正那么species.rate也会为正。因此,任何与此SinkOrSource实例连接的RegionalPopulation其种群大小都会受到损失。相反,当species.rate为负时,growth变量就会为正。在这种情况下,所连接的RegionalPopulation其种群大小将会增加。

通过定义并继承SinkOrSource模型可以隐藏大部分复杂性。这样一来,像Reproduction这样的模型可以将方程写得更直观,更好地反映模型的行为。例如:growth = alpha*species.population

虽然先前没有显示,IconReproduction模型如下所示:

饥饿

正如先前介绍的Reproduction模型一样,Starvation模型也继承了SinkOrSource模型。不过,模型与decline变量有关的行为可以描述如下:

within ModelicaByExample.Components.LotkaVolterra.Components;
model Starvation "Model of starvation"
  extends Interfaces.SinkOrSource;
  parameter Real gamma "Starvation coefficient";
equation
  decline = gamma*species.population
    "Decline is proporational to population (competition)";
end Starvation;

猎食

建立经典猎食者猎物系统模型前,我们需要考虑最后一个效应是猎食模型。

回想此前的SinkOrSource模型以及与正负号约定相关的潜在问题等等的讨论。SinkOrSource模型设计时仅仅考虑了与一个RegionalPopulation进行交互(因为模型只有一个Species连接器)。为了解决可能的正负号规定混淆问题,对于涉及两个不同地区种群之间相互效应,我们定义了下面的partial模型Interaction

within ModelicaByExample.Components.LotkaVolterra.Interfaces;
partial model Interaction "Used to describe interactions between two species"
  Species a "Species A"
    annotation (Placement(transformation(extent={{-110,-10},{-90,10}})));
  Species b "Species B"
    annotation (Placement(transformation(extent={{90,-10},{110,10}})));
protected
  Real a_growth "Growth in population of species A (if positive)";
  Real a_decline "Decline in population of species A (if positive)";
  Real b_growth "Growth in population of species B (if positive)";
  Real b_decline "Decline in population of species B (if positive)";
equation
  a_decline = -a_growth;
  a.rate = a_decline;
  b_decline = -b_growth;
  b.rate = b_decline;
end Interaction;

再一次,我们定义了有增长变量以及减少变量。但是这一次,有每个变量均有两个版本。其中一个与a连接器相关联,另一个则与b连接器相关联。

使用上述定义,我们很容易可以定义下面的Predation

within ModelicaByExample.Components.LotkaVolterra.Components;
model Predation "Model of predation"
  extends Interfaces.Interaction;
  parameter Real beta "Prey (Species A) consumed";
  parameter Real delta "Predators (Species B) fed";
equation
  b_growth = delta*a.population*b.population;
  a_decline = beta*a.population*b.population;
end Predation;

此模型描述了,“B”(掠食者)种群大小增长与捕食者以及猎物种群大小两者乘积成正比。同样,“A”(猎物)种群大小的减少也与捕食者以及猎物的种群大小这两者乘积成正比(虽然使用的是另一个比例常数)。

虽然先前没有显示,Predation模型的Icon如下所示:

需要注意的是Predation模型并不对称的。b连接器应连接到猎食者种群。而a连接器应连接到猎物种群。模型的图像以及不对称的图标也强调了这一点。

经典系统模型

有了上述组件,我们就可以很容易地构建面向组件版本的经典猎食者猎物系统。通过拖放组件,可以得到以下系统配置:

Component-oriented version of the classic Lotka-Volterra model

这里我们看到Starvation模型连接到foxes种群上。而Reproduction模型则连接到了rabbits种群。Predation模型连接则同时连接到两个种群上。a(猎物)连接器连接到rabbits处。而b(捕食)连接器则连接到foxes处。

我们可以从下面的图中看出,这个系统的行为等同于之前讨论的猎食者猎物系统的结果:

../../../_images/CLV.png

引入第三个物种

我们将一而再地看到,建立组件模型对比直接写出其封装的公式而言,需要一定的初始投入。但这个过程也有显著的“回报”。因为作为结果,我们就可以基于概要地去组合系统。例如对于猎食者猎物系统,我们可以在经典的猎食者猎物系统添加第三个物种狼。

引入狼群

要建立带有第三个物种模型,并不需要定义任何其他组件模型。相反,我们不但可以重用现有的PredationStarvationRegionalPopulation模型,我们还可以重用ClassicLotkaVolterra模型本身:

within ModelicaByExample.Components.LotkaVolterra.Examples;
model ThirdSpecies "Adding a third species to Lotka-Volterra"
  import ModelicaByExample.Components.LotkaVolterra.Components.RegionalPopulation.InitializationOptions.FixedPopulation;
  extends ClassicLotkaVolterra(rabbits(initial_population=25), foxes(initial_population=2));
  Components.RegionalPopulation wolves(init=FixedPopulation, initial_population=4)
    annotation (Placement(transformation(extent={{30,40},{50,60}})));
  Components.Starvation wolf_starvation(gamma=0.4)
    annotation (Placement(transformation(extent={{70,40},{90,60}})));
  Components.Predation wolf_predation(beta=0.04, delta=0.08) "Wolves eating Foxes"
    annotation (Placement(transformation(extent={{-10,-10},{10,10}}, rotation=90, origin={60,20})));
  Components.Predation wolf_rabbit_predation(beta=0.02, delta=0.01) "Wolves eating rabbits"
    annotation (Placement(transformation(extent={{-10,30},{10,50}})));
equation
  connect(wolf_predation.b, wolves.species) annotation (Line(
      points={{60,30},{60,80},{40,80},{40,60},{40,60}},
      color={0,127,0},
      smooth=Smooth.None));
  connect(wolf_rabbit_predation.a, rabbits.species) annotation (Line(
      points={{-10,40},{-40,40},{-40,-20}},
      color={0,127,0},
      smooth=Smooth.None));
  connect(wolf_predation.a, foxes.species) annotation (Line(
      points={{60,10},{60,0},{40,0},{40,-20}},
      color={0,127,0},
      smooth=Smooth.None));
  connect(wolf_starvation.species, wolves.species) annotation (Line(
      points={{80,60},{80,80},{40,80},{40,60}},
      color={0,127,0},
      smooth=Smooth.None));
  connect(wolves.species, wolf_rabbit_predation.b) annotation (Line(
      points={{40,60},{40,80},{20,80},{20,40},{10,40}},
      color={0,127,0},
      smooth=Smooth.None));
  annotation (experiment(StopTime=100, Tolerance=1e-006));
end ThirdSpecies;

你通常不会被通过键入上述源代码去建立这样的模型。相反,在一个图形化开发环境里,通过扩大现有的ClassicLotkaVolterra模型去组装这样的系统,将需要不到一分钟。可视化后,生成系统的示意图为如下所示:

The Classic Lotka-Volterra model augmented with an additional predatory species

变更后的动态

通过建立这样一个模型,我们可以很容易观察经典模型与三种群模型在系统动力学之间的差异。下面的图显示了这三个物种是进行如何相互作用的:

../../../_images/ThirdS.png

通过在不同的RegionalPopulation实例内的init参数,我们还可以快速创建一个模型来求解三个种群的在平衡态的大小。

within ModelicaByExample.Components.LotkaVolterra.Examples;
model ThreeSpecies_Quiescent "Three species in a quiescent state"
  import ModelicaByExample.Components.LotkaVolterra.Components.RegionalPopulation.InitializationOptions.SteadyState;
  extends ThirdSpecies(
    rabbits(init=SteadyState, population(start=30)),
    foxes(init=SteadyState, population(start=2)),
    wolves(init=SteadyState, population(start=4)));
  annotation (experiment(StopTime=100, Tolerance=1e-006));
end ThreeSpecies_Quiescent;

需要做的仅仅是以ThirdSpecies模型为基础进行扩展,然后在包含的每个种群模型里修改init参数。对模型进行仿真就可以得到每个物种在平衡态的种群大小

../../../_images/ThreeS.png

从生态学的角度来看,我们已可以对这个系统做出一个有趣的观察。如果我们从一个非平衡条件启动仿真,系统将迅速变为不稳定。换句话说,引入狼群对原本稳定、只涉及兔子狐狸的生态系统内的种群动态产生了显著影响。