具有跨物種多級傳播特征的包蟲病優化算法*

黃光球，陸秋琴

西安建筑科技大學 管理學院，西安710055

1 引言

工程中存在一類非線性優化問題，此類問題的數學表達式中存在不可導的數學表達式。基于梯度的優化方法無法求解此類問題。為了求解此類問題，人們提出了啟發式搜索方法[1]，群智能優化算法[2]就是其中的一種。因群智能優化算法對工程優化問題沒有特殊要求，故具有廣泛的適應性。近幾年，群智能優化算法得到快速發展。

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

隨著互聯網的高速發展，社會已進入大數據時代，人們所遇到的問題變得越來越復雜，故需要提出高智能性的算法來解決這些復雜問題。對于群智能算法來說，如何增加群智能算法的智能特征，顯然是需要解決的關鍵問題。解決該問題的一種策略是：精心選擇出某種具有特殊性質的生物進化場景，巧妙解決目前群智能算法所存在問題，即可設計出高智能性的群智能算法。

若一類生物進化場景具有豐富的演化內涵，其中個體之間的相互作用關系十分豐富，則依據此類生物進化場景就可以設計出具有算子眾多、邏輯結構清晰的高智能性群智能優化算法。能夠跨三個物種實施多級傳播的包蟲傳染病攻擊家犬、牛羊和人類的過程，就是這樣一個場景。

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

由于包蟲病可同時危害人類和家畜養殖業，近年來全世界研究人員對該傳染病給予了高度關注，主要研究成果表現為如下幾方面：

（1）包蟲病流行病學分析。為了了解包蟲病的流行病學特征，通過對包蟲病病例進行流行病學特征調查研究，探討病例可能的感染來源，為制定包蟲病防控措施提供依據[15]。

（2）包蟲病免疫學研究。通過對家畜和人的包蟲病免疫機理進行研究，為包蟲病的防治和新型疫苗的生產提供科學參考[16]。

（3）包蟲病診斷方法研究。提出包蟲病的診斷方法，為包蟲病的治療提供科學、快速和準確的依據[17]。

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

包蟲病具有跨多物種傳播特征，能夠對人類造成很大痛苦，但在這悲慘的自然現象后面，卻隱藏著一個性能優良的群智能優化算法，本文稱其為包蟲傳染病優化算法（hydatid disease optimization，HDO）。本文著重解決如下五個問題：

（1）如何將包蟲傳染病模型轉化為能求解復雜優化問題的HDO 算法。

（2）如何使得HDO 算法中的算子能充分反映動物個體之間的相互作用關系，以便體現包蟲傳染病模型的思想。

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

（4）如何確定HDO 算法的最佳參數設置。

（5）如何用HDO 算法求解關聯區域VOCs（volatile organic compounds）聯防聯控最優減排優化問題。

2 包蟲傳染病模型的優化算法設計

設要求解的優化問題為：

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

2.1 算法場景設計

在一個草原牧區生活有一群牧民，其數量有N3名，其編號分別是1,2,…,N3；牧民們以飼養藏羊為生，藏羊的數量有N2只，其編號分別是1,2,…,N2；與此同時，牧民們還飼養了一群狗用于控制和保護羊群，狗的數量有N1只，其編號分別是1,2,…,N1。一條狗、一只羊或一名牧民統稱為一個個體，每個個體均由n個特征（器官）來表征。令u表示個體類型，u=1表示狗類，u=2 表示羊類，u=3 表示牧民類，即對類型為u的個體i來說，其表征特征為。

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

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

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

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

將上述場景映射到對優化問題式（1）全局最優解的搜索過程中，含義如下所述。

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

2.2 可跨多物種多級傳播的包蟲傳染病模型

在時期t，對于類型為u個體來說，共存在5 種狀態：易感（未感染蟲卵）的個體，其在時期t的占比為Su(t)，處于此狀態的個體用Su表示；已暴露（已感染蟲卵但未發病）的個體，其在時期t的占比為Eu(t)，處于此狀態的個體用Eu表示；已發病（感染蟲卵后已發病）的個體，其在時期t的占比為Iu(t)，處于此狀態的個體用Iu表示；已治愈的個體，其在時期t的占比為Ru(t)，處于此狀態的個體用Ru表示；已死亡的個體，其在時期t的占比為Du(t)，處于此狀態的用Du表示。

