太原理工大学数学学院，山西 太原

收稿日期：2023年4月28日；录用日期：2023年5月21日；发布日期：2023年5月31日

摘要

随着社会经济发展，羊肉制品的需求量与日俱增，羊养殖规模不断扩大，预防羊疫病工作的难度系数也日渐严峻。目前我国羊养殖业依然受到布鲁氏菌病的威胁，由于羊群养殖量和羊肉价格有直接关系，羊肉价格和羊肉的供需有直接关系，本文在考虑布鲁氏菌病对羊群养殖的影响时加入了羊肉价格和羊肉市场库存量这两个变量进而建立数学模型，给出模型的基本再生数，无病平衡点和正平衡点并分别讨论了稳定性。随后在传染病模型基础上加入三个控制措施，分别为减少易感羊和染病羊的接触；增大对感染羊的捕杀力度；减少布鲁氏菌对环境的排放；以此建立控制模型。以养殖户损失最小为目的建立目标函数，寻求对养殖户最有利的控制策略，利用Pontryagin极小值原理得出最优控制策略的理论结果。以河南羊的活重价格和人间布鲁氏菌病的累计病例数为例，对模型进行数值模拟。最后对上述三个控制措施进行组合，给出七种不同的策略。经研究表明，相比起单一控制策略和双重控制策略，三项措施同时实施对养殖户来说损失最小，控制成本、羊患病数和环境中布病浓度达到最低，故本文给养殖户的建议是综合实施上述三项措施最有利。

关键词

布鲁氏菌病，价格，最优控制，Pontryagin极小值原理

Modeling of Brucellosis Dynamics and Its Control Measures Based on Price Factors

Yuxin Wang, Yuzhen Chai^*, Mingtao Li, Junjun Liu

School of Mathematics, Taiyuan University of Technology, Taiyuan Shanxi

Received: Apr. 28^th, 2023; accepted: May 21^st, 2023; published: May 31^st, 2023

ABSTRACT

With the socio-economic development, the demand for sheep meat products is increasing day by day, the scale of sheep farming is expanding, and the difficulty factor of sheep disease prevention work is getting more and more severe. At present, the sheep farming industry in China is still threatened by brucellosis. Since the amount of sheep farming and the price of sheep meat are directly related, and the price of sheep meat and the supply and demand of sheep meat are directly related, this paper adds two variables, the price of sheep meat and the amount of sheep meat market inventory, when considering the influence of brucellosis on sheep farming and then establishes a mathematical model, gives the basic regeneration number of the model, the disease-free equilibrium point and the positive equilibrium point and discusses respectively. The stability of the model was discussed. Three control measures were then added to the infectious disease model, namely, reducing contact between susceptible and infected sheep; increasing the culling of infected sheep; and reducing brucellosis emissions to the environment; thus establishing the control model. The objective function was established with the aim of minimizing farmers’ losses, and the most favorable control strategy for farmers was sought. The theoretical results of the optimal control strategy were derived using the Pontryagin principle of minimal values. The model is numerically simulated with the live weight price of Henan sheep and the cumulative number of human brucellosis cases. Finally, the three controls mentioned above were combined to give seven different strategies. It was shown that compared to the single control strategy and dual control strategy, the simultaneous implementation of the three measures was the least lossy for farmers, with the lowest control cost, number of sheep diseases and brucellosis concentration in the environment, so the recommendation to farmers in this paper is that the combined implementation of the above three measures is the most beneficial.

Keywords:Brucellosis, Price, Optimal Control, Pontryagin Minimal Value Principle

Copyright © 2023 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

布鲁氏菌病是由布鲁氏杆菌引起的一种人畜共患的传染病，它能导致家畜的流产、不育，人感染后会出现发热、全身无力等症状，严重危害到畜牧业的生产和人类的健康。因此，要加强对这一问题的研究 [1] 。世界动物卫生组织将布病列为必须报告的动物疫病，我国将其列为二类动物疫病 [2] 。必须要加强对牛羊布鲁氏菌病的研究，减少养殖户的经济损失，促进畜牧业长远发展 [3] 。最优控制理论作为现代控制理论的重点内容，被广泛地应用到传染病模型中 [4] [5] 。最优控制为研究预防疾病的传播提供了一种新的思路为解决相关问题提供了良好策略，对整个地区控制疾病具有积极意义。结合羊布鲁氏菌病的发病机理建立一类相关的控制模型寻求最优控制策略从而控制染病类数量，对减小疾病传播有现实意义。

Aȉnseba等人 [6] 提出了绵羊布鲁氏菌病的两种传播途径：直接接触感染个体和间接暴露于污染环境。Li等人 [7] 提出了一个确定性模型来描述羊布鲁氏菌病在内蒙古兴安盟的传播。Hou等人 [8] 通过分析内蒙古人羊群布病发病情况，建立动力学模型，结果发现结合环境消毒和易感羊群免疫措施能有效控制内蒙古布鲁氏菌病。Zhou等人 [9] 建立人羊耦合的动力学模型，提出了一个多目标优化问题，利用加权和方法将其转化为一个关于最小化控制总成本的标量优化问题，利用庞特里亚金的极大值原理，得到最优控制策略。Nannyonga等人 [10] 利用数学模型研究了羊布鲁氏菌病直接从受感染个体、通过接触受污染环境或通过母婴垂直传播时的动态变化。Nyerere等人 [11] 对控制措施进行了成本效益分析，结果表明，牲畜疫苗接种、通过屠宰逐步扑杀，环境卫生和人类个人防护的结合具有较高的影响和较低的预防成本。Nie等 [12] 认为输入、扑杀和消毒是影响吉林省布鲁氏菌病的重要因素。后来，他们发现布鲁氏菌病的爆发具有明显的季节性和周期性繁殖机制 [13] [14] 。Beauvais等人使用了一个根据季节的过渡模型来研究牛羊混合养殖场的传播。考虑到年龄结构 [15] [16] [17] [18] ，Li和Hou等提出了关于青年羊和成羊的多阶段动态模型。

根据经济理论，当需求大于供给时，市场价格就会上涨，因此根据需求和供给对价格的影响， [19] [20] [21] 结合线性需求函数，非线性需求函数和实际供给函数，研究了各种形式的价格差异模型并进行了一系列动力学研究和实际应用。结合农户行为考虑价格与养殖数量的交互作用，Wang等人 [22] 在建立模型时考虑了动态价格因素。上述学者考虑了各种各样的传染病模型，没有在模型中加入价格影响，以及随着价格这一变量的出现，还会出现供需，库存等的影响，这些因素对布病的控制都有着不可忽视的影响。实施控制措施时需要一定的成本，所以使得养殖户在整个过程中所造成的损失最小也是一个关键性的问题。

