最小比率旅行商問題的陰陽平衡優(yōu)化算法

許秋艷，馬 良，劉 勇

(1. 上海理工大學管理學院，上海 200093；2. 鹽城工學院信息工程學院，江蘇 鹽城224051)

1 引言

最小比率旅行商問題(minimum ratio traveling salesman problem，MRTSP)是經(jīng)典旅行商問題的一種擴展形式：除考慮行程距離外，推銷員從一個城市到另一個城市將獲得一定的收益，該問題的目標就是選擇最佳旅行線路，最小化總路程和總收益的比值。該模型同時考慮距離和收益，和只考慮距離最小相比，該目標函數(shù)等價于最小化費用效益比，更符合實際情況。

旅行商問題屬于NP難問題，在其基礎(chǔ)上構(gòu)建的最小比率旅行商問題MRTSP也屬于NP難問題[1]。目前，用于求解MRTSP的算法主要有模擬退火算法、競爭決策算法、大洪水算法、蝙蝠算法、蟻群優(yōu)化算法、引力搜索算法和中心引力優(yōu)化算法等[2，3]。這些方法為求解MRTSP提供了思路，但是計算效果需要進一步優(yōu)化。為有效處理MRTSP問題，本文將基于陰陽平衡優(yōu)化(Yin-Yang-pair optimization，YYPO)算法[4]設(shè)計求解方法。

陰陽平衡優(yōu)化算法是在2016年由Punnathanam等提出的，其優(yōu)化策略源于中國傳統(tǒng)文化——陰陽學說。在智能計算領(lǐng)域，算法的局部開發(fā)與全局探索能力既對立又互補，如何實現(xiàn)兩者的平衡是算法設(shè)計的關(guān)鍵。YYPO算法將局部開發(fā)與全局探索分別視為陰和陽兩個方面，期望通過陰陽平衡實現(xiàn)局部開發(fā)與全局探索的均衡[5]。

因獨特的優(yōu)化機制，陰陽平衡優(yōu)化算法備受關(guān)注。目前，該算法已經(jīng)成功用于發(fā)動機系統(tǒng)設(shè)計[6]、光伏逆變器PID調(diào)節(jié)[7]、風力渦輪機設(shè)計[8]和連續(xù)型選址[9]等。總體來看，陰陽平衡優(yōu)化算法的研究處于起步階段，主要用于求解連續(xù)優(yōu)化問題，還沒有涉及到組合優(yōu)化問題。本文將陰陽平衡優(yōu)化算法用于求解最小比率旅行商問題，擴大算法的應(yīng)用范圍，也為該問題的求解提供可行有效的方法。為此，本文采用佳點集法產(chǎn)生初始群體，提高初始解質(zhì)量；基于超球體和歸檔集對解不斷進行更新；并采用相對位置索引法直接產(chǎn)生可行解，避免多個約束條件的處理；基于綜卦變換引入反轉(zhuǎn)操作，提高算法的局部搜索性能。數(shù)值實驗驗證了算法求解最小比率旅行商問題的可行性和有效性。

2 最小比率旅行商問題MRTSP的數(shù)學模型

最小比率旅行商問題的具體模型可描述如下[1-3]

本文考慮的網(wǎng)絡(luò)屬于完全圖G=(V，E) (否則，可通過計算任兩個點間最短路轉(zhuǎn)換成等價的完全圖)。令V={1，2，…，n} 代表頂點集，表示城市的位置；E={(i，j)：i≠j，i，j∈V}代表邊的集合；D=(dij)n×n代表距離矩陣，表示任兩個城市間的距離，其中?i，j∈V，dij=dji，dii=+∞，dij>0；P=(pij)n×n代表收入矩陣，表示旅行商(推銷員)完成從一個城市到另一個城市的銷售可得到的收入，其中?i，j∈V，pij=pji，pii=0，pij>0。設(shè)那么最小比率旅行商問題的模型為

(1)

