具有跨物種多級(jí)傳播特征的包蟲(chóng)病優(yōu)化算法*

黃光球，陸秋琴

西安建筑科技大學(xué) 管理學(xué)院，西安710055

1 引言

工程中存在一類非線性優(yōu)化問(wèn)題，此類問(wèn)題的數(shù)學(xué)表達(dá)式中存在不可導(dǎo)的數(shù)學(xué)表達(dá)式?；谔荻鹊膬?yōu)化方法無(wú)法求解此類問(wèn)題。為了求解此類問(wèn)題，人們提出了啟發(fā)式搜索方法[1]，群智能優(yōu)化算法[2]就是其中的一種。因群智能優(yōu)化算法對(duì)工程優(yōu)化問(wèn)題沒(méi)有特殊要求，故具有廣泛的適應(yīng)性。近幾年，群智能優(yōu)化算法得到快速發(fā)展。

群智能優(yōu)化算法是依據(jù)特定的生物進(jìn)化場(chǎng)景而構(gòu)建出來(lái)的，其算子依據(jù)個(gè)體間的相互作用而構(gòu)建，其邏輯結(jié)構(gòu)是依據(jù)生物進(jìn)化場(chǎng)景的內(nèi)涵而構(gòu)建[3]。早期的群智能算法典型代表有粒子群算法[4-7]、蟻群算法[8-9]、蜂群算法[10]、布谷鳥(niǎo)算法[11]、蝙蝠算法[12]、遺傳算法[13]等。目前，很多新型群智能算法先后被提出，如細(xì)菌覓食優(yōu)化、蛙跳算法、猴群算法、螢火蟲(chóng)算法、鴿群算法、果蠅算法和頭腦風(fēng)暴算法等[14]。然而，目前所提出的群智能優(yōu)化算法的共同缺陷包括如下幾方面：（1）算法所依賴生物進(jìn)化場(chǎng)景均非常簡(jiǎn)單，因而依據(jù)此類場(chǎng)景設(shè)計(jì)出來(lái)的群智能算法的邏輯結(jié)構(gòu)非常簡(jiǎn)單；（2）所能開(kāi)發(fā)出的算子非常少，從而導(dǎo)致個(gè)體之間的信息交換不充分，搜索易限于局部陷阱；（3）算法的全局收斂性難以保證。

隨著互聯(lián)網(wǎng)的高速發(fā)展，社會(huì)已進(jìn)入大數(shù)據(jù)時(shí)代，人們所遇到的問(wèn)題變得越來(lái)越復(fù)雜，故需要提出高智能性的算法來(lái)解決這些復(fù)雜問(wèn)題。對(duì)于群智能算法來(lái)說(shuō)，如何增加群智能算法的智能特征，顯然是需要解決的關(guān)鍵問(wèn)題。解決該問(wèn)題的一種策略是：精心選擇出某種具有特殊性質(zhì)的生物進(jìn)化場(chǎng)景，巧妙解決目前群智能算法所存在問(wèn)題，即可設(shè)計(jì)出高智能性的群智能算法。

若一類生物進(jìn)化場(chǎng)景具有豐富的演化內(nèi)涵，其中個(gè)體之間的相互作用關(guān)系十分豐富，則依據(jù)此類生物進(jìn)化場(chǎng)景就可以設(shè)計(jì)出具有算子眾多、邏輯結(jié)構(gòu)清晰的高智能性群智能優(yōu)化算法。能夠跨三個(gè)物種實(shí)施多級(jí)傳播的包蟲(chóng)傳染病攻擊家犬、牛羊和人類的過(guò)程，就是這樣一個(gè)場(chǎng)景。

包蟲(chóng)病是棘球蚴寄生在人體所致的一種人獸共患寄生蟲(chóng)病。我國(guó)有囊型和泡型包蟲(chóng)病兩種，分別由棘球蚴和泡球蚴寄生在人體組織器官中所致。囊型包蟲(chóng)病呈世界性分布，畜牧業(yè)發(fā)達(dá)的國(guó)家和地區(qū)多見(jiàn)。我國(guó)包蟲(chóng)病高發(fā)流行區(qū)主要集中在高山草甸地區(qū)及氣候寒冷、干旱少雨的牧區(qū)及半農(nóng)半牧區(qū)，以新疆、青海、甘肅、寧夏、西藏、內(nèi)蒙古、陜西、河北、山西和四川等地較為嚴(yán)重。包蟲(chóng)來(lái)源于動(dòng)物的排泄物，沒(méi)有什么有效的治療方式，在西藏被視為西藏第一癌癥。

由于包蟲(chóng)病可同時(shí)危害人類和家畜養(yǎng)殖業(yè)，近年來(lái)全世界研究人員對(duì)該傳染病給予了高度關(guān)注，主要研究成果表現(xiàn)為如下幾方面：

（1）包蟲(chóng)病流行病學(xué)分析。為了了解包蟲(chóng)病的流行病學(xué)特征，通過(guò)對(duì)包蟲(chóng)病病例進(jìn)行流行病學(xué)特征調(diào)查研究，探討病例可能的感染來(lái)源，為制定包蟲(chóng)病防控措施提供依據(jù)[15]。

（2）包蟲(chóng)病免疫學(xué)研究。通過(guò)對(duì)家畜和人的包蟲(chóng)病免疫機(jī)理進(jìn)行研究，為包蟲(chóng)病的防治和新型疫苗的生產(chǎn)提供科學(xué)參考[16]。

（3）包蟲(chóng)病診斷方法研究。提出包蟲(chóng)病的診斷方法，為包蟲(chóng)病的治療提供科學(xué)、快速和準(zhǔn)確的依據(jù)[17]。

（4）包蟲(chóng)病的治療方法研究。從手術(shù)和藥物等角度探討包蟲(chóng)病的各種治療方法[18]。

