年龄-环境耦合的两阶段结核病传播动力学模型的稳定性与预设目标控制

The Stability and Preset Target Control of a Two-Stage Tuberculosis Transmission Dynamical Model Coupled With Age and Environment

摘要

构建了一类兼顾易感人群年龄结构与环境中结核分枝杆菌分布的动力学模型，系统分析了模型的稳定性，并设计了预设目标控制策略。首先，推导得到基本再生数{\mathrm{ℛ}}_{\mathrm{0}}的显式表达式， 并基于李雅普诺夫（Lyapunov）稳定性理论，证明了无病平衡点与地方病平衡点的全局渐近稳定性。接着，引入疫苗接种与直接督导下短程化疗（DOTS）策略作为控制变量，提出预设目标控制问题并完成求解。最后，通过数值模拟发现，同时实施疫苗接种和DOTS策略可将活动性结核病患者的人数控制并稳定在预设防控目标之下，感染避免率达到\mathrm{45.01}\mathrm{\%}， 显著优于单一控制措施， 为结核病防控策略的优化提供了理论依据。

Abstract

By incorporating both the age structure of the susceptible population and Mycobacterium tuberculosis in environment， we construct a two-stage dynamical model of tuberculosis transmission and systematically analyze its stability. A preset target control strategy is further developed. First， we derive an explicit expression for the basic reproduction number and prove the globally asymptotic stability of both the disease-free and endemic equilibria using Lyapunov theory. Next， we formulate and solve a preset-target control problem that treats vaccination and the DOTS （directly observed treatment， short-course） strategy as control variables. Numerical simulations show that the combined implementation of vaccination and DOTS can maintain the number of active TB cases below the preset target and yields a significantly higher infection-avoidance rate than either measure used alone， providing a theoretical basis for optimizing tuberculosis prevention and control strategies.

Graphical abstract

关键词

结核病 / 动力学模型 / 基本再生数 / 稳定性 / 预设目标控制

Key words

tuberculosis / dynamical model / basic reproduction number / stability / preset target control

结核病 （TB） 是由结核分枝杆菌 （MTB） 引起的高传染性慢性呼吸道疾病， 主要侵犯肺部。该病既可通过人际密切接触直接传播， 也能随气流在环境中扩散， 进而导致广泛流行^［1］。 结核病的治疗需采用多种药物联合方案，治疗周期长、成本高，这使得结核病的临床治疗面临巨大挑战^［2］。 近年来， 全球每年有120 万~160万人死于结核病， 这与 2030 年终止结核病流行的目标仍相去甚远^［3］。 因此，相关部门正加大资金投入，多策并举遏制结核病传播，降低发病率与死亡率。这一目标的实现，亟须研究者厘清其传播机制，在有限资源下寻求最优干预措施。

基于结核病的传播机理，构建数学模型对结核病流行趋势进行预测，并探寻有效控制策略，是生物数学及公共卫生领域的研究热点。现有研究多聚焦人际直接传播^［4-6］， 如Zhang等^［5］提出并分析了一个考虑媒体影响的结核病预防性治疗六维模型，采用中国 4 个地区 2009—2019 年的新报告结核病病例数据对模型进行拟合与参数估计，研究发现，适当提高感染者的及时治疗比例以及潜伏结核感染人群的预防性治疗寻求比例，可实现结核病消除的目标；但媒体影响仅能在有限程度上减少活动性感染者数量，无法改变结核病的流行程度。

以上研究忽视了一个关键事实： 活动性结核病患者排出的结核分枝杆菌可在空气中存活超9 h， 附着尘埃后感染力可持续 8~10 d^［7-8］， 这意味着环境间接传播同样会助推疾病的传播。目前，关于间接传播对结核病流行影响的研究仍较为匮乏^［9-13］。

张正斌等^［9］检索了国内外近 5 a关于结核病季节分布特征的文献，发现若干季节性因素可促进结核分枝杆菌在人群中的传播，例如冬季人们户外活动减少、室内空气污染较重（病原菌浓度较高），以及部分国家和地区冬季雾霾天气频发等。 Cai 等^［10］运用数学模型探讨了环境中结核分枝杆菌载量对结核病传播的影响，并基于江苏省的实际数据估计模型参数，发现其基本再生数{\mathrm{ℛ}}_{\mathrm{0}}>\mathrm{1}，江苏省的结核病呈地方病流行态势；为有效控制该地区结核病，提出降低病菌排放率、提高患者康复率与环境病原体清除率的综合控制方案。Li 等^［11］建立了一个考虑环境传播的复发性结核病SVEIRB模型， 利用Lyapunov稳定性理论和LaSalle不变集原理，证明了模型的无病平衡点和地方病平衡点具有全局渐近稳定性；数值模拟结果表明，环境传播会加剧结核病的扩散， 同时，提出疫苗接种联合降低环境病原体负荷是控制疾病传播的有效策略。

此外， 考虑到年龄差异会直接导致人群免疫力水平的异质性， 年龄也是影响结核病传播动力学特征的关键因素^［14-18］。 杨应周等^［14］指出，结核病传播受年龄和环境协同影响显著：学生免疫系统尚未完善，加之学校空间拥挤、通风不良，易发生聚集性疫情，应针对性落实环境优化等防控措施。 张立兴等^［15］对比1980—2002年北京市肺结核发病率、感染率变化及DOTS实施效果，发现年龄小于30岁组的发病率下降与人群感染率降低密切相关；DOTS 通过控制传染源显著降低该年龄组感染与发病风险，但对已感染者后续发病影响有限。 Jing 等^［16］构建含环境因素的年龄结构肺结核传播模型，结合江苏省监测数据，采用马尔可夫链蒙特卡洛（MCMC）方法估计平均感染周期为 44.3 d；模拟曲线与实际新发病例数高度吻合，为区域防控提供了量化依据，并指出诊断+疫苗接种的联合干预效果显著。考虑到从潜伏感染个体发展为活动性感染个体存在一定的时间延迟，Gao 等^［17］建立了具有年龄结构和复发特性的结核病传播动力学模型，得到了关于模型解稳定性的相关结论，并将该模型应用于描述中国结核病的传播态势，结果显示，模型预测的总人口数和年度新报告结核病病例数均与统计数据高度契合。 Xue 等^［18］考虑到中国不同年龄段结核病患病率差异显著，提出了一个具有年龄结构和季节性传播率的非自治微分方程模型， 证明了当时，唯一的无病周期解P_{\mathrm{0}}全局渐近稳定；反之，疾病将均匀持续存在，且至少存在一个正的周期解。研究结果表明，为65岁以上及20~24岁易感人群接种疫苗，在降低结核病患病率方面效果显著。

在对模型稳定性等相关动力学性质开展系统分析的基础上， 已有部分学者将最优控制策略引入结核病传播模型的研究框架中^{［13，19-20］}。 Khoda 等^［19］建立了一个结核病传播模型， 并应用Pontryagin极大值原理（Pontryagin’s maximum principle）确定了媒体和教育提升公众意识、提升检测率和改善治疗服务3种干预措施的最佳控制水平。 Obsu 等^［20］将最优控制理论应用于由非线性常微分方程系统描述的结核病模型，探讨了治疗失败对结核病流行的影响， 通过Pontryagin极大值原理推导了最优路径的特征， 设计了不同的模拟案例，并与分析结果进行对比。结果表明，通过媒体报道提高公众意识且在治疗期间进行持续监督的综合干预效果，有助于降低治疗失败率，进而减少社区内结核病的流行。

综上所述， 当前结核病领域的相关探索多聚焦于环境、年龄及复发等关键影响因素^［9-18］， 并在此基础上进一步拓展至最优控制问题的研究^{［13，19-20］}。 然而， 当控制目标接近于零时， 所需投入的控制成本往往呈非线性增长态势， 且实际防控资源存在客观约束， 导致相关控制措施难以有效落地实施。 因此， 为更精准解析结核病的传播机制、科学制定防控策略， 文中构建年龄-环境耦合的两阶段结核病传播动力学模型， 并将疫苗接种与直接督导下短程化疗（DOTS）两项核心控制措施纳入模型框架， 研究模型的预设目标控制问题， 旨在以有限的控制成本实现更优的防控效果。

1 模型的建立和基本再生数的计算

1.1 模型的建立

结核病的传播是传染源、环境条件和宿主免疫等多重因素相互交织的动态过程。 为了建模需要，强调如下事实：

（ⅰ） 结核病的核心传播途径为呼吸道飞沫传播。即活动性结核病患者咳嗽、打喷嚏时产生的含菌气溶胶， 经易感人群吸入后引发感染， 而接触传播仅在特殊场景 （如含菌污染物直接接触破损黏膜） 下偶发， 并非主要传播形式^［7-8］。 因此，间接传染率大于直接传染率。基于此，同时纳入直接传染与间接传染两种传播方式。

（ⅱ） 结核病的感染率与年龄存在一定相关性。在中国，仅新生儿接种卡介苗，且接种后该疫苗的保护效果会逐年下降 （尤其在接种15 a后） ^［21］。据此， 可将易感人群分为 “青少年及幼儿” 与 “成人” 两类， 以体现不同免疫屏障水平对疾病传播的影响。

（ⅲ） 处于潜伏期的人群虽携带病原体，但因菌量极低而不具备传播能力。此外，感染者即便被成功治愈，肺部遗留的不可逆病灶仍使其终生面临复发风险。因此，假设康复者具备一定防护意识与免疫记忆，若再次发病，仍需经历新的潜伏期，方可重新进入活动期^［22］。

