李雅芝，劉利利

（1.黔南民族師范學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，貴州 都勻 558000；2.黔南州復(fù)雜系統(tǒng)與智能優(yōu)化重點(diǎn)實(shí)驗(yàn)室，貴州 都勻 558000；3.山西大學(xué) 復(fù)雜系統(tǒng)研究所，太原 030006）

引 言

近些年來(lái)，登革熱在世界各地迅速流行.據(jù)估計(jì)，每年約1 億有癥狀登革熱病例，另有約3 億無(wú)癥狀感染病例，負(fù)擔(dān)最大的是亞洲(75%)[1]，這不僅威脅到了人們的身體健康，而且造成了嚴(yán)重的經(jīng)濟(jì)負(fù)擔(dān)，因此需要尋找有效的防控措施以緩解登革熱在世界范圍內(nèi)的傳播.

目前還沒有針對(duì)登革熱的特定疫苗或抗病毒治療，故媒介控制是防控登革熱病毒傳播的重要策略.這些策略有：通過(guò)環(huán)境治理防止蚊子獲得產(chǎn)卵地、妥善處理固體廢物，去除人為制造的蚊蟲棲息地；對(duì)室外儲(chǔ)水容器施用適當(dāng)殺蟲劑；使用個(gè)人家庭防護(hù)，如紗窗、長(zhǎng)袖衣服、經(jīng)殺蟲劑處理的材料、蚊香和噴霧式殺蟲劑；不育技術(shù)和基因控制等[2].在實(shí)際中，往往會(huì)同時(shí)實(shí)施多種控制措施以達(dá)到更好的控制效果.但是任何一種控制措施都要耗費(fèi)一定的人力、物力和財(cái)力，如何在花費(fèi)盡可能少的前提下達(dá)到最好的控制效果是一個(gè)值得研究的問(wèn)題.

關(guān)于登革熱的最優(yōu)控制已有一些科研工作者進(jìn)行了研究，如Agusto 和Khan[3]考慮了使用殺蟲劑和疫苗的最優(yōu)控制.Prasetyo 等[4]建立了一個(gè)SVIR 模型研究最優(yōu)疫苗控制策略.Xue 等[5]建立了一個(gè)多菌株的登革熱模型并考慮了最優(yōu)疫苗控制.Fister 等[6]研究了投放絕育蚊子和生活環(huán)境改造對(duì)蚊媒傳染病的最優(yōu)控制.Masud 等[7]建立了一個(gè)SIR-SI 模型，加入兩種控制措施，研究其在降低蚊子對(duì)人的叮咬率和蚊子被人感染的有效性，并考慮了最優(yōu)控制.從現(xiàn)有文獻(xiàn)中發(fā)現(xiàn)，在模型中加入基因控制研究最優(yōu)控制策略的研究還較少，關(guān)于基因控制目前主要是投放攜帶Wolbachia 的蚊子.有實(shí)驗(yàn)發(fā)現(xiàn)Wolbachia 可以抑制登革熱病毒在蚊子體內(nèi)的復(fù)制，所以讓蚊子攜帶Wolbachia 有助于抑制登革熱的傳播.而只要成年雌蚊感染了Wolbachia，其后代也會(huì)感染此菌，并且被感染的雄性與正常雌性的后代會(huì)因?yàn)榧?xì)胞質(zhì)不相容(CI)而死亡[8]，這些特性都有助于Wolbachia 在蚊子種群的傳播.相比傳統(tǒng)的控制措施，使用Wolbachia 是一種可持續(xù)且環(huán)境友好的控制方式.

為了研究登革熱的最優(yōu)控制策略，本文將建立登革熱在蚊子和人之間的傳播模型，加入自我保護(hù)、殺蟲劑和Wolbacia 三種控制措施，研究最優(yōu)綜合控制策略.

1 模型建立

關(guān)于登革熱傳播模型的建立，有以下幾點(diǎn)解釋：

(Ⅰ)對(duì)于人采用SIS 倉(cāng)室模型，因?yàn)榈歉餆岬臐摲谙鄬?duì)人的壽命是很短的，在此忽略.而登革熱有四種類型，人感染了一種類型的登革熱并恢復(fù)后，仍可能感染其他類型的登革熱，所以SIS 模型是合理的.

(Ⅱ)對(duì)于蚊子采用SEI 倉(cāng)室模型，因?yàn)榈歉餆岵《驹谖米芋w內(nèi)的潛伏期(約7 d)相對(duì)蚊子的整個(gè)生命周期(30～100 d)是不能被忽略的.

