Packing问题属于NP难的组合优化问题，不等圆Packing问题是其中的一个子问题，简单描述就是把n个半径不等的小圆放入一个圆形容器中，在给定小圆个数和小圆半径的情况下使得圆形容器的半径值最小。Packing问题在集装箱装货、下料排样、卫星舱布局等领域具有重要的理论意义和应用价值，例如计算一个集装箱的装载容量，可以估算出需要多少个集装箱来装货；再例如卫星舱布局需要确定零件的安装位置，可以先把类似圆的零件用最小包络圆替代，进而计算零件的摆放位置。

使用传统的数学方法无法对Packing问题进行求解，因为这是NP难问题。因此，人们提出了使用启发式算法来求解此类问题。针对此问题使用的启发式算法有遗传算法、退火算法以及一些混合启发式算法。黄文奇等^[1-2]提出了一个拟物算法，这是一种基于弹性势能模型的算法，根据此算法得到了多个最优结果。刘朝霞和刘景发^[3]提出了一种模拟退火算法，得到了更快的求解速度。陆一平和查建中^[4]提出了一种膨胀方法，对装填物进行合适的膨胀，通过排斥操作来求解，对多种几何形状的装填物有较好的适应性。康燕^[5]基于弹性力学模型使用禁忌搜索的启发式算法求解不等圆Packing问题，该算法首先把小圆按给定的优先级进行分组，然后逐组地使用拟物拟人法来放置小圆，在整个过程中通过禁忌搜索的思想即禁止重复之前的工作来使得搜索逃离局部极小值的陷阱，提高了算法效率。

本文采用遗传模拟退火算法来解决不等圆Packing问题，遗传算法的优点是全局搜索能力强，缺点是容易早熟收敛，局部搜索能力差；模拟退火算法可以以一定概率选择次优解，从而避免陷入局部最优解，其优点是局部搜索能力强，缺点是全局搜索能力差，这样一来模拟退火算法和遗传算法的优缺点恰好互补，因此采用遗传算法和模拟退火算法的混合算法来求解不等圆Packing问题。

1 数学模型

为解决不等圆Packing问题，建立其数学模型，该问题的数学描述为：在二维平面内，设有n个半径为r_i的小圆(r_i为变量)，把这n个小圆放入圆形容器内，在满足约束条件的情况下使圆形容器的半径值R尽量小，求解小圆的圆心坐标。

约束条件为：①小圆与圆形容器没有嵌入；②任意2个已放入圆形容器的小圆彼此互不嵌入。

求解该问题的数学模型为：将二维笛卡尔坐标系的原点取在圆形容器的圆心处，令其半径为R(半径值R可以转化为所有小圆中圆心到圆形容器圆心的距离值加该小圆半径值中的最大值)。小圆的圆心坐标为(x_i, y_i), i=1, 2, …, n, 问是否存在2n个实数x₁, y₁, x₂, y₂, …, x_n, y_n, 使得：

(1)

(2)

(3)

2 一种改进的遗传模拟退火算法

模拟退火算法是一种概率搜索算法^[6-7]，适合于在一个较大的搜索空间内寻找最优解。该算法的原理是对热力学中固体退火过程的模拟，开始给定一个较高的初始温度，然后缓慢下降温度，使得算法能够在多项式时间里给出一个近优解。模拟退火算法每次从当前解的临近解空间中选择一个最优解作为当前解，但该算法也以一定的概率接受一个比当前解更差的解，因此就有一定概率跳出局部最优解而达到全局最优解。模拟退火算法具有良好的局部搜索能力，可以快速搜索到局部的最优解，但是模拟退火算法是单点搜索，不能保证搜索过程进入最优解所在区域进行搜索，这样就导致模拟退火算法在前期搜索效率不高。

遗传算法是一种具有全局搜索能力的启发式算法，具有很广泛的应用范围，但是遗传算法也有依赖初始种群、容易早熟收敛的缺点。遗传算法是群体进化，优点恰恰在于全局搜索能力强，而局部搜索能力弱，如此一来，结合遗传算法和模拟退火算法的混合算法应运而生^[8-10]。

2.1 算法流程图

本文改进遗传模拟退火算法的流程是：首先随机产生一组初始解，然后对初始解进行选择、交叉、变异操作来产生一组新的个体，然后分别对每个个体进行模拟退火操作，使用得到的每个新个体来组成下一代群体。具体流程图如图 1所示。

算法描述：

1) 根据公式(4)可以得到算法的初始值R；

2) 随机生成满足条件的个体，组成初始种群；

3) 对种群中的每一个个体适应度值进行计算。计算当代种群的平均适应度值，和上一代种群的平均适应度值做比较，若两者差值小于ε，则满足收敛条件，算法运算结束，转步骤11)，否则进行迭代。

4) 计算模拟退火算法的初始温度；

