基于动态价格下的家畜养殖模型及其分析

Livestock Breeding Model Based on Dynamic Price and Its Analysis

摘要: 在任何生物经济学研究中，市场价格都是最重要的工具之一。养殖户的每一步都取决于其市场价格。羊的市场价格与消费者的需求量，羊的市场供给量有密切的关系。本文以羊养殖为例，以Logistic增长模型为基础，根据市场价格变化对养殖户行为的影响建立了价格影响下的动态单种群养殖模型，利用定性理论证明了平衡点的存在性、局部稳定性、极限环的不存在性、出现的跨临界点分支，利用Poincar-Bendixson准则证明平衡点的全局稳定性。该模型与内蒙古地区每年年底羊的存栏量数据相吻合，通过重要参数变化下的分支现象进行现实解释，可以发现养殖成本、出栏率、市场需求量对当前养殖量的影响，反映了动态价格影响下，养殖量的变化过程，为养殖户的长期决策提供了一定的建议。

Abstract: In any Bioeconomic research, market price is one of the most important tools. Every step of farmers depends on their market price. The market price of sheep is closely related to the demand of con-sumers and the market supply of sheep. Taking sheep breeding as an example, based on the logistic growth model, a dynamic single population breeding model under the influence of market price change on Farmers’ behavior is established. The existence, local stability, nonexistence of limit cycle and cross critical point branches of the equilibrium point are proved by using qualitative theory, and the global stability of the equilibrium point is proved by using Poincar Bendixson criterion. The model is consistent with the data of sheep stock at the end of each year in Inner Mongolia. Through the realistic explanation of the branching phenomenon under the change of important parameters, we can find the impact of breeding cost, marketing rate and market demand on the current breed-ing quantity, reflect the change process of breeding quantity under the influence of dynamic price, and provide some suggestions for farmers’ long-term decision-making.

1. 引言

在自然界中存在各种生物资源，如海洋中的鱼群、养殖场中羊群等，如何将这些有限生物资源可持续的开发和利用，是当今经济管理学家和生态学者都关心的问题。动力学的方法对研究种群生态学的相关问题，特别是研究种群个体数量和结构的变化规律问题有着重要作用。

单种群生物模型是种群生态学中发展最早的数学模型，文献 [1] [2] [3] 给出了一些单种群增长模型，文献 [4] [5] [6] 描述了许多被开发的单种群生物模型，对模型进行了定性分析之后，确定了最大捕捞量和可持续量。文献 [7] [8] 研究了自治单种群捕获模型，以最大可持续捕获量为决策目标，得到了最优开发策略。但是捕获量不能直接影响经济效益，换句话说，捕获量大未必经济效益就好，还与劳动成本、市场价格等因素有关。文献 [9] 中Gordon H S假设生产的单位生物量价格和生产成本都是常数，并把价格和成本因素引入生物模型中，确定了最优生产是将收入和成本之间的差异最大化。Clark C W [10] 研究数学生物经济学，致力于求解一些单种群捕获模型的最优捕获策略。文献 [11] [12] 讨论了以最大可持续捕获量和最大净收益为决策目标，研究了单种群生物模型，得到最优捕获策略。

这些模型都给出了开放获取的概念，而忽视了市场自由概念，忽视了供求之间的平衡，价格是一个长期浮动变化的变量，因此Clark C W [13] 从经济学的角度出发，将价格作为变量处理，对此前的模型进行了补充。文献 [14] [15] 中对价格随供求变化的开放式单种群捕获模型进行研究，建立了一般生产函数，用微分方程定性与稳定性理论分析了模型的动力学行为，并从生物经济学的角度出发解释了相关结果。高海音 [16] 建立了具有Gompertz增长的单种群模型，建立了价格成本变化的单种群经济捕捞模型。他们的研究都是通过微分方程的定性理论得到单种群生物资源的最优收获策略，目前还有一些研究是通过最优控制理论得到最优捕获策略。比如，文献 [17] 利用最大化原理考虑了在可变价格与可变成本下，开放式单种群捕获模型的最优策略，并针对Logistic模型进行了具体讨论，得到了最优控制的隐式表达。文献 [18] 也利用最大化原理，讨论了养殖式单种群生物模型，在贴现率为0和价格为常数情况下，饲养场最优生产方案，但是这与市场规律不符合，在实际过程中，价格会随着供应的改变而改变。大多数的单种群模型都是开放式单种群捕捞模型，关于养殖式单种群捕捞模型的讨论很少，养殖式单种群模型具有一些特殊性，所以还需要很多工作。

在任何生物经济学研究中，市场价格都是最重要的工具之一。目前的研究虽然考虑了价格受市场供求的影响，但是市场的影响是双向的，反过来说，价格也会影响市场供应，进而影响资源的开发。养殖户的每一步都取决于其市场价格。羊的市场价格与消费者的需求量，羊的市场供给量有密切的关系。如何将其互相影响的关系建立起来，是本文的研究目的之一。

