求解最小包容圆问题的一种有效算法

2021-09-13王成露

四川师范大学学报（自然科学版） 2021年5期

蔡园，蒋毅，王成露

（四川师范大学数学科学学院可视化计算与虚拟现实四川省重点实验室，四川成都610066）

1 背景介绍

最小包容圆问题是求解平面上包容所有给定圆的最小的圆，该问题简写为SEC［1］.这个问题最初是由Sylvester［2］在1857年提出的，并且可以在线性时间内求解.最小包容圆问题有很多的应用背景，比如规划共享设施的位置、环境科学、模式识别、机械工程以及计算机制图等，具体参考文献［3-6］.许多文献都涉及了最小包容圆问题的算法，例如切平面法［7］、二次规划法［3］以及基于Voronoi图的算法［8］.最小包容圆问题的数学模型描述为如下形式

其中，（x，y）和R分别代表最优圆的圆心和半径，¯=｛o1，o2，…，om｝表示m个给定的圆，它们的圆心和半径分别为｛（a1，b1），（a2，b2），…，（am，bm）｝和｛r1，r2…，rm｝.根据半径的不同，最小包容圆问题又可以分为以下的2种情况：ri=0以及ri＞0，其中i=1，2，…，m.对于ri=0，已经有许多文献报道，比如文献［9-10］.本文主要研究ri＞0的情况，可参考文献［3，7，11］等.在文献［3］中，一共研究了4种著名的算法来解最小包容圆问题，包括二次规划法、次梯度法，随机增量法以及二阶锥规划法.在他们的计算实验中，表明了最好的算法是二次规划法.文献［3］中二次规划法最主要的思想就是把问题（1）转化为如下的一个含有线性约束的二次规划问题

（x0，y0）为任意给定的初始值，此时问题（2）是凸优化［12］.在本文数值实验中，与文献［3］中最好的二次规划法进行比较，新算法处理的数据越大越有效.值得提出的是，本文思想来源于文献［13-14］中的方法.

2 算法描述

下面介绍本文提出的新算法.根据文献［3］中的定理2.1，知道最小包容圆问题（1）与非凸二次规划问题（2）在一些条件下是等价的.对于问题（2），是含有非凸约束的二次规划问题.引入一个新变量V来替换平方项R2，并结合文献［15-21］的有关线性松弛的研究，可以得到问题（2）的如下松弛问题

定理2.1问题（3）是问题（2）的线性松弛.

证明设（x，y，R，z）是问题（2）的一个可行解，则满足如下不等式：

根据最小包容圆问题（2）中R和¯R的定义，知道

由（5）和（7）式得，（x，y，R，z，V）（V=R2）是问题（3）的一个可行解.因此，问题（3）是问题（2）的线性松弛.

由定理2.1的证明，可以得到以下引理.

引理2.2当V=R2时，问题（3）与问题（2）是等价的.

问题（3）是含有线性约束的二次规划问题，可以有效地求解.如果V=R2，则问题（3）的最优解是问题（2）的最优解；否则，根据文献［13］中提到的切平面方法，在问题（3）中加入有效的切平面并求解新的二次规划问题.基于这种思想，给出以下求解最小包容圆问题的算法.

算法1

步骤1 初始化：k=0.给定初始点（x0，y0），并计算

再用MATLAB求解问题（9），更新k：=k+1，得到问题（9）的最优解（xk，yk，Rk，zk，Vk）.

引理2.3［3］设（x*，y*，R*）是问题（1）的一个最优解当且仅当（x*，y*，0）是问题（2）的最优解，且有R=R*时，问题（2）与问题（1）是等价的.

有关算法1中使用的割平面的研究，可参见文献［13］.根据文献［13］，可以得到（8）式是有效的切平面且算法1是收敛的.又由引理2.2和2.3，停机准则｜xk2+yk2-zk｜≤ε和｜Rk2-Vk｜≤ε是适定的.

定理2.4算法1产生的序列的极限点是问题（1）的最优解.

证明设（x*，y*，R*，z*，V*）是算法1产生的序列的极限点，则有x*2+y*2-z*=0且V*=R*2.由V*=R*2，可以得到V*-2αR*+α2≥0是成立的.因此，问题（3）与问题（2）是等价的，即（x*，y*，R*，z*，V*）是问题（2）的最优解.由引理2.2和2.3可知，（x*，y*，R*）是问题（1）的最优解.

3 数值实验

本节主要比较算法1和文献［3］中二次规划法在MATLAB中的数值表现.所有的测试数据都是随机产生的，每一组圆的产生均满足独立的正态分布N（0，16），每组圆的半径均满足均匀分布U（0，1），并且保证这些圆都是不相交的.在数值实验中，取初始值（x0，y0）=（0，0），ε为10-6.测试了不同大小的随机例子，m的范围是从18 000到30 000.对于不同的m，分别取50组随机产生的数据进行计算，再求这50组数据得到的平均值.产生的结果是在内存4 GB、2.5 GHz英特尔酷睿处理器的个人电脑中得到.数值结果概括到表1～4中.在这些表中，m代表圆的个数，k代表算法1产生的迭代次数，nA是平均迭代次数，nmax是最大迭代次数，nmin是最小迭代次数，t是平均CPU时间（单位用s表示），QP［3］表示文献［3］中的二次规划法，Delta-T表示2种算法的时间差.最终的数值结果表明算法1处理的数据越大速度比二次规划法越快.

表1 m＝30 000时，算法1每次迭代的计算结果Tab.1 Computational results for 30 000 circles by the Algorithm 1

表1给出了当m=30 000时，算法1求解最小包容圆问题的迭代次数及每步迭代所用的CPU时间.值得注意的是在k=3和k=4时，圆心已经相等.

表2给出了算法1的迭代次数.已知圆的个数m是从18 000到30 000之间取了5组值，针对同一组m分别取50组随机数据所得的平均值结果.其中算法1的平均迭代次数约是4，最大迭代次数是6，最小迭代次数是3.不难发现随着已知圆的个数m的逐渐增加，算法1的迭代次数几乎没有变化.因此，可以得出算法1的迭代次数与圆的个数之间没有联系.

表2 算法1的平均迭代次数Tab.2 The number of iterations on Algorithm 1

表3给出了QP［3］与算法1求出的最优值.对于同一组m，可以看到2种算法求出的最优值是相同的.

表3 2种算法的最优目标值Tab.3 Objective function value of two methods

表4给出了QP［3］与算法1的平均CPU时间差.数据结果表明，当已知圆个数m取18 000到22 000时，算法1与QP［3］的计算时间几乎相等.但随着已知圆个数m（大于22 000）的逐渐增加，算法1的计算速度比QP［3］更快.当m=30000时，QP［3］的平均CPU时间约是算法1的1.6倍.因此，数值结果验证了算法1的有效性.