一类捕食者-食饵模型的敏感性分析和最优控制

Sensitivity Analysis and Optimal Control of a Class of Predator-Prey Model

摘要

对一类具有Allee效应和食饵收获的捕食者-食饵模型进行研究。选取拉丁超立方抽样法，通过多次抽样计算偏秩相关系数，分析影响种群数量的关键因素。研究了在维持系统稳定性的前提下获得最大收获问题，运用Pontryagin 极大值原理， 以收获系数作为控制参数，得到了最优收获策略。研究结果表明，食饵的出生率和收获强度是影响种群数量的主要因素。在合理的收获策略下，两种群的数量都有所增加，这对保持生态稳定性具有一定的参考价值。

Abstract

This paper studies a predator-prey model with Allee effect and prey harvesting. We select the Latin hypercube sampling method, calculate the partial rank correlation coefficient through multiple sampling, and analyze the key factors affecting population density. At the same time, we study the problem of obtaining the maximum harvesting under the premise of maintaining the stability of the system, using the Pontryagin principle of maximum value and using the harvesting coefficient as the control parameter to obtain the optimal harvesting strategy. The research results show that the birth rate and harvesting intensity of prey are the main factors affecting population density. Under a reasonable harvesting strategy, the number of both populations will increase, which has some reference value for maintaining ecological stability.

Graphical abstract

关键词

捕食者⁃食饵模型 / 敏感性分析 / 最优控制

Key words

predator-prey model / sensitivity analysis / optimal control

0　 引 言

自Lotka和Volterra在关于种群动力学的开创性工作以来，许多生态学家和应用数学家通过构建数学模型的方法来研究食饵种群和捕食者种群之间的关系^[1],这种类型的模型通常称为捕食者-食饵模型。对捕食者-食饵模型的研究有助于了解不同生态效应对相互作用的种群的具体影响。由于这种模型在生态学中具有重要意义，它将继续作为生态领域的主题之一被深入地研究下去^[2~4]。

种群动力学的经典理论预言，种群越小(或其密度越低)，遭受种内竞争约束的个体就越少;个体越少，可用资源就越多，种群稳定增长。尽管表面上这种增加的资源可用性使种群受益，但该种群将由于缺乏同种物种，丧失了同种物种数量较多时所具备的优势(其优势包括合作捕食或防御入侵以及社会温度调节)^[5]，这可能导致种群的健康状况下降，这种现象称为Allee效应。Allee效应会导致平均增长率不是单调的，甚至在种群密度特别低的情况下(强Allee效应)会变为负值。Allee效应具有临界阈值，低于该阈值，就有可能种群灭绝。出现Allee效应的机制很丰富，包括近交衰退，无法抵御天敌，以及两性生殖成熟的时间异步等^[5]。

另一方面，由于人类的生活、生产需求对渔业、林业和野生动植物资源进行开发和收获，势必对生态系统产生巨大影响。数学模型研究结果表明，收获对种群的动力学行为演变有很大的影响^[6]。考虑到Allee效应在动植物种群中的广泛存在和收获对人类生活生产的重要性，将Allee效应和收获考虑在内一定会增强模型的实际意义。因此，首先考虑广义的Rosenzweig和MacArthur模型

