考虑扩散的SVEITR肺结核传播模型的最优控制策略

图1 肺结核在人群中的传播流程图

又因为人口运动和空间分布的交互作用在结核病传播中发挥着重要作用, 因此每个仓室加入扩散项.而在实际中, 为了达到更好的防控效果, 相关部门付出了很大努力, 例如: 开展肺结核预防和控制的宣传教育以提高个人防护意识 ($u_{1}$),及时的检测和治疗 ($u_{2}$), 提高治疗水平 ($u_{3}$), 以及降低复发率 ($u_{4}$) 等.将扩散和控制措施引入图 1, 可以得到下面模型

(2.1)$\begin{equation} \left\{\begin{array}{l} \left.\begin{array}{l} \frac{\partial S}{\partial t}-D_{1}\Delta S=(1-q_{1})\Lambda-\beta(1-u_{1}) S(I+\rho T)-\mu S,\\ \frac{\partial V}{\partial t}-D_{2}\Delta V=q_{1}\Lambda-\sigma\beta(1-u_{1}) V(I+\rho T)-\mu V,\\ \frac{\partial E}{\partial t}-D_{3}\Delta E=\eta\beta(1-u_{1}) S(I+\rho T)+\sigma\beta(1-u_{1}) V(I+\rho T)\\ \quad\quad\quad\quad\quad\quad\quad+\delta_{1}I+\theta_{1}T-(\alpha+\mu)E,\\ \frac{\partial I}{\partial t}-D_{4}\Delta I=(1-\eta)\beta(1-u_{1}) S(I+\rho T)+\alpha E+\gamma(1-u_{4}) R\\ \quad\quad\quad\quad\quad\quad\quad-(\delta_{1}+\delta_{2}(1+u_{2})+\mu+d_{1})I,\\ \frac{\partial T}{\partial t}-D_{5}\Delta T=\delta_{2}(1+u_{2}) I-(\theta_{1}+\theta_{2}(1+u_{3})+\mu+d_{2})T,\\ \frac{\partial R}{\partial t}-D_{6}\Delta R=\theta_{2}(1+u_{3}) T-(\gamma(1-u_{4})+\mu)R,\\ \end{array}\right\}(x,t)\in Q=\Omega\times(0,T),\\ \left.\begin{array}{l} S(x,0)=S_{0}(x)\geq0, V(x,0)=V_{0}(x)\geq0, E(x,0)=E_{0}(x)\geq0, \\ I(x,0)=I_{0}(x)\geq0, T(x,0)=T_{0}(x)\geq0, R(x,0)=R_{0}(x)\geq0,\\ \end{array}\right\} x\in\Omega,\\ \frac{\partial S}{\partial \textbf{n}}=0, \frac{\partial V}{\partial \textbf{n}}=0, \frac{\partial E}{\partial \textbf{n}}=0, \frac{\partial I}{\partial \textbf{n}}=0, \frac{\partial T}{\partial \textbf{n}}=0, \frac{\partial R}{\partial \textbf{n}}=0, (x,t)\in\partial\Omega\times(0,T), \end{array}\right. \end{equation}$

其中$\Omega$表示$\mathbb{R}^{n}$中的一个有界空间区域, 且具有光滑边界$\partial\Omega.$$\Delta$表示 Laplacian 算子. $\textbf{n}$表示$\partial\Omega$的外向单位法向量. 模型 (2.1) 中各参数含义见表 1.

表1 模型 (2.1) 各参数含义

注: 表 1 中各参数取值范围和来源主要参考文献 [7].

新窗口打开| 下载CSV

对于最优控制问题, 我们希望以最少的花费实现最少的感染人数, 因此建立如下目标泛函

(2.2)$\begin{equation} \begin{array}{l} J(u_{1},u_{2},u_{3},u_{4})=\int_{0}^{T}\int_{\Omega}\big[A_{1}E+A_{2}I+\sum\limits_{i=1}^{4}\frac{1}{2}B_{i}u_{i}^{2}\big]\text{d}x\text{d}t\\ \quad\quad\quad\quad\quad\quad\quad =\int_{Q}\big[A_{1}E+A_{2}I+\sum\limits_{i=1}^{4}\frac{1}{2}B_{i}u_{i}^{2}\big]\text{d}x\text{d}t, \end{array} \end{equation}$

其中$T$为终端时刻,$A_{i}(i=1,2)$分别为$E$和$I$的权重常数,$B_{i}(i=1,2,3,4)$分别为四种控制的权重常数.

控制容许集为

$U=\{(u_{1},u_{2},u_{3},u_{4})\in L^{2}(Q):0\leq u_{i}\leq M_{i}, i=1,2,3,4\}.$

其中$M_{i}$表示第$i$个控制变量$u_{i}$的上限值.

我们需要寻找最优控制$(u_{1}^{*},u_{2}^{*},u_{3}^{*},u_{4}^{*})\in U,$使得

