梁利東，何東，朱良恒

（安徽工程大學(xué) 機(jī)械工程學(xué)院，安徽蕪湖 241000）

不等圓Packing問題是一類典型的NP-hard問題[1]，具有組合優(yōu)化問題的典型特征，廣泛應(yīng)用于紡織、造船、皮革工業(yè)、無線通信、航天器設(shè)計(jì)等領(lǐng)域[2-4]。作為圖形Packing問題一個(gè)分支，研究和探索更優(yōu)的排列模式和算法對(duì)于眾多行業(yè)都具有重要的科研和經(jīng)濟(jì)效益[5-6]。隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展，計(jì)算機(jī)數(shù)值算法尤其是啟發(fā)式算法來求解圖形Packing問題成為了這一領(lǐng)域研究的主要方向。

針對(duì)圓Packing問題，國(guó)內(nèi)外學(xué)者從各種角度提出了許多不同的求解策略。黃文奇等[7]將擬物算法與禁忌搜索策略相結(jié)合提出帶全局變換禁忌搜索算法GT-TS，通過禁忌規(guī)則和特赦規(guī)則的約束實(shí)現(xiàn)當(dāng)前局部最優(yōu)布局的全局變化策略。Hifi等[8]基于左下角優(yōu)先的策略求解矩形容器的不等圓Packing問題，提出了一種求解約束和無約束圓形Packing問題的啟發(fā)式算法。Birgin等[9]首次提出了使用鄰接矩陣來加速底層算法的概念，可實(shí)現(xiàn)將底層算法的計(jì)算復(fù)雜度從O(N2) 降 低到O(N)，大大縮短計(jì)算時(shí)間。王英聰?shù)萚10]利用群智能勞動(dòng)分工的任務(wù)分配來實(shí)現(xiàn)不等圓Packing問題的空間分配，從分配的角度分析不等圓Packing問題和群智能勞動(dòng)分工方法；之后又提出一種基于位置選擇的構(gòu)造法[11]，即將圓放置過程看成位置選擇過程，借鑒群智能勞動(dòng)分工中的刺激-響應(yīng)原理，將未布局空間的完整度和已布局空間的緊密度看作圓形物體選擇位置時(shí)的刺激和閾值，其自適應(yīng)閾值調(diào)整策略具有較好的可行性。Specht[12]設(shè)計(jì)并實(shí)現(xiàn)了一種基于圖論方法分析容器中未被填充區(qū)域的算法，并將未被填充的區(qū)域稱之為“孔”，通過探測(cè)“孔”的位置及其大小，結(jié)合跳轉(zhuǎn)策略，有效地提高了填充密度。余麗娟等[13]通過改進(jìn)基于Power圖的區(qū)域劃分，提出一種收斂速度更快的圓Packing算法，該算法是對(duì)局部Power圖算法（Local power diagram，LPD）的一種改進(jìn)。劉景法等[14]針對(duì)正三角容器內(nèi)的等圓Packing問題，提出了將啟發(fā)式格局更新策略與基于梯度法的局部搜索策略相融合的模擬退火算法，并與二分搜索結(jié)合，實(shí)現(xiàn)能量更低的最優(yōu)化格局，算例證明了其有效性。張維等[15]提出一種改進(jìn)遺傳模擬退火算法，通過計(jì)算生成一個(gè)合適大小的初始圓形容器來指導(dǎo)初始種群的生成，采用最優(yōu)保存策略來保證歷代的最優(yōu)解，并結(jié)合遺傳算法和模擬退火算法分階段來提升局部搜索和全局搜索性能。何琨等[16]提出了基于動(dòng)作空間的擬物求解算法，借鑒了求解矩形Packing問題中動(dòng)作空間的概念，通過化“圓”為“方”，將不規(guī)則的空閑空間近似表示為一系列矩形空間，有效地解決了跳出局部最優(yōu)的難點(diǎn)。

本文針對(duì)不同類型的擬物算法進(jìn)行驗(yàn)算和性能分析，以擬物算法為基礎(chǔ)，綜合分支搜索策略、多重退火策略提出多重退火分支搜索算法（Multiple annealing branch search，MABS）。針對(duì)連續(xù)優(yōu)化和組合優(yōu)化兩個(gè)層面描述不等圓Packing的最優(yōu)化策略和算法思想，并通過國(guó)際算例分析驗(yàn)證了不同形狀容器下的實(shí)驗(yàn)數(shù)據(jù)，證明該算法是有效可行的。