(Ⅲ)考慮Wolbachia 在蚊子種群中的垂直傳播和CI 影響.

(Ⅳ)控制變量u1(t)表示通過(guò)使用蚊帳、穿長(zhǎng)袖衣服等保護(hù)人類免受蚊蟲叮咬的效果，控制變量u2(t)表示為減少蚊子數(shù)量而進(jìn)行的捕殺工作.

綜合以上幾種因素，可得到如下模型：

其中Nmi=Smi+Emi+Imi,Nmu=Smu+Emu+Imu,Nm=Nmi+Nmu,Nh=Sh+Ih.各狀態(tài)變量和參數(shù)的含義分別見表1和表2，“W”表示“Wolbachia”，“D”表示“登革熱”.由于有Wolbachia 菌株的垂直傳播率和CI 發(fā)生可能性均接近1，故在模型(1)中，將兩種影響的可能性視為1，即只要成年雌蚊感染了Wolbachia，其后代也一定會(huì)感染，感染雄蚊和正常雌蚊的后代一定會(huì)因?yàn)镃I 影響而死亡.模型(1)第二個(gè)方程中的項(xiàng)表示正常雌性產(chǎn)生正常后代.另外，由于對(duì)蚊子的捕殺不可能超過(guò)蚊子的出生，否則蚊子種群將會(huì)由于捕殺而滅絕，故假設(shè)bm>u2.

表1 模型(1)中各狀態(tài)變量的含義Table 1 The meanings of variables in model (1)

表2 模型(1)中各參數(shù)的含義Table 2 The meanings of parameters in model (1)

定義

關(guān)于模型(1)的適定性，有如下結(jié)果.

定理1模型(1)的解非負(fù)且在D內(nèi)是有界的.進(jìn)一步地，D是模型(1)的正不變集.

證明由于系統(tǒng)(1)的初值非負(fù)，易知模型(1)的解非負(fù)，故只需證明Nm和Nh是有界的.

分別將系統(tǒng)(1)中的前六個(gè)方程和后兩個(gè)方程相加，可以得到

積分可得

得證.

2 常數(shù)控制分析

令u1(t)=u1,u2(t)=u2，其中u1,u2為常數(shù).下考慮常數(shù)控制對(duì)模型(1)基本再生數(shù)的影響.

在系統(tǒng)(1)中令Emi=0,Imi=0,Emu=0,Imu=0,Ih=0，則得到如下子系統(tǒng)：

定義g(0,0)=0，則系統(tǒng)(2)被連續(xù)延拓到 (0,0).分別令f(Smi,Smu)=0,g(Smi,Smu)=0,v(Sh)=0，通過(guò)計(jì)算模型(2)的平衡態(tài)可知模型(1)最多存在三個(gè)無(wú)病平衡態(tài) (,,0,0,0,0,0)，分別為

由解對(duì)初值的連續(xù)依賴性可知E00是不穩(wěn)定的.記E01和E02對(duì)應(yīng)的基本再生數(shù)分別為R01和R02，使用下一代矩陣方法可以計(jì)算得到

通過(guò)計(jì)算系統(tǒng)(1)在E0i(i=1,2)處 的Jacobi 矩陣對(duì)應(yīng)特征多項(xiàng)式的特征根，可以得到：當(dāng)R01<1時(shí)，E01是局部漸近穩(wěn)定的；當(dāng)R02<1時(shí)，E02是局部漸近穩(wěn)定的.

比較R01和R02的表達(dá)式，其區(qū)別在于參數(shù)和 βh，由表2可知，<βh，故有R01

下面分析R0i(i=1,2)對(duì) 控制參數(shù)u1和u2的 敏感性，以R01為例.

在此采用Chitnis 等[9]定義的標(biāo)準(zhǔn)正向靈敏度指數(shù)衡量變量P對(duì)參數(shù)u的敏感性.結(jié)合R01的表達(dá)式，通過(guò)計(jì)算可以得到

可以看出，R01與u1和u2均為負(fù)相關(guān)的，即增大常數(shù)控制參數(shù)u1和u2，均會(huì)使得R0i(i=1,2)減小，進(jìn)而減少登革熱的感染人數(shù).又? ||/?u1=1/(1?u1)2>0,?||/?u2=bm/2(bm?u2)2>0，所以u(píng)1和u2越 大，R01與它們的相關(guān)性越強(qiáng).

