基于反应扩散SIR流行病模型的治疗时间最优控制问题

《Mathematics and Computers in Simulation》：An optimal control problem of treatment time from reaction diffusion epidemic model SIR

　　

摘要

本研究聚焦于带有出生率和死亡率的时空????????传染病模型中的最优治疗时间控制问题。状态方程采用扩散型偏微分方程（PDE）进行建模，以分析疾病在连续时空中的传播过程，这对于理解医院和交通系统等人烟稠密区域的疫情传播，以及COVID-19等呼吸道传染病的传播至关重要。我们考虑了一个依赖于自由时间变量的泛函目标，该变量代表优化时的治疗时间，因为治疗时长是每次治疗成功与否的关键因素。我们获得了关于治疗浓度和时长的第一阶和第二阶最优性条件。为说明理论结果的有效性，我们针对疾病从中心或角落开始的几种情景提供了数值模拟。数值模拟采用迭代离散方案进行，并按照类似前向-后向扫描法（FBSM）的适当检验方法收敛。采用了一种算法方法来确定最优治疗时间。

引言

室模型在流行病学中对于描述、控制和分析传染病的传播至关重要。这些模型根据感染阶段将总人口划分为不同的类别。经典模型记为????????，其中??类包括易感疾病的个体且从未被感染。一旦感染，他们转变为??类，该类包括具有传染性且可将疾病传播给易感个体的个体。康复后，个体进入??类，假设他们终身免疫。许多疾病已使用????????模型进行了研究（见，例如，[1]、[2]、[3]、[4]、[5]、[6]、[7]、[8]、[9]、[10]、[11]、[12]）。在[13]中，Gomes等人假设所有新出生人口进入易感类S，自然死亡率与疾病无关，因此各室间相等。标准????????演化系统可表示为：

??? ?S/?t = μN ? βSI ? γS,
?I/?t = βSI ? (γ + δ)I, t ∈ [0, T]
?R/?t = δI ? γR.

其中，N = S + I + R。参数μ表示出生率，β表示感染性参数，项βS(t)I(t)描述了易感个体被感染的速率。参数δ是康复率，γ表示自然死亡率。

众多研究调查了各种类型的????????模型以防止和控制传染病。一些模型按年龄构建，并使用耦合积分微分方程表述[14]。其他模型使用随机过程[3]、[15]进行表述，而一些模型包含时滞[6]、[16]、[17]或可变感染性[18]、[19]。为考虑空间异质性的影响，研究人员在几项工作中[20]、[21]、[22]利用了反应-扩散系统。沿此思路，我们考虑????????反应-扩散模型，这是对????????传染病模型(1)的扩展，以纳入疫情的空间传播。所得模型如下：

??? ?S(t, r)/?t ? α?ΔS(t, r) = μ(r)N(t, r) ? β(r)S(t, r)I(t, r) ? γ(r)S(t, r),
?I(t, r)/?t ? α?ΔI(t, r) = β(r)S(t, r)I(t, r) ? (δ(r) + γ(r))I(t, r), (t, r) ∈ [0, T] × Ω
?R(t, r)/?t ? α?ΔR(t, r) = δ(r)I(t, r) ? γ(r)R(t, r).

通常的Laplace算子，Δ = ?2/?x2 + ?2/?y2，用于描述空间移动性。基于生物学原因，系统(2)的初始数据假设满足：S(0, r) = S? > 0，I(0, r) = I? ≥ 0，且R(0, r) = R? ≥ 0，r ∈ Ω。

本研究的贡献在于使用时空反应-扩散模型，该方法能够准确描述封闭地理区域内疾病的演变。与受限于纯时间动态并施加全局均匀观点的常微分方程（ODE）模型不同，我们的研究克服了这一抽象，以反映呼吸道传染病的物理现实。对于COVID-19等病毒，传染本质上与空间接近度和气溶胶的运动相关联。因此，扩散项Δ对于建模局部密度变化以及感染波前的物理传播至关重要。纯粹的时间分析将无法反映局部感染簇的形成。

模型的一个关键方面是齐次Neumann边界条件的应用：

?S(t, r)/?η = ?I(t, r)/?η = ?R(t, r)/?η = 0, (t, r) ∈ [0, T] × ?Ω

这在数学上是代表地理区域隔离所必需的。通过在边界处施加零通量，这些条件确保人群和病毒都不会越过研究区域的边界。这使得能够分析病原体在封闭栖息地内的生存情况。与Dirichlet边界条件（后者任意施加外部值，仿佛外部环境决定了边界状态）不同，Neumann条件保留了领地的内在动态。它们为孤立城市或国家内的疫情自主性提供了逼真的模拟，从而支持领地遏制策略。