其中，式(1)表示目標函數(shù)，在經(jīng)過所有點的回路中求總行程與總收入之比最?。皇?2)和(3)要求在一次行程中每個城市只訪問一次；式(4)要求沒有任何子回路生成；|S|代表子圖S中包含的圖G中頂點個數(shù)；式(5)表示決策變量取值范圍。同時滿足上述約束條件的解就形成了一條Hamilton 回路。

3 最小比率旅行商問題的陰陽平衡優(yōu)化算法

受到太極圖的啟示，陰陽平衡優(yōu)化算法利用半徑為1的超球體表示其搜索空間。算法首先在超球體內(nèi)隨機生成兩個點，并比較這兩個點的適應(yīng)度值。令P1表示值較優(yōu)的點，而令P2表示值較差的點。然后分別以P1和P2為中心，δ1和δ2為步長在超球體中進行優(yōu)化搜索。在尋優(yōu)過程中，δ1會逐步減小而δ2會逐步變大。因此，點P1的搜索區(qū)域?qū)⒅饾u收縮而進行小范圍尋優(yōu)，在算法中被稱為局部開發(fā)，有收縮的含義，和陰對應(yīng)；點P2的搜索區(qū)域?qū)⒅饾u擴展而進行大范圍尋優(yōu)，在算法中被稱為全局探索，有擴張的含義，和陽對應(yīng)。利用P1和P2進行的局部開發(fā)和全局探索缺一不可，有陰陽互根的含義。點P1的搜索范圍越來越小，而點P2的搜索范圍越來越大，有陰陽對立的含義。在算法尋優(yōu)中，不斷增強基于點P1的局部搜索和點P2的全局搜索能力，有陰陽皆長的含義。算法期望通過基于點P1和點P2的兩種搜索以實現(xiàn)局部開發(fā)和全局探索的權(quán)衡，有陰陽平衡的含義。以下給出陰陽平衡優(yōu)化算法求解最小比率旅行商問題的具體實現(xiàn)過程。

3.1 基于佳點集的初始化

考慮到最小比率旅行商問題NP-難的特征，設(shè)置合理的初始解，可改善算法的搜索效率。此外，雖然智能優(yōu)化算法不依賴于初始解，但是提高初始解質(zhì)量能夠加快優(yōu)化速度。而陰陽平衡優(yōu)化算法僅通過P1和P2兩個點進行優(yōu)化搜索，提高P1和P2對應(yīng)初始點質(zhì)量，可提高算法優(yōu)化效率。本文不再采用基本陰陽平衡優(yōu)化算法對P1和P2隨機初始化的方法，而采用一種性能較好的初始化方法——佳點集法。在佳點集中選擇最好的兩個解賦值給P1和P2。

佳點集的概念由華羅庚等首先提出，具體方法如下[10]：

