基于动态价格因素下的新疆布鲁氏菌病最优控制策略研究

Research on Optimal Brucellosis Control Strategy in Xinjiang Based on Dynamic Price Factors

摘要: 近年来，新疆地区动物感染布鲁氏菌病的发病率以及新发人间病例呈持续上升趋势，这不仅对畜牧业发展构成了显著影响，而且给公共卫生安全带来了严峻挑战。因此，加强布鲁氏菌病风险评估与防控工作，对于保障畜牧业的可持续发展具有关键意义。然而畜牧业的发展不仅受到动物疫病因素的影响，市场价格与供需关系等经济因素同样发挥重要作用。基于上述背景，本文将市场价格因素与布鲁氏菌病传播特点相结合，建立了受动态价格影响的布鲁氏菌病传播动力学模型。首先对此模型的无病平衡点及基本再生数进行了求解和验证。其次通过对传染病模型施加三种控制措施，并以实现最小传染规模和最低防控成本为目标，构建了目标函数。最后在数值模拟中，选取新疆维吾尔自治区2011至2023年的活羊出栏价格数据以及人间布鲁氏菌病病例数，对模型进行参数估计，验证了其合理性。同时，通过参数敏感性分析和数值实验，得出结论：在保持牲畜种群规模不扩张的前提下，选择出栏牲畜的策略有利于传染病的控制；而当三种控制策略联合施行时，实现疾病防控效果最优并且疾病控制成本最低。

Abstract: In recent years, the incidence of brucellosis infection in animals and new human cases in Xinjiang has been steadily increasing, which not only has a significant impact on the development of animal husbandry, but also poses a serious challenge to public health safety. Therefore, strengthening the risk assessment, prevention and control of brucellosis is essential to ensure the sustainable development of animal husbandry. However, the development of animal husbandry is not only influenced by animal diseases, but also by economic factors such as market prices and supply and demand relations, which also play an important role. Based on the above background, this paper combines the market price factor with the characteristics of brucellosis transmission and establishes a model of brucellosis transmission dynamics affected by dynamic price. Firstly, the disease-free equilibrium point and the basic regeneration number of this model are solved and verified. Secondly, an objective function was constructed by imposing three control measures on the transmission model and aiming to achieve the minimum level of transmission and the minimum cost of prevention and control. Finally, in the numerical simulation, the live sheep slaughter price data from 2011 to 2023 and the number of human brucellosis cases in Xinjiang Uygur Autonomous Region were selected to estimate the parameters of the model and verify its reasonableness. At the same time, through parameter sensitivity analysis and numerical experiments, it was concluded that the strategy of choosing livestock selling is conducive to the control of infectious diseases under the premise that the size of the livestock population does not increase; and when the three control strategies are implemented together, the optimal effect of disease prevention and control and the lowest cost of disease control can be achieved.

1. 引言

布鲁氏菌病是一种由布鲁氏菌引起的急性人畜共患病[1]，作为一种古老而顽固的人畜共患病，它在全球范围内对公共卫生安全和畜牧业发展构成威胁[2] [3]。新疆作为我国五大牧区之一，流行牛羊布鲁氏菌，且近年来疫情趋势严峻[4]，2021~2023年间新疆人间布鲁菌病疫情波及14个地市，95个县区，累计报告病例18,094例，年均发病率达到23.33/10万，对人民群众身心健康与畜牧业发展造成严重威胁[5]。对于畜牧业而言，布鲁氏菌病不仅造成牲畜个体健康损伤，还对养殖效益产生不利影响。考虑到畜牧业是新疆地区的传统支柱产业以及严峻的疫情形势，加强新疆地区的布鲁氏菌病研究对于推进畜牧业的长久发展具有重要的现实意义。

近年来，许多学者对布鲁氏菌病保持高度关注，根据布病传播机理，采用动力学方法建立模型深入研究。Liang [6]等对布鲁氏菌病的潜伏期进行分析，证明平衡点的稳定性，并进一步求出系统的阈值。Hou [7]等基于布鲁氏菌病的传播特点，提出SEIB模型，结果表明，平衡点的全局稳定性完全取决于基本再生数，时滞对平衡点的稳定性无影响。Hu [8]等结合布鲁氏菌病非局部传播以及季节性扩散的特点建立季节布鲁氏菌病的传播模型，筛查随访了多名牧羊人，测量出布鲁氏菌病流行区高暴露人群的血清流行率，并报告血清阳性无症状受试者的结果和持续时间。Lolika [9]等结合两个离散的时滞以及对感染动物的扑杀建立了一个布鲁氏菌病的传播模型。其中第一个时滞表示潜伏期，第二个时滞表示检测和扑杀感染性动物所需时间。通过分析和数值方法确定了模型稳态的可行性和稳定性。结果表明，无论是分析还是数值，都表明这两个延迟会使系统不稳定，且可以通过Hopf分叉产生周期解。Zhang [10]等结合浙江省布鲁氏菌病的病理学、报告的阳性数据和饲养方法等建立了动态模型，拟合奶牛实际疾病情况，预测布鲁氏菌感染趋势。结果表明，除了已采取的措施以外，感染地区的消毒时间定为每周两次并管理奶牛出生率，该措施可以更有力地控制布鲁氏菌病的传播。Feng [11]等人基于染病人数建立动力学模型，得出疾病的最优控制解，分析检测行为对布鲁氏菌病的传播影响，发现当参数满足一定条件时，系统会出现周期解。Li [12] [13]等人考虑到年龄结构，建立了包含多年龄阶段的布鲁氏菌病的传播模型。Nie [14]等人结合吉林省布鲁氏菌病感染情况，提出了一种具有外部转移量的动态模型SEIV描述布鲁氏菌病在奶牛中的传播，并结合吉林省近20年的数据进行参数估计以及数值模拟。

同时，许多学者也针对布鲁氏菌病在新疆地区产生的严重威胁展开了相关研究。Liu [15]等人利用血清学检测，分子生物学检测和问卷调查的方法，统计新疆布病流行状况，通过Logistics回归法分析验证定期监测以及制定合理的免疫防控的重要性。Zheng [16]等人以青河县乡镇与养羊场作为调查单位，对2016~2021年采集的牛羊血清进行布鲁氏菌病监测。结果显示，不同区域以及不同年龄段的牛羊布鲁氏菌病检测阳性率不同。Liao [17]等人通过数据分析发现，2011~2021年期间，牛布氏杆菌感染阳性率为0.77%~6.76%，平均感染率为2.43%；羊布氏杆菌感染阳性率为0.12%~9.50%，平均感染率为1.88%。

值得注意的是，养殖收益不止受到动物疫病的影响，同样受到市场供需关系影响。牲畜的产量及需求量通过对市场价格产生直接影响，进而影响养殖决策以及养殖利润。当市场价格显著偏高时，消费者的购买意愿会随之下降，从而导致需求减少；反之，若价格过低，则可能削弱生产者的盈利能力，进而阻碍养殖业的发展。基于此，深入探讨商品价格波动的内在机制具有重要的理论和实践意义。Rahman [18]等人使用区间参数法，开发了一个非确定性库存模型，其中的市场需求率被视作为区间值，且取决于销售价格。许多学者[19] [20]通过建立线性需求函数、非线性需求函数等建立微分方程模型，探究市场供需对价格的影响。Wang [21]等考虑动态价格因素建立模型，Chen [22]等考虑最大养殖量建立模型。但以往学者并未考虑将动态的市场价格与养殖场的承载能力进行有机结合，并基于新疆地区养殖情景进行具体分析，而这些因素对于制定新疆地区的布病控制策略是至关重要的。

本文以新疆羊群养殖为例，根据市场价格以及养殖规模的动态变化构建动力学模型，考虑到布鲁氏菌病的传播特点，加入控制策略寻求使得养殖户成本损失最小的策略。本文的结构如下：在文章的第2节中阐述模型构建过程。第3节中对模型进行动力学分析，求解出无病平衡点以及正平衡点，并验证其局部稳定性。第4节中提出最优控制模型，利用庞特里亚金极值原理给出最优控制的理论解。在第5节中，选取新疆2011~2023年的活羊出栏价格以及人间布鲁氏菌病的累计病例进行参数估计。最优控制策略的分析及结论将在第6节中重点讨论。最后，在第7节中对文章的主要工作以及内容进行阐述。

2. 模型建立

2.1. 基础模型

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

