蔡 園， 蔣 毅， 王成露

（四川師范大學(xué) 數(shù)學(xué)科學(xué)學(xué)院 可視化計算與虛擬現(xiàn)實四川省重點實驗室，四川 成都610066）

1 背景介紹

最小包容圓問題是求解平面上包容所有給定圓的最小的圓，該問題簡寫為SEC［1］.這個問題最初是由Sylvester［2］在1857年提出的，并且可以在線性時間內(nèi)求解.最小包容圓問題有很多的應(yīng)用背景，比如規(guī)劃共享設(shè)施的位置、環(huán)境科學(xué)、模式識別、機械工程以及計算機制圖等，具體參考文獻［3-6］.許多文獻都涉及了最小包容圓問題的算法，例如切平面法［7］、二次規(guī)劃法［3］以及基于Voronoi圖的算法［8］.最小包容圓問題的數(shù)學(xué)模型描述為如下形式

其中，（x，y）和R分別代表最優(yōu)圓的圓心和半徑，ˉ=｛o1，o2，…，om｝表示m個給定的圓，它們的圓心和半徑分別為｛（a1，b1），（a2，b2），…，（am，bm）｝和｛r1，r2…，rm｝.根據(jù)半徑的不同，最小包容圓問題又可以分為以下的2種情況：ri=0以及ri＞0，其中i=1，2，…，m.對于ri=0，已經(jīng)有許多文獻報道，比如文獻［9-10］.本文主要研究ri＞0的情況，可參考文獻［3，7，11］等.在文獻［3］中，一共研究了4種著名的算法來解最小包容圓問題，包括二次規(guī)劃法、次梯度法，隨機增量法以及二階錐規(guī)劃法.在他們的計算實驗中，表明了最好的算法是二次規(guī)劃法.文獻［3］中二次規(guī)劃法最主要的思想就是把問題（1）轉(zhuǎn)化為如下的一個含有線性約束的二次規(guī)劃問題

（x0，y0）為任意給定的初始值，此時問題（2）是凸優(yōu)化［12］.在本文數(shù)值實驗中，與文獻［3］中最好的二次規(guī)劃法進行比較，新算法處理的數(shù)據(jù)越大越有效.值得提出的是，本文思想來源于文獻［13-14］中的方法.

2 算法描述

下面介紹本文提出的新算法.根據(jù)文獻［3］中的定理2.1，知道最小包容圓問題（1）與非凸二次規(guī)劃問題（2）在一些條件下是等價的.對于問題（2），是含有非凸約束的二次規(guī)劃問題.引入一個新變量V來替換平方項R2，并結(jié)合文獻［15-21］的有關(guān)線性松弛的研究，可以得到問題（2）的如下松弛問題

定理2.1問題（3）是問題（2）的線性松弛.

證明設(shè)（x，y，R，z）是問題（2）的一個可行解，則滿足如下不等式：

根據(jù)最小包容圓問題（2）中R和ˉR的定義，知道

由（5）和（7）式得，（x，y，R，z，V）（V=R2）是問題（3）的一個可行解.因此，問題（3）是問題（2）的線性松弛.

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

引理2.2當V=R2時，問題（3）與問題（2）是等價的.

問題（3）是含有線性約束的二次規(guī)劃問題，可以有效地求解.如果V=R2，則問題（3）的最優(yōu)解是問題（2）的最優(yōu)解；否則，根據(jù)文獻［13］中提到的切平面方法，在問題（3）中加入有效的切平面并求解新的二次規(guī)劃問題.基于這種思想，給出以下求解最小包容圓問題的算法.

算法1

步驟1 初始化：k=0.給定初始點（x0，y0），并計算

再用MATLAB求解問題（9），更新k：=k+1，得到問題（9）的最優(yōu)解（xk，yk，Rk，zk，Vk）.

引理2.3［3］設(shè)（x*，y*，R*）是問題（1）的一個最優(yōu)解當且僅當（x*，y*，0）是問題（2）的最優(yōu)解，且有R=R*時，問題（2）與問題（1）是等價的.