包蟲(chóng)病具有跨多物種傳播特征，能夠?qū)θ祟愒斐珊艽笸纯啵谶@悲慘的自然現(xiàn)象后面，卻隱藏著一個(gè)性能優(yōu)良的群智能優(yōu)化算法，本文稱其為包蟲(chóng)傳染病優(yōu)化算法（hydatid disease optimization，HDO）。本文著重解決如下五個(gè)問(wèn)題：

（1）如何將包蟲(chóng)傳染病模型轉(zhuǎn)化為能求解復(fù)雜優(yōu)化問(wèn)題的HDO 算法。

（2）如何使得HDO 算法中的算子能充分反映動(dòng)物個(gè)體之間的相互作用關(guān)系，以便體現(xiàn)包蟲(chóng)傳染病模型的思想。

（3）如何證明HDO 算法的全局收斂性。

（4）如何確定HDO 算法的最佳參數(shù)設(shè)置。

（5）如何用HDO 算法求解關(guān)聯(lián)區(qū)域VOCs（volatile organic compounds）聯(lián)防聯(lián)控最優(yōu)減排優(yōu)化問(wèn)題。

2 包蟲(chóng)傳染病模型的優(yōu)化算法設(shè)計(jì)

設(shè)要求解的優(yōu)化問(wèn)題為：

式中，Rn是n維歐氏空間；X=(x1,x2,…,xn)是一個(gè)n維決策向量；H為搜索空間；F(X)為目標(biāo)函數(shù)；gi(X)≥0 為第i個(gè)不等式約束條件，I為不等式約束條件個(gè)數(shù)。目標(biāo)函數(shù)F(X)和約束條件gi(X)沒(méi)有限制。

2.1 算法場(chǎng)景設(shè)計(jì)

在一個(gè)草原牧區(qū)生活有一群牧民，其數(shù)量有N3名，其編號(hào)分別是1,2,…,N3；牧民們以飼養(yǎng)藏羊?yàn)樯?，藏羊的?shù)量有N2只，其編號(hào)分別是1,2,…,N2；與此同時(shí)，牧民們還飼養(yǎng)了一群狗用于控制和保護(hù)羊群，狗的數(shù)量有N1只，其編號(hào)分別是1,2,…,N1。一條狗、一只羊或一名牧民統(tǒng)稱為一個(gè)個(gè)體，每個(gè)個(gè)體均由n個(gè)特征（器官）來(lái)表征。令u表示個(gè)體類型，u=1表示狗類，u=2 表示羊類，u=3 表示牧民類，即對(duì)類型為u的個(gè)體i來(lái)說(shuō)，其表征特征為。

該草原牧區(qū)有一種稱為棘球絳蟲(chóng)的寄生蟲(chóng)病在狗群中流行。棘球絳蟲(chóng)寄生于狗的小腸中，蟲(chóng)卵隨狗的糞便排出。由于狗隨羊群到處游動(dòng)，其排泄的糞便廣泛地污染了牧區(qū)的水源、土壤和草場(chǎng)。棘球絳蟲(chóng)蟲(chóng)卵在自然界有非常強(qiáng)的適應(yīng)能力，它能在自然狀態(tài)下保持持續(xù)感染力。羊在吃草、喝水的過(guò)程中，就會(huì)把附著在草葉上或水中的蟲(chóng)卵一同吃入體內(nèi)。生活在羊胃中的蟲(chóng)卵會(huì)在胃酸的作用下大量繁殖，一部分蟲(chóng)卵會(huì)在羊的體內(nèi)發(fā)育成包囊，另一部分蟲(chóng)卵會(huì)隨羊的糞便又排泄到牧區(qū)的水源、土壤和草場(chǎng)中。

牧民不但與羊群有密切接觸，而且與狗也有更加密切的接觸。牧民與羊群密切接觸或吃羊肉、喝羊血時(shí)，就會(huì)把蟲(chóng)卵吃入口中；狗有舔舐其肛門和身體的習(xí)慣，會(huì)把蟲(chóng)卵從肛門帶到其身體上，牧民與狗密切接觸時(shí)，就會(huì)把蟲(chóng)卵吃入口中；牧民喝被污染的水或食物時(shí)，更會(huì)將蟲(chóng)卵吃入口中。

生活在人胃中的棘球絳蟲(chóng)蟲(chóng)卵在胃酸的作用下加速繁殖、數(shù)量大增。一部分蟲(chóng)卵會(huì)隨血液鉆入人體的某些器官內(nèi)，如肝、肺、腦中，然后發(fā)育成巨型包囊；另一部分蟲(chóng)卵會(huì)隨人的糞便又排泄到牧區(qū)的水源、土壤和草場(chǎng)中。

棘球絳蟲(chóng)蟲(chóng)卵的傳播特征是糞口傳播，具有很強(qiáng)的傳染能力，既能夠在同物種類內(nèi)傳播，又能夠跨物種傳播。當(dāng)狗、羊和牧民傳染上包蟲(chóng)病后，一部分能夠治愈，另一部分則死亡。棘球絳蟲(chóng)蟲(chóng)卵攻擊的是狗、羊和牧民的部分特征（器官）。

將上述場(chǎng)景映射到對(duì)優(yōu)化問(wèn)題式（1）全局最優(yōu)解的搜索過(guò)程中，含義如下所述。

優(yōu)化問(wèn)題式（1）的搜索空間與草原牧區(qū)相對(duì)應(yīng)，該草原牧區(qū)中的一條狗、一只羊雞或一名牧民分別對(duì)應(yīng)于優(yōu)化問(wèn)題式（1）的一個(gè)試探解，即對(duì)類型為u的Nu個(gè)個(gè)體所對(duì)應(yīng)的試探解集就是且類型為u的個(gè)體i的特征j與試探解的變量相對(duì)應(yīng)。對(duì)于優(yōu)化問(wèn)題式（1），類型為u的個(gè)體i的適應(yīng)度Fit按下式計(jì)算：