設(shè)Gn為n維空間的單位立方體，r=(r1，r2，…，rn) 為Gn內(nèi)的點。令Pm(k)是Gn內(nèi)的包含m個點的集合，如果形為Pm(k)={({r1k}，{r2k}，…，{rsk}，k=1，2，…，m}的偏差φ(m)滿足：φ(m)=C(r，ε)m-1+ε，那么就稱r為佳點，Pm(k)為佳點集。其中，ε代表無窮?。籆(r，ε)代表僅和r與ε有關(guān)的常量。

在具體應(yīng)用時，可以設(shè)置rk={2cos 2πk/q}，1≤k≤n，q為符合(q-3)/2≥n的最小素數(shù)，那么r即為佳點。此外，也可以設(shè)置rk={ek}，1≤k≤n那么r也是佳點。其中，{t}代表t的小數(shù)部分。研究表明，相對隨機點集而言，佳點集分布更加均勻[2，11]。

3.2 基于超球體的解更新

算法分別采用P1和P2作中心點，δ1和δ2作步長，并在超球體中進行優(yōu)化搜索。因為P1和P2兩點采用相同的更新方法，為方便起見，用P表示P1或P2，并用δ表示δ1或δ2。首先將P進行復(fù)制，生成2D(D表示變量維數(shù))個同樣的點，并分別用NP1，NP2，…，NP2D進行表示；然后對這些點進行更新。在算法中，可以對點的一個分量或者所有分量進行更新。對點的一個分量進行更新時，具體方法如下：

(6)

(7)

在對點的所有分量進行更新時，具體方法如下：

(8)

(9)

在算法中，通過等概率來選擇上述兩種方法的一種進行點的更新。此外，在利用點P1生成的2D個新解中，挑選適應(yīng)度值最優(yōu)的點，直接替換P1。對P2生成的新解也進行同樣的操作。

3.3 基于歸檔集的解更新

在計算過程中，算法采取精英保留策略，在進行基于超球體的解更新前，會將P1與P2兩個點存儲到歸檔集內(nèi)。每隔I代將利用歸檔集中的2I個點對當前的P1與P2兩點做更新操作。其中，I表示歸檔集的更新頻率，并且I是在Imin和Imax之間隨機產(chǎn)生的整數(shù)。

該階段的解更新方法如下：第一步，從歸檔集內(nèi)挑選適應(yīng)度值最優(yōu)的點并和P1進行比較，如果該點優(yōu)于P1，則這兩個點進行交換；第二步，再次從歸檔集內(nèi)中挑選適應(yīng)度值最優(yōu)的點并和P2進行比較，如果該點優(yōu)于P2，則這兩個點也進行交換；第三步，清空歸檔集，并生成的新的歸檔集更新頻率參數(shù)I。

此外，算法在該階段還對兩個搜索步長δ1和δ2進行更新，具體方法如下：

δ1=δ1-(δ1/α)

(10)

δ2=δ2+(δ2/α)

(11)

其中，α表示收縮/擴張因子；δ1和δ2分別表示是點P1和P2的搜索步長。

3.4 候選解的生成

為充分發(fā)揮陰陽平衡優(yōu)化算法求解連續(xù)優(yōu)化問題的優(yōu)勢，仍然采用連續(xù)空間表示算法的搜索空間。但最小比率旅行商問題屬于組合優(yōu)化問題，可行解表示城市的訪問順序。如何實現(xiàn)算法搜索個體的位置到問題解的轉(zhuǎn)換是必須要解決的問題?？紤]到最小比率旅行商問題的解是離散的特征，本文采用基于相對位置索引法的解生成方法[12]。

在相對位置索引法中，首先將最小實數(shù)轉(zhuǎn)換為最小整數(shù)1，然后將第二小的實數(shù)轉(zhuǎn)換為下一個高位整數(shù)2，以此類推，可以對所有實數(shù)進行編號。最后利用編號構(gòu)成城市的訪問順序。該方法和隨機鍵編碼方法類似，以下通過一個算例作進一步說明。例如，一個有四個城市的最小比率旅行商問題，陰陽平衡優(yōu)化算法的一個搜索個體位置為(0.52，0.68，0.15，0.97)，按照上述方法，得到的城市的訪問順序為(2，3，1，4)。

此外，如果在轉(zhuǎn)換前向量中有多個相同的分量，這種方法可能會產(chǎn)生不可行解。本文設(shè)計一種修復(fù)機制，即根據(jù)向量位置的先后順序進行編碼。以下以相對位置索引法舉例說明，轉(zhuǎn)換前的向量(0.02，0.58，0.78，0.98，0.78)，第三個分量和第五個分量取值相同，根據(jù)上述修復(fù)機制，其編號分別為3和4，轉(zhuǎn)換后的向量為(1，2，3，5，4)。

因此，采用上述候選解的生成方法，可以直接產(chǎn)生可行解，每個城市只被訪問一次，并且在任何一個子集中不會構(gòu)成完整的Hamilton 回路。在算法搜索過程中，可以直接計算目標函數(shù)值，不再考慮模型中四個約束條件。