(2.3)$J(u_{1}^{*},u_{2}^{*},u_{3}^{*},u_{4}^{*})=\inf\limits_{U} J(u_{1},u_{2},u_{3},u_{4}).$

由文献 [16] 第 25 章 25.1 部分的内容可知: 模型 (2.1) 中的任意状态解在$\Omega\times(0,T)$均为非负的;而在$L^{2}((0,T);H^{1}(\Omega))$(状态函数所属空间) 和 $L^{2}((0,T);H^{1}(\Omega)^{*})$(偏导函数所属空间) 内,系统 (2.1) 的弱解是存在唯一的; 并且最优控制问题 (2.1), (2.2), (2.3) 存在唯一的最优控制.

3 问题求解

下面求解最优控制问题 (2.1), (2.2), (2.3).由文献 [16] 知道可以使用一组伴随偏微分方程来求解最优控制和相应的状态系统: 伴随函数类似于连接状态系统和成本目标泛函的 Lagrange 乘子. 通过求解由状态偏微分方程、伴随偏微分方程和最优控制特征组成的最优性系统, 就可以获得最优控制问题的解.

3.1 灵敏度系统

构造 “控制-状态” 映射$(u_{1},u_{2},u_{3},u_{4})\rightarrow \textbf{X}(u_{1},u_{2},u_{3},u_{4}),$其中$\textbf{X}=(S,V,E,I,T,R).$对于控制$u_{i}(i=1,2,3,4),$考虑另一个控制$u_{i}^{\epsilon}=u_{i}+\epsilon l_{i},$其中$l_{i}$是一个变分函数且$\epsilon>0.$令$X_{j}=X_{j}(u_{1},u_{2},u_{3},u_{4}),$$X_{j}^{\epsilon}=X_{j}(u_{1}^{\epsilon},u_{2}^{\epsilon},u_{3}^{\epsilon},u_{4}^{\epsilon}),$其中$X_{j}$是$\textbf{X}$的第$j(j=1,2,\cdots,6)$个分量.

下面定义灵敏度函数{$\psi_{i}~(i=1,2,\cdots,6)$,} 它们是 “控制-状态” 映射的 Gateaux 导数. 为此, 构造差商

(3.1)$\frac{X_{j}^{\epsilon}-X_{j}}{\epsilon},$

满足$\epsilon\rightarrow0$ 时,$X_{j}^{\epsilon}\rightarrow X_{j}$且$\frac{X_{j}^{\epsilon}-X_{j}}{\epsilon}\rightarrow \psi_{j} (j=1,2,\cdots,6).$

差商(3.1)的收敛性以及灵敏度函数及其对应偏微分方程的存在性可参考文献 [16]. 灵敏度函数满足与状态方程的线性化形式相对应的系统, 下面计算$\psi_{1}$满足的偏微分方程.对于差商$\frac{S^{\epsilon}-S}{\epsilon},$由状态系统(2.1)的第一个方程可知

$\frac{\partial S^{\epsilon}}{\partial t}-D_{1}\Delta S^{\epsilon}=(1-q_{1})\Lambda-\beta(1-(u_{1}+\epsilon l_{1})) S^{\epsilon}(I^{\epsilon}+\rho T^{\epsilon})-\mu S^{\epsilon},$

$\frac{\partial S}{\partial t}-D_{1}\Delta S=(1-q_{1})\Lambda-\beta(1-u_{1}) S(I+\rho T)-\mu S.$

所以

(3.2)$\begin{aligned} &~~~~\frac{\frac{\partial S^{\epsilon}}{\partial t}-\frac{\partial S}{\partial t}}{\epsilon}-D_{1}\frac{\Delta S^{\epsilon}-\Delta S}{\epsilon}\\ &=\frac{\partial \left(\frac{S^{\epsilon}-S}{\epsilon}\right)}{\partial t}-D_{1}\Delta\left(\frac{S^{\epsilon}-S}{\epsilon}\right)\\ &=-\beta(1-u_{1})\left[\frac{(S^{\epsilon}-S)I^{\epsilon}}{\epsilon}+\frac{S(I^{\epsilon}-I)}{\epsilon} +\rho\frac{(S^{\epsilon}-S)T^{\epsilon}}{\epsilon}+\rho\frac{S(T^{\epsilon}-T)}{\epsilon}\right]\\ & \quad +\beta l_{1}S^{\epsilon}(I^{\epsilon}+\rho T^{\epsilon})-\mu\frac{S^{\epsilon}-S}{\epsilon}. \end{aligned}$

令$\epsilon\rightarrow0,$对 (3.2) 式取极限可得