定义$N\left(t\right)$ 为时刻$t$ 时羊类种群的数量。种群在未开发，按照自然状态进行生长时，选择采用Logistic增长模型进行描述。在羊群资源不断进行开发时，羊群将按照出栏率$\theta$ 投入到市场当中，羊群的存栏量增长趋势为：

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

对上式观察易知，若有$\theta >r$ 时，羊群的数量逐渐减少，趋于灭绝；反之羊群的数量将趋于稳定在

$N^{*}=k\left(1-\frac{\theta }{r}\right)$ 。为考虑更切合实际的情形，本文中假设考虑模型中始终成立：$r>\theta$ 。要注意的是，表达

式中的最大养殖量$k$ 为一个随着时间进行变化的量，而非一个常数。

根据经济学理论，供需平衡时市场价格趋于稳定；而当供需失衡时，市场价格便会发生波动，具体表现为：供应量超过需求量时价格呈下降趋势，而供应不足以满足需求时价格则倾向于上升。假设需求函数为$D\left(p\right)$ ，基于羊肉将根据需求下降的背景，将$D\left(p\right)$ 设立为关于价格的单调递减函数，得出的价格模型为$D\left(p\right)=D-bep$ ，市场需求不可能为负数，因此始终成立$D\left(p\right)=D-bep>0$ 。与出栏量$\theta N\left(t\right)$ 进行结合后得到动态价格模型：

$\frac{\text{d}p}{\text{d}t}=\phi \left(D-bep-\theta N\left(t\right)\right)$

市场价格波动对养殖业的生产决策具有显著影响。当市场价格处于较高水平时，较高的盈利预期将促使养殖业者增加投入并扩大养殖规模以实现利润最大化；而在价格下行阶段，由于盈利空间受到挤压，养殖业者则趋向于缩减投入并缩小规模，从而减少潜在损失。

本文设定了一个用于刻画市场价格水平的价值阈值。考虑到该阈值的上限不得超过市场价格的最高值，因此，其数值应限定在一个特定区间内浮动。以下将刻画最大养殖量与市场价格变化之间的关系，在某一时间段内最大养殖量将随着市场价格与价格阈值之差进行变动，其中比例系数$\epsilon >0$ ，表述形式为：

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

考虑到自然资源等限制因素，采用参数$k_m$ 刻画养殖规模增长的上限；变量$p$ 表示随时间变化的价格。

将上述种群增长、动态价格以及最大养殖量的变动结合起来进行描述，可以得到以下模型：

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

2.2. 布鲁氏菌病传播下的养殖模型

在无传染病情形下，已建立基础养殖模型(1)。考虑到布鲁氏菌病在人群与羊群之间的传播机制，进一步构建综合考虑羊群、人群、价格及养殖规模的动力学模型。具体而言，将羊群细分为易感羊群$S$ 和感染羊群$I$ ，引入变量$W$ 以描述环境中布鲁氏菌的浓度，并用$p$ 表示活羊出栏时的市场价格；同时，将人群划分为易感人群$S_h$ 和患病人群$I_h$ 。

根据布鲁氏菌的传播特点，做出以下假设：

(1) 假设染病羊群不再进行繁殖，故只考虑布鲁氏菌病接触感染者或接触受污环境时的情形[23]。

(2) 假设羊群增长满足Logistic增长，且新增长的羊只进入易感羊群仓室。

(3) 检验结果为阳性的羊将立即进行扑杀处理，只有易感仓室中的羊群进行出栏以供应市场。

基于以上模型假设，我们构造了如下形式的满足Logistic增长的布病传播养殖模型：

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=rS\left(1-\frac{S+I}{k}\right)-{\beta }_{1}SI-{\beta }_{2}SW-\theta S\\ \frac{\text{d}I}{\text{d}t}={\beta }_{1}SI+{\beta }_{2}SW-dI\\ \frac{\text{d}W}{\text{d}t}=\eta I-\delta W\\ \frac{\text{d}k}{\text{d}t}=\epsilon \left(p-p_0-cdI\right)k\left(1-\frac{k}{k_m}\right)\\ \frac{\text{d}p}{\text{d}t}=\phi \left(D-bep-\theta S\right)\end{array}$ (2)

其中，模型中的$r$ 表示羊群的内禀增长率，考虑到布鲁氏菌病在牲畜中的传播途径为接触染病个体与接触存在病菌体的环境两种方式，采用双线性发生率${\beta }_{1}SI+{\beta }_{2}SW$ 来描述传播过程，其中${\beta }_1$ 和${\beta }_2$ 分别为两种感染方式的感染率。此外，采用$\theta$ 表示易染个体进入到市场中的出栏率，故牲畜出栏量为$\theta S$ 。当易染个体被感染时，牲畜将从易感仓室$S$ 进入感染仓室$I$ ，在对羊群进行抽检的过程中如果发现染病个体，将采用扑杀率为$d$ 的措施对患病羊群进行扑杀处理。由于环境中的布鲁氏菌病来源为患病个体的排放，并且其浓度随着时间的推移而减弱，此处采用$\eta$ 表示受感染牲畜单位时间内排放布鲁氏菌到环境中的排放率，$\delta$ 表示环境中的布鲁氏菌的衰减率。在养殖过程中，市场价格受到多方因素影响，进而又反过来影响养殖规模。参数$c$ 表示扑杀感染牲畜时所造成的损失，进而对价格产生影响。

由于布鲁氏菌病的传播具有家畜单向传播给人且人间不传播的特点，并且人间的布鲁氏菌病感染方式同样主要考虑接触受感染牲畜以及环境中的布鲁氏菌，采用双线性发生率${\beta }_h{S}_{h}I+{\beta }_w{S}_{h}W$ 来描述传播过程，其中${\beta }_h$ 和${\beta }_w$ 分别为两种感染方式的感染率。以常数$A_h$ 的输入表示人口出生，$d_h$ 为人的自然死亡率，$m$ 为疾病康复率，得到方程组：

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=rS\left(1-\frac{S+I}{k}\right)-{\beta }_{1}SI-{\beta }_{2}SW-\theta S\\ \frac{\text{d}I}{\text{d}t}={\beta }_{1}SI+{\beta }_{2}SW-dI\\ \frac{\text{d}W}{\text{d}t}=\eta I-\delta W\\ \frac{\text{d}k}{\text{d}t}=\epsilon \left(p-p_0-cdI\right)k\left(1-\frac{k}{k_m}\right)\\ \frac{\text{d}p}{\text{d}t}=\phi \left(D-bep-\theta S\right)\\ \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\\ \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\end{array}$ (3)

由于整个羊群被划分为易感个体与染病个体两个仓室，故满足$N=S+I$ ，系统的前两个方程相加后可以得到：

$\frac{\text{d}\left(S+I\right)}{\text{d}t}=rS\left(1-\frac{S+I}{k}\right)-\theta S-dI<rS\left(1-\frac{S+I}{k}\right)$

总人数为$N_h=S_h+I_h$ ，系统的最后两个方程相加得到：

$\frac{\text{d}\left(S_h+I_h\right)}{\text{d}t}=A_h-d_h{N}_h$

得到布鲁氏菌病传播系统(2)的正不变集：

$\Gamma =\left\{\left(S,I,W,k,p,S_h,I_h\right)\in R_7^{+},0\le S+I\le k,0\le W\le \frac{\eta k}{\delta },0<k<k_m,0<p<\frac{D}{be},0\le \left(S_h+I_h\right)\le \frac{A_h}{d_h}\right\}$ 。

3. 模型动力学分析

由于系统中的后两个方程独立于前五个方程，因此对系统的动力学分析可以直接分析系统(2)。

3.1. 无病平衡点

设系统(2)中所有方程的右端为0，发现系统总是存在着两个不同的无病平衡点：

$E_1=\left(\overline{S_1^0},0,0,\overline{k_1^0},\overline{p_1^0}\right)$ ，$E_2=\left(\overline{S_2^0},0,0,\overline{k_2^0},\overline{p_2^0}\right)$ 。

其中：$\overline{S_1^0}=\frac{\left(r-\theta \right)k_m}{r}$ ，$\overline{k_1^0}=k_m$ ，$\overline{p_1^0}=\frac{Dr-\theta \left(r-\theta \right)k_m}{ber}$ ；$\overline{S_2^0}=\frac{D-be{p}_0}{\theta }$ ，$\overline{k_2^0}=\frac{r\left(D-be{p}_0\right)}{\theta \left(r-\theta \right)}$ ，$\overline{p_2^0}=p_0$ 。

3.2. 基本再生数

根据[24] [25]中基本再生数和下一代矩阵的定义，考虑辅助系统：

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