令S_{\mathrm{1}}(t)\mathrm{、}{S}_{\mathrm{2}}(t)\mathrm{、}E(t)\mathrm{、}I(t)\mathrm{、}R(t)和W(t)表示t时刻易感青少年及幼儿、易感成年人、潜伏者、活动性结核病感染者、康复者的数量， 以及环境中的结核分枝杆菌载量， 且N(t)=S_{\mathrm{1}}(t)+S_{\mathrm{2}}(t)+E(t)+I(t)+R(t)。为表达简便， 将S_{\mathrm{1}}(t)\mathrm{、}{S}_{\mathrm{2}}(t)\mathrm{、}E(t)\mathrm{、}I(t)\mathrm{、}R(t)\mathrm{、}W(t)\mathrm{、}N(t)记为S_{\mathrm{1}}\mathrm{、}{S}_{\mathrm{2}}\mathrm{、}E\mathrm{、}I\mathrm{、}R\mathrm{、}W\mathrm{、}N。年龄-环境耦合的两阶段结核病传播机制如图1所示，其动力学模型为

\left\{\begin{array}{l}\frac{\mathrm{d}{S}_{\mathrm{1}}(t)}{\mathrm{d}t}=A-S_{\mathrm{1}}(t)({\beta }_{\mathrm{1}}I(t)+{\beta }_{\mathrm{2}}W(t))-\gamma S_{\mathrm{1}}(t)-\delta S_{\mathrm{1}}(t),\\ \frac{\mathrm{d}{S}_{\mathrm{2}}(t)}{\mathrm{d}t}=\gamma S_{\mathrm{1}}(t)-S_{\mathrm{2}}(t)({\beta }_{\mathrm{3}}I(t)+{\beta }_{\mathrm{4}}W(t))-\delta S_{\mathrm{2}}(t),\\ \frac{\mathrm{d}E(t)}{\mathrm{d}t}=c{S}_{\mathrm{1}}(t)({\beta }_{\mathrm{1}}I(t)+{\beta }_{\mathrm{2}}W(t))+c{S}_{\mathrm{2}}(t)({\beta }_{\mathrm{3}}I(t)+{\beta }_{\mathrm{4}}W(t))+kR(t)-\phi E(t)-\delta E(t)\\ \frac{\mathrm{d}I(t)}{\mathrm{d}t}=(\mathrm{1}-c)S_{\mathrm{1}}(t)({\beta }_{\mathrm{1}}I(t)+{\beta }_{\mathrm{2}}W(t))+(\mathrm{1}-c)S_{\mathrm{2}}(t)({\beta }_{\mathrm{3}}I(t)+{\beta }_{\mathrm{4}}W(t))+\phi E(t)-hI(t)-(\delta +\mu )I(t),\\ \frac{\mathrm{d}R(t)}{\mathrm{d}t}=hI(t)-kR(t)-\delta R(t),\\ \frac{\mathrm{d}W(t)}{\mathrm{d}t}=\sigma I(t)-dW(t)\mathrm{。}\end{array}\right.（1）

根据实际意义， 所有参数均非负， 其生物学意义如表1所示。

由模型（1）可得

\frac{\mathrm{d}N}{\mathrm{d}t}=\frac{\mathrm{d}(S_{\mathrm{1}}+S_{\mathrm{2}}+E+I+R)}{\mathrm{d}t}=A-\delta (S_{\mathrm{1}}+S_{\mathrm{2}}+E+I+R)-\mu I\le A-\delta N，

从而有

\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}\mathrm{s}\mathrm{u}\mathrm{p}}N(t)=\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}\mathrm{s}\mathrm{u}\mathrm{p}}(S_{\mathrm{1}}+S_{\mathrm{2}}+E+I+R)\le \frac{A}{\delta }\mathrm{。}

类似地， 可得

\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}\mathrm{s}\mathrm{u}\mathrm{p}}W(t)=\frac{\sigma A}{\delta d}\mathrm{。}

所以， 集合

\Omega =\left\{\left(S_{\mathrm{1}},S_{\mathrm{2}},E,I,R,W\right)\in {\mathbf{R}}_{+}^{\mathrm{6}}:S_{\mathrm{1}}+S_{\mathrm{2}}+E+I+R\le \frac{A}{\delta },W\le \frac{\sigma A}{\delta d}\right\}

为模型（1）的正不变集。

1.2 模型的基本再生数

模型（1）的无病平衡点为P^{\mathrm{0}}=(S_{\mathrm{1}}^{\mathrm{0}},S_{\mathrm{2}}^{\mathrm{0}},\mathrm{0,0},\mathrm{0,0})， 其中

S_{\mathrm{1}}^{\mathrm{0}}=\frac{A}{\gamma +\delta },S_{\mathrm{2}}^{\mathrm{0}}=\frac{\gamma A}{\delta (\gamma +\delta )}\mathrm{。}

由模型（1）可得