我们的模型结构基于偏微分方程（PDE）反应-扩散系统，旨在将呼吸道感染的动力学与领地的人口动态耦合。与简化的流行病学模型不同，该系统不假设恒定人口或均匀的空间混合。反应动力学由质量作用发病率（βSI）控制。从生物学角度看，这意味着传播取决于局部接触强度，而不是在全人群中均匀稀释。对于气溶胶传播疾病，这一假设尤其合理，因为在近距离下传染性风险显著增加。

人口更新通过比例招募μN建模，其中N = S + I + R。该项表示自然出生过程，并确保易感个体的持续流入，这可能维持病原体的长期持续存在。

参数α?、α?和α?（扩散系数）：分别代表易感、感染和康复个体的空间移动性。扩散算子将地理传播建模为由运动和局部相互作用驱动的波浪状过程。此外，?S/?η表示Ω ? ?2上的外法向导数。

许多研究以各种形式，如治疗、疫苗接种等，在流行病学领域研究了????????反应-扩散模型的最优控制[23]、[24]、[25]、[26]、[27]。目前，最小化治疗时间对于管理COVID-19和埃博拉等传染病至关重要，因为较短的治疗时长会导致较高的治愈率。在本研究中，我们旨在调查针对感染人群的控制方法，形式为最优治疗，同时考虑最优治疗时长。具体来说，我们将以最优治疗的形式研究时空控制，以减少患者数量和治疗的总成本。为实现这一目标，我们考虑以下目标函数：

J((S, I, R), u, ξ) = σ?∫?^ξ ∫_Ω I(t, r) dr dt + θ?‖u‖2_{L2([0,ξ]×Ω)} + ρξ.

在方程(2)中，我们纳入项u(t, r)以表示正在接受治疗的比例。我们假设所有接受治疗的感染者都被转移到康复类，记为R。控制系统的动力学如下：

??? ?S/?t ? α?ΔS = μN ? βSI ? γS,
?I/?t ? α?ΔI = βSI ? (γ + δ)I ? u(t, r)I, (t, r) ∈ [0, T] × Ω
?R/?t ? α?ΔR = δI ? γR + u(t, r)I,
?S/?η = ?R/?η = ?I/?η = 0, (t, r) ∈ [0, T] × ?Ω
S(0, r) = S?, R(0, r) = R? 且 I(0, r) = I?, r ∈ Ω.

在本研究中，疫情动力学由一组反应-扩散方程建模，其中状态变量，即人口密度S(x, t)、I(x, t)和R(x, t)，是空间位置x（以公里km表示）和时间t（以天day表示）的显式函数，单位为人/km2。传播过程由异质感染率β(x)控制，其单位固定为km2/(人·天)。这一特定量纲对于局部发病率项β(x)S(x, t)I(x, t)保持一致性是必需的：将β(x)与两个局部密度的乘积相乘，得到单位面积每天的新感染流量。

公共卫生干预通过时空控制函数u(x, t)整合，定义为t时刻x点的有效治疗率。以天?1表示，u(x, t)代表即时医疗护理的能力。这一表述确保了与环境中的异质人口参数完全一致，如出生率μN、康复率δ和自然死亡率γ，均以天?1定义。最后，移动性由恒定的扩散系数α?、α?、α?控制，单位为km2/天，其与Laplace算子Δ的相互作用保持了系统时空流动的整体连贯性。

目标泛函(5)旨在平衡三个目标：最小化感染者数量、限制干预成本和确保快速的公共卫生响应。其表达如下：

通过引入加权系数σ?、θ和ρ，确保此泛函的量纲一致性，它们作为转换和归一化因子。

最小化感染者数量的优先级（σ?）：系数σ?将感染者的密度（人/km2）转换为社会成本和公共卫生价值。其单位为Cost/(人·天)，确保在时空积分后，该项代表与疾病存在相关的主要惩罚。

控制成本（θ）：参数θ对二次治疗努力进行加权。由于控制u(x, t)是一个速率（天?1），θ的单位为Cost·km2/天2，确保与感染人数项在积分后的量纲平衡。

时长最小化（ρ）：系数ρ（Cost/天）与最终自由时间ξ相关联。它允许将疫情总时长的缩短明确地作为优化目标纳入。

这些单位的平衡使模型能够识别一种优先减少感染者数量的策略u(x, t)，同时确定尽可能短的治疗时间ξ。控制变量u的平方反映了治疗副作用的严重程度。我们的目标是最小化目标泛函，在尽可能短的治疗时间内通过减少感染者数量和最小化控制变量u来实现，其中ξ ∈ [0; T]。我们定义控制集U_ad为：