5) 保存初始种群中的最优个体，避免最优个体没有被选中；

6) 采用轮盘赌选择方法进行选择；

7) 进行基因重组操作，保存重组后的最优个体，避免最优个体在变异操作中丢失；

8) 进行基因变异操作，保存变异后的最优个体，避免最优个体在退火操作中丢失；

9) 进行模拟退火操作；

10)转步骤3)；

11)对计算结果进行判断。若结果满足约束条件，则给R赋新值R-a，a为强制迭代的步长，然后再次运行算法；

12)若算法运算失败，则说明在当前R值的约束下已无法求解，应输出上一步结果，说明上一步结果是本算法求到的最好解。

2.2 圆形容器初始R值的计算

根据给定的小圆个数，可以把小圆按照图 2的方式排列，然后用一个大圆把图中的正方形包络起来。以图 2为例，图中有4行4列共16个半径为1的小圆，当小圆个数n满足条件9 < n≤16时，可以算出的圆形容器的半径即为图中所示，。

根据这种方法，计算出一个整数K，当半径为1的小圆个数n满足条件(K-1)² < n≤K²时，可以得到K行K列的小圆组成的正方形边长为2K，进一步得到大圆直径为，圆形容器的半径值为R=。同理，在不等圆中，假设有n个不等圆，半径值分别为r_i，找到半径最大的小圆r_max给出圆形容器初始半径值：

2.3 初始种群函数

种群的大小用N表示，即该种群中包含N个个体，种群的大小在每一次迭代中都不变。当N取值较小时，可以提高遗传算法运行速度，但却降低了种群的多样性；当N取值较大时，会使得遗传算法运行缓慢。针对本问题，每一个圆的坐标由x和y构成，所以该问题的变量个数为小圆个数的2倍，N的取值定为变量个数的10倍，当N值超过100时，限定N值最高为100。

由于该问题要求精度高，搜索空间大，因此适合采用浮点数编码。浮点数编码是指个体的每个基因值用某一范围内的一个浮点数来表示，而个体的编码长度等于其决策变量的个数。任一个体T_i的染色体是由x₁, y₁, x₂, y₂, …, x_n, y_n的浮点数形式组成。

初始种群是随机产生的，但是要满足问题描述中的2个约束条件，即小圆与圆形容器没有嵌入；任意两个已放入圆形容器的小圆彼此互不嵌入。

使用随机函数分别赋予x₁, y₁, x₂, y₂, …, x_n, y_n满足约束条件的值，组成个体T_i记为：

该种群中共有N个这样的个体，记为：

式中，S表示种群，T_i表示种群中的第i个个体。

2.4 适应度函数

不等圆装箱问题的目标函数是f(x), f(x)的值越小越好，在使用遗传模拟退火算法过程中，适应度函数的值越大越好，所以定义适应度函数为：

(4)

同时保存种群中最大的适应度值max_-fitness，最小适应度值min_-fitness，平均适应度值mean_-fitness。

2.5 计算模拟退火算法的初始温度

模拟退火算法的初始温度越高，获得高质量解的几率越大，但是花费的时间也越长。在此通过适应度函数计算种群中每个个体的适应度值，使用经验公式(5)可以给出一个合适的初始温度。

(5)

式中，pr为初始接受概率，本文定义pr的值为0.4。

2.6 保存最优个体函数

在遗传算法和模拟退火算法中，在寻优的每一步中，交叉环节、变异环节、退火环节等，只要是产生新解的环节中，都有可能破坏掉当前的最优解，导致最后的结果反而不如过程中的结果。所以在每次产生新解之前，都使种群中的最优个体保存下来不参与改变，然后用这个最优个体把新种群中适应度值最低的个体替换，从而保证产生新解过程中的最优解不被破坏。

2.7 交叉函数和变异函数

遗传算法中交叉概率的选取影响算法的性能和最终结果，交叉概率越大，产生的新个体就越多，但是优良的个体被破坏的可能性也越大；交叉概率过小，则产生的新个体过少，使得整个搜索过程变得缓慢。在此采用一种自适应调整交叉概率的方法。

在迭代过程中，如果当前种群的个体适应度值低于平均适应度值，则提高交叉概率；如果当前种群的个体适应度值高于平均适应度值，则降低交叉概率。交叉概率的自适应调整公式为：

(6)

式中，f_max为种群中最大的适应度值，f_mean为每代群体中的平均适应度值, f为要交叉的2个个体中较大的适应度值，本文中P_c1=0.9, P_c2=0.5。

交叉操作具体步骤为：

步骤1  假设种群大小为N，对群体中的个体进行随机配对，可以组成[N/2]对配对个体组；

步骤2  对于每一个配对个体组，产生一个介于0和1之间服从均匀分布的随机数u；

步骤3  比较配对个体组中较大的适应度值和种群的平均适应度值，利用公式(6)计算得到该配对个体组的交叉概率；