$\begin{align*} \frac{\partial \psi_{1}}{\partial t}-D_{1}\Delta\psi_{1} &=-\beta(1-u_{1})(\psi_{1}I+S\psi_{4}+\rho\psi_{1}T+\rho S\psi_{5})+\beta l_{1}S(I+\rho T)-\mu \psi_{1}\\ &=[-\beta(1-u_{1})(I\!+\!\rho T)-\mu]\psi_{1}-\beta(1-u_{1})S\psi_{4}-\beta(1-u_{1})\rho S\psi_{5}\!+\!l_{1}\beta S(I\!+\!\rho T). \end{align*}$

对于$\psi_{k} (k=2,3,4,5,6)$的偏微分方程可类似计算. 则可相应得到关于 {${\psi}=(\psi_{1},\psi_{2},\psi_{3},\psi_{4},$$\psi_{5},\psi_{6})$} 的偏微分方程组, 其向量形式如下

(3.3)$\begin{equation} \mathcal{L}\left( \begin{array}{c} \psi_{1} \\ \psi_{2} \\ \psi_{3} \\ \psi_{4} \\ \psi_{5} \\ \psi_{6} \\ \end{array} \right)=\left( \begin{array}{c} l_{1}\beta S(I+\rho T) \\ l_{1}\sigma\beta V(I+\rho T) \\ -l_{1}\eta\beta S(I+\rho T)-l_{1}\sigma\beta V(I+\rho T) \\ -l_{1}(1-\eta)\beta S(I+\rho T)-l_{2}\delta_{2}I-l_{4}\gamma R \\ l_{2}\delta_{2}I-l_{3}\theta_{2}T \\ l_{3}\theta_{2}T+l_{4}\gamma R \\ \end{array} \right),(x,t)\in Q=\Omega\times(0,T), \end{equation}$

在 (3.3) 式中

$\mathcal{L}\left( \begin{array}{c} \psi_{1} \\ \psi_{2} \\ \psi_{3} \\ \psi_{4} \\ \psi_{5} \\ \psi_{6} \\ \end{array} \right)= \left( \begin{array}{c} \frac{\partial \psi_{1}}{\partial t}-D_{1}\Delta \psi_{1} \\ \frac{\partial \psi_{2}}{\partial t}-D_{2}\Delta \psi_{2} \\ \frac{\partial \psi_{3}}{\partial t}-D_{3}\Delta \psi_{3} \\ \frac{\partial \psi_{4}}{\partial t}-D_{4}\Delta \psi_{4} \\ \frac{\partial \psi_{5}}{\partial t}-D_{5}\Delta \psi_{5} \\ \frac{\partial \psi_{6}}{\partial t}-D_{6}\Delta \psi_{6} \\ \end{array} \right)+\textbf{M}\left( \begin{array}{c} \psi_{1} \\ \psi_{2} \\ \psi_{3} \\ \psi_{4} \\ \psi_{5} \\ \psi_{6} \\ \end{array} \right),$

且

$\textbf{M}=(\textbf{M}_{1},\textbf{M}_{2}),$

其中

$\textbf{M}_{1}=\left( \begin{array}{ccc} \beta(1-u_{1})(I+\rho T)-\mu & 0 & 0 \\ 0 & \sigma\beta(1-u_{1})(I+\rho T)-\mu & 0 \\ -\eta\beta(1-u_{1})(I+\rho T) & -\sigma\beta(1-u_{1})(I+\rho T) & \alpha+\mu \\ -(1-\eta)\beta(1-u_{1})(I+\rho T) & 0 & -\alpha \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ \end{array} \right),$

$\textbf{M}_{2}=\left( \begin{array}{ccc} \beta(1-u_{1})S & \beta(1-u_{1})\rho S & 0 \\ \sigma\beta(1-u_{1})V & \sigma\beta(1-u_{1})\rho V & 0 \\ -\eta\beta(1-u_{1})S & -\eta\beta(1-u_{1})\rho S & 0 \\ -\sigma\beta(1-u_{1})V-\delta_{1}& -\sigma\beta(1-u_{1})\rho V-\theta_{1}&\\ -(1-\eta)\beta(1-u_{1})S & -(1-\eta)\beta(1-u_{1})\rho S & -\gamma(1-u_{4}) \\ +(\delta_{1}+\delta_{2}(1+u_{2})+\mu+d_{1})& &\\ -\delta_{2}(1+u_{2}) & \theta_{1}+\theta_{2}(1+u_{3})+\mu+d_{2} & 0 \\ 0 & -\theta_{2}(1+u_{3}) & \gamma(1-u_{4})+\mu \\ \end{array} \right).$

可以发现灵敏度函数满足状态方程(2.1)的线性化形式相对应的系统.

初值条件为

(3.4)$\psi_{1}(x,0)=\psi_{2}(x,0)=\psi_{3}(x,0)=\psi_{4}(x,0)=\psi_{5}(x,0)=\psi_{6}(x,0)=0, x\in\Omega.$

边界条件为

