梁利東，何東，朱良恒

（安徽工程大學 機械工程學院，安徽蕪湖 241000）

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

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

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

1 問題描述及算法模型

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

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

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

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

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

1）相對勢能

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

圖3 相對勢能與絕對勢能對比

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

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

圖4 使用鄰接矩陣前后對比

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

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

3）定步長序列梯度下降算法

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

圖5 3種算法的性能比較

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

圖6 實際運算場景下兩種擬物算法勢能收斂性對比

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

2.1 鄰域算子

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

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

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

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

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

2.2 變長分支搜索策略

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

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

圖8 變長分支搜索策略流程圖

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

圖9 不同長度分支搜索性能對比

2.3 模擬退火接受準則及多重退火策略

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

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

圖10 模擬退火算法流程圖

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

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

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

圖12 多重退火勢能溫度變化曲線

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

3 實例分析

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

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

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

表1 圓形容器中26個國際公開算例的實驗結(jié)果

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

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

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

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

4 結(jié)束語

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