2.2 可跨多物種多級(jí)傳播的包蟲(chóng)傳染病模型

在時(shí)期t，對(duì)于類型為u個(gè)體來(lái)說(shuō)，共存在5 種狀態(tài)：易感（未感染蟲(chóng)卵）的個(gè)體，其在時(shí)期t的占比為Su(t)，處于此狀態(tài)的個(gè)體用Su表示；已暴露（已感染蟲(chóng)卵但未發(fā)病）的個(gè)體，其在時(shí)期t的占比為Eu(t)，處于此狀態(tài)的個(gè)體用Eu表示；已發(fā)?。ǜ腥鞠x(chóng)卵后已發(fā)?。┑膫€(gè)體，其在時(shí)期t的占比為Iu(t)，處于此狀態(tài)的個(gè)體用Iu表示；已治愈的個(gè)體，其在時(shí)期t的占比為Ru(t)，處于此狀態(tài)的個(gè)體用Ru表示；已死亡的個(gè)體，其在時(shí)期t的占比為Du(t)，處于此狀態(tài)的用Du表示。

為了簡(jiǎn)單起見(jiàn)，對(duì)上述場(chǎng)景進(jìn)行一些簡(jiǎn)化：不考慮個(gè)體的自然死亡；當(dāng)一個(gè)個(gè)體因染上包蟲(chóng)病而死亡后，立即就有一個(gè)新個(gè)體添加到該草原牧區(qū)中，從而確保各類個(gè)體的總數(shù)為常數(shù)；棘球絳蟲(chóng)蟲(chóng)卵可以在同一物種內(nèi)傳播；不同物種間的傳播規(guī)律是：狗可以將蟲(chóng)卵傳播給羊和牧民，羊只能將蟲(chóng)卵傳播給牧民，牧民不會(huì)將蟲(chóng)卵傳播給狗和羊。上述場(chǎng)景可采用傳染病傳播倉(cāng)室模型[19-20]來(lái)描述，如圖1 所示。根據(jù)圖1，可以寫出其相應(yīng)的動(dòng)力學(xué)方程組如下：

式中，βu表示包蟲(chóng)病在類型為u的個(gè)體中的傳染率，0 ＜βu＜1;λu表示類型為u的個(gè)體的復(fù)原率，0 ＜λu＜1；γu表示類型為u的個(gè)體的發(fā)病率，0 ＜γu＜1；δu表示類型為u的個(gè)體的治愈率，0 ＜δu＜1；θu表示類型為u的個(gè)體的死亡率，0 ＜θu＜1。

時(shí)期t，對(duì)于類型為u的任意一個(gè)個(gè)體來(lái)說(shuō)，Su(t)、Eu(t)、Iu(t)、Ru(t)、Du(t)分別表示該個(gè)體處于狀態(tài)Su、Eu、Iu、Ru、Du的概率。必須指出，在同一時(shí)期，任意一個(gè)個(gè)體只能處在狀態(tài)Su、Eu、Iu、Ru、Du中的某一個(gè)。通常情況下模型式（3）中的參數(shù)βu、γu、δu、λu、θu的取值是隨時(shí)間變化的，可將式（3）應(yīng)用到類型為u的任意一個(gè)個(gè)體上，并將式（3）改寫成如式（4）所示離散遞推形式。

2.3 個(gè)體演化狀態(tài)識(shí)別

Fig.1 Flow chart of infectious disease transmission in interaction of dogs,sheep and herdsmen圖1 狗、羊與牧民相互作用的傳染病傳播流程圖

Fig.2 Legal types of state transition圖2 合法的狀態(tài)轉(zhuǎn)移類型

在時(shí)期t，采用模型式（4）計(jì)算類型為u的個(gè)體i的易感概率、暴露概率、發(fā)病概率、治愈概率和死亡概率；類型為u的個(gè)體i在時(shí)期t處于狀態(tài)Su、Eu、Iu、Ru、Du中的哪個(gè)狀態(tài)，由它們之中的最大者確定。依據(jù)圖1，可以識(shí)別出所有合法的狀態(tài)轉(zhuǎn)移類型，如圖2 所示，圖2 中所描述的狀態(tài)轉(zhuǎn)移類型的含義及其所對(duì)應(yīng)的算子如表1 所示。

Table 1 Legal state transitions表1 合法狀態(tài)轉(zhuǎn)換

從表1 可知，圖2 所確定的算子有3×9=27 個(gè)。較多的算子可提升HDO 算法的智能性。

2.4 演化算子設(shè)計(jì)

2.4.1 特征集合生成方法

設(shè)當(dāng)前個(gè)體類型u，個(gè)體編號(hào)為i，個(gè)體狀態(tài)s∈{Su,Eu,Iu,Ru,Du}，則：

2.4.2 狀態(tài)轉(zhuǎn)移算子設(shè)計(jì)方法

對(duì)圖2 進(jìn)行分解，可得如圖3 所示的三種狀態(tài)轉(zhuǎn)移，存在下列三種情況：

（1）個(gè)體從時(shí)期t的狀態(tài)A轉(zhuǎn)移到時(shí)期t+1 的狀態(tài)B，如圖3（a）所示，其中A、B∈{Su,Eu,Iu,Ru,Du|u=1,2,3}但A≠B。有如下兩種情形：

Fig.3 Three situations of individual state transition圖3 個(gè)體三種狀態(tài)轉(zhuǎn)移情形