1 問題描述及算法模型

不等圓Packing問題是將一組已知個(gè)數(shù)的大小不等的圓不重疊放置在給定形狀的容器中，求容器的面積最小值或容器的填充密度最大值。以圓形容器為例，設(shè)容器半徑為r0，其中放置的第i個(gè)小圓的半徑為ri， 坐標(biāo)為 (xi,yi)。假定有N個(gè)圓需要放置，則對(duì)于某一布局，數(shù)學(xué)模型表示為：

擬物算法的求解思想是將需要放置的各圓看作為實(shí)心光滑彈性體，而容器看作是在填滿整個(gè)平面的無限彈性體中挖去相應(yīng)部分而形成的光滑空腔。通過模擬各圓在容器中受彈性力作用產(chǎn)生擠壓運(yùn)動(dòng)來實(shí)現(xiàn)連續(xù)優(yōu)化，其中各圓所受的彈性力合力方向也就是使目標(biāo)函數(shù)下降最快的負(fù)梯度方向，算法在迭代中通過調(diào)整容器大小以滿足布局中絕對(duì)彈性勢(shì)能最小化的約束要求，從而求的最優(yōu)解，算法的直觀幾何圖像如圖1所示。

本文的MABS算法將相對(duì)勢(shì)能作為約束條件并通過變鄰接系數(shù)加速方法構(gòu)建算法模型，以定步長(zhǎng)序列梯度下降算法為連續(xù)優(yōu)化方法，采用領(lǐng)域算子、變長(zhǎng)分支搜索以及多重模擬退火策略為組合優(yōu)化策略實(shí)現(xiàn)不等圓排列優(yōu)化布局，算法思想及總體結(jié)構(gòu)如圖2所示。

圖2 本文算法的總體結(jié)構(gòu)

相關(guān)描述和定義如下：

1）相對(duì)勢(shì)能

在大量實(shí)驗(yàn)過程中發(fā)現(xiàn)，對(duì)于不同尺寸大小的布局圖形，降低其絕對(duì)勢(shì)能的難易程度存在很大區(qū)別，主要體現(xiàn)在：較小的圓半徑和容器尺寸更容易得到較小的勢(shì)能，而對(duì)于擁有較大圓半徑和容器尺寸的算例，使勢(shì)能降至同一水平的難度明顯增大。其原因在于絕對(duì)勢(shì)能是建立在圖形間的絕對(duì)擠壓深度上，當(dāng)圖形整體擴(kuò)大至原來的k倍時(shí)，則所有擠壓深度也擴(kuò)大至原來的k倍，根據(jù)絕對(duì)勢(shì)能的定義，絕對(duì)勢(shì)能將擴(kuò)大至原來的k2倍。因此為了消除圖形絕對(duì)尺寸的影響，本文將容器的尺寸進(jìn)行歸一化預(yù)處理：將尺寸為L(zhǎng)的容器大小縮小到原來的1/L，其勢(shì)能相應(yīng)縮小到原來的1/L2。考慮到相對(duì)勢(shì)能的數(shù)值很小，不利于對(duì)比觀察，可將其擴(kuò)大106倍。以圓形容器為例，對(duì)于任一布局，相對(duì)勢(shì)能為約束的問題模型可表示為

圖3 相對(duì)勢(shì)能與絕對(duì)勢(shì)能對(duì)比

2）變鄰接系數(shù)加速方法

對(duì)于n個(gè)不等圓，擬物算法在每輪迭代都需計(jì)算n（n-1）/2個(gè)距離和每個(gè)圓所受的矢量力以及彈性勢(shì)能。然而從直觀的幾何圖像來看，某個(gè)圓所受的矢量力顯然只與相鄰圓有關(guān)，因此可通過構(gòu)建鄰接矩陣來加速搜索并降低計(jì)算復(fù)雜度，圖4為使用鄰接矩陣前后的所需計(jì)算距離表示。

圖4 使用鄰接矩陣前后對(duì)比

鄰接矩陣是表達(dá)圖形間接觸關(guān)系的矩陣，lij為各圓的鄰接關(guān)系屬性值，dij表示為各圓之間的距離，定義如下：