可以得到：

$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}{d}& 0\\ \frac{\eta }{\delta d}& \frac{1}{\delta }\end{array}\right]$

系统的基本再生数为：

$R_0=\rho \left(F{V}^{-1}\right)=\frac{\left({\beta }_1\delta +{\beta }_2\eta \right)S_0}{\delta d}$ 。

3.3. 正平衡点

设系统(2)中所有方程的右端为0，当$R_0>1$ 时，可以得到系统(2)的正平衡点。下面将根据最大养殖量的不同情况对正平衡点进行讨论。

当有$k=k_m$ 时，此时有$S_0=\overline{S_1^0}$ ，可得

$R_{01}=\frac{{\beta }_1\delta +{\beta }_2\eta }{\delta d}\frac{\left(r-\theta \right)k_m}{r}$

可得系统的正平衡点$E_3=\left(S_1^{*},I_1^{*},W_1^{*},k_m,p_1^{*}\right)$ ，其中：

$W_1^{*}=\frac{\eta I_1^{*}}{\delta }$ ，$S_1^{*}=\frac{\delta d}{{\beta }_1\delta +{\beta }_2\eta }$ ，$p_1^{*}=\frac{D\left({\beta }_1\delta +{\beta }_2\eta \right)-\theta \delta d}{be\left({\beta }_1\delta +{\beta }_2\eta \right)}$ 。

$I$ 满足方程式：

$f\left(I\right)=-\frac{r\delta +{\beta }_1{k}_m\delta +{\beta }_2{k}_m\eta }{k_m\delta }I+\frac{\left(R_{01}-1\right)\delta dr}{k_m\left({\beta }_1\delta +{\beta }_2\eta \right)}$ ,

易知当$k=k_m$ 且满足$R_{01}>1$ 时，方程存在唯一的正解

$I_1^{*}=\frac{\left(R_{01}-1\right){\delta }^{2}dr}{\left({\beta }_1\delta +{\beta }_2\eta \right)\left(r\delta +{\beta }_1{k}_m\delta +{\beta }_2{k}_m\eta \right)}$ 。

此时有：

$W_1^{*}=\frac{\left(R_{01}-1\right)\delta dr\eta }{\left({\beta }_1\delta +{\beta }_2\eta \right)\left(r\delta +{\beta }_1{k}_m\delta +{\beta }_2{k}_m\eta \right)}$ 。

当$k\ne k_m$ 时，有$S_0=\overline{S_2^0}$ ，可得：

$R_{02}=\frac{{\beta }_1\delta +{\beta }_2\eta }{\delta d}\frac{D-be{p}_0}{\theta }$ 。

可得系统(2)的正平衡点$E_4=\left(S_2^{*},I_2^{*},W_2^{*},k_2^{*},p_0\right)$ ，其中：

$W_1^{*}=\frac{\eta I_1^{*}}{\delta }$ ，$S_1^{*}=\frac{\delta d}{{\beta }_1\delta +{\beta }_2\eta }$ 。

$I$ 满足方程式：

$g\left(I\right)=-cdI+\frac{\left(R_{02}-1\right)\theta \delta d}{be\left({\beta }_1\delta +{\beta }_2\eta \right)}$

易知当$k\ne k_m$ 且满足$R_{02}>1$ 时，方程存在唯一的正解：

$I_2^{*}=\frac{\left(R_{02}-1\right)\theta \delta }{cbe\left({\beta }_1\delta +{\beta }_2\eta \right)}$ 。

此时有

$W_2^{*}=\frac{\left(R_{02}-1\right)\theta \eta }{cbe\left({\beta }_1\delta +{\beta }_2\eta \right)}$ ，$k_2^{*}=\frac{r\delta \left[cbed+\left(R_{02}-1\right)\theta \right]}{\left({\beta }_1\delta +{\beta }_2\eta \right)\left[cbe\left(r-\theta \right)-\left(R_{02}-1\right)\theta \right]}$ 。

3.4. 无病平衡点的局部稳定性

定理3.1：价格阈值满足$p_0<f\left(\theta \right)$ ，$R_{01}<1$ 时，无病平衡点$E_1$ 是局部渐近稳定的，其中 $f\left(\theta \right)=\frac{k_m}{ber}{\theta }^2-\frac{k_m}{be}\theta +\frac{D}{be}$ 。

证明：无病平衡点$E_1$ 的雅可比矩阵为

$J\left(E_1\right)=\left[\begin{array}{ccccc}-r+\theta & \frac{\left(\theta -r\right)\left(r+{\beta }_1{k}_m\right)}{r}& \frac{-{\beta }_2\left(r-\theta \right)k_m}{r}& \frac{{\left(r-\theta \right)}^2}{r}& 0\\ 0& \frac{{\beta }_1\left(r-\theta \right)k_m}{r}-d& \frac{{\beta }_2\left(r-\theta \right)k_m}{r}& 0& 0\\ 0& \eta & -\delta & 0& 0\\ 0& 0& 0& \epsilon \left(p_0-\frac{Dr-\theta \left(r-\theta \right)k_m}{ber}\right)& 0\\ -\phi \theta & 0& 0& 0& -\phi be\end{array}\right]$

其中

$D_1=-r+\theta <0$

$D_2=\left|\begin{array}{cc}-r+\theta & \frac{\left(\theta -r\right)\left(r+{\beta }_1{k}_m\right)}{r}\\ 0& \frac{{\beta }_1\left(r-\theta \right)k_m}{r}-d\end{array}\right|=\left(\frac{{\beta }_1\left(r-\theta \right)k_m}{r}-d\right)\left(-r+\theta \right)>0$

$D_3=\left|\begin{array}{ccc}-r+\theta & \frac{\left(\theta -r\right)\left(r+{\beta }_1{k}_m\right)}{r}& \frac{-{\beta }_2\left(r-\theta \right)k_m}{r}\\ 0& \frac{{\beta }_1\left(r-\theta \right)k_m}{r}-d& \frac{{\beta }_2\left(r-\theta \right)k_m}{r}\\ 0& \eta & -\delta \end{array}\right|=\left(-r+\theta \right)\left[-\delta \left(\frac{{\beta }_1\left(r-\theta \right)k_m}{r}-d\right)-\eta \frac{{\beta }_2\left(r-\theta \right)k_m}{r}\right]<0$

$\begin{array}{c}{D}_4=\left|\begin{array}{cccc}-r+\theta & \frac{\left(\theta -r\right)\left(r+{\beta }_1{k}_m\right)}{r}& \frac{-{\beta }_2\left(r-\theta \right)k_m}{r}& \frac{{\left(r-\theta \right)}^2}{r}\\ 0& \frac{{\beta }_1\left(r-\theta \right)k_m}{r}-d& \frac{{\beta }_2\left(r-\theta \right)k_m}{r}& 0\\ 0& \eta & -\delta & 0\\ 0& 0& 0& \epsilon \left(p_0-\frac{Dr-\theta \left(r-\theta \right)k_m}{ber}\right)\end{array}\right|\\ =\left(-r+\theta \right)\epsilon \left(p_0-\frac{Dr-\theta \left(r-\theta \right)k_m}{ber}\right)\left[-\delta \left(\frac{{\beta }_1\left(r-\theta \right)k_m}{r}-d\right)-\eta \frac{{\beta }_2\left(r-\theta \right)k_m}{r}\right]\end{array}$

$D_5=\left|\begin{array}{ccccc}-r+\theta & \frac{\left(\theta -r\right)\left(r+{\beta }_1{k}_m\right)}{r}& \frac{-{\beta }_2\left(r-\theta \right)k_m}{r}& \frac{{\left(r-\theta \right)}^2}{r}& 0\\ 0& \frac{{\beta }_1\left(r-\theta \right)k_m}{r}-d& \frac{{\beta }_2\left(r-\theta \right)k_m}{r}& 0& 0\\ 0& \eta & -\delta & 0& 0\\ 0& 0& 0& \epsilon \left(p_0-\frac{Dr-\theta \left(r-\theta \right)k_m}{ber}\right)& 0\\ -\phi \theta & 0& 0& 0& -\phi be\end{array}\right|=\left(-\phi be\right)D_4$

可知，当$p_0<f\left(\theta \right)$ ，$R_{01}<1$ 时，有$D_4>0$ ，$D_5<0$ ，所有的特征根都具有负实部，无病平衡点$E_1$ 是局部渐近稳定的，此时$f\left(\theta \right)=\frac{k_m}{ber}{\theta }^2-\frac{k_m}{be}\theta +\frac{D}{be}$ 。