步骤4  如果交叉概率P_c大于u, 则对该配对个体组进行交叉操作；

同样的，变异概率P_m也影响算法的性能和结果，P_m过小就不易产生新的个体，P_m过大则变成了纯粹的随即搜索算法。在此也采用自适应的变异概率。变异概率的自适应调整公式为：

(7)

式中，f_max、f_mean定义同上，f为要变异个体的适应度值，本文中P_m1=0.1, P_m2=0.01。

2.8 退火环节

模拟退火算法的核心操作就是使用Metropolis接受准则来选择是否接受新的状态。Metropolis接受准则是：在温度T时，从之前状态i到新状态j，若E_j < E_i，则接受j状态为当前状态；若E_j>E，计算概率，生成一个[0, 1)区间的随机数rand，如果p>rand，则仍接受j状态为当前状态，反之保留i状态，其中K为玻尔兹曼常数。该步骤可以实现在高温下，接受与当前状态能量相差较大的新状态，当温度降低到低温时，只接受与当前状态能量相差较小的新状态。

模拟退火算法中用E表示不同状态的能量，可根据问题中的目标函数进行定义，算法过程中E的值不断降低；对于本问题，算法中用fitness表示个体的适应度值，算法过程中fitness的值不断变大，因此本问题的退火环节具体步骤与模拟退火算法有所不同，不同之处在退火环节的具体步骤中体现。

退火环节具体步骤为：

步骤1  把步骤2~5对种群中的每一个个体执行一次，总共N次，N为种群大小；

步骤2  计算当前个体的适应度值为fitness_-A;

步骤3  因为是浮点数编码，所以对每一个个体加上一个可正可负的随机数，来产生新解，该随机数为round[(rand-0.5)×2]×0.03，其中round是一个四舍五入的函数，rand是一个可以产生[0, 1)区间随机数的函数；

步骤4  计算新的个体的适应度值为fitness_-B；

步骤5  判断fitness_-A和fitness_-B的大小，若fitness_-B较大，则接受新的个体代替原来的个体；若fitness_-A较大，计算概率p=exp[-(fitness_-A-fitness_-B)/KT]，如果p大于随机数rand，则仍然接受新的个体代替原来的个体，否则保持原来的个体；

2.9 停止条件和约束条件的判断

遗传算法的停止条件多种多样，有使用运算时间作为停止条件的，有使用迭代次数作为停止条件的，有使用预期适应度值作为停止条件的等等。在此使用函数容限作为停止条件，即当代种群的平均适应度值相较上一代种群的平均适应度值没有明显变化，则算法停止。假设上一代种群的平均适应度值为avg_-F(x)，当代种群的平均适应度值为avg_-F(x)，如果满足公式(8)，则停止算法。

(8)

式中，ε=0.001。

圆形容器的R值不断缩小，总会有缩小到装不下这么多小圆的时候。计算结果表明，当圆形容器装不下这么多小圆时，小圆会重叠在一起，即违反“任意2个已放入圆形容器的小圆彼此互不嵌入”的约束条件。算法求解的结果是一组坐标值，可以根据公式(2)和公式(3)来判别计算结果是否满足约束条件，若满足约束，则说明满足条件，R值减a，继续迭代；否则输出上一步结果，说明上一步结果是该算法能求到的最小R值。

3 实验结果与分析

本文的目标是使圆形容器R值越小越好，在结果分析中把面积利用率作为指标之一，是因为面积利用率反映的结果更直观，圆形容器R值越小，面积利用率则越高。面积利用率η=(π·r₁²·n₁+π·r₂²·n₂)/π·R²。

本文采用的改进遗传模拟退火算法，根据工程中的实际要求计算了6个算例，对于每个算例分别调整了2种小圆的大小和数量，最终求解得到的圆形容器的面积利用率在72%左右，满足了工程应用的要求，说明了该算法的有效性。

表 2中展示了遗传模拟退火算法改进前后的求解结果，可以看出，未改进算法的求解结果出现了早熟收敛的现象，求得圆形容器的面积利用率仅在50%上下。

4 结论

本文提出了一种改进遗传模拟退火算法来求解不等圆Packing问题，算法通过给定初始圆半径来指导初始种群的生成，并通过一步步缩小圆形容器的半径值来优化初始种群的生成，使得遗传算法依赖初始种群的缺点成为了可以利用的优点；在算法中加入保存最优个体的策略有效避免了当代最优个体被破坏，使得算法终止时得到的结果是历代出现过的适应度值最高的个体；自适应的交叉概率和变异概率有效调整了较优个体和较差个体产生新解的概率，从而保证群体向着有利的方向进化；退火环节克服了遗传算法容易陷入局部最优的缺点，帮助种群中的个体跳出局部最优，进而得到更好的解。通过算例验证了算法的有效性，其计算结果满足了工程应用的需求。