本文以羊养殖为例，以Logistic增长模型为基础，根据市场价格变化的动态影响建立了相应的改进的价格影响下的动态养殖单种群模型，研究了其内在机理。本文的研究结构如下：在第2节中给出了价格影响下的动态养殖模型及其参数的描述，第3节利用定性理论证明了平衡点的存在性、局部稳定性、出现的跨临界点分支、极限环的不存在性以及平衡点的全局稳定性。第4节得出模型的解与内蒙古地区年末羊存栏量数据进行拟合情况和一些重要参数的分支结果，并且给出了长期养殖状态。最后在第5节中讨论主要发现和进一步研究范围。

2. 模型的建立

定义 $N\left(t\right)$ 为时间t种群的数量。种群在没有开发状态下，按照自然状态下进行生长，可以用经典的logistic增长模型来描述 [19]：

$\frac{\text{d}N}{\text{d}t}=rN\left(1-\frac{N}{k}\right)$ (1)

该模型中k代表环境的最大承载量即代表着栖息地可以容纳的最大种群数量，是一个常数。

为了获得经济收益，养殖户会对羊群进行开发和利用，假定羊群的生产函数是由当前种群数 $N\left(t\right)$ 和出栏率 $\theta$ 决定，进一步可得，生产函数如下： $h\left(t\right)=\theta N(\; t\; )$

因此生物资源在受人们开发时，其增长曲线为：

$\frac{\text{d}N}{\text{d}t}=rN\left(1-\frac{N}{k}\right)-h\left(t\right)$ (2)

具体来说，对于羊养殖过程， $N\left(t\right)$ 为养殖场在t时刻羊的存栏量，如果生产的羊全部投入市场，那么生产函数 $h\left(t\right)$ 表示t时刻转化为羊商品投放到市场的羊数量，当前存栏量的增长的趋势与(2)式一致。系统

(2)的一个重要特征是，如果出栏率 $\theta >r$，羊群可能会灭绝，相反羊群的数量将会稳定到 $N^{*}=k\left(1-\frac{\theta }{r}\right)$ 其

中 $\theta ,r,k$ 都为正数。无开发下的种群模型系统(1)，羊群永远不会灭绝，在 $N^{*}=k$ 时，趋于稳定。

在实际的养殖过程中，羊的市场价格显然不是一成不变的，价格受供求关系的影响，在需求一定时，供给量越小价格越高，供给量越大价格越低。市场价格为关于市场投放量h的函数。本研究参考Clark C W的价格函数：

$p\left(h\right)=a\left(1-\frac{h}{b}\right)$ (3)

在这里，a表示消费者愿意为羊支付的最高价格，b表示以任何价格都可以出售的羊的最大数量可以反映

市场羊的需求量，a和b是独立的，都为正常数。可以发现， $\underset{h\to 0}{\lim }p\left(h\right)=a$， $\underset{h\to b}{\lim }p\left(h\right)=0$。当h趋近于0

值时，羊的市场价格趋于最大的市场价格a，h趋近于b值时，羊的市场价格趋近于0。在实际情况下，羊的供应和价格都不可能为0，故而市场投放量h的取值范围为(0, b)。

对养殖户来说，羊的市场价格变化对他们有一定的影响，当市场价格高时，养殖的收益增多，会吸引更多人投资养羊，相应的扩大养殖场，以增加当前养殖量。反之，当市场价格低时，养殖羊不收益，或者收益不足，那么现存的养殖户，会离开这个行业，进而很多养殖场关闭，养殖规模缩小。由此，发现了羊商品的市场价格对养殖规模的影响过程。市场价格对当前养殖量的动态影响如图1所示。

Figure 1. Dynamic diagram of sheep breeding process

图1. 羊养殖过程动态图

如果将羊市场价格的高低用一个价格阈值 $p_0$ 来描述，由于价格阈值不可能超过市场的最大价格，故而 $p_0\in \left(0,a\right)$。根据上述分析某一特定时间的最大养殖量随着该时间段羊市场价格和价格阈值 $p_0$ 之差而变化，比例常数 $\epsilon >0$。我们都知道，受自然资源条件的限制，最大养殖量不可能无限制的扩张，把最大养殖量的限制定义为 $k_m$。基于这些假设，建立了价格影响下的最大养殖量变化曲线：

$\frac{\text{d}k}{\text{d}t}=\epsilon \left(p-p_0\right)\left(1-\frac{k}{k_m}\right)$ (4)

在模型(3)中，最大养殖量被认为是一个参数，模型(4)可以得到最大养殖量不是一个常数，而是一个随价格变化的时变变量，由此可以得到市场价格动态影响下的养殖模型可以表示为：

