基于數(shù)據(jù)驅(qū)動(dòng)的SEIQRW新冠肺炎模型及最優(yōu)控制

苗 慧，李淑萍

(中北大學(xué) 理學(xué)院，山西 太原 030051)

0 引 言

從2019年起，由嚴(yán)重急性呼吸綜合征冠狀病毒2(SARS-CoV-2)引起的新冠肺炎在世界各地大規(guī)模流行.據(jù)世界衛(wèi)生組織WHO報(bào)告數(shù)據(jù)來(lái)看，截至2022年1月17日,全世界累計(jì)確診329 109 769例，死亡555 3630例.全球新冠肺炎病例不斷增加，全面了解它的傳播途徑，采取對(duì)應(yīng)的防控措施非常重要.研究表明，呼吸道飛沫和直接接觸是新冠病毒傳播的主要途徑[1-2]，且新冠病毒在物體表面也能存活數(shù)天[3]，因此，人類(lèi)通過(guò)間接接觸被污染的物體表面也可能會(huì)導(dǎo)致感染[4-5].

專(zhuān)家學(xué)者對(duì)不同地區(qū)新冠肺炎的傳播進(jìn)行了數(shù)學(xué)建模.Wang Liyan等在揚(yáng)州疫情的基礎(chǔ)上考慮了潛伏期和治愈期所服從的新密度函數(shù)，以及感染模型中的疫苗接種效應(yīng)，研究了長(zhǎng)潛伏期新冠肺炎的傳播機(jī)制，認(rèn)為通過(guò)足夠的隔離和高接種疫苗措施可促使疫情消失[6].Kolebaje O T等提出了一種結(jié)合非藥物干預(yù)措施的新冠肺炎SEIIAQHR數(shù)學(xué)模型，其中S代表易感者，E代表潛伏者，I代表感染者，IA代表無(wú)癥狀感染者，Q代表隔離檢疫者，H代表住院者，R代表恢復(fù)者，并將此模型用于研究一些非洲國(guó)家的累計(jì)確診病例數(shù)的時(shí)間序列演變[7].Shen Zhonghua等考慮了巴基斯坦的新冠肺炎報(bào)告病例，建立了帶有疫苗接種V倉(cāng)室的SVEAIR模型(A為無(wú)癥狀感染倉(cāng)室)，并進(jìn)行最優(yōu)控制分析[8].Abioye A I等建立了帶有隔離倉(cāng)室的SEIQR模型，并結(jié)合戴口罩、積極篩查檢測(cè)、控制社交距離等嚴(yán)格的控制策略，預(yù)測(cè)尼日利亞的新冠肺炎將被根除[9].

上述模型有的考慮了隔離檢疫，有的考慮了疫苗接種，但沒(méi)有考慮環(huán)境病毒對(duì)疫情傳播的影響.自新冠疫情爆發(fā)以來(lái)，我國(guó)政府采取了一系列隔離防控措施，如通過(guò)建立火神山醫(yī)院集中隔離收治患者[10]，并取得了很大的成效.因此，本文建立了一個(gè)隔離倉(cāng)室Q的SEIQR模型.另外，近期的深圳疫情很可能是由國(guó)外入境的新冠病毒污染物品引起的[11]，本文將研究環(huán)境病毒W(wǎng)對(duì)新冠肺炎傳播的影響，以期為政府部門(mén)制定防控策略提供理論依據(jù).

1 模型的建立

S(t)為t時(shí)刻易感者數(shù)量，E(t)為t時(shí)刻潛伏者數(shù)量，I(t)為t時(shí)刻感染者數(shù)量，Q(t)為t時(shí)刻隔離檢疫者數(shù)量，R(t)為t時(shí)刻恢復(fù)者數(shù)量，W(t) 為t時(shí)刻環(huán)境中的新冠病毒數(shù)量，模型倉(cāng)室圖如圖1 所示.圖1 中，Λ為總?cè)丝诘妮斎肼?，?為潛伏者對(duì)易感者的傳染率，β2為染病者對(duì)易感者的傳染率，β3為環(huán)境中的新冠病毒對(duì)易感者的傳染率，μ為自然死亡率，k為潛伏者到染病者的轉(zhuǎn)化率，b為染病者隔離檢疫的比例系數(shù)，α為染病者的因病死亡率，γ為染病者的恢復(fù)率，δ為隔離后的康復(fù)率，m為潛伏者向外界環(huán)境釋放病毒的速率，μ1為環(huán)境中新冠病毒的自然死亡率.模型為

圖1 SEIQRW流程圖Fig.1 Flow chart of SEIQRW transmission

(1)

由N=S+E+I+Q+R得

(2)