定理3.2：价格阈值满足$p_0>f\left(\theta \right)$ ，$R_{02}<1$ 时，无病平衡点$E_2$ 是局部渐近稳定的，其中 $f\left(\theta \right)=\frac{k_m}{ber}{\theta }^2-\frac{k_m}{be}\theta +\frac{D}{be}$ 。

证明：定义$s\left(M\right)$ 为矩阵$M$ 的最大特征值，其中$M=F-V$ 。

根据[24]中定理2，易知无病平衡点$E_2$ 是局部稳定的所需满足的五个条件中，显然满足前四个条件；且可知：$R_{02}<1\iff s\left(M\right)<0$ ，$R_{02}>1\iff s\left(M\right)>0$ 。下证无病平衡点$E_2$ 的雅可比矩阵的所有特征值为具有负实部。对系统(2)的各项以$I$ ，$W$ ，$k$ ，$p$ ，$S$ 重新排列，得到无病平衡点$E_2$ 的雅可比矩阵为

${J|}_{E_2}=\left[\begin{array}{cc}M& 0\\ J_1& J_2\end{array}\right]$

其中

$J_2=\left[\begin{array}{ccc}0& \epsilon \frac{r\left(D-be{p}_0\right)}{\theta \left(r-\theta \right)}\left(1-\frac{r\left(D-be{p}_0\right)}{k_m\theta \left(r-\theta \right)}\right)& 0\\ 0& -\phi be& -\phi \theta \\ \frac{{\left(r-\theta \right)}^2}{r}& 0& \theta -r\end{array}\right]$

当$R_{02}<1$ 时，有$s\left(M\right)<0$ 。因此只需要证明$J_2$ 的所有特征值都具有负实部，$J_2$ 的特征方程为：

${\lambda }^3+m_1{\lambda }^2+m_2\lambda +m_3=0$ 。

其中

$m_1=\left(r-\theta \right)+\phi be$ ，$m_2=\left(r-\theta \right)\phi \theta$ ，

$m_3=\left(\theta -r\right)\phi be\epsilon \frac{r\left(D-be{p}_0\right)}{\theta \left(r-\theta \right)}\left(1-\frac{r\left(D-be{p}_0\right)}{k_m\theta \left(r-\theta \right)}\right)$ 。

可知当$p_0>f\left(\theta \right)$ ，$R_{02}<1$ 时，有$m_1>0$ ，$m_2>0$ ，$m_3>0$ ，其中$f\left(\theta \right)=\frac{k_m}{ber}{\theta }^2-\frac{k_m}{be}\theta +\frac{D}{be}$ 。

$H_1=m_1>0$ ，$H_2=\left[\begin{array}{cc}{m}_1& m_3\\ 1& m_2\end{array}\right]=m_1{m}_2-m_3>0$ ，$H_3=\left[\begin{array}{ccc}{m}_1& m_3& 0\\ 1& m_2& 0\\ 0& m_1& m_3\end{array}\right]=m_3{H}_2>0$ ，

根据Hurwitz准则，方程的所有根都具有负实部，无病平衡点$E_2$ 是局部渐近稳定的。

3.5. 正平衡点的局部稳定性

定理3.3：价格阈值满足$p_0<g\left(\theta \right)$ ，$R_{01}>1$ 时，正平衡点$E_3=\left(S_1^{*},I_1^{*},W_1^{*},k_m,p_1^{*}\right)$ 是局部渐近稳定的，其中$g\left(\theta \right)=\left(\frac{cd\delta k_m}{B}-\frac{d\delta }{beA}\right)\theta +\frac{D}{be}-\frac{cr{k}_{m}d\delta }{B}+\frac{cr{\left(d\delta \right)}^2}{AB}$ ，$A={\beta }_1\delta +{\beta }_2\eta$ ，$B=r\delta +{\beta }_1{k}_m\delta +{\beta }_2{k}_m\eta$ 。

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

$J\left(E_3\right)=\left[\begin{array}{ccccc}r-\theta -\frac{2r{S}_1^{*}-r{I}_1^{*}}{k_m}-{\beta }_1{I}_1^{*}-{\beta }_2{W}_1^{*}& -\frac{r{S}_1^{*}}{k_m}-{\beta }_1{S}_1^{*}& -{\beta }_2{S}_1^{*}& \frac{r{S}_1^{*}\left(S_1^{*}+I_1^{*}\right)}{{\left(k_m\right)}^2}& 0\\ {\beta }_1{I}_1^{*}& {\beta }_1{S}_1^{*}-d& {\beta }_2{S}_1^{*}& 0& 0\\ 0& \eta & -\delta & 0& 0\\ 0& 0& 0& -\epsilon \left(p_1^{*}-p_0-cd{I}_1^{*}\right)& 0\\ -\phi \theta & 0& 0& 0& -\phi be\end{array}\right]$ ，

可得特征方程为

$\left(\lambda +\phi be\right)\left(\lambda +\epsilon \left(p_1^{*}-p_0-cd{I}_1^{*}\right)\right)\left({\lambda }^3+a_1{\lambda }^2+a_2\lambda +a_3\right)=0$ ，

可得特征根${\lambda }_1=-\phi be<0$ ，${\lambda }_2=-\epsilon \left(p_1^{*}-p_0-cd{I}_1^{*}\right)$ 。

下面分析

${\lambda }^3+a_1{\lambda }^2+a_2\lambda +a_3=0$ ，

其中$a_1=g_1+\xi$ ，$a_2=g_1\xi +g_2$ ，$a_3=g_2\xi +g_3$ ，

$g_1=d-{\beta }_1{S}_1^{*}+\delta$ ，

$g_2={\beta }_1{I}_1^{*}\left(\frac{r{S}_1^{*}}{k_m}+{\beta }_1{S}_1^{*}\right)-\delta \left({\beta }_1{S}_1^{*}-d\right)-\eta {\beta }_2{S}_1^{*}$ ，

$g_3=-\delta {\beta }_1{I}_1^{*}\left(\frac{r{S}_1^{*}}{k_m}+{\beta }_1{S}_1^{*}\right)+\eta {\beta }_2{S}_1^{*}$ ，

$\xi =\frac{2r{S}_1^{*}-r{I}_1^{*}}{k_m}+{\beta }_1{I}_1^{*}+{\beta }_2{W}_1^{*}-\left(r-\theta \right)$ 。

同时对于特征根${\lambda }_2=-\epsilon \left(p_1^{*}-p_0-cd{I}_1^{*}\right)$ ，可知有$p_0<g\left(\theta \right)$ ，$R_{01}>1$ 时，有

${\lambda }_2=-\epsilon \left(p_1^{*}-p_0-cd{I}_1^{*}\right)<0$ ，

$\xi >0$ ，$g_1>0$ ，$g_2>0$ ，$g_3>0$ ，

故有$a_1>0$ ，$a_2>0$ ，$a_3>0$ 。

$H_1=a_1>0$ ，$H_2=\left[\begin{array}{cc}{a}_1& a_3\\ 1& a_2\end{array}\right]=a_1{a}_2-a_3>0$ ，

$H_3=\left[\begin{array}{ccc}{m}_1& m_3& 0\\ 1& a_2& 0\\ 0& a_1& a_3\end{array}\right]=a_3{H}_2>0$ ，

根据Hurwitz准则，方程的所有根都具有负实部。

因此$p_0<g\left(\theta \right)$ ，$R_{01}>1$ 时，正平衡点$E_3$ 是局部渐近稳定的，其中

$g\left(\theta \right)=\left(\frac{cd\delta k_m}{B}-\frac{d\delta }{beA}\right)\theta +\frac{D}{be}-\frac{cr{k}_{m}d\delta }{B}+\frac{cr{\left(d\delta \right)}^2}{AB}$ ，

且有$A={\beta }_1\delta +{\beta }_2\eta$ ，$B=r\delta +{\beta }_1{k}_m\delta +{\beta }_2{k}_m\eta$ 。

定理3.4：价格阈值满足$p_0>g\left(\theta \right)$ ，$R_{02}>1$ 时，正平衡点$E_4=\left(S_2^{*},I_2^{*},W_2^{*},k_2^{*},p_0\right)$ 是局部渐近稳定的，其中$g\left(\theta \right)=\left(\frac{cd\delta k_m}{B}-\frac{d\delta }{beA}\right)\theta +\frac{D}{be}-\frac{cr{k}_{m}d\delta }{B}+\frac{cr{\left(d\delta \right)}^2}{AB}$ ，$A={\beta }_1\delta +{\beta }_2\eta$ ，$B=r\delta +{\beta }_1{k}_m\delta +{\beta }_2{k}_m\eta$ 。

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

