不等圆Packing问题的多策略优化方法

梁利东，何东，朱良恒

（安徽工程大学 机械工程学院，安徽芜湖 241000）

不等圆Packing问题是一类典型的NP-hard问题[1]，具有组合优化问题的典型特征，广泛应用于纺织、造船、皮革工业、无线通信、航天器设计等领域[2-4]。作为图形Packing问题一个分支，研究和探索更优的排列模式和算法对于众多行业都具有重要的科研和经济效益[5-6]。随着计算机技术的飞速发展，计算机数值算法尤其是启发式算法来求解图形Packing问题成为了这一领域研究的主要方向。

针对圆Packing问题，国内外学者从各种角度提出了许多不同的求解策略。黄文奇等[7]将拟物算法与禁忌搜索策略相结合提出带全局变换禁忌搜索算法GT-TS，通过禁忌规则和特赦规则的约束实现当前局部最优布局的全局变化策略。Hifi等[8]基于左下角优先的策略求解矩形容器的不等圆Packing问题，提出了一种求解约束和无约束圆形Packing问题的启发式算法。Birgin等[9]首次提出了使用邻接矩阵来加速底层算法的概念，可实现将底层算法的计算复杂度从O(N2) 降 低到O(N)，大大缩短计算时间。王英聪等[10]利用群智能劳动分工的任务分配来实现不等圆Packing问题的空间分配，从分配的角度分析不等圆Packing问题和群智能劳动分工方法；之后又提出一种基于位置选择的构造法[11]，即将圆放置过程看成位置选择过程，借鉴群智能劳动分工中的刺激-响应原理，将未布局空间的完整度和已布局空间的紧密度看作圆形物体选择位置时的刺激和阈值，其自适应阈值调整策略具有较好的可行性。Specht[12]设计并实现了一种基于图论方法分析容器中未被填充区域的算法，并将未被填充的区域称之为“孔”，通过探测“孔”的位置及其大小，结合跳转策略，有效地提高了填充密度。余丽娟等[13]通过改进基于Power图的区域划分，提出一种收敛速度更快的圆Packing算法，该算法是对局部Power图算法（Local power diagram，LPD）的一种改进。刘景法等[14]针对正三角容器内的等圆Packing问题，提出了将启发式格局更新策略与基于梯度法的局部搜索策略相融合的模拟退火算法，并与二分搜索结合，实现能量更低的最优化格局，算例证明了其有效性。张维等[15]提出一种改进遗传模拟退火算法，通过计算生成一个合适大小的初始圆形容器来指导初始种群的生成，采用最优保存策略来保证历代的最优解，并结合遗传算法和模拟退火算法分阶段来提升局部搜索和全局搜索性能。何琨等[16]提出了基于动作空间的拟物求解算法，借鉴了求解矩形Packing问题中动作空间的概念，通过化“圆”为“方”，将不规则的空闲空间近似表示为一系列矩形空间，有效地解决了跳出局部最优的难点。

本文针对不同类型的拟物算法进行验算和性能分析，以拟物算法为基础，综合分支搜索策略、多重退火策略提出多重退火分支搜索算法（Multiple annealing branch search，MABS）。针对连续优化和组合优化两个层面描述不等圆Packing的最优化策略和算法思想，并通过国际算例分析验证了不同形状容器下的实验数据，证明该算法是有效可行的。

1 问题描述及算法模型

不等圆Packing问题是将一组已知个数的大小不等的圆不重叠放置在给定形状的容器中，求容器的面积最小值或容器的填充密度最大值。以圆形容器为例，设容器半径为r0，其中放置的第i个小圆的半径为ri， 坐标为 (xi,yi)。假定有N个圆需要放置，则对于某一布局，数学模型表示为：

拟物算法的求解思想是将需要放置的各圆看作为实心光滑弹性体，而容器看作是在填满整个平面的无限弹性体中挖去相应部分而形成的光滑空腔。通过模拟各圆在容器中受弹性力作用产生挤压运动来实现连续优化，其中各圆所受的弹性力合力方向也就是使目标函数下降最快的负梯度方向，算法在迭代中通过调整容器大小以满足布局中绝对弹性势能最小化的约束要求，从而求的最优解，算法的直观几何图像如图1所示。

本文的MABS算法将相对势能作为约束条件并通过变邻接系数加速方法构建算法模型，以定步长序列梯度下降算法为连续优化方法，采用领域算子、变长分支搜索以及多重模拟退火策略为组合优化策略实现不等圆排列优化布局，算法思想及总体结构如图2所示。

图2 本文算法的总体结构

相关描述和定义如下：

1）相对势能

