1引言
将个不重叠的物体（例如圆形、球体）放入一个规则的凸容器（例如在二维或更高维的欧氏空间中的正方形、正六面体），使得包装密度最大化的问题是一类数量繁多且已被广泛研究的问题。对于二维欧式空间，对应的一类问题称为圆包装问题(circle packing problems)。其中最具代表性的问题是将相同大小的圆放入一个正方形中[packing equal circles in a square(PECS)]的问题。相应地，对于三维欧式空间，对应的一类问题称为球体包装(sphere packing)问题，其中研究最多的则是将相同大小的球体放入一个正六面体中[packing equal spheres ina cube(PESC)]。
2问题描述
PECS：给定个单位圆，我们需要将它们不重叠地放入一个正方形容器中，同时使得容器的边长最小。
其数学模型如下：
其中，正方形的中心为二维笛卡尔坐标系的原点；和分别表示单位圆和单位圆的圆心坐标。约束(2)称为不重叠约束，保证两个单位圆之间不发生重叠（圆心之间距离不小于2）；约束(3)(4)称为容纳约束(containment constraints)，保证每一个单位圆都完全容纳于正方形容器内（圆心与边界距离不小于1）。
将此模型拓展到三维欧氏空间的PESC问题中：
其中，正六面体的中心为三维笛卡尔坐标系的原点；和分别表示单位球体和单位球体的球心坐标。约束(6)保证两个单位球体之间不发生重叠，约束(7)~(9)保证每一个单位球体都完全容纳于正六面体容器内。
3PECS和PESC的子问题建模中心思想
由于PECS和PESC的目标函数是关于边长单调的，因此我们可以不断修正的大小，并通过使用惩罚函数的方式，将原来带约束的优化问题转化为将逐渐变小的边长作为参数的无约束优化子问题，进而通过无约束优化方法来解决这些子问题。
固定边长的PECS(PECS-L)
前面说到的无约束子问题(用PECS-L表示)目标是尝试在固定边长的正方形容器中找到一个将个单位圆不重叠的放置方法。我们用重叠部分的深度的平方和作为目标函数（惩罚函数）。
给定正方形容器的固定边长，以及候选的放置方式，表示为个圆的坐标向量，PECS-L的目标函数定义为：
其中：
其中，表示单位圆和之间重叠的深度，表示单位圆和垂直边界之间重叠的深度，表示单位圆和水平边界之间重叠的深度，如下图所示：
目标函数刻画了放置方式在边长为的正方形容器中违背原问题约束的程度。也就是说，当且仅当时，为可行解。
固定边长的PESC(PESC-L)
PESC是PECS在三位欧氏空间上的自然拓展，相应地，其无约束子问题PESC-L定义如下：
给定正六面体容器的固定边长，放置方式，表示单位球的球心坐标，目标函数定义为：
其中表示单位球体和之间重叠的深度，表示单位球体和平面或之间重叠的深度，表示单位球体和平面或之间重叠的深度，表示单位球体和平面或之间重叠的深度。同样地，当且仅当时，为可行解。
4PBTS算法
基于扰动的阈值搜索算法[ perturbation-based thresholding search algorithm(PBTS)]是一个用于解决PTCS和PTSC的随机优化算法。
算法主框架
伪代码如下：
其中，和表示当前搜索到的最优可行放置方法及其对应边长，表示当前解。
我们可以看到，整个搜索过程分为两个阶段：第一阶段旨在通过逐步缩小边长以快速找到高密度的放置方法；第二阶段则通过多起点方法(multistart approach)来寻找全局最优解。为了保证搜索的高效性，算法使用了容器调整法(container adjustment method)和结合了混合扰动策略与局部最优的阈值搜索法(a thresholding search procedure integrating a mixed perturbation strategy and local optimization)。
在第一阶段中，我们从经验估计得到的包装密度（圆（球）的总面（体）积占容器总面（体）积的比例，即对PECS，）开始，算出初始容器大小，并随机生成初始放置方法（一般不可行）。接着，我们通过阈值搜索法不断降低目标函数的值，从而找到一个当前边长下的可行解。然后，我们在保证可行的前提下，逐渐提升密度（相应地，缩短边长），继续通过阈值搜索法尝试找到可行解，直到无法找到可行解为止（最后得到的解仍不可行）。最后，通过容器调整方法对前面得到的解进行局部优化调整，得到的可行解保存在和中。
在第二阶段中，每次迭代开始时，首先根据随机生成一个放置方式，接着通过阈值搜索法尝试找到可行解。如果找到了一个可行解，则在保证可行的前提下，通过容器调整法最小化边长。如果得到的边长优于，则更新和。为了提升算法的鲁棒性，如果当次迭代没能找到可行解或改进最优边长，则在下一次迭代时换一种扰动方式（伪代码中以Flag表示）。当到达时间限制时迭代停止。
下文中，我们将分别介绍算法主框架中的阈值搜索法、容器调整法；在此之前，我们先介绍阈值搜索法中用到的局部优化方法和两个扰动算子：
局部优化
给定放置方法，为了使其达到关于目标函数和容器调整法中目标函数的局部最优解，阈值搜索法和容器调整法都用到了一个流行且高效的拟牛顿法——限制变尺度法(L-BFGS)。L-BFGS使用线性搜索方法，在多数算例中可以使梯度的最大范数达到高达的精度。
关于什么是L-BFGS：
https://cloud.tencent.com/developer/article/1573935
关于什么是范数：
https://blog.csdn.net/qq_24819773/article/details/85718335
扰动算子
为了确保有效性和泛用性，阈值搜索方法使用了两种互补且通用的扰动方式以避免陷入局部最优：均匀随机扰动[uniformly random perturbation(URP)]和连续随机扰动[ sequential random perturbation (SRP)]。如前文中提到的那样，优化过程是自适应地使用这两种扰动方式，从而提高鲁棒性的。
均匀随机扰动的伪代码如下：
其中，为含个变量的候选解，为扰动强度，默认值为0.8（根据Grosso et al. (2010)）。URP算子对每个变量做出（均匀分布）的改变。
均匀随机扰动的伪代码如下：
SRP算子旨在通过连续地小幅扰动和局部优化从而跳过低质量的局部最优，使得优化过程可以克服目标函数中的小障碍，但不改变全局最优解，如下图所示：
其中，黑色实线是目标函数的曲线，红色虚线是通过SRP过程后呈现出的预期效果。可以看到，SRP过程通过跳出目标函数较小的下凸部分，使得目标函数变得平滑了。通过这个效果，SRP算子帮助算法找到更小并且更有希望产生全局最优的样本集合（平滑后函数的极小点），同时在许多情况下不改变全局最优解的同时帮助目标函数达到全局最优。
值得注意的是，SRP具有较强的泛用性，对于任意一阶导函数非凸的函数都可以使用，而且有较好的效果。
介绍完局部优化方法和两个扰动算子后，我们接下来就认识一下阈值搜索法：
阈值搜索法
为了搜索PECS-L和PESC-L的全局最优解，我们需要用到阈值搜索法，其包含了一个局部搜索过程以及上节提到的两种扰动算子，伪代码如下：
可以看到，进入函数的初始解首先进行一次局部优化，并作为函数中找到的当前最优解；接下来进入循环当中继续对其进行优化，当迭代超过一定次数（以表示，代表搜索深度）没有改进或找到可行解时（需要注意，解决的问题存在固定的边长），迭代停止。
而在每次的迭代过程当中，当前解首先经过一次URP或SRP扰动，由决定。上节中提到，URP的扰动强度值恒为0.8，而为了保证强度与分化程度的平衡，我们使用下面的周期函数来计算SRP的扰动强度：
其中，为当前迭代次数；为周期，在这里设置为72；和分别表示的最小值和最大值（由实验得出）。随的周期性变化如下图所示：
接着，我们对扰动后的解进行局部优化，然后将其目标函数值与优化前的目标函数值进行比较，记其差值为。若小于阈值且（一般意味着没有改进或产生了对称的无效改进），则接受成为新的当前解。和分别记录了接受与拒绝更新当前解的次数。若，即接受次数过多时，需要将阈值降低至，否则提高至。其中代表预期接受率，系数。
值得注意的是，与SRP类似，阈值搜索法同样具有较强的泛用性，可以用于其它非凸可微函数的全局优化当中。
容器调整法
容器调整法通过对边长和单位圆（球）心坐标的小幅调整，使解变得可行，并且局部优化边长。
这里用到的是惩罚函数法(SUMT)，又称外点法、顺序无约束极小化技术。这种方法将前文中PECS和PESC的带约束优化问题模型转化为一系列的无约束问题，进而按顺序解决它们，直到达到局部可行最小值。转换的无约束问题表示如下：
其中，为惩罚因子，它的每个值对应一个无约束优化函数；的定义与基本相同（对PECS忽略项）：
关于惩罚函数法：
https://zhuanlan.zhihu.com/p/376587367
其伪代码如下：
其中，惩罚因子初始值为，并在次局部优化中每次增大5倍以增强惩罚力度，使得最终在的情况下获得足以看作可行的解。
以上即为PBTS的全部内容，你看明白了吗？
-The End-
⬇️
https://github.com/INFORMSJoC/2022.0004
参考文献：
原论文：
[1]Xiangjing Lai, Jin-Kao Hao, Renbin Xiao, Fred Glover (2023) Perturbation-Based Thresholding Search for Packing Equal Circles and Spheres. INFORMS Journal on Computing
文案&排版：吕其乐（华中科技大学管理学院本科二年级）
指导老师：秦虎老师 （华中科技大学管理学院）
如对文中内容有疑问，欢迎交流。PS：部分资料来自网络。
如有需求，可以联系：
秦虎老师（华中科技大学管理学院，微信号tigerqin_china）
吕其乐（华中科技大学管理学院本科二年级：[email protected]）
欢迎大家加入数据魔术师粉丝群，我们的活动将会通过粉丝群优先发布， 学习资料将通过粉丝群分享。
转载自tigerqin_china（微信公众号，久章智能于 2023-09-11 转载），内容保持原貌，版权归原作者及原账号所有。阅读原文 ↗