\mathit{V}=\left[\begin{array}{cccc}\phi +\delta & \mathrm{0}& -k& \mathrm{0}\\ -\phi & h+\delta +\mu & \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -h& k+\delta & \mathrm{0}\\ \mathrm{0}& -\sigma & \mathrm{0}& d\end{array}\right]，\mathit{F}=\left[\begin{array}{cccc}\mathrm{0}& c{S}_{\mathrm{1}}^{\mathrm{0}}{\beta }_{\mathrm{1}}+c{S}_{\mathrm{2}}^{\mathrm{0}}{\beta }_{\mathrm{3}}& \mathrm{0}& c{S}_{\mathrm{1}}^{\mathrm{0}}{\beta }_{\mathrm{2}}+c{S}_{\mathrm{2}}^{\mathrm{0}}{\beta }_{\mathrm{4}}\\ \mathrm{0}& a_{\mathrm{1}}& \mathrm{0}& a_{\mathrm{2}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]\mathrm{。}

其中：a_{\mathrm{1}}=(\mathrm{1}-c)S_{\mathrm{1}}^{\mathrm{0}}{\beta }_{\mathrm{1}}+(\mathrm{1}-c)S_{\mathrm{2}}^{\mathrm{0}}{\beta }_{\mathrm{3}}， a_{\mathrm{2}}=(\mathrm{1}-c)S_{\mathrm{1}}^{\mathrm{0}}{\beta }_{\mathrm{2}}+(\mathrm{1}-c)S_{\mathrm{2}}^{\mathrm{0}}{\beta }_{\mathrm{4}}。由文献［13］可得模型（1）的基本再生数为

{\mathrm{ℛ}}_{\mathrm{0}}=\rho \left(F{V}^{-\mathrm{1}}\right)=\frac{A(k+\delta )(\delta {\beta }_{\mathrm{1}}+\gamma {\beta }_{\mathrm{3}})(\phi +\delta -c\delta )}{\delta (\gamma +\delta )[(k+\delta )(\phi +\delta )(h+\delta +\mu )-kh\phi ]}\left(\mathrm{1}+\frac{\sigma }{d}\right)\mathrm{。}

定理1 当{\mathrm{ℛ}}_{\mathrm{0}}>\mathrm{1}时， 模型（1）存在唯一的地方病平衡点P^{\mathrm{*}}=\left(S_{\mathrm{1}}^{\mathrm{*}},S_{\mathrm{2}}^{\mathrm{*}},E^{\mathrm{*}},I^{\mathrm{*}},R^{\mathrm{*}},W^{\mathrm{*}}\right), 且这些分量满足：

\left\{\begin{array}{l}A=S_{\mathrm{1}}^{\mathrm{*}}\left({\beta }_{\mathrm{1}}{I}^{\mathrm{*}}+{\beta }_{\mathrm{2}}{W}^{\mathrm{*}}\right)+\gamma S_{\mathrm{1}}^{\mathrm{*}}+\delta S_{\mathrm{1}}^{\mathrm{*}},\\ \gamma S_{\mathrm{1}}^{\mathrm{*}}=S_{\mathrm{2}}^{\mathrm{*}}\left({\beta }_{\mathrm{3}}{I}^{\mathrm{*}}+{\beta }_{\mathrm{4}}{W}^{\mathrm{*}}\right)+\delta S_{\mathrm{2}}^{\mathrm{*}},\\ c{S}_{\mathrm{1}}^{\mathrm{*}}\left({\beta }_{\mathrm{1}}{I}^{\mathrm{*}}+{\beta }_{\mathrm{2}}{W}^{\mathrm{*}}\right)+c{S}_{\mathrm{2}}^{\mathrm{*}}\left({\beta }_{\mathrm{3}}{I}^{\mathrm{*}}+{\beta }_{\mathrm{4}}{W}^{\mathrm{*}}\right)+k{R}_{\mathrm{1}}^{\mathrm{*}}=\phi E^{\mathrm{*}}+\delta E^{\mathrm{*}},\\ h{I}^{\mathrm{*}}+\left(\delta +\mu \right)I^{\mathrm{*}}=\left[S_{\mathrm{1}}^{\mathrm{*}}\left({\beta }_{\mathrm{1}}{I}^{\mathrm{*}}+{\beta }_{\mathrm{2}}{W}^{\mathrm{*}}\right)\right.+\left.S_{\mathrm{2}}^{\mathrm{*}}\left({\beta }_{\mathrm{3}}{I}^{\mathrm{*}}+{\beta }_{\mathrm{4}}{W}^{\mathrm{*}}\right)\right]\left(\mathrm{1}-c\right)+\phi E^{\mathrm{*}},\\ h{I}^{\mathrm{*}}=k{R}^{\mathrm{*}}+\delta R^{\mathrm{*}},\\ \sigma I^{\mathrm{*}}=d{W}^{\mathrm{*}}\mathrm{。}\end{array}\right.（2）

证明 由式（2）可知

\begin{array}{l}{W}^{\mathrm{*}}=\frac{\sigma }{d}{I}^{\mathrm{*}},R^{\mathrm{*}}=\frac{h}{k+\delta }{I}^{\mathrm{*}},S_{\mathrm{1}}^{\mathrm{*}}=\frac{dA}{d{\beta }_{\mathrm{1}}{I}^{\mathrm{*}}+\sigma {\beta }_{\mathrm{2}}{I}^{\mathrm{*}}+d\gamma +d\delta },\\ S_{\mathrm{2}}^{\mathrm{*}}=\frac{dA\gamma }{\left(d{\beta }_{\mathrm{1}}{I}^{\mathrm{*}}+\sigma {\beta }_{\mathrm{2}}{I}^{\mathrm{*}}+d\gamma +d\delta \right)\left(d{\beta }_{\mathrm{3}}{I}^{\mathrm{*}}+\sigma {\beta }_{\mathrm{4}}{I}^{\mathrm{*}}+d\delta \right)},\\ E^{\mathrm{*}}=\frac{cA(d{\beta }_{\mathrm{1}}{I}^{\mathrm{*}}+\sigma {\beta }_{\mathrm{2}}{I}^{\mathrm{*}})}{m}+\frac{kh{I}^{\mathrm{*}}}{(\phi +\delta )(k+\delta )}+\frac{cdA\gamma (d{\beta }_{\mathrm{3}}{I}^{\mathrm{*}}+\sigma {\beta }_{\mathrm{4}}{I}^{\mathrm{*}})}{m(d{\beta }_{\mathrm{3}}{I}^{\mathrm{*}}+\sigma {\beta }_{\mathrm{4}}{I}^{\mathrm{*}}+d\delta )}\mathrm{。}\end{array}

其中m=(\phi +\delta )(d{\beta }_{\mathrm{1}}{I}^{\mathrm{*}}+\sigma {\beta }_{\mathrm{2}}{I}^{\mathrm{*}}+d\gamma +d\delta )。

将S_{\mathrm{1}}^{\mathrm{*}}\mathrm{、}{S}_{\mathrm{2}}^{\mathrm{*}}\mathrm{、}{E}^{\mathrm{*}}\mathrm{、}{R}^{\mathrm{*}}\mathrm{、}{W}^{\mathrm{*}}代入式（2）的第3个式子， 可得

b_{\mathrm{1}}{I}^{\mathrm{*}\mathrm{2}}+b_{\mathrm{2}}{I}^{\mathrm{*}}+b_{\mathrm{3}}=\mathrm{0}。

其中

\begin{array}{l}{b}_{\mathrm{1}}=\frac{d(d{\beta }_{\mathrm{1}}+\sigma {\beta }_{\mathrm{2}})(d{\beta }_{\mathrm{3}}+\sigma {\beta }_{\mathrm{4}})}{k+\delta }[kh\phi -(\phi +\delta )(k+\delta )(h+\delta +\mu )],\\ b_{\mathrm{2}}=dA(\phi +\delta )(d{\beta }_{\mathrm{1}}+\sigma {\beta }_{\mathrm{2}})(d{\beta }_{\mathrm{3}}+\sigma {\beta }_{\mathrm{4}})-c\delta (d{\beta }_{\mathrm{1}}+\sigma {\beta }_{\mathrm{2}})(d{\beta }_{\mathrm{3}}+\sigma {\beta }_{\mathrm{4}})+\\ \frac{d^{\mathrm{2}}\delta kh\phi I^{\mathrm{*}}(d{\beta }_{\mathrm{1}}+\sigma {\beta }_{\mathrm{2}})}{k+\delta }+\frac{d^{\mathrm{2}}kh\phi I^{\mathrm{*}}(d{\beta }_{\mathrm{3}}+\sigma {\beta }_{\mathrm{4}})(\gamma +\delta )}{k+\delta },\\ b_{\mathrm{3}}=\frac{d^{\mathrm{3}}\delta (\gamma +\delta )}{k+\delta }\left(\mathrm{1}-{\mathrm{ℛ}}_{\mathrm{0}}\right)[hk\phi -(k+\delta )(\phi +\delta )(h+\delta +\mu )]\mathrm{。}\end{array}

由于， 当{\mathrm{ℛ}}_{\mathrm{0}}>\mathrm{1}时， ， b_{\mathrm{3}}>\mathrm{0}。因此， 模型（1）存在唯一的地方病平衡点P^{\mathrm{*}}。

2 模型(1)平衡点的稳定性分析

2.1 无病平衡点P^{\mathrm{0}}的稳定性

定理2 当时， 模型（1）的无病平衡点P^{\mathrm{0}}局部渐近稳定； 当{\mathrm{ℛ}}_{\mathrm{0}}>\mathrm{1}时， 模型（1）的无病平衡点P^{\mathrm{0}}不稳定。

证明 定义s(\mathit{M})=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\mathrm{R}\mathrm{e}\lambda :\lambda \right.是\mathit{M}的特征根\}， 其中\mathit{M}=\mathit{F}-\mathit{V}。 点P^{\mathrm{0}}处的Jacobi矩阵为

\mathit{J}=\left[\begin{array}{cc}{\mathit{J}}_{\mathrm{1}}& \mathrm{*}\\ \mathrm{0}& \mathit{M}\end{array}\right]。

其中：\mathrm{*}表示\mathrm{2}\times \mathrm{4}维非零矩阵，

{\mathit{J}}_{\mathrm{1}}=\left[\begin{array}{cc}-\gamma -\delta & \mathrm{0}\\ \gamma & -\delta \end{array}\right]。

因此， 矩阵\mathit{J}的特征方程为

|\lambda \mathit{E}-\mathit{J}|=(\lambda +\gamma +\delta )(\lambda +\delta )|\lambda \mathit{E}-\mathit{M}|=\mathrm{0}。（3）

显然， -\gamma -\delta和-\delta为特征方程（3）的两个负实根。对|\lambda \mathit{E}-\mathit{M}|=\mathrm{0}， 由文献 ［23］ 中定理2可知，当时， ， 故特征方程（3）所有特征值均有负实部。因此， 模型（1）的无病平衡点P^{\mathrm{0}}局部渐近稳定。当{\mathrm{ℛ}}_{\mathrm{0}}>\mathrm{1}时， s(\mathit{M})>\mathrm{0}， 这意味着特征方程（3）至少存在一个实部为正的特征值。因此， 模型（1）的无病平衡点P^{\mathrm{0}}是不稳定的。

定理3 当时， 模型（1）的无病平衡点P^{\mathrm{0}}全局渐近稳定。

证明 定义Lyapunov函数为

V_{\mathrm{1}}=E+\frac{\phi +\delta }{\phi }I+\frac{k}{k+\delta }R+\frac{c\phi ({\beta }_{\mathrm{2}}{S}_{\mathrm{1}}^{\mathrm{0}}+{\beta }_{\mathrm{4}}{S}_{\mathrm{2}}^{\mathrm{0}})}{d\phi }W+\frac{(\mathrm{1}-c)(\phi +\delta )({\beta }_{\mathrm{2}}{S}_{\mathrm{1}}^{\mathrm{0}}+{\beta }_{\mathrm{4}}{S}_{\mathrm{2}}^{\mathrm{0}})}{d\phi }W.

计算V_{\mathrm{1}}沿着模型（1）轨线的全导数， 得

\begin{array}{l}\frac{\mathrm{d}{V}_{\mathrm{1}}}{\mathrm{d}t}=\frac{\mathrm{d}E}{\mathrm{d}t}+\frac{\phi +\delta }{\phi }\frac{\mathrm{d}I}{\mathrm{d}t}+\frac{k}{k+\delta }\frac{\mathrm{d}R}{\mathrm{d}t}+\frac{c\phi ({\beta }_{\mathrm{2}}{S}_{\mathrm{1}}^{\mathrm{0}}+{\beta }_{\mathrm{4}}{S}_{\mathrm{2}}^{\mathrm{0}})}{d\phi }\frac{\mathrm{d}W}{\mathrm{d}t}+\frac{(\mathrm{1}-c)(\phi +\delta )({\beta }_{\mathrm{2}}{S}_{\mathrm{1}}^{\mathrm{0}}+{\beta }_{\mathrm{4}}{S}_{\mathrm{2}}^{\mathrm{0}})}{d\phi }\frac{\mathrm{d}W}{\mathrm{d}t}=\\ c{S}_{\mathrm{1}}^{\mathrm{0}}({\beta }_{\mathrm{1}}I+{\beta }_{\mathrm{2}}W)+c{S}_{\mathrm{2}}^{\mathrm{0}}({\beta }_{\mathrm{3}}I+{\beta }_{\mathrm{4}}W)+kR-\phi E-\delta E+\frac{k}{k+\delta }(hI-kR-\delta R)+\\ \frac{\phi +\delta }{\phi }[(\mathrm{1}-c)S_{\mathrm{1}}^{\mathrm{0}}({\beta }_{\mathrm{1}}I+{\beta }_{\mathrm{2}}W)+\phi E-hI+(\mathrm{1}-c)S_{\mathrm{2}}^{\mathrm{0}}({\beta }_{\mathrm{3}}I+{\beta }_{\mathrm{4}}W)-(\delta +\mu )I]+\\ \frac{c\phi ({\beta }_{\mathrm{2}}{S}_{\mathrm{1}}^{\mathrm{0}}+{\beta }_{\mathrm{4}}{S}_{\mathrm{2}}^{\mathrm{0}})}{d\phi }(\sigma I+dW)+\frac{(\mathrm{1}-c)(\phi +\delta )({\beta }_{\mathrm{2}}{S}_{\mathrm{1}}^{\mathrm{0}}+{\beta }_{\mathrm{4}}{S}_{\mathrm{2}}^{\mathrm{0}})}{d\phi }(\sigma I+dW)=\\ \frac{(k+\delta )(\phi +\delta )(h+\delta +\mu )-hk\phi }{\phi (\phi +\delta )}({\mathrm{ℛ}}_{\mathrm{0}}-\mathrm{1})I\mathrm{。}\end{array}

由LaSalle不变集原理可知， 当时， 模型（1）的无病平衡点P^{\mathrm{0}}全局渐近稳定。

2.2 地方病平衡点P^{\mathrm{*}}的稳定性

定理4 当{\mathrm{ℛ}}_{\mathrm{0}}>\mathrm{1}时， 模型（1）的地方病平衡点P^{\mathrm{*}}=(S_{\mathrm{1}}^{\mathrm{*}},S_{\mathrm{2}}^{\mathrm{*}},E^{\mathrm{*}},I^{\mathrm{*}},R^{\mathrm{*}},W^{\mathrm{*}})全局渐近稳定。

证明 令

s_{\mathrm{1}}=\frac{S_{\mathrm{1}}}{S_{\mathrm{1}}^{\mathrm{*}}},s_{\mathrm{2}}=\frac{S_{\mathrm{2}}}{S_{\mathrm{2}}^{\mathrm{*}}},e=\frac{E}{E^{\mathrm{*}}},i=\frac{I}{I^{\mathrm{*}}},r=\frac{R}{R^{\mathrm{*}}},w=\frac{W}{W^{\mathrm{*}}}，

则模型（1）可变换为

\left\{\begin{array}{l}\frac{\mathrm{d}{s}_{\mathrm{1}}}{\mathrm{d}t}=s_{\mathrm{1}}\left[\frac{A}{S_{\mathrm{1}}^{\mathrm{*}}}\left(\frac{\mathrm{1}}{s_{\mathrm{1}}}-\mathrm{1}\right)-{\beta }_{\mathrm{1}}{I}^{\mathrm{*}}(i-\mathrm{1})\right]-s_{\mathrm{1}}{\beta }_{\mathrm{2}}{W}^{\mathrm{*}}(w-\mathrm{1}),\\ \frac{\mathrm{d}{s}_{\mathrm{2}}}{\mathrm{d}t}=s_{\mathrm{2}}\left[\frac{\gamma S_{\mathrm{1}}^{\mathrm{*}}}{S_{\mathrm{2}}^{\mathrm{*}}}\left(\frac{s_{\mathrm{1}}}{s_{\mathrm{2}}}-\mathrm{1}\right)-{\beta }_{\mathrm{3}}{I}^{\mathrm{*}}(i-\mathrm{1})\right]-s_{\mathrm{2}}{\beta }_{\mathrm{4}}{W}^{\mathrm{*}}(w-\mathrm{1}),\\ \frac{\mathrm{d}e}{\mathrm{d}t}=ec{\beta }_{\mathrm{1}}\frac{S_{\mathrm{1}}^{\mathrm{*}}{I}^{\mathrm{*}}}{E^{\mathrm{*}}}\left(\frac{s_{\mathrm{1}}i}{e}-\mathrm{1}\right)+ec{\beta }_{\mathrm{2}}\frac{S_{\mathrm{1}}^{\mathrm{*}}{W}^{\mathrm{*}}}{E^{\mathrm{*}}}\left(\frac{s_{\mathrm{1}}w}{e}-\mathrm{1}\right)+e\frac{k{R}^{\mathrm{*}}}{E^{\mathrm{*}}}\left(\frac{r}{e}-\mathrm{1}\right)+\\ ec{\beta }_{\mathrm{4}}\frac{S_{\mathrm{2}}^{\mathrm{*}}{W}^{\mathrm{*}}}{E^{\mathrm{*}}}\left(\frac{s_{\mathrm{2}}w}{e}-\mathrm{1}\right)+ec{\beta }_{\mathrm{3}}\frac{S_{\mathrm{2}}^{\mathrm{*}}{I}^{\mathrm{*}}}{E^{\mathrm{*}}}\left(\frac{s_{\mathrm{2}}i}{e}-\mathrm{1}\right),\\ \frac{\mathrm{d}i}{\mathrm{d}t}=i(\mathrm{1}-c){\beta }_{\mathrm{1}}{S}_{\mathrm{1}}^{\mathrm{*}}(s_{\mathrm{1}}-\mathrm{1})+i(\mathrm{1}-c){\beta }_{\mathrm{2}}\frac{S_{\mathrm{2}}^{\mathrm{*}}{W}^{\mathrm{*}}}{I^{\mathrm{*}}}\left(\frac{s_{\mathrm{2}}w}{i}-\mathrm{1}\right)+i\frac{\phi E^{\mathrm{*}}}{I^{\mathrm{*}}}\left(\frac{e}{i}-\mathrm{1}\right)+\\ i(\mathrm{1}-c){\beta }_{\mathrm{4}}\frac{S_{\mathrm{2}}^{\mathrm{*}}{W}^{\mathrm{*}}}{I^{\mathrm{*}}}\left(\frac{s_{\mathrm{2}}w}{i}-\mathrm{1}\right)+i(\mathrm{1}-c){\beta }_{\mathrm{3}}{S}_{\mathrm{2}}^{\mathrm{*}}(s_{\mathrm{2}}-\mathrm{1}),\\ \frac{\mathrm{d}r}{\mathrm{d}t}=r\frac{h{I}^{\mathrm{*}}}{R^{\mathrm{*}}}\left(\frac{i}{r}-\mathrm{1}\right),\\ \frac{\mathrm{d}w}{\mathrm{d}t}=w\frac{\sigma I^{\mathrm{*}}}{W^{\mathrm{*}}}\left(\frac{i}{w}-\mathrm{1}\right)\mathrm{。}\end{array}\right.（4）

模型（4）具有唯一的地方病平衡点P_{\mathrm{1}}^{\mathrm{*}}=(\mathrm{1,1},\mathrm{1,1},\mathrm{1,1})， 并且P_{\mathrm{1}}^{\mathrm{*}}与P^{\mathrm{*}}的全局渐近稳定性等价， 所以只需证明P_{\mathrm{1}}^{\mathrm{*}}是全局渐近稳定的。

定义Lyapunov函数为

V_{\mathrm{2}}=S_{\mathrm{1}}^{\mathrm{*}}g(s_{\mathrm{1}})+S_{\mathrm{2}}^{\mathrm{*}}g(s_{\mathrm{2}})+E^{\mathrm{*}}g(e)+I^{\mathrm{*}}g(i)+\frac{k{R}^{\mathrm{*}\mathrm{2}}}{h{I}^{\mathrm{*}}}g(r)+\frac{({\beta }_{\mathrm{2}}{S}_{\mathrm{1}}^{\mathrm{*}}+{\beta }_{\mathrm{4}}{S}_{\mathrm{2}}^{\mathrm{*}})W^{\mathrm{*}\mathrm{2}}}{\sigma I^{\mathrm{*}}}g(w)。

其中：g(x)=x-\mathrm{1}-\mathrm{l}\mathrm{n}x(x>\mathrm{0})。显然， g(x)\ge \mathrm{0}， 并且g(x)=\mathrm{0}当且仅当x=\mathrm{1}， 故V_{\mathrm{2}}是正定的。计算V_{\mathrm{2}}沿着模型（4）轨线的全导数， 得

\begin{array}{l}\frac{\mathrm{d}{V}_{\mathrm{2}}}{\mathrm{d}t}=S_{\mathrm{1}}^{\mathrm{*}}\left(\mathrm{1}-\frac{\mathrm{1}}{s_{\mathrm{1}}}\right)\frac{\mathrm{d}{s}_{\mathrm{1}}}{\mathrm{d}t}+S_{\mathrm{2}}^{\mathrm{*}}\left(\mathrm{1}-\frac{\mathrm{1}}{s_{\mathrm{2}}}\right)\frac{\mathrm{d}{s}_{\mathrm{2}}}{\mathrm{d}t}+E^{\mathrm{*}}\left(\mathrm{1}-\frac{\mathrm{1}}{e}\right)\frac{\mathrm{d}e}{\mathrm{d}t}+I^{\mathrm{*}}\left(\mathrm{1}-\frac{\mathrm{1}}{i}\right)\frac{\mathrm{d}i}{\mathrm{d}t}+\\ \frac{({\beta }_{\mathrm{2}}{S}_{\mathrm{1}}^{\mathrm{*}}+{\beta }_{\mathrm{4}}{S}_{\mathrm{2}}^{\mathrm{*}})W^{\mathrm{*}\mathrm{2}}}{\sigma I^{\mathrm{*}}}\left(\mathrm{1}-\frac{\mathrm{1}}{w}\right)\frac{\mathrm{d}w}{\mathrm{d}t}+\frac{k{R}^{\mathrm{*}\mathrm{2}}}{h{I}^{\mathrm{*}}}\left(\mathrm{1}-\frac{\mathrm{1}}{r}\right)\frac{\mathrm{d}r}{\mathrm{d}t}\mathrm{。}\end{array}

由模型（4）得

\begin{array}{l}\frac{\mathrm{d}{V}_{\mathrm{2}}}{\mathrm{d}t}=-g\left(\frac{\mathrm{1}}{s_{\mathrm{1}}}\right)A-g\left(\frac{s_{\mathrm{1}}}{s_{\mathrm{2}}}\right)S_{\mathrm{1}}^{\mathrm{*}}\gamma -\left(\frac{s_{\mathrm{1}}i}{e}\right)c{\beta }_{\mathrm{1}}{S}_{\mathrm{1}}^{\mathrm{*}}{I}^{\mathrm{*}}-g\left(\frac{s_{\mathrm{1}}w}{i}\right)c{\beta }_{\mathrm{2}}{S}_{\mathrm{1}}^{\mathrm{*}}{W}^{\mathrm{*}}-g\left(\frac{s_{\mathrm{2}}i}{e}\right)c{\beta }_{\mathrm{3}}{S}_{\mathrm{2}}^{\mathrm{*}}{I}^{\mathrm{*}}-\\ g\left(\frac{s_{\mathrm{2}}w}{i}\right)(\mathrm{1}-c){\beta }_{\mathrm{4}}{S}_{\mathrm{2}}^{\mathrm{*}}{W}^{\mathrm{*}}-g\left(\frac{s_{\mathrm{2}}w}{e}\right)c{\beta }_{\mathrm{4}}{S}_{\mathrm{2}}^{\mathrm{*}}{W}^{\mathrm{*}}-g\left(\frac{e}{i}\right)E^{\mathrm{*}}\phi -g\left(\frac{i}{r}\right)h{I}^{\mathrm{*}}-g\left(\frac{r}{e}\right)k{R}^{\mathrm{*}}-\\ g\left(\frac{s_{\mathrm{1}}w}{e}\right)[(\mathrm{1}-c)\beta +\mathrm{2}{S}_{\mathrm{1}}^{\mathrm{*}}{W}^{\mathrm{*}}]\le \mathrm{0}\mathrm{。}\end{array}

由LaSalle不变集原理可知， 模型（4）的地方病平衡点P_{\mathrm{1}}^{\mathrm{*}}=(\mathrm{1,1},\mathrm{1,1},\mathrm{1,1})是全局渐近稳定的， 故模型（1）的地方病平衡点P^{\mathrm{*}}也是全局渐近稳定的。

3 模型的预设目标控制问题

本节运用最优控制的方法研究模型（1）的预设目标控制问题。

3.1 控制问题的建立与最优控制的存在性

引入两个控制变量u_{\mathrm{1}}(t)(\mathrm{0}\le u_{\mathrm{1}}(t)\le \mathrm{1})和u_{\mathrm{2}}(t)(\mathrm{0}\le u_{\mathrm{2}}(t)\le \mathrm{1})， 分别表示疫苗接种和实施DOTS治疗方案， 建立控制模型：

\left\{\begin{array}{l}\frac{\mathrm{d}{S}_{\mathrm{1}}}{\mathrm{d}t}=A-S_{\mathrm{1}}\left({\beta }_{\mathrm{1}}I+{\beta }_{\mathrm{2}}W\right)-\gamma S_{\mathrm{1}}-\delta S_{\mathrm{1}},\\ \frac{\mathrm{d}{S}_{\mathrm{2}}}{\mathrm{d}t}=\gamma S_{\mathrm{1}}-S_{\mathrm{2}}\left({\beta }_{\mathrm{3}}I+{\beta }_{\mathrm{4}}W\right)-\delta S_{\mathrm{2}}-u_{\mathrm{1}}(t)\upsilon S_{\mathrm{2}},\\ \frac{\mathrm{d}E}{\mathrm{d}t}=c{S}_{\mathrm{1}}({\beta }_{\mathrm{1}}I+{\beta }_{\mathrm{2}}W)+c{S}_{\mathrm{2}}({\beta }_{\mathrm{3}}I+{\beta }_{\mathrm{4}}W)+kR-\phi E-\delta E,\\ \frac{\mathrm{d}I}{\mathrm{d}t}=(\mathrm{1}-c)S_{\mathrm{1}}({\beta }_{\mathrm{1}}I+{\beta }_{\mathrm{2}}W)+\phi E+(\mathrm{1}-c)S_{\mathrm{2}}({\beta }_{\mathrm{3}}I+{\beta }_{\mathrm{4}}W)-hI-(\delta +\mu )I-u_{\mathrm{2}}(t)I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}=hI+u_{\mathrm{2}}(t)I-kR-\delta R,\\ \frac{\mathrm{d}W}{\mathrm{d}t}=\sigma I-dW\mathrm{。}\end{array}\right.（5）

其中：u_{\mathrm{1}}(t)=\mathrm{0}表示对易感成年人不接种疫苗， u_{\mathrm{1}}(t)=\mathrm{1}表示对易感成年人完全接种疫苗； u_{\mathrm{2}}(t)=\mathrm{0}表示对感染者不实施DOTS措施， u_{\mathrm{2}}(t)=\mathrm{1}表示对感染者完全实施DOTS措施； \upsilon表示疫苗的有效率， \upsilon u_{\mathrm{1}}(t)表示疫苗的有效接种率。初始条件为

S_{\mathrm{1}}(\mathrm{0})\ge \mathrm{0},S_{\mathrm{2}}(\mathrm{0})\ge \mathrm{0},E(\mathrm{0})\ge \mathrm{0},I(\mathrm{0})\ge \mathrm{0},R(\mathrm{0})\ge \mathrm{0},W(\mathrm{0})\ge \mathrm{0}。

令\mathit{x}={\left(S_{\mathrm{1}},S_{\mathrm{2}},E,I,R,W\right)}^{\mathrm{T}}， \mathit{u}={\left(u_{\mathrm{1}}(t),u_{\mathrm{2}}(t)\right)}^{\mathrm{T}}。 构造目标函数为

J(u_{\mathrm{1}}(t),u_{\mathrm{2}}(t))={\int }_{\mathrm{0}}^{T}L(\mathit{x},\mathit{u})\mathrm{d}t。（6）

其中

L=A_{\mathrm{1}}{S}_{\mathrm{2}}+A_{\mathrm{2}}{\left(I-I^{\mathrm{1}}\right)}^{\mathrm{2}}+B_{\mathrm{1}}{u}_{\mathrm{1}}(t)S_{\mathrm{2}}+B_{\mathrm{2}}{u}_{\mathrm{2}}(t)I+\frac{\mathrm{1}}{\mathrm{2}}\left(C_{\mathrm{1}}{u}_{\mathrm{1}}^{\mathrm{2}}(t)+C_{\mathrm{2}}{u}_{\mathrm{2}}^{\mathrm{2}}(t)\right)。

这里\left[\mathrm{0},T\right]表示实施两种不同控制措施的时间区间， 正常数A_i\mathrm{、}{B}_i和C_i(i=\mathrm{1},\mathrm{2})表示S_{\mathrm{2}}(t)\mathrm{、}I(t)\mathrm{、}{u}_{\mathrm{1}}(t)和u_{\mathrm{2}}(t)的加权系数。实施这两种控制措施的总目标是以最小的成本在有限的时间间隔\left[\mathrm{0},T\right]内使结核病感染者I(t)的数量稳定在预设目标I^{\mathrm{1}}之下， 也就是说， 要找到一个最优控制对{\left(u_{\mathrm{1}}^{\mathrm{*}}(t),u_{\mathrm{2}}^{\mathrm{*}}(t)\right)}^{\mathrm{T}}，使得J\left(u_{\mathrm{1}}^{\mathrm{*}}(t),u_{\mathrm{2}}^{\mathrm{*}}(t)\right)=\mathrm{m}\mathrm{i}\mathrm{n}\left\{J(u_{\mathrm{1}}(t),u_{\mathrm{2}}(t))|(u_{\mathrm{1}}(t),u_{\mathrm{2}}{(t))}^{\mathrm{T}}\in U\right\}， 其中U=\left\{(u_{\mathrm{1}}(t),u_{\mathrm{2}}{(t))}^{\mathrm{T}}|u_i(t)(i=\mathrm{1,2})\right.是Lebesgue可测的， u_i(t)\in [\mathrm{0,1}],\left.t\in [\mathrm{0},T]\right\}。

定理5 对控制模型（5）， 存在最优控制对{\mathit{u}}^{\mathrm{*}}={\left(u_{\mathrm{1}}^{\mathrm{*}}(t),u_{\mathrm{2}}^{\mathrm{*}}(t)\right)}^{\mathrm{T}}\in U以及对应的最优状态变量(S_{\mathrm{1}}^{\mathrm{1}\mathrm{*}},S_{\mathrm{2}}^{\mathrm{1}\mathrm{*}},E^{\mathrm{1}\mathrm{*}},I^{\mathrm{1}\mathrm{*}},R^{\mathrm{1}\mathrm{*}},W^{\mathrm{1}\mathrm{*}}{)}^{\mathrm{T}}， 使得J(u_{\mathrm{1}}^{\mathrm{*}}(t),u_{\mathrm{2}}^{\mathrm{*}}(t))=\mathrm{m}\mathrm{i}{\mathrm{n}}_{\left(u_{\mathrm{1}}(t),u_{\mathrm{2}}(t)\right)\in U}J(u_{\mathrm{1}}(t),u_{\mathrm{2}}(t))。

证明 令\mathit{p}=(p_{\mathrm{1}}(t),p_{\mathrm{2}}{(t))}^{\mathrm{T}},\mathit{q}=(q_{\mathrm{1}}(t),q_{\mathrm{2}}{(t))}^{\mathrm{T}}\in U， 对任意\kappa \in [\mathrm{0,1}]， 有

\kappa \mathit{p}+(\mathrm{1}-\kappa )\mathit{q}=(\kappa p_{\mathrm{1}}(t)+(\mathrm{1}-\kappa )q_{\mathrm{1}}(t),\kappa p_{\mathrm{2}}(t)+(\mathrm{1}-\kappa )q_{\mathrm{2}}{(t))}^{\mathrm{T}}\in U。

因为状态变量和控制变量都是非负的， 且U是封闭有界的， 所以U是凸集。根据正不变集\Omega可以推断， 对于每一个有界的\mathit{u}\in U， 模型（5）的解是有界的。由式（6）得

\begin{array}{l}L\left(\mathit{x},\kappa \mathit{p}+(\mathrm{1}-\kappa )\mathit{q}\right)-\kappa L(\mathit{x},\mathit{p})-(\mathrm{1}-\kappa )L(\mathit{x},\mathit{q})=\\ \frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum _{j=\mathrm{1}}^{\mathrm{2}}}\left\{B_j{\left[\kappa p_j(t)+(\mathrm{1}-\kappa )q_j(t)\right]}^{\mathrm{2}}\right\}-\frac{\mathrm{1}}{\mathrm{2}}\kappa {\displaystyle \sum _{j=\mathrm{1}}^{\mathrm{2}}}{B}_j{p}_j^{\mathrm{2}}(t)-\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{1}-\kappa ){\displaystyle \sum _{j=\mathrm{1}}^{\mathrm{2}}}{B}_j{q}_j^{\mathrm{2}}(t)=\\ \frac{\mathrm{1}}{\mathrm{2}}\kappa (\kappa -\mathrm{1}){\displaystyle \sum _{j=\mathrm{1}}^{\mathrm{2}}}{B}_j{\left(p_j(t)-q_j(t)\right)}^{\mathrm{2}}\le \mathrm{0},\end{array}

