具有意识影响的传染病模型分析及优化控制

Model Analysis and Optimal Control of an Epidemic Model with the Effects of Awareness

摘要

文章首先建立了一类具有高、低风险的易感者和有、无意识的感染者的传染病模型。其次，通过下一代矩阵法求得基本再生数，利用Hurwitz判据和LaSalle不变集原理分析了平衡点的稳定性，结合中心流形定理证明系统存在向后分岔。然后，运用庞特里亚金（Pontryagin）最大值原理研究了模型对应的优化控制问题。最后，利用数值模拟对上述理论进行了验证，同时对参数进行了灵敏度分析，讨论了参数对基本再生数的影响。结果表明，增强感染者的意识和对高危易感者接种疫苗，能够有效控制传染病的传播。

Abstract

This paper first establishes a class of infectious disease model with high and low risk susceptible persons and infected persons who are unconscious. Secondly, the basic regeneration number is obtained by the next generation matrix method, the stability of the equilibrium point is analyzed by Hurwitz criterion and LaSalle invariant set principle, and the backward bifurcation is proved by the central prevalence theorem. Then, the optimal control problem corresponding to the model is studied using the Pontryagin maximum principle. Finally, the above theory is verified by numerical simulation, and the sensitivity analysis of parameters is carried out, and the influence of parameters on the basic regeneration number is discussed. The results showed that enhancing the awareness of infected people and vaccinating high-risk susceptible people could effectively control the spread of infectious diseases.

Graphical abstract

关键词

意识计划 / 稳定性 / 高风险人群 / 灵敏度分析

Key words

awareness programs / stability / high-risk population / sensitivity analysis

0 引言

传染病的传播严重危害人类健康，给社会造成了巨大的生态、经济损失。例如，1918年流感大流行在全世界造成数千万人死亡^［1］。SARS（严重急性呼吸系统综合征）的暴发导致了大规模的公共卫生危机^［2-3］。新型冠状病毒感染仍然是世界上最具破坏性的传染病之一，给社会和人民带来了深远的影响^［4-5］。

对于传染病，易感者中会存在一些特定的高风险人群，比如医护人员。同时，研究表明，流感的高危人群主要包括妊娠期妇女、儿童、65岁以上老年人、慢性基础性疾病患者等^［6］。Tylicki等^［7］表明，老年人、男性、高血压患者和并发症患者构成了新型冠状病毒感染严重表现最典型的高危人群。控制措施在缓解传染病传播方面发挥着重要作用。关于控制效果的研究经常使用最优控制方法，该方法旨在以最低的成本将感染人数降至最低^［8］。

人类对传染病的个人意识以及他们采取有效的预防措施以减少疾病传播的意愿也是传染病控制措施的一个重要理念^［9］。当传染病在人群中暴发时，有防范意识的感染个体产生警觉，采取自我保护措施（如戴口罩、居家隔离、取消行程等），减少与易感者接触，从而有效控制疾病的传播。许多学者使用动力学模型来描述个人意识对疾病传播的影响。Ndii等^［8］将易感人群分为有意识易感者和无意识易感者，使用数学模型和最优控制方法研究了意识计划和病媒控制对疟疾传播动力学的影响。模型结果表明，个人采取预防措施的意识可以有效减少与蚊子的接触，从而间接控制疾病的发生。Kar等^［10］研究了媒体意识对疾病传播动态的数学模型，该模型考虑了由于媒体报道的影响，有意识的易感个体转移到恢复者类。研究发现，通过媒体宣传活动增强意识可以有效地降低疾病的发生率，从而有效抑制疾病的传播。Musa等^［11］提出并分析了一个关于尼日利亚病毒传播动态的数学模型，该模型结合了针对轻度和重度病例的意识计划和不同的住院策略，以评估个人意识对新型冠状病毒感染传播动态的影响。研究结果表明，增强意识计划可以有效控制和缓解流行病。Aldila^［12］提出的模型考虑了通过媒体活动的干预以增强人类的意识，以及通过快速测试以跟踪现场未发现的病例。研究结果表明媒体宣传和快速检测策略相结合，能有效减少了感染者的数量。Kumar等^［13］提出了一个非线性数学模型来评估社交媒体广告在提高人们对传染病传播的认识和限制治疗方面的效果。研究结果表明，增加媒体和治疗的强度对控制疾病传播至关重要。Al Basir等^［14］调查了在没有适当疫苗的情况下，增强意识运动和最佳控制技术对传染病的影响。结果表明，宣传活动对控制传染病至关重要，而最优控制理论与媒体意识相结合，是控制传染病的必要策略。上述研究意识影响的数学模型是将易感人群分为有意识易感者和无意识易感者或将有意识的易感个体转移到恢复者类来制定的。

感染者的意识水平直接影响其对预防措施的采纳程度。意识到自己可能携带病毒的个体更可能采取诸如佩戴口罩、保持社交距离、勤洗手等个人防护措施，从而降低传染给他人的风险。具有较高健康意识的感染者在出现症状时可能会更快地寻求医疗帮助，并及时报告自己的病情，这对于早期识别疫情、追踪接触者以及控制疾病的进一步扩散至关重要。所以在传染病模型中纳入感染者的意识因素，有助于更加准确地模拟疾病传播过程，并为科学防控决策提供支持。例如：Ibrahim等^［15］建立了疟疾的数学模型，并研究了意识对疟疾传播动力学的影响。该模型中将受感染者分为有意识和无意识的受感染者。研究发现，这种认识是有助于减少疟疾发病率的一个重要因素。所以，本文考虑了感染者的意识，将感染人群分为有意识感染者和无意识感染者。同时考虑到易感者中存在一些特定的高风险人群，将易感人群分为高风险易感人群和低风险易感人群。识别并区分出高危易感人群有助于公共卫生部门制定针对性地预防和干预策略。例如，在疫苗分配上，可以优先保障高危群体。因此建立了一个具有意识影响的S₁S₂EI₁I₂R模型，并利用微分方程定性理论分析了传染病模型的动态问题，此外，通过最优控制方法研究了疫苗接种和增强意识计划对传染病传播动力学的影响。通过实施有效的控制措施，可以最大限度地减少感染人数，从而有效控制传染病的传播。

1 模型构建

将传染病传播区域的总人数分为高危易感者（S_{\mathrm{1}}）、低危易感者（S_{\mathrm{2}}）、潜伏者（E）、无意识感染者（I_{\mathrm{1}}）、有意识感染者（I_{\mathrm{2}}）、康复者（R），其中N=S_{\mathrm{1}}+S_{\mathrm{2}}+E+I_{\mathrm{1}}+I_{\mathrm{2}}+R。根据病毒的传播机理，构建如下具有意识影响的S₁S₂EI₁I₂R传染病传播动力学模型：