$J\left(E_4\right)=\left[\begin{array}{ccccc}a& b^{*}& c& d& 0\\ e^{*}& f& g& 0& 0\\ 0& h& i& 0& 0\\ 0& j& 0& k& n\\ m& 0& 0& 0& x\end{array}\right]$ ，

其中

$a=r-\theta -\frac{2r{S}_2^{*}-r{I}_2^{*}}{k_2^{*}}-{\beta }_1{I}_2^{*}-{\beta }_2{W}_2^{*}$ ，$b^{*}=-\frac{r{S}_2^{*}}{k_2^{*}}-{\beta }_1{S}_2^{*}$ ，$c=-{\beta }_2{S}_2^{*}$ ，$d=\frac{r{S}_2^{*}\left(S_2^{*}+I_2^{*}\right)}{{\left(k_2^{*}\right)}^2}$ ，

$e^{*}={\beta }_1{I}_2^{*}$ ，$f={\beta }_1{S}_2^{*}-d$ ，$g={\beta }_2{S}_2^{*}$ ，$h=\eta$ ，$g=-\delta$ ，

$j=-\left(\epsilon cd\right)k_2^{*}\left(1-\frac{k_2^{*}}{k_m}\right)$ ，$k=-\left(\epsilon cd{I}_2^{*}\right)\left(1-\frac{2k_2^{*}}{k_m}\right)$ ，$n=\epsilon k_2^{*}\left(1-\frac{k_2^{*}}{k_m}\right)$ ，$m=-\phi \theta$ ，$x=-\phi be$ 。

$J\left(E_4\right)$ 的特征方程为：

${\lambda }^5+b_1{\lambda }^4+b_2{\lambda }^3+b_3{\lambda }^2+b_4\lambda +b_5=0$ ，

$b_1=h_1+\psi$ ，$b_2=h_1\psi +h_2$ ，$b_3=h_2\psi +h_3$ ，$b_4=h_3\psi +h_4$ ，$b_5=h_4\psi +h_5$ ，

其中

$\psi =-x$ ，$h_1=-\left(a+k+f+i\right)$ ，$h_2=k\left(a+f+i\right)+\left(a+f\right)i+fa-eb-hg$ ，

$h_3=\left(hg-\left(a+f\right)i-fa+eb\right)k+\left(-fa+eb\right)i+hga-\left(ch+dj\right)e+mnf$ ，

$h_4=-mndx+\left(i\left(fa-eb\right)-h\left(ag-ec\right)\right)k+d\left(\left(ej+mn\right)i+mnf\right)$ ，

$h_5=mnd\left(-x^2+\left(f+i\right)x-\left(if-hg\right)\right)$ 。

可知当$p_0>g\left(\theta \right)$ ，$R_{02}>1$ 时，有$\psi >0$ ，$h_i>0\left(i=1,2,3,4,5\right)$ ，故有$b_i>0\left(i=1,2,3,4,5\right)$ 。

$H_1=b_1>0$ ，$H_2=\left[\begin{array}{cc}{b}_1& b_3\\ 1& b_2\end{array}\right]=b_1{b}_2-b_3>0$ ，$H_3=\left[\begin{array}{ccc}{b}_1& b_3& 0\\ 1& b_2& b_4\\ 0& b_1& b_3\end{array}\right]=b_3{H}_2>0$ ，

$H_4=\left[\begin{array}{cccc}{b}_1& b_3& b_5& 0\\ 1& b_2& b_4& 0\\ 0& b_1& b_3& b_5\\ 0& 1& b_2& b_2\end{array}\right]=b_4{H}_3>0$ ，$H_5=\left[\begin{array}{ccccc}{b}_1& b_3& b_5& 0& 0\\ 1& b_2& b_4& 0& 0\\ 0& b_1& b_3& b_5& 0\\ 0& 1& b_2& b_4& 0\\ 0& 0& b_1& b_3& b_5\end{array}\right]=b_5{H}_4>0$ 。

根据Hurwitz准则，方程的所有根都具有负实部。

因此$p_0>g\left(\theta \right)$ ，$R_{02}>1$ 时，正平衡点$E_4=\left(S_2^{*},I_2^{*},W_2^{*},k_2^{*},p_0\right)$ 是局部渐近稳定的，其中

$g\left(\theta \right)=\left(\frac{cd\delta k_m}{B}-\frac{d\delta }{beA}\right)\theta +\frac{D}{be}-\frac{cr{k}_{m}d\delta }{B}+\frac{cr{\left(d\delta \right)}^2}{AB}$ ，

且有$A={\beta }_1\delta +{\beta }_2\eta$ ，$B=r\delta +{\beta }_1{k}_m\delta +{\beta }_2{k}_m\eta$ 。

4. 最优控制策略

4.1. 建立控制模型

为降低布鲁氏菌病爆发带来的公共卫生风险与经济损失，应在疾病传播过程中采取相应控制措施。当动物疫病在牲畜种群中爆发时，为抑制其传染规模所采取的控制措施将导致大量成本支出，同时市场供需及价格也会受到冲击，从而引发额外的养殖成本与损失，对养殖户收益产生负面影响。基于此，本文从养殖户视角出发，以实现疾病传染规模最小化和控制成本最小化为目标，构建相应的目标函数。假设最优控制策略的终端时间$T$ 已知，本文中控制策略的实施时间设定为2023年至2033年。

本节将基于布鲁氏菌病的传播特点，对传染病模型添加控制措施，以求抑制布鲁氏菌病的扩散以及感染人数的激增。在系统(2)的基础上，增加了三种控制措施$u_i\left(t\right)\left(i=1,2,3\right)$ 以减少甚至根除这种疾病。

控制项$u_1\left(t\right)\left(0\le u_1\left(t\right)\le 1\right)$ 用于减少易感动物与受感染动物的接触，通过将动物进行分区隔离来实现管理。控制项$u_2\left(t\right)\left(0\le u_2\left(t\right)\le 1\right)$ 表示对感染牲畜进行扑杀。养殖户应及时对动物种群进行检测，及时扑杀检测结果为阳性的感染个体。控制项$u_3\left(t\right)\left(0\le u_3\left(t\right)\le 1\right)$ 表示为减少环境中的布鲁氏菌，通过定期清洁设备与动物分泌物等努力实现。

得到最优控制模型为：

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=rS\left(1-\frac{S+I}{k}\right)-\left(1-u_1\left(t\right)\right){\beta }_{1}SI-{\beta }_{2}SW-\theta S\\ \frac{\text{d}I}{\text{d}t}=\left(1-u_1\left(t\right)\right){\beta }_{1}SI+{\beta }_{2}SW-\left(1-u_2\left(t\right)\right)dI\\ \frac{\text{d}W}{\text{d}t}=\left(1-u_3\left(t\right)\right)\left(\eta I-\delta W\right)\\ \frac{\text{d}k}{\text{d}t}=\epsilon \left(p-p_0-c\left(1+u_2\left(t\right)\right)dI\right)k\left(1-\frac{k}{k_m}\right)\\ \frac{\text{d}p}{\text{d}t}=\phi \left(D-bep-\theta S\right)\end{array}$ (4)

其中$S>0$ ，$I>0$ ，$W>0$ ，$k>0$ ，$p>0$ ，并且$\left(u_1,u_2,u_3\right)\in U$ 。$U$ 为控制集，满足

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

目标函数的目标为最大限度的控制疾病的传播以及降低养殖户控制疾病的成本。目标函数为：

$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}$ (5)

存在最优控制$\left(u_1^{*},u_2^{*},u_3^{*}\right)$ ，使得$J\left(u_1^{*},u_2^{*},u_3^{*}\right)=\min J\left(u_1,u_2,u_3\right)$ ，$\left(u_1^{*},u_2^{*},u_3^{*}\right)\in U$ 。

4.2. 控制解

4.2.1. 最优控制解的存在性

对于最优控制问题，下证明最优控制解$\left(u_1^{*},u_2^{*},u_3^{*}\right)$ 的存在性。

定理4.1：使得目标函数最小化的控制组合解$\left(u_1^{*},u_2^{*},u_3^{*}\right)$ 存在，如果满足以下条件：

(1) 控制集和状态变量集合均是非空的；

(2) 控制集$U$ 是凸闭集；

(3) 目标函数的被积函数是关于控制集$U$ 的凸函数；