情形1當(dāng)A≠Du，B≠Du時(shí)，大量的個(gè)體狀態(tài)轉(zhuǎn)移都屬于這種情況。例如Su→Eu、Eu→Iu等。為了實(shí)現(xiàn)圖3（a）所示的狀態(tài)轉(zhuǎn)移，可將已處于狀態(tài)B的若干個(gè)體的某些特征的狀態(tài)值經(jīng)加工處理后傳給處于狀態(tài)A的當(dāng)前個(gè)體的對(duì)應(yīng)特征，從而使其具有處于狀態(tài)B的個(gè)體的特征。例如，對(duì)于狀態(tài)轉(zhuǎn)移Su→Eu，將已處于暴露狀態(tài)(Eu)的若干個(gè)體的某些特征的狀態(tài)值經(jīng)加工處理后傳給處于易感狀態(tài)(Su)的當(dāng)前個(gè)體，即可使其感染上包蟲(chóng)病。此相當(dāng)于已感染棘球絳蟲(chóng)蟲(chóng)卵但未發(fā)病的暴露個(gè)體將其帶有蟲(chóng)卵的東西傳給當(dāng)前易感者，使其也感染上蟲(chóng)卵。其他轉(zhuǎn)移的含義可類似解釋。

情形2當(dāng)A=Iu，B=Du時(shí)，當(dāng)前個(gè)體從時(shí)期t的發(fā)病狀態(tài)Iu轉(zhuǎn)移到時(shí)期t+1 的死亡狀態(tài)Du，此時(shí)可將適應(yīng)度極低的虛弱個(gè)體用適應(yīng)度較高的強(qiáng)壯個(gè)體來(lái)取代，從而實(shí)現(xiàn)虛弱個(gè)體的死亡。此情形可隨機(jī)淘汰極差個(gè)體。

（2）當(dāng)前個(gè)體在時(shí)期t處于某個(gè)狀態(tài)A時(shí)，A∈{Su,Eu,Iu,Ru|u=1,2,3}，在時(shí)期t+1 沒(méi)有發(fā)生狀態(tài)轉(zhuǎn)移，即相當(dāng)于A→A，如圖3（b）所示。圖2 中的每個(gè)節(jié)點(diǎn)包含了圖3（b）所示的情形。例如，Su→Su、Eu→Eu、Iu→Iu和Ru→Ru。注意，不存在Du→Du。

生物個(gè)體在生存競(jìng)爭(zhēng)過(guò)程中總是盡量使自身強(qiáng)壯，以便更好地生存發(fā)展。為達(dá)到此目的，可將處于相同狀態(tài)的若干個(gè)強(qiáng)壯個(gè)體的某些特征的狀態(tài)值經(jīng)加工處理后傳給當(dāng)前個(gè)體的對(duì)應(yīng)特征，使得當(dāng)前個(gè)體獲得強(qiáng)壯個(gè)體的特征，從而向更好的方向發(fā)展。

（3）在時(shí)期t，當(dāng)前個(gè)體在處于狀態(tài)C的個(gè)體的作用下從狀態(tài)A轉(zhuǎn)移到狀態(tài)B，如圖3（c）所示，其中A∈{S2,S3},B∈{E2,E3},C∈{E1∪I1,E2∪I2}。例如，對(duì)于圖3（c）中的節(jié)點(diǎn)S2，有。

設(shè)當(dāng)前個(gè)體類型u∈{1,2,3}，個(gè)體編號(hào)為i，則HDO 算法中各個(gè)算子的數(shù)學(xué)表達(dá)式如表2 所示。

對(duì)于優(yōu)化問(wèn)題式（1），其生長(zhǎng)算子可以描述為：

2.5 HDO 算法構(gòu)造方法

（1）初始化：①令時(shí)期t=0，按第4 章介紹的方法設(shè)置本算法中的參數(shù)：演化時(shí)期數(shù)G，全局最優(yōu)解誤差ε，個(gè)體特征被包蟲(chóng)病攻擊的最大概率E0、L、N1、N2、N3；②分別初始化個(gè)體，u∈{1,2,3}。

（3）計(jì)算類型為u的個(gè)體i的Su、Eu、Iu、Ru、Du狀態(tài)，u∈{1,2,3}。函數(shù)GetSEIRD()用于確定類型為u的個(gè)體i將處于Su、Eu、Iu、Ru、Du狀態(tài)中的哪一個(gè)狀態(tài)。

Table 2 Mathematical expressions of operators in HDO algorithm表2 HDO 算法中各個(gè)算子的數(shù)學(xué)表達(dá)式

（4）執(zhí)行下列操作：

（5）結(jié)束。

函數(shù)GetSEIRD(pS,pE,pI,pR,pD)的定義如下：

FunctionGetSEIRD(pS,pE,pI,pR,pD)//pS、pE、pI、pR、pD分別表示易感狀態(tài)S、暴露狀態(tài)E、發(fā)病狀態(tài)I、治愈狀態(tài)R、死亡狀態(tài)D 的概率

2.6 HDO 算法的特性

（1）HDO 算法演化過(guò)程具有Markov 特性。從Su-Su、Su-Eu、Eu-Eu、Eu-Iu、Iu-Iu、Iu-Ru、Iu-Du、Ru-Ru、Ru-Su的定義知，當(dāng)i、k為種群編號(hào)，j=1～n,u∈{1,2,3}，時(shí)，這些算子均有如下特征：

從式（6）和式（7）可知，時(shí)期t+1（新一代）的試探解的產(chǎn)生只與該試探解在時(shí)期t（當(dāng)前代）的狀態(tài)有關(guān)，而與該試探解在時(shí)期t以前是如何演變到時(shí)期t所對(duì)應(yīng)的狀態(tài)的歷程無(wú)關(guān)。因此，依據(jù)Markov 特性的定義[21]，HDO 算法演化過(guò)程具有Markov 特性。

（2）HDO 算法演化過(guò)程具有“步步不差”特性。從生長(zhǎng)算子的定義可知，因總有，故新一代試探解的適應(yīng)度不會(huì)劣于其老一代試探解的適應(yīng)度。