本文以羊养殖为例，根据市场价格变化的动态影响建立相应动力学模型，研究其内在机理，随后加入控制策略寻求对养殖户最有利的措施。本文的研究结构如下：在第2节中给出不加疾病时的基础模型和传染病模型，并给出了各个参数的含义。第3节利用定性理论证明了无病平衡点的全局稳定性和正平衡点的局部稳定性。第4节给出最优控制模型以及给出三个控制项，建立目标函数，然后利用Pontryagin最大化原理得出最优控制策略的理论结果。第5节以河南省羊的活重价格和人间布鲁氏菌病的累计病例数为例，进行参数估计。第6节给出最优控制策略的分析，对给出三个控制措施进行组合，一共可以给出七种不同的策略。最后在第7节中阐述了本文的主要结论。

2. 模型建立

2.1. 基础模型

首先考虑无传染病时羊群养殖的模型，此模型共包括羊群养殖量N，羊肉市场价格p和羊肉的市场库存量M三个变量。

养殖户在养殖羊群时，羊群的入栏和出栏都是会受到羊肉市场价格的影响。假设羊群的入栏分为两个部分，一部分每周期内羊群规模以常数A进行补充，另一部分受羊肉市场价格的影响，当出栏价格大于养殖成本价格 $p_0$ 时，养殖户为了获取利润会扩大养殖规模，增大羊群输入，相反当出栏价格小于初始养殖价格时，会缩小养殖规模，减少羊群输入，因此种群规模受价格影响的这部分设为 $\epsilon e\left(p-p_0\right)$ ，ε为价格调节系数，e为羊出栏时的平均重量。羊群的出栏与价格以及养殖量都有正比的关系，α为出栏调节系数，αpN为羊群出栏量。

商品价格会受到供需变化的调节：需求越大价格越高，销售供应越多价格越低。假设存在最大需求量，由于羊肉为非必需品，则需求函数关于价格为反比例关系，当价格越高，需求量越小，需求调节率为b，需求函数为 $D-bep$ ，且需求不可能为负值，因此 $D-bep>0$ 。μ是销售供应系数，M为库存量。销售供应与库存量成正比，价格上涨库存量也会增加，但是随着价格的升高，销售供应不会无限的增大，

故销售供应函数定义为 $\frac{\mu Mp}{1+p}$ ， $\varphi$ 为供需调节价格系数。

羊群出栏量αpN定义为生产供应量， $\frac{\mu Mp}{1+p}$ 为销售供应量，这两个供应量的差值即为市场的库存量。

综上分析的结果，得到羊群养殖量，价格，库存的模型：