(4) 由控制系统等式右端所组成的向量函数$f\left(t,x,u\right)$ ，满足$f\left(t,x,u\right)=f_1\left(t,x\right)+f_2\left(t,x\right)u$ ，且$\left|f\left(t,x,u\right)\right|\le a\left(1+\left|x\right|+\left|u\right|\right)$ ，其中常数$a>0$ ；

(5) 存在常数$a_1,a_2,\theta >0$ ，使得目标函数的被积函数$L\left(t,x,u_1,u_2,u_2\right)$ 满足

$L\left(t,x,u_1,u_2,u_2\right)\ge a_2-a_1{\left({\left|u_1\right|}^2+{\left|u_2\right|}^2+{\left|u_3\right|}^2\right)}^{\frac{\theta }{2}}$

证明：

(1) 由于控制系统中的系数均为正数且系统的解有界，因此系统中的解存在，并且满足初始条件，因此可得：对于给定的初值$X_0$ ，控制$u$ 构成的集合$\left\{\left(X_0,u\right)\right\}$ 是非空的；控制集和状态变量集合均是非空；

(2) 控制集$U$ 是凸闭集；

(3) 目标被积函数$L\left(t,x,u\right)$ 为

$\begin{array}{l}L\left(t,x,u\right)=r_1\left(t,x,u\right)+r_2\left(t,x,u\right)\\ L\left(t,x,\left(1-k\right)v+kw\right)\ge \left(1-k\right)L\left(t,x,v\right)+k\end{array}$

其中

$r_1\left(t,x,u\right)=A_{1}I+A_{2}W$ ，$r_2\left(t,x,u\right)=\frac{B_1}{2}{u}_1^2+\frac{B_2}{2}{u}_2^2+\frac{B_3}{2}{u}_3^2$

取任意两点$v=\left(v_1,v_2,v_3\right)\in U$ ，$w=\left(w_1,w_2,w_3\right)\in U$ ，且$k\in \left[0,1\right]$ 。下证：

$L\left(t,x,\left(1-k\right)v+kw\right)\le \left(1-k\right)L\left(t,x,v\right)+kL\left(t,x,w\right)$ ，

即证

$r_2\left(t,x,\left(1-k\right)v+kw\right)\le \left(1-k\right)r_2\left(t,x,v\right)+k{r}_2\left(t,x,w\right)$ 。

根据定义有

$r_2\left(t,x,\left(1-k\right)v+kw\right)=\frac{1}{2}{\displaystyle \sum _{j=1}^3\left[C_j{\left(\left(1-k\right)v_j+k{w}_j\right)}^2\right]}$

$\left(1-k\right)r_2\left(t,x,v\right)+k{r}_2\left(t,x,w\right)=\frac{1}{2}\left(1-k\right){\displaystyle \sum _{j=1}^3{C}_j{v}_j^2-\frac{1}{2}k}{\displaystyle \sum _{j=1}^3{C}_j{w}_j^2}$

则

$r_2\left(t,x,\left(1-k\right)v+kw\right)-\left(\left(1-k\right)r_2\left(t,x,v\right)+k{r}_2\left(t,x,w\right)\right)=\frac{1}{2}\left(k^2-k\right){\displaystyle \sum _{j=1}^3{C}_j{\left(v_j-w_j\right)}^2}\le 0$ ，

因此目标函数的被积函数是关于控制集$U$ 的凸函数得以验证。

(4) 控制系统等式右端所组成的向量函数为

$g\left(t,x,u\right)=\left[\begin{array}{c}rS\left(1-\frac{S+I}{k}\right)-\left(1-u_1\left(t\right)\right){\beta }_{1}SI-{\beta }_{2}SW-\theta S\\ \left(1-u_1\left(t\right)\right){\beta }_{1}SI+{\beta }_{2}SW-\left(1+u_2\left(t\right)\right)dI\\ \left(1-u_3\left(t\right)\right)\left(\eta I-\delta W\right)\\ \epsilon \left(p-p_0-c\left(1+u_2\left(t\right)\right)dI\right)k\left(1-\frac{k}{k_m}\right)\\ \phi \left(D-bep-\theta S\right)\end{array}\right]$

显然函数$g\left(t,x,u\right)$ 可以表示为$g\left(t,x,u\right)=g_1\left(t,x\right)+g_2\left(t,x\right)u$ ，其中

$g_1\left(t,x\right)=\left[\begin{array}{c}rS\left(1-\frac{S+I}{k}\right)-{\beta }_{1}SI-{\beta }_{2}SW-\theta S\\ {\beta }_{1}SI+{\beta }_{2}SW-dI\\ \eta I-\delta W\\ \epsilon \left(p-p_0-cdI\right)k\left(1-\frac{k}{k_m}\right)\\ \phi \left(D-bep-\theta S\right)\end{array}\right]$ ，

$g_2\left(t,x\right)=\left[\begin{array}{ccc}{\beta }_{1}SI& 0& 0\\ -{\beta }_{1}SI& dI& 0\\ 0& 0& -\left(\eta I-\delta W\right)\\ 0& -\epsilon cdIk\left(1-\frac{k}{k_m}\right)& 0\\ 0& 0& 0\end{array}\right]$ 。

由于控制系统的有界性，有

$\begin{array}{c}\left|g\left(t,x,u\right)\right|=\left|g_1\left(t,x\right)+g_2\left(t,x\right)u\right|\\ \le \left|g_1\left(t,x\right)\right|+\left|g_2\left(t,x\right)\right|\left|u\right|\\ \le \sqrt{r^2+{\phi }^2{D}^2}+N_1\left|x\right|+N_2\left|x\right|\\ \le m\left(1+\left|x\right|+\left|u\right|\right)\end{array}$

其中$m=\max \left\{\frac{N_1}{\sqrt{r^2+{\phi }^2{D}^2}},\frac{N_1}{\sqrt{r^2+{\phi }^2{D}^2}}\right\}$ 。

(5) 由于$I\left(0\right)\ge 0$ ，$W\left(0\right)\ge 0$ ，有

$\begin{array}{c}L\left(t,x,u\right)=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\\ \ge \frac{B_1}{2}{u}_1^2+\frac{B_2}{2}{u}_2^2+\frac{B_3}{2}{u}_3^2\\ \ge \frac{B_{\min }}{2}\left(u_1^2+u_2^2+u_3^2\right)=\frac{B_{\min }}{2}{\left|u\right|}^2\end{array}$

其中$B_{\min }=\min \left\{B_1,B_2,B_3\right\}>0$ 。取$a_1=\frac{B_{\min }}{2}>0$ ，$n=3$ ，$a_2\ge 0$ ，于是有$L\left(t,x,u\right)\ge a_2-a_1{\left|u\right|}^2$ 。

综上所述，控制系统存在最优控制$u^{*}=\left(u_1^{*},u_2^{*},u_3^{*}\right)$ 使得$J\left(u_1^{*},u_2^{*},u_3^{*}\right)=\min J\left(u_1,u_2,u_3\right)$ ，$u^{*}\in U$ 。

4.2.2. 最优控制的理论解

给出控制解$\left(u_1^{*},u_2^{*},u_3^{*}\right)$ 满足的存在的必要条件，将控制系统下的目标函数通过庞特里亚金极小值原理转化为控制集$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(rS\left(1-\frac{S+I}{k}\right)-\left(1-u_1\left(t\right)\right){\beta }_{1}SI-{\beta }_{2}SW-\theta S\right)\\ +{\lambda }_2\left(\left(1-u_1\left(t\right)\right){\beta }_{1}SI+{\beta }_{2}SW-\left(1+u_2\left(t\right)\right)dI\right)\\ +{\lambda }_3\left(\left(1-u_3\left(t\right)\right)\left(\eta I-\delta W\right)\right)\\ +{\lambda }_4\left(\epsilon \left(p-p_0-cd\left(1+u_2\left(t\right)\right)I\right)k\left(1-\frac{k}{k_m}\right)\right)\\ +{\lambda }_5\left(\phi \left(D-bep-\theta S\right)\right)\end{array}$

定理4.2：控制系统(4)下使得目标函数(5)极小化的控制组合解$u^{*}=\left(u_1^{*},u_2^{*},u_3^{*}\right)$ 和对应的状态变量$\left(S^{\ast },I^{\ast },W^{\ast },k^{\ast },p^{\ast }\right)$ 存在，且最优控制变量的解为