（3）HDO 算法性能優(yōu)良。由于HDO 算法擁有27個(gè)算子，是迄今為止擁有算子最多的群智能算法，由文獻(xiàn)[22]已經(jīng)證明，群智能算法的性能優(yōu)良程度與算子個(gè)數(shù)成正比，因此HDO 算法具有性能優(yōu)良的特征。

（4）HDO 算法收斂速度快。HDO 算法利用包蟲(chóng)病病毒只攻擊個(gè)體的極少部分特征這一優(yōu)勢(shì)而獲得每次只需要處理n/1 000～n/100 個(gè)變量這一能力。因被處理變量大幅減少，故當(dāng)求解復(fù)雜優(yōu)化問(wèn)題，特別是高維優(yōu)化問(wèn)題時(shí)，能夠顯著提升收斂速度。

2.7 時(shí)間復(fù)雜度

令N=(N1+N2+N3)/3，HDO 算法的時(shí)間復(fù)雜度計(jì)算過(guò)程如表3 所示。

3 HDO 算法的收斂性證明

由HDO 算法知，草原牧區(qū)與搜索空間H等價(jià)，將N1只狗只羊和N3個(gè)牧民重新排列成M=N1+N2+N3個(gè)個(gè)體，形成新的個(gè)體序列為。每個(gè)個(gè)體是優(yōu)化問(wèn)題式（1）的一個(gè)試探解，其目標(biāo)函數(shù)值按式（1）計(jì)算為F(Xti)，所有個(gè)體的狀態(tài)所形成的集合為：

Table 3 Time complexity table of HDO algorithm表3 HDO 算法的時(shí)間復(fù)雜度計(jì)算表

進(jìn)一步令：

不失一般性，令F1即為所求的全局最優(yōu)解。由式（8）的下標(biāo)形成的集合為：

U={1,2,…,M}

?i∈U,i就是個(gè)體i執(zhí)行隨機(jī)搜索時(shí)可能處的狀態(tài)。假設(shè)在時(shí)期t個(gè)體i搜索到的最好目標(biāo)函數(shù)值為Fi，其對(duì)應(yīng)的狀態(tài)為i。顯然，由式（8）知，在時(shí)期t+1 個(gè)體i進(jìn)行搜索時(shí)，若向更優(yōu)的狀態(tài)k轉(zhuǎn)移，則應(yīng)滿足k＜i；相反，若向更差的狀態(tài)k轉(zhuǎn)移，則應(yīng)滿足k＞i，如圖4 所示。

?X∈H，有F1≤F(X)≤FM，將H劃分為如下非空子集：

另外，顯然有：

Fig.4 State transition in random search of HDO algorithm圖4 HDO 算法隨機(jī)搜索時(shí)的狀態(tài)轉(zhuǎn)移

引理1在HDO 算法中，i=1,2,…,M，滿足：

證明（1）引理式（10）的證明：設(shè)狀態(tài)i為時(shí)期t個(gè)體i的狀態(tài)，其在H 中對(duì)應(yīng)的位置為Xt，由2.6 節(jié)知，HDO 算法的隨機(jī)搜索過(guò)程具有步步不差的特性，故在時(shí)期t+1，個(gè)體i不會(huì)轉(zhuǎn)移到任何更差的狀態(tài)上去，如圖4 所示，故有：

（2）引理式（11）的證明：設(shè)狀態(tài)i為時(shí)期t個(gè)體i的狀態(tài)，在時(shí)期t+1，個(gè)體i隨機(jī)選擇算子Su-Su、Su-Eu、Eu-Eu、Eu-Iu、Iu-Iu、Iu-Ru、Iu-Du、Ru-Ru、Ru-Su進(jìn)行演化以便轉(zhuǎn)移到更好的狀態(tài)k上。此時(shí)，存在有如下兩種情況：

①若狀態(tài)i就是全局最優(yōu)狀態(tài)，也即i=1，則下一步轉(zhuǎn)移仍留在原狀態(tài)，即k=i=1，這是因?yàn)橛?.6節(jié)知，個(gè)體i不會(huì)轉(zhuǎn)移到比原狀態(tài)i更差的其他狀態(tài)上去，故必以概率p1,1=1 留在原狀態(tài)i上。因p1,1=1 ＞0,命題得證。

②若狀態(tài)i不是全局最優(yōu)狀態(tài)，則在狀態(tài)1 和當(dāng)前狀態(tài)i之間必至少存在一個(gè)中間狀態(tài)k，如圖4 所示，使得F1≤Fk＜Fi，也就是1 ≤k＜i，此時(shí)個(gè)體i可以從當(dāng)前狀態(tài)i轉(zhuǎn)移到更優(yōu)的新?tīng)顟B(tài)k上去，也即pi,k＞0。命題得證。

綜上所述，可得?k＜i,pi,k＞0。 □

定理1[23]設(shè)P′是一n階可歸約隨機(jī)矩陣，即通過(guò)相同的行和列變換后可得到，其中C是m階本原隨機(jī)矩陣，且T≠0,R≠0，則有：

上述矩陣是一個(gè)穩(wěn)定隨機(jī)矩陣，P′∞=1′P′∞,P′∞=P′0P′∞唯一確定且與初始分布無(wú)關(guān)，P′∞滿足條件：

定理2HDO 算法具有全局收斂性。

證明由2.6 節(jié)知，HDO 算法的搜索過(guò)程具有Markov 特性。對(duì)于每個(gè),i=1,2,…,M可看作是有限Markov 鏈上的一個(gè)狀態(tài)，根據(jù)引理中式（10）的結(jié)論可得該Markov 鏈的狀態(tài)轉(zhuǎn)移矩陣為：