為了簡單起見，對上述場景進行一些簡化：不考慮個體的自然死亡；當一個個體因染上包蟲病而死亡后，立即就有一個新個體添加到該草原牧區中，從而確保各類個體的總數為常數；棘球絳蟲蟲卵可以在同一物種內傳播；不同物種間的傳播規律是：狗可以將蟲卵傳播給羊和牧民，羊只能將蟲卵傳播給牧民，牧民不會將蟲卵傳播給狗和羊。上述場景可采用傳染病傳播倉室模型[19-20]來描述，如圖1 所示。根據圖1，可以寫出其相應的動力學方程組如下：

式中，βu表示包蟲病在類型為u的個體中的傳染率，0 ＜βu＜1;λu表示類型為u的個體的復原率，0 ＜λu＜1；γu表示類型為u的個體的發病率，0 ＜γu＜1；δu表示類型為u的個體的治愈率，0 ＜δu＜1；θu表示類型為u的個體的死亡率，0 ＜θu＜1。

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

2.3 個體演化狀態識別

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，采用模型式（4）計算類型為u的個體i的易感概率、暴露概率、發病概率、治愈概率和死亡概率；類型為u的個體i在時期t處于狀態Su、Eu、Iu、Ru、Du中的哪個狀態，由它們之中的最大者確定。依據圖1，可以識別出所有合法的狀態轉移類型，如圖2 所示，圖2 中所描述的狀態轉移類型的含義及其所對應的算子如表1 所示。

Table 1 Legal state transitions表1 合法狀態轉換

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

2.4 演化算子設計

2.4.1 特征集合生成方法

設當前個體類型u，個體編號為i，個體狀態s∈{Su,Eu,Iu,Ru,Du}，則：

2.4.2 狀態轉移算子設計方法

對圖2 進行分解，可得如圖3 所示的三種狀態轉移，存在下列三種情況：

（1）個體從時期t的狀態A轉移到時期t+1 的狀態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 個體三種狀態轉移情形

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

情形2當A=Iu，B=Du時，當前個體從時期t的發病狀態Iu轉移到時期t+1 的死亡狀態Du，此時可將適應度極低的虛弱個體用適應度較高的強壯個體來取代，從而實現虛弱個體的死亡。此情形可隨機淘汰極差個體。

（2）當前個體在時期t處于某個狀態A時，A∈{Su,Eu,Iu,Ru|u=1,2,3}，在時期t+1 沒有發生狀態轉移，即相當于A→A，如圖3（b）所示。圖2 中的每個節點包含了圖3（b）所示的情形。例如，Su→Su、Eu→Eu、Iu→Iu和Ru→Ru。注意，不存在Du→Du。

生物個體在生存競爭過程中總是盡量使自身強壯，以便更好地生存發展。為達到此目的，可將處于相同狀態的若干個強壯個體的某些特征的狀態值經加工處理后傳給當前個體的對應特征，使得當前個體獲得強壯個體的特征，從而向更好的方向發展。

（3）在時期t，當前個體在處于狀態C的個體的作用下從狀態A轉移到狀態B，如圖3（c）所示，其中A∈{S2,S3},B∈{E2,E3},C∈{E1∪I1,E2∪I2}。例如，對于圖3（c）中的節點S2，有。

設當前個體類型u∈{1,2,3}，個體編號為i，則HDO 算法中各個算子的數學表達式如表2 所示。

對于優化問題式（1），其生長算子可以描述為：

2.5 HDO 算法構造方法

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

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

Table 2 Mathematical expressions of operators in HDO algorithm表2 HDO 算法中各個算子的數學表達式

（4）執行下列操作：

（5）結束。

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

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

2.6 HDO 算法的特性

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

從式（6）和式（7）可知，時期t+1（新一代）的試探解的產生只與該試探解在時期t（當前代）的狀態有關，而與該試探解在時期t以前是如何演變到時期t所對應的狀態的歷程無關。因此，依據Markov 特性的定義[21]，HDO 算法演化過程具有Markov 特性。

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

（3）HDO 算法性能優良。由于HDO 算法擁有27個算子，是迄今為止擁有算子最多的群智能算法，由文獻[22]已經證明，群智能算法的性能優良程度與算子個數成正比，因此HDO 算法具有性能優良的特征。

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

2.7 時間復雜度

令N=(N1+N2+N3)/3，HDO 算法的時間復雜度計算過程如表3 所示。

3 HDO 算法的收斂性證明