即L\left(\mathit{x},\kappa \mathit{p}+(\mathrm{1}-\kappa )\mathit{q}\right)\le \kappa L(\mathit{x},\mathit{p})+(\mathrm{1}-\kappa )L(\mathit{x},\mathit{q}),\forall \kappa \in [\mathrm{0,1}]。 因此， 目标函数的被积函数在控制集U上是凸函数。又因为存在常数\theta =\mathrm{2}， {\xi }_{\mathrm{1}}=\mathrm{0.5}\mathrm{m}\mathrm{i}\mathrm{n}\{B_{\mathrm{1}},B_{\mathrm{2}}\}>\mathrm{0}与{\xi }_{\mathrm{2}}=\mathrm{0}， 使得

A_{\mathrm{1}}{S}_{\mathrm{2}}+A_{\mathrm{2}}(I-I^{\mathrm{1}}{)}^{\mathrm{2}}+B_{\mathrm{1}}{u}_{\mathrm{1}}{S}_{\mathrm{2}}+B_{\mathrm{2}}{u}_{\mathrm{2}}I+\frac{\mathrm{1}}{\mathrm{2}}(C_{\mathrm{1}}{u}_{\mathrm{1}}^{\mathrm{2}}(t)+C_{\mathrm{2}}{u}_{\mathrm{2}}^{\mathrm{2}}(t))\ge {\xi }_{\mathrm{1}}{\left(|u_{\mathrm{1}}(t){|}^{\mathrm{2}}+|u_{\mathrm{2}}(t){|}^{\mathrm{2}}\right)}^{\frac{\theta }{\mathrm{2}}}-{\xi }_{\mathrm{2}}。

