拟人拟物法求解不等圆Packing问题

NP难度这门课还是比较有意思的，老师布置了一道作业，写一个用拟人拟物法求解不等圆Packing问题的小程序。

问题描述：在一个已知的容器中希望能放下N个不同形状大小的物体，其中界限容器的封闭边境以及各个物体都是不可入的刚性实体，如果客观上放不下，我们要求做出放不下的判断；如果客观上放得下，则要求给出每个物体的位置和方向。

这就是所谓的Packing问题，我们将问题简化了，只考虑圆形的容器和物体，给定容器的半径、物体的数量、物体的半径，求解是否存在一种布局，使得容器和物体两两之间都不相交。Packing问题是一个NP难度问题，因为我们要确定不存在满足条件的布局，必须穷举所有的布局，这是不可能的。

拟物方法是解决Packing问题的比较经典的算法。我们将这N个圆形物体和容器想象成光滑弹性的气球，容器是固定不动的，当这N个物体放入容器后，就会受到挤压而在弹力的作用下开始运动，最终有可能各个物体和容器都恢复自己的大小和形状，问题就得到了解决。因此我们引入以下概念：

• 任一时刻，称N个物体的圆心坐标x1，y1，x2，y2……，xn，yn为系统在此刻的格局；
  
• 物体间的距离Lij，如果i、j物体相切，Lij=0；如果i、j物体相离，Lij>0；如果i、j物体相交，Lij<0，且|Lij| = Ri + Rj - 圆心距离；
• 物体间的弹性势能Uij，当Lij>=0时，Uij=0；当Lij<0时，Uij=Lij^2；
• 物体和容器所构成系统的总弹性势能U，等于各个物体弹性势能之和；
• 物体移动的步长h，物体向当前时刻所受合力方向移动的距离；
• 梯度t，h梯度下降的速度。

很显然，U等于0代表着问题得到了解决。具体的拟物方法：1.初始化阶段，为物体随机产生一个格局；2.计算当前格局的总弹性势能，如果U=0退出算法；3.比较当前总势能和上一格局的总势能，如果势能上升了，h=h*t；4.计算每一个物体所受的合力向量，朝其合力方向移动h长度，获得新的格局；5.返回第二步。第三步的意思是梯度下降，当发现总势能上升时，说明格局一步跳大了越过了平衡点，这时应减小步长h。

但是纯粹的拟物方法有可能掉入局部最小值的陷阱，即各个物体和容器的受力并不为零，但是合力为零，物体将不会继续运动，这时就需要引入拟人方法。想象一下拥挤的公交车，车里面最想挪位置的是挤得最厉害的人，因此我们将当前挤压最严重的物体（姑且称其为“难民”吧）拿出来放到另外一个位置上，走出局部最小值陷阱。先引入以下概念：

• 任一格局之下，物体i的绝对痛苦指数DPi为此时自身具有的弹性势能；
• 任一格局之下，物体i的相对痛苦指数RDPi为此时其自身具有的弹性势能除以其半径的平方。

当h很小时，物体的移动速度很慢，我们就认为系统掉入了局部最小值陷阱。此时选出RDP最大的物体，随机放到五个位置（容器外围的四角和圆心）中去，如果这次的难民连续两次中标，就改为调整DP最小的物体。

算法的拟物部分较死板一些，而拟人部分有不少地方值得思考。包括如何判断局部最小值陷阱、如何选择难民以及如何重新放置难民，这几个点的设计直接会影响到算法的效率。关于如何放置难民，我测试了三种方案。

1. 随机取一个坐标，这是书上推荐的办法，但是我的测试却并不理想，解题的时间很看运气，往往很长时间找不到最优解。
2. 只从固定的点中选择，我选择的候选点是能包含容器的最小矩形的四个角，然后将难民分配到最远的角去，但是这种办法可能导致死循环，难民在两个点来回奔波，而且如果大量的空白区域被几个物体围在中心，会导致其它物体长时间不能进入。
3. 在四角的基础上，增加容器中心的位置，然后从五个位置随机选择一个点安放难民，试验表明这种方法是三种方法里面最好的。

可以看出，随机性的算法可能很长时间得不到解，而非随机性的算法就可能进入死循环，这是不是说明了“过于理性的思考可能使人发疯”呢。关于如何判断局部最小值陷阱，我觉得是比较难的，试了两种方法，理论上都有缺陷，但测试都还没发现问题。

1. 当h很小时，就判定进入了局部最小值陷阱，这是书上推荐的方法，实现很简单，但是限制了物体移动步长的最小值，可能恰好只需要一个更小的步长就能到达最优解；
2. 当h/U很小时，就判定进入局部最小值陷阱，这种方法可以判断出所有的局部最小值陷阱，同1一样，也存在将最优解误判为陷阱的可能性，但是我认为概率更小。
   

在实现的时候需要注意，h应该是和容器半径成比例的。用于测试的实验输入数据：

1. 容器的半径为6，小圆个数是7，它们的半径都是2；
2. 容器的半径为2.4143，小圆的个数是9，其中四个小圆的半径是1，五个小圆的半径是0.41415；
3. 容器的半径为2.4143，小圆的个数是17，其中四个小圆的半径是1，五个小圆的半径是0.41415，八个小圆的半径是0.2。

第一种输入最简单，可以用来检验程序是否存在大的逻辑问题，程序正确的话应很快就能计算出来。

第二种输入较复杂些，如果按照书上的算法，当落入局部最小值陷阱时，将难民随机放置到一个位置，运气不好可能较长时间才能解出来。所以我对其进行了优化，在圆的四角取四个点，加上圆心，将难民随机放到这五个位置之一，通常几秒钟就能找到答案。

第三种输入是最难的，在第二种输入的情况下增加了8个小圆，解题需要非常长的时间。大概跑了20分钟得到了结果。

程序是用MFC写的，遇到了在别的机器不能运行的问题，原因是工程选择了使用MFC的动态链接库，需要在vs里将工程的properties->General->Use of MFC改为使用静态链接库，同时properties->C/C++->Runtime Library也要改成MT。

CirclePacking小程序下载地址

CirclePacking源代码下载地址

参考文献

[1] 黄文奇，许如初。近世计算理论导引——NP难度问题的背景、前景以及求解算法研究。 科学出版社。