Ω={(S,E,I,Q,R,W)|S,E,I,Q,R,W≥0,

(3)

2 平衡點(diǎn)和基本再生數(shù)的計(jì)算

利用下一代矩陣法[12]，系統(tǒng)(1)中受感染的相關(guān)倉(cāng)室可寫(xiě)為

式中：X=(E,I,Q,W)T；F為新感染的矩陣；V為倉(cāng)室個(gè)體轉(zhuǎn)移矩陣，則

在無(wú)病平衡點(diǎn)E0處對(duì)(E,I,Q,W)求偏導(dǎo)，可得

計(jì)算得

因此，基本再生數(shù)為

為求解正平衡點(diǎn)，令模型(1)右端為0，則

(4)

當(dāng)R0>1時(shí)，正平衡點(diǎn)為E1=(S*,E*,I*,Q*,R*,W*).其中，

3 平衡點(diǎn)的穩(wěn)定性分析

3.1 無(wú)病平衡點(diǎn)的穩(wěn)定性

定理 1當(dāng)R0<1時(shí)，系統(tǒng)(1)的無(wú)病平衡點(diǎn)E0在Ω內(nèi)是局部漸近穩(wěn)定的；當(dāng)R0>1時(shí),無(wú)病平衡點(diǎn)不穩(wěn)定.

證明系統(tǒng)(1)在E0的Jacobian矩陣為

J|E0的特征方程為

(λ+μ)2(λ+μ+α+δ)(λ3+a1λ2+

a2λ+a3)=0.

(5)

顯而易見(jiàn)，λ1,λ2=-μ,λ3=-(μ+α+δ).

其中，

μ1+b+μ+α+γ>0,

a2=(μ+k)(b+μ+α+γ)+μ1(μ+k)+

μ1(b+μ+α+γ)>0,

a3=μ1(μ+k)(b+μ+α+γ)(1-R0)>0,

μ1(b+μ+α+γ)2>0.

由Hurwitz判據(jù)[13]可知，無(wú)病平衡點(diǎn)E0在可行域內(nèi)局部漸近穩(wěn)定.

定理 2當(dāng)R0<1時(shí)，系統(tǒng)(1)的無(wú)病平衡點(diǎn)E0在Ω內(nèi)是全局漸近穩(wěn)定的.

證明可以將原系統(tǒng)(1)的染病倉(cāng)室重新表示為

x′=(F-V)x-f(x,y),

(6)

式中：x=(E,I,Q,W)T為染病倉(cāng)室中的人群；y=(S,R)T為非染病倉(cāng)室中的人群，f(x,y)=(0,0,0,0)T.

運(yùn)用矩陣?yán)碚撝械腜erron特征向量[14]，非負(fù)矩陣V-1F的左特征向量ωT對(duì)應(yīng)的特征值為ρ(V-1F)=ρ(FV-1)=R0，可求得ωT=(β1,β2,0,β3).若f(x,y)≥0,F≥0,V-1≥0，則定義Lyapunov函數(shù)為

Q=ωTV-1x=

當(dāng)R0<1時(shí)，

Q′=ωTV-1x′=

(R0-1)ωTx-ωTV-1f(x,y)=

(R0-1)(β1E+β2I+β3W)≤0.

所以，當(dāng)R0<1時(shí)，無(wú)病平衡點(diǎn)X0在Ω內(nèi)全局漸近穩(wěn)定.

3.2 正平衡點(diǎn)的穩(wěn)定性

b1=A+B+C+μ+μ1,

b2=A(B+C)+(A+C)μ1+(B+C)μ+μμ1+

b3=ABC+(AB+AC)μ1+Cμμ1+Aμ1β1S*+

b4=ABCμ1+ACμ1β1S*,

A=β1E*+β2I*+β3W*,

B=μ+k-β1S*=[kμ1β2(μ+k)+

mβ3C(μ+k)]/Ζ,

C=b+μ+α+γ,

Z=μ1β1C+kμ1β2+mβ3C.

證明系統(tǒng)(1)在E1處的Jacobian矩陣為

J|E1的特征方程為

(λ+μ)(λ+μ+α+δ)(λ4+b1λ3+b2λ2+

b3λ+b4)=0.

(7)

顯然，λ1=-μ,λ2=-(μ+α+δ).b1,b2,b3,b4>0，且有

b1b2-b3=(A+B)a2+ABμ+ACμ+AC2+

μ1(A+C)(C+μ+μ1)+μ(B+C)(C+μ+

μ1)+μμ1(μ+μ1)+Aμβ1S*+(C+μ1)·

4 最優(yōu)控制策略

為找到合理有效的控制策略，本文采用最優(yōu)控制理論[15]，并運(yùn)用龐特里亞金極大值原理[16]找到模型中最優(yōu)控制所需要的條件.

在系統(tǒng)(1)中引入3個(gè)控制項(xiàng)u1,u2,u3，將控制系統(tǒng)表示為

(8)

式中：u1為減少易感者與新冠病毒的直接接觸和間接接觸的隔離策略；u2為增強(qiáng)感染者的治療恢復(fù)率的治療策略；u3為消滅環(huán)境中病毒的環(huán)境衛(wèi)生策略.

控制函數(shù)集為

U={U(·)∈L1([0,T]；R3)|0

umax<1,?t∈[0,T]}.