其中當(dāng)鄰接系數(shù)kn= 1時(shí)，鄰接矩陣為嚴(yán)格鄰接矩陣，客觀表達(dá)各圓的鄰接關(guān)系；當(dāng)kn＞ 1時(shí)，鄰接矩陣定義為過定鄰接矩陣，即不僅包含已有的鄰接關(guān)系，而且存在潛在鄰接關(guān)系，極有可能在迭代后期產(chǎn)生鄰接關(guān)系。為保證鄰接關(guān)系變化的相對(duì)穩(wěn)定性，嘗試在迭代的初期使用較大的kn以包含更多潛在的鄰接關(guān)系，同時(shí)可減小鄰接矩陣的計(jì)算間隔以應(yīng)對(duì)鄰接關(guān)系的頻繁變化。這樣，就可以依據(jù)迭代不同階段的特點(diǎn)通過調(diào)整kn與間隔代數(shù)以實(shí)現(xiàn)全計(jì)算周期的鄰接矩陣加速。實(shí)驗(yàn)也表明使用變鄰接系數(shù)鄰接矩陣可將算法的計(jì)算復(fù)雜度從O（N2）降低至O（N），大大縮短計(jì)算時(shí)間。

3）定步長(zhǎng)序列梯度下降算法

根據(jù)實(shí)現(xiàn)機(jī)制和策略的不同組合，可衍生出眾多擬物算法[17-18]。本文以定步長(zhǎng)序列梯度下降算法、定步長(zhǎng)同步BFGS算法和變步長(zhǎng)同步梯度下降算法為例進(jìn)行測(cè)試對(duì)比和分析。為了避免單次求解的隨機(jī)性干擾，所有數(shù)據(jù)都是基于隨機(jī)生成的初始布局運(yùn)行10次所取的平均值。圖5為3種算法的運(yùn)算時(shí)間和收斂勢(shì)能對(duì)比圖，結(jié)果表明：定步長(zhǎng)序列梯度下降擬物算法具有最高的運(yùn)算速度，其次是定步長(zhǎng)同步BFGS擬物算法、變步長(zhǎng)同步梯度下降算法，勢(shì)能收斂?jī)?yōu)劣情況則正好相反。

圖5 3種算法的性能比較

實(shí)驗(yàn)表明，盡管變步長(zhǎng)同步梯度下降擬物算法在對(duì)隨機(jī)初始布局的勢(shì)能收斂性方面表現(xiàn)最優(yōu)，但考慮到后期的組合優(yōu)化對(duì)布局?jǐn)_動(dòng)的穩(wěn)定性因素，本文又特別模擬實(shí)際運(yùn)算場(chǎng)景測(cè)試了變步長(zhǎng)同步梯度下降擬物算法與定步長(zhǎng)序列梯度下降算法對(duì)經(jīng)過擾動(dòng)的穩(wěn)定布局的勢(shì)能收斂性情況，結(jié)果如圖6所示?？梢姡ú介L(zhǎng)序列梯度下降算法在實(shí)際運(yùn)算場(chǎng)景的多數(shù)情況下勢(shì)能收斂性接近變步長(zhǎng)同步梯度下降擬物算法，而且又擁有最快的運(yùn)算速度，故采用定步長(zhǎng)序列梯度下降算法作為MABS算法的連續(xù)優(yōu)化模塊。

圖6 實(shí)際運(yùn)算場(chǎng)景下兩種擬物算法勢(shì)能收斂性對(duì)比

2 算法的組合優(yōu)化策略

2.1 鄰域算子

不等圓Packing的最優(yōu)布局在很大程度上取決于各圓的排列位置。多起點(diǎn)策略由于所生成的各個(gè)布局之間沒有任何聯(lián)系，在生成新布局時(shí)也無法利用之前生成的布局信息，搜索盲目而低效。因此運(yùn)用鄰域算子在不大幅改變布局的前提下，可通過對(duì)布局進(jìn)行小的變動(dòng)而得到更優(yōu)布局。本文采用3種鄰域算子：1）拋擲算子：隨機(jī)挑選一個(gè)圓并隨機(jī)“拋擲”到容器內(nèi)某處；2）交換算子：交換當(dāng)前布局中任意兩個(gè)不等圓的位置；3）振蕩算子：將所有圓在各自的隨機(jī)方向上做小幅度移動(dòng)。圖7為拋擲算子與交換算子生成新鄰域的過程示意。

圖7 利用交換算子與拋擲算子生成新鄰域的過程

相關(guān)定義如下：

定義2（鄰域算子）：對(duì)穩(wěn)定布局進(jìn)行一定變換，得到與之不同的穩(wěn)定布局，這一變換稱之為鄰域算子。