$\{\begin{array}{l}{u}_1^{*}=\min \left\{1,\max \left\{0,\frac{\left({\lambda }_2-{\lambda }_1\right){\beta }_1{S}^{*}{I}^{*}}{B_1}\right\}\right\}\\ u_2^{*}=\min \left\{1,\max \left\{0,\frac{{\lambda }_2{I}^{*}d+{\lambda }_4\epsilon cd{k}^{*}I\left(1-\frac{k^{*}}{k_m}\right)}{B_2}\right\}\right\}\\ u_3^{*}=\min \left\{1,\max \left\{0,\frac{{\lambda }_3\left(\eta I-\delta W\right)}{B_3}\right\}\right\}\end{array}$

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

$\frac{\text{d}{\lambda }_1}{\text{d}t}=-\frac{\partial H}{\partial S}$ ，$\frac{\text{d}{\lambda }_2}{\text{d}t}=-\frac{\partial H}{\partial I}$ ，$\frac{\text{d}{\lambda }_3}{\text{d}t}=-\frac{\partial H}{\partial W}$ ，$\frac{\text{d}{\lambda }_4}{\text{d}t}=-\frac{\partial H}{\partial k}$ ，$\frac{\text{d}{\lambda }_5}{\text{d}t}=-\frac{\partial H}{\partial p}$ ，

且有${\lambda }_i\left(T\right)=0,i=1,2,3,4,5$ ；

则伴随变量满足

$\{\begin{array}{l}\frac{\text{d}{\lambda }_1}{\text{d}t}=-{\lambda }_1\left[r\left(1-\frac{2S+I}{k}\right)-\left(1-u_1\left(t\right)\right){\beta }_{1}I-{\beta }_{2}W-\theta \right]-{\lambda }_2\left[\left(1-u_2\left(t\right)\right){\beta }_{1}I+{\beta }_{2}W\right]+{\lambda }_5\phi \theta \\ \frac{\text{d}{\lambda }_2}{\text{d}t}=-A_1+{\lambda }_1\left(\frac{rS}{k}+\left(1-u_1\left(t\right)\right){\beta }_{1}S\right)+{\lambda }_{2}d\left(u_2\left(t\right)+1\right)-{\lambda }_3\left(1-u_3\left(t\right)\right)\eta +{\lambda }_4\epsilon cd\left(u_2\left(t\right)+1\right)k\left(1-\frac{k}{k_m}\right)\\ \frac{\text{d}{\lambda }_3}{\text{d}t}=-A_2+{\lambda }_1{\beta }_{2}S-{\lambda }_2{\beta }_{2}S+{\lambda }_3\delta \left(1-u_3\left(t\right)\right)\\ \frac{\text{d}{\lambda }_4}{\text{d}t}=\frac{-{\lambda }_{1}rS\left(S+I\right)}{k^2}-{\lambda }_4\epsilon \left(p-p_0-cd\left(u_2\left(t\right)+1\right)I\right)\left(1-\frac{2k}{k_m}\right)\\ \frac{\text{d}{\lambda }_5}{\text{d}t}=\phi be{\lambda }_5-{\lambda }_4\epsilon k\left(1-\frac{k}{k_m}\right)\end{array}$

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

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

$\{\begin{array}{l}{u}_1^{*}=\min \left\{1,\max \left\{0,\frac{\left({\lambda }_2-{\lambda }_1\right){\beta }_1{S}^{*}{I}^{*}}{B_1}\right\}\right\}\\ u_2^{*}=\min \left\{1,\max \left\{0,\frac{{\lambda }_2{I}^{*}d+{\lambda }_4\epsilon cd{k}^{*}I\left(1-\frac{k^{*}}{k_m}\right)}{B_2}\right\}\right\}\\ u_3^{*}=\min \left\{1,\max \left\{0,\frac{{\lambda }_3\left(\eta I-\delta W\right)}{B_3}\right\}\right\}\end{array}$

5. 数值模拟

为验证模型的合理性，本文选取新疆维吾尔自治区2011~2023年的活羊出栏价格数据以及人间布鲁氏菌病的累计患病数据进行参数估计，进而预测布鲁氏菌病的传播形式，寻求最优控制策略。

5.1. 参数估计

首先对无传染的价格养殖模型进行参数估计，通过查找[26] [27]，得到2011~2023年新疆地区的各项详细数据见表1，由于羊群的平均规模为19,260,000头，2011年羊的活重价格为每公斤为22.2202元，所以$N\left(0\right)=19260000$ ，$p\left(0\right)=21.2202$ 。假设$X\left(t\right)$ 为羊的实际出栏价格，$x\left(t\right)$ 为羊的理论出栏价格，并

且$\frac{\text{d}p}{\text{d}t}=\phi \left(D-bep-\theta S\right)$ ，采用最小二乘法进行参数估计，构建目标函数为：

$\min {\displaystyle \sum _t{\left[X\left(t\right)-x\left(t\right)\right]}^2}$

得到数据拟合结果以及参数估计值$b=2133$ ，$\phi =8.68\times {10}^{-7}$ ，$\theta =0.177$ 。数据拟合结果见图1(a)，图1(a)中的散点表示为实际价格数据值，区间范围为95%的置信区间。

Table 1. Live sheep slaughter prices and the number of human brucellosis cases in Xinjiang from 2011 to 2023

表1. 2011~2023年新疆地区的活羊出栏价格以及人间布鲁氏菌病病例数

年份	羊出栏的价格(元/kg)	人间布鲁氏菌病新增病例数	人间布鲁氏菌病累计病例数
2011	21.2202	1361	1361
2012	25.2128	2251	3612
2013	26.3232	3930	7542
2014	26.842	7493	15,035
2015	21.538	8997	24,032
2016	22.1288	7858	31,890
2017	23.0012	5447	37,337
2018	25.0772	3543	40,880
2019	27.4496	3627	44,507
2020	31.8692	2570	47,077
2021	31.6362	4152	51,229
2022	31.1984	5665	56,894
2023	27.6244	8277	65,171

(a) (b)

Figure 1. Fitting results of live sheep prices and human brucellosis data in Xinjiang

图1. 新疆地区活羊价格与人间布鲁氏菌病数据拟合结果

在此基础上，继续选用人间布鲁氏菌病累计病例采用最小二乘法进行参数估计，获取${\beta }_1$ 、${\beta }_2$ 、${\beta }_h$ 、${\beta }_w$ 的数值。假设初始值为$S\left(0\right)=1927\times {10}^4$ ，$I\left(0\right)=171534$ ，$W\left(0\right)=5\times {10}^5$ ，$k\left(0\right)=7\times {10}^7$ ，$p\left(0\right)=21.2202$ ，$S_h\left(0\right)=1361$ ，$I_h\left(0\right)=1361$ 。假设$Y\left(t\right)$ 为人间布鲁氏菌病的实际累计案例，$y\left(t\right)$ 为人间布鲁氏菌病的理

论累计案例，且满足$\frac{\text{d}y}{\text{d}t}={\beta }_h{S}_{h}I+{\beta }_w{S}_{w}W$ ，2011~2023年的实际患病数据见表1，其中人间布鲁氏菌患病

数据来自国家疾病通报监测系统(NNDSS) [28]，构建目标函数为：

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

得到数据拟合效果与参数值${\beta }_1=2.26\times {10}^{-10}$ ，${\beta }_2=3.05\times {10}^{-10}$ ，${\beta }_h=2.62\times {10}^{-10}$ ，${\beta }_w=3.23\times {10}^{-10}$ 。数据拟合结果见图1(b)，其中的散点表示为实际患病数据值，着色区间范围为95%的置信区间。

模型(3)中的所有参数详见表2。

Table 2. Parameter values and sources

表2. 参数值和来源

参数	数值	来源
$r$	0.449	参考
$\theta$	0.177	模拟
$\epsilon$	3.563 × 10^−^3	假设
$k_m$	6.5 × 10^7	假设
$D$	7.8 × 10^6	假设
$p_0$	18.8132	参考
$b$	2133	模拟
$e$	54.95	参考
$\phi$	8.68 × 10^−^7	模拟
$\eta$	0.000005	假设
$\delta$	0.0003	假设
$d$	0.00028	假设
$c$	1.45 × 10^−^3	假设
$A_h$	333,901	参考
$d_h$	4.56 × 10^−^7	参考
$m$	0.5	参考
${\beta }_1$	2.26 × 10^−^10	模拟
${\beta }_2$	3.05 × 10^−^10	模拟
${\beta }_h$	2.62 × 10^−^10	模拟
${\beta }_w$	3.23 × 10^−^10	模拟

5.2. 基本再生数分析

为深入探讨模型参数对基本再生数的影响，本节选取关键参数进行分析。保持模型中其他参数不变，预定参数区间，系统评估各参数对基本再生数的作用。通过理论分析求解得出模型具有两个不同的正平衡点，对于基本再生数$R_{01}$ ，选取参数$r$ 和$k_m$ 进行参数分析，此时$r$ 的变化阈值为[0.4, 0.9]，$k_m$ 的变化阈值为$\left[5.5\times {10}^7,7.5\times {10}^7\right]$ ；对于基本再生数$R_{02}$ ，选取参数$p_0$ 和$\theta$ 进行参数分析，此时$p_0$ 的变化阈值为[13, 22]，$\theta$ 的变化阈值为[0.4, 0.9]。具体的参数分析结果详见图2。

根据图2中的(a)，(b)，对于固定的$r$ ，当$k_m$ 越大时，基本再生数$R_{01}$ 的数值越大；对于固定的$k_m$ ，当$r$ 越大时，基本再生数$R_{01}$ 的数值越大，传染病的传播能力越强。此时对于种群规模大以及养殖密度高的种群来说，传染病蔓延的速度更快；并且种群内部的快速繁殖也增加了疾病传播风险。图2中的(c)，(d)表明，对于固定的$\theta$ ，基本再生数$R_{02}$ 的大小不会随着$p_0$ 的增大或减小而改变；而对于固定的$p_0$ ，基本再生数$R_{02}$ 将随着$\theta$ 的增大而减小。在现实养殖情境中，高出栏率往往带来牲畜养殖周期缩短，降低了种群内部接触频率，从而有利于控制种群内部疾病传播。

(a) (b)

(c) (d)

Figure 2. Basic reproduction number analysis results: (a) $R_{01}$ variation diagram regarding $r$ and $k_m$ , (b) $R_{01}$ stereoscopic results regarding $r$ and $k_m$ , (c) $R_{02}$ variation diagram regarding $p_0$ and $\theta$ , (d) $R_{02}$ stereoscopic results regarding $p_0$ and $\theta$

图2. 基本再生数分析结果：(a) $R_{01}$ 关于$r$ 和$k_m$ 的变化情况平面图，(b) $R_{01}$ 关于$r$ 和$k_m$ 的变化情况立体图，(c) $R_{02}$ 关于$p_0$ 和$\theta$ 的变化情况平面图，(d) $R_{02}$ 关于$p_0$ 和$\theta$ 的变化情况立体图

6. 最优控制策略分析

本节在三种控制措施的基础上进行自由组合，共形成七种控制策略，并据此探究布鲁氏菌病的最优防控方案。假设控制策略的实施期为2023年至2033年。具体模拟情景如下：(1) 使用单项控制$u_1$ ，对感染动物进行隔离；(2) 使用单项控制$u_2$ ，扑杀感染动物；(3) 使用单项控制$u_3$ ，清理环境中的布鲁氏菌病；(4) 使用双项控制$u_1$ 和$u_2$ ；(5) 使用双项控制$u_1$ 和$u_3$ ；(6) 使用双项控制$u_2$ 和$u_3$ ；(7) 三项控制同时使用。

6.1. 单项控制措施

(a) (b)

(c)

Figure 3. Numerical results of single control measures. (a) graph of farmer costs, (b) graph of changes in the number of infected sheep, and (c) graph of changes in control terms

图3. 单项控制措施的数值结果。(a) 控制成本图，(b) 感染羊的数量变化图，(c) 控制项变化图

图3展示了在实施单项控制措施时的数值结果。由图3(a)可见，各控制措施均使传染病防控成本相较无控制情景显著下降，其中策略3的防控成本最低。图3(b)则反映了不同情境下受感染羊群数量的变化，均呈现不同程度的减少，而策略2对抑制疾病传播的效果最为显著。图3(c)综合了不同控制措施下的具体实施情况。综合分析几项单一控制措施，发现在单项控制策略中，策略2可使感染规模达到最小，而策略3则可实现最低防控成本。

6.2. 双项控制措施

图4展示了在应用双项控制措施时的数值结果。由图4(a)可见，任意两种控制措施的联合应用均可使传染病防控成本相较于无控制情形明显降低，其中策略6产生的成本最低。图4(b)则呈现了各控制策略下受感染羊群数量的变化，结果表明实施策略6时，传染范围达到最小。图4(c)~(e)分别刻画了在分别实施控制策略4、5和6时，各项控制措施随时间的动态变化。综合上述结果，不论从防控效果还是成本效益的角度考虑，两项控制措施的联合应用中，策略6都是最佳选择。

6.3. 三项控制措施

图5展示了同时施加三种控制措施的数值结果。图5(a)显示，在三种控制措施联合实施时，防控成本相较于无控制情形及其他控制情形均显著降低；图5(b)则反映出，在此策略下受感染羊群数量出现了明显减少。图5(c)则刻画了三种控制措施随时间的动态变化。综合来看，无论从疾病防控效果还是防控成本的角度考虑，同时施加三种控制措施的策略(策略7)均表现出优良效果。

(a) (b)

(c) (d)

(e)

Figure 4. Numerical results of dual control measures. (a) graph of farmer costs, (b) graph of changes in the number of infected sheep, and (c) plot of changes in control terms for strategy 4, (d) plot of changes in control terms for strategy 5, and (e) plot of changes in control terms for strategy 6

图4. 双项控制措施的数值结果。(a) 控制成本图，(b) 感染羊的数量变化图，(c) 策略4的控制项变化图，(d) 策略5的控制项变化图，(e) 策略6的控制项变化图

(a) (b)

(c)

Figure5. Numerical results of three control measures. (a) graph of farmer losses, (b) graph of changes in the number of infected sheep, (c) graph of changes in control terms for strategy 7

图5. 三项控制措施的数值结果。(a) 控制成本图，(b) 感染羊的数量变化图，(c) 策略7的控制项变化图

综上所述，若只考虑施加一种控制措施，则选择策略2将可以达到疾病最优防控效果，选择策略3将可以达到防控成本最优。若同时施加两种控制措施，则选择策略六(同时施加控制项$u_2$ 和$u_3$ )可以同时达到疾病防控和成本防控最优效果。但相较于其余的控制策略效果，同时施加三种控制措施取得的防控效果最好，同时控制成本最低。

7. 结论

近年来，新疆布鲁氏菌病发病率和发病数均呈现较高态势，公共卫生领域以及畜牧业发展面临着严峻挑战。为有效防控布病，农业农村部印发了《畜间布鲁氏菌病防控五年行动方案(2022~2026年)》，但自2021年以来新疆牧区动物发病率和新发人间病例依旧不断上升，危害畜牧业的发展[5] [15]。因此，加强新疆地区的布病防控至关重要。

本文基于市场动态价格因素与布病传播特点，构建了一个受价格因素影响的传染病模型。首先对模型进行动力学分析，求解了该模型的平衡点及基本再生数，并分别验证其稳定性。为防控布病传播，对模型施加三个传染病控制措施，并从养殖户视角出发，以疾病传播范围最小与防控成本最低为目的建立目标函数，利用Pontryagin极值原理推导出控制项的理论最优解。接着选取新疆维吾尔自治区2011至2023年的活羊价格数据及人间布鲁氏菌病例数据进行参数估计，通过分析基本再生数随特定参数变动的情况，得出：养殖规模较大的种群中疾病传播效力更强，同时在出栏率较高的种群中疾病传播速率较低。最后对控制策略进行数值模拟，得出结论：(1) 在单项控制情形下，策略2实现了最佳疾病防控效果，而策略3对应的防控成本最低；(2) 在双重控制情形下，策略6能同时达到疾病防控与成本控制的最优效果；(3) 当三项控制措施同时施行时，在缩小传染规模和降低防控成本方面均优于其他控制策略。

不同于以往传染病控制问题，本文在传染病模型中创新性加入动态市场价格因素，综合考虑养殖户自身养殖行为以及传染病防控措施等因素，从多方面入手研究最优养殖策略以及传染病最优控制策略，对养殖户的实际养殖与疾病防控具有指导意义。在下一步工作中，我们将尝试制定动态变量刻画市场需求，以求更贴近现实场景。

基金项目

本文受到国家自然科学基金(批准号：12371493，11801398，12101443)，山西省基础研究计划自然科学研究面上项目(批准号：202303021221024，202103021224095)和山西省回国留学人员科研教研资助项目(2022-73)的资助。

NOTES

^*通讯作者。