猎食者猎物系统¶

目前为止，我们已经介绍了热学、电学和机械方面的实例。实际上，这些都是工程实例。然而，Modelica语言并不仅限于工程科学。为了强化这一点，本节将介绍基于猎食者和猎物之间关系的共同生态系统动力学模型。我们将使用[Lotka]-[Volterra]方程去建立上述模型。

经典猎食者猎物系统¶

经典的猎食者猎物模型涉及两个物种。一个物种叫做“猎物”。本节中，其物种种群用\(x\)表示。另一个物种称为“猎食者”，其物种种群用\(y\)表示。

有三个因素对猎食者猎物系统影响较大。首先是“猎物”物种的繁殖率。假设其繁殖率和种群成正比。如果你对化学反应比较熟悉，质量作用定律 [Guldberg] 和其在概念上是完全相同的。如果你对质量作用定律不了解，可以理解为生存环境中潜在的配偶越多，种群的繁殖率就会越高。我们可以用数学公式表达如下：

\[\dot{x}_{r} = \alpha x\]

其中，\(x\)表示猎物的种群，\(\alpha\)表示种群繁殖率的比例常数，\(\dot{x}_r\)表示由于繁殖引起的种群变化。

其次影响猎食者猎物系统的因素是猎食者的饥饿程度。如果没有足够的“猎物”来充饥，部分猎食者就会饿死。在建模饥饿程度时，要考虑竞争关系对其的重要影响。我们也是用比例关系描述上述模型。不同于繁殖率模型，此次是反比例关系，因为越多的猎食者越容易引发饥饿。在数学表达式上，这和繁殖率模型具有相同的表达形式：

\[\dot{y}_{s} = -\gamma y\]

其中，\(y\)表示猎食者种群。\(\gamma\)表示种群饥饿的比例常数。\(\dot{y}_s\)表示由于饥饿引起的猎食者种群变化。

最后一个我们要考虑的因素是“捕食”，即猎食者对猎物的消耗。如果没有天敌，猎物会指数增长（至少从数学观点）。因此，捕食是保持猎物种群的重要因素。同样的，没有猎物，猎食者也会全部死光。因此，捕食平衡了这种效应，并且防止猎食者种群的消失。再次的，我们得到一个比例关系。但是，这个比例关系实际上是一个双线性关系。而且，概念上也类似于质量作用定律。这种比例关系在数学上很容易得到。事实上，猎食者发现并捕获猎物的机会与猎物种群、猎食者种群两者都存在着正比关系。由于这种特殊的影响因素，需要两个物种都参与。因此，这里的数学表达式与前述的繁殖率及猎食者饥饿程度的公式都在结构上有所不同，即：

\[\begin{split}\dot{x}_p &= -\beta x y \\ \dot{y}_p &= \delta x y\end{split}\]

其中，\(\dot{x}_p\)表示由于捕食而导致的猎物种群减少。\(\dot{y}_p\)表示由于捕食而导致猎食者种群的增加。\(\beta\)表示捕获猎物概率的比例常数。\(\delta\)表示猎食者由于足够的营养以提高繁殖率可能性的比例常数。

全面考虑各种因素，可以得到每个种群的整体变化。用以下两个方程表示：

\[\begin{split}\dot{x} &= \dot{x}_r + \dot{x}_p \\ \dot{y} &= \dot{y}_p + \dot{y}_s\end{split}\]

根据上述两个方程，通过数学运算，可以重新整合方程等式的右侧，得到如下方程组：

\[\begin{split}\dot{x} &= x (\alpha - \beta y) \\ \dot{y} &= y (\delta x - \gamma)\end{split}\]

根据前几章我们学到的知识，将上述方程转化成Modelica语言应该相当简单：

model ClassicModel "This is the typical equation-oriented model"
  parameter Real alpha=0.1 "Reproduction rate of prey";
  parameter Real beta=0.02 "Mortality rate of predator per prey";
  parameter Real gamma=0.4 "Mortality rate of predator";
  parameter Real delta=0.02 "Reproduction rate of predator per prey";
  parameter Real x0=10 "Start value of prey population";
  parameter Real y0=10 "Start value of predator population";
  Real x(start=x0) "Prey population";
  Real y(start=y0) "Predator population";