有關(guān)算法1中使用的割平面的研究，可參見文獻［13］.根據(jù)文獻［13］，可以得到（8）式是有效的切平面且算法1是收斂的.又由引理2.2和2.3，停機準則｜xk2+yk2-zk｜≤ε和｜Rk2-Vk｜≤ε是適定的.

定理2.4算法1產(chǎn)生的序列的極限點是問題（1）的最優(yōu)解.

證明設(shè)（x*，y*，R*，z*，V*）是算法1產(chǎn)生的序列的極限點，則有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）的最優(yōu)解.由引理2.2和2.3可知，（x*，y*，R*）是問題（1）的最優(yōu)解.

3 數(shù)值實驗

本節(jié)主要比較算法1和文獻［3］中二次規(guī)劃法在MATLAB中的數(shù)值表現(xiàn).所有的測試數(shù)據(jù)都是隨機產(chǎn)生的，每一組圓的產(chǎn)生均滿足獨立的正態(tài)分布N（0，16），每組圓的半徑均滿足均勻分布U（0，1），并且保證這些圓都是不相交的.在數(shù)值實驗中，取初始值（x0，y0）=（0，0），ε為10-6.測試了不同大小的隨機例子，m的范圍是從18 000到30 000.對于不同的m，分別取50組隨機產(chǎn)生的數(shù)據(jù)進行計算，再求這50組數(shù)據(jù)得到的平均值.產(chǎn)生的結(jié)果是在內(nèi)存4 GB、2.5 GHz英特爾酷睿處理器的個人電腦中得到.數(shù)值結(jié)果概括到表1～4中.在這些表中，m代表圓的個數(shù)，k代表算法1產(chǎn)生的迭代次數(shù)，nA是平均迭代次數(shù)，nmax是最大迭代次數(shù)，nmin是最小迭代次數(shù)，t是平均CPU時間（單位用s表示），QP［3］表示文獻［3］中的二次規(guī)劃法，Delta-T表示2種算法的時間差.最終的數(shù)值結(jié)果表明算法1處理的數(shù)據(jù)越大速度比二次規(guī)劃法越快.

表1 m＝30 000時，算法1每次迭代的計算結(jié)果Tab.1 Computational results for 30 000 circles by the Algorithm 1

表1給出了當m=30 000時，算法1求解最小包容圓問題的迭代次數(shù)及每步迭代所用的CPU時間.值得注意的是在k=3和k=4時，圓心已經(jīng)相等.

表2給出了算法1的迭代次數(shù).已知圓的個數(shù)m是從18 000到30 000之間取了5組值，針對同一組m分別取50組隨機數(shù)據(jù)所得的平均值結(jié)果.其中算法1的平均迭代次數(shù)約是4，最大迭代次數(shù)是6，最小迭代次數(shù)是3.不難發(fā)現(xiàn)隨著已知圓的個數(shù)m的逐漸增加，算法1的迭代次數(shù)幾乎沒有變化.因此，可以得出算法1的迭代次數(shù)與圓的個數(shù)之間沒有聯(lián)系.

表2 算法1的平均迭代次數(shù)Tab.2 The number of iterations on Algorithm 1

表3給出了QP［3］與算法1求出的最優(yōu)值.對于同一組m，可以看到2種算法求出的最優(yōu)值是相同的.

表3 2種算法的最優(yōu)目標值Tab.3 Objective function value of two methods

表4給出了QP［3］與算法1的平均CPU時間差.數(shù)據(jù)結(jié)果表明，當已知圓個數(shù)m取18 000到22 000時，算法1與QP［3］的計算時間幾乎相等.但隨著已知圓個數(shù)m（大于22 000）的逐漸增加，算法1的計算速度比QP［3］更快.當m=30000時，QP［3］的平均CPU時間約是算法1的1.6倍.因此，數(shù)值結(jié)果驗證了算法1的有效性.

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

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

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