$\{\begin{array}{l}\frac{\text{d}N}{\text{d}t}=rN\left(1-\frac{N}{k}\right)-\theta N\hfill \\ \frac{\text{d}k}{\text{d}t}=\epsilon \left(p-p_0\right)\left(1-\frac{k}{k_m}\right)\hfill \end{array}$ (5)

将价格函数(3)代入系统(5)，我们得到了以下模型：

$\{\begin{array}{l}\frac{\text{d}N}{\text{d}t}=rN\left(1-\frac{N}{k}\right)-\theta N\hfill \\ \frac{\text{d}k}{\text{d}t}=\epsilon \left(a\left(1-\frac{\theta N\left(t\right)}{b}\right)-p_0\right)\left(1-\frac{k}{k_m}\right)\hfill \end{array}$ (6)

系统(6)初始条件 $\left(N_0,k_0\right)\in R^{+}$，由价格函数的定义及其模型假设，模型(6)中使用的变量参数解释和取值范围如表1所示。显而易见，模型(6)是logistic模型的推广，例如，如果 $\epsilon =0$，则该模型就是一个logistic模型。

Table 1. Description of parameters in model (6)

表1. 模型(6)中参数的描述

3. 模型的动力学分析

3.1. 解的正性和有界性

在本节中，我们首先试图找到系统(6)的解是正的和有界的条件。

定理3.1：如果k(t)总是非负的，且系统(6)的初始条件N(0)为正的，则系统(6)的所有解都为正的，并且有界。系统(6)正不变集和可行域为 $\Gamma =\left\{\left(N,k\right)\in R_2^{+}|0<N\le k_m,0<k\le k_m\right\}$。

证明：1) 正性：系统(6)的第一个方程可以改写为

$\frac{\text{d}}{\text{d}t}\left[N\left(t\right){\rho }_1\left(t\right)\right]=0$ (7)

这里 ${\rho }_{\text{1}}\left(t\right)={\text{e}}^{{\displaystyle {\int }_0^t\left[-r\left(1-\frac{N\left(s\right)}{k\left(s\right)}\right)+\theta \right]\text{d}s}}>0$

因此，

$N\left(t\right){\rho }_1\left(t\right)-N\left(0\right)=0$ (8)

$N\left(t\right)={\rho }_1^{-1}\left(t\right)N\left(0\right)>0$ (9)

根据假设 $\forall t,k\left(t\right)>0$，可以得出结论，系统(6)的所有解都是正的。

2) 有界性：

从系统(6)的第二个方程可知，

$\frac{\text{d}k}{\text{d}t}=\epsilon \left(a\left(1-\frac{\theta N}{b}\right)-p_0\right)\left(1-\frac{k}{k_m}\right)\le ϵ\left(a-p_0\right)+\frac{ϵ\left(p_0-a\right)k}{k_m},$ (10)

所以有 $\underset{h\to \infty }{\lim }\sup k\left(t\right)\le k_m$。

由系统(6)的第一个方程显而易见， $N\le k_m$

综上所示，可得系统(6)的可行域为 $\Gamma =\left\{\left(N,k\right)\in R_2^{+}|0<N\le k_m,0<k\le k_m\right\}$。

3.2. 平衡点的存在性

接下来我们将讨论系统(6)中参数在定义范围下平衡点的存在性，显然，系统(6)总有一个零平衡点 $E_1=\left(0,k_m\right)$，系统(6)的正平衡点可通过解下列代数方程得到

将价格函数(3)代入系统(5)，我们得到了以下模型：