equation
  der(x) = x*(alpha-beta*y);
  der(y) = y*(delta*x-gamma);
end ClassicModel;

在上述代码中，只有一点我们还没有讨论，即变量x和y的start属性。正如我们在前面引入物理章节里对NewtonCoolingWithUnits示例所介绍的那样，我们可以为变量指定各种的属性（在后面的章节内建类型中将详尽地讨论变量的属性）。我们已经在前面章节讨论了变量的unit属性。这是我们第一次看到start属性。

细心的读者可能已经注意到程序中对变量x0和y0的声明以及其所代表的种群初始值。根据前面实例的讲解，有人可能认为这些初始化条件是以如下的方式被模型获取的：

model ClassicModelInitialEquations "This is the typical equation-oriented model"
  parameter Real alpha=0.1 "Reproduction rate of prey";
  parameter Real beta=0.02 "Mortality rate of predator per prey";
  parameter Real gamma=0.4 "Mortality rate of predator";
  parameter Real delta=0.02 "Reproduction rate of predator per prey";
  parameter Real x0=10 "Initial prey population";
  parameter Real y0=10 "Initial predator population";
  Real x(start=x0) "Prey population";
  Real y(start=y0) "Predator population";
initial equation
  x = x0;
  y = y0;
equation
  der(x) = x*(alpha-beta*y);
  der(y) = y*(delta*x-gamma);
end ClassicModelInitialEquations;

然而，对于ClassicModel示例，我们采用了更简洁的方法。我们可以在声明变量的start属性时直接指定其初始条件。在初始化章节中将对这一功能做简单的介绍

值得注意的是，这种方法既有优点也有缺点。方法的优点是比较灵活。其实，start属性更多的是一种绑定关系的暗示。如果Modelica编译器将某个特定变量识别为状态（即某个变量需要一个初始化条件），并且模型的initial equation部分没有提供非充足的初始条件，这时start属性就可以作为替代，声明变量的初始条件。换句话说，你可以认为start属性是在必须具备初始条件的情况下的“后备初始条件”。

你需要注意start属性所具备的一些缺点。首先，此属性只是一个完全可以忽略的提示工具。其次，其是否被忽略也很难预测。因为哪些变量作为状态在不同的工具里选择也会不同。

避免上述缺点我们可以使用fixed属性（在内建类型章节中作进一步讨论）。变量的fixed属性被用来通知编译器start属性必须作为初始条件来使用。换句话说，如下的initial equation部分：

  Real x;
initial equation
  x = 5;

等同于下述变量声明中关于fixed和start属性的使用：

Real x(start=5, fixed=true);

最后，start属性的另一个复杂之处是其功能被“重载”了。这就意味着，此属性有两种不同的用途。如果，所讨论的变量不是一个状态，而是一个“迭代变量”（即变量的解依赖于非线性方程组），那么start属性可以通过编译器作为初始的猜测值（即非线性求解方程首次迭代时使用的变量值）。

是否指定一个变量的start属性取决于你希望初始条件会得到多严格的。决定这点对于需要对语言有一定的经验。而这已经超出了本书的范围。然而，至少值得在这里指出不同选择的存在，以及简单解释一下权衡要素。

不管使用哪种初始化的方法，这些模型的仿真结果都是相同的。猎食者猎物系统的典型结果如下图所示：

通过上图，可以看到每个种群的周期性行为。最初，没有充足的食物来支撑较多的猎食者。存在的猎食者捕食任何可以找到的猎物。但即便如此，饥饿还是会发生并且导致猎食者种群数量的下降。在此期间，猎食者消耗猎物种群的速率如此之高，以致于猎物种群的繁殖速率不足以弥补因捕食而减少的种群数。因此，猎物种群的数量也在下降。

在某些点上，猎食者种群的数量变的如此之低，导致猎物种群的繁殖数量高于因捕食而减少的种群数。猎物物种的数量开始反弹。因为猎食者种群大小的回升需要较长的时间。所以在此期间，猎物种群的增长速度几乎不受猎食者种群的影响。最终，由于猎物的丰盛，猎食者种群的数量开始回升，之至系统恢复到猎食者和猎物种群的原始状态，然后整个循环再重演，无穷无尽。

