拟人拟物法求解不等圆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难度问题的背景、前景以及求解算法研究。 科学出版社。

©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页