具有跨物种多级传播特征的包虫病优化算法*

黄光球，陆秋琴

西安建筑科技大学 管理学院，西安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 算法的研究更加深入。