(3.5)$\frac{\partial \psi_{1}}{\partial \textbf{n}}=\frac{\partial \psi_{2}}{\partial \textbf{n}}= \frac{\partial \psi_{3}}{\partial \textbf{n}}=\frac{\partial \psi_{4}}{\partial \textbf{n}}= \frac{\partial \psi_{5}}{\partial \textbf{n}}=\frac{\partial \psi_{6}}{\partial \textbf{n}}=0, (x,t)\in\partial\Omega\times(0,T).$

3.2 伴随系统

对于伴随变量${\lambda}=(\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4},\lambda_{5},\lambda_{6}),$在合适的弱$L^{2}(Q)$意义下, 伴随算子可由下面关系得到

(3.6)$\langle{\lambda},\mathcal{L}{\psi}\rangle=\langle\mathcal{L}^{*}{\lambda},{\psi}\rangle.$

伴随系统中各方程等号的右端表达式由目标泛函 (2.2) 中的被积函数分别关于状态变量求偏导得到.下面使用表达式 (3.6) 来导出伴随系统中的算子, 通过在$\langle{\lambda},\mathcal{L}{\psi}\rangle$上使用分部积分, 再结合灵敏度系统的边值条件(3.5), 可得到如下向量形式的伴随系统

(3.7)$\begin{equation} \mathcal{L}^{*}\left( \begin{array}{c} \lambda_{1} \\ \lambda_{2} \\ \lambda_{3} \\ \lambda_{4} \\ \lambda_{5} \\ \lambda_{6} \\ \end{array} \right)=\left( \begin{array}{c} 0 \\ 0 \\ A_{1} \\ A_{2} \\ 0 \\ 0 \\ \end{array} \right),(x,t)\in Q=\Omega\times(0,T), \end{equation}$

其中

$\mathcal{L}^{*}\left( \begin{array}{c} \lambda_{1} \\ \lambda_{2} \\ \lambda_{3} \\ \lambda_{4} \\ \lambda_{5} \\ \lambda_{6} \\ \end{array} \right)=\left( \begin{array}{c} -\frac{\partial \lambda_{1}}{\partial t}-D_{1}\Delta \lambda_{1} \\ -\frac{\partial \lambda_{2}}{\partial t}-D_{2}\Delta \lambda_{2} \\ -\frac{\partial \lambda_{3}}{\partial t}-D_{3}\Delta \lambda_{3} \\ -\frac{\partial \lambda_{4}}{\partial t}-D_{4}\Delta \lambda_{4} \\ -\frac{\partial \lambda_{5}}{\partial t}-D_{5}\Delta \lambda_{5} \\ -\frac{\partial \lambda_{6}}{\partial t}-D_{6}\Delta \lambda_{6} \\ \end{array} \right)+\textbf{M}^{T}\left( \begin{array}{c} \lambda_{1} \\ \lambda_{2} \\ \lambda_{3} \\ \lambda_{4} \\ \lambda_{5} \\ \lambda_{6} \\ \end{array} \right).$

其初值条件为

(3.8)$\lambda_{1}(x,T)=\lambda_{2}(x,T)=\lambda_{3}(x,T)=\lambda_{4}(x,T)=\lambda_{5}(x,T)=\lambda_{6}(x,T)=0, x\in\Omega.$

边界条件为

(3.9)$\begin{equation} \frac{\partial \lambda_{1}}{\partial \textbf{n}}=\frac{\partial \lambda_{2}}{\partial \textbf{n}}= \frac{\partial \lambda_{3}}{\partial \textbf{n}}=\frac{\partial \lambda_{4}}{\partial \textbf{n}}= \frac{\partial \lambda_{5}}{\partial \textbf{n}}=\frac{\partial \lambda_{6}}{\partial \textbf{n}}=0, (x,t)\in\partial\Omega\times(0,T). \end{equation}$

3.3 最优控制的特征

记$\textbf{u}=(u_{1},u_{2},u_{3},u_{4}),$$\textbf{l}=(l_{1},l_{2},l_{3},l_{4}).$映射$\textbf{u}\rightarrow J(\textbf{u})$在$\textbf{l}$方向上关于$\textbf{u}^{*}$的导数为

$\lim\limits_{\epsilon\rightarrow0^{+}}\frac{J(\textbf{u}^{*}+\epsilon \textbf{l})-J(\textbf{u}^{*})}{\epsilon}\geq0,$

因为$J(u^{*})$是目标泛函$J(u)$的极小值, 所以在最优控制${u}^{*}$处, 上述导数是非负的.

易知

$J(\textbf{u}^{*}+\epsilon \textbf{l})=\int_{Q}\big[A_{1}E^{\epsilon}+A_{2}I^{\epsilon}+\sum\limits_{i=1}^{4}\frac{1}{2}B_{i}(u_{i}^{*}+\epsilon l_{i})^{2}\big]\text{d}x\text{d}t,$