由式（9）知，P′中每行的概率之和等于1。又根據(jù)引理中式（11）的結(jié)論可得：

由以上可知，P′是一個(gè)M階可歸約隨機(jī)矩陣，即Markov 狀態(tài)轉(zhuǎn)移矩陣，滿足定理1 的條件，故有：

因C∞=C=(1)，T∞=0，故必有R∞=(1,1,…,1)T，這是因?yàn)镻′中任意一行的概率之和總為1，故有：

上式表明，當(dāng)k→∞時(shí)，pi,1=1,i=1,2,…,M，即無(wú)論各個(gè)體初始狀態(tài)如何，最后都能以概率1 轉(zhuǎn)移到全局最優(yōu)狀態(tài)1 上去。于是有：

因此，HDO 算法具有全局收斂性。 □

4 HDO 算法的參數(shù)選擇

HDO 算法參數(shù)包括兩部分：一部分是包蟲(chóng)病傳染病模型參數(shù)，該部分參數(shù)為算法內(nèi)置參數(shù)，無(wú)需用戶再進(jìn)行設(shè)置；另一部分是算法運(yùn)行控制參數(shù)，此類參數(shù)需要用戶根據(jù)情況進(jìn)行設(shè)置。

（1）包蟲(chóng)病傳染病模型參數(shù)確定方法。包蟲(chóng)病傳染病模型參數(shù)的選擇依據(jù)是確保具有足夠的隨機(jī)性，u∈{1,2,3}。依據(jù)文獻(xiàn)[19]介紹的參數(shù)取值方法并經(jīng)隨機(jī)化后，可得βu=Rand(0.42,0.93),γu=Rand(0.12,0.24),δu=Rand(0.31,0.57),λu=Rand(0.32,0.53),θu=Rand(0.18,0.43)。應(yīng)用此取值策略，任取測(cè)試情況如圖5 所示。從圖5 可知，具有極好的隨機(jī)性。參數(shù)βu、γu、δu、λu、θu的取值方法可作為算法內(nèi)置參數(shù)進(jìn)行設(shè)置，無(wú)須用戶干預(yù)。

Fig.5 Randomicity of state E2圖5 狀態(tài)E2出現(xiàn)的隨機(jī)性

（2）算法運(yùn)行控制參數(shù)設(shè)置方法。HDO 算法的運(yùn)行控制參數(shù)有G、ε、E0、L、N1、N2、N3。G和ε是兩個(gè)互補(bǔ)參數(shù)，只要滿足一個(gè)即可，ε由所求解的工程問(wèn)題決定，通?？扇ˇ?10-5～10-8即可；G可取充分大，不妨設(shè)G=108。HDO 算法關(guān)鍵參數(shù)只有E0、L、N1、N2和N3，可令N=N1=N2=N3。下面主要討論關(guān)鍵參數(shù)E0、L、N的取值方法。由于BUMP 優(yōu)化問(wèn)題與本文所要求解的工程問(wèn)題形態(tài)相似，且BUMP優(yōu)化問(wèn)題極難求解，故下面以BUMP 優(yōu)化問(wèn)題為例來(lái)探明E0、L、N的取值方法。BUMP 優(yōu)化問(wèn)題如下所示。

設(shè)平均最優(yōu)目標(biāo)函數(shù)值用Y表示，平均CPU 時(shí)間用T表示。當(dāng)L取不同值時(shí)，采用HDO 算法求解BUMP 優(yōu)化問(wèn)題，令n=50,G=108,N=100,E0=0.01,運(yùn)行100 次，表4 描述了L與Y和T之間的關(guān)系。結(jié)果表明，當(dāng)L=3～7 時(shí)，Y的精度達(dá)到最高，而T增加較低。由此可見(jiàn)，L=3～7 為L(zhǎng)的最佳取值區(qū)間。

Table 4 Relationship of L with Y and T表4 L 與Y 和T 之間的關(guān)系

令n=50,G=108,N=100,E0=0.01,L=3，HDO 算法運(yùn)行100 次。表5 描述了參數(shù)E0與Y和T之間的關(guān)系。結(jié)果表明，當(dāng)E0=0.006～0.100 時(shí)，Y的精度較高，且T較少；當(dāng)E0＞0.100 時(shí)，T增加很大，且Y的精度也大大降低；特別是當(dāng)E0=1.000 時(shí)，無(wú)法獲得最佳解。由此可見(jiàn)，E0的最佳取值區(qū)間為E0=0.006～0.100。

令G=108,N=100,L=3，HDO 算法運(yùn)行100 次。表6 描述了Y、n、N和T之間的關(guān)系。從表6 可以看到：

（1）當(dāng)n增加時(shí)，T大大增加；

（2）對(duì)于給定的n，如果N增加，T大大增加；

Table 5 Rrelationship of E0 with Y and T表5 E0與Y 和T 之間的關(guān)系

（3）對(duì)于給定的n和N，如果E0增加，Y的精度會(huì)降低，但T增加。

因此，如果n＞500,N=50～100 就足夠了；如果n＜500，N=100 就足夠了。

5 應(yīng)用實(shí)例研究

5.1 VOCs最優(yōu)減排模型