该系统一次次的返回到相同的初始条件处（当然要忽略数值计算错误）。这一现象是系统最有趣的地方之一。特别是考虑猎食者猎物系统方程实际上是非线性的，这也就更加令人注目了。

稳定状态初始化¶

让我们想象一下，这些物种种群的极端波动可能会导致一些不良的生态后果。在这种情况下，理解减少或消除这些波动的方法对整改系统将会很有帮助。一个简单的方法就是使两个种群维持在相对平衡的状态。但是，如何使用这些模型来帮助我们确定这样的“稳定”状态？

上述问题的答案就在于初始化条件。我们可以用能描述系统处于平衡状态的方程来初始化系统（相应的设置方法可以从前面讲述的实例FirstOrderSteady中获得），而不是直接给猎食者和猎物种群指定初始值。幸运的是，Modelica包含足够丰富的初始化方法，允许我们指定上述（或者其他）类型的初始化条件。

为了确保我们的系统开始的时候处于平衡状态，我们只需要简单的定义什么是平衡。从数学描述上来说，系统如果满足以下两个条件，就可以说该系统是处于平衡状态的：

\[\begin{split}\dot{x} &= 0 \\ \dot{y} &= 0\end{split}\]

为了在Modelica模型中获取上述特征，我们只需在initial equation区域增加以下方程，如下所示：

model QuiescentModel "Find steady state solutions to LotkaVolterra equations"
  parameter Real alpha=0.1 "Reproduction rate of prey";
  parameter Real beta=0.02 "Mortality rate of predator per prey";
  parameter Real gamma=0.4 "Mortality rate of predator";
  parameter Real delta=0.02 "Reproduction rate of predator per prey";
  Real x "Prey population";
  Real y "Predator population";
initial equation
  der(x) = 0;
  der(y) = 0;
equation
  der(x) = x*(alpha-beta*y);
  der(y) = y*(delta*x-gamma);
end QuiescentModel;

上述模型和之前模型最大的区别就是包含突出显示的初始化方程。回到上述模型，你可能想知道那些初始化方程到底意味着什么？毕竟，我们需要求解的是变量x和y。但是这些变量甚至都没有出现在我们的初始化方程内。系统是如何对这些变量求解的？

问题的答案在于理解函数\(x(t)\)和\(y(t)\)均是通过对带初始条件的微分方程进行积分来求解的。在仿真过程中，我们可以看到\(x\)和\(\dot{x}\)是通过下面的方程进行“耦合”的：

\[x(t) = \int_{t_0}^{t_f} \dot{x}\ \mathrm{d}x + x(t_0)\]

（当然了，在\(y\)和\(\dot{y}\)之间也存在着类似的关系）

但是，在系统初始化的过程中（即在计算初始条件时），上述关系是不成立的。在这种情况下，\(x\)和\(\dot{x}\)之间并不存在“耦合”关系（对\(y\):和\(\dot{y}\)并不适用）。换句话说，即便知道变量\(x\)和\(y\)是如何定义的也不能提供求解方程\(\dot{x}\)或\(\dot{y}\)的线索。在初始化问题上，我们可以认为\(x\)、\(y\)、\(\dot{x}\)和\(\dot{y}\)是四个相互独立的变量。

换一种说法，即在仿真的过程中，我们通过对\(\dot{x}\)进行积分来求解\(x\)。因此，积分方程是用于求解\(x\)的方程。但是在初始化过程中，我们不能用这个等式，所以（对每个在仿真过程中求解的微分方程），我们需要一个额外的方程。

任何情况下，我们要明白是在初始化过程中需要四个不同的方程来得到唯一解。在我们的QuiescentModel模型中，这四个方程如下所示：

\[\begin{split}\dot{x} &= 0 \\ \dot{x} &= x \ (\alpha - \beta y) \\ \dot{y} &= 0 \\ \dot{y} &= y \ (\delta x - \gamma)\end{split}\]

这些方程之间并不相互矛盾，理解这一点非常重要。如果你有编程背景，对前面两个方程可能会有疑问。“到底\(\dot{x}\)取值是多少？是零？还是\(x (\alpha - \beta y)\)？”问题的答案是都是。没有理由认为这两个方程不可能同时为真。

