求解最小包容圓問題的一種有效算法

蔡 園， 蔣 毅， 王成露

（四川師范大學 數學科學學院 可視化計算與虛擬現實四川省重點實驗室，四川 成都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的有效性.

表4 2種算法的平均CPU時間Tab.4 Average running time of two methods

圖1給出了m=9時，由算法1產生的最小包容圓.可以清楚地看到，算法1求出的最優圓是非常精確的.

圖1 給定9個圓，由算法1產生的最小包容圓Fig.1 The smallest enclosing circle for nine circles by Algorithm 1