U_ad = {L∝(0, T; L∝(Ω))/0 < a ≤ u(t, r) ≤ b < 1, a.e. (t, r) ∈ [0, T] × Ω}.

其中，a和b分别代表药物的最小和最大疗效。

在流行病学领域，研究人员使用控制理论方法为各种模型推导自由终端时间最优控制问题。Taesoo Jang等人（2011）讨论了基于共轭梯度的HIV模型自由终端时间最优控制问题[28]，而Gaurav Pachpute和Kriti Saxena探索了HIV病毒的自由最终时间最优控制[29]。Amine Hamdache等人考虑了HIV感染治疗的自由终端时间最优控制问题[30]。这些研究基于Pontryagin最大值原理调查了由常微分方程控制的系统的最优自由时间控制。然而，关于偏微分方程（PDE）传染病模型，特别是反应-扩散系统的研究很少。各种理论研究[31]、[32]、[33]、[34]、[35]、[36]、[37]、[38]、[39]已经探索了从反应-扩散系统出发的最优时间控制，但不是在流行病学的背景下。在本研究中，我们旨在在流行病学背景下，针对????????反应-扩散系统调查时间自由的、最优治疗时间的最优控制。该问题在空间区域内对抗不同类型疾病疫情传播方面具有重大实际意义，特别是对于COVID-19疾病，因为最小化治疗时间是挽救患者生命的主要因素。为实现我们的目标，我们的工作基于控制理论和变分方法。我们首先证明正则性条件以建立状态系统的强解的存在性，然后使用极小序列技术证明控制的存在性。尽管????????反应-扩散系统的书写带来了困难，我们将证明系统关于控制以及目标泛函相对于控制对偶变量、治疗形式和时间的Fréchet可微性，以给出控制的特征描述。

在本文中，我们将第2节专门用于介绍后续各节将使用的符号和预备知识。接下来在第3节，我们建立状态方程的存在性和正则性。随后，在第4节，我们证明最优控制的存在性。第5节专门用于建立线性化状态方程的可微性和唯一可解性。在第6节中，我们提供了所考虑的最优控制问题的必要最优性条件。第7节专门用于第二阶最优性条件。作为应用，第8节呈现了与控制问题相关的数值结果。最后，我们在第9节中结束本文。

符号和预备知识

为了简化本工作余下部分的符号表示，时空柱体[0, T] × Ω记为Q，[0, T] × ?Ω记为Σ。我们令y = (y?, y?, y?) = (S, I, R)，y? = (y??, y??, y??) = (S?, I?, R?)。

本文研究的最优控制问题可以表述为如下形式：

minimize $J(y, u, \xi) = \sigma_1 \int_{\xi_0} \int_{\Omega} y^2(t, r)\, dr\, dt + \theta_2 \|u\|_{L^2([0,\xi]\times\Omega)}^2 + \rho\,\xi,$

$$
\begin{cases}
\frac{\partial y_1}{\partial t} = \alpha_1 \Delta y_1 + \mu(y_1 + y_2 + y_3) - \beta\, y_1 y_2 - \gamma\, y_1, \\
\frac{\partial y_2}{\partial t} = \alpha_2 \Delta y_2 + \beta\, y_1 y_2 - (\gamma + \delta)\, y_2 - u\, y_2, & (t, r) \in Q, \\
\frac{\partial y_3}{\partial t} = \alpha_3 \Delta y_3 + \delta\, y_2 - \gamma\, y_3 + u\, y_2, \\
\frac{\partial y_1}{\partial \eta} = \frac{\partial y_2}{\partial \eta} = \frac{\partial y_3}{\partial \eta} = 0, & (t, r) \in \Sigma, \\
y_1(0, r) = y_{01},\quad y_2(0, r) = y_{02},\quad y_3(0, r) = y_{03}, & r \in \Omega,
\end{cases}
$$

其中 $u \in U_{ad}$。

**状态方程的存在性与正则性**

**定理1** (i) 设 $\Omega$ 为 $\mathbb{R}^2$ 中的有界区域，其边界光滑性足够好，（假设1）和（假设2）成立，则问题（9）–（10）存在唯一的（全局）强解 $y \in W^{1,2}([0, T])$，使得 $y_1, y_2, y_3 \in E(0, T) \cap L^\infty$。此外，$y_1, y_2, y_3$ 均非负。并且，存在 $C > 0$，独立于 $u$，对所有 $t \in [0, T]$ 有
$$
\left\|\frac{\partial y_i}{\partial t}\right\|_{L^2} + \|y_i\|_{L^2, H^2} + \|y_i\|_{H^1(\Omega)} + \|y_i\|_{L^\infty} \leq C, \quad i = 1, 2, 3.
$$