3.5 基于綜卦變換的局部搜索

基本陰陽平衡優(yōu)化算法的全局搜索能力較強，而局部搜索能力較弱[13，14]。為求解最小比率旅行商問題，需要提高算法的局部尋優(yōu)能力。和陰陽學說一樣，《易經(jīng)》也屬于中國傳統(tǒng)文化。從算法設(shè)計角度來看，卦的變化蘊含了優(yōu)化思想。例如，綜卦是把本卦倒過來看，例如圖1所示，姤卦的綜卦是夬卦。綜卦的哲理是從另一個角度來看同一問題[15，16]。

圖1 綜卦變化示意圖

受此啟發(fā)，將當前解視為本卦，對其進行綜卦變化，具體將采用反轉(zhuǎn)操作。在反轉(zhuǎn)操作中，旋轉(zhuǎn)變化的城市位置不再是起點和終點，而是隨機選擇兩個城市位置作為起點和終點，然后將隨機選出的兩個城市之間的訪問順序顛倒過來。圖2給出一個示例，其中T表示反轉(zhuǎn)操作前的解，T′ 表示反轉(zhuǎn)操作后的解。

圖2 反轉(zhuǎn)操作示意圖

3.6 算法主要計算步驟

綜上所述，求解最小比率旅行商問題的陰陽平衡優(yōu)化算法具體步驟如下：

步驟1：根據(jù)佳點集法產(chǎn)生初始點，并選擇最好的兩個點分別賦值給P1和P2

步驟2：在基于超球體的解更新階段，按等概率原則分別利用P1與P2進行優(yōu)化搜索；

步驟3：每間隔I次迭代，則進入基于歸檔集的解更新階段，否則，轉(zhuǎn)步驟4；

步驟4：若迭代次數(shù)達到預(yù)設(shè)最大值，轉(zhuǎn)步驟6；否則，轉(zhuǎn)步驟5。

步驟5：若點P2適應(yīng)度值優(yōu)于點P1適應(yīng)度值，則交換這兩個點的取值與搜索步長，轉(zhuǎn)步驟2；

步驟6：對當前最好解進行基于綜卦的反轉(zhuǎn)操作，算法停止，輸出最好結(jié)果。

4 數(shù)值實驗

實驗共分為三個部分：首先采用并行程序設(shè)計陰陽平衡優(yōu)化算法；然后比較陰陽平衡優(yōu)化算法和微粒群優(yōu)化算法、引力搜索算法、生物地理學優(yōu)化算法以及最有價值球員算法的求解性能；最后對距離矩陣和收益矩陣的元素進行參數(shù)靈敏度分析。在本文中，采用最小比率旅行商問題的3個典型問題進行仿真。這些算例定義如下[1-3]：

算例1： 設(shè)給定對稱賦權(quán)完全圖G={V，E}，距離矩陣D和收益矩陣P定義如下

算例2： 設(shè)給定對稱賦權(quán)完全圖G={V，E}，距離矩陣D和收益矩陣P定義如下

算例3： 設(shè)給定對稱賦權(quán)完全圖G={V，E}，距離矩陣D和收益矩陣P定義如下：

4.1 算法并行程序設(shè)計

迄今為止，在智能優(yōu)化算法領(lǐng)域，使用并行程序設(shè)計算法是一個重要方向，現(xiàn)有工作主要是將計算任務(wù)分配給多臺計算機完成，而在一臺計算機上實現(xiàn)并行智能優(yōu)化算法的研究工作極其罕見。當前計算機至少配置了雙核或者四核的CPU，在體系結(jié)構(gòu)上具有實現(xiàn)并行計算的硬件條件。此外，還有許多支持并行計算的軟件環(huán)境，例如，Windows多線程編程、MPI并行編程和Matlab并行計算工具箱等[17]。