某個(gè)關(guān)聯(lián)區(qū)域由n個(gè)區(qū)域組成，在氣象因素作用下，每個(gè)區(qū)域所排放的VOCs（揮發(fā)性有機(jī)廢氣）一部分留存在本區(qū)域，另一部分會(huì)擴(kuò)散到其他區(qū)域。關(guān)聯(lián)區(qū)域VOCs 聯(lián)防聯(lián)控減排方案的含義是，如何控制關(guān)聯(lián)區(qū)域內(nèi)每個(gè)區(qū)域的VOCs 排放量，才能使關(guān)聯(lián)區(qū)域內(nèi)VOCs 對(duì)大氣環(huán)境的影響降到最低。假設(shè)VOCs減排工作已進(jìn)行了t-1 期，在當(dāng)前時(shí)期t，VOCs 在n個(gè)區(qū)域的預(yù)測(cè)總產(chǎn)生量分別為Q1(t),Q2(t),…,Qn(t)；對(duì)區(qū)域i來(lái)說(shuō)，其他區(qū)域遷入到區(qū)域i的VOCs 分別為R1i(t),R2i(t),…,Rni(t),Rji(t)≥0,j=1～n,i=1～n；而從區(qū)域i遷出的VOCs 分別為Si1(t),Si2(t),…,Sin(t),Sij(t)≥0,j=1～n,i=1～n。在時(shí)期t，n個(gè)區(qū)域的減排量分別為x1(t),x2(t),…,xn(t)，它們是待求的變量，{x1(t),x2(t),…,xn(t)}的一種取值方案構(gòu)成一個(gè)時(shí)期t的減排方案。由于該方案既考慮到了前t-1 期的減排情況，又考慮了n個(gè)區(qū)域間的相互影響，因此該減排方案具有跨時(shí)間和跨空間的特征，是一種跨時(shí)空協(xié)同減排方案。關(guān)聯(lián)區(qū)域VOCs聯(lián)防聯(lián)控最優(yōu)減排模型如式（12）所示。

Table 6 Relationship of E0,n,N with Y and T表6 E0、n、N 與Y 和T 之間的關(guān)系

式中，F(xiàn)(X(t))為減排方案的總效用；減排總量控制目標(biāo)；減排任務(wù)控制目標(biāo)；政府補(bǔ)貼效用控制目標(biāo)；總成本控制目標(biāo)f4i(xi(t))=-ci(t)xi(t)；X(t)=(x1(t),x2(t),…,xn(t))；ci(t)為區(qū)域i的單位VOCs 減排成本；Pi(s)為時(shí)期s區(qū)域i的VOCs 實(shí)際減排量；Vi為受區(qū)域i影響的其他區(qū)域集合；V0為受區(qū)域間相互影響的最大允許值，簡(jiǎn)稱區(qū)域影響系數(shù)；wi為區(qū)域i的補(bǔ)貼系數(shù)，區(qū)域越重要，補(bǔ)貼系數(shù)越大；O1～O4為4 個(gè)目標(biāo)函數(shù)的優(yōu)先級(jí)，優(yōu)先級(jí)次序要求滿足O1＞O2＞O3＞O4，因此可設(shè)O1=1 000，O2=100,O3=10,O4=1;qt為前t-1 期未完成的減排量在時(shí)期t所分?jǐn)偟谋壤?jiǎn)稱分?jǐn)偙壤?/p>

計(jì)算時(shí)，減排方案X(t)也稱為試探解；若減排方案X(t) 不滿足約束條件，則令f(X(t))=-∞,f(X(t))∈{f1(X(t)),f2(X(t)),f3(X(t)),f4(X(t))}。優(yōu)化模型式（12）是一個(gè)非線性優(yōu)化問(wèn)題，傳統(tǒng)的基于函數(shù)連續(xù)性和可導(dǎo)性的數(shù)學(xué)優(yōu)化方法無(wú)法求解該優(yōu)化模型。本文采用HDO 算法對(duì)其求解。

5.2 優(yōu)化模型求解過(guò)程

本文以西安市為例來(lái)說(shuō)明本文所提出方法的使用過(guò)程。結(jié)合西安市各區(qū)縣2018 年10 月至12 月VOCs 排放情況來(lái)制定2019 年1 月份的VOCs 最佳 減排方案，以此來(lái)說(shuō)明算法HDO 算法的使用方法。表7 給 出了該市2018 年10 月至12 月VOCs 實(shí)際減 排量，西安市區(qū)縣數(shù)n=13。

表8為該市2018年10月至12月VOCs產(chǎn)生量。利用EViews8軟件可以預(yù)測(cè)出該市各區(qū)縣在2019年1月份的VOCs產(chǎn)生量，如表8 的最后一列所示。

依據(jù)該市的氣象規(guī)律，每年10 月至12 月的氣象特性較為相似。因此，依據(jù)該市各區(qū)縣2018年10月至12 月的氣象資料，采用HYSPLIT4 模式軟件可以計(jì)算出2018 年10 月至12 月各區(qū)縣VOCs 遷入和遷出量，如表9 所示。表9 最后一列是對(duì)2019 年1 月各區(qū)縣VOCs 遷入和遷出百分比的預(yù)測(cè)值。從表9 也可以確定出不同時(shí)期受某個(gè)區(qū)縣影響的其他區(qū)縣的集合。

Table 7 Actual emission reduction of VOCs in various counties from October to December 2018表7 2018 年10 月至12 月各區(qū)縣VOCs實(shí)際減排情況 億m3

Table 8 VOCs production amount in various counties from October to December 2018表8 2018年10月至12月各區(qū)縣VOCs產(chǎn)生量 億m3

表10 給出了該市各區(qū)縣VOCs 減排成本和減排補(bǔ)貼系數(shù)。前期未完成的減排量在當(dāng)前時(shí)期所分?jǐn)偟谋壤畹蜑閝t=7%；區(qū)域間相互影響的最大允許值V0=5%。

