生物界经过数以亿万年的进化发展，形成了无数精妙的形态和行为机制，当科学家耗费几十年、甚至穷尽一生才得到解决难题的设计方案时，往往发现这样精巧的设计在大自然早已存在，模仿大自然成为产生新方法的有效途径之一，并且取得了丰硕的成果[1-2]．仿生算法是仿生学在算法领域应用的结果，即通过模仿生物的行为机制，构造能解决现实问题的算法[3-4]．与传统优化方法相比，仿生算法在求解非确定性多项式难题(NP-hard问题)方面表现出了显著的优越性，大量仿生算法相继被提出，例如模拟生物进化机制的遗传算法[5-6]，模拟蚂蚁觅食机制的蚁群算法[7-8]，模拟鸟群觅食机制的粒子群算法[9-10]，模拟人类大脑神经系统的人工神经网络和深度学习算法[11-12]，等等．当然，虽然仿生算法已经取得了长足的发展，但迄今为止，在复杂优化问题的求解方面，依然还存在很多难题，这也使得不断地探索新的仿生算法，以发现更多的求解问题的新途径非常必要．1 牵制平衡算法原理1.1　生态平衡的优化机制生物界虽然纷繁复杂，但数以千万计的生物之间能奇妙地达到生态平衡，从数学复杂度来讲，由于涉及的设计变量太多、约束条件太复杂，这几乎是一个无法求解的问题．而生物界以一个简单的机制，实现了这种精妙的平衡，即生物种群数量之间的牵制机制．生物之间的数量通过在一定容差范围内的相互牵制，实现了生物界的动态平衡，即生态平衡，从而实现了整个生物界种群规模的相对稳定．从生态平衡机制可以看出相互之间的牵制约束，也可以形成优化的驱动力，促进系统达到优化的稳态．生态平衡机制中，生物间的牵制关系是形成生态平衡的驱动力：各种生物具备自繁衍能力，能扩大种群规模，而各种生物的生存要依赖于其他生物，这种依赖性控制了种群规模，使得各种群之间的规模存在相互牵制的关系，从而达到一种动态平衡，而动态平衡点的位置取决于自然资源的规模．例如由狼群、羊群、青草所组成的简单生态系统，最终达到生态平衡时各种群的规模取决于这个生态系统中土地的面积、青草的覆盖面积及自然生长速度，即自然资源约束决定了最终的平衡点位置．可见生态平衡机制是一种以种群规模为设计变量，以自然资源为约束条件，以系统达到稳态平衡为优化目标的优化机制．1.2　牵制平衡算法机制牵制平衡算法的基本思路是将设计变量对应于生态系统中各种群的规模，构造单一设计变量的自成长函数和设计变量之间的牵制函数，以各设计变量的值达到容差范围内的稳态为目标，形成迭代机制．在此定义如下概念．自成长函数：在无约束的情况下，设计变量随着时间的变化规律函数．以动物种群规模的成长为例，设无约束条件下种群增长率为r，设t时刻时种群规模为x(t)，则在无约束的条件下，该种群在Δt时间内增长的个体数量为x(t+Δt)=x(t)+rx(t)Δt．(1)牵制函数：表达设计变量之间相互影响和牵制的依存关系函数．依存关系存在于种群之间，也存在于种群和资源之间，不但包括食物链上游对本环节的捕食消耗作用，还包括食物链下游对本环节限制或促进作用，二者共同构成了牵制关系．成长函数：综合考虑生物自成长和生物间的牵制关系而形成的函数，该函数既包含自成长函数，又包含牵制函数，二者共同组成了描述种群规模随着时间变化的函数．为了将上述概念及求解机制表述清楚，本研究以一个极简的生态系统(草、羊、狼组成的生态系统)为例来进行阐述．1.2.1　草草在没有其他物种食用的自然环境中，按照一定的速度自由生长，当生长到一定阶段时，由于对土壤资源的竞争而形成自阻滞作用，使得自成长规模也不是无限的，其自成长函数为PG(t+Δt)=PG(t)+[1-G(t)/NG]r1AΔt，(2)式中：PG(t+Δt)为草在自成长条件下t时刻的种群规模；G(t)为草在t时刻的实际种群规模；NG为草的环境容纳上限；r1为草的固有增长率；A为草生长的面积；G(t)/NG反映随着草规模的增长，由于土地资源的限制导致的自阻滞作用，即增长速度随着规模的增加而放缓．以上模型为假设没有食草生物存在的情况，而在“狼-羊-草”系统中，由于羊-草之间的食物链关系，草的规模将随着羊群扩张而相应减少．草随着羊群规模增加而减少的函数即牵制函数为QG(Δt)=r2S(t)Δt，(3)式中：QG(Δt)为草的规模在Δt时段的损失量；r2为单位规模羊群对草消耗速度；S(t)为t时刻羊群规模．因此，考虑到种群内自阻滞、物种间捕食关系及草作为生产者的特性，构造在t时刻草的种群规模函数，即成长函数为G(t+Δt)=G(t)+[1-G(t)/NG]r1AΔt-r2S(t)Δt. (4)1.2.2　羊羊在资源有限且没有天敌的自然环境中，其规模增长趋势可用如下所示的自成长函数进行说明，PS(t+Δt)=PS(t)+[1-S(t)/NS]r3S(t)Δt，(5)式中：PS(t+Δt)为羊在自成长条件下t时刻的种群规模；NS为羊群的环境容纳上限；r3为羊群的固有增长率；S(t)/NS反映随着羊群规模的增长，由于食物等环境资源有限而导致的羊群内部的自阻滞作用．显然，当羊群规模不断增加，逐渐接近其环境容纳上限时，增长曲线趋于平缓．区别于假设自然资源无限的J型增长模型，式(5)描述的是资源有限条件下生物增长趋势，如图1的S型曲线所示即逻辑斯蒂增长模型．10.13245/j.hust.239409.F001图1生物增长两类模型由于草-羊-狼三级捕食关系的存在，羊作为中间环节，受到草和狼两种生物的影响：一方面，草为羊群提供能量，因此草一定程度上促进羊群规模的增长；另一方面，狼通过捕食羊群获得能量，因此狼对羊群规模存在消耗作用．二者牵制关系为QS1(Δt)=r4G(t)Δt；(6)QS2(Δt)=r5L(t)Δt，(7)式中：QS1(Δt)为羊群规模在Δt时段的增量；QS2(Δt)为羊群规模在Δt时段的损失量；r4和r5分别为单位规模草对羊群的补给速度和单位规模狼群对羊群的消耗速度；L(t)为狼在t时刻的种群规模．考虑到种群内自阻滞、羊作为食物链中间环节受到草的促进和狼群消耗三种因素的影响，构造羊群t时刻的种群规模即成长函数为      S(t+Δt)=S(t)+[1-S(t)/NS]∙r3S(t)Δt+r4G(t)Δt-r5L(t)Δt. (8)1.2.3　狼狼在食物源充分的自然环境中，其种群规模主要与自阻滞作用相关．由于狼作为“狼-羊-草”食物链中最高一级生物，没有对其捕食、消耗的天敌，因此，除了自阻滞作用外，还应将狼群的自然死亡率作为影响种群规模自成长的因素，可用如下自成长函数表示PL(t+Δt)=PL(t)+[1-L(t)/NL]∙r6L(t)Δt-aL(t)Δt, (9)式中：PL(t+Δt)为在自成长条件下t时刻的狼群规模；NL为狼群环境容纳上限；r6为狼的固有增长率；L(t)/NL为狼群内部自阻滞作用；a为狼群自然死亡率．在“狼-羊-草”组成的捕食关系链中，羊群对狼群规模增长具有促进作用，可用如下的牵制函数表示为QL(Δt)=r7S(t)Δt，(10)式中：QL(Δt)为在Δt时段狼群规模增量；r7为单位规模羊群对狼的补给速度．因此，考虑狼群内部自阻滞作用、自然死亡率及羊群的促进作用，狼的成长函数构造为     L(t+Δt)=L(t)+[1-L(t)/NL]⋅r6L(t)Δt+r7S(t)Δt-aL(t)Δt. (11)1.2.4　迭代规则通过模拟自然界生物间巧妙的平衡机制，本研究设计了式(2)～(11)用于描述狼、羊、草三类物种间的牵制关系及生长规律，从而综合得出三类物种的成长函数，即物种经过Δt时间的种群规模增长量．算法通用的求解流程与之相同：首先，输入环境容纳上限等相关参数；然后，根据问题设计各种群规模的自成长函数和牵制函数，并计算当前实际种群规模；每经历一次迭代，各类种群规模都将依据成长函数进行更新；经过一定时间的迭代，当各种群规模的变化都在某一容差范围ε内时，算法停止，将小幅波动的中心值输出即为问题最终的解．2 算法测试2.1　收敛性实验为验证本文所提算法的收敛性，分别对草、羊、狼三类种群取不同的初始规模，考察最终所得结果是否收敛，如表1~3所示．在本组实验中，通过查阅相关生物种群的生长特性，设置各参数的取值如下：NG，NS，NL分别取400，200，100；草的面积A取300；r1～r7分别取0.3，0.2，0.2，0.1，0.3，0.1，0.02；d取0.05；最大迭代次数取100．此处参数取值仅改变各生物种群规模最终收敛结果值的大小，对于算法的收敛性并不影响，即在“草-羊-狼”生态系统中，两个核心函数的参数取值可以自行设定，这并不影响算法的求解能力和收敛性能．10.13245/j.hust.239409.T001表1各物种不同初始规模对应收敛结果(草/100/50)实验编号草的初始种群收敛值草羊狼1100238.2182.090.32300238.2182.090.310.13245/j.hust.239409.T002表2各物种不同初始规模对应收敛结果(200/羊/50)实验编号羊的初始种群收敛值草羊狼350238.2182.090.34150238.2182.090.310.13245/j.hust.239409.T003表3各物种不同初始规模对应收敛结果(200/100/狼)实验编号狼的初始种群收敛值草羊狼520238.2182.090.3660238.2182.090.3在此次收敛性实验中，通过控制变量，每次实验确保其中的一类物种规模变化，而另外两类物种的种群规模不变，以此研究各物种的初始规模大小对结果的影响．如表1所示，“(草/100/50)”表示每次实验中草的初始种群规模取值按照本列所示的规律变化，而羊群和狼群的初始规模固定，分别取100和50．显然，由表1~3可知：种类相异、数量不同的初始种群规模在实验中皆收敛至相同结果，即最终达到一致平衡状态，表明算法具有较好的收敛性．图2为实验编号1，3，5对应的收敛曲线．10.13245/j.hust.239409.F002图2迭代过程曲线在图2中，虽然迭代后期各种群规模看似是一条平滑的曲线，但实际上种群规模一直处于小幅波动中，只因迭代后期波动幅度小，故看似为平滑曲线，即在容差范围内形成了动态平衡．2.2　算法求解性能测试2.2.1　不同基础资源测试在狼-羊-草生态系统模型中，草的面积A为基础资源，是系统的约束条件．本研究通过设置不同A值，以探究基础资源约束对于优化结果的影响．算法的其他相关参数设置见2.1节，在此分别设置生长面积A为100，200，300，最大迭代次数为300，求解结果见表4．由表4可知：求解过程收敛且收敛结果相异，基础资源越充分，稳态收敛状态下各种群规模越大．可见基础资源决定了最终的求解结果，草的生长面积作为基础资源，不但对草的种群规模产生直接影响，且通过各物种间的牵制关系，将影响传递给了整个食物链，从而影响了最终的求解结果．因此，基于基础资源会对食物链中各种群规模产生影响这一特征，草的面积A可视作问题求解中的约束条件，决定最终的平衡值．不同的面积取值A收敛趋势如图3所示．10.13245/j.hust.239409.T004表4不同生长面积对应的收敛结果A草羊狼100129.6101.476.5200191.1156.786.3300238.2182.090.310.13245/j.hust.239409.F003图3不同基础资源条件下算法收敛图2.2.2　多物种求解测试“狼-羊-草”捕食系统中，物种数为3，在此基础上进一步引入“虎-狼”捕食关系，使得食物链中的营养级达到4级，以探究当更多生物种类存在时，该模型能否获得较好的收敛效果．实验相关参数见2.1节，另外在“草-羊-狼-虎”实验中补充虎的环境容纳上限为80，虎对狼的消耗速度和狼对虎的食物供给速度两个参数分别为0.2和0.02，最大迭代次数为100，求解结果见表5．10.13245/j.hust.239409.T005表54级捕食系统的求解结果(200/100/50/虎)虎的初始种群求解结果草羊狼虎10248.9170.0100.064.720248.9170.0100.064.730248.9170.0100.064.7求解迭代过程如图4所示，实验结果显示：即使增加了食物链的营养级，当对多物种捕食系统进行求解时，牵制平衡算法仍有较好的收敛效果．10.13245/j.hust.239409.F004图4物种数为4的捕食系统迭代过程2.3　算法复杂度分析在牵制平衡算法研究的多个物种构成的生态系统中，根据物种间的相互牵制作用，每进行一次迭代，则对物种的种群规模完成一次更新．假设在某一生态系统中存在M个物种，算法迭代N次停止，则牵制平衡算法的复杂度计算为：经过物种间的牵制关系对M个物种的种群规模进行更新，更新操作共须迭代N次；因此，算法的时间复杂度为O(MN)，空间复杂度为O(M)，皆在可接受范围内．3 案例分析为了验证本文所提牵制平衡算法的性能，选取了3个工程设计问题进行仿真实验，通过与现有文献的比对，阐明牵制平衡算法的有效性和竞争力．3.1　案例1(三杆桁架设计问题)三杆桁架在工程中较为常见，经常被用于测试优化方法的有效性和优越性[13-17]，其结构如图5所示[18]．在三杆桁架结构中，三杆截面积分别为X1，X2，X3，其中X1=X3．节点处受力为P，杆的应力约束为σ，在满足题设约束的条件下，确定横截面积使得三杆体积最小的数学模型为10.13245/j.hust.239409.F005图5三杆桁架结构图min  f(x)=(22X1+X2)l, (12)s.t.     g1=[(2X1+X2)/(2X12+2X1X2)]P-σ≤0;    g2=[X2/(2X12+2X1X2)]P-σ≤0;g3=[1/(X1+2X2)]P-σ≤0, (13)式中：l=100cm;P=2 kN/cm2；σ=2 kN/cm2；0≤Xi≤1 (i=1,2)．3.1.1　求解方法根据牵制平衡算法的求解机制，在此对三杆桁架问题中的各项参数与算法的对应关系进行说明．a. 环境容纳上限：即生物种群规模扩张的上限．因两个设计变量的取值范围为[0，1]，故设计变量相对应的生物种群规模取值上限为1，即环境容纳上限为1．b. 面积：即影响生物种群规模扩张的基础资源．问题中包含了三个约束条件(见式(13))，因为每个约束条件都要被设计变量满足，而约束条件的边界将影响设计变量的最终取值，故可将三个约束条件视作三种环境资源，所有的生物种群规模都受到三类资源的牵制．c. 物种数量：即生物种群的个数．物种个数相当于问题的设计变量个数(本案例中包含2个生物种群)，物种规模大小相当于设计变量的取值．d. 增长率、死亡率、消耗速度等细化参数：主要用于描述种群规模的变化规律，然而结合工程设计问题的特性，此处自成长函数并非照搬式(1)～(11)的核心函数，而是在牵制平衡算法机制的基础上引入了随机扰动因子，故以随机扰动代替增长率、死亡率、消耗速度等细化参数．结合上述参数设置规则，对自成长函数和牵制函数的构造说明如下．a. 三杆桁架问题中设计变量X1和X2的取值对应于牵制平衡算法中两个物种的种群规模，由于二者取值范围皆为[0，1]，因此当构造自成长函数时，采用二进制编码方式表示设计变量的取值以满足环境容纳上限约束．为达到1×10-6的编码精度，确定编码长度为20位，两个物种的自成长函数以如下所示的随机扰动完成对种群规模的更新，X1=[x1,x2,…,x20];X2=[y1,y2,…,y20],xi,yi=0  (rand(0,1)∈(0.0,0.5));1  (rand(0,1)∈[0.5,1.0)),i=1，2，…，20．b. 与三杆桁架问题相对应，种群间的牵制关系主要由约束条件体现，在此将题中包含的3个约束条件视为环境中的3类基础资源，在此基础上2个种群互相牵制和影响，争夺有限的环境基础资源．每次迭代须保证当前种群规模取值满足约束条件；否则，应根据牵制函数对两个物种的种群规模进行调整，使之满足案例约束．在这个过程中两个物种并非独立变化，二者相互牵制和影响，一个种群规模的改变往往导致另一种群规模随之变化，直至满足约束条件．两个物种的牵制函数见式(13)．3.1.2　结果分析牵制平衡算法在20次独立运行中，对于三杆桁架问题获得的最优设计方案如下：X1=0.788 99，X2=0.407 37，g1=-8.0×10-6，g2=-1.465 10，g3=-0.534 91，fmin=263.897 02．所得结果与基于邻域的连续优化一致性(neighborhood-based consensus for continuous optimization，NCCO)算法[13]、正弦余弦算法(sine cosine algorithm，SCA)[14]、乌鸦搜索算法(crow search algorithm，CSA)[15]、矿井爆破算法(mine blast algorithm，MBA)[16]和社会与文明算法(society and civilization algorithm，SC)[17]等算法求解结果进行对比见表6．10.13245/j.hust.239409.T006表6多种算法求解三杆桁架设计问题结果对比算法最优解最差解平均值标准差NCCO279.724 50279.758 00279.734 101.06×10-2SCA279.731 20282.842 70281.405 901.56CSA263.895 84263.895 84263.895 841.01MBA263.895 85263.915 98263.898 003.93×10-3SC263.895 85263.969 76263.903 361.30×10-2牵制平衡263.897 02264.737 01264.103 060.26由表6可知：与其他算法相比，本文提出的牵制平衡算法在三杆桁架设计问题中取得较好求解效果，其优化结果明显优于基于邻域的连续优化一致性算法和正弦余弦算法，在最优解、最差解和平均值三个方面的优化程度分别超过5.66%，5.37%和5.59%．当然，牵制平衡算法与乌鸦搜索算法，矿井爆破算法和社会与文明算法求解结果相比有微小差距，但在最优解、最差解和平均值三个方面的优化程度与最优方案的差距不足4.464×10-6，0.318%和0.080%．实验结果表明：牵制平衡算法在三杆桁架设计问题的求解中，能获得较满意的求解效果，具有较好的有效性和可行性．由图6所示的求解三杆桁架问题的收敛曲线可知，X1和X2在求解过程中，二者之间相互牵制以满足题设约束条件，同时使得目标函数不断向最优方向靠近，最终在300代以内达到收敛状态，获得最优解，显示本算法具有较快的收敛速度和较高的求解效率．10.13245/j.hust.239409.F006图6牵制平衡算法求解三杆桁架问题收敛曲线3.2　案例2(管柱设计问题)管柱结构优化设计问题[19]如图7所示．采用牵制平衡算法对管柱结构设计问题进行求解，并将10.13245/j.hust.239409.F007图7管柱结构图所得结果与文献[20-22]中布谷鸟搜索算法(Cuckoo search algorithm，CS)、标准萤火虫算法(BVWFA)等算法求解结果进行对比(见表7)．10.13245/j.hust.239409.T007表7多种算法求解管柱设计问题结果对比参数文献[22]文献[19]CSBVWFA牵制平衡算法d5.450 705.440 005.451 395.451 165.460 07δ0.292 000.293 000.291 960.291 970.291 49g1-7.800×10-5-0.857 90-0.024 10-1.618×105-3.247×10-3g20.131 700.002 60-0.109 50-1.765×105-3.650×10-6fmin25.531 6026.532 3026.532 1726.499 5126.517 38由表7可知：本文提出的牵制平衡算法在管柱设计问题中能获得较好求解效果，所得结果优于布谷鸟搜索算法及文献[19]所提算法，稍劣于标准萤火虫算法，但与其差距不足0.067%，在可接受范围内．另外，虽然文献[22]求解结果优于当前文献记载数据，但其所得结果超出约束条件g2，因此很难将其视为完美解．实验结果表明牵制平衡算法对于管柱设计问题的求解具有较好的有效性和可行性，其收敛曲线如图8所示．10.13245/j.hust.239409.F008图8牵制平衡算法求解管柱设计问题收敛曲线3.3　案例3(压力容器设计问题)压力容器结构设计问题[23]如图9所示．使用牵制平衡算法对压力容器设计问题进行求解，结果如下：Ts=0.814 33，Th=0.406 76，R=42.122 83，L=176.737 76，g1=-0.001 36，g2=-0.004 91，g3=10.13245/j.hust.239409.F009图9压力容器结构图-2 247.18，g4=-63.262 2，fmin=5 981.841．与粒子群优化算法(particle swarm optimization algorithm，PSOA)[15]、状态转换算法(state transition algorithm，STA)[24]、混沌灰狼优化算法(chaotic grey wolf optimization algorithm，CGWOA)[25]、乌鸦搜索算法[15]及改进的乌鸦搜索算法(modified crow search algorithm，MCSA)[15]等算法求解结果对比见表8．10.13245/j.hust.239409.T008表8多种算法求解压力容器设计问题结果对比测试算法最优值最差值平均值标准差PSOA6 693.721 2014 076.324 008 756.680 301 492.567 00STA6 287.037 206 377.878 256 332.254 50150.883 80CGWOA6 134.180 006 388.110 906 159.320 00254.500 00CSA6 059.714 307 332.841 606 342.499 10384.945 00MCSA6 059.714 306 420.714 306 086.214 30104.407 10牵制平衡算法5 981.841 067 036.462 416 432.879 80326.433 00由表8可知：牵制平衡算法对于压力容器设计问题的求解取得最优结果，不仅优于传统粒子群优化算法和布谷鸟搜索算法，相较于状态转换算法、混沌灰狼优化算法和改进的乌鸦搜索算法等，其求解优势亦十分显著．实验结果表明牵制平衡算法对于压力容器设计问题的求解具有显著优越性和竞争力，其收敛曲线如图10所示．10.13245/j.hust.239409.F010图10牵制平衡算法求解压力容器问题收敛曲线4 结语仿生算法是求解复杂NP-Hard问题的有效途径，由于复杂问题的千差万别，目前没有公认的对各种问题都有效的算法，仿生算法也一直处于改进和扩充之中．本研究模拟大自然精妙简洁的生态平衡机制，提出了一种新的仿生算法——牵制平衡算法，并阐述了这一算法的两个核心函数(自成长函数和牵制函数)的构造原理与方法，提出了模拟生态平衡机制的牵制平衡算法机制．通过收敛性测试、不同基础资源测试、多物种求解测试等不同角度的算法性能测试，证明了算法的有效性和可靠性．以三个工程问题为比对实验案例，将牵制平衡算法的求解结果与现有文献中多种算法的求解结果进行比对，结果表明：与现有仿生算法相比，牵制平衡算法在求解工程设计问题中具有显著的优越性和竞争力，该算法既扩充了仿生算法的种类，又为求解复杂问题提供了新思路和新途径．当然，牵制平衡算法作为新提出的算法机制，其求解其他各类问题的有效性和优越性尚有待长期的实践检验．相信和其他仿生算法一样，牵制平衡算法也将在未来的应用实践中得以发展和完善．