為充分利用計算資源，設(shè)計者需要認識到底層多核的存在，采用并行程序設(shè)計智能優(yōu)化算法，提高算法性能。在陰陽平衡優(yōu)化算法中，所有個體基于超球體的解更新操作具有相互獨立性，具有天然的并行特征。此外，和其它智能優(yōu)化算法一樣，陰陽平衡優(yōu)化算法采用群體搜索策略，具有隱含并行性，例如各個體適應(yīng)度函數(shù)值的評價可以并行化等。在本文中，采用Matlab的并行計算工具箱(Parallel Computing Toolbox)對陰陽平衡優(yōu)化算法進行并行程序設(shè)計。

以算例1為例進行實驗，對比算法并行和串行的計算時間。在Windows10操作系統(tǒng)下，CPU：Intel(R) Xeon(R) E-2186M、32G內(nèi)存、2.90GHz主頻的計算機上運行，調(diào)用六核進行并行計算。算法參數(shù)設(shè)置為：Imin=5，Imax=10，α=25，a=2。此外，δ1與δ2初始值都設(shè)為0.5；為防止步長δ2無限制的增大，超出算法的搜索區(qū)域，將其最大值設(shè)為0.75；最大迭代次數(shù)T=500。

和大多數(shù)智能優(yōu)化方法一樣，陰陽平衡優(yōu)化算法屬于隨機方法，需要多次運行同一算法，統(tǒng)計其優(yōu)化性能。在實驗中，循環(huán)次數(shù)分別為1，100，200，300，400，500，比較并行和串行計算的耗時，結(jié)果如圖3所示。

圖3 算法并行和串行計算時間對比

此外，還計算了這5次的加速比(串行算法的執(zhí)行時間除以并行算法的執(zhí)行時間)，最大加速比為4.5231。采用并行程序設(shè)計算法，能夠充分利用空閑的CPU內(nèi)核，能夠減少運行時間。因此，本文所有算法均采用并行程序設(shè)計，以提高算法運行速度。

4.2 算法優(yōu)化性能測試

為測試陰陽平衡優(yōu)化算法YYPO性能，采用微粒群優(yōu)化算法(particle swarm optimization，PSO)[18]、引力搜索算法(gravitational search algorithm，GSA)[19]、生物地理學優(yōu)化算法(biogeography-based optimization，BBO)[20]和最有價值球員算法(most valuable player algorithm，MVPA)[21]進行比較。

為保證實驗公平性，這5種算法設(shè)置相同的函數(shù)評價次數(shù)，即2000n，其中n為問題規(guī)模。YYPO算法其它參數(shù)設(shè)置保持不變，其它算法根據(jù)文獻[18-21]進行參數(shù)設(shè)置。為避免算法單次運行結(jié)果的偶然性對分析優(yōu)化性能的影響，每種算法各自獨立運行30次，分別統(tǒng)計最優(yōu)值、最差值、平均值和標準差等指標。此外，統(tǒng)計結(jié)果還考慮每個算法在30次實驗中達到已知最優(yōu)解的成功次數(shù)，并進行排序。具體計算結(jié)果如表1所示。

表1 不同算法計算結(jié)果比較

圖4 算法優(yōu)化過程對比

通過分析表1中的結(jié)果不難發(fā)現(xiàn)：對于最小規(guī)模的算例1，在相同的函數(shù)評價次數(shù)下，這5種算法都能收斂到已知最優(yōu)解。但是對更大規(guī)模的算例2和3，雖然微粒群優(yōu)化算法PSO、引力搜索算法GSA、生物地理學優(yōu)化算法BBO和最有價值球員算法MVPA都可以找到已知最優(yōu)值，但是在最差值、平均值、標準差、成功次數(shù)和排序這五個指標上，這些算法都明顯劣于陰陽平衡優(yōu)化算法YYPO。此外，陰陽平衡優(yōu)化算法YYPO所獲得的標準差明顯優(yōu)于其它4種算法，說明本文算法的優(yōu)化性能更加穩(wěn)定。上述結(jié)果說明陰陽平衡優(yōu)化算法YYPO在計算精度方面明顯優(yōu)于其它4種優(yōu)化算法。