選擇7 種優(yōu)化算法與HDO 算法進(jìn)行比較，這些算法包括NP-PSO（non-parametric particle swarm optimization）[7]、DASA（differential ant-stigmergy algor-ithm）[8]、ABC（artificial bee colony）[10]、RCGA（realcoded genetic algorithm）[13]、MBBO（metropolis biogeography-based optimization）[24]、DE（differential evolution）[25]和SaDE（self-adaptive differential evolution）[26]。RCGA 是一種實(shí)數(shù)編碼型遺傳算法，搜索效率高；DASA 是一種仿蟻群算法設(shè)計(jì)思路的新型算法，收斂速度很快；NP-PSO 是一種新型粒子群算法，不需要參數(shù)設(shè)置，收斂速度很快；MBBO 是BBO 算法的改進(jìn)版，強(qiáng)化了島嶼的都市特征，搜索效率極佳；DE 通過(guò)在基本DE 算法中引入了局部誘導(dǎo)遺傳算子，大幅提升了其收斂速度；SaDE 是在傳統(tǒng)自適應(yīng)差分算法中引入了高效自調(diào)整參數(shù)演化算子，使得該算法的搜索效率大幅提升；ABC 是一種基于基因組合參數(shù)控制策略的蜂群算法，收斂速度很快。

Table 9 Emigration percentage of VOCs in counties from October 2018 to January 2019表9 2018 年10 月至2019 年1 月各區(qū)縣VOCs遷出百分比

Table 10 Cost of VOCs emission reduction and subsidies coefficient for emission reduction in counties表10 各區(qū)縣VOCs減排成本和減排補(bǔ)貼系數(shù)

計(jì)算時(shí)，7 種優(yōu)化算法的參數(shù)按表11 進(jìn)行初始化；HDO 算法的參數(shù)設(shè)置為G=108,N=100,L=3,E0=0.01。

采用HDO 算法和7 個(gè)比較算法進(jìn)行求解，每個(gè)算法共運(yùn)行50 次，得到2019 年1 月西安市各區(qū)縣VOCs最佳減排方案平均值如表12 所示。

從表12 可以看出，HDO 算法所獲得的目標(biāo)函數(shù)值要優(yōu)于其他7 種算法。對(duì)新城區(qū)來(lái)說(shuō)，HDO 算法所獲得的VOCs 最佳減排方案最好，SaDE 和DE 算法所獲得的VOCs 最佳減排方案略遜于HDO 算法，但好于RCGA、DASA 和ABC 算法，而NP-PSO、MBBO所獲得的VOCs 最佳減排方案最差；對(duì)其他區(qū)縣來(lái)說(shuō)，各算法所獲得的VOCs 最佳減排方案對(duì)比可作類似的分析。與其他算法相比，HDO 算法發(fā)現(xiàn)VOCs最佳減排方案的平均計(jì)算時(shí)間只有43 s，略優(yōu)于RCGA(66 s)、ABC(74 s)、SaDE(79 s)和DE(87 s)，明顯優(yōu)于MBBO(102 s)、NP-PSO(243 s)和DASA(253 s)；HDO 算法發(fā)現(xiàn)VOCs 最佳減排方案的平均適應(yīng)度計(jì)算次數(shù)只有58 473 次，明顯優(yōu)于RCGA(202 193 次)和ABC(243 192 次)，大幅優(yōu)于SaDE(338 794 次)和DE(357 597 次)，極優(yōu)于MBBO(464 381 次)、NP-PSO(523 675 次)和DASA(572 247 次)。圖6 給出了各算法求解優(yōu)化模型式（12）時(shí)的樣本收斂曲線，從該圖也可以看出，HDO 算法的收斂過(guò)程要明顯快于其他7種參與比較的算法。

Table 11 Parameters of 7 optimization algorithms表11 7 種優(yōu)化算法的參數(shù)

Table 12 Average of VOCs optimum emission reduction schemes for counties in January 2019表12 2019 年1 月各區(qū)縣VOCs最佳減排方案平均值 億m3

Fig.6 Sample convergence curve圖6 樣本收斂曲線

6 結(jié)束語(yǔ)

與其他算法相比，HDO 算法具有如下優(yōu)點(diǎn)：

（1）HDO 算法中包括形態(tài)為Su-Su、Su-Eu、Eu-Eu、Eu-Iu、Iu-Iu、Iu-Ru、Iu-Du、Ru-Ru、Ru-Su的27 個(gè)算子，擁有3 個(gè)不同物種的個(gè)體，可顯著地提升算法的搜索能力。

（2）在HDO 算法中，Su-Su、Eu-Eu、Iu-Iu、Ru-Ru算子可利用強(qiáng)壯個(gè)體的特征來(lái)改善虛弱個(gè)體的特征，從而提升算法的求精（exploitation）能力；Su-Eu、Eu-Iu、Iu-Ru、Ru-Su算子可改良個(gè)體的適應(yīng)度分布特征，從而提升算法的探索（exploration）能力；Iu-Du算子可使極虛弱個(gè)體得到有效清除，從而降低算法陷入局部陷阱的概率。

（3）HDO 算法進(jìn)行演化計(jì)算時(shí)，每次只需要處理極少部分?jǐn)?shù)量，具有自動(dòng)降維特性，具有求解高維優(yōu)化問(wèn)題能力。

（4）HDO 算法具有全局收斂性。

HDO 算法今后的改進(jìn)方向：

（1）已經(jīng)發(fā)現(xiàn)，某些傳染病能夠跨不低于4 個(gè)物種，可以利用HDO 算法的設(shè)計(jì)思路，提出跨多物種的傳染病優(yōu)化算法。

（2）將HDO 算法的狀態(tài)數(shù)從當(dāng)前的S（易感）、E（暴露）、I（染?。?、R（治愈）、D（死亡）等5個(gè)狀態(tài)擴(kuò)展到S（易感）、E（暴露）、I（發(fā)?。（免疫）、R（治愈）、D（死亡）等6個(gè)狀態(tài)，從而使HDO算法擁有更多的算子。

（3）深入分析Su-Su、Su-Eu、Eu-Eu、Eu-Iu、Iu-Iu、Iu-Ru、Iu-Du、Ru-Ru、Ru-Su的性能。

（4）在HDO 算法中納入DNA 機(jī)制、免疫機(jī)制，可使HDO 算法的研究更加深入。