(ii) 设 $Y = \bar{y} - y$，其中 $y$ 和 $\bar{y}$ 分别表示对应于 $u$ 和 $u + h$ 的满足（14）的（9）的两个解。则

**最优解的存在性**

本节致力于证明最优解的存在性。本节的主要结果如下：

**定理2** 在定理1的相同假设条件下，最优控制问题（8）–（10）存在最优解 $(y^*, u^*, \xi^*)$。

**证明** 设 $J^* = \inf\{J(y, u, \xi)\}$，其中 $0 < \xi \leq T$，$u \in U_{ad}$，$y$ 为（9）的对应解。显然 $J^*$ 为有限值，且可行点集非空，我们可以选取一个可行极小化序列 $(y_n, u_n, \xi_n)$，其中 $\xi_n \in [0, T]$，$u_n \in U_{ad}$，$y_n = (y_{n1}, y_{n2}, y_{n3}) \in W^{1,2}([0, T])$，使得

**可微性**

**定义1** 解映射 $l: u \mapsto y(u)$，从 $L^2$ 到 $F(0, T)$，若在 $L^2$ 上Frechét可微，则对任意 $u \in L^2$，存在 $S(u) \in \mathcal{L}(L^2, F(0, T))$，使得对任意 $u \in L^2$，
$$
\lim_{\|h\|_{L^2} \to 0} \frac{\|l(u+h) - l(u) - S(u)\, h\|_{F(0,T)}}{\|h\|_{L^2}} = 0.
$$
算子 $S(u)$ 称为 $y$ 在 $u$ 处的Frechét导数，记为 $D_u y(u)$。则 $D_u y(u)\, h$ 称为 $y$ 在 $u$ 处沿 $h \in L^2$ 方向的Frechét导数。

我们旨在证明状态系统关于控制变量的Frechét可微性，我们研究如下线性化系统：
$$
\begin{cases}
\frac{\partial z}{\partial t} = \alpha \Delta z + A z + F_y(y)\, z + (G z)\, N
\end{cases}
$$

**必要最优性条件**

在本段中，为了刻画我们的控制策略（以治疗形式）和最优治疗时间（作为自由时间），我们受到以下文献 [34]、[35] 中思想的启发。

设 $p = (p_1, p_2, p_3)$ 为伴随变量，用 $X_T$ 表示算子 $X$ 的伴随。我们可以将系统（8）–（10）的伴随系统写为：
$$
\begin{cases}
-\frac{\partial p}{\partial t} - \alpha \Delta p - A^T p - (F_y(y^*))^T p - (G^T p)\, u = \sigma\, B\, y^*, & \text{on } [0, \xi^*] \\
p(\xi^*, r) = 0
\end{cases}
$$

其中 $u^*$ 为最优控制，$y^* = (y_1^*, y_2^*, y_3^*)$ 为

**二阶最优性条件**

在本节中，我们为本控制问题建立二阶最优性条件。为建立时间最优问题（8）–（9）的条件，我们需要进行变量替换，将其转化为具有固定终端时间的问题。

我们以与文献 [43] 类似的方式进行处理。对任意 $\xi > 0$，定义函数 $\psi: [0, 1] \to [0, \xi]$，设 $\psi(s) = \xi\, s$。令 $(\bar{y}_1, \bar{y}_2, \bar{y}_3)(x, s) = (y_1, y_2, y_3)(x, \psi(s))$，$v(x, s) = u(x, \psi(s))$。

通过对（8）–（9）进行变量替换 $t = \psi(s) = \xi\, s$，我们得到如下最优控制问题：

**数值结果**

在本节中，我们展示数值结果以说明前几节中获得的理论结论。具体而言，我们检验了以治疗形式确定控制策略的方法，并确定最优治疗时间。最优治疗时长 $T^*$ 的确定基于时间最优性判据。在数学上，$T^*$ 被确定为目标泛函对时间的导数为零的点：$\frac{dJ}{dT}\bigg|_{T=T^*} = 0$。

**结论**

我们研究了自由时间最优控制问题，该问题建模了流行病反应扩散模型 SIR 中的治疗时长，该模型描述了在人口密集地区（如口头传播的 COVID-19）的传染病传播。有效治疗传染病需要最小化治疗时间，因为较短的治疗时长与较高的治愈率相关。在本研究中，我们建立了包含自由时间变量的目标函数

**CRediT 作者贡献声明**

Adil El Alami Laaroussi：数据整理、概念构思。Sefian Lamarti Mohammed：撰写初稿。Mostafa Rachik：监督指导。

**致谢**

作者感谢编辑和审稿人提出的宝贵意见和建议，帮助我们改进了本文的研究成果。

Adil El Alami Laaroussi | Sefian Lamarti Mohammed | Mostafa Rachik