$J(\textbf{u}^{*})=\int_{Q}\big[A_{1}E+A_{2}I+\sum\limits_{i=1}^{4}\frac{1}{2}B_{i}(u_{i}^{*})^{2}\big]\text{d}x\text{d}t.$

记 $W=B_{1}l_{1}u_{1}^{*}+B_{2}l_{2}u_{2}^{*}+B_{3}l_{3}u_{3}^{*}+B_{4}l_{4}u_{4}^{*},$ 则

(3.10)$\begin{array}{l} \quad \lim\limits_{\epsilon\rightarrow0^{+}}\frac{J(\textbf{u}^{*}+\epsilon \textbf{l})-J(\textbf{u}^{*})}{\epsilon}\\ =\lim\limits_{\epsilon\rightarrow0^{+}}\int_{Q}\big[A_{1}(\frac{E^{\epsilon}-E}{\epsilon})+A_{2}(\frac{I^{\epsilon}-I}{\epsilon}) +\sum\limits_{i=1}^{4}\frac{1}{2}B_{i}[\frac{(u_{i}^{*}+\epsilon l_{i})^{2}-(u_{i}^{*})^{2}}{\epsilon}]\big]\text{d}x\text{d}t\\ =\int_{Q}\big[A_{1}\psi_{3}+A_{2}\psi_{4}+W\big]\text{d}x\text{d}t\\ =\int_{Q}\big[(\psi_{1},\psi_{2},\psi_{3},\psi_{4},\psi_{5},\psi_{6})(0,0,A_{1},A_{2},0,0)^{T}+W\big]\text{d}x\text{d}t\\ =\int_{Q}\big[(\psi_{1},\psi_{2},\psi_{3},\psi_{4},\psi_{5},\psi_{6})\mathcal{L}^{*}(\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4},\lambda_{5},\lambda_{6})^{T}+W\big]\text{d}x\text{d}t\\ =\int_{Q}\big[(\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4},\lambda_{5},\lambda_{6})\mathcal{L}(\psi_{1},\psi_{2},\psi_{3},\psi_{4},\psi_{5},\psi_{6})^{T}+W\big]\text{d}x\text{d}t\\ =\int_{Q}\big\{l_{1}\big[\beta[{\lambda_{1}S+\lambda_{2}\sigma V-\lambda_{3}\eta S-\lambda_{3}\sigma V-\lambda_{4}(1-\eta) S ](I+\rho T)}+B_{1}u_{1}^{*}\big]\big. \\ \quad +l_{2}[{(\lambda_{5}-\lambda_{4})\delta_{2}I+B_{2}u_{2}^{*}]+l_{3}[(\lambda_{6}-\lambda_{5})\theta_{2}T+B_{3}u_{3}^{*}] +l_{4}[(\lambda_{6}-\lambda_{4})\gamma R}+B_{4}u_{4}^{*}]\big\}\text{d}x\text{d}t.\\ \end{array}$

关于最优控制的特征可以得到如下定理

定理 3.1 在控制容许集 $U$ 上使得 $J(u_{1},u_{2},u_{3},u_{4})$ 取得最小值的$(u_{1}^{*},u_{2}^{*},u_{3}^{*},u_{4}^{*})$为

$u_{1}^{*}=\min\left(M_{1},\max\left(\frac{\beta[\lambda_{3}^*\eta S^*+\lambda_{3}^*\sigma V^*+\lambda_{4}^*(1-\eta) S^* -\lambda_{1}^*S^*-\lambda_{2}^*\sigma V^*](I^*+\rho T^*)}{B_{1}},0\right)\right),$

$u_{2}^{*}=\min\left(M_{2},\max\left(\frac{(\lambda_{4}^*-\lambda_{5}^*)\delta_{2}I^*}{B_{2}},0\right)\right),$

$u_{3}^{*}=\min\left(M_{3},\max\left(\frac{{(\lambda_{5}^*-\lambda_{6}^*)\theta_{2}T^*}}{B_{3}},0\right)\right),$

$u_{4}^{*}=\min\left(M_{4},\max\left(\frac{{(\lambda_{4}^*-\lambda_{6}^*)\gamma R^*}}{B_{4}},0\right)\right),$

其中 $(S^*,V^*,E^*,I^*,T^*,R^*)$ 分别是相应状态变量的最优解, $(\lambda_1^*,\lambda_2^*,\lambda_3^*,\lambda_4^*,\lambda_5^*,\lambda_6^*)$ 是系统 (3.7) 在条件 (3.8)-(3.9) 下的解.

证由文献 [16] 可知, 因为在控制容许集$U$上, 最优控制变量可以在其有界限内调整以优化目标泛函 (2.2). 因此, 变分$\textbf{l}$可以具有任意符号.所以在集合$U$上, 有

$\beta[\lambda_{1}S+\lambda_{2}\sigma V-\lambda_{3}\eta S-\lambda_{3}\sigma V-\lambda_{4}(1-\eta) S ](I+\rho T)+B_{1}u_{1}^{*}=0,$