由HDO 算法知，草原牧區與搜索空間H等價，將N1只狗只羊和N3個牧民重新排列成M=N1+N2+N3個個體，形成新的個體序列為。每個個體是優化問題式（1）的一個試探解，其目標函數值按式（1）計算為F(Xti)，所有個體的狀態所形成的集合為：

Table 3 Time complexity table of HDO algorithm表3 HDO 算法的時間復雜度計算表

進一步令：

不失一般性，令F1即為所求的全局最優解。由式（8）的下標形成的集合為：

U={1,2,…,M}

?i∈U,i就是個體i執行隨機搜索時可能處的狀態。假設在時期t個體i搜索到的最好目標函數值為Fi，其對應的狀態為i。顯然，由式（8）知，在時期t+1 個體i進行搜索時，若向更優的狀態k轉移，則應滿足k＜i；相反，若向更差的狀態k轉移，則應滿足k＞i，如圖4 所示。

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

另外，顯然有：

Fig.4 State transition in random search of HDO algorithm圖4 HDO 算法隨機搜索時的狀態轉移

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

證明（1）引理式（10）的證明：設狀態i為時期t個體i的狀態，其在H 中對應的位置為Xt，由2.6 節知，HDO 算法的隨機搜索過程具有步步不差的特性，故在時期t+1，個體i不會轉移到任何更差的狀態上去，如圖4 所示，故有：

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

①若狀態i就是全局最優狀態，也即i=1，則下一步轉移仍留在原狀態，即k=i=1，這是因為由2.6節知，個體i不會轉移到比原狀態i更差的其他狀態上去，故必以概率p1,1=1 留在原狀態i上。因p1,1=1 ＞0,命題得證。

②若狀態i不是全局最優狀態，則在狀態1 和當前狀態i之間必至少存在一個中間狀態k，如圖4 所示，使得F1≤Fk＜Fi，也就是1 ≤k＜i，此時個體i可以從當前狀態i轉移到更優的新狀態k上去，也即pi,k＞0。命題得證。

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

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

上述矩陣是一個穩定隨機矩陣，P′∞=1′P′∞,P′∞=P′0P′∞唯一確定且與初始分布無關，P′∞滿足條件：

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

證明由2.6 節知，HDO 算法的搜索過程具有Markov 特性。對于每個,i=1,2,…,M可看作是有限Markov 鏈上的一個狀態，根據引理中式（10）的結論可得該Markov 鏈的狀態轉移矩陣為：

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

由以上可知，P′是一個M階可歸約隨機矩陣，即Markov 狀態轉移矩陣，滿足定理1 的條件，故有：

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

上式表明，當k→∞時，pi,1=1,i=1,2,…,M，即無論各個體初始狀態如何，最后都能以概率1 轉移到全局最優狀態1 上去。于是有：

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

4 HDO 算法的參數選擇

HDO 算法參數包括兩部分：一部分是包蟲病傳染病模型參數，該部分參數為算法內置參數，無需用戶再進行設置；另一部分是算法運行控制參數，此類參數需要用戶根據情況進行設置。

（1）包蟲病傳染病模型參數確定方法。包蟲病傳染病模型參數的選擇依據是確保具有足夠的隨機性，u∈{1,2,3}。依據文獻[19]介紹的參數取值方法并經隨機化后，可得β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)。應用此取值策略，任取測試情況如圖5 所示。從圖5 可知，具有極好的隨機性。參數βu、γu、δu、λu、θu的取值方法可作為算法內置參數進行設置，無須用戶干預。

Fig.5 Randomicity of state E2圖5 狀態E2出現的隨機性

（2）算法運行控制參數設置方法。HDO 算法的運行控制參數有G、ε、E0、L、N1、N2、N3。G和ε是兩個互補參數，只要滿足一個即可，ε由所求解的工程問題決定，通常可取ε=10-5～10-8即可；G可取充分大，不妨設G=108。HDO 算法關鍵參數只有E0、L、N1、N2和N3，可令N=N1=N2=N3。下面主要討論關鍵參數E0、L、N的取值方法。由于BUMP 優化問題與本文所要求解的工程問題形態相似，且BUMP優化問題極難求解，故下面以BUMP 優化問題為例來探明E0、L、N的取值方法。BUMP 優化問題如下所示。

設平均最優目標函數值用Y表示，平均CPU 時間用T表示。當L取不同值時，采用HDO 算法求解BUMP 優化問題，令n=50,G=108,N=100,E0=0.01,運行100 次，表4 描述了L與Y和T之間的關系。結果表明，當L=3～7 時，Y的精度達到最高，而T增加較低。由此可見，L=3～7 為L的最佳取值區間。