重要的是要记住，这些方程不具备赋值语句的功能。下面系统的方程与上述系统在数学表达上是完全一样的，并且清楚的展现了\(x\)和\(y\)是如何求解的：

\[\begin{split}\dot{x} &= 0 \\ \dot{y} &= 0 \\ x (\alpha - \beta y) &= \dot{x} \\ y (\delta x - \gamma) &= \dot{y}\end{split}\]

在这种格式下，我们更容易理解如何对\(x\)和\(y\)值进行求解。首先要注意的是，我们不能显式求解变量\(x\)和\(y\)的值。换句话说，如果没有变量\(x\)出现在等式的右边，我们就不能将这些方程的形式变换为\(x=...\)。所以我们不得不接受，这个方程组同时包含变量\(x\)和\(y\)。

但是，更加复杂的是该仿真系统是非线性的（这也是为什么我们不能用线性代数将其变换为显式方程组）。事实上，如果我们仔细研究这些方程，可以发现存在两种可能解。一个解是平凡解（\(x=0;y=0\)），另一个则是非零解。

如果我们试着仿真建立的QuiescentModel模型，会出现什么样的结果呢？仿真结果如下图所示：

根据第一种解，猎食者和猎物的种群数都变为了零。这种情况下系统没有繁殖、捕食以及饥饿。因为上述影响因素都和物种种群数成比例，而物种种群数都为零。所以系统没有变化。但这并不是一个很有趣的解决方案。

因为该系统是非线性的，因此系统方程有两个解。我们怎样才能使非线性求解器远离这个零根呢？如果你比较关注经典猎食者猎物系统模型的讨论，那么已经给你了答案的暗示。

还记得，前面提到过start属性是被重载了？在讨论经典猎食者猎物系统模型时，我们曾经指出，如果具有start属性的变量被选为迭代变量，start属性的其中一个作用是提供初始化猜想值。在我们的QuiescentModel模型中，恰巧变量x和y就是迭代变量。因为该变量必须通过系统的非线性方程组来求解。这也就意味着，我们要对变量x和y的start属性值进行指定，以尽量“避开”系统的零解（或者说至少接近我们期望的非零解）。例如：

model QuiescentModelUsingStart "Find steady state solutions to LotkaVolterra equations"
  parameter Real alpha=0.1 "Reproduction rate of prey";
  parameter Real beta=0.02 "Mortality rate of predator per prey";
  parameter Real gamma=0.4 "Mortality rate of predator";
  parameter Real delta=0.02 "Reproduction rate of predator per prey";
  Real x(start=10) "Prey population";
  Real y(start=10) "Predator population";
initial equation
  der(x) = 0;
  der(y) = 0;
equation
  der(x) = x*(alpha-beta*y);
  der(y) = y*(delta*x-gamma);
end QuiescentModelUsingStart;

上述模型引导我们设置一组初始条件，使得仿真结果更加符合我们最初的目标（即猎食者和猎物种群大小都得到了非零解）。

值得指出的是（在内建类型章节中将会讨论），start属性的默认值是零。这也就是为什么当我们仿真QuiescentModel模型时，我们会恰巧准确的得到了系统的零解。因为这是我们初始的猜测，并且恰巧也是系统的精确解，因此系统就无需进行迭代或求解其他根了。

避免重复¶

基于猎食者猎物方程，我们已经讨论过几个不同的模型（ClassicModel，QuiescentModel以及QuiescentModelUsingStart模型）。你有没有注意到这些模型间的共同点？如果你仔细观察，你会发现它们几乎所有的内容都相同。实际上模型之间几乎没有任何差别！

在软件工程中，有一种说法是 “冗余是一切罪恶的根源” 。这里的情况也不例外（其实很大程度上如此。因为Modelica代码其实也是属于软件）。目前为止，我们所写的代码维护起来将会特别恼人。这是因为我们发现的任何错误都必须在每个模型中修改。此外，我们所作的任何改进也必须应用到每一个模型上。目前，我们处理的模型数量相对较少，但这种“复制粘贴”的模型开发方法会导致大量只有轻微差异的模型。