在大量实验过程中发现，对于不同尺寸大小的布局图形，降低其绝对势能的难易程度存在很大区别，主要体现在：较小的圆半径和容器尺寸更容易得到较小的势能，而对于拥有较大圆半径和容器尺寸的算例，使势能降至同一水平的难度明显增大。其原因在于绝对势能是建立在图形间的绝对挤压深度上，当图形整体扩大至原来的k倍时，则所有挤压深度也扩大至原来的k倍，根据绝对势能的定义，绝对势能将扩大至原来的k2倍。因此为了消除图形绝对尺寸的影响，本文将容器的尺寸进行归一化预处理：将尺寸为L的容器大小缩小到原来的1/L，其势能相应缩小到原来的1/L2。考虑到相对势能的数值很小，不利于对比观察，可将其扩大106倍。以圆形容器为例，对于任一布局，相对势能为约束的问题模型可表示为

图3 相对势能与绝对势能对比

2）变邻接系数加速方法

对于n个不等圆，拟物算法在每轮迭代都需计算n（n-1）/2个距离和每个圆所受的矢量力以及弹性势能。然而从直观的几何图像来看，某个圆所受的矢量力显然只与相邻圆有关，因此可通过构建邻接矩阵来加速搜索并降低计算复杂度，图4为使用邻接矩阵前后的所需计算距离表示。

图4 使用邻接矩阵前后对比

邻接矩阵是表达图形间接触关系的矩阵，lij为各圆的邻接关系属性值，dij表示为各圆之间的距离，定义如下：

其中当邻接系数kn= 1时，邻接矩阵为严格邻接矩阵，客观表达各圆的邻接关系；当kn＞ 1时，邻接矩阵定义为过定邻接矩阵，即不仅包含已有的邻接关系，而且存在潜在邻接关系，极有可能在迭代后期产生邻接关系。为保证邻接关系变化的相对稳定性，尝试在迭代的初期使用较大的kn以包含更多潜在的邻接关系，同时可减小邻接矩阵的计算间隔以应对邻接关系的频繁变化。这样，就可以依据迭代不同阶段的特点通过调整kn与间隔代数以实现全计算周期的邻接矩阵加速。实验也表明使用变邻接系数邻接矩阵可将算法的计算复杂度从O（N2）降低至O（N），大大缩短计算时间。

3）定步长序列梯度下降算法

根据实现机制和策略的不同组合，可衍生出众多拟物算法[17-18]。本文以定步长序列梯度下降算法、定步长同步BFGS算法和变步长同步梯度下降算法为例进行测试对比和分析。为了避免单次求解的随机性干扰，所有数据都是基于随机生成的初始布局运行10次所取的平均值。图5为3种算法的运算时间和收敛势能对比图，结果表明：定步长序列梯度下降拟物算法具有最高的运算速度，其次是定步长同步BFGS拟物算法、变步长同步梯度下降算法，势能收敛优劣情况则正好相反。

图5 3种算法的性能比较

实验表明，尽管变步长同步梯度下降拟物算法在对随机初始布局的势能收敛性方面表现最优，但考虑到后期的组合优化对布局扰动的稳定性因素，本文又特别模拟实际运算场景测试了变步长同步梯度下降拟物算法与定步长序列梯度下降算法对经过扰动的稳定布局的势能收敛性情况，结果如图6所示。可见，定步长序列梯度下降算法在实际运算场景的多数情况下势能收敛性接近变步长同步梯度下降拟物算法，而且又拥有最快的运算速度，故采用定步长序列梯度下降算法作为MABS算法的连续优化模块。

图6 实际运算场景下两种拟物算法势能收敛性对比

2 算法的组合优化策略

2.1 邻域算子

不等圆Packing的最优布局在很大程度上取决于各圆的排列位置。多起点策略由于所生成的各个布局之间没有任何联系，在生成新布局时也无法利用之前生成的布局信息，搜索盲目而低效。因此运用邻域算子在不大幅改变布局的前提下，可通过对布局进行小的变动而得到更优布局。本文采用3种邻域算子：1）抛掷算子：随机挑选一个圆并随机“抛掷”到容器内某处；2）交换算子：交换当前布局中任意两个不等圆的位置；3）振荡算子：将所有圆在各自的随机方向上做小幅度移动。图7为抛掷算子与交换算子生成新邻域的过程示意。

图7 利用交换算子与抛掷算子生成新邻域的过程

相关定义如下：

定义2（邻域算子）：对稳定布局进行一定变换，得到与之不同的稳定布局，这一变换称之为邻域算子。

邻域表征的是稳定布局之间的关系，借助于邻域算子，有可能对布局进行累计的小幅度改进以获得最优布局。此外，抛掷算子、交换算子与振荡算子的多次叠加使用也可以看作是一个邻域算子。显然，叠加次数越多，所产生的布局就可能更不同于原有布局。换言之，可以认为在一定范围内，3种算子的叠加次数与初末布局之间的相似度呈负相关。

2.2 变长分支搜索策略

分支搜索策略源于优胜劣汰思想，即从当前最优的父节点出发，使用邻域算子生成一系列子节点，并比较子节点与父节点的势能。若父节点势能低于子节点势能，则继续从父节点生成新的子节点；否则将此子节点更新为父节点，并从新的父节点出发继续搜索其更优子节点。最优布局的定义如下：