\left\{\begin{array}{l}\frac{\mathrm{d}u}{\mathrm{d}t}=au(\mathrm{1}-\frac{u}{k})-\frac{buv}{c+u}\\ \frac{\mathrm{d}v}{\mathrm{d}t}=\frac{muv}{c+u}-h(v)\end{array}\right.（1）

其中，u和v分别代表食饵和捕食者的种群密度，a是食饵出生率，k是环境容纳量，b是食饵的最大吸收率，c是捕食者具有最大捕杀率的食饵密度，m是捕食者出生率，函数h(v)是捕食者的死亡率。进行无量纲化，u\to kU,v\to \frac{ac}{b}V,\frac{k}{c}\to \beta ,\frac{a}{m}\to \alpha ,t\to \frac{T}{m}，得到^[7]

\left\{\begin{array}{l}\frac{\mathrm{d}u}{\mathrm{d}t}=\alpha u(\mathrm{1}-u)-\frac{\alpha uv}{\mathrm{1}+\beta u}\\ \frac{\mathrm{d}v}{\mathrm{d}t}=\frac{\beta uv}{\mathrm{1}+\beta u}-h(v)\end{array}\right.（2）

引入弱Allee效应和线性食饵收获得到下面的模型

\left\{\begin{array}{l}\frac{\mathrm{d}u}{\mathrm{d}t}=\alpha u(\mathrm{1}-u)(u+A)-\frac{\alpha uv}{\mathrm{1}+\beta u}-qEu\\ \frac{\mathrm{d}v}{\mathrm{d}t}=\frac{\beta uv}{\mathrm{1}+\beta u}-h(v)\end{array}\right.（3）

这里考虑h(v)=\frac{\gamma v^{\mathrm{2}}}{e+\eta v}的双曲死亡率情况，e和\eta是浮游生物死亡情况下水和自阴影造成的光衰减系数，\gamma是捕食者的死亡率，(u+A)是弱Allee效应项^[8]，q是捕捉系数，通过设置q=\mathrm{1}使收获项单位化，E是食饵收获强度。参数取值分别为：\alpha =\mathrm{0.65}，\beta =\mathrm{6}，\gamma =\mathrm{1.22}，e=\mathrm{1}，\eta =\mathrm{1.22}，A=\mathrm{0.2}，E=\mathrm{0.05}时,我们给出模型（3）唯一渐进稳定正平衡点(0.037 94,0.186 6)的相图(图1）且本文后续模拟的参数取值均固定为这组参数。

通常学者们对于捕食者-食饵模型的研究关注点在稳定性分析、分支计算等方面。本文重点关注参数变化对种群数量的影响，探讨决定种群数量的关键因素。拟通过拉丁超立方体抽样(Latin hypercube sampling,LHS)对参数进行抽样，计算偏秩等级相关系数(partial rank correlation coefficient,PRCC)，找到对种群数量影响较为明显的参数。人类对自然资源的开发过程中，经济发展与生态平衡之间总是存在着矛盾，不破坏生态系统稳定性的前提下最大限度的收获是本文的一个重要的研究内容，拟将食饵收获强度作为控制变量进行最优控制研究。

1　 种群的敏感性因素分析

拉丁超立方体抽样也称为无重复层级抽样，是蒙特卡罗(Monte Carlo，MC)抽样法中的一种，于1979年首次由Mckay提出^[9]。拉丁超立方体抽样的优点是避免了抽样的聚集现象，与其他随机抽样方法相比其抽样迭代次数较少。在执行参数抽样时选择均匀分布作为先验分布(均匀分布是无信息分布)，以捕食者种群v为例，抽样次数取2 000次，各参数取值范围见表1。\mathrm{P}\mathrm{R}\mathrm{C}\mathrm{C}值的计算方法参考文献[10]。请注意，当P值大于0.05时，则认为该参数对捕食者种群v的影响不显著，\mathrm{P}\mathrm{R}\mathrm{C}\mathrm{C}的绝对值越大就表示该参数对种群v的影响越明显。

由表1可知，参数\alpha和\gamma对捕食者种群的影响为正影响，其他参数为负影响。以影响较为明显的参数食饵出生率\alpha和收获强度参数E为例进行说明。食饵出生率增大会使食饵种群数量增大，由于捕食者需要捕食食饵才能稳定种群数量生存下去，食饵出生率间接影响了捕食者的数量，且为正影响。对食饵的收获变相等于和捕食者竞争食物，提高收获强度会使捕食者食物匮乏，所以食饵的收获强度参数E对捕食者种群的影响为负影响。图2给出了各个参数的PRCC散点图，其中横坐标表示所研究参数的秩变换值与其他参数秩变换值之间的线性回归残差，纵坐标表示输出的秩变换值与所有研究参数的秩变换值之间的线性回归残差。根据文献[10]，两个残差的散点图可以表征输入和输出之间是否存在显著的线性关系，同时散点图的倾斜程度(斜率)越大也对应PRCC的绝对值越大，散点图呈现递增状态时对应的PRCC值为正，呈递减状态时对应 PRCC 值为负。在同样的输入和输出下与传统的相关性指标的相较，偏秩相关系数（PRCC）可以使得参数和输出变化之间的线性关系变得更加明显。图3研究了连续时间段0到50内各参数的时变敏感性，曲线越靠近0代表着参数的变化对变量的影响越小。由图3可见，参数\alpha和E明显对变量有显著影响且分别为正影响和负影响。

2　 收获率的优化控制

在具有收获项的捕食者-食饵模型中，以收获强度参数作为控制变量，运用Pontryagin极大值原理进行优化，通过该方式得到的最优控制策略能在维持系统稳定和利益最大化之间得到平衡。本文运用Biswas等人^[11]的方法建立最优控制问题。

将E记作x来表示u(t)的收获率，因生态系统中的食饵数量和收获的强度都是有限的，结合实际情况可以得到下面的不等式

其中,\Omega (t)是单位时间所能收获的食饵数量随时间变化的上限;因食饵遭受Allee效应，捕食者的捕杀和人类的收获，我们引入约束条件使种群数量保持在一个较低水平。令X={[u,v]}^{\mathrm{T}}，收获强度x(t)作为控制变量，考虑到食饵种群数量和收获强度，把控制变量和状态变量的二次函数作为目标函数，可以得到关于收获强度的最优控制问题(P)

\mathrm{m}\mathrm{i}\mathrm{n}J={\int }_{\mathrm{0}}^{t_f}(B{u}^{\mathrm{2}}+x^{\mathrm{2}})\mathrm{d}t

\mathrm{s}.\mathrm{t}.

\dot{u}=\alpha u(\mathrm{1}-u)(u+A)-\frac{\alpha uv}{\mathrm{1}+\beta u}-xu

\dot{v}=\frac{\beta uv}{\mathrm{1}+\beta u}-\frac{\gamma v^{\mathrm{2}}}{e+\eta v}

u(\mathrm{0})=u_{\mathrm{0}},v(\mathrm{0})=v_{\mathrm{0}}

其中，B\in {\mathbb{R}}^{+}为食饵种群的权重系数，t_f\in {\mathbb{R}}^{+}是最终时间。问题(P)的最优解记为x_{\mathrm{*}}。利用下面的非负参变量{\delta }_i(i=\mathrm{1,2},\mathrm{3,4})将不等式约束转变为等式约束以便利用Pontryagin极大值原理考虑问题(P) ^[12]

\left\{\begin{array}{l}-x+{\delta }_{\mathrm{1}}=\mathrm{0}\\ x-x_{\mathrm{m}\mathrm{a}\mathrm{x}}+{\delta }_{\mathrm{2}}=\mathrm{0}\\ u-u_{\mathrm{m}\mathrm{a}\mathrm{x}}+{\delta }_{\mathrm{3}}=\mathrm{0}\\ xu-\Omega +{\delta }_{\mathrm{4}}=\mathrm{0}\end{array}\right.

令\delta =[{\delta }_{\mathrm{1}},{\delta }_{\mathrm{2}},{\delta }_{\mathrm{3}},{\delta }_{\mathrm{4}}{]}^{\mathrm{T}}，得到(P)的Hamiltonian函数

\begin{array}{l}H=B{u}^{\mathrm{2}}+x^{\mathrm{2}}+{\phi }_u(\alpha u(\mathrm{1}-u)(u+A)-\\ \frac{\alpha uv}{\mathrm{1}+\beta u}-xu)+{\phi }_v(\frac{\beta uv}{\mathrm{1}+\beta u}-\frac{\gamma v^{\mathrm{2}}}{e+\eta v})+\\ {\varphi }_{\mathrm{1}}(-x+{\delta }_{\mathrm{1}})+{\varphi }_{\mathrm{2}}(x-x_{\mathrm{m}\mathrm{a}\mathrm{x}}+{\delta }_{\mathrm{2}})+\\ {\varphi }_{\mathrm{3}}(u-u_{\mathrm{m}\mathrm{a}\mathrm{x}}+{\delta }_{\mathrm{3}})+{\varphi }_{\mathrm{4}}(xu-\Omega +{\delta }_{\mathrm{4}})\end{array}

\mathit{\phi }=[{\phi }_u,{\phi }_v{]}^{\mathrm{T}}是伴随变量，\mathit{\varphi }=[{\varphi }_{\mathrm{1}},{\varphi }_{\mathrm{2}},{\varphi }_{\mathrm{3}},{\varphi }_{\mathrm{4}}{]}^{\mathrm{T}}是非负参变量。根据Pontryagin极大值原理可以得到一阶必要条件

\dot{x}=\frac{\partial H}{\partial \mathit{\phi }}，\dot{\phi }=-\frac{\partial H}{\partial X}

其中，\mathit{\delta }\ge \mathrm{0}，\mathit{\varphi }\ge \mathrm{0}，{\mathit{\delta }}^{\mathrm{T}}\mathit{\varphi }=\mathrm{0}。

控制变量的最优条件为

\frac{\partial H}{\partial x}=\mathrm{0}

则由一阶必要条件可得下列伴随系统

\left\{\begin{array}{l}\begin{array}{l}{\dot{\phi }}_u=\frac{({\phi }_u\alpha -{\phi }_v\beta )v_{\mathrm{*}}}{(\mathrm{1}+\beta u_{\mathrm{*}}{)}^{\mathrm{2}}}+({\phi }_u-{\varphi }_{\mathrm{4}}x)-\\ {\phi }_u\alpha (A-u_{\mathrm{*}}^{\mathrm{2}}-\mathrm{2}A{u}_{\mathrm{*}})-{\varphi }_{\mathrm{3}}-\mathrm{2}A{u}_{\mathrm{*}}\end{array}\\ {\dot{\phi }}_v=\frac{({\phi }_u\alpha -{\phi }_v\beta )u_{\mathrm{*}}}{\mathrm{1}+\beta u_{\mathrm{*}}}+\frac{\mathrm{2}{\phi }_v\gamma e{v}_{\mathrm{*}}}{(e+\eta v_{\mathrm{*}}{)}^{\mathrm{2}}}\end{array}\right.

横截条件为{\phi }_u(t_f)={\phi }_v(t_f)=\mathrm{0}。

由\mathrm{0}\le x\le x_{\mathrm{m}\mathrm{a}\mathrm{x}}和xu\le \Omega，可得

\mathrm{0}\le x\le \mathrm{m}\mathrm{i}\mathrm{n}(x_{\mathrm{m}\mathrm{a}\mathrm{x}},\frac{\Omega }{u})

令u_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{1},又因为\Omega \le u_{\mathrm{m}\mathrm{a}\mathrm{x}},可得\frac{\Omega }{u}\le x_{\mathrm{m}\mathrm{a}\mathrm{x}}。

通过解最优条件可得

x_{\mathrm{*}}=\frac{({\phi }_u-{\varphi }_{\mathrm{4}})u_{\mathrm{*}}+{\varphi }_{\mathrm{1}}-{\varphi }_{\mathrm{2}}}{\mathrm{2}}

首先考虑\frac{\Omega }{u}>x_{\mathrm{m}\mathrm{a}\mathrm{x}}的情况。此时成立，所以{\varphi }_{\mathrm{4}}=\mathrm{0}\mathrm{。}为了使最优控制表达式为不含{\varphi }_{\mathrm{1}}和{\varphi }_{\mathrm{2}}的显式表达式，考虑下列三种假设：

1)　当时，{\varphi }_{\mathrm{1}}={\varphi }_{\mathrm{2}}=\mathrm{0},这时x_{\mathrm{*}}=\frac{{\phi }_u{u}_{\mathrm{*}}}{\mathrm{2}}；

2)　当x_{\mathrm{*}}=\mathrm{0}时，{\varphi }_{\mathrm{2}}=\mathrm{0},所以可得\mathrm{0}=x_{\mathrm{*}}=\frac{{\phi }_u{u}_{\mathrm{*}}+{\varphi }_{\mathrm{1}}}{\mathrm{2}},因为{\varphi }_{\mathrm{1}}\ge \mathrm{0},所以有{\phi }_u{u}_{\mathrm{*}}\ge \mathrm{0}；

3)　当x_{\mathrm{*}}=x_{\mathrm{m}\mathrm{a}\mathrm{x}}时，有{\varphi }_{\mathrm{1}}=\mathrm{0},可得x_{\mathrm{*}}=x_{\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{{\phi }_u{u}_{\mathrm{*}}-{\varphi }_{\mathrm{2}}}{\mathrm{2}},此外{\varphi }_{\mathrm{2}}\ge \mathrm{0},所以存在\frac{{\phi }_u{u}_{\mathrm{*}}}{\mathrm{2}}\ge u_{\mathrm{m}\mathrm{a}\mathrm{x}}\mathrm{。}

结合以上情况，有

x_{\mathrm{*}}=\mathrm{m}\mathrm{a}\mathrm{x}\{\mathrm{0},\mathrm{m}\mathrm{i}\mathrm{n}\{\frac{{\phi }_u{u}_{\mathrm{*}}}{\mathrm{2}},x_{\mathrm{m}\mathrm{a}\mathrm{x}}\}\}

同时考虑\frac{\Omega }{u}\le x_{\mathrm{m}\mathrm{a}\mathrm{x}}的情况可得

x_{\mathrm{*}}=\mathrm{m}\mathrm{a}\mathrm{x}\{\mathrm{0},\mathrm{m}\mathrm{i}\mathrm{n}\{\frac{{\phi }_u{u}_{\mathrm{*}}}{\mathrm{2}},\frac{\Omega }{u_{\mathrm{*}}}\}\}

综上所述，最优收获率的解析表达式为

x_{\mathrm{*}}=\mathrm{m}\mathrm{a}\mathrm{x}\{\mathrm{0},\mathrm{m}\mathrm{i}\mathrm{n}\{\frac{{\phi }_u{u}_{\mathrm{*}}}{\mathrm{2}},x_{\mathrm{m}\mathrm{a}\mathrm{x}},\frac{\Omega }{u_{\mathrm{*}}}\}\}

运用保辛伪谱方法(SPM)进行数值模拟^[13]，从生物背景出发,若以平凡平衡点和轴平衡点为初始值进行控制研究不具有现实意义。以正平衡点(0.037 94,0.186 6)作为食饵和捕食者种群的初始值,设定食饵的收获强度权重B=\mathrm{0.01},食饵种群数最大值u_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{3}\mathrm{000},收获强度最大值x_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{1},模拟最终时间t_f=\mathrm{20} a,单位时间收获食饵数量最大值\Omega =\mathrm{200},得到控制前后的对比图，如图4和图5。

由图4和图5可以看出，对收获进行控制后，食饵和捕食者两种群的数量都有所上升。通过控制收获的强度使得生态资源可持续发展，满足了人类的需求又促进了生态种群的增长。这符合我们的期望。

3　 结 语

本文研究了一类具有食饵收获和Allee效应的捕食者-食饵模型。首先，从各种生态因素对种群数量的影响大小出发，找到对种群数量影响的关键因素，通过敏感性分析研究参数对种群数量的影响，我们发现收获强度对种群数量的影响明显，且收获强度是一个可以人为控制干扰的参数。其次，我们将收获强度作为控制参数，通过计算和数值模拟得到了优化结果，即控制收获后种群数量呈现增加。本文的研究结果对于解决维持生态平衡下得到最优控制策略的问题提供了理论依据，对稀缺生物的保护，以及人类可持续发展具有一定的借鉴意义。

参考文献

[1] MURRAY J D. Mathematical Biology ［M］. 3rd Ed. Berlin： Springer， 2001.

[2] BAURMANN M， GROSS T， FEUDEL U. Instabilities in spatially extended predator-prey systems： Spatio-temporal patterns in the neighborhood of Turing-Hopf bifurcations ［J］. Journal of Theoretical Biology， 2007，245（2）：220-229. DOI： 10.1016/j.jtbi.2006.09.036 .

[3] 伏升茂， 苏发儒. 带恐惧因子和强Allee效应的捕食者-食饵扩散模型的Hopf分支［J］.西北师范大学学报（自然科学版）， 2019，55（3）：14-20.DOI：10.16783/j.cnki.nwnuz.2019.03.004 .

[4] FU S M， SU F R. Hopf bifurcation of a diffusive predator-prey model with fear factors and strong Allee effects［J］. Journal of Northwest Normal University（Natural Science）， 2019，55（3）：14-20. DOI：10.16783/j.cnki.nwnuz.2019.03.004（Ch ）.

[5] 王欢， 邢慧. 一类带有Allee效应及恐惧因子的捕食-食饵模型［J］.纺织高校基础科学学报， 2020，33（4）：85-90.DOI：10.13338/j.issn.1006-8341.2020.04.014 .

[6] WANG H， XING H. A class of predator-prey model with Allee effect and fear factor ［J］. Basic Sciences Journal of Textile Universities， 2020，33（4）：85-90. DOI：10.13338/j.issn.1006-8341.2020.04.014（Ch ）.

[7] STEPHENS P A， SUTHERLAND W J， FRECKLETON R P. What is the Allee effect ［J］. Oikos， 1999， 87（1）：185-190. DOI： 10.2307/3547011 .

[8] KAR T K， PAHARI U K ． Modelling and analysis of a prey-predator system with stage-structure and harvesting ［J］. Nonlinear Analysis： Real World Applications， 2007，8（2）：601-609. DOI： 10.1016/j.nonrwa.2006.01.004 .

[9] NAGANO S， MAEDA Y ． Phase transitions in predator-prey systems ［J］. Physical Review E， 2012， 85（1）：011915. DOI：10.1103/PhysRevE.85.011915 .

[10] BANERJEE M， TAKEUCHI Y. Maturation delay for the predators can enhance stable coexistence for a class of prey-predator models ［J］. Journal of Theoretical Biology， 2017，412：154-171. DOI：dx.doi.org/10.1016/j.jtbi.2016.10.016 .

[11] MCKAY M D， BECKMAN R J， CONOVER W J ． Comparison of three methods for selecting values of input variables in the analysis of output from a computer code［J］. Technometrics， 1979，21（2）： 239-245. DOI：dx.doi.org/10.1080/00401706.1979.10489755 .

[12] MARINO S， HOGUE I B， RAY C J， et al ．A methodology for performing global uncertainty and sensitivity analysis in systems biology ［J］. Journal of Theoretical Biology， 2008， 254（1）： 178-196. DOI：10.1016/j.jtbi.2008.04.011 .

[13] BISWAS M H A， PAIVA L T， DE PINHO M D ． A SEIR model for control of infectious diseases with constraints［J］. Mathematical Biosciences & Engineering， 2014，11（4）： 761-784. DOI： 10.3934/mbe.2014.11.761 .

[14] WANG X W， PENG H J， SHI B Y， et al ．Optimal vaccination strategy of a constrained time-varying SEIR epidemic model［J］. Communications in Nonlinear Science and Numerical Simulation， 2018， 67： 37-48. DOI： 10.1016/j.cnsns.2018.07.003 .

[15] PENG H J， WANG X W， SHI B Y， et al ． Stabilizing constrained chaotic system using a symplectic psuedospectral method［J］. Communications in Nonlinear Science and Numerical Simulation， 2019， 56： 77-92. DOI： 10.1016/j.cnsns.2017.07.028 .

基金资助

甘肃省科技计划项目(20JR5RA512)

教育部人文社会科学研究基金项目(20XJAZH006)

中央高校基本科研业务费专项资金资助(31920220066)

国家民委领军人才支持计划项目