那么，我们能做些什么来避免上述情况的发生呢？在面向对象的编程语言中，基本上存在两种机制来减少冗余代码。他们是组合（在后面的组件章节中进行讨论）和继承。这里我们只对继承进行简单的介绍。

如果我们仔细观察QuiescentModelUsingStart模型，会发现和原来的ClassicModel模型之间几乎没有差别。事实上，唯一的不同之处如下所示：

model ClassicModel "This is the typical equation-oriented model"
  parameter Real alpha=0.1 "Reproduction rate of prey";
  parameter Real beta=0.02 "Mortality rate of predator per prey";
  parameter Real gamma=0.4 "Mortality rate of predator";
  parameter Real delta=0.02 "Reproduction rate of predator per prey";
  parameter Real x0=10 "Start value of prey population";
  parameter Real y0=10 "Start value of predator population";
  Real x(start=x0) "Prey population";
  Real y(start=y0) "Predator population";
equation
  der(x) = x*(alpha-beta*y);
  der(y) = y*(delta*x-gamma);
end ClassicModel;

model QuiescentModelUsingStart "Find steady state solutions to LotkaVolterra equations"
  parameter Real alpha=0.1 "Reproduction rate of prey";
  parameter Real beta=0.02 "Mortality rate of predator per prey";
  parameter Real gamma=0.4 "Mortality rate of predator";
  parameter Real delta=0.02 "Reproduction rate of predator per prey";
  Real x(start=10) "Prey population";
  Real y(start=10) "Predator population";
initial equation
  der(x) = 0;
  der(y) = 0;
equation
  der(x) = x*(alpha-beta*y);
  der(y) = y*(delta*x-gamma);
end QuiescentModelUsingStart;

换句话说，唯一的区别是增加了initial equation区域（原始的ClassicModel模型已经包含了变量x和y的非零start属性值）。理想情况下，我们可以通过简单的定义一个模型与另外模型间的差异，就可以由此避免任何冗余代码。事实证明，这也正是extends关键字的作用。考虑QuiescentModelUsingStart模型的以下版本：

model QuiescentModelWithInheritance "Steady state model with inheritance"
  extends ClassicModel;
initial equation
  der(x) = 0;
  der(y) = 0;
end QuiescentModelWithInheritance;

注意extends关键字的使用。从概念上讲，“扩展子句”只是简单的要求编译器将另外的模型（本例是ClassicModel模型）插入到所定义的模型中。通过这种方式，我们将从ClassicModel模型中复制（或“继承”）其包含的所有内容，而无需重复定义。因此，除了新加入的初始化方程外，QuiescentModelWithInheritance模型和ClassicModel模型其他部分完全一样。

但是，有一种情况，如果我们不想完全精确地继承模型中的内容，该怎么办？例如，如果我们想要改变参数gamma和delta的值？

当使用extends功能时，Modelica语言允许我们对模型加入相应的“修改语句”。模型对应的修改语句紧跟在所继承模型的名字后面，如下所示：

model QuiescentModelWithModifications "Steady state model with modifications"
  extends QuiescentModelWithInheritance(gamma=0.3, delta=0.01);
end QuiescentModelWithModifications;

还要注意的是，我们可以从ClassicModel中继承。但是，随后我们还必须重新定义初始化方程组以获得静止的初始条件。但是，我们不是通过继承QuiescentModelWithModifications模型，而是通过重复使用两个不同模型的内容来避免重复。

更多的种群动态

这就是本章所包含的所有实例。如果你想更深入的探索猎食者猎物方程，在下面再探猎食者猎物模型章节中将会重新讲述如何通过示意图并连接在一起的图形化组件来构建复杂的种群动态模型。

[Lotka]

Lotka, A.J., “Contribution to the Theory of Periodic Reaction”, J. Phys. Chem., 14 (3), pp 271–274 (1910)

[Volterra]

Volterra, V., Variations and fluctuations of the number of individuals in animal species living together in Animal Ecology, Chapman, R.N. (ed), McGraw–Hill, (1931)

[Guldberg]

C.M. Guldberg and P. Waage,”Studies Concerning Affinity” C. M. Forhandlinger: Videnskabs-Selskabet i Christiana (1864), 35

Table Of Contents

猎食者猎物系统¶

经典猎食者猎物系统¶

稳定状态初始化¶

避免重复¶