分支搜索策略总是从最优节点生成新的节点，而生成的新节点只能经过一次邻域算子的处理，即分支长度为1的分支生长过程。由于分支较短易于陷入局部最优布局，通过延长分支长度来扩大搜索范围，变长分支搜索策略的算法流程如图8所示。

图8 变长分支搜索策略流程图

实验证明延长分支长度确实可以有效提高搜索质量，但随着分支长度的增长，计算耗费的时间也同比增长，因此需要在计算时间和求解质量之间做出权衡。本文对分支长度从1至3的分支搜索进行了实验对比，实验结果表明分支长度为2的分支搜索在耗时和求解质量两项指标上均取得了较好的表现，如图9所示。

图9 不同长度分支搜索性能对比

2.3 模拟退火接受准则及多重退火策略

为进一步提高布局多样性，本文基于模拟退火算法，并以一定概率接受更差解的策略来达到全局最优的目的。在退火过程中，假定某布局Xi，经过邻域算子扰动生成新布局Xj，两布局的势能分别是fr(Xi)和fr(Xj)。根据Metropolis准则，接受新布局Xj的概率表示为：

式中T为温度，当迭代次数增加时取Tk+1=α·Tk，α=0．99。 当fr(Xj)＞fr(Xi) 时 ，有可知此时Pj∈(0,1)， 且温度越低Pj越接近于0。这表明随着迭代进行，模拟退火接受更差解的概率逐渐趋于0，算法流程如图10所示。

图10 模拟退火算法流程图

实现发现，200 ℃以上的初始温度对解的质量增长不大，但迭代次数和运算时间却随着初始温度升高而大幅增加。因此考虑到求解质量和运算时间之间的平衡控制，本文依据实验数据采用400 ℃的退火初温，并与贪婪算法进行比较，实验数据如图11所示。

图11 不同初始温度下模拟退火算法的求解质量对比

由于模拟退火算法跳出局部最优的能力主要体现在搜索前期，当进入势能较为稳定的搜索中后期，找到更优解或跳出当前局部最优解的可能都很小。对此，提出多重退火的改进思路，即当判定模拟退火陷入局部最优布局或缺乏找到较优解的潜力时，及时跳出本次退火进程，通过升高温度开始新一轮退火过程。当温度重新升高，接受差解的概率将大大增加，搜索随机性增强，相当于对当前布局进行大幅重组，从而可跳到新的搜索域，以期望获得可能的更优布局。当达到限定时间或限定重退火次数即可结束搜索，图12为多重退火策略实施的势能和温度的变化曲线。

图12 多重退火势能温度变化曲线

可见，在时间相同的情况下，通过增加多重退火判断准则（势能下降量准则），多轮模拟退火搜索在局部最优布局上或缺少搜索前景的布局上极大降低了计算资源，同时实现了在不同状态空间的寻优过程，既扩展了搜索广度，也提高了最终求解的质量。

3 实例分析

本文的MABS算法基于MATLAB软件进行编写，用于实验的硬件环境为：笔记本电脑，CPU为四核Intel酷睿i5-7400，主频3.0 GHz，运行内存16 GB。算法在6种不同形状的容器中进行不等圆排列布局优化，取N=20，ri=（i=1,2,···,N），获得的最优布局如图13所示，其中φ 为填充密度。

图13 MABS算法在不同形状容器中的最优布局

其中，针对packomania官网中圆形容器与正方形容器的对应算例，本文算法所得最优解与其最佳记录的排样密度持平。其他4种形状容器的不等圆Packing算例未见记录，所得填充密度为本算法独立得出。

表1 圆形容器中26个国际公开算例的实验结果

由此可知，MABS算法在近一半的算例上与之前最优解持平，10个算例的解略差于记录最优解，改进了4个算例的此前最优解。可见MABS是具有较高性能的不等圆Packing算法。

本文算法基于最小包络圆规则也可将不规则图形Packing问题转化为不等圆问题进行求解。图14表示为给定的4个不规则图形及其形成的最小包络圆，图15为MABS算法所得的最优布局。

图14 不规则图形及其最小包络圆

图15 不规则图形Packing求解应用

4 结束语

本文基于拟物算法模型提出了分支搜索策略、多重退火策略的高效不等圆Packing算法，在连续优化和组合优化两个层次对算法进行了优化和改进。以定步长序列梯度下降拟物算法作为底层排列方法，并在运算迭代的不同阶段使用不同的邻接系数及计算间隔计算邻接矩阵，改进了邻接矩阵加速方法；采用分支搜索策略，通过增加分支长度和领域算子来扩展搜索空间以提升解的质量；针对模拟退火算法迭代后期难以跳出局部最优的缺陷提出了多重退火策略，并设置接受概率保证全局优化过程。通过对国际公开算例集的实验结果以及在不规则图形排列问题的求解应用表明，该算法具有一定的有效性。