$(\lambda_{5}-\lambda_{4})\delta_{2}I+B_{2}u_{2}^{*}=0,$

$(\lambda_{6}-\lambda_{5})\theta_{2}T+B_{3}u_{3}^{*}=0,$

$(\lambda_{6}-\lambda_{4})\gamma R+B_{4}u_{4}^{*}=0.$

所以当$B_{i}\neq0$时, 即可得到最优控制$u_{i}^{*}~(i=1,2,3,4)$的表达式.

4 数值模拟

本部分对理论结果进行可视化展示, 以下所有图中同一仓室在空间上不同位置处初值一样,分别取为:$S(:,1)=501000,V(:,1)=79300,E(:,1)=126000,I(:,1)=36,T(:,1)=10,R(:,1)=5,$以上初值单位为 “个”.参数取值 (无单位) 分别为: $\mu=1/77.11,~\Lambda=1,q_{1}=0.9,~\beta=9.0000e-06,~\sigma=0.0011,~\rho=0.0102, \eta=0.9380,~\delta_{1}=0.7,~\delta_{2}=0.3,~\theta_{1}=0.5,~\theta_{2}=0.5, \alpha=0.5002;~\gamma=0.0100,d_{1}=7.0600e-04, d_{2}=7.8300e-04, D_{1}=0.3,D_{2}=0.3,D_{3}=0.3,D_{4}=0.1,D_{5}=0.2,D_{6}=0.3.$ 目标泛函中权重系数取值分别为: $A_1=0.2,A_2=0.2,B_1=1,B_2=1,B_3=1,B_4=1.$ 控制容许集中 $M_{i}=0.8~(i=1,2,3,4).$ 空间总长取为$x=10$ (公里), 步长取为$0.2.$ 时间总长取为$t=30$ (天), 步长取为$0.5.$

在最优控制下, 四个控制变量$u_{i}~(i=1,2,3,4)$的时间序列图见图 2. 为了更直观地看出最优控制变量的变化情况, 可在图 2 中截取$x=8$时的时间序列剖面图, 具体见图 3. 从图 3 可以看出, 对于控制变量$u_{1}$, 在一开始选取最大控制值, 持续 2 天后即可持续下降至 0.对于控制变量$u_{2}$, 大约半天内从 0 增到最大控制, 持续约 2 天后, 持续下降至 0.对于控制变量$u_{3}$, 一开始即需要维持在最高控制水平, 持续约 28 天后才可持续下降.对于控制变量$u_{4}$, 首先需要在最大控制水平维持约 2.5 天, 随后即可持续下降.

图2

图2 最优控制变量的时间序列图.Time (单位: 天), $x$ (单位: 公里)

图3

图3 最优控制变量在 $x=8$ 时的时间序列剖面图.Time (单位: 天)

在上述最优控制水平下, 对于各状态变量的控制效果见图 4. 可以看出: 在最优控制下, 潜伏者、TB 活跃者和接受治疗者的数量均少于无控制的情形,而恢复者的数量却多于无控制的情形, 这说明最优综合控制策略对于肺结核的传播起到了较好的控制效果.

图4

图4 在 $x=8$ 时最优控制和无控制情况下状态变量值的比较.Time (单位: 天)

对于单一控制和综合控制, 其控制效果见图 5. 该图显示了两种控制策略下染病者数量在一个月内的变化情况, 可以看出:相比只有 $u_{1}$ 存在的情况, 综合控制下染病者的峰值要低, 且最终感染数量也要少;相比只有 $u_{2}$ 存在的情况, 综合控制下虽然染病者峰值与仅 $u_{2}$ 存在时相差不大, 但最终感染数量要少;相比只有 $u_{3} $ 存在的情况, 综合控制下虽然染病者最终感染数量与仅 $u_{3}$ 存在时相差不大, 但峰值要低, 达峰时间要晚;相比只有 $u_{4}$ 存在的情况, 综合控制下染病者的峰值要低,达峰时间要晚, 且最终感染数量也要少.这充分说明最优综合控制策略的效果要优于单一控制策略.

图5

图5 在 $x=8$ 时单一控制和综合控制情况下染病者数量的比较.Time (单位: 天)

不同控制组合下的控制效果对比见图 6. 该图中给出了三种控制措施同时存在的情况下染病者的数量.共有四组不同的控制组合, 分别为: 组合 1: $u_{1},u_{2},u_{3}\neq0, u_{4}=0;$组合 2: $u_{1},u_{2},u_{4}\neq0, u_{3}=0;$ 组合 3: $u_{1},u_{3},u_{4}\neq0, u_{2}=0;$ 组合 4: $u_{2},u_{3},u_{4}\neq0, u_{1}=0.$从图中可以看出: 组合 2 的控制效果最好, 其次为组合 1, 而组合 3 和组合 4 控制效果差异不大.