所以， 存在最优控制对{\mathit{u}}^{\mathrm{*}}=(u_{\mathrm{1}}^{\mathrm{*}}(t),u_{\mathrm{2}}^{\mathrm{*}}{(t))}^{\mathrm{T}}\in U， 满足J(u_{\mathrm{1}}^{\mathrm{*}}(t),u_{\mathrm{2}}^{\mathrm{*}}(t))=\mathrm{m}\mathrm{i}{\mathrm{n}}_{(u_{\mathrm{1}}(t),u_{\mathrm{2}}(t))\in U}J(u_{\mathrm{1}}(t),u_{\mathrm{2}}(t))。

3.2 控制问题的表征

定理6 设u_{\mathrm{1}}^{\mathrm{*}}(t)和u_{\mathrm{2}}^{\mathrm{*}}(t)为最优控制变量， S_{\mathrm{1}}^{\mathrm{1}\mathrm{*}}\mathrm{、}{S}_{\mathrm{2}}^{\mathrm{1}\mathrm{*}}\mathrm{、}E^{\mathrm{1}\mathrm{*}}\mathrm{、}{I}^{\mathrm{1}\mathrm{*}}\mathrm{、}{R}^{\mathrm{1}\mathrm{*}}\mathrm{、}{W}^{\mathrm{1}\mathrm{*}}是模型（5）在初始条件下对应的最优状态变量。则存在伴随变量\mathit{\lambda }(t)=({\lambda }_{\mathrm{1}}(t),{\lambda }_{\mathrm{2}}(t),{\lambda }_{\mathrm{3}}(t),{\lambda }_{\mathrm{4}}(t),{\lambda }_{\mathrm{5}}(t),{\lambda }_{\mathrm{6}}{(t))}^{\mathrm{T}}\in {\mathbf{R}}^{\mathrm{6}}满足下列伴随方程。即