為進一步分析算法優(yōu)化性能，圖4給出了這5種算法優(yōu)化速度對比示意圖。從圖中可以發(fā)現(xiàn)，陰陽平衡優(yōu)化算法YYPO不僅具有更高的計算精度，而且具有更快的收斂速度。該算法為最小比率旅行商問題的求解提供一個可行有效的算法。

4.3 靈敏度分析

在最小比率旅行商問題中，兩個城市間的距離和旅行商的收入直接影響總行程和總收益之比。這里，分析距離矩陣D和收益矩陣P發(fā)生變化時對結(jié)果的影響。為方便起見，考慮距離矩陣和收益矩陣中元素同時增加或同時減少的情況。當這兩個矩陣元素取值反向變化時，總行程和總收益之比容易判斷變化趨勢。例如，當距離矩陣元素取值都變大而收益矩陣元素取值都變小時，總行程和總收益之比變大。

這里，以算例1為例，分析這兩個矩陣元素取值同向變化時，總行程和總收益之比變化情況。原總行程為199，總收益為234，總行程和總收益之比為0.8504。

當距離矩陣元素值都增加5，收益矩陣元素都增加1時，總行程變?yōu)?24，總收益為239，兩者比值為0.9372。此時，目標函數(shù)值變差。當距離矩陣元素值都增加10，收益矩陣元素都增加20時，總行程變?yōu)?49，總收益為334，兩者比值為0.7455。此時，目標函數(shù)值變好。

當距離矩陣元素值都減少6，收益矩陣元素都減少2時，總行程變?yōu)?69，總收益為219，兩者比值為0.7717。此時，目標函數(shù)值變好。當距離矩陣元素值都減少8，收益矩陣元素都減少16時，總行程變?yōu)?59，總收益為154，兩者比值為1.0325。此時，目標函數(shù)值變差。

距離矩陣元素值都增加相當于旅行商通勤距離變大，服務(wù)對象擴展到郊區(qū)。此時，旅行商的收益也會增加，但不表示費用效益比也會增加，即總行程和總收益之比可能變大也可能變小。與此類似，距離矩陣元素值都減小相當于交通條件的改善，通勤距離變短。此時，旅行商的收益會減少，但不表示費用效益比也會減少，即總行程和總收益之比可能變小也可能變大。

最小比率旅行商問題的目標是考慮總行程和總收益之比，類似于經(jīng)濟學中費用效益比，與僅關(guān)注總距離相比，既考慮總行程，又考慮旅行商的收益，更具有現(xiàn)實意義。在滿足約束的情況下，如何兼顧成本和旅行商收益需要根據(jù)具體問題進行定量計算，以實現(xiàn)最佳的調(diào)度方案。

5 結(jié)論

為有效求解最小比率旅行商問題，構(gòu)建了陰陽平衡優(yōu)化算法。首先，引入佳點集法來初始化搜索群體，改善算法初始解的性能；其次，采用基于超球體和歸檔集的兩種策略對不同階段的解進行更新，并根據(jù)相對位置索引法進行解的編碼來產(chǎn)生可行解，以滿足模型中所有的約束條件；最后，根據(jù)綜卦變換采用反轉(zhuǎn)操作，避免算法早熟收斂。此外，算法采用并行程序設(shè)計，以充分利用多核處理器計算資源。與其它四種算法相比，在計算精度和優(yōu)化速度兩方面，所設(shè)計的算法均排名第一。通過距離矩陣和收益矩陣的靈敏度分析發(fā)現(xiàn)，配送距離和收益同時增加或者減小時，需要根據(jù)具體問題進行定量計算，才可以確定最小比率(費用效益比)的變化情況。將算法用于車輛路徑問題是進一步的研究工作。