图6

图6 在 $x=8$ 时四种不同控制策略组合下的染病者数量比较.Time (单位: 天)

根据上述结果分析, 我们可以提出以下实际操作建议, 以优化肺结核的防控措施: (1) 个人防护意识提升方面: 在疫情初期, 迅速加强公众的防护意识教育, 通过宣传和教育活动使个人防护措施 (如佩戴口罩、保持社交距离等) 得到广泛普及. 建议在疾病暴发初期投入最大宣传力度, 然后根据防控效果逐渐减弱宣传强度. (2) 检测和治疗能力的提升方面: 在疫情初期, 快速提高检测和治疗的能力, 确保更多潜在病例能被及时发现并进行治疗. 此措施在初期保持高强度, 然后根据疫情的发展逐步减少检测频率, 但仍需保持一定的检测能力以防止疫情反弹. (3) 治疗水平的持续改进方面: 在整个疫情控制阶段, 保持高水平的治疗质量和能力, 这意味着需要足够的医疗资源和人员培训, 以确保患者能够得到及时有效的治疗. (4) 降低复发率的措施方面: 在疾病暴发初期, 集中采取措施降低复发率, 这可能包括针对康复者的后续监控和必要的医学干预, 以降低复发可能性. 总之, 应结合各方面的措施以达到更好的综合防控效果, 综合策略应根据疫情变化进行灵活调整, 确保各个环节的协调和优化. 另外, 还应建立持续的监测机制, 对防控措施的效果进行评估和反馈, 及时调整策略以适应新的挑战.

5 结束语

数学建模在肺结核研究中的应用是多方面的, 包括评估疫情发展、制定和评价干预策略、优化资源配置等. 通过对模型的不断完善和数据的更新, 数学建模将继续在控制肺结核传播和防治策略的制定中发挥关键作用. 本文对一个 SVEITR 模型加入了扩散和 4 个控制变量, 研究扩散存在下的最优控制策略. 理论分析得到了最优控制的特征, 并数值展示了最优控制变量的变化情况以及控制效果. 所得结果可为肺结核防控提供理论建议. 然而, 对于肺结核的研究仍有许多需要解决的问题, 例如, 在突发呼吸道传染病的各类非药物干预策略下, 我国肺结核在预防、传播、治疗和控制等各阶段将受到怎样的影响, 其发病率和死亡率如何变化, 如何建立动力学模型, 并定性研究其内在传播机制等均有待进一步研究.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Global Tuberculosis report 2023. Geneva: World Health Organization; 2023. Licence:CC BY-NC-SA 3.O IGO

[2]

Mao

K X

, Zhen

C X

, Hong

Yan L U

, et al.

Protective eff ect of vaccination of Bacille Calmette-Gnerin on children

Chin J Contemp Pediatr, 2003, 5(4): 325-327

[3]

Abubakar

, Pimpin

, et al.

Systematic review and meta-analysis of the current evidence on the duration of protection by bacillus Calmette-Guerin vaccination against tuberculosis

Health Technol Assess, 2013, 17(37): 1-372

DOI:10.3310/hta17370 PMID:24021245 [本文引用: 2]