$\{\begin{array}{l}\frac{\text{d}N}{\text{d}t}=rN\left(1-\frac{N}{k}\right)-\theta N=0\hfill \\ \frac{\text{d}k}{\text{d}t}=\epsilon \left(a\left(1-\frac{\theta N\left(t\right)}{b}\right)-p_0\right)\left(1-\frac{k}{k_m}\right)=0\hfill \end{array}$ (11)

对方程(11)进行转换可得，

$\left(k-k_m\right)\left(k-\frac{br\left(a-p_0\right)}{a\theta \left(r-\theta \right)}\right)=0$ (12)

易知(12)式的两个根为 $k_1=k_m$， $k_2=\frac{br\left(a-p_0\right)}{a\theta \left(r-\theta \right)}$

对应代入(11)式第一个方程，对应求得

$N_1=\frac{\left(r-\theta \right)k_m}{r},N_2=\frac{b\left(a-p_0\right)}{a\theta }$

然后可得系统(6)的两个正平衡点 $E_2=\left(N_1,k_1\right)$， $E_3=\left(N_2,k_2\right)$。我们知道当(12)式的判别式 $\Delta =B_1^2-4B_0{B}_1={\left[\epsilon a\theta k_m\left(r-\theta \right)+r\epsilon b\left(a-p_0\right)\right]}^2=0$，方程有两个相同的根，故而系统(6)有两个相同的正平衡点。

接下来为了便于下面讨论，现在定义

$\begin{array}{l}{D}_1=\left\{\left(\theta ,p_0\right)|0<p_0<\psi \left(\theta \right)\right\},\\ D_2=\left\{\left(\theta ,p_0\right)|\psi \left(\theta \right)<p_0<a\right\},\\ C_1=\left\{\left(\theta ,p_0\right)|p_0=\psi \left(\theta \right)\right\},\end{array}$

此时， $\psi \left(\theta \right)\triangleq \frac{k_{m}a}{rb}{\theta }^2-\frac{k_{m}a}{b}\theta +a$。图2描绘了 $D_1,D_2,C_1$ 在 $\theta -p_0$ 平面内的位置。

根据上面的分析我们很容易得到，(1)如果 $\left(\theta ,p_0\right)\in C_1$，系统(6)有两个相同的正平衡点 $E_2$。(2)如果 $\left(\theta ,p_0\right)\in D_1\cup D_2$，系统(6)有两个正平衡点 $E_2,E_3$。

Figure 2. Location of $D_1,D_2,C_1$ in the $\theta -p_0$ plan

图2. $D_1,D_2,C_1$ 在 $\theta -p_0$ 平面内的位置

3.3. 局部稳定性

本文利用线性化技术给出了平衡点的不同稳定性。在这种方法中，如果在平衡点处计算的雅可比矩阵的所有特征值都有负实部，则称平衡点是局部渐近稳定的(参见参考文献 [20])。在任意平衡点 $E_i$ 处，系统(6)的雅可比矩阵如下所示：

$J\left(E_i\right)=\left(\begin{array}{cc}r-\theta -\frac{2rN}{k}& \frac{r{N}^2}{k^2}\\ -\frac{\epsilon a\theta }{b}+\frac{\epsilon ak\theta }{k_{m}b}& -\frac{\epsilon \left(a-p_0\right)}{k_m}+\frac{\epsilon ak\theta }{k_{m}b}\end{array}\right)$ (13)

现在证明下列关于不同平衡点稳定性的定理，首先先来讨论零平衡点的局部稳定性，

定理3.2：系统(6)的零平衡点 $E_1$ 总是不稳定的。

证明：系统(6)在无病平衡点 $E_1$ 处的雅可比矩阵为

$J\left(E_1\right)=\left(\begin{array}{cc}r-\theta & 0\\ 0& -\frac{\epsilon \left(a-p_0\right)}{k_m}\end{array}\right)$ (14)

特征值为 $r-\theta >0$， $-\frac{\epsilon \left(a-p_0\right)}{k_m}<0$，不满足所有的特征值都为负的，故而，零平衡点 $E_1$ 总是不稳定的。

定理3.3：如果 $p_0<\psi \left(\theta \right)$，系统(6)的正平衡点 $E_2$ 是局部渐进稳定的。

证明：正平衡点 $E_2$ 的雅可比矩阵，

$J\left(E_2\right)=\left(\begin{array}{cc}-r+\theta & \frac{{\left(r-\theta \right)}^2}{r}\\ 0& -\frac{\epsilon \left(a-p_0\right)}{k_m}+\frac{\epsilon a\theta \left(r-\theta \right)}{br}\end{array}\right)$ (15)

然后， $J\left(E_2\right)$ 的特征方程为：

$\left(\lambda +r-\theta \right)\left(\lambda -\frac{\epsilon \left(a-p_0\right)}{k_m}-\frac{\epsilon a\theta \left(r-\theta \right)}{br}\right)=0$ (16)

在这种情况下，两个特征值分别为 ${\lambda }_1=-r+\theta$， ${\lambda }_2=\frac{\epsilon \left(a-p_0\right)}{k_m}-\frac{\epsilon a\theta \left(r-\theta \right)}{br}$，

当且仅当 $\theta <r,p_0<\psi \left(\theta \right)$，可得 ${\lambda }_1<0,{\lambda }_2<0$。很容易得到，只要 $p_0<\psi \left(\theta \right)$，系统(6)的正平衡点 $E_2$ 是局部渐进稳定的。

定理3.4：如果 $\psi \left(\theta \right)<p_0<a$，系统(6)的正平衡点 $E_3$ 是局部渐进稳定的。

证明：不难求得正平衡点 $E_3$ 的雅可比矩阵，

$J\left(E_3\right)=\left(\begin{array}{cc}r-\theta & \frac{{\left(r-\theta \right)}^2}{r}\\ -\frac{\epsilon d}{b}+\frac{\epsilon r\left(a-p_0\right)}{k_m\left(r-\theta \right)}& 0\end{array}\right)$ (17)

然后， $J\left(E_3\right)$ 的特征方程为：

$C_0{\lambda }^2+C_1\lambda +C_2=0$ (18)

其中，

$C_0=1,C_1=r-\theta ,C_2=\frac{\epsilon d{\left(r-\theta \right)}^2}{br}-\frac{\epsilon \left(a-p_0\right)\left(r-\theta \right)}{k_m}$

当 $\theta <r$， $\psi \left(\theta \right)<p_0<a$ 时，

$\begin{array}{l}{\Delta }_1={\lambda }_1+{\lambda }_2=-C_1<0,\\ {\Delta }_2={\lambda }_1{\lambda }_2=C_2>0,\end{array}$

可以看出，二次方程(18)的两个解将是负的或具有负实部， $E_3$ 的两个特征值都将是负的或具有负实部，因此，当 $\psi \left(\theta \right)<p_0<a$ 时，系统(6)的正平衡点 $E_3$ 是局部渐进稳定的。

在图2中，零平衡点 $E_1$ 存在时即为一个鞍点。如果 $\left(\theta ,p_0\right)\in D_1$，平衡点 $E_2$ 是局部渐进稳定的。如果 $\left(\theta ,p_0\right)\in D_2$，平衡点 $E_3$ 是局部渐进稳定的。可以发现，参数 $p_0=\psi \left(\theta \right)$，两个正平衡点交换彼此稳定性，表明跨临界分支出现的情况 [21]，此时，其中一个特征值为零，接下来将利用正规形和中心流形理论来证明的分支行为。

定理3.5：假设 $\theta <r$，当参数 $p_0$ 从 $p_0<\psi \left(\theta \right)$ 变化到 $\psi \left(\theta \right)<p_0<a$，平衡点 $E_2$ 发生了跨临界点分支，使得稳定结点 $E_2$ 变成鞍点 $E_2$，鞍点 $E_3$ 变为结点。

证明：首先将平衡点坐标转换为(0, 0)，我们通过引入新变量进行坐标变换，如下所示：

$\tilde{N}=N-\frac{\left(r-\theta \right)k_m}{r},\tilde{k}=k-k_m$

然后系统(6)变化如下：

$\left(\begin{array}{c}\frac{\text{d}\tilde{N}}{\text{d}t}\\ \frac{\text{d}\tilde{k}}{\text{d}t}\end{array}\right)=J\left(E_2\right)\left(\begin{array}{c}\tilde{N}\\ \tilde{k}\end{array}\right)+\left(\begin{array}{c}\mu \left(\tilde{N},\tilde{k}\right)\\ \nu \left(\tilde{N},\tilde{k}\right)\end{array}\right)$ (19)

其中，

$J\left(E_2\right)=\left(\begin{array}{cc}-r+\theta & \frac{{\left(r-\theta \right)}^2}{r}\\ 0& 0\end{array}\right)$

$\mu \left(\tilde{N},\tilde{k}\right)=-\frac{r{\tilde{N}}^2}{k_m}+\frac{2\left(r-\theta \right)\tilde{N}\tilde{k}}{k_m}+O\left({\left|\tilde{N},\tilde{k}\right|}^3\right)$

$\nu \left(\tilde{N},\tilde{k}\right)=-\frac{\epsilon a\theta \tilde{N}\tilde{k}}{b{k}_m}+O\left({\left|\tilde{N},\tilde{k}\right|}^3\right)$

当 $p_0=\psi \left(\theta \right)$， $J\left(E_2\right)$ 的特征值 ${\lambda }_1=\theta -r$， ${\lambda }_2=0$，

变化矩阵T

$T=\left(\begin{array}{cc}1& \frac{r-\theta }{r}\\ 0& 1\end{array}\right)$

坐标变化

$\left(\begin{array}{c}x\\ y\end{array}\right)=T^{-1}\left(\begin{array}{c}\tilde{N}\\ \tilde{k}\end{array}\right),$ (20)

令 $\left(\theta -r\right)\text{d}t=\text{d}\tau$ 将系统(19)化为标准型

$\left(\begin{array}{c}\frac{\text{d}x}{\text{d}\tau }\\ \frac{\text{d}y}{\text{d}\tau }\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)\left(\begin{array}{c}\tilde{N}\\ \tilde{k}\end{array}\right)+\left(\begin{array}{c}\tilde{\mu }\left(x,y\right)\\ \tilde{\nu }\left(x,y\right)\end{array}\right)$ (21)

其中

$\tilde{\mu }\left(x,y\right)=\frac{r{x}^2}{k_m\left(r-\theta \right)}-\frac{\left(br+\epsilon a\theta \right)\left(r-\theta \right)}{b{r}^2{k}_m}-\frac{\epsilon a\theta xy}{br{k}_m}+O\left({\left|x,y\right|}^3\right)$

$\tilde{\nu }\left(x,y\right)=\frac{\epsilon a\theta xy}{b{k}_m\left(r-\theta \right)}+\frac{\epsilon a\theta y^2}{br{k}_m}+O\left({\left|x,y\right|}^3\right)$

通过局部中心流形定理，使得存在一个中心流形 $x=h\left(y\right)=h{y}^2+o\left({\left|y\right|}^3\right)$，满足

$\frac{\text{d}x}{\text{d}\tau }=\left(2hy\right)\frac{\text{d}y}{\text{d}\tau }+O\left({\left|y\right|}^3\right)$ (22)

可得 $\left[h-\frac{\left(br+\epsilon a\theta \right)\left(r-\theta \right)}{b{r}^2{k}_m}\right]y^2+O\left({\left|y\right|}^3\right)$

所以， $h=\frac{\left(br+\epsilon a\theta \right)\left(r-\theta \right)}{b{r}^2{k}_m}$ 将 $x=\frac{\left(br+\epsilon a\theta \right)\left(r-\theta \right)}{b{r}^2{k}_m}{y}^2+o\left({\left|y\right|}^3\right)$，代入系统(21)的第一个方程，可得

$\frac{\text{d}y}{\text{d}\tau }=\frac{\epsilon a\theta y^2}{br{k}_m}{y}^2+O\left({\left|y\right|}^3\right)$ (23)

通过 $\left(\theta -r\right)\text{d}t=\text{d}\tau$ 逆缩

$\frac{\text{d}y}{\text{d}\tau }=-\frac{\epsilon a\theta \left(r-\theta \right)y^2}{br{k}_m}{y}^2+O\left({\left|y\right|}^3\right)$ (24)

由此，根据文献 [22] 定理7.1可得平衡点 $E_2$ 是一个结点，假设 $\theta <r$，若从 $p_0<\psi \left(\theta \right)$ 变化到 $\psi \left(\theta \right)<p_0<a$，稳定结点 $E_2$ 变成鞍点，鞍点 $E_3$ 变成结点，证毕。

3.4. 全局稳定性

本节我们考虑正平衡点 $E_2$， $E_3$ 的全局渐近稳定。此类研究最常用的方法是构造Lyapunov函数，但是对于系统(6)，构造Lyapunov函数相对困难，因此，我们采用Poincar-Bendixson定理(参考文献 [23])，首先证明原系统中不存在极限环。

定理3.6：系统(6)不存在极限环。

证明：在大部分文献中，证明的方法是Dulac准则。这个准则是针对于光滑向量场的，显而易见系统(6)是一个光滑的向量场。

若系统(6)存在极限环，那么它必存在于可行域内， $\Gamma =\left\{\left(N,k\right)\in R_2^{+}|0<N\le k_m,0<k\le k_m\right\}$

选择

$B\left(N,k\right)=\frac{1}{\left(k_m-k\right)N}$ (25)

设

$f_1\left(N,k\right)=rN\left(1-\frac{N}{k}\right)-\theta N\left(t\right)$ (26)

$f_2\left(N,k\right)=\epsilon \left(a\left(1-\frac{\theta N}{b}\right)-p_0\right)\left(1-\frac{k}{k_m}\right)$ (27)

然后可得，对于 $0<N\le k_m$， $0<k<k_m$

$\begin{array}{c}\frac{\partial \left(B{f}_1\right)}{\partial N}+\frac{\partial \left(B{f}_2\right)}{\partial k}=\frac{\partial }{\partial N}\left[\frac{r-\theta }{k_m-k}-\frac{rN}{\left(k_m-k\right)k}\right]+\frac{\partial }{\partial k}\left[\frac{\epsilon \left(a-p_0\right)}{k_{m}N}-\frac{\epsilon a\theta }{b{k}_m}\right]\\ =-\frac{r}{\left(k_m-k\right)k}<0\end{array}$ (28)

对于任意的 $0<k<k_m$，都有 $\frac{\partial \left(B{f}_1\right)}{\partial N}+\frac{\partial \left(B{f}_2\right)}{\partial k}$ 总是负的，证毕。

定理3.7：如果 $\left(\theta ,p_0\right)\in D_1$， $E_2$ 是全局渐进稳定的。反之，如果 $\left(\theta ,p_0\right)\in D_2$， $E_3$ 是全局渐进稳定的

证明：由定理3.2，如果 $\left(\theta ,p_0\right)\in D_1$， $E_2$ 是局部渐进稳定的，根据Poincar-Bendixson定理，在可行域 $\Gamma$ 内每个轨迹都接近一个平衡点，故而 $E_2$ 是全局渐进稳定的。

同理如果 $\left(\theta ,p_0\right)\in D_2$， $E_3$ 也是全局渐进稳定的。

4. 模拟结果

在这一部分中，我们通过将模型解与2010年至2020年报告的内蒙古年底羊存栏量数量(见表2)进行拟合，用最小二乘法拟合来确定模型的合理性，并估计参数 $p_0$ 和 $k_m$。

Table 2. Annual sheep stock in Inner Mongolia from 2010 to 2020

表2. 2010~2020内蒙古地区年度羊存栏量

4.1. 参数描述

在进行数值模拟之前，我们首先确定模型(6)中参数值。然而，涉及肉羊的数据不容易获得，我们必须依靠在线新闻、假设估计或数据拟合。在模型(6)中，时间单位为年，表3给出了参数值和变量初始值。

(A) 根据价格函数 $p\left(h\right)=a\left(1-h/b\right)$，根据羊的市场价格和市场供应量，采用最小二乘法进行线性拟合，可以得到参数a和b的值，可以见表3。

(B) 在新疆地区，羊群的净生殖率约为153%，世代时间约为15~18个月，由此可得羊群的內禀增长率为 $r=\ln \left(1.53\right)/0.67\approx 0.83$ [24]。羊的存栏时间大约为两到三年，因此估计出羊的出栏率为 $\theta =0.48$。

(C) 当羊肉市场价格超过价格阈值时，养殖户有两种选择，一种是扩张养殖规模，一种是保持现有的养殖规模，故而比例系数可以确定为0.5。在该模型中，时间单位为年。参数值见表3。

Table 3. Values of parameters in model (6)

表3. 模型(6)中参数的取值

4.2. 数值模拟

利用系统(6)，模拟了内蒙古地区2010年到2020年年底的羊存栏量的数据，并预测未来肉羊存栏量的发展趋势。 $\widetilde{N\left(t\right)}$ 是实际上内蒙古地区每年年底的存栏量， $N\left(t\right)$ 是理论上的第t年每年年底的存栏量，满足，

$\frac{\text{d}N}{\text{d}t}=rN\left(1-\frac{N}{k}\right)-\theta N$

根据表2中2010年到2020年肉羊的存栏数据，最小二乘拟合用于获取参数值 $\theta$ 和 $k_0$，满足：

$\min {\displaystyle \underset{t=2010}{\overset{2020}{\sum }}{\left[\widetilde{N\left(t\right)}-N\left(t\right)\right]}^2}$

由于信息错误和实践中的漏报，肉羊存栏的数据存在一定的随机误差。因此，假设实际数据 $\widetilde{N\left(t\right)}$ 存在泊松分布误差。进行了1000次随机模拟，数据拟合结果如图3所示，数值模拟表明，具有合理参数值的模型与报告的实际数据具有良好的匹配性，估计的参数值如表2所示。

Figure 3. The fitting results between the dynamic model and the stock of sheep, in which the blue dotted line is the 95% confidence interval, the blue solid line is the average value, and the red star point is the actual data

图3. 动力学模型与羊存栏量的拟合结果，其中蓝色虚线为95%置信区间，蓝色实线为平均值，红色星点为实际数据

4.3. 单参数变化下的分支描述

本节我们探究模型中重要参数 $p_0,\theta ,b$ 对当前养殖量长久发展的潜在影响。根据数值模拟我们预测了参数值，接下来考虑以下两组假设参数值：

(a) $r=0.83,\theta =0.4,\epsilon =0.5,a=4506,b=4.692\times {10}^7,k_m=7.338\times {10}^7,p_0=2000$

(b) $r=0.83,\theta =0.4,\epsilon =0.5,a=4506,b=4.692\times {10}^7,k_m=7.338\times {10}^7,p_0=4000$

通过第3节的分析，确保了这两组参数下，平衡点都为正的，代入具体数据得到 $\psi \left(\theta \right)=3046$。

图4(a)为(a)组参数下状态变量的稳定状态，羊的存栏量约为3.8 × 10⁷，最大养殖量约为7.338 × 10⁷，注意这里养殖户的最大存栏量为可以扩张的最大存栏量，在图4(b) (b)组参数下羊的最大存栏量比(a)组减少，相应的羊的存栏量也对应缩小为1.346 × 10⁷，通过(a)组和(b)组参数差别，主要是价格阈值带来的系统平衡点稳定性的变化。接下来通过分支图来直观的感受参数变化下，系统(6)平衡点的稳态变化。

(a) (b)

Figure 4. The time evolution of model (6) when the initial condition is 3.6 × 10⁷, 5.4 × 10⁷. (a) Figure E₂ is local asymptotic stability under (a) group parameters, (b) figure E₃ local asymptotic stability under (b) group parameters

图4. 初始条件为(3.6 × 10⁷, 5.4 × 10⁷)时，模型(6)的时间演化。(a) 图(a)组参数下E₂局部渐进稳定，(b) 图(b)组参数下E₃局部渐进稳定

图5展示在确定其他的参数下，关于参数 $p_0$ 的单参数分支图，其中蓝色实线代表的是稳定平衡点的变化趋势，红点表示分支点。可见，正平衡点随着 $p_0$ 的变化发生了两个分支点，当 $p_0=3046$ 时，发生了跨临界分岔，两个平衡点相交并交换其稳定性。在 $p_0<3046$ 时，平衡点 $E_2$ 稳定，在 $p_0>3046$ 时，转化为平衡点 $E_3$ 稳定。

需要指出的是，该系统 $p_0$ 是养殖户选择扩张的价格阈值，表示环境承载量增加的一个阈值，这在实际情况中，养殖的成本是价格阈值的一个重要影响因素，价格阈值大代表养殖户盈利需要的羊的市场价格高，即相应的成本高。因此在 $p_0>3046$ 时，养殖场最大的养殖量，即养殖规模减少。反之，当 $p_0<3046$ 时，养殖户盈利需要的羊市场价格底，养殖户会选择用最大的养殖规模进行养殖，即选择用最大的养殖量 $k_m$ 如图5(b)所示，相应的也可以发现在 $p_0=3046$ 左右当前养殖量的变化(图5(a))，再次观察图5也可以发现，当 $p_0<3046$ 存栏量和最大养殖量受价格阈值的影响不大，即价格阈值不影响养殖发展的最终状态，此时的养殖成本是养殖户可以接受的范围，不会影响养殖户的长期决策，但是当 $p_0>3046$ 价格阈值越高，存栏量随着最大养殖量的减少而减少，成本高带来的养羊的收益小，促使存栏量的减少。

(a) (b)

Figure 5. Single parameter branching diagram of parameters, p₀ in which the blue solid line represents the change trend of the stable equilibrium point and the red dot represents the branching point. Other parameters are set to r = 0.83, θ = 0.4, a = 4506, b = 4.692 × 10⁷, k_m = 7.338 × 10⁷

图5. 参数p₀的单参数分支图，其中蓝色实线代表的是稳定平衡点的变化趋势，红点表示分支点。其中其他参数设定为r = 0.83，θ = 0.4，a = 4506，b = 4.692 × 10⁷，k_m = 7.338 × 10⁷

对于其他重要参数 $\theta ,b$ 也绘制了相同的分支图，对于参数 $\theta$ 的分支图(图6)可以清楚的看出这里存在两个跨临界分支点，在 $\theta =0.08$ 、 $\theta =0.72$ 发生了跨临界分岔，两个平衡点相交并交换其稳定性 $E_2\to E_3\to E_2$。由此可以看出随着参数 $\theta$ 的增长，当前存栏量 $N\left(t\right)$ 逐渐降低，直至为零，但是在0.08~0.72时降低速率降低，接近于双曲线趋势降低。由此可见对于养殖场来说，想要长久养殖，其出栏量未必是越大越好，出栏量越大，意味着养殖场消失得越快，这与实际情况类似，这也就反映了，养殖场养殖过程中将需要采用合理的出栏率，如何决定这个出栏率，是我们未来需要考虑的一个方向。

再次通过观察图6中(c)和(d)可以发现，当 $b=4.691\times {10}^7$ 时，发生了跨临界点分支，随着参数b的增加，养殖量随之增加，当增加到一定程度之后，养殖量保持不变，此时的最大养殖量是该地所能扩张的最大养殖量。随着社会的发展市场需求也不是稳定不变的，是波动变化的，故而可以根据市场需求量确定当前的养殖量，以避免造成社会资源的浪费。

(a) (b)(c) (d)

Figure 6. (a) (b) for the single parameter branching graph with parameters θ, two transcritical branches occur. (c) (d) On the single parameter bifurcation graph of parameters b, when b = 4.691 × 10⁷ a transcritical bifurcation occurs.

图6. (a) (b)关于参数θ的单参数分支图，发生了两个跨临界分支。(c) (d)关于参数b的单参数分支图当b = 4.691 × 10⁷发生了一个跨临界分支

5. 结论

本文以羊养殖厂为研究对象，考虑市场价格影响下的羊养殖模型，建立了受价格动态影响的养殖模型，主要考虑了价格变化对最大养殖量的影响，通过调节最大养殖量的变化，进而来影响养殖量的变化。当羊市场价格超过价格阈值时，养殖户会选择扩大养殖规模，当羊市场价格低于价格阈值时，养殖户会选择缩减养殖规模，进而形成了一个价格对养殖过程的影响状态，如图1所示，根据分析建立了系统(5)价格影响下的动态养殖模型。

该模型有三个可行平衡点，它们是与时间无关的状态变量值。这三个平衡点，其中零平衡点对应与不存在该物种，显然这是不现实的，它也只会在出栏率超过内柄增长率的时候才会出现。其他两个正平衡点是在不同的参数条件下是稳定的，在临界点处发生了跨临界点分支。本文选择了三个重要的参数，观察他们的变化是如何导致稳定性的切换。通过单参数分支图可以发现，在 $p_0<\psi \left(\theta \right)$ 时，平衡点 $E_2$ 稳定，在 $p_0>\psi \left(\theta \right)$ 时，转化为平衡点 $E_3$ 稳定。一旦价格阈值超过某一个界限，平衡点的稳定状态立马发生了改变。在羊价格阈值较低时，系统的稳定状态不变的，可以发现此时，价格对养殖户长期的决策没有影响。我们知道养殖成本是价格阈值的一个重要决定因素，因此价格阈值的高低可以评价成本的高低，根据单参数分支图(图5)，当养殖成本低时，羊的价格是不会影响养殖户的长期决策。需要注意的是，一旦价格阈值超过某一个限定值，养殖量和最大养殖量会随着价格阈值的增加而减少，意味着如果养殖成本高，那么养殖将会逐渐缩小，以致该行业的消失。实际中影响养殖量成本由多种条件决定，受其他行业的或直接或间接的影响，由此养殖户可以通过当前的成本价格，决定当前的养殖量和养殖规模，以来降低自己的养殖成本。

描述市场需求的参数b与养殖量关系的单参数分支图也有着有趣的结果，如果b增加，市场需求增加，养殖量相应的也会随着市场需求的增加而增加，如果增加到一定程度，即增加到地区所能承受的最大范围后不再增加，这种情况下，一直保持着地区的最大养殖。故而养殖户可以根据市场需求确定自己的养殖量，避免盲目的出栏扩张，导致资源的浪费。由此可见对于养殖场来说，想要长久养殖，其出栏量未必是越大越好，出栏量越大，意味着养殖场消失得越快，这与实际情况类似，这也就反映了，养殖场养殖过程中将需要采用合理的出栏率，如何决定这个出栏率，是我们未来需要考虑的一个方向。如前所述，该模型可以反映市场价格对养殖过程的动态影响，通过单参数分支分析可以发现随着成本、需求、出栏的变化，养殖状态的最终状态也会随之变化，养殖户可以根据变化调整自己的养殖状态，以达到长久发展的目的。

本文相较于之前模型的区别在于，它将价格的影响考虑到养殖模型中去，考虑价格影响下的动态养殖方程，虽然这个工作反映了现实的一个生态经济问题，可以更加准确的反映养殖状态，为养殖户的养殖提供一定的建议，但是仍有一些问题需要解决，例如，不同状态下的收入水平，养殖户养殖的主要目的还是收益，养殖户该如何控制参数以达到自己收入最高。

总之，考虑价格动态影响下的种群模型，讨论重要参数对系统稳定性的影响，对养殖式单种群模型的研究有着重要的意义，也为养殖户提供了一定的建议，本文讨论的模型也同样适用与其他单种群开发模型，对其研究生物经济提供了一定的帮助。

基金项目

国家自然科学基金(12101443, 11801398)；山西省自然科学基金(20210302124260, 202103021224095)。

NOTES

^*通讯作者。