\left\{\begin{array}{l}\frac{\mathrm{d}{S}_{\mathrm{1}}}{\mathrm{d}t}=(\mathrm{1}-\tau )\Lambda -{\beta }_{\mathrm{1}}{S}_{\mathrm{1}}(I_{\mathrm{1}}+\sigma I_{\mathrm{2}})-\mu S_{\mathrm{1}},\\ \frac{\mathrm{d}{S}_{\mathrm{2}}}{\mathrm{d}t}=\tau \Lambda -{\beta }_{\mathrm{2}}{S}_{\mathrm{2}}(I_{\mathrm{1}}+\sigma I_{\mathrm{2}})-\mu S_{\mathrm{2}}+r_{\mathrm{3}}R,\\ \frac{\mathrm{d}E}{\mathrm{d}t}=({\beta }_{\mathrm{1}}{S}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{S}_{\mathrm{2}})(I_{\mathrm{1}}+\sigma I_{\mathrm{2}})-\mu E-\alpha E,\\ \frac{\mathrm{d}{I}_{\mathrm{1}}}{\mathrm{d}t}=p\alpha E-r_{\mathrm{1}}{I}_{\mathrm{1}}-\mu I_{\mathrm{1}},\\ \frac{\mathrm{d}{I}_{\mathrm{2}}}{\mathrm{d}t}=(\mathrm{1}-p)\alpha E-r_{\mathrm{2}}{I}_{\mathrm{2}}-\mu I_{\mathrm{2}}-d{I}_{\mathrm{2}},\\ \frac{\mathrm{d}R}{\mathrm{d}t}=r_{\mathrm{1}}{I}_{\mathrm{1}}+r_{\mathrm{2}}{I}_{\mathrm{2}}-\mu R-r_{\mathrm{3}}R,\end{array}\right.（1）

其中： \Lambda表示易感者的常数输入率，\tau表示低危易感者人群的比例，{\beta }_{\mathrm{1}}表示高危易感者与感染者的接触感染率，{\beta }_{\mathrm{2}}表示低危易感者与感染者的接触感染率，\sigma表示有意识感染者的修正系数，p表示无意识感染者人群的比例，\alpha表示从暴露到感染级的进展率，r_{\mathrm{1}}和r_{\mathrm{2}}分别表示无意识感染者和有意识感染者的恢复率，r_{\mathrm{3}}表示康复者重返低危易感者人群的比例，\mu表示自然死亡率，d表示因病死亡率。

假设模型（1）满足的初始条件为S_{\mathrm{1}}(\mathrm{0})\ge \mathrm{0}，S_{\mathrm{2}}(\mathrm{0})\ge \mathrm{0}，E(\mathrm{0})\ge \mathrm{0}，I_{\mathrm{1}}(\mathrm{0})\ge \mathrm{0}，I_{\mathrm{2}}(\mathrm{0})\ge \mathrm{0}，R_{\mathrm{0}}\ge \mathrm{0}，其中N=S_{\mathrm{1}}+S_{\mathrm{2}}+E+I_{\mathrm{1}}+I_{\mathrm{2}}+R。由于\frac{\mathrm{d}N}{\mathrm{d}t}=\Lambda -\mu N-d{I}_{\mathrm{2}}，则模型（1）的可行域为

\Omega =\left\{(S_{\mathrm{1}},S_{\mathrm{2}},E,I_{\mathrm{1}},I_{\mathrm{2}},R)\in {\mathfrak{R}}_{+}^{\mathrm{6}}|S_{\mathrm{1}}+S_{\mathrm{2}}+E+I_{\mathrm{1}}+I_{\mathrm{2}}+R\le \frac{\Lambda }{\mu }\right\}\mathrm{。}

2 基本再生数和平衡点的存在性

易知模型（1）的一个无病平衡点P_{\mathrm{0}}=\left(\frac{(\mathrm{1}-\tau )\Lambda }{\mu },\frac{\tau \Lambda }{\mu },\mathrm{0,0},\mathrm{0,0}\right)。利用下一代矩阵法 ^［16］，计算基本再生数。

\mathit{F}=\left(\begin{array}{ccc}\mathrm{0}& \frac{{\beta }_{\mathrm{1}}(\mathrm{1}-\tau )\Lambda }{\mu }+\frac{{\beta }_{\mathrm{2}}\tau \Lambda }{\mu }& (\frac{{\beta }_{\mathrm{1}}(\mathrm{1}-\tau )\Lambda }{\mu }+\frac{{\beta }_{\mathrm{2}}\tau \Lambda }{\mu })\sigma \\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)， \mathit{V}=\left(\begin{array}{ccc}\mu +\alpha & \mathrm{0}& \mathrm{0}\\ -p\alpha & r_{\mathrm{1}}+\mu & \mathrm{0}\\ -(\mathrm{1}-p)\alpha & \mathrm{0}& r_{\mathrm{2}}+\mu +d\end{array}\right)，

经过推导，基本再生数为：

R_{\mathrm{0}}=\rho (\mathit{F}{\mathit{V}}^{-\mathrm{1}})=\frac{\left[{\beta }_{\mathrm{1}}(\mathrm{1}-\tau )\Lambda +{\beta }_{\mathrm{2}}\tau \Lambda \right]\left[p\alpha (r_{\mathrm{2}}+\mu +d)+\sigma (\mathrm{1}-p)\alpha (r_{\mathrm{1}}+\mu )\right]}{\mu (\mu +\alpha )(r_{\mathrm{1}}+\mu )(r_{\mathrm{2}}+\mu +d)}。

设模型（1）存在正平衡点P^{\mathrm{*}}=(S_{\mathrm{1}}^{\mathrm{*}},S_{\mathrm{2}}^{\mathrm{*}},E^{\mathrm{*}},I_{\mathrm{1}}^{\mathrm{*}},I_{\mathrm{2}}^{\mathrm{*}},R^{\mathrm{*}})，令模型（1）中右面等于0，求解得

E^{\mathrm{*}}=\frac{m}{p\alpha }{I}_{\mathrm{1}}^{\mathrm{*}}，I_{\mathrm{2}}^{\mathrm{*}}=\frac{m(\mathrm{1}-p)}{np}{I}_{\mathrm{1}}^{\mathrm{*}}，R^{\mathrm{*}}=\frac{r_{\mathrm{1}}pn+r_{\mathrm{2}}(\mathrm{1}-p)m}{pqn}{I}_{\mathrm{1}}^{\mathrm{*}}，S_{\mathrm{1}}^{\mathrm{*}}=\frac{(\mathrm{1}-\tau )\Lambda pn\alpha }{{\beta }_{\mathrm{1}}Q{I}_{\mathrm{1}}^{\mathrm{*}}+\mu np\alpha }，

S_{\mathrm{2}}^{\mathrm{*}}=\frac{gmn({\beta }_{\mathrm{1}}Q{I}_{\mathrm{1}}^{\mathrm{*}}+\mu np\alpha )-{\beta }_{\mathrm{1}}(\mathrm{1}-\tau )\Lambda pn\alpha Q}{{\beta }_{\mathrm{2}}Q({\beta }_{\mathrm{1}}Q{I}_{\mathrm{1}}^{\mathrm{*}}+\mu np\alpha )}=\frac{\alpha npq\tau \Lambda +r_{\mathrm{3}}M{I}_{\mathrm{1}}^{\mathrm{*}}}{{\beta }_{\mathrm{2}}Q{I}_{\mathrm{1}}^{\mathrm{*}}+\mu npq\alpha }，

其中：

g=\mu +\alpha，m=r_{\mathrm{1}}+\mu，n=r_{\mathrm{2}}+\mu +d，q=\mu +r_{\mathrm{3}}，

Q=\alpha pn+\sigma (\mathrm{1}-p)m\alpha，M=r_{\mathrm{1}}pn\alpha +r_{\mathrm{2}}(\mathrm{1}-p)m\alpha，T=\frac{{\beta }_{\mathrm{1}}(\mathrm{1}-\tau )\Lambda +{\beta }_{\mathrm{2}}\tau \Lambda }{\mu }。

I_{\mathrm{1}}^{\mathrm{*}}的值由方程（2）的解给出，

c_{\mathrm{2}}{I_{\mathrm{1}}^{\mathrm{*}}}^{\mathrm{2}}+c_{\mathrm{1}}{I}_{\mathrm{1}}^{\mathrm{*}}+c_{\mathrm{0}}=\mathrm{0},（2）

其中：

c_{\mathrm{2}}=Q^{\mathrm{2}}{\beta }_{\mathrm{1}}{\beta }_{\mathrm{2}}p\alpha (r_{\mathrm{3}}M-qgmn)，c_{\mathrm{1}}=q{\beta }_{\mathrm{1}}{\beta }_{\mathrm{2}}{Q}^{\mathrm{2}}{p}^{\mathrm{2}}{\alpha }^{\mathrm{2}}\Lambda n\left(\mathrm{1}-\frac{{\beta }_{\mathrm{1}}(\mathrm{1}-\tau )+{\beta }_{\mathrm{2}}\tau }{{\beta }_{\mathrm{2}}{R}_{\mathrm{0}}}\right)+\frac{c_{\mathrm{2}}p\alpha \mu n}{{\beta }_{\mathrm{1}}Q}，

c_{\mathrm{0}}=p^{\mathrm{3}}{\alpha }^{\mathrm{3}}{n}^{\mathrm{3}}{\mu }^{\mathrm{2}}gmq(R_{\mathrm{0}}-\mathrm{1})。

通过分析，有以下定理：

定理1 1）当R_{\mathrm{0}}>\mathrm{1}且qgmn>r_{\mathrm{3}}M时，模型（1）有一个正平衡点。

2）当且qgmn>r_{\mathrm{3}}M时，（i）当时，模型（1）无正平衡点；（ii）当c_{\mathrm{1}}>\mathrm{0}时，若c_{\mathrm{1}}^{\mathrm{2}}=\mathrm{4}{c}_{\mathrm{2}}{c}_{\mathrm{0}}，则模型（1）有一个正平衡点，若c_{\mathrm{1}}^{\mathrm{2}}>\mathrm{4}{c}_{\mathrm{2}}{c}_{\mathrm{0}}，则模型（1）有两个正平衡点。

3 平衡点的稳定性

定理2 如果且m=n，则模型（1）的无病平衡点P_{\mathrm{0}}是局部渐近稳定的。

证明 模型（1）在无病平衡点P_{\mathrm{0}}处的Jacobian矩阵为

\mathit{J}(P_{\mathrm{0}})=\left(\begin{array}{cccccc}-\mu & \mathrm{0}& \mathrm{0}& -\frac{{\beta }_{\mathrm{1}}(\mathrm{1}-\tau )\Lambda }{\mu }& -\frac{{\beta }_{\mathrm{1}}(\mathrm{1}-\tau )\Lambda \sigma }{\mu }& \mathrm{0}\\ \mathrm{0}& -\mu & \mathrm{0}& -\frac{{\beta }_{\mathrm{2}}\tau \Lambda }{\mu }& -\frac{{\beta }_{\mathrm{2}}\tau \Lambda \sigma }{\mu }& r_{\mathrm{3}}\\ \mathrm{0}& \mathrm{0}& -(\alpha +\mu )& \frac{{\beta }_{\mathrm{1}}(\mathrm{1}-\tau )\Lambda +{\beta }_{\mathrm{2}}\tau \Lambda }{\mu }& \frac{\sigma [{\beta }_{\mathrm{1}}(\mathrm{1}-\tau )\Lambda +{\beta }_{\mathrm{2}}\tau \Lambda ]}{\mu }& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& p\alpha & -(r_{\mathrm{1}}+\mu )& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& (\mathrm{1}-p)\alpha & \mathrm{0}& -(r_{\mathrm{2}}+\mu +d)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& r_{\mathrm{1}}& r_{\mathrm{2}}& -(\mu +r_{\mathrm{3}})\end{array}\right)。

矩阵\mathit{J}(P_{\mathrm{0}})的特征方程为

{(\lambda +\mu )}^{\mathrm{2}}(\lambda +\mu +r_{\mathrm{3}})(\lambda +m)\left[{\lambda }^{\mathrm{2}}+\lambda (g+m)+\frac{TQ}{m}\left(\frac{\mathrm{1}}{R_{\mathrm{0}}}-\mathrm{1}\right)\right]=\mathrm{0}。（3）

（3）式有四个负特征根{\lambda }_{\mathrm{1}}={\lambda }_{\mathrm{2}}=-\mu，{\lambda }_{\mathrm{3}}=-(\mu +r_{\mathrm{3}})，{\lambda }_{\mathrm{4}}=-m，其余两个特征根满足下面方程

{\lambda }^{\mathrm{2}}+\lambda (g+m)+\frac{TQ}{m}\left(\frac{\mathrm{1}}{R_{\mathrm{0}}}-\mathrm{1}\right)=\mathrm{0}。（4）

其中k_{\mathrm{1}}=g+m，k_{\mathrm{2}}=\frac{TQ}{m}\left(\frac{\mathrm{1}}{R_{\mathrm{0}}}-\mathrm{1}\right)。

当且m=n时，k_{\mathrm{1}}>\mathrm{0}，k_{\mathrm{2}}>\mathrm{0}，则{\lambda }_{\mathrm{5}}{\lambda }_{\mathrm{6}}>\mathrm{0}，，所以{\lambda }_{\mathrm{5}}，{\lambda }_{\mathrm{6}}均小于0，所以方程（3）的所有根都具有负实部。由Hurwitz判据^［17］可得，无病平衡点P_{\mathrm{0}}是局部渐近稳定的。

定理3 如果，则模型（1）的无病平衡点P_{\mathrm{0}}是全局渐近稳定的。

证明 构造Lyapunov函数：

V=E+\frac{{\beta }_{\mathrm{1}}{S}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{S}_{\mathrm{2}}}{m}{I}_{\mathrm{1}}+\frac{({\beta }_{\mathrm{1}}{S}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{S}_{\mathrm{2}})\sigma }{n}{I}_{\mathrm{2}}。

沿着模型（1）的解对V进行求导，可得

\begin{array}{c}\dot{V}=\dot{E}+\frac{{\beta }_{\mathrm{1}}{S}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{S}_{\mathrm{2}}}{m}{\dot{I}}_{\mathrm{1}}+\frac{({\beta }_{\mathrm{1}}{S}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{S}_{\mathrm{2}})\sigma }{n}{\dot{I}}_{\mathrm{2}}=\\ ({\beta }_{\mathrm{1}}{S}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{S}_{\mathrm{2}})(I_{\mathrm{1}}+\sigma I_{\mathrm{2}})-gE+\frac{{\beta }_{\mathrm{1}}{S}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{S}_{\mathrm{2}}}{m}(p\alpha E-m{I}_{\mathrm{1}})+\frac{({\beta }_{\mathrm{1}}{S}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{S}_{\mathrm{2}})\sigma }{n}[(\mathrm{1}-p)\alpha E-n{I}_{\mathrm{2}}]\le \\ \left(\frac{{\beta }_{\mathrm{1}}(\mathrm{1}-\tau )\Lambda }{\mu }+\frac{{\beta }_{\mathrm{2}}\tau \Lambda }{\mu }\right)(I_{\mathrm{1}}+\sigma I_{\mathrm{2}})-gE+\frac{{\beta }_{\mathrm{1}}(\mathrm{1}-\tau )\Lambda +{\beta }_{\mathrm{2}}\tau \Lambda }{\mu m}(p\alpha E-m{I}_{\mathrm{1}})+\\ \frac{[{\beta }_{\mathrm{1}}(\mathrm{1}-\tau )\Lambda +{\beta }_{\mathrm{2}}\tau \Lambda ]\sigma }{\mu n}[(\mathrm{1}-p)\alpha E-n{I}_{\mathrm{2}}]\le \end{array}

\begin{array}{c}T(I_{\mathrm{1}}+\sigma I_{\mathrm{2}})-gE+\frac{Tp\alpha E}{m}-T{I}_{\mathrm{1}}+\frac{T\sigma (\mathrm{1}-p)\alpha E}{n}-T\sigma I_{\mathrm{2}}=\\ -gE+\frac{Tp\alpha E}{m}+\frac{T\sigma (\mathrm{1}-p)\alpha E}{n}=g(R_{\mathrm{0}}-\mathrm{1})E\mathrm{。}\end{array}（5）

当R_{\mathrm{0}}\le \mathrm{1}时，\dot{V}\le \mathrm{0}。因此，由LaSalle的不变性原理^［18］，定理3得证。

定理4 如果R_{\mathrm{0}}>\mathrm{1}且qgmn>r_{\mathrm{3}}M，模型（1）的地方病平衡点P^{\mathrm{*}}是全局渐近稳定的。

证明 构造Lyapunov函数：

\begin{array}{c}L=\left(S_{\mathrm{1}}-S_{\mathrm{1}}^{*}-S_{\mathrm{1}}^{*}\mathrm{l}\mathrm{n}\frac{S_{\mathrm{1}}}{S_{\mathrm{1}}^{\mathrm{*}}}\right)+\left(S_{\mathrm{2}}-S_{\mathrm{2}}^{*}-S_{\mathrm{2}}^{*}\mathrm{l}\mathrm{n}\frac{S_{\mathrm{2}}}{S_{\mathrm{2}}^{\mathrm{*}}}\right)+\left(E-E^{*}-E^{*}\mathrm{l}\mathrm{n}\frac{E}{E^{\mathrm{*}}}\right)+\\ \left(I_{\mathrm{1}}-I_{\mathrm{1}}^{*}-I_{\mathrm{1}}^{*}\mathrm{l}\mathrm{n}\frac{I_{\mathrm{1}}}{I_{\mathrm{1}}^{\mathrm{*}}}\right)+\left(I_{\mathrm{2}}-I_{\mathrm{2}}^{*}-I_{\mathrm{2}}^{*}\mathrm{l}\mathrm{n}\frac{I_{\mathrm{2}}}{I_{\mathrm{2}}^{\mathrm{*}}}\right)+\left(R-R^{*}-R^{*}\mathrm{l}\mathrm{n}\frac{R}{R^{\mathrm{*}}}\right)\mathrm{。}\end{array}（6）

沿着模型（1）的解对L进行求导，可得

\dot{L}=\left(\mathrm{1}-\frac{S_{\mathrm{1}}^{\mathrm{*}}}{S_{\mathrm{1}}}\right){\dot{S}}_{\mathrm{1}}+\left(\mathrm{1}-\frac{S_{\mathrm{2}}^{\mathrm{*}}}{S_{\mathrm{2}}}\right){\dot{S}}_{\mathrm{2}}+\left(\mathrm{1}-\frac{E^{\mathrm{*}}}{E}\right)\dot{E}+\left(\mathrm{1}-\frac{I_{\mathrm{1}}^{\mathrm{*}}}{I_{\mathrm{1}}}\right){\dot{I}}_{\mathrm{1}}+\left(\mathrm{1}-\frac{I_{\mathrm{2}}^{\mathrm{*}}}{I_{\mathrm{2}}}\right){\dot{I}}_{\mathrm{2}}+\left(\mathrm{1}-\frac{R^{\mathrm{*}}}{R}\right)\dot{R}\mathrm{。}

令x_{\mathrm{1}}=\frac{S_{\mathrm{1}}}{S_{\mathrm{1}}^{\mathrm{*}}},x_{\mathrm{2}}=\frac{S_{\mathrm{2}}}{S_{\mathrm{2}}^{\mathrm{*}}},x_{\mathrm{3}}=\frac{E}{E^{\mathrm{*}}},x_{\mathrm{4}}=\frac{I_{\mathrm{1}}}{I_{\mathrm{1}}^{\mathrm{*}}},x_{\mathrm{5}}=\frac{I_{\mathrm{2}}}{I_{\mathrm{2}}^{\mathrm{*}}},x_{\mathrm{6}}=\frac{R}{R^{\mathrm{*}}},化简上述等式

\begin{array}{c}\dot{L}=-(\mathrm{1}-\tau )\Lambda x_{\mathrm{1}}{\left(\mathrm{1}-\frac{\mathrm{1}}{x_{\mathrm{1}}}\right)}^{\mathrm{2}}-\tau \Lambda x_{\mathrm{2}}{\left(\mathrm{1}-\frac{\mathrm{1}}{x_{\mathrm{2}}}\right)}^{\mathrm{2}}-{\beta }_{\mathrm{1}}{I}_{\mathrm{1}}^{\mathrm{*}}{S}_{\mathrm{1}}\left(x_{\mathrm{4}}-\frac{x_{\mathrm{4}}}{x_{\mathrm{1}}}+\frac{\mathrm{1}}{x_{\mathrm{1}}}-\mathrm{1}\right)-\\ \sigma {\beta }_{\mathrm{1}}{I}_{\mathrm{2}}^{\mathrm{*}}{S}_{\mathrm{1}}\left(x_{\mathrm{5}}-\frac{x_{\mathrm{5}}}{x_{\mathrm{1}}}+\frac{\mathrm{1}}{x_{\mathrm{1}}}-\mathrm{1}\right)-{\beta }_{\mathrm{2}}{I}_{\mathrm{1}}^{\mathrm{*}}{S}_{\mathrm{2}}\left(x_{\mathrm{4}}-\frac{x_{\mathrm{4}}}{x_{\mathrm{2}}}+\frac{\mathrm{1}}{x_{\mathrm{2}}}-\mathrm{1}\right)-\\ \sigma {\beta }_{\mathrm{2}}{I}_{\mathrm{2}}^{\mathrm{*}}{S}_{\mathrm{2}}\left(x_{\mathrm{5}}-\frac{x_{\mathrm{5}}}{x_{\mathrm{2}}}+\frac{\mathrm{1}}{x_{\mathrm{2}}}-\mathrm{1}\right)-r_{\mathrm{3}}{R}^{\mathrm{*}}\left(-x_{\mathrm{6}}+x_{\mathrm{2}}+\frac{x_{\mathrm{6}}}{x_{\mathrm{2}}}-\mathrm{1}\right)-\\ \left(-x_{\mathrm{1}}{x}_{\mathrm{4}}+\frac{x_{\mathrm{1}}{x}_{\mathrm{4}}}{x_{\mathrm{3}}}+x_{\mathrm{3}}-\mathrm{1}\right){\beta }_{\mathrm{1}}{S}_{\mathrm{1}}^{\mathrm{*}}{I}_{\mathrm{1}}^{\mathrm{*}}-\left(-x_{\mathrm{1}}{x}_{\mathrm{5}}+\frac{x_{\mathrm{1}}{x}_{\mathrm{5}}}{x_{\mathrm{3}}}+x_{\mathrm{3}}-\mathrm{1}\right)\sigma {\beta }_{\mathrm{1}}{S}_{\mathrm{1}}^{\mathrm{*}}{I}_{\mathrm{2}}^{\mathrm{*}}-\\ \left(-x_{\mathrm{2}}{x}_{\mathrm{4}}+\frac{x_{\mathrm{2}}{x}_{\mathrm{4}}}{x_{\mathrm{3}}}+x_{\mathrm{3}}-\mathrm{1}\right){\beta }_{\mathrm{2}}{S}_{\mathrm{2}}^{\mathrm{*}}{I}_{\mathrm{1}}^{\mathrm{*}}-\left(-x_{\mathrm{2}}{x}_{\mathrm{5}}+\frac{x_{\mathrm{2}}{x}_{\mathrm{5}}}{x_{\mathrm{3}}}+x_{\mathrm{3}}-\mathrm{1}\right)\sigma {\beta }_{\mathrm{2}}{S}_{\mathrm{2}}^{\mathrm{*}}{I}_{\mathrm{2}}^{\mathrm{*}}-\\ p\alpha E^{\mathrm{*}}\left(-x_{\mathrm{3}}+x_{\mathrm{4}}+\frac{x_{\mathrm{3}}}{x_{\mathrm{4}}}-\mathrm{1}\right)-(\mathrm{1}-p)\alpha E^{\mathrm{*}}\left(-x_{\mathrm{3}}+x_{\mathrm{5}}+\frac{x_{\mathrm{3}}}{x_{\mathrm{5}}}-\mathrm{1}\right)-\\ r_{\mathrm{1}}{I}_{\mathrm{1}}^{\mathrm{*}}\left(-x_{\mathrm{4}}+x_{\mathrm{6}}+\frac{x_{\mathrm{4}}}{x_{\mathrm{6}}}-\mathrm{1}\right)-r_{\mathrm{2}}{I}_{\mathrm{2}}^{\mathrm{*}}\left(-x_{\mathrm{5}}+x_{\mathrm{6}}+\frac{x_{\mathrm{5}}}{x_{\mathrm{6}}}-\mathrm{1}\right)\mathrm{。}\end{array}（7）

因此，当且仅当x_{\mathrm{1}}=\mathrm{1},x_{\mathrm{2}}=\mathrm{1},x_{\mathrm{3}}=\mathrm{1},x_{\mathrm{4}}=\mathrm{1},x_{\mathrm{5}}=\mathrm{1},x_{\mathrm{6}}=\mathrm{1},则\dot{L}\le \mathrm{0}。所以，由LaSalle的不变性原理^［18］，定理4得证。

定理5 当m=n时，模型（1）在{\beta }_{\mathrm{1}}={\beta }_{\mathrm{1}}^{\mathrm{*}}（即在R_{\mathrm{0}}=\mathrm{1}）处表现出后向分岔。

证明 设x=(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}},x_{\mathrm{4}},x_{\mathrm{5}},x_{\mathrm{6}}{)}^{\mathrm{{\rm T}}}=(S_{\mathrm{1}},S_{\mathrm{2}},E,I_{\mathrm{1}},I_{\mathrm{2}}{,R)}^{\mathrm{{\rm T}}}，可以将模型（1）重写为\frac{\mathrm{d}x}{\mathrm{d}t}=f(x)，其中f(x)=(f(x_{\mathrm{1}}),f(x_{\mathrm{2}}),f(x_{\mathrm{3}}),f(x_{\mathrm{4}}),f(x_{\mathrm{5}}),f(x_{\mathrm{6}}{))}^{\mathrm{{\rm T}}}，即