鄰域表征的是穩(wěn)定布局之間的關(guān)系，借助于鄰域算子，有可能對(duì)布局進(jìn)行累計(jì)的小幅度改進(jìn)以獲得最優(yōu)布局。此外，拋擲算子、交換算子與振蕩算子的多次疊加使用也可以看作是一個(gè)鄰域算子。顯然，疊加次數(shù)越多，所產(chǎn)生的布局就可能更不同于原有布局。換言之，可以認(rèn)為在一定范圍內(nèi)，3種算子的疊加次數(shù)與初末布局之間的相似度呈負(fù)相關(guān)。

2.2 變長(zhǎng)分支搜索策略

分支搜索策略源于優(yōu)勝劣汰思想，即從當(dāng)前最優(yōu)的父節(jié)點(diǎn)出發(fā)，使用鄰域算子生成一系列子節(jié)點(diǎn)，并比較子節(jié)點(diǎn)與父節(jié)點(diǎn)的勢(shì)能。若父節(jié)點(diǎn)勢(shì)能低于子節(jié)點(diǎn)勢(shì)能，則繼續(xù)從父節(jié)點(diǎn)生成新的子節(jié)點(diǎn)；否則將此子節(jié)點(diǎn)更新為父節(jié)點(diǎn)，并從新的父節(jié)點(diǎn)出發(fā)繼續(xù)搜索其更優(yōu)子節(jié)點(diǎn)。最優(yōu)布局的定義如下：

分支搜索策略總是從最優(yōu)節(jié)點(diǎn)生成新的節(jié)點(diǎn)，而生成的新節(jié)點(diǎn)只能經(jīng)過一次鄰域算子的處理，即分支長(zhǎng)度為1的分支生長(zhǎng)過程。由于分支較短易于陷入局部最優(yōu)布局，通過延長(zhǎng)分支長(zhǎng)度來擴(kuò)大搜索范圍，變長(zhǎng)分支搜索策略的算法流程如圖8所示。

圖8 變長(zhǎng)分支搜索策略流程圖

實(shí)驗(yàn)證明延長(zhǎng)分支長(zhǎng)度確實(shí)可以有效提高搜索質(zhì)量，但隨著分支長(zhǎng)度的增長(zhǎng)，計(jì)算耗費(fèi)的時(shí)間也同比增長(zhǎng)，因此需要在計(jì)算時(shí)間和求解質(zhì)量之間做出權(quán)衡。本文對(duì)分支長(zhǎng)度從1至3的分支搜索進(jìn)行了實(shí)驗(yàn)對(duì)比，實(shí)驗(yàn)結(jié)果表明分支長(zhǎng)度為2的分支搜索在耗時(shí)和求解質(zhì)量?jī)身?xiàng)指標(biāo)上均取得了較好的表現(xiàn)，如圖9所示。

圖9 不同長(zhǎng)度分支搜索性能對(duì)比

2.3 模擬退火接受準(zhǔn)則及多重退火策略

為進(jìn)一步提高布局多樣性，本文基于模擬退火算法，并以一定概率接受更差解的策略來達(dá)到全局最優(yōu)的目的。在退火過程中，假定某布局Xi，經(jīng)過鄰域算子擾動(dòng)生成新布局Xj，兩布局的勢(shì)能分別是fr(Xi)和fr(Xj)。根據(jù)Metropolis準(zhǔn)則，接受新布局Xj的概率表示為：

式中T為溫度，當(dāng)?shù)螖?shù)增加時(shí)取Tk+1=α·Tk，α=0．99。 當(dāng)fr(Xj)＞fr(Xi) 時(shí) ，有可知此時(shí)Pj∈(0,1)， 且溫度越低Pj越接近于0。這表明隨著迭代進(jìn)行，模擬退火接受更差解的概率逐漸趨于0，算法流程如圖10所示。

圖10 模擬退火算法流程圖

實(shí)現(xiàn)發(fā)現(xiàn)，200 ℃以上的初始溫度對(duì)解的質(zhì)量增長(zhǎng)不大，但迭代次數(shù)和運(yùn)算時(shí)間卻隨著初始溫度升高而大幅增加。因此考慮到求解質(zhì)量和運(yùn)算時(shí)間之間的平衡控制，本文依據(jù)實(shí)驗(yàn)數(shù)據(jù)采用400 ℃的退火初溫，并與貪婪算法進(jìn)行比較，實(shí)驗(yàn)數(shù)據(jù)如圖11所示。

圖11 不同初始溫度下模擬退火算法的求解質(zhì)量對(duì)比