為考慮感染者的治療成本和控制成本，將目標(biāo)函數(shù)定義為

(9)

龐特里亞金極大值原理是將式(8)和(9)轉(zhuǎn)化為一個(gè)關(guān)于控制變量的最小化的哈密頓問(wèn)題，構(gòu)造哈密頓函數(shù)，即

(10)

式中：λi(i=1,2,…,6)為伴隨變量.

(11)

那么，最優(yōu)控制為

證明由伴隨系統(tǒng)(11)得

根據(jù)最優(yōu)化條件

得最優(yōu)控制值為

5 數(shù)值模擬

5.1 平衡點(diǎn)的穩(wěn)定性

假設(shè)Λ=100,β1=0.000 2,β2=0.000 1,β3=0.000 05,μ=0.069,μ1=0.06,γ=0.04,k=0.69,b=0.7,α=0.006 6,m=0.009,δ=0.83，可得R0=0.56<1，如圖2(a)所示，驗(yàn)證了無(wú)病平衡點(diǎn)是全局漸近穩(wěn)定的.

(a) R0<1

假設(shè)β1=0.000 9,β2=0.000 2,β3=0.000 05,m=0.09,b=0.07,δ=0.63，可得R0=3.28>1，如圖2(b)所示，故地方病平衡點(diǎn)是全局漸近穩(wěn)定的.

5.2 數(shù)據(jù)擬合

基于政府在新冠肺炎疫情前期采取隔離防控措施時(shí)感染者的數(shù)據(jù)以及部分參數(shù)估計(jì)[17]，得到k=1/5.2,γ=1/7,α=0.043,δ=1/10,并且假設(shè)Λ=3 000,β1=0.65×10-5,β2=0.85×10-5,β3=0.59×10-6,μ=0.06,μ1=0.056,m=0.89,b=0.07.運(yùn)用Matlab模擬對(duì)比了有無(wú)考慮環(huán)境病毒時(shí)估計(jì)的染病者人數(shù)曲線，如圖3 所示，發(fā)現(xiàn)帶病毒項(xiàng)的系統(tǒng)與實(shí)際感染人數(shù)更相符.

圖3 帶病毒、不帶病毒的模型與患病人數(shù)擬合Fig.3 Fitting model with and without virus to number of sick people

5.3 敏感性分析

運(yùn)用PRCC法進(jìn)行敏感性分析，可以更好地觀察各個(gè)參數(shù)與基本再生數(shù)R0的相關(guān)性，如圖4 所示.可知參數(shù)Λ,β1,β2,β3,m與R0呈正相關(guān)，μ,μ1,k,b,α,γ與R0呈負(fù)相關(guān).這也說(shuō)明了直接和間接接觸都加速了病毒傳播，而控制接觸和消滅環(huán)境病毒可降低基本再生數(shù)，從而減少感染人數(shù).

圖4 R0與參數(shù)的相關(guān)性Fig.4 Correlation between R0 and parameters

5.4 最優(yōu)控制

利用前后掃描法和四階龍格庫(kù)塔公式進(jìn)行模擬，為控制成本，節(jié)約醫(yī)療資源，令u1=0.8,u2=0.7,u3=0.6，得到最優(yōu)控制函數(shù)圖，如圖5(a) 所示.令Λ=100,β1=0.000 2,β2=0.000 1,β3=0.000 05,μ=0.069,μ1=0.06,m=0.009,k=0.69,b=0.7,α=0.006 6,γ=0.04,δ=0.83，計(jì)算得R0=0.56<1.如圖5(b)和(c)，當(dāng)加入控制時(shí)，感染人數(shù)和環(huán)境中病毒數(shù)量在較短時(shí)間內(nèi)迅速減少為0，由此可知，最優(yōu)控制策略可以阻斷新冠肺炎的傳播.

(a) 控制變量u(t)與t的函數(shù)關(guān)系

6 結(jié) 論

本文建立了一個(gè)包括隔離和環(huán)境中新冠病毒倉(cāng)室的模型，通過(guò)動(dòng)力學(xué)分析，得到基本再生數(shù)R0、無(wú)病平衡點(diǎn)和正平衡點(diǎn)的穩(wěn)定性及最優(yōu)控制策略.基于數(shù)值模擬對(duì)理論分析進(jìn)行了驗(yàn)證，參數(shù)的敏感性分析和數(shù)據(jù)擬合進(jìn)一步證實(shí)隔離是疫情防控必不可少的措施.將環(huán)境病毒考慮在內(nèi)可使得模型曲線與患病人數(shù)更加接近，強(qiáng)調(diào)了環(huán)境病毒在疫情傳播過(guò)程是不可忽視的.如果不在模型中納入環(huán)境間接傳播，會(huì)導(dǎo)致低估疫情的嚴(yán)重程度.研究結(jié)果為世衛(wèi)組織建議的洗手和表面清潔和消毒的做法提供了強(qiáng)有力的支撐，也為我國(guó)采取及時(shí)隔離的政策提供了理論依據(jù).