\left\{\begin{array}{l}\frac{\mathrm{d}{x}_{\mathrm{1}}}{\mathrm{d}t}=f_{\mathrm{1}}=(\mathrm{1}-\tau )\Lambda -{\beta }_{\mathrm{1}}{x}_{\mathrm{1}}(x_{\mathrm{4}}+\sigma x_{\mathrm{5}})-\mu x_{\mathrm{1}},\\ \frac{\mathrm{d}{x}_{\mathrm{2}}}{\mathrm{d}t}=f_{\mathrm{2}}=\tau \Lambda -{\beta }_{\mathrm{2}}{x}_{\mathrm{2}}(x_{\mathrm{4}}+\sigma x_{\mathrm{5}})-\mu x_{\mathrm{2}}+r_{\mathrm{3}}{x}_{\mathrm{6}},\\ \frac{\mathrm{d}{x}_{\mathrm{3}}}{\mathrm{d}t}=f_{\mathrm{3}}=({\beta }_{\mathrm{1}}{x}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{x}_{\mathrm{2}})(x_{\mathrm{4}}+\sigma x_{\mathrm{5}})-\mu x_{\mathrm{3}}-\alpha x_{\mathrm{3}},\\ \frac{\mathrm{d}{x}_{\mathrm{4}}}{\mathrm{d}t}=f_{\mathrm{4}}=p\alpha x_{\mathrm{3}}-r_{\mathrm{1}}{x}_{\mathrm{4}}-\mu x_{\mathrm{4}},\\ \frac{\mathrm{d}{x}_{\mathrm{5}}}{\mathrm{d}t}=f_{\mathrm{5}}=(\mathrm{1}-p)\alpha x_{\mathrm{3}}-r_{\mathrm{2}}{x}_{\mathrm{5}}-\mu x_{\mathrm{5}}-d{x}_{\mathrm{5}},\\ \frac{\mathrm{d}{x}_{\mathrm{6}}}{\mathrm{d}t}=f_{\mathrm{6}}=r_{\mathrm{1}}{x}_{\mathrm{4}}+r_{\mathrm{2}}{x}_{\mathrm{5}}-\mu x_{\mathrm{6}}-r_{\mathrm{3}}{x}_{\mathrm{6}}\mathrm{。}\end{array}\right.（8）

设{\beta }_{\mathrm{1}}={\beta }_{\mathrm{1}}^{\mathrm{*}}=\frac{\mu (\mu +\alpha )(r_{\mathrm{1}}+\mu )(r_{\mathrm{2}}+\mu +d)-{\beta }_{\mathrm{2}}\tau \Lambda [p\alpha n+\sigma (\mathrm{1}-p)\alpha m]}{(\mathrm{1}-\tau )\Lambda [p\alpha n+\sigma (\mathrm{1}-p)\alpha m]}。当m=n且{\beta }_{\mathrm{1}}={\beta }_{\mathrm{1}}^{\mathrm{*}}时，模型（8）在无病平衡点P_{\mathrm{0}}处的Jacobian矩阵为\mathit{J}(P_{\mathrm{0}})，其对应的特征多项式为\lambda {(\lambda +\mu )}^{\mathrm{2}}(\lambda +\mu +r_{\mathrm{3}})(\lambda +m)[(\lambda +(g+m)]=\mathrm{0}\mathrm{。}特征值为{\lambda }_{\mathrm{1}}={\lambda }_{\mathrm{2}}=-\mu，{\lambda }_{\mathrm{3}}=-(\mu +r_{\mathrm{3}})，{\lambda }_{\mathrm{4}}=-m，{\lambda }_{\mathrm{5}}=-(g+m)，{\lambda }_{\mathrm{6}}=\mathrm{0}。因此与零特征值相关的右特征向量为\mathit{\omega }=({\omega }_{\mathrm{1}},{\omega }_{\mathrm{2}},{\omega }_{\mathrm{3}},{\omega }_{\mathrm{4}},{\omega }_{\mathrm{5}},{\omega }_{\mathrm{6}}{)}^{\mathrm{{\rm T}}}，左特征向量为\mathit{v}=(v_{\mathrm{1}},v_{\mathrm{2}},v_{\mathrm{3}},v_{\mathrm{4}},v_{\mathrm{5}},v_{\mathrm{6}})，其中{\omega }_{\mathrm{1}}=\frac{mn(\mathrm{1}-\tau )\Lambda }{{\beta }_{\mathrm{1}}Q{\omega }_{\mathrm{3}}+\mu mn}，{\omega }_{\mathrm{2}}=\frac{qmn\tau \Lambda +r_{\mathrm{3}}[r_{\mathrm{1}}p\alpha n+r_{\mathrm{2}}(\mathrm{1}-p)\alpha m]{\omega }_{\mathrm{3}}}{q{\beta }_{\mathrm{2}}Q{\omega }_{\mathrm{3}}+\mu mnq}，{\omega }_{\mathrm{3}}={\omega }_{\mathrm{3}}>\mathrm{0}，{\omega }_{\mathrm{4}}=\frac{p\alpha {\omega }_{\mathrm{3}}}{m}，{\omega }_{\mathrm{5}}=\frac{(\mathrm{1}-p)\alpha {\omega }_{\mathrm{3}}}{n}，{\omega }_{\mathrm{6}}=\frac{[r_{\mathrm{1}}p\alpha n+r_{\mathrm{2}}(\mathrm{1}-p)\alpha m]{\omega }_{\mathrm{3}}}{mnq}，v_{\mathrm{1}}=v_{\mathrm{2}}=v_{\mathrm{6}}=\mathrm{0}，v_{\mathrm{3}}=v_{\mathrm{3}}>\mathrm{0}，v_{\mathrm{4}}=\frac{[{\beta }_{\mathrm{1}}(\mathrm{1}-\tau )\Lambda +{\beta }_{\mathrm{2}}\tau \Lambda ]v_{\mathrm{3}}}{\mu m}，v_{\mathrm{5}}=\frac{\sigma [{\beta }_{\mathrm{1}}(\mathrm{1}-\tau )\Lambda +{\beta }_{\mathrm{2}}\tau \Lambda ]v_{\mathrm{3}}}{\mu n}，并且\mathit{\omega }和\mathit{v}满足\mathit{v}\cdot \mathit{\omega }=\mathrm{1}。基于中心流形理论^［19］，需要计算a={\displaystyle \sum _{k,i,j=\mathrm{1}}^{\mathrm{6}}}{v}_k{\omega }_i{\omega }_j\frac{{\partial }^{\mathrm{2}}{f}_k}{\partial x_i\partial x_j}(P_{\mathrm{0}},{\beta }_{\mathrm{1}}^{\mathrm{*}})，b={\displaystyle \sum _{k,i=\mathrm{1}}^{\mathrm{6}}}{v}_k{\omega }_i\frac{{\partial }^{\mathrm{2}}{f}_k}{\partial x_i\partial \varphi }(P_{\mathrm{0}},{\beta }_{\mathrm{1}}^{\mathrm{*}})。因此，

\begin{array}{c}a=v_{\mathrm{3}}{\displaystyle \sum _{i,j=\mathrm{1}}^{\mathrm{6}}}{\omega }_i{\omega }_j\frac{{\partial }^{\mathrm{2}}{f}_{\mathrm{3}}}{\partial x_i\partial x_j}+v_{\mathrm{4}}{\displaystyle \sum _{i,j=\mathrm{1}}^{\mathrm{6}}}{\omega }_i{\omega }_j\frac{{\partial }^{\mathrm{2}}{f}_{\mathrm{4}}}{\partial x_i\partial x_j}+v_{\mathrm{5}}{\displaystyle \sum _{i,j=\mathrm{1}}^{\mathrm{6}}}{\omega }_i{\omega }_j\frac{{\partial }^{\mathrm{2}}{f}_{\mathrm{5}}}{\partial x_i\partial x_j}=\\ \frac{\mathrm{2}{v}_{\mathrm{3}}{\omega }_{\mathrm{3}}Q}{mn}\left[\frac{mn(\mathrm{1}-\tau )\Lambda {\beta }_{\mathrm{1}}}{{\beta }_{\mathrm{1}}Q{\omega }_{\mathrm{3}}+\mu mn}+\frac{qmn\tau \Lambda +{\omega }_{\mathrm{3}}{r}_{\mathrm{3}}(r_{\mathrm{1}}p\alpha n+r_{\mathrm{2}}(\mathrm{1}-p)\alpha m)}{qQ{\omega }_{\mathrm{3}}+\mu mnq/{\beta }_{\mathrm{2}}}\right]>\mathrm{0},\end{array}

b=\frac{(\mathrm{1}-\tau )\Lambda Q{\omega }_{\mathrm{3}}{v}_{\mathrm{3}}}{\mu mn}>\mathrm{0}\mathrm{。}

通过计算可知a>\mathrm{0}和b>\mathrm{0}。因此，在R_{\mathrm{0}}=\mathrm{1}时，a>\mathrm{0}且b>\mathrm{0}，模型（1）经历后向分岔。

4 数值模拟和敏感性分析

下面通过数值模拟，验证上述传染病传播模型得到的理论分析结果。选择初始数值\left(S_{\mathrm{1}},S_{\mathrm{2}},E,I_{\mathrm{1}},I_{\mathrm{2}},R\right)=\left(\mathrm{0.296}\mathrm{4,0.697}\mathrm{6,0.000}\mathrm{049}\mathrm{7,0.000}\mathrm{381}\mathrm{3,0},\mathrm{0.005}\mathrm{569}\right)。其他参数值如表1所示。

为了说明平衡点的局部稳定性，我们选择了不同的初始条件来模拟模型。图1证明了当基本再生数R_{\mathrm{0}}>\mathrm{1}时，地方病平衡点P^{\mathrm{*}}的局部渐近稳定性。图2模拟了无病平衡点和地方病平衡点的全局稳定性。

从前面的分析中，我们可以看到，基本再生数决定了传染病模型的定性行为。因此，分析对R_{\mathrm{0}}值影响较大的参数是很重要的。为此，我们用公式{\epsilon }_{R_{\mathrm{0}}}^p=\frac{\partial R_{\mathrm{0}}}{\partial p}\frac{p}{R_{\mathrm{0}}}来研究表1中给出的各参数对R_{\mathrm{0}}的影响^［20］。

从表1中可以看出，参数\mathrm{\Lambda },{\beta }_{\mathrm{1}},{\beta }_{\mathrm{2}},\sigma ,p和\alpha与R_{\mathrm{0}}呈正相关，这些参数每增长或减少10%，R_{\mathrm{0}}分别增加或减少10%，4.6%，5.3%，4.2%，2.8%和0.000 206 9%。随着这些参数值的增加，病毒的感染程度将变得更加严重。而\tau ,r_{\mathrm{1}},r_{\mathrm{2}},\mu和d每增加10%，R_{\mathrm{0}}则分别减少5.3%，5.7%，4.4%，10%，和0.056%。r_{\mathrm{3}}对R_{\mathrm{0}}没有影响。

5 模型的最优控制

传染病的最优控制是一个重要而复杂的问题。基于上述的敏感性分析，我们扩展了模型（1），将高危易感人群的疫苗接种和增强意识计划（如社交平台、广播、电视和其他社交媒体上的媒体活动）包括在内。假设高危易感人群的疫苗接种成功率为u_{\mathrm{1}}。此外，增强意识方案可以增加具有保护意识的感染者的数量。增强意识项目的比率为u_{\mathrm{2}}。控制模型由以下微分方程组表示：

\left\{\begin{array}{l}\frac{\mathrm{d}{S}_{\mathrm{1}}}{\mathrm{d}t}=(\mathrm{1}-\tau )\Lambda -{\beta }_{\mathrm{1}}{S}_{\mathrm{1}}(I_{\mathrm{1}}+\sigma I_{\mathrm{2}})-\mu S_{\mathrm{1}}-u_{\mathrm{1}}{S}_{\mathrm{1}},\\ \frac{\mathrm{d}{S}_{\mathrm{2}}}{\mathrm{d}t}=\tau \Lambda -{\beta }_{\mathrm{2}}{S}_{\mathrm{2}}(I_{\mathrm{1}}+\sigma I_{\mathrm{2}})-\mu S_{\mathrm{2}}+r_{\mathrm{3}}R+u_{\mathrm{1}}{S}_{\mathrm{1}},\\ \frac{\mathrm{d}E}{\mathrm{d}t}=({\beta }_{\mathrm{1}}{S}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{S}_{\mathrm{2}})(I_{\mathrm{1}}+\sigma I_{\mathrm{2}})-\mu E-\alpha E,\\ \frac{\mathrm{d}{I}_{\mathrm{1}}}{\mathrm{d}t}=p\alpha E-r_{\mathrm{1}}{I}_{\mathrm{1}}-\mu I_{\mathrm{1}}-u_{\mathrm{2}}{I}_{\mathrm{1}},\\ \frac{\mathrm{d}{I}_{\mathrm{2}}}{\mathrm{d}t}=(\mathrm{1}-p)\alpha E-r_{\mathrm{2}}{I}_{\mathrm{2}}-\mu I_{\mathrm{2}}-d{I}_{\mathrm{2}}+u_{\mathrm{2}}{I}_{\mathrm{1}},\\ \frac{\mathrm{d}R}{\mathrm{d}t}=r_{\mathrm{1}}{I}_{\mathrm{1}}+r_{\mathrm{2}}{I}_{\mathrm{2}}-\mu R-r_{\mathrm{3}}R\mathrm{。}\end{array}\right.（9）

定义目标函数

J(u_{\mathrm{1}},u_{\mathrm{2}})=\left(A_{\mathrm{1}}E+A_{\mathrm{2}}{I}_{\mathrm{1}}+A_{\mathrm{3}}{I}_{\mathrm{2}}+\frac{B_{\mathrm{1}}}{\mathrm{2}}{u}_{\mathrm{1}}^{\mathrm{2}}+\frac{B_{\mathrm{2}}}{\mathrm{2}}{u}_{\mathrm{2}}^{\mathrm{2}}\right)\mathrm{d}t,（10）

其中T是实施控制策略的最终时间。系数A_{\mathrm{1}},A_{\mathrm{2}}和A_{\mathrm{3}}分别是E,I_{\mathrm{1}}和I_{\mathrm{2}}正权重常数，B_{\mathrm{1}}和B_{\mathrm{2}}表示收益和成本的权重常数。定义控制集U=\left\{u=(u_{\mathrm{1}}(t),u_{\mathrm{2}}(t):\mathrm{0}\le u_i(t)\le \mathrm{1},t\in [\mathrm{0},T],i=\mathrm{1,2}\right\}是Lebesgue可测控制量。目标是找到一个最优控制对(u_{\mathrm{1}}^{\mathrm{*}},u_{\mathrm{2}}^{\mathrm{*}})，使得

J(u_{\mathrm{1}}^{\mathrm{*}},u_{\mathrm{2}}^{\mathrm{*}})=\mathrm{m}\mathrm{i}\mathrm{n}\left\{J\left(u_{\mathrm{1}}(t),u_{\mathrm{2}}(t)\right)\left|\left(u_{\mathrm{1}}(t),u_{\mathrm{2}}(t)\right)\right.\in U\right\}。

由模型（9）右端的有界性及式（10）的被积函数关于u_i(t)(i=\mathrm{1,2})的凸性^［22-23］，得到最优控制对(u_{\mathrm{1}}^{\mathrm{*}},u_{\mathrm{2}}^{\mathrm{*}})的存在性。

使用Pontryagin的最大原理^［23］导出最优控制函数的必要条件，将最优控制问题转化为哈密顿函数最小化问题。定义如下Hamilton函数

\begin{array}{c}H=A_{\mathrm{1}}E+A_{\mathrm{2}}{I}_{\mathrm{1}}+A_{\mathrm{3}}{I}_{\mathrm{2}}+\frac{B_{\mathrm{1}}}{\mathrm{2}}{u}_{\mathrm{1}}^{\mathrm{2}}+\frac{B_{\mathrm{2}}}{\mathrm{2}}{u}_{\mathrm{2}}^{\mathrm{2}}+{\lambda }_{\mathrm{1}}[(\mathrm{1}-\tau )\mathrm{\Lambda }-{\beta }_{\mathrm{1}}{S}_{\mathrm{1}}(I_{\mathrm{1}}+\sigma I_{\mathrm{2}})-\mu S_{\mathrm{1}}-u_{\mathrm{1}}{S}_{\mathrm{1}}]+\\ {\lambda }_{\mathrm{2}}[\tau \Lambda -{\beta }_{\mathrm{2}}{S}_{\mathrm{2}}(I_{\mathrm{1}}+\sigma I_{\mathrm{2}})-\mu S_{\mathrm{2}}+r_{\mathrm{3}}R+u_{\mathrm{1}}{S}_{\mathrm{1}}]+{\lambda }_{\mathrm{3}}[({\beta }_{\mathrm{1}}{S}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{S}_{\mathrm{2}})(I_{\mathrm{1}}+\sigma I_{\mathrm{2}})-\mu E-\alpha E]+\\ {\lambda }_{\mathrm{4}}[p\alpha E-r_{\mathrm{1}}{I}_{\mathrm{1}}-\mu I_{\mathrm{1}}-u_{\mathrm{2}}{I}_{\mathrm{1}}]+{\lambda }_{\mathrm{5}}[(\mathrm{1}-p)\alpha E-r_{\mathrm{2}}{I}_{\mathrm{2}}-\mu I_{\mathrm{2}}-d{I}_{\mathrm{2}}+u_{\mathrm{2}}{I}_{\mathrm{1}}]+{\lambda }_{\mathrm{6}}[r_{\mathrm{1}}{I}_{\mathrm{1}}+r_{\mathrm{2}}{I}_{\mathrm{2}}-\mu R-r_{\mathrm{3}}R],\end{array}（11）

其中\lambda =({\lambda }_{\mathrm{1}},{\lambda }_{\mathrm{2}},{\lambda }_{\mathrm{3}},{\lambda }_{\mathrm{4}},{\lambda }_{\mathrm{5}},{\lambda }_{\mathrm{6}}{)}^{\mathrm{{\rm T}}}，{\lambda }_i是与S_{\mathrm{1}},S_{\mathrm{2}},E,I_{\mathrm{1}},I_{\mathrm{2}}和R相关的伴随变量。

定理6 设u_{\mathrm{1}}^{\mathrm{*}}和u_{\mathrm{2}}^{\mathrm{*}}是（9）式中给出的控制问题的控制函数，S_{\mathrm{1}}^{\mathrm{*}},S_{\mathrm{2}}^{\mathrm{*}},E^{\mathrm{*}},I_{\mathrm{1}}^{\mathrm{*}},I_{\mathrm{2}}^{\mathrm{*}}和R^{\mathrm{*}}是状态变量的解。则存在满足以下方程的伴随变量{\lambda }_{\mathrm{1}},{\lambda }_{\mathrm{2}},{\lambda }_{\mathrm{3}},{\lambda }_{\mathrm{4}},{\lambda }_{\mathrm{5}}和{\lambda }_{\mathrm{6}}。

\left\{\begin{array}{l}\frac{\mathrm{d}{\lambda }_{\mathrm{1}}}{\mathrm{d}t}={\beta }_{\mathrm{1}}(I_{\mathrm{1}}+\sigma I_{\mathrm{2}})({\lambda }_{\mathrm{1}}-{\lambda }_{\mathrm{3}})+u_{\mathrm{1}}({\lambda }_{\mathrm{1}}-{\lambda }_{\mathrm{2}})+\mu {\lambda }_{\mathrm{1}},\\ \frac{\mathrm{d}{\lambda }_{\mathrm{2}}}{\mathrm{d}t}={\beta }_{\mathrm{2}}(I_{\mathrm{1}}+\sigma I_{\mathrm{2}})({\lambda }_{\mathrm{2}}-{\lambda }_{\mathrm{3}})+\mu {\lambda }_{\mathrm{2}},\\ \frac{\mathrm{d}{\lambda }_{\mathrm{3}}}{\mathrm{d}t}={\lambda }_{\mathrm{3}}(\mu +\alpha )-{\lambda }_{\mathrm{4}}p\alpha -{\lambda }_{\mathrm{5}}(\mathrm{1}-p)\alpha -A_{\mathrm{1}},\\ \frac{\mathrm{d}{\lambda }_{\mathrm{4}}}{\mathrm{d}t}=({\lambda }_{\mathrm{1}}-{\lambda }_{\mathrm{3}}){\beta }_{\mathrm{1}}{S}_{\mathrm{1}}+({\lambda }_{\mathrm{2}}-{\lambda }_{\mathrm{3}}){\beta }_{\mathrm{2}}{S}_{\mathrm{2}}+{\lambda }_{\mathrm{4}}(r_{\mathrm{1}}+\mu +u_{\mathrm{2}})-{\lambda }_{\mathrm{5}}{u}_{\mathrm{2}}-{\lambda }_{\mathrm{6}}{r}_{\mathrm{1}}-A_{\mathrm{2}},\\ \frac{\mathrm{d}{\lambda }_{\mathrm{5}}}{\mathrm{d}t}=({\lambda }_{\mathrm{1}}-{\lambda }_{\mathrm{3}})\sigma {\beta }_{\mathrm{1}}{S}_{\mathrm{1}}+({\lambda }_{\mathrm{2}}-{\lambda }_{\mathrm{3}})\sigma {\beta }_{\mathrm{2}}{S}_{\mathrm{2}}+{\lambda }_{\mathrm{5}}(r_{\mathrm{2}}+\mu +d)-{\lambda }_{\mathrm{6}}{r}_{\mathrm{2}}-A_{\mathrm{3}},\\ \frac{\mathrm{d}{\lambda }_{\mathrm{6}}}{\mathrm{d}t}=-{\lambda }_{\mathrm{2}}{r}_{\mathrm{3}}+{\lambda }_{\mathrm{6}}(\mu +r_{\mathrm{3}}),\end{array}\right.（12）

具有边界条件{\lambda }_i(T)=\mathrm{0},i=\mathrm{1,2},\mathrm{3,4},\mathrm{5,6}。此外，最优控制策略u_{\mathrm{1}}^{\mathrm{*}}和u_{\mathrm{2}}^{\mathrm{*}}为

u_{\mathrm{1}}^{\mathrm{*}}=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\mathrm{0},\mathrm{m}\mathrm{i}\mathrm{n}\left\{\mathrm{1},\frac{({\lambda }_{\mathrm{1}}-{\lambda }_{\mathrm{2}})S_{\mathrm{1}}}{B_{\mathrm{1}}}\right\}\right\},u_{\mathrm{2}}^{\mathrm{*}}=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\mathrm{0},\mathrm{m}\mathrm{i}\mathrm{n}\left\{\mathrm{1},\frac{({\lambda }_{\mathrm{4}}-{\lambda }_{\mathrm{5}})I_{\mathrm{1}}}{B_{\mathrm{2}}}\right\}\right\}\mathrm{。}（13）

证明 利用Pontryagin的最大原理，对每个状态的哈密顿量进行导数运算，得到以下方程

\frac{\mathrm{d}{\lambda }_i}{\mathrm{d}t}=-\frac{\partial H}{\partial x_i},（14）

其中{\lambda }_i，i=\mathrm{1,2},\mathrm{3,4},\mathrm{5,6}分别表示S_{\mathrm{1}},S_{\mathrm{2}},E,I_{\mathrm{1}},I_{\mathrm{2}},R。

最优控制u_{\mathrm{1}}^{\mathrm{*}}和u_{\mathrm{2}}^{\mathrm{*}}满足\frac{\partial H}{\partial u_i^{\mathrm{*}}}=\mathrm{0}，i=\mathrm{1,2}。因此

u_{\mathrm{1}}^{\mathrm{*}}=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\mathrm{0},\mathrm{m}\mathrm{i}\mathrm{n}\left\{\mathrm{1},\frac{({\lambda }_{\mathrm{1}}-{\lambda }_{\mathrm{2}})S_{\mathrm{1}}}{B_{\mathrm{1}}}\right\}\right\},u_{\mathrm{2}}^{\mathrm{*}}=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\mathrm{0},\mathrm{m}\mathrm{i}\mathrm{n}\left\{\mathrm{1},\frac{({\lambda }_{\mathrm{4}}-{\lambda }_{\mathrm{5}})I_{\mathrm{1}}}{B_{\mathrm{2}}}\right\}\right\}\mathrm{。}

在Matlab^［24］中利用前后扫描法^［25］对最优系统进行了数值仿真。当实施（高危易感人群的疫苗接种和增强意识计划）组合策略时，从图3（a），图3（c），图3（d），图3（e）和图3（f）中可以观察到，高危易感人群，无意识感染人群，有意识感染人群，潜伏人群和恢复人群的数量都减少了。从图3（b）中观察到，低危易感人群数量逐渐增加到0.93。这意味着，社会上低危易感人群居多。控制结果表明，通过对高危易感人群接种疫苗以及提高感染者的意识，可以有效降低感染率，从而有效控制传染病的传播。

6 结论

在本文中，我们提出了一个S₁S₂EI₁I₂R模型，该模型考虑了高、低风险的易感者和有、无意识的感染者对传染病传播动力学的影响。首先，对这个数学模型进行了定性分析，求解了基本再生数R_{\mathrm{0}}，并进行了稳定性分析和分岔分析。然后，对模型进行了数值模拟和灵敏度分析，讨论了参数对基本再生数的影响。最后，为了更好地控制疫情，最大限度地减少感染者的数量，本文以高危易感人群的疫苗接种率以及提高感染者的意识计划为控制变量，建立了最优控制系统。数值模拟验证了理论分析的正确性。因此，为了减少和消除传染病的感染人数，应大力支持高危易感者接种疫苗和增强感染者的意识计划。

参考文献

[1] BEACH B, CLAY K, SAAVEDRA M. The 1918 Influenza Pandemic and Its Lessons for COVID-19[J]. J Econ Lit, 2022, 60(1): 41-84. DOI: 10.1257/jel.20201641 .

[2] WANG W D, RUAN S G. Simulating the SARS Outbreak in Beijing with Limited Data[J]. J Theor Biol, 2004, 227(3): 369-379. DOI: 10.1016/j.jtbi.2003.11.014 .

[3] YAN X F, ZOU Y. Optimal and Sub-optimal Quarantine and Isolation Control in SARS Epidemics[J]. Math Comput Model, 2008, 47(1): 235-245. DOI: 10.1016/j.mcm.2007.04.003 .

[4] RONG X M, YANG L, CHU H D, et al. Effect of Delay in Diagnosis on Transmission of COVID-19[J]. Math Biosci Eng, 2020, 17(3): 2725-2740. DOI: 10.3934/mbe.2020149 .

[5] TANG B, WANG X, LI Q, et al. Estimation of the Transmission Risk of the 2019-nCoV and Its Implication for Public Health Interventions[J]. J Clin Med, 2020, 9(2): 462. DOI: 10.3390/jcm9020462 .

[6] LOUBET P, LOULERGUE P, GALTIER F, et al. Seasonal Influenza Vaccination of High-risk Adults[J]. Expert Rev Vaccines, 2016, 15(12): 1507-1518. DOI: 10.1080/14760584.2016.1188696 .

[7] TYLICKI L, PUCHALSKA-REGLIŃSKA E, TYLICKI P, et al. Predictors of Mortality in Hemodialyzed Patients after SARS-COV-2 Infection[J]. J Clin Med, 2022, 11(2): 285. DOI: 10.3390/jcm11020285 .

[8] NDII M Z, ADI Y A. Understanding the Effects of Individual Awareness and Vector Controls on Malaria Transmission Dynamics Using Multiple Optimal Control[J]. Chaos Soliton Fract, 2021, 153:(1) 111476. DOI: 10.1016/j.chaos.2021.111476 .

[9] YUAN Y R, LI N. Optimal Control and Cost-effectiveness Analysis for a COVID-19 Model with Individual Protection Awareness[J]. Physica A, 2022, 603: 127804. DOI: 10.1016/j.physa.2022.127804 .

[10] KAR T K, NANDI S K, JANA S, et al. Stability and Bifurcation Analysis of an Epidemic Model with the Effect of Media[J]. Chao Soliton Fract, 2019, 120: 188-199. DOI: 10.1016/j.chaos.2019.01.025 .

[11] MUSA S S, QURESHI S, ZHAO S, et al. Mathematical Modeling of COVID-19 Epidemic with Effect of Awareness Programs[J]. Infect Dis Model, 2021, 6: 448-460. DOI: 10.1016/j.idm.2021.01.012 .

[12] ALDILA D. Analyzing the Impact of the Media Campaign and Rapid Testing for COVID-19 as an Optimal Control Problem in East Java, Indonesia[J]. Chao Soliton Fract, 2020, 141: 110364. DOI: 10.1016/j.chaos.2020.110364 .

[13] KUMAR A, DUBEY U S, DUBEY B. The Impact of Social Media Advertisements and Treatments on the Dynamics of Infectious Diseases with Optimal Control Strategies[J]. Math Comput Simul, 2024, 219: 50-86. DOI: 10.1016/j.matcom.2023.12.015 .

[14] BASIR F AL, RAJAK B, RAHMAN B, et al. Hopf Bifurcation Analysis and Optimal Control of an Infectious Disease with Awareness Campaign and Treatment[J]. Axioms, 2023, 12(6): 608. DOI: 10.3390/axioms12060608 .

[15] IBRAHIM M M, AHMAD KAMRAN M, NAEEM MANNAN M M, et al. Impact of Awareness to Control Malaria Disease: A Mathematical Modeling Approach[J]. Complexity, 2020, 2020: 8657410. DOI: 10.1155/2020/8657410 .

[16] VAN DEN DRIESSCHE P, WATMOUGH J. Reproduction Numbers and Sub-threshold Endemic Equilibria for Compartmental Models of Disease Transmission[J]. Math Biosci, 2002, 180(1/2): 29-48. DOI: 10.1016/s0025-5564(02)00108-6 .

[17] 马知恩, 周义仓. 常微分方程定性与稳定性方法[M]. 北京: 科学出版社, 2001.

[18] MA Z E, ZHOU Y C. Qualitative and Stability Methods of Ordinary Differential Equations[M]. Beijing: Science Press, 2001.

[19] LASALLE J P. The Stability of Dynamical Systems[M]. Philadelphia: Society for Industrial and Applied Mathematics, 1976.

[20] CASTILLO-CHAVEZ C, SONG B J. Dynamical Models of Tuberculosis and Their Applications[J]. Math Biosci Eng, 2004, 1(2): 361-404. DOI: 10.3934/mbe.2004.1.361 .

[21] MAKINDE O D, OKOSUN K O. Impact of Chemo-therapy on Optimal Control of Malaria Disease with Infected Immigrants[J]. Biosystems, 2011, 104(1): 32-41. DOI: 10.1016/j.biosystems.2010.12.010 .

[22] IBRAHIM A, HUMPHRIES U W, KHAN A, et al. COVID-19 Model with High- and Low-risk Susceptible Population Incorporating the Effect of Vaccines[J]. Vaccines, 2022, 11(1): 3. DOI: 10.3390/vaccines11010003 .

[23] ALEXANDER M E, MOGHADAS S M, RÖST G, et al. A Delay Differential Model for Pandemic Influenza with Antiviral Treatment[J]. Bull Math Biol, 2008, 70(2): 382-397. DOI: 10.1007/s11538-007-9257-2 .

[24] LUKES D L. Differential Equations: Classical to Controlled[M]. New York: Academic Press, 1982.

[25] MathWorks The, Inc. MATLAB Version 8.2.0.701 (R2013b)[CP]. Natick, Massachusetts: The MathWorks, Inc., 13 August 2013.

[26] LI K Z, ZHU G H, MA Z J, et al. Dynamic Stability of an SIQS Epidemic Network and Its Optimal Control[J]. Commun Nonlinear Sci Numer Simul, 2019, 66: 84-95. DOI: 10.1016/j.cnsns.2018.06.020 .

基金资助

国家自然科学基金(11301491)

山西省自然科学基金(2018010221040)