圖1為當(dāng)R0i(i=1,2)>1時(shí)，控制參數(shù)u1和u2對(duì) 登革熱感染人數(shù)的影響.參數(shù)取值為：bm=0.8,dm=4E?2,βm=0.001 5,em=0.05,w=2,dh=4E?5,γh=0.002,βh=0.000 1,q3=0.9.從圖中可以看出：控制參數(shù)取值越大，登革熱感染人數(shù)越少，且病例數(shù)呈顯著上升趨勢(shì)所需的時(shí)間越長(zhǎng).

圖1 常數(shù)控制對(duì)登革熱感染人數(shù)的影響Fig.1 Effects of constant control on the number of dengue infections

3 最優(yōu)控制分析

由于Wolbachia 主要是通過(guò)抑制登革熱病毒在蚊子體內(nèi)的復(fù)制，進(jìn)而達(dá)到抑制登革熱傳播的目的，對(duì)于蚊子種群屬于內(nèi)部控制，控制效果與最初的投放量有關(guān)，故在此僅研究自我保護(hù)和殺蟲劑兩種外部控制措施的最優(yōu)綜合控制策略.

首先建立如下目標(biāo)泛函：

其中(u1(t),u2(t))∈K，這里的K為控制集，定義為A表示由攜帶登革熱病毒的病人造成的人均損失，B1和B2分別表示人的自我保護(hù)和滅蚊工作相關(guān)的成本，T表示控制的終端時(shí)間.假設(shè)時(shí)間區(qū)間為 [0,T]，我們需要尋找最優(yōu)的u1(t)和u2(t)使 得目標(biāo)泛函J(u1(t),u2(t))達(dá)到最小.

為了研究系統(tǒng)(1)在目標(biāo)泛函(3)下的最優(yōu)控制策略，現(xiàn)在給出最優(yōu)控制的存在性結(jié)論如下.

定理2對(duì)于模型(1)，存在最優(yōu)控制(t),(t)，使得

證明易知目標(biāo)函數(shù)J(u1(t),u2(t))中 的被積函數(shù)關(guān)于u1(t)和u2(t)是凸函數(shù).由定理1 知系統(tǒng)(1)的解是有界的，故系統(tǒng)關(guān)于變量Smi(t),Smu(t),Emi(t),Emu(t),Imi(t),Imu(t),Sh(t)，Ih(t)滿足Lipschitz 性質(zhì)，進(jìn)而存在最優(yōu)的二元組((t),(t)).得證.

關(guān)于動(dòng)力系統(tǒng)的最優(yōu)控制方法很多，動(dòng)力系統(tǒng)類型不同，求解最優(yōu)控制的方法也不同[10-11].模型(1)為一個(gè)高維ODE 系統(tǒng)，可以采用Pontryagin 極值原理[12]尋找具體的最優(yōu)控制((t),(t)).

首先構(gòu)造Hamilton 函數(shù)H如下：

其中λi(t)(i=1,2,···,8)稱為協(xié)態(tài)變量，滿足下面的伴隨方程：

橫截條件為

利用Hamilton 函數(shù)，我們將狀態(tài)方程(1)、目標(biāo)泛函(3)和伴隨方程(4)聯(lián)立起來(lái)構(gòu)成一個(gè)最優(yōu)控制系統(tǒng).由最優(yōu)條件可得

綜上，關(guān)于最優(yōu)控制的具體形式有如下結(jié)果.

定理3在控制集K上使得J(u1(t),u2(t))取 得最小值的最優(yōu)控制((t),(t))為

4 數(shù)值模擬

本部分采用四階Runge-Kutta 前推回代法計(jì)算最優(yōu)控制和狀態(tài)值[13].記狀態(tài)變量向量為

x=(Smi,Smu,Emi,Emu,Imi,Imu,Sh,Ih),

協(xié)態(tài)變量向量為

λ=(λ1,λ2,λ3,λ4,λ5,λ6,λ7,λ8),

控制變量向量為

u=(u1,u2).

其基本算法步驟如下：

步1 在時(shí)間范圍[0,T]內(nèi)對(duì)u進(jìn)行初步估計(jì).

步2 使用初始條件和估計(jì)的u值，依據(jù)系統(tǒng)(1)關(guān)于時(shí)間向前求解x.

步3 使用橫截條件(5)和u，x的值, 依據(jù)系統(tǒng)(4)關(guān)于時(shí)間向后求解λ.

步4 通過(guò)在最優(yōu)控制的表征(6)中輸入新的x和λ 值來(lái)更新u.

步5 檢查收斂性.重復(fù)前面的步驟，直到當(dāng)前的狀態(tài)值、伴隨值和控制值充分收斂.