Recent evidence suggests that the duration of protection by bacillus Calmette-Guérin (BCG) may exceed previous estimates with potential implications for estimating clinical and cost-efficacy.To estimate the protection and duration of protection provided by BCG vaccination against tuberculosis, explore how this protection changes with time since vaccination, and examine the reasons behind the variation in protection and the rate of waning of protection.Electronic databases including MEDLINE, Excerpta Medica Database (EMBASE), Cochrane Databases, NHS Economic Evaluation Database (NHS EED), Database of Abstracts of Reviews of Effects (DARE), Web of Knowledge, Biosciences Information Service (BIOSIS), Latin American and Caribbean Health Sciences Literature (LILACs), MEDCARIB Database, Cumulative Index to Nursing and Allied Health Literature (CINAHL) were searched from inception to May 2009. Index to Theses, System for Information on Grey Literature in Europe (SIGLE), Centre for Agricultural Bioscience International (CABI) Abstracts, Scopus, Article First, Academic Complete, Africa-Wide Information, Google Scholar, Global Health, British National Bibliography for Report Literature, and clinical trial registration websites were searched from inception to October 2009.Electronic databases searches, screening of identified studies, data extraction and analysis were undertaken. Meta-analysis was used to present numerical and graphical summaries of clinical efficacy and efficacy by time since vaccination. Evidence of heterogeneity was assessed using the tau-squared statistic. Meta-regression allowed the investigation of observed heterogeneity. Factors investigated included BCG strain, latitude, stringency of pre-BCG vaccination tuberculin testing, age at vaccination, site of disease, study design and vulnerability to biases. Rate of waning of protection was estimated using the ratio of the measure of efficacy after 10 years compared with the efficacy in the first 10 years of a study.Study selection. A total of 21,030 references were identified, providing data on 132 studies after abstract and full-text review. Efficacy. Protection against pulmonary tuberculosis in adults is variable, ranging from substantial protection in the UK MRC trial {rate ratio 0.22 [95% confidence interval (CI) 0.16 to 0.31]}, to absence of clinically important benefit, as in the large Chingleput trial [rate ratio 1.05 (95% CI 0.88 to 1.25)] and greater in latitudes further away from the equator. BCG vaccination efficacy was usually high, and varied little by form of disease (with higher protection against meningeal and miliary tuberculosis) or study design when BCG vaccination was given only to infants or to children after strict screening for tuberculin sensitivity. High levels of protection against death were observed from both trials and observational studies. The observed protective effect of BCG vaccination did not differ by the strain of BCG vaccine used in trials.Reviewed studies showed that BCG vaccination protects against pulmonary and extrapulmonary tuberculosis for up to 10 years. Most studies either did not follow up participants for long enough or had very few cases after 15 years. This should not be taken to indicate an absence of effect: five studies (one trial and four observational studies) provided evidence of measurable protection at least 15 years after vaccination. Efficacy declined with time. The rate of decline was variable, with faster decline in latitudes further from the equator and in situations where BCG vaccination was given to tuberculin-sensitive participants after stringent tuberculin testing.The main limitation of this review relates to quality of included trials, most of which were conducted before current standards for reporting were formulated. In addition, data were lacking in some areas and the review had to rely on evidence from observational studies.BCG vaccination protection against tuberculosis varies between populations, to an extent that cannot be attributed to chance alone. Failure to exclude those already sensitised to mycobacteria and study latitude closer to the equator were associated with lower efficacy. These factors explained most of the observed variation. There is good evidence that BCG vaccination protection declines with time and that protection can last for up to 10 years. Data on protection beyond 15 years are limited; however, a small number of trials and observational studies suggest that BCG vaccination may protect for longer. Further studies are required to investigate the duration of protection by BCG vaccination.The National Institute for Health Research Health Technology Assessment programme.

[4]

徐瑞, 周凯娟, 白宁.

一类基于游离病毒感染和细胞-细胞传播的宿主体内 HIV-1 感染动力学模型

数学物理学报, 2024, 44A(3): 771-782

, Zhou

K J

, Bai

A with-in host HIV-1 infection dynamics model based on virus-to-cell infection and cell-to-cell transmission

Acta Math Sci, 2024, 44A(3): 771-782

[5]

杨俊元, 张晨琳, 杨丽.

年龄结构流感模型综合控制策略研究

数学物理学报, 2024, 44A(3): 804-814

Yang

J Y

, Zhang

C L

, Yang

Study on comprehensive control strategies of an age structured influenza model

Acta Math Sci, 44A (3): 804-814

[6]

Liu

J L

, Zhang

T L

Global stability for a tuberculosis model

Math Comput Model, 2011, 54: 836-845

[7]

, Liu

X N

, Yuan

Y Y

, et al.

Global analysis of tuberculosis dynamical model and optimal control strategies based on case data in the United States

Appl Math Comput, 2022, 422: 126983

[本文引用: 4]

[8]

Yang

Y L

, Tang

S Y

, Ren

X H

, et al.

Global stability and optimal control for a tuberculosis model with vaccination and treatment

Discrete Contin Dyn Syst Ser B, 2017, 21(3): 1009-1022

[9]

Gao

D P

, Huang

N J

A note on global stability for a tuberculosis model

Appl Math Lett, 2017, 73: 163-168

[10]

Liu

S Y

, Bi

Y J

, Liu

Y W

Modeling and dynamic analysis of tuberculosis in mainland China from 1998 to 2017: the effect of DOTS strategy and further control

Theor Biol Med Model, 2020, 17(6): 1-10

[11]

Gao

D P

, Huang

N J

Optimal control analysis of a tuberculosis model

Appl Math Model, 2018, 58: 47-64

[12]

Egonmwan

A O

, Okuonghae

Mathematical analysis of a tuberculosis model with imperfect vaccine

Int J Biomath, 2019, 12(7): 1950073

[13]

Bhih

A E

, Laaroussi

A E A

, Ghazzali

, et al.

An optimal chemoprophylaxis and treatment control for a spatiotemporal tuberculosis model

Commun Math Biol Neurosci, 2021, 2021: Article 40

[14]

Wang

Z P

, Xu

Stability and traveling waves of an epidemic model with relapse and spatial diffusion

J Appl Anal Comput, 2014, 4(3): 307-322

[15]

Zhang

, Liu

L L

, Feng

X M

, et al.

Existence of traveling wave solutions for a diffusive tuberculosis model with fast and slow progression

Appl Math Lett, 2021, 112: 106848