Table 4 Relationship of L with Y and T表4 L 與Y 和T 之間的關系

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

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

（1）當n增加時，T大大增加；

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

Table 5 Rrelationship of E0 with Y and T表5 E0與Y 和T 之間的關系

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

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

5 應用實例研究

5.1 VOCs最優減排模型

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

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

式中，F(X(t))為減排方案的總效用；減排總量控制目標；減排任務控制目標；政府補貼效用控制目標；總成本控制目標f4i(xi(t))=-ci(t)xi(t)；X(t)=(x1(t),x2(t),…,xn(t))；ci(t)為區域i的單位VOCs 減排成本；Pi(s)為時期s區域i的VOCs 實際減排量；Vi為受區域i影響的其他區域集合；V0為受區域間相互影響的最大允許值，簡稱區域影響系數；wi為區域i的補貼系數，區域越重要，補貼系數越大；O1～O4為4 個目標函數的優先級，優先級次序要求滿足O1＞O2＞O3＞O4，因此可設O1=1 000，O2=100,O3=10,O4=1;qt為前t-1 期未完成的減排量在時期t所分攤的比例，簡稱分攤比例。

計算時，減排方案X(t)也稱為試探解；若減排方案X(t) 不滿足約束條件，則令f(X(t))=-∞,f(X(t))∈{f1(X(t)),f2(X(t)),f3(X(t)),f4(X(t))}。優化模型式（12）是一個非線性優化問題，傳統的基于函數連續性和可導性的數學優化方法無法求解該優化模型。本文采用HDO 算法對其求解。

5.2 優化模型求解過程

本文以西安市為例來說明本文所提出方法的使用過程。結合西安市各區縣2018 年10 月至12 月VOCs 排放情況來制定2019 年1 月份的VOCs 最佳 減排方案，以此來說明算法HDO 算法的使用方法。表7 給 出了該市2018 年10 月至12 月VOCs 實際減 排量，西安市區縣數n=13。

表8為該市2018年10月至12月VOCs產生量。利用EViews8軟件可以預測出該市各區縣在2019年1月份的VOCs產生量，如表8 的最后一列所示。

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

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

Table 8 VOCs production amount in various counties from October to December 2018表8 2018年10月至12月各區縣VOCs產生量 億m3

表10 給出了該市各區縣VOCs 減排成本和減排補貼系數。前期未完成的減排量在當前時期所分攤的比例最低為qt=7%；區域間相互影響的最大允許值V0=5%。

選擇7 種優化算法與HDO 算法進行比較，這些算法包括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 是一種實數編碼型遺傳算法，搜索效率高；DASA 是一種仿蟻群算法設計思路的新型算法，收斂速度很快；NP-PSO 是一種新型粒子群算法，不需要參數設置，收斂速度很快；MBBO 是BBO 算法的改進版，強化了島嶼的都市特征，搜索效率極佳；DE 通過在基本DE 算法中引入了局部誘導遺傳算子，大幅提升了其收斂速度；SaDE 是在傳統自適應差分算法中引入了高效自調整參數演化算子，使得該算法的搜索效率大幅提升；ABC 是一種基于基因組合參數控制策略的蜂群算法，收斂速度很快。

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

Table 10 Cost of VOCs emission reduction and subsidies coefficient for emission reduction in counties表10 各區縣VOCs減排成本和減排補貼系數

計算時，7 種優化算法的參數按表11 進行初始化；HDO 算法的參數設置為G=108,N=100,L=3,E0=0.01。

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

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

Table 11 Parameters of 7 optimization algorithms表11 7 種優化算法的參數

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

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

6 結束語

與其他算法相比，HDO 算法具有如下優點：

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

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

（3）HDO 算法進行演化計算時，每次只需要處理極少部分數量，具有自動降維特性，具有求解高維優化問題能力。

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

HDO 算法今后的改進方向：

（1）已經發現，某些傳染病能夠跨不低于4 個物種，可以利用HDO 算法的設計思路，提出跨多物種的傳染病優化算法。

（2）將HDO 算法的狀態數從當前的S（易感）、E（暴露）、I（染病）、R（治愈）、D（死亡）等5個狀態擴展到S（易感）、E（暴露）、I（發病）、V（免疫）、R（治愈）、D（死亡）等6個狀態，從而使HDO算法擁有更多的算子。

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

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