取參數(shù)值為：bm=0.8,dm=4E?5,βm=0.001 5,em=0.05,w=2,dh=4E?5,γh=0.002,βh=0.001,q3=0.9.本組參數(shù)值使得無(wú)控制系統(tǒng)的基本再生數(shù)為232.5204，也即登革熱會(huì)在人群中流行.控制的終端時(shí)間T=30.初值取為：Smi(0)=1 000,Smu=1 000,Emi=100,Emu=100,Imi=20,Imu=20,Sh=100,Ih=5.

控制策略對(duì)登革熱感染人數(shù)的影響見圖2，權(quán)重取值為：A=1，B1=1，B2=1.可以看出，通過(guò)采取控制措施，在一個(gè)月內(nèi)被感染人數(shù)有了顯著降低，且最優(yōu)綜合控制策略的效果最好，蚊子捕殺的控制效果優(yōu)于人的自我保護(hù)的控制效果.這充分說(shuō)明直接對(duì)傳播媒介本身采取例如捕殺等控制措施是非常有必要的.

圖2 控制策略對(duì)登革熱感染人數(shù)的影響Fig.2 Effects of the control strategy on the number of dengue infections

最優(yōu)綜合控制與單一控制的比較見圖3，權(quán)重取值為：A=1,B1=1,B2=1.由圖3(a)可知，對(duì)于人的自我保護(hù)，如果同時(shí)還對(duì)蚊子進(jìn)行捕殺，在前2 d，人的自我保護(hù)在控制下界，之后大約需要在上界持續(xù)12 d 左右，隨后逐漸下降到控制下界；如果只有人的自我保護(hù)，沒有捕殺蚊子，可以發(fā)現(xiàn)幾乎在整個(gè)模擬周期內(nèi)人的自我保護(hù)都要維持在上界.由圖3(b)可知，對(duì)于蚊子的捕殺，如果同時(shí)還有人的自我保護(hù)，對(duì)蚊子的捕殺最開始就要在控制上界，而且要持續(xù)24 d 左右，之后再緩慢下降到控制下界；如果只捕殺蚊子，沒有人的自我保護(hù)，對(duì)蚊子的捕殺要在控制上界持續(xù)27 d 左右.總之，相比綜合控制策略，單一控制策略在上界的時(shí)間要長(zhǎng).

圖3 最優(yōu)綜合控制與單一控制的比較：(a)兩種控制策略下u 1的 變化圖；(b)兩種控制策略下u 2的變化圖Fig.3 Comparison between the optimal integrated control and the single control: (a)the variation diagram of u 1 under 2 control strategies;(b)the variation diagram of u 2 under 2 control strategies

不同權(quán)重對(duì)控制變量的影響見圖4.B1和B2取值越大表示對(duì)應(yīng)控制的成本越高，A取值越大表示感染登革熱造成的損失越大.從圖4(a)、(b)可以發(fā)現(xiàn)，隨著控制成本增大，兩種控制措施的控制力度都在減小.從圖4(c)、(d)可以發(fā)現(xiàn)，隨著疾病造成的損失增大，人的自我保護(hù)力度要增大到最大.

圖4 不同權(quán)重對(duì)控制變量的影響：(a)改變 B1對(duì)u 1的 影響，A =1,B2=1； (b)改變 B2對(duì)u 2的 影響，A =1,B1=1；(c)A=1 時(shí)，u 1,u2 的變化圖，B 1=1,B2=1； (d)A =5 時(shí)，u 1,u2 的 變化圖，B1=1,B2=1Fig.4 The influences of different weights on control variables: (a)the influence of B1 on u 1, A =1,B2=1; (b)the influence of B2 on u 2 , A =1,B1=1;(c)the variation diagram of u 1 and u2, A =1,B1=1,B2=1; (d)the variation diagram of u 1 and u2, A =5 , B1=1,B2=1

5 結(jié) 語(yǔ)

本文建立了一個(gè)SEI-SIS 模型描述登革熱在人和蚊子之間的傳播，加入了Wolbachia、人的自我保護(hù)和蚊子捕殺三種控制措施，研究最優(yōu)綜合控制策略，并采用數(shù)值模擬探討控制策略對(duì)登革熱感染人數(shù)的影響，對(duì)比綜合控制策略與單一控制策略，對(duì)比權(quán)重對(duì)控制策略的影響.在實(shí)際中，登革熱的控制方法有多種，可以根據(jù)實(shí)際情況，對(duì)不同的控制方法進(jìn)行組合來(lái)研究最優(yōu)綜合控制策略.