由于模擬退火算法跳出局部最優(yōu)的能力主要體現(xiàn)在搜索前期，當(dāng)進(jìn)入勢(shì)能較為穩(wěn)定的搜索中后期，找到更優(yōu)解或跳出當(dāng)前局部最優(yōu)解的可能都很小。對(duì)此，提出多重退火的改進(jìn)思路，即當(dāng)判定模擬退火陷入局部最優(yōu)布局或缺乏找到較優(yōu)解的潛力時(shí)，及時(shí)跳出本次退火進(jìn)程，通過升高溫度開始新一輪退火過程。當(dāng)溫度重新升高，接受差解的概率將大大增加，搜索隨機(jī)性增強(qiáng)，相當(dāng)于對(duì)當(dāng)前布局進(jìn)行大幅重組，從而可跳到新的搜索域，以期望獲得可能的更優(yōu)布局。當(dāng)達(dá)到限定時(shí)間或限定重退火次數(shù)即可結(jié)束搜索，圖12為多重退火策略實(shí)施的勢(shì)能和溫度的變化曲線。

圖12 多重退火勢(shì)能溫度變化曲線

可見，在時(shí)間相同的情況下，通過增加多重退火判斷準(zhǔn)則（勢(shì)能下降量準(zhǔn)則），多輪模擬退火搜索在局部最優(yōu)布局上或缺少搜索前景的布局上極大降低了計(jì)算資源，同時(shí)實(shí)現(xiàn)了在不同狀態(tài)空間的尋優(yōu)過程，既擴(kuò)展了搜索廣度，也提高了最終求解的質(zhì)量。

3 實(shí)例分析

本文的MABS算法基于MATLAB軟件進(jìn)行編寫，用于實(shí)驗(yàn)的硬件環(huán)境為：筆記本電腦，CPU為四核Intel酷睿i5-7400，主頻3.0 GHz，運(yùn)行內(nèi)存16 GB。算法在6種不同形狀的容器中進(jìn)行不等圓排列布局優(yōu)化，取N=20，ri=（i=1,2,···,N），獲得的最優(yōu)布局如圖13所示，其中φ 為填充密度。

圖13 MABS算法在不同形狀容器中的最優(yōu)布局

其中，針對(duì)packomania官網(wǎng)中圓形容器與正方形容器的對(duì)應(yīng)算例，本文算法所得最優(yōu)解與其最佳記錄的排樣密度持平。其他4種形狀容器的不等圓Packing算例未見記錄，所得填充密度為本算法獨(dú)立得出。

表1 圓形容器中26個(gè)國(guó)際公開算例的實(shí)驗(yàn)結(jié)果

由此可知，MABS算法在近一半的算例上與之前最優(yōu)解持平，10個(gè)算例的解略差于記錄最優(yōu)解，改進(jìn)了4個(gè)算例的此前最優(yōu)解。可見MABS是具有較高性能的不等圓Packing算法。

本文算法基于最小包絡(luò)圓規(guī)則也可將不規(guī)則圖形Packing問題轉(zhuǎn)化為不等圓問題進(jìn)行求解。圖14表示為給定的4個(gè)不規(guī)則圖形及其形成的最小包絡(luò)圓，圖15為MABS算法所得的最優(yōu)布局。

圖14 不規(guī)則圖形及其最小包絡(luò)圓

圖15 不規(guī)則圖形Packing求解應(yīng)用

4 結(jié)束語(yǔ)

本文基于擬物算法模型提出了分支搜索策略、多重退火策略的高效不等圓Packing算法，在連續(xù)優(yōu)化和組合優(yōu)化兩個(gè)層次對(duì)算法進(jìn)行了優(yōu)化和改進(jìn)。以定步長(zhǎng)序列梯度下降擬物算法作為底層排列方法，并在運(yùn)算迭代的不同階段使用不同的鄰接系數(shù)及計(jì)算間隔計(jì)算鄰接矩陣，改進(jìn)了鄰接矩陣加速方法；采用分支搜索策略，通過增加分支長(zhǎng)度和領(lǐng)域算子來擴(kuò)展搜索空間以提升解的質(zhì)量；針對(duì)模擬退火算法迭代后期難以跳出局部最優(yōu)的缺陷提出了多重退火策略，并設(shè)置接受概率保證全局優(yōu)化過程。通過對(duì)國(guó)際公開算例集的實(shí)驗(yàn)結(jié)果以及在不規(guī)則圖形排列問題的求解應(yīng)用表明，該算法具有一定的有效性。