\left\{\begin{array}{l}\begin{array}{l}\frac{\mathrm{d}{\lambda }_{\mathrm{1}}(t)}{\mathrm{d}t}={\lambda }_{\mathrm{1}}(t)({\beta }_{\mathrm{1}}I+{\beta }_{\mathrm{2}}W+\gamma +\delta )-{\lambda }_{\mathrm{2}}(t)\gamma -{\lambda }_{\mathrm{3}}(t)c({\beta }_{\mathrm{1}}I+{\beta }_{\mathrm{2}}W)-{\lambda }_{\mathrm{4}}(t)(\mathrm{1}-c)({\beta }_{\mathrm{1}}I+{\beta }_{\mathrm{2}}W),\\ \frac{\mathrm{d}{\lambda }_{\mathrm{2}}(t)}{\mathrm{d}t}=-A_{\mathrm{1}}-B_{\mathrm{1}}{u}_{\mathrm{1}}(t)+{\lambda }_{\mathrm{2}}(t)({\beta }_{\mathrm{3}}I+{\beta }_{\mathrm{4}}W+\delta +u_{\mathrm{1}}(t)\upsilon )-{\lambda }_{\mathrm{3}}(t)c({\beta }_{\mathrm{3}}I+{\beta }_{\mathrm{4}}W)-\\ {\lambda }_{\mathrm{4}}(t)(\mathrm{1}-c)({\beta }_{\mathrm{3}}I+{\beta }_{\mathrm{4}}W),\\ \frac{\mathrm{d}{\lambda }_{\mathrm{3}}(t)}{\mathrm{d}t}={\lambda }_{\mathrm{3}}(t)(\phi +\delta )-{\lambda }_{\mathrm{4}}(t)\phi ,\end{array}\\ \begin{array}{l}\frac{\mathrm{d}{\lambda }_{\mathrm{4}}(t)}{\mathrm{d}t}=-\mathrm{2}{A}_{\mathrm{2}}(I-I^{\mathrm{1}})-B_{\mathrm{2}}{u}_{\mathrm{2}}(t)+{\lambda }_{\mathrm{1}}(t)S_{\mathrm{1}}{\beta }_{\mathrm{1}}+{\lambda }_{\mathrm{2}}(t)S_{\mathrm{2}}{\beta }_{\mathrm{3}}-{\lambda }_{\mathrm{3}}(t)c(S_{\mathrm{1}}{\beta }_{\mathrm{1}}+S_{\mathrm{2}}{\beta }_{\mathrm{3}})-{\lambda }_{\mathrm{4}}(t)[(\mathrm{1}-c)S_{\mathrm{1}}{\beta }_{\mathrm{1}}+\\ (\mathrm{1}-c)S_{\mathrm{2}}{\beta }_{\mathrm{3}}-h-\delta -\mu -u_{\mathrm{2}}(t)]-{\lambda }_{\mathrm{5}}(t)(h+u_{\mathrm{2}}(t))-{\lambda }_{\mathrm{6}}(t)\sigma ,\\ \frac{\mathrm{d}{\lambda }_{\mathrm{5}}(t)}{\mathrm{d}t}=-{\lambda }_{\mathrm{3}}(t)k+{\lambda }_{\mathrm{5}}(t)(k+\delta ),\\ \frac{\mathrm{d}{\lambda }_{\mathrm{6}}(t)}{\mathrm{d}t}={\lambda }_{\mathrm{1}}(t)S_{\mathrm{1}}{\beta }_{\mathrm{2}}+{\lambda }_{\mathrm{2}}(t)S_{\mathrm{2}}{\beta }_{\mathrm{4}}-{\lambda }_{\mathrm{3}}(t)c(S_{\mathrm{1}}{\beta }_{\mathrm{2}}+S_{\mathrm{2}}{\beta }_{\mathrm{4}})-{\lambda }_{\mathrm{4}}(t)(\mathrm{1}-c)(S_{\mathrm{1}}{\beta }_{\mathrm{2}}+S_{\mathrm{2}}{\beta }_{\mathrm{4}})+{\lambda }_{\mathrm{6}}(t)d,\end{array}\end{array}\right.

其横截条件为{\lambda }_i(T)=\mathrm{0}(i=\mathrm{1,2},\cdots ,\mathrm{6})， 给出最优控制，即

u_i^{\mathrm{*}}(t)=\mathrm{m}\mathrm{i}\mathrm{n}\{\mathrm{m}\mathrm{a}\mathrm{x}\{D_i,\mathrm{0}\},\mathrm{1}\},i=\mathrm{1,2}。

其中：

D_{\mathrm{1}}=\frac{\upsilon S_{\mathrm{2}}^{\mathrm{1}\mathrm{*}}}{C_{\mathrm{1}}}({\lambda }_{\mathrm{2}}(t)-B_{\mathrm{1}}),D_{\mathrm{2}}=\frac{I^{\mathrm{1}\mathrm{*}}}{C_{\mathrm{2}}}({\lambda }_{\mathrm{4}}(t)-{\lambda }_{\mathrm{5}}(t)-B_{\mathrm{2}})。（7）

证明 为了方便， 将\mathit{\lambda }(t)=({\lambda }_{\mathrm{1}}(t),{\lambda }_{\mathrm{2}}(t),{\lambda }_{\mathrm{3}}(t),{\lambda }_{\mathrm{4}}(t),{\lambda }_{\mathrm{5}}(t),{\lambda }_{\mathrm{6}}{(t))}^{\mathrm{T}}简记为\mathit{\lambda }=({\lambda }_{\mathrm{1}},{\lambda }_{\mathrm{2}},{\lambda }_{\mathrm{3}},{\lambda }_{\mathrm{4}},{\lambda }_{\mathrm{5}},{\lambda }_{\mathrm{6}}{)}^{\mathrm{T}}，则模型 （5）的Hamiltonian函数为

\begin{array}{l}H=A_{\mathrm{1}}{S}_{\mathrm{2}}+A_{\mathrm{2}}(I-I^{\mathrm{1}}{)}^{\mathrm{2}}+B_{\mathrm{1}}{u}_{\mathrm{1}}{S}_{\mathrm{2}}+B_{\mathrm{2}}{u}_{\mathrm{2}}I+\frac{\mathrm{1}}{\mathrm{2}}(C_{\mathrm{1}}{u}_{\mathrm{1}}^{\mathrm{2}}+C_{\mathrm{2}}{u}_{\mathrm{2}}^{\mathrm{2}})+{\lambda }_{\mathrm{1}}A-{\lambda }_{\mathrm{1}}\delta S_{\mathrm{1}}-\\ {\lambda }_{\mathrm{1}}{S}_{\mathrm{1}}({\beta }_{\mathrm{1}}I+{\beta }_{\mathrm{2}}W)-{\lambda }_{\mathrm{1}}\gamma S_{\mathrm{1}}+{\lambda }_{\mathrm{2}}\gamma S_{\mathrm{1}}-{\lambda }_{\mathrm{2}}{S}_{\mathrm{2}}({\beta }_{\mathrm{3}}I+{\beta }_{\mathrm{4}}W)-{\lambda }_{\mathrm{2}}\delta S_{\mathrm{2}}-{\lambda }_{\mathrm{2}}{u}_{\mathrm{1}}(t)\upsilon S_{\mathrm{2}}+\\ {\lambda }_{\mathrm{3}}c{S}_{\mathrm{1}}({\beta }_{\mathrm{1}}I+{\beta }_{\mathrm{2}}W)+{\lambda }_{\mathrm{3}}c{S}_{\mathrm{2}}({\beta }_{\mathrm{3}}I+{\beta }_{\mathrm{4}}W)+{\lambda }_{\mathrm{3}}kR-{\lambda }_{\mathrm{3}}\phi E+{\lambda }_{\mathrm{4}}\phi E-{\lambda }_{\mathrm{4}}(\delta +\mu )I-\\ {\lambda }_{\mathrm{3}}\delta E+{\lambda }_{\mathrm{4}}(\mathrm{1}-c)S_{\mathrm{1}}({\beta }_{\mathrm{1}}I+{\beta }_{\mathrm{2}}W)+{\lambda }_{\mathrm{4}}(\mathrm{1}-c)S_{\mathrm{2}}({\beta }_{\mathrm{3}}I+{\beta }_{\mathrm{4}}W)-{\lambda }_{\mathrm{4}}hI-{\lambda }_{\mathrm{4}}{u}_{\mathrm{2}}(t)I+\\ {\lambda }_{\mathrm{5}}(hI+u_{\mathrm{2}}(t)I-kR-\delta R)+{\lambda }_{\mathrm{6}}(\sigma I-dW)\mathrm{。}\end{array}

依据Pontryagin极大值原理， 利用如下正则方程

\frac{\mathrm{d}{\lambda }_{\mathrm{1}}(t)}{\mathrm{d}t}=-\frac{\partial H}{\partial S_{\mathrm{1}}},\frac{\mathrm{d}{\lambda }_{\mathrm{2}}(t)}{\mathrm{d}t}=-\frac{\partial H}{\partial S_{\mathrm{2}}},\frac{\mathrm{d}{\lambda }_{\mathrm{3}}(t)}{\mathrm{d}\mathrm{t}}=-\frac{\partial H}{\partial E},\frac{\mathrm{d}{\lambda }_{\mathrm{4}}(t)}{\mathrm{d}t}=-\frac{\partial H}{\partial I},\frac{\mathrm{d}{\lambda }_{\mathrm{5}}(t)}{\mathrm{d}t}=-\frac{\partial H}{\partial R},\frac{\mathrm{d}{\lambda }_{\mathrm{6}}(t)}{\mathrm{d}t}=-\frac{\partial H}{\partial W},

可得最优控制u_{\mathrm{1}}^{\mathrm{*}}(t)和u_{\mathrm{2}}^{\mathrm{*}}(t)如（7）式所示。

4 数值模拟

本节将运用MATLAB进行数值模拟^［24-26］， 对前述理论结果的准确性进行验证。 进一步， 通过计算不同控制策略对应的感染避免率 （IAR）， 评估各种策略对结核病传播的抑制效果。

4.1 模型(1)平衡点的稳定性

模型（1）的初始值为S_{\mathrm{1}}(\mathrm{0})=\mathrm{200},S_{\mathrm{2}}(\mathrm{0})=\mathrm{500},E(\mathrm{0})=\mathrm{10},I(\mathrm{0})=\mathrm{3},R(\mathrm{0})=\mathrm{2},W(\mathrm{0})=\mathrm{40}， 参数取值为A=\mathrm{15},\delta =\mathrm{0.1},c=\mathrm{0.9},\phi =\mathrm{0.15},{\beta }_{\mathrm{1}}=\mathrm{0.001}\mathrm{6},{\beta }_{\mathrm{2}}=\mathrm{0.000}\mathrm{15},{\beta }_{\mathrm{3}}=\mathrm{0.003}\mathrm{6},{\beta }_{\mathrm{4}}=\mathrm{0.000}\mathrm{8},h=\mathrm{0.4},\sigma =\mathrm{1.8},d=\mathrm{0.8},k=\mathrm{0.15},\mu =\mathrm{0.3},\upsilon=\mathrm{0.8}^［5，27］。模型（1）随着时间变化的曲线图如图2所示。由图2（a）可以明显看出， 当时， 随着时间的推移， (S_{\mathrm{1}},S_{\mathrm{2}},E,I,R,W)收敛至无病平衡点P^{\mathrm{0}}=(\mathrm{120.8,27.3,0},\mathrm{0,0},\mathrm{0}), 这表明该疾病最终将在人群中消亡。由图2（b）可以看出， 当{\mathrm{ℛ}}_{\mathrm{0}}=\mathrm{2.830}\mathrm{6}>\mathrm{1}时， 随着时间的推移， (S_{\mathrm{1}},S_{\mathrm{2}},E,I,R,W)将收敛至地方病平衡点P^{\mathrm{*}}=(\mathrm{103.4,147.7,167.3,41.3,66.2,82.7})， 这表明该疾病最终会在人群中流行。数值模拟结果验证了定理3和定理4的准确性。

4.2 敏感性分析

为了分析模型（1）中各参数对结核病传播的影响， 采用偏秩相关法 （PRCC） 对基本再生数{\mathrm{ℛ}}_{\mathrm{0}}进行参数敏感性分析。仿真结果如图 3 所示， 据此可得，参数A （相对变化率为0.950）、σ（相对变化率为0.857）、{\beta }_{\mathrm{4}}（相对变化率为0.745）与{\mathrm{ℛ}}_{\mathrm{0}}呈强正相关， 是加剧结核病传播风险的关键驱动因素；φ（相对变化率为0.113）、k（相对变化率为0.089）等参数对结核病传播的影响较小。参数δ （相对变化率为-0.682）、d（相对变化率为-0.440）、h（相对变化率为-0.414）与{\mathrm{ℛ}}_{\mathrm{0}}呈较强负相关， 表明此类参数取值增大可抑制{\mathrm{ℛ}}_{\mathrm{0}}升高，进而降低疾病的传播能力。

4.3 不同控制策略对结核病传播的影响

基于两种干预措施u_{\mathrm{1}}(t)和u_{\mathrm{2}}(t)， 设计了4种差异化控制策略。策略1： u_{\mathrm{1}}=\mathrm{0},u_{\mathrm{2}}=\mathrm{0}； 策略2： u_{\mathrm{1}}\ne \mathrm{0},u_{\mathrm{2}}=\mathrm{0}； 策略3： u_{\mathrm{1}}=\mathrm{0},u_{\mathrm{2}}\ne \mathrm{0}； 策略4： u_{\mathrm{1}}\ne\mathrm{0},u_{\mathrm{2}}\ne \mathrm{0}。为了比较4种策略的控制效果， 将活动性结核病患者的控制目标设置为I^{\mathrm{1}}=\mathrm{20} （以不施加控制时活动性结核病患者数量的50%为目标）。

图 4 分别刻画了不同控制策略对模型（1）各状态变量的影响。从图4（a）可以观察到， 策略2、策略3和策略4对S_{\mathrm{1}}的影响差异较小； 而图4（d） 显示， 策略2导致康复者数量减少， 策略3则使康复者数量呈增长趋势。 这两种截然不同的结果， 源于各控制策略的干预措施针对特定人群实施，干预对象的差异导致了数量变化趋势的不同。在图4（b）中可见， 策略2可以减少易感人群S_{\mathrm{2}}的数量， 而策略3和策略4可以使易感人群S_{\mathrm{2}}的数量有所增加， 该结果符合实际情况。进一步地， 从图 4（c）、4（e） 和 4（f） 中可以看到， 实施策略2、策略3和策略4时， 不仅能降低各状态变量的峰值， 还可以延缓疫情高峰的到来， 这为传染病防控争取了关键的准备时间。在图4 （f）中还可看到，实施策略1时，活动性结核病患者数量先快速上升至峰值，随后缓慢下降并趋于稳定， 而实施策略4时， 活动性结核病患者数量相较于策略1上升更缓慢， 最终稳定在控制目标I^{\mathrm{1}}之下。在图4（e）中可见，策略4的控制效果显著优于策略2和策略3， 原因在于实施策略4时，活动性结核病患者人数减少，使得环境中的结核杆菌含量低于实施其他策略时的水平。图 5展示了实施策略2、策略3和策略4时，两种控制措施强度随着时间的变化情况。

4.4 控制措施效益分析

感染避免率 （IAR） 作为评估防控策略有效性的关键指标，其定义为

\mathrm{感染}\mathrm{避免}\mathrm{率}=\frac{\mathrm{避免}\mathrm{感染}\mathrm{的人}\mathrm{数}}{\mathrm{康复}\mathrm{人数}}\times \mathrm{100}\mathrm{\%}。

其中 “避免感染的人数” 指不实施控制措施时总感染人数与实施控制措施时总感染人数的差值。 该指标的比值越大，表明对应防控策略越有效。各策略的感染避免率计算结果（表2）显示：策略1的IAR为0， 策略2的IAR为3.32%， 策略3的IAR为36.03%， 策略4的IAR达到40.01%。这表明， 实施策略4时， 活动性结核病患者人数最少。

根据第4.3节不同控制策略的干预效果分析及第4.4节的效益评估结果可知，在资源有限的情况下，应优先选择以疫苗接种联合DOTS措施的综合防控策略。该策略凭借最高的感染避免率 IAR， 能够将活动性结核病患者数量有效控制在预设防控目标I^{\mathrm{1}}范围内。 在此基础上， 可结合阶段性防控成效与现有资源， 动态调整并制定下一阶段的量化控制指标， 逐步压缩疾病的传播空间。这种精准适配资源承载力的阶梯式防控模式， 既能确保各项控制措施有效落地， 又能稳步推进 “终结结核病流行” 终极目标的实现， 为资源约束条件下的结核病防控提供了兼具可行性与科学性的实施框架。

5 结论

通过分析结核病的传播特点，首先构建年龄-环境耦合的两阶段结核病传播动力学模型，并分析了模型无病平衡点与地方病平衡点的稳定性。其次，引入疫苗接种与DOTS控制措施，建立最优控制模型，运用Pontryagin极大值原理，开展模型的预设目标控制研究。最后，通过数值模拟验证了不同控制策略下的防控效果。结果表明，两种控制措施联合实施对结核病传播的防控效果最优。此外，结核病的耐药性、时滞等因素同样是影响疾病传播的重要因素，将此类因素纳入模型可进一步拓展模型的适用范围，这也是本研究未来的重点方向。

参考文献

[1] BLOWER S M， DALEY C L. Problems and solutions for the stop TB partnership［J］. The Lancet Infectious Diseases， 2002， 2（6）： 374-376.

[2] JACKSON S， SLEIGH A C， WANG Guojie， et al. Poverty and the economic effects of TB in rural China［J］. The International Journal of Tuberculosis and Lung Disease， 2006， 10（10）：1104-1110.

[3] BAGCCHI S. WHO’s global tuberculosis report 2022［J］. The Lancet Microbe， 2023， 4（1）： e20. DOI： 10.1016/S2666- 5247（22）00359-7 .

[4] WANG Lei， TENG Zhidong， RIFHAT R， et al. Modelling of a drug resistant tuberculosis for the contribution of resistance and relapse in Xinjiang， China［J］. Discrete and Continuous Dynamical Systems：B， 2023， 28（7）： 4167-4189.

[5] ZHANG Jun， TAKEUCHI Y， DONG Yueping， et al. Modelling the preventive treatment under media impact on tuberculosis： A comparison in four regions of China［J］. Infectious Disease Modelling， 2024， 9（2）： 483-500.

[6] SONG Pengfei， XIAO Yanni. Analysis of an epidemic system with two response delays in media impact function［J］. Bulletin of Mathematical Biology，2019，81（5）：1582-1612.

[7] ERNST J D. The immunological life cycle of tuberculosis［J］. Nature Reviews Immunology， 2012， 12（8）： 581-591.

[8] DING Zuqin， LI Yaxiao， WANG Xiameng， et al. The impact of air pollution on the transmission of pulmonary tuberculosis［J］. Mathematical Biosciences and Engineering， 2020， 17（4）： 4317-4327.

[9] 张正斌，鲁周琴，谢红，等.结核病季节性分布特征及影响因素［J］.中华流行病学杂志，2016，37（8）：1183-1186.

[10] CAI Yongli， ZHAO Shi， NIU Yun， et al. Modelling the effects of the contaminated environments on tuberculosis in Jiangsu， China［J］. Journal of Theoretical Biology，2021，508：110453.DOI：10.1016/j.jtbi.2020.110453 .

[11] LI Qiuyun， WANG Fengna. An epidemiological model for tuberculosis considering environmental transmission and reinfection［J］.Mathematics，2023，11（11）： 2423. DOI：10.3390/math11112423 .

[12] AGNELLI J P， BUFFA B， KNOPOFF D， et al. A spatial kinetic model of crowd evacuation dynamics with infectious disease contagion［J］. Bulletin of Mathematical Biology， 2023， 85（4）： 23.DOI：10.1007/s11538-023- 01127-6 .

[13] SHI Lei， QI Longxing. Dynamic analysis and optimal control of a class of SISP respiratory diseases［J］. Journal of Biological Dynamics， 2022， 16（1）： 64-97.

[14] 杨应周.关注脆弱人群的结核病防控［J］.中国防痨杂志，2013，35（11）：868-870.

[15] 张立兴，屠德华，安燕生，等.北京市结核病发病趋势研究［J］.中国防痨杂志，2003，25（4）：204-208.

[16] JING Shuanglin， XUE Ling， WANG Hao，et al. Global analysis of an age-structured tuberculosis model with an application to Jiangsu， China［J］. Journal of Mathematical Biology， 2024， 88（5）： 52. DOI：10.1007/s00285- 024-02066-z .

[17] GAO Chunjie， ZHANG Tao， LIAO Ying， et al. Modelling of tuberculosis dynamics incorporating indirect transmission of contaminated environment and infectivity of smear-negative individuals： A case study for Xinjiang， China［J］. Acta Tropica， 2024， 254： 107130.DOI：10.1016/j.actatropica.2024.107130 .

[18] XUE Ling， JING Shuanglin， WANG Hao.Evaluating strategies for tuberculosis to achieve the goals of WHO in China： A seasonal age-structured model study［J］. Bulletin of Mathematical Biology， 2022， 84（6）： 61. DOI： 10.1007/s11538-022- 01019-1 .

[19] KHODA P， BAJIYA V PAL， PRASAD S N. Effective strategies toward controlling tuberculosis： Optimal control and cost-effectiveness analysis［J］. The European Physical Journal Plus，2025，140（1）：14.DOI：10.1140/ epjp/s13360-025-05978-x .

[20] OBSU L L. Optimal control analysis of a tuberculosis model［J］. Journal of Biological Systems，2022，30（4）： 837-855.

[21] HUANG Wei， FANG Zhixiong， LUO Si， et al. The effect of BCG vaccination and risk factors for latent tuberculosis infection among college freshmen in China［J］. International Journal of Infectious Diseases， 2022， 122： 321-326.

[22] IMPERIAL M Z， NAHID P， PHILLIPS P P J， et al. A patient-level pooled analysis of treatment-shortening regimens for drug-susceptible pulmonary tuberculosis［J］. Nature Medicine， 2018， 24（11）： 1708-1715.

[23] VAN DEN DRIESSCHE P， WATMOUGH J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission［J］. Mathematical Biosciences， 2002， 180（1/2）： 29-48.

[24] 亢婷.随机年龄结构固定资产系统倒向Euler法的p阶矩耗散性［J］.宁夏大学学报（自然科学版），2024，45（1）：9-15.

[25] 曹博强，亢婷.两阶段分数阶羊布鲁氏菌病传播模型的非线性自适应控制［J］.应用数学，2025，38（3）：703-710.

[26] 张月蕾，朱磊，朱家明.基于DEA的单车共享经济的计量分析［J］.宁夏大学学报（自然科学版），2018，39（2）：115-120.

[27] RONOH M， JAROUDI R， FOTSO P， et al. A mathematical model of tuberculosis with drug resistance effects［J］. Applied Mathematics，2016，7（12）： 1303-1316.

基金资助

国家自然科学基金资助项目(12302017)

国家自然科学基金资助项目(12201330)

宁夏自然科学基金资助项目(2024AAC03001)

宁夏自然科学基金资助项目(2025AAC030001)