$\{\begin{array}{l}\frac{\text{d}N}{\text{d}t}=A+\epsilon e\left(p-p_0\right)-\alpha pN\hfill \\ \frac{\text{d}p}{\text{d}t}=\varphi \left(D-bep-\frac{\mu Mp}{1+p}\right)\hfill \\ \frac{\text{d}M}{\text{d}t}=\alpha pN-\frac{\mu Mp}{1+p}\hfill \end{array}$ (1)

根据定义可知需求函数 $D-bep>0$ ，从而 $Ab+D\epsilon -be{p}_0\epsilon >0$

此外中国羊肉供不应求的背景下 $D-bep>A$ ，进而可以得到 $D>A$

因此模型有唯一的正平衡点 $E^{\ast }=\left(N^{\ast },p^{\ast },M^{\ast }\right)$

$N^{\ast }=\frac{e\left(Ab+D\epsilon -b\epsilon e{p}_0\right)}{\alpha \left(D+\epsilon e{p}_0-A\right)}$ ， $p^{\ast }=\frac{D+\epsilon e{p}_0-A}{e\left(\epsilon +b\right)}$

$M^{\ast }=\frac{\left(Ab+D\epsilon -b\epsilon e{p}_0\right)\left(\epsilon e+be+D+\epsilon e{p}_0-A\right)}{\mu \left(D+\epsilon e{p}_0-A\right)\left(\epsilon +b\right)}$

和正不变集

$\Gamma =\left\{\left(N,P,M\right)|N,P,M\ge 0,0\le N\le \frac{Ab+D\epsilon -b\epsilon e{p}_0}{b},0\le p\le \frac{D}{be},0\le M\le \frac{\alpha D\left(Ab+D\epsilon -b\epsilon e{p}_0\right)}{e{b}^2}\right\}$

定理2.1：系统(1)的平衡点是局部渐近稳定的。

证明：该平衡点处的雅可比矩阵为

${J|}_{E^{\ast }}=\left(\begin{array}{ccc}-\alpha p^{\ast }& \epsilon e-\alpha N^{\ast }& 0\\ 0& -\varphi \left[be+\frac{\mu M^{\ast }}{{\left(1+p^{\ast }\right)}^2}\right]& -\frac{\varphi \mu p^{\ast }}{1+p^{\ast }}\\ \alpha p^{\ast }& \alpha N^{\ast }-\frac{\mu M^{\ast }}{{\left(1+p^{\ast }\right)}^2}& -\frac{\mu p^{\ast }}{1+p^{\ast }}\end{array}\right)$

其中：

$D_1=-\alpha p^{\ast }<0$

$D_2=\left|\begin{array}{cc}-\alpha p^{\ast }& \epsilon e-\alpha N^{\ast }\\ 0& -\varphi \left[be+\frac{\mu M^{\ast }}{{\left(1+p^{\ast }\right)}^2}\right]\end{array}\right|=\alpha p^{\ast }\varphi \left[be+\frac{\mu M^{\ast }}{{\left(1+p^{\ast }\right)}^2}\right]>0$

$D_3=\left|\begin{array}{ccc}-\alpha p^{\ast }& \epsilon e-\alpha N^{\ast }& 0\\ 0& -\varphi \left[be+\frac{\mu M^{\ast }}{{\left(1+p^{\ast }\right)}^2}\right]& -\frac{\varphi \mu p^{\ast }}{1+p^{\ast }}\\ \alpha p^{\ast }& \alpha N^{\ast }-\frac{\mu M^{\ast }}{{\left(1+p^{\ast }\right)}^2}& -\frac{\mu p^{\ast }}{1+p^{\ast }}\end{array}\right|=-\alpha p^{\ast }\frac{\varphi \mu p^{\ast }}{1+p^{\ast }}\left[\left(\epsilon +b\right)e+\frac{2\mu M^{\ast }}{{\left(1+p^{\ast }\right)}^2}\right]<0$

可得出矩阵 ${J|}_{E^{\ast }}$ 是负定矩阵，其特征根都小于0，则平衡点是局部渐近稳定的。

2.2. 传染病模型

基于无传染病时的养殖模型(1)，结合人间和羊群布鲁氏菌病的传播特征，建立了羊，人，价格和库存的动态模型。羊群包括易感羊群S和已感染羊群I，人群包括易感人群 $S_h$ 和染病人群 $I_h$ 。W表示环境中布鲁氏菌浓度，p表示市场上卖出羊的价格。

此处说明一下和模型(1)的不同之处。羊中布鲁氏菌病的感染途径包括与受感染羊的接触和与环境中的病菌接触，因此传播以 ${\beta }_{1}SI$ 和 ${\beta }_{1}SW$ 的双线性作用的形式发生，其中 ${\beta }_1$ 和 ${\beta }_2$ 易感分别为两种感染方式的感染率。当易感羊被感染时，就会进入被感染的仓室I，羊群会定期进行抽检，一旦抽检发现被感染的羊就会将其捕杀且在这个抽检过程中的最大感染羊的捕杀率设置为γ。环境中布鲁氏菌的来源是受感染的羊，受感染羊单位时间排放到环境中的布鲁氏菌量记为kI，环境中布鲁氏菌的浓度单位时间衰减率记为δW。布鲁氏菌病具有家畜传播人而人不传播家畜，人间不传播的特点。故易感人群主要通过接触受感染的羊和环境中布鲁氏菌进行感染，感染发病率为 ${\beta }_h{S}_{h}I+{\beta }_W{S}_{h}W$ ，以常数 $A_h$ 为人口出生率， $d_h$ 为人口死亡率，m为疾病康复率。

因此，该数学模型被描述为以下常微分方程：

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=A+\epsilon e\left(p-p_0\right)-{\beta }_{1}SI-{\beta }_{2}SW-\alpha pS\hfill \\ \frac{\text{d}I}{\text{d}t}={\beta }_{1}SI+{\beta }_{2}SW-\gamma I\hfill \\ \frac{\text{d}W}{\text{d}t}=kI-\delta W\hfill \\ \frac{\text{d}p}{\text{d}t}=\varphi \left(D-bep-\frac{\mu Mp}{1+p}\right)\hfill \\ \frac{\text{d}M}{\text{d}t}=\alpha pS-\frac{\mu Mp}{1+p}\hfill \\ \frac{\text{d}{S}_h}{\text{d}t}=A_h-{\beta }_h{S}_{h}I-{\beta }_w{S}_{h}W-d_h{S}_h+m{I}_h\hfill \\ \frac{\text{d}{I}_h}{\text{d}t}={\beta }_h{S}_{h}I+{\beta }_w{S}_{h}W-d_h{I}_h-m{I}_h\hfill \end{array}$ (2)

整个羊群为 $N=S+I$ ，将系统(2)的前两个方程相加得到

$\frac{\text{d}N}{\text{d}t}=A+\epsilon e\left(p-p_0\right)-\alpha pS-\gamma I$

整个人口数为 $N_h=S_h+I_h$ ，将系统(2)的后两个方程相加得到

$\frac{\text{d}{N}_h}{\text{d}t}=A_h-d_h{N}_h$

系统(2)的正不变集为

$\begin{array}{l}X=\{\left(S,I,W,P,M\right)|S,I,W,P,M\ge 0,0\le S+I\le \frac{Ab+D\epsilon -b\epsilon e{p}_0}{b},0\le W\frac{k\left(Ab+D\epsilon -b\epsilon e{p}_0\right)}{b\delta },\\ 0\le p\le \frac{D}{be},0\le M\le \frac{\alpha D\left(Ab+D\epsilon -b\epsilon e{p}_0\right)}{e{b}^2},0\le S_h+I_h\le \frac{A_h}{d_h}\}\end{array}$

3. 模型的动力学分析

由于系统(2)中最后两个方程独立于前五个方程，因此对(2)的动力学分析等同于对下列系统分析：

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=A+\epsilon e\left(p-p_0\right)-{\beta }_{1}SI-{\beta }_{2}SW-\alpha pS\hfill \\ \frac{\text{d}I}{\text{d}t}={\beta }_{1}SI+{\beta }_{2}SW-\gamma I\hfill \\ \frac{\text{d}W}{\text{d}t}=kI-\delta W\hfill \\ \frac{\text{d}p}{\text{d}t}=\varphi \left(D-bep-\frac{\mu Mp}{1+p}\right)\hfill \\ \frac{\text{d}M}{\text{d}t}=\alpha pS-\frac{\mu Mp}{1+p}\hfill \end{array}$ (3)

3.1. 无病平衡点

设系统(3)中所有方程的右端为0，可以得到一个无病平衡点 $E_0=\left(S^0,0,0,p^0,M^0\right)$ ，其中

$S^0=\frac{e\left(Ab+D\epsilon -b\epsilon e{p}_0\right)}{\alpha \left(D+\epsilon e{p}_0-A\right)}$ ， $p^0=\frac{D+\epsilon e{p}_0-A}{e\left(\epsilon +b\right)}$

$M^0=\frac{\left(Ab+D\epsilon -b\epsilon e{p}_0\right)\left(\epsilon e+be+D+\epsilon e{p}_0-A\right)}{\mu \left(D+\epsilon e{p}_0-A\right)\left(\epsilon +b\right)}$

3.2. 基本再生数

根据 [23] 和 [24] 基本再生数以及下一代矩阵的定义，考虑以下包含致病项的辅助系统：

$\{\begin{array}{l}\frac{\text{d}I}{\text{d}t}={\beta }_{1}SI+{\beta }_{2}SW-\gamma I\hfill \\ \frac{\text{d}W}{\text{d}t}=kI-\delta W\hfill \end{array}$

从文章 [23] 中的方法，得到

$F=\left(\begin{array}{cc}{\beta }_1{S}^0& {\beta }_2{S}^0\\ 0& 0\end{array}\right),V^{-1}=\left(\begin{array}{cc}\frac{1}{\gamma }& 0\\ \frac{k}{\delta \gamma }& \frac{1}{\delta }\end{array}\right)$

则系统(3)的基本再生数为

$R_0=\rho \left(F{V}^{-1}\right)=\frac{e\left(Ab+D\epsilon -b\epsilon e{p}_0\right)\left({\beta }_1\delta +{\beta }_{2}k\right)}{\alpha \delta \gamma \left(D+\epsilon e{p}_0-A\right)}$

3.3. 正平衡点

设系统(3)中所有方程的右端为0，当 $R_0>1$ 时，可以得到系统(3)的正平衡点 $E^1=\left(S^1,I^1,W^1,p^1,M^1\right)$

$S^1=\frac{\delta \gamma }{{\beta }_1\delta +{\beta }_{2}k}>0$

$I^1=\frac{e\left(Ab+D\epsilon -be\epsilon p_0\right)\left({\beta }_1\delta +{\beta }_{2}k\right)-\alpha \delta \gamma \left(D+\epsilon e{p}_0-A\right)}{\alpha \delta {\gamma }^2+be\gamma \left({\beta }_1\delta +{\beta }_{2}k\right)}>0$

$W^1=\frac{ke\left(Ab+D\epsilon -be\epsilon p_0\right)\left({\beta }_1\delta +{\beta }_{2}k\right)-k\alpha \delta \gamma \left(D+\epsilon e{p}_0-A\right)}{\alpha {\delta }^2{\gamma }^2+be\delta \gamma \left({\beta }_1\delta e+{\beta }_{2}k\right)}>0$

$p^1=\frac{D\left({\beta }_1\delta +{\beta }_{2}k\right)}{\alpha \delta \gamma +be\left({\beta }_1\delta +{\beta }_{2}k\right)}>0$

$M^1=\frac{beu{\left({\beta }_1\delta +{\beta }_{2}k\right)}^2+\alpha \gamma u\delta \left({\beta }_1\delta +{\beta }_{2}k\right)}{\alpha \delta \gamma \left(\alpha \delta \gamma +\left(D+be\right)\left({\beta }_1\delta +{\beta }_{2}k\right)\right)}>0$

3.4. 无病平衡点的全局稳定

定理3.1：当 $R_0<1$ ，系统的无病平衡点 $E_0=\left(S^0,0,0,p^0,M^0\right)$ 是全局稳定的。

证明：首先证明无病平衡点的局部稳定性，根据 [23] 中无病平衡点是局部稳定的所需满足的五个条件中显然前四个条件是满足的，下证无病平衡点的雅可比矩阵的所有特征值为具有负实部。

根据 [25] 中定义 $s\left(M\right)$ 为矩阵 $M$ 的最大特征值，其中 $M=F-V$ 。 [23] 中的定理2可得： $R_0<1\iff s\left(M\right)<0$ ， $R_0>1\iff s\left(M\right)>0$

对系统(3)的各项以 $I,W,p,M,S$ 重排，得到无病平衡点的雅可比矩阵：

${J|}_{E_0}=\left(\begin{array}{cc}M& 0\\ J_3& J_4\end{array}\right)$

其中

$J_3=\left(\begin{array}{cc}0& 0\\ 0& 0\\ -{\beta }_1{S}^0& -{\beta }_2{S}^0\end{array}\right),J_4=\left(\begin{array}{ccc}-\varphi be-\frac{\varphi \mu M^0}{{\left(1+p^0\right)}^2}& -\frac{\varphi \mu p^0}{1+p^0}& 0\\ \alpha S^0-\frac{\mu M^0}{{\left(1+p^0\right)}^2}& -\frac{\mu p^0}{1+p^0}& \alpha p^0\\ \epsilon e-\alpha S^0& 0& -\alpha p^0\end{array}\right)$

当 $R_0<1$ 时，有 $s\left(M\right)<0$ ，因此只需证明 $J_4$ 的特征根都为负，其中：

$D_1=-\varphi be-\frac{\varphi \mu M^0}{{\left(1+p^0\right)}^2}<0$

$D_2=\left|\begin{array}{cc}-\varphi be-\frac{\varphi \mu M^0}{{\left(1+p^0\right)}^2}& -\frac{\varphi \mu p^0}{1+p^0}\\ \alpha S^0-\frac{\mu M^0}{{\left(1+p^0\right)}^2}& -\frac{\mu p^0}{1+p^0}\end{array}\right|=\frac{\varphi \mu p^0}{1+p^0}\left(be+\alpha S^0\right)>0$

$D_3=\left|\begin{array}{ccc}-\varphi be-\frac{\varphi \mu M^0}{{\left(1+p^0\right)}^2}& -\frac{\varphi \mu p^0}{1+p^0}& 0\\ \alpha S^0-\frac{\mu M^0}{{\left(1+p^0\right)}^2}& -\frac{\mu p^0}{1+p^0}& \alpha p^0\\ \epsilon e-\alpha S^0& 0& -\alpha p^0\end{array}\right|=-\alpha p^0\frac{\varphi \mu p^0}{1+p^0}\left(be+\epsilon e\right)<0$

可得出矩阵 $J_4$ 是负定矩阵，其特征根都小于0，则平衡点是局部稳定的。下证明其全局稳定：

从系统(3)的第一个方程可知

$\frac{\text{d}S}{\text{d}t}=A+\epsilon e\left(p-p_0\right)-{\beta }_{1}SI-{\beta }_{2}SW-\alpha pS\le A+\epsilon e\left(\frac{D}{be}-p_0\right)-\alpha \frac{D}{be}S$

因无传染病时的模型(1)的唯一平衡点 $\left(N^{\ast },p^{\ast },M^{\ast }\right)$ 是局部稳定的，则

$\underset{t\to \infty }{\lim }N=\frac{Abe+De\epsilon -b\epsilon e^2{p}_0}{\alpha \left(D+\epsilon e{p}_0-A\right)},\underset{t\to \infty }{\lim }p=\frac{D+\epsilon e{p}_0-A}{\epsilon e+be}$

存在 $t_i>0$ 使得对任意的 $t>t_i$ 时，有下列式子成立：

$\underset{t\to \infty }{\lim }\sup S=\frac{Abe+De\epsilon -b\epsilon e^2{p}_0}{\alpha D}=S^0+{\epsilon }_1$

其中：

${\epsilon }_1=\frac{\left(Abe+De\epsilon -b\epsilon e^2{p}_0\right)\left(\epsilon e{p}_0-A\right)}{\alpha D\left(D+\epsilon e{p}_0-A\right)}$

当 $t>t_i$ 时，考虑以下辅助系统

$\{\begin{array}{l}\frac{\text{d}I}{\text{d}t}\le \left(S^0+{\epsilon }_1\right)\left({\beta }_{1}I+{\beta }_{2}W\right)-\gamma I\hfill \\ \frac{\text{d}W}{\text{d}t}=kI-\delta W\hfill \end{array}$

相对应的线性系统为：

$\{\begin{array}{l}\frac{\text{d}{I}^{\prime }}{\text{d}t}=\left(S^0+{\epsilon }_1\right)\left({\beta }_1{I}^{\prime }+{\beta }_2{W}^{\prime }\right)-\gamma I^{\prime }\hfill \\ \frac{\text{d}{W}^{\prime }}{\text{d}t}=k{I}^{\prime }-\delta W^{\prime }\hfill \end{array}$ (4)

令 $M_1=M+{\epsilon }_2{M}_0$ ，其中 $M_0$ 的定义如下：

$M_0=\left(\begin{array}{cc}{\beta }_1& {\beta }_2\\ 0& 0\end{array}\right)$

根据文章 [23] 中的定理2可知，存在充分小的正数 ${\epsilon }_2$ ，使得当 $R_0<1$ 时有 $s\left(M_1\right)<0$ 。由Perron-Frobenius定理可知，当 $s\left(M_1\right)<0$ 时矩阵 $M_1$ 的所有特征根都具有负实部。因此系统(4)满足下列式子 $\left(I^{\prime }\left(t\right),W^{\prime }\left(t\right)\right)\to \left(0,0\right),t\to \infty$

利用比较定理 [25] ，可知 $\left(I\left(t\right),W\left(t\right)\right)\to \left(0,0\right),t\to \infty$

再利用渐近自治系统理论 [26] ，可以得出 $\left(S\left(t\right),p\left(t\right),M\left(t\right)\right)\to \left(S^0,p^0,M^0\right),t\to \infty$

所以当 $R_0<1$ 时可知无病平衡点 $P_0$ 是全局吸引的。结合 $P_0$ 的局部稳定性可知，当 $R_0<1$ 时可知系统(3)的无病平衡点 $P_0$ 是全局渐近稳定的。

3.5. 正平衡点的局部稳定

定理3.2：当 $R_0>1$ ，系统的正平衡点 $E^1=\left(S^1,I^1,W^1,p^1,M^1\right)$ 是局部稳定的。

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

${J|}_{E^1}=\left(\begin{array}{ccccc}-{\beta }_1{I}^1-{\beta }_2{W}^1-\alpha p^1& -{\beta }_1{S}^1& -{\beta }_2{S}^1& \epsilon e-\alpha S^1& 0\\ {\beta }_1{I}^1+{\beta }_2{W}^1& {\beta }_1{S}^1-\gamma & {\beta }_2{S}^1& 0& 0\\ 0& k& -\delta & 0& 0\\ 0& 0& 0& -\varphi \left[be+\frac{\mu M^1}{{\left(1+p^1\right)}^2}\right]& -\frac{\varphi \mu p^1}{1+p^1}\\ \alpha p^{\ast }& 0& 0& \alpha S^1-\frac{\mu M^1}{{\left(1+p^1\right)}^2}& -\frac{\varphi \mu p^1}{1+p^1}\end{array}\right)$

其中

$D_1=-{\beta }_1{I}^1-{\beta }_2{W}^1-\alpha p^1<0$

$D_2=\left|\begin{array}{cc}-{\beta }_1{I}^1-{\beta }_2{W}^1-\alpha p^1& -{\beta }_1{S}^1\\ {\beta }_1{I}^1+{\beta }_2{W}^1& {\beta }_1{S}^1-\gamma \end{array}\right|=\gamma \left({\beta }_1{I}^1+{\beta }_2{W}^1\right)-\alpha p^1\left({\beta }_1{S}^1-\gamma \right)>0$

其中： ${\beta }_1{S}^1-\gamma =\frac{-{\beta }_{2}k\gamma }{{\beta }_1\delta +{\beta }_{2}k}<0$

$\begin{array}{c}{D}_3=\left|\begin{array}{ccc}-{\beta }_1{I}^1-{\beta }_2{W}^1-\alpha p^1& -{\beta }_1{S}^1& -{\beta }_2{S}^1\\ {\beta }_1{I}^1+{\beta }_2{W}^1& {\beta }_1{S}^1-\gamma & {\beta }_2{S}^1\\ 0& k& -\delta \end{array}\right|\\ =\left(\delta {\beta }_1{S}^1-k{\beta }_2{S}^1-\delta \gamma \right)\alpha p^1-\left(\delta {\beta }_1{S}^1+\delta \gamma +k{\beta }_2{S}^1\right)\left({\beta }_1{I}^1+{\beta }_2{W}^1\right)<0\end{array}$

其中： $\delta {\beta }_1{S}^1-k{\beta }_2{S}^1-\delta \gamma =\frac{-2{\beta }_{2}k\delta \gamma }{{\beta }_1\delta +{\beta }_{2}k}<0$

$\begin{array}{c}{D}_4=\left|\begin{array}{cccc}-{\beta }_1{I}^1-{\beta }_2{W}^1-\alpha p^1& -{\beta }_1{S}^1& -{\beta }_2{S}^1& \epsilon e-\alpha S^1\\ {\beta }_1{I}^1+{\beta }_2{W}^1& {\beta }_1{S}^1-\gamma & {\beta }_2{S}^1& 0\\ 0& k& -\delta & 0\\ 0& 0& 0& -\varphi \left[be+\frac{\mu M^1}{{\left(1+p^{\ast }\right)}^2}\right]\end{array}\right|\\ =-\varphi \left[be+\frac{\mu M^1}{{\left(1+p^1\right)}^2}\right]D_3>0\end{array}$

$\begin{array}{c}{D}_5=\left|\begin{array}{ccccc}-{\beta }_1{I}^1-{\beta }_2{W}^1-\alpha p^1& -{\beta }_1{S}^1& -{\beta }_2{S}^1& \epsilon e-\alpha S^1& 0\\ {\beta }_1{I}^1+{\beta }_2{W}^1& {\beta }_1{S}^1-\gamma & {\beta }_2{S}^1& 0& 0\\ 0& k& -\delta & 0& 0\\ 0& 0& 0& -\varphi \left[be+\frac{\mu M^1}{{\left(1+p^1\right)}^2}\right]& -\frac{\varphi \mu p^1}{1+p^1}\\ \alpha p^{\ast }& 0& 0& \alpha S^1-\frac{\mu M^1}{{\left(1+p^1\right)}^2}& -\frac{\varphi \mu p^1}{1+p^1}\end{array}\right|\\ =\delta \left({\beta }_1{S}^1-\gamma \right)\left({\beta }_1{I}^1+{\beta }_2{W}^1+\alpha p^1\right)\varphi \left[be+\frac{\mu M^1}{{\left(1+p^1\right)}^2}\right]\frac{\varphi \mu p^1}{1+p^1}<0\end{array}$

其中： ${\beta }_1{S}^1-\gamma =\frac{-{\beta }_{2}k\gamma }{{\beta }_1\delta +{\beta }_{2}k}<0$

得出矩阵 ${J|}_{E^1}$ 是负定矩阵，其特征根都小于0，则正平衡点局部稳定。

4. 最优控制策略

4.1. 建立控制模型

布鲁氏菌严重危害到畜牧业的生产和人类的健康，并且一旦爆发布病，感染羊的被扑杀和流产将直接影响到养殖户的收益问题。因此，必须对布鲁氏菌病加以实施相应的控制策略以抑制疫病的蔓延。对布病的检测可以有效控制其流行，但一般需要很大的控制力度，这关系到检测扑杀处理的费用，包括检测试剂、检测过程、扑杀带来的损失等。为了在达到控制布鲁氏菌病目的的同时降低控制成本，利用最优控制理论对疫病进行时变控制。假设最优控制策略的终端时间T已知，表明该策略是在0~T的时间段内实现的。

在模型(3)中，为了分析控制布鲁氏菌病所需的最佳策略，引入依赖时间的控制量 $u_1\left(t\right)$ ， $u_2\left(t\right)$ ， $u_3\left(t\right)$ ，如下式所示：

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=A+\epsilon e\left(p-p_0\right)-\left(1-u_1\left(t\right)\right){\beta }_{1}SI-{\beta }_{2}SW-\alpha pS\hfill \\ \frac{\text{d}I}{\text{d}t}=\left(1-u_1\left(t\right)\right){\beta }_{1}SI+{\beta }_{2}SW-u_2\left(t\right)\gamma I\hfill \\ \frac{\text{d}W}{\text{d}t}=\left(1-u_3\left(t\right)\right)kI-\delta W\hfill \\ \frac{\text{d}p}{\text{d}t}=\varphi \left(D-bep-\frac{\mu Mp}{1+p}\right)\hfill \\ \frac{\text{d}M}{\text{d}t}=\alpha pS-\frac{\mu Mp}{1+p}\hfill \end{array}$ (5)

其中 $S,I,W,p,M>0$ 且 $\left(u_1\left(t\right),u_2\left(t\right),u_3\left(t\right)\right)\in U$ ，U是控制集，满足

$U=\left\{\left(u_1,u_2,u_3\right)|u_i\left(t\right)\in L^{\infty },0\le u_i\left(t\right)\le 1,t\in \left[0,T\right],i=1,2,3\right\}$

控制项 $u_1\left(t\right)$ ( $0\le u_1\left(t\right)\le 1$ )用于减少易感动物与受感染动物的接触。这可以通过合理化养殖方式，合理建房屋，确保生活区，饲养区和饲料区相互隔离，以及避免来自不同畜群或畜群的动物混杂在一起，并限制动物从受感染地区向其他地区移动来实现。与这种控制相关的成本包括使用零放牧或避免游牧，以及建设市场。

控制项 $u_2\left(t\right)$ ( $0\le u_2\left(t\right)\le 1$ )是养殖户屠杀感染羊做出的努力，为了通过扑杀成功控制布鲁氏菌病，需要努力识别受感染的动物并坚持进行这项工作。与扑杀相关的成本可能与立即宰杀所有测试呈阳性的动物有关。这种方法需要业主的合作，补偿通常是必要的。在检测不太可靠的地方，对绵羊和山羊来说，检测和屠宰政策的成本可能更高。为了成功控制布鲁氏菌病，政府部门必须适当规划、协调和提供资源。

控制项 $u_3\left(t\right)$ ( $0\le u_3\left(t\right)\le 1$ )是减少布鲁氏菌病对环境排放所需的努力。这可以通过焚烧被感染动物的流产胎儿、组织和排泄物来实现。与这种控制相关的成本包括清洁和消毒。

建立目标函数使得养殖户损失最小，即将羊群中受感染羊的数量和环境中布鲁氏菌浓度降到最低，同时将与三个控制量相关的成本尽可能降到最低。

目标函数

$J\left(u_1,u_2,u_3\right)={\displaystyle \underset{0}{\overset{T}{\int }}\left[A_{1}I+A_{2}W+\frac{B_1}{2}{u}_1^2+\frac{B_2}{2}{u}_2^2+\frac{B_3}{2}{u}_3^2\right]\text{d}t}$ (6)

其中 $A_1$ 和 $A_2$ 是对应受感染动物的数量和环境中布鲁氏菌浓度的权重因子， $B_1,B_2,B_3$ 是分别平衡与三个控制策略 $u_1,u_2,u_3$ 相关的成本因素的权重因子。

4.2. 控制解

给出控制解 $\left(u_1^{*},u_2^{*},u_3^{*}\right)$ 满足的必要条件，控制系统(5)下的目标函数(6)极小化问题通过庞特里亚金极小值原理转化为控制集U下的哈密顿量H极小化问题，构造哈密顿量为：

$\begin{array}{c}H=A_{1}I+A_{2}W+\frac{B_1}{2}{u}_1^2+\frac{B_2}{2}{u}_2^2+\frac{B_3}{2}{u}_3^2\\ +{\lambda }_1\left[A+\epsilon e\left(p-p_0\right)-\left(1-u_1\right){\beta }_{1}SI-{\beta }_{2}SW-\alpha pS\right]\\ +{\lambda }_2\left[\left(1-u_1\right){\beta }_{1}SI+{\beta }_{2}SW-u_2\gamma I\right]+{\lambda }_3\left[\left(1-u_3\right)kI-\delta W\right]\\ +{\lambda }_4\left[\varphi \left(D-bep-\frac{\mu Mp}{1+p}\right)\right]+{\lambda }_5\left[\alpha pS-\frac{\mu Mp}{1+p}\right]\end{array}$ (7)

其中 ${\lambda }_i\left(i=1,2,3,4,5\right)$ 为状态变量 $S,I,W,p,M$ 所对应的伴随变量。

定理4.1：控制系统(5)下使得目标函数(6)极小化的控制组合解 $\left(u_1^{*},u_2^{*},u_3^{*}\right)$ 和对应的状态变量 $\left(S^{\ast },I^{\ast },W^{\ast },p^{\ast },M^{\ast }\right)$ 存在，则伴随变量满足：

$\{\begin{array}{l}\frac{\text{d}{\lambda }_1}{\text{d}t}={\lambda }_1\left(1-u_1^{\ast }\right){\beta }_1{I}^{\ast }+{\lambda }_1{\beta }_2{W}^{\ast }+{\lambda }_1\alpha p^{\ast }-{\lambda }_2\left(1-u_1^{\ast }\right){\beta }_1{I}^{\ast }-{\lambda }_2{\beta }_2{W}^{\ast }-{\lambda }_5\alpha p^{\ast }\hfill \\ \frac{\text{d}{\lambda }_2}{\text{d}t}=-A_1+{\lambda }_1\left(1-u_1^{\ast }\right){\beta }_1{S}^{\ast }-{\lambda }_2\left(1-u_1^{\ast }\right){\beta }_1{S}^{\ast }+{\lambda }_2{u}_2^{\ast }\gamma -{\lambda }_3\left(1-u_3^{\ast }\right)k\hfill \\ \frac{\text{d}{\lambda }_3}{\text{d}t}=-A_2+{\lambda }_1{\beta }_2{S}^{\ast }-{\lambda }_2{\beta }_2{S}^{\ast }+{\lambda }_3\delta \hfill \\ \frac{\text{d}{\lambda }_4}{\text{d}t}=-{\lambda }_1\epsilon e+{\lambda }_1\alpha S^{\ast }+{\lambda }_4\varphi be+{\lambda }_4\varphi \frac{\mu M^{\ast }}{{\left(1+p^{\ast }\right)}^2}-{\lambda }_5\alpha S^{\ast }+{\lambda }_5\frac{\mu M^{\ast }}{{\left(1+p^{\ast }\right)}^2}\hfill \\ \frac{\text{d}{\lambda }_5}{\text{d}t}={\lambda }_4\varphi \frac{\mu p^{\ast }}{1+p^{\ast }}+{\lambda }_5\frac{\mu p^{\ast }}{1+p^{\ast }}\hfill \end{array}$

并且具有横截性条件 ${\lambda }_i\left(T\right)=0,i=1,2,3,4,5$

最优控制变量的解为

$\{\begin{array}{l}{u}_1^{\ast }=\min \left\{1,\max \left\{0,\frac{\left({\lambda }_2-{\lambda }_1\right){\beta }_1{S}^{\ast }{I}^{\ast }}{B_1}\right\}\right\}\hfill \\ u_2^{\ast }=\min \left\{1,\max \left\{0,\frac{{\lambda }_2\gamma I^{\ast }}{B_2}\right\}\right\}\hfill \\ u_3^{\ast }=\min \left\{1,\max \left\{0,\frac{{\lambda }_{3}k{I}^{\ast }}{B_3}\right\}\right\}\hfill \end{array}$

证明：根据Pontryagin极小值原理可以得到伴随系统和横截性条件。由哈密顿量函数(7)对状态变量 $S,I,W,p,M$ 求导得到：

$\frac{\text{d}{\lambda }_1}{\text{d}t}=-\frac{\partial H}{\partial S},{\lambda }_1\left(T\right)=0$ ， $\frac{\text{d}{\lambda }_2}{\text{d}t}=-\frac{\partial H}{\partial I},{\lambda }_2\left(T\right)=0$ ， $\frac{\text{d}{\lambda }_3}{\text{d}t}=-\frac{\partial H}{\partial W},{\lambda }_3\left(T\right)=0$

$\frac{\text{d}{\lambda }_4}{\text{d}t}=-\frac{\partial H}{\partial p},{\lambda }_4\left(T\right)=0$ ， $\frac{\text{d}{\lambda }_5}{\text{d}t}=-\frac{\partial H}{\partial M},{\lambda }_5\left(T\right)=0$

为了求解控制组 $\left(u_1^{*},u_2^{*},u_3^{*}\right)$ ，得到哈密顿量H的极值解，最优控制方程为 $\frac{\partial H}{\partial u_i}=0,i=1,2,3$

由于控制变量 $u_i\in \left[0,1\right]$ ，因此得到：

$\{\begin{array}{l}{u}_1^{\ast }=\min \left\{1,\max \left\{0,\frac{\left({\lambda }_2-{\lambda }_1\right){\beta }_1{S}^{\ast }{I}^{\ast }}{B_1}\right\}\right\}\hfill \\ u_2^{\ast }=\min \left\{1,\max \left\{0,\frac{{\lambda }_2\gamma I^{\ast }}{B_2}\right\}\right\}\hfill \\ u_3^{\ast }=\min \left\{1,\max \left\{0,\frac{{\lambda }_{3}k{I}^{\ast }}{B_3}\right\}\right\}\hfill \end{array}$

5. 数值模拟

以河南省为例进行模拟，分别对2016~2021年河南羊的活重价格和2016~2020人间布鲁氏菌病的累计病例数进行了模拟，以验证模型(3)的合理性，进而预测布鲁氏菌病的发展。

5.1. 模拟价格数据

模型(1)用于模拟羊的活重价格和估计参数 $\epsilon ,\varphi ,\mu$ ，通过查找中国统计年鉴 [27] 和《国家农产品成本效益数据汇编》 [28] ，得到2016~2021年河南省羊群的平均规模为19,260,000头，2016年每公斤价格为27.6278元，所以 $N\left(0\right)=19260000$ ， $p\left(0\right)=27.6278$ 。定义 $P\left(t\right)$ 为出栏羊的实际价格，定义 $P\left(t\right)$ 为出栏羊的理论

价格，并且 $\frac{\text{d}p}{\text{d}t}=\varphi \left(D-bep-\frac{\mu Mp}{1+p}\right)$ ，则目标函数为：

$\min {\displaystyle \underset{t}{\sum }{\left[P\left(t\right)-p\left(t\right)\right]}^2}$

2016~2021年的实际价格数据见表1，运用MATLAB，采用最小二乘估计，得到数据拟合结果和参数值 $\epsilon =70.28,\varphi =1.27\times {10}^{-6},\mu =1.242$ ，见图1中浅绿色阴影区域利用MCMC方法显示模拟10000次的95%置信区间，蓝色符号标记了实际价格数据。

表1. 河南2016~2021年的实际数据，羊价格数据来自中国统计出版社出版的《国家农产品成本效益数据汇编》

图1. 河南地区羊价格的拟合结果

5.2. 模拟人患布鲁氏菌病累计病例数据

模型(2)用于模拟人患布鲁氏菌病的累积病例和估计参数 ${\beta }_1$ ， ${\beta }_2$ ， ${\beta }_h$ ， ${\beta }_w$ ，假设初始值为 $S\left(0\right)=\text{18238868}$ ， $I\left(0\right)=1021132$ ， $W\left(0\right)=78240$ ， $p\left(0\right)=27.6278$ ， $M\left(0\right)=8296000$ ， $S_h\left(0\right)=8852559$ ， $I_h\left(0\right)=3993$ ，目标函数为：

$\min {\displaystyle \underset{t}{\sum }{\left[Y\left(t\right)-y\left(t\right)\right]}^2}$

其中 $Y\left(t\right)$ 为人患布鲁氏菌病的实际累积病例， $y\left(t\right)$ 为人患布鲁氏菌病的理论累积病例。设 $y\left(0\right)=3993$ 并

且 $y\left(t\right)$ 满足 $\frac{\text{d}y}{\text{d}t}={\beta }_h{S}_{h}I+{\beta }_w{S}_{w}W$ 。

与上述模拟相同，2016~2020年的实际人患病数据见表2，运用MATLAB，采用最小二乘估计，得到了数据拟合结果和参数值 ${\beta }_1=7.4\times {10}^{-11}$ ， ${\beta }_2=9.2\times {10}^{-11}$ ， ${\beta }_h=1.64\times {10}^{-10}$ ， ${\beta }_w=1.72\times {10}^{-10}$ 。

见图2中浅绿色的阴影区域利用MCMC方法显示模拟10,000次结果的95%置信区间，蓝色记号标记了人类布鲁氏菌病实际累积病例。所以的参数值见表3。

模型(3)的所有参数都是已知见表3。因此可以得到系统的基本再生数为 $R_0=1.6338>1$ ，此时布病时呈上升趋势的，如果不及时改变控制措施，将会导致布鲁氏菌病在羊群和人类中出现更大规模的爆发。

表2. 河南2016~2020年的实际数据，人患布鲁氏菌病数据来自国家疾病通报监测系统(NNDSS)

图2. 河南地区累计人患病的拟合结果

表3. 参数值和来源

6. 最优控制策略分析

最优控制策略的理论分析已在第4节中给出。本节根据三个控制项目的七种不同组合给出数值分析结果，为河南省提供最佳策略。假设目标函数(6)中的终点时间为 $T=10$ ，表示该段控制策略的实施周期为2020年至2030年这十年，目标函数涉及的权重因子在此做一个假设 $A=100$ ， $B_1=5000$ ， $B_2=10000$ ， $B_3=3000$ 。面对最优控制问题，使用四阶龙格–库塔的向前向后扫描法来求解。七种可能的组合策略分为三种情景，具体详情如下：

6.1. 情形A (使用单项控制策略)

策略一 $\left(u_1\ne 0,u_2=0,u_3=0\right)$ ：采用单一对照 $u_1\left(t\right)$ ，即减少易感动物与受感染动物的接触；

策略二 $\left(u_1=0,u_2\ne 0,u_3=0\right)$ ：采用单一对照 $u_2\left(t\right)$ ，即加大屠杀感染羊的力度；

策略三 $\left(u_1=0,u_2=0,u_3\ne 0\right)$ ：采用单一对照 $u_3\left(t\right)$ ，即降低布鲁氏菌病对环境的排放；

见图3，其显示了情形A中三种策略的数值结果。如图3(a)所示，三种策略与目前的控制相比受感染的羊只数量都大大减少，其中策略一的效果最佳。如图3(b)所示，三种策略与目前的控制相比环境中

图3. 情形A的数值结果。(a) 感染羊的数量变化图，(b) 环境中布鲁氏菌的数量变化图，(c) 养殖户损失图，(d) 控制项变化图

布鲁氏菌浓度的都大大减少，其中策略三的效果最佳。对养殖户来说，最重要的减少损失，如图3(c)所示，三种策略与目前的控制相比养殖户的损失都有所减少，其中策略三的效果最佳。图3(d)是控制项的变化图。考虑到所有因素，如果养殖户采用单一控制策略，最合理的策略是采用策略三。

6.2. 情形B (使用双项控制策略)

策略四 $\left(u_1\ne 0,u_2\ne 0,u_3=0\right)$ ：采取双重控制 $u_1\left(t\right)$ 和 $u_2\left(t\right)$ ：1) 减少易感动物与受感染动物的接触；2) 加大屠杀感染羊的力度；

策略五 $\left(u_1=0,u_2\ne 0,u_3\ne 0\right)$ ：采取双重对照 $u_2\left(t\right)$ 和 $u_3\left(t\right)$ ：1) 加大屠杀感染羊的力度；2) 降低布鲁氏菌病对环境的排放；

策略六 $\left(u_1\ne 0,u_2=0,u_3\ne 0\right)$ ：采取双重对照 $u_1\left(t\right)$ 和 $u_3\left(t\right)$ ：1) 减少易感动物与受感染动物的接触；2) 降低布鲁氏菌病对环境的排放；

见图4，其显示了情形B中三种策略的数值结果。如图4(a)所示，三种策略与目前的控制相比受感染的羊只数量都大大减少，其中策略四的效果最佳。如图4(b)所示，三种策略与目前的控制相比环境中布鲁氏菌浓度的都大大减少，其中策略六的效果最佳。如图4(c)所示，三种策略与目前的控制相比养殖户的损失都有所减少，其中策略六的效果最佳。图4(d)是控制项的变化图。考虑到所有因素，如果养殖户采用两种控制策略的组合，最合理的策略是采用策略六。

(a) (b) (c) (d) (e) (f)

Figure 4. Numerical results for scenario B. (a) Plot of changes in the number of infected sheep, (b) Plot of changes in the number of Brucella abortus in the environment, (c) Plot of farmer losses, (d) Plot of changes in control terms for strategy 4, (e) Plot of changes in control terms for strategy 5, and (f) Plot of changes in control terms for strategy 6