求解旅行商問題的多尺度量子自由粒子優化算法

楊云亭，王 鵬

（1.中國科學院成都計算機應用研究所，成都610041； 2.中國科學院大學，北京100049；3.西南民族大學計算機科學與技術學院，成都610225）

（?通信作者電子郵箱wp002005@163.com）

0 引言

近些年來，元啟發式優化算法［1］發展迅猛，不斷地被挖掘應用到現實中的大規模問題中，如神經網絡模型及參數優化、強化學習中的策略選擇等熱門技術中，相比傳統算法中的精確算法中的動態規劃、分支限定等，更加靈活和通用。元啟發式優化算法求解優化問題的過程中，不受限于目標函數的可導性質和優化問題的對偶規則，就可以對問題進行算法的求解，同時可以在非常大的候選解空間中進行迭代搜索［2］。在通用的啟發式優化算法中，模擬退火（Simulate Anneal，SA）算法［3-5］將優化問題看作是物理退火的過程，根據Metropolis準則選擇是否接受搜索到的新解，直至退火過程結束；遺傳算法（Genetic Algorithm，GA）［6-8］在當前搜索到解的基礎上進行自然法則的變化：遺傳、變異、交叉產生新解，在新解中應用“物競天擇，適者生存”法則進行最優解的替換；粒子群優化（Particle Swarm Optimization，PSO）算法［9］則根據鳥類的飛行規律進行解空間中的搜索，通過不同位置的速度和位移變化進行新解的搜索和替換。這些不同的搜索方式會使解沿著某種方向趨近于全局最優解，當搜索足夠充分的情況下，這類啟發式算法會以概率1收斂于全局最優。本文模擬量子體系下自由粒子的波函數的概率解釋進行優化算法搜索行為的指導，提出了多尺度量子自由粒子優化算法（Multi-scale Free Particle Optimization Algorithm，MFPOA）。

組合優化問題是現實生活中很多問題的抽象，求解此類問題可以借助元啟發式優化算法進行離散狀態空間中的搜索。本文使用MFPOA對組合優化問題中的旅行商問題（Travelling Salesman Problem，TSP）［10-12］進行求解，首先對MFPOA進行了物理模型和基本流程進行了闡述，并采用了光骨頭煙花算法（Bare Bones Fireworks Algorithm，BBFWA）［13］中的尺度迭代策略對算法進行了改進，將算法應用到TSP數據集上進行優化問題的求解。另外，通過與目前較優的啟發式優化算法中的混合粒子群優化（Hybrid Particle Swarm Optimization，HPSO）算法［14］、模擬退火算法、遺傳算法和蟻群優化（Ant Colony Optimization，ACO）算法［15］進行TSP求解能力的對比，進一步驗證自由粒子模型下優化算法的性能。使用量子概率驅動解的搜索過程是一種創新的方式，同時完善的量子理論可以支撐算法的求解可行性。

1 基于量子概率的MFPOA

1.1 MFPOA的物理模型描述

多尺度量子諧振子優化算法（Multi-Scale Quantum Harmonic oscillator Optimization Algorithm，MQHOA）［16-18］是一種基于一維量子諧振子波函數原理提出的優化算法，模擬量子系統下的諧振子運動規律，將優化問題的求解通過量子諧振子過程和多尺度過程實現全局區域的逐步搜索定位和局部區域的充分搜索。在該算法的基礎之上，對量子微觀系統下的自由粒子的運動方式進行了分析，并將自由粒子的運動規律映射到函數優化過程中。物理學中的薛定諤方程可以通過給定的勢阱函數求解出粒子的波函數，進一步通過波函數平方計算得到粒子在勢阱下的概率密度函數。而復雜的勢阱較難通過薛定諤方程求解出波函數，在此通過勢阱函數的泰勒展開公式的近似求得近似波函數，在波函數的指導下進行可行解空間的采樣搜索。在優化問題中將優化問題的目標函數近似薛定諤方程中的勢阱，采用勢阱函數的零階泰勒公式近似，可以得到簡潔的勢阱表達，然后通過薛定諤方程求解出近似波函數形式。借助物理學中波函數的概率意義，指導優化算法在空間中進行采樣。

多尺度量子自由粒子優化算法是從微觀系統中的自由粒子衍生出來，自由粒子是量子系統下簡單的粒子之一，在空間中不受力的作用，并且沒有能級結構，轉化到優化求解過程中，可以實現簡潔的搜索結構。目標函數的零階泰勒展開式為常數，在薛定諤方程中使用0表示勢阱，這與自由粒子在空間中不受力等效。首先將自由粒子的勢阱0帶到薛定諤方程中得到波函數的表達，然后通過波函數的平方表征粒子在空間中出現的概率，該概率形式指導算法在當前搜索的最優解和當前的搜索尺度下依此概率進行迭代搜索。將量子空間中簡單的粒子形式映射到優化算法中，通過均勻分布采樣和尺度收縮在求解過程中實現當前最優解的替換，當搜索足夠多次后實現最優解的收斂。

量子理論中的不受力的自由粒子的薛定諤方程寫作如下形式：

其中：?為普朗克常數，m為粒子的質量，ψ(x)為粒子的波函數，E為粒子的能量。

通過方程求解的波函數如下所示。

粒子在空間中的存在是不確定的，但是可以通過波函數模的平方表示粒子在空間某位置上的概率，如式（3）：

式（3）可以看出自由粒子的波函數的概率意義表示粒子在空間中任何位置存在的概率是一個常數，也即在空間任何位置中存在的概率相同，在算法求解過程中可以使用均勻分布進行近似。MFPOA模擬粒子的運動過程，采用均勻分布在可行域中進行采樣，根據采樣點在優化函數中的函數值的信息選取下次迭代過程中的采樣中心，并根據多次采樣信息決定尺度的收縮，直到搜索過程中達到結束條件，退出迭代。

1.2 MFPOA的基本工作流程

MFPOA在算法搜索過程中是逐漸減小搜索區域的，直至搜索區域減小到給定搜索區域的最小值，說明算法求解已經收斂到接近全局最優的程度，算法結束。算法基于概率的形式進行采樣，在某一個尺度是否達到充分的搜索只能通過一些預設條件的判斷，但是這種判斷一般都是存在一定誤差，在判定充分搜索后可能還存在跳過的區域?；谶@種判斷的誤差性，本節介紹了一種自適應的尺度變化策略，這種策略不同于預先設定迭代充分搜索的次數，而是使用搜索過程中的信息進行尺度的持續自適應變化。模擬量子自由粒子在空間中的存在概率，MFPOA通過概率采樣在可行解空間中進行搜索，并實現優化問題的求解，其偽代碼如下所示。

算法1 MFPOA。

算法初始化自由粒子的個數k，搜索空間的最小值min，最大值max，最后停止迭代的搜索空間的最小值δmin，初始的搜索區間在[min，max]，和充分搜索的判定條件c值，實驗一般設定c值為10，f(x)表示目標函數。2）～3）行表示隨機生成k個解空間中粒子，并計算每個粒子的在目標函數中的函數值，并記錄其中最優解f(x)best；4）～12）行表示算法迭代搜索的過程，分別以k個粒子為中心進行均勻分布采樣產生新的粒子，采樣的尺度跟當前的δ有關。采樣后的粒子在目標函數中表現更優的情況下進行粒子的替換，并將當前粒子中的最差用當前k個粒子中的最好粒子更新。若此次迭代過程中的最優值相較上次迭代過程中的最優值沒有發生變化，則最優值的優越性系數為2，下次迭代若最優值依然不變，優越系數增加1，如此累計，直到優越系數等于c時，表示此搜索區域已搜索充分，搜索區域減半，進行更精細的搜索。否則，在當前尺度下進行最優值替換，進一步地迭代搜索，尺度不減小。

模擬自由粒子在空間中存在的概率分布的形式，MFPOA通過均勻分布采樣，探索可行域中最優解的存在位置，不斷地通過迭代進行中心替換和搜索尺度變化，逼近可行域中的最優解。這種方式是將最終解的形式看作是一種在特定區域的概率分布，而不再是一個特定的值。當無限次迭代后，可行解會慢慢逼近到一個極小的區域中，回到最終的值的形式。從偽代碼中也可以看出，算法中控制搜索區域變化的參數c決定整個算法求解在當前搜索尺度下收斂的情況。c值設置得過小，在當前尺度下進行采樣搜索不夠充分，進行尺度減半會丟失大部分解空間，導致算法求解性能下降。c值設置得過大，會使算法在解空間中進行大量采樣和最優值判斷及替換，大幅增加算法求解的時間，而使算法求解速度不可忍受。c值的設定使算法增加了可調參數的設置，c值設置的不合理性會嚴重影響算法的收斂能力和求解精度，一般設置為10。

基于此，本文參考了2018年提出的BBFWA中搜索區域變化的策略。BBFWA在煙花算法的基礎上簡化了算法求解的結構，對煙花爆炸的半徑距離進行了動態的調整來提高算法的求解效果。煙花算法模型同樣采用了均勻分布的采樣方式進行可行解的探索，不同的是MFPOA在搜索過程中通過判斷當前尺度搜索的充分性進行尺度的縮小，而BBFWA在搜索過程中根據采樣前后最優值的變化情況進行尺度的更新。這種尺度的自適應更新實質是對尺度持續衰減的一種改變，在每次采樣完后進行在原尺度的更新，在最優值變得更優的情況下，下一次采樣尺度在當前尺度上擴大至原來的1.2倍；否則，在當前尺度上縮小至原來的0.9倍。這種尺度的變化更加緩慢，相對于尺度減半策略不會丟失解空間中的大部分信息。

為了減少MFPOA中可調參數c值的影響和避免尺度減半導致的解空間的信息的丟失，采用了BBFWA中煙花爆炸半徑的改變策略對算法的搜索尺度進行了一種改進，使用迭代過程中的信息指導尺度的減小，去掉c值的判定條件，采用量子自由粒子的概率解釋在搜索空間中尋優的方式是一種新型搜索算法的框架。BBFWA實質上是對自由粒子模型下的算法進行了搜索尺度的改進，在MFPOA簡單結構的基礎上進一步地改進策略的指導，可以產生性能更優的算法，同樣自由粒子算法可以對同樣機理搜索的算法進行可行性的解釋。光骨頭煙花算法在自由粒子優化算法的解釋下的偽代碼可以進一步改寫成如下形式，并記為MFPOA-DS（Multi-scale Free Particle Optimization Algorithm-Dynamic Search）。

從算法1和算法2中的兩段偽代碼中可以明顯看出，MFPOA-DS版本中尺度的更新使用了公式δ=δ*λ，其中λ的取值有兩種0.9和1.2：只有迭代中產生的新粒子比上一次迭代中的粒子的最優值好，設置λ=1.2，進行搜索區域的擴大；否則進行搜索區域的衰減。MFPOA減半的搜索方式使得算法搜索過程中會忽視搜索中心外圍區域，過于關注內部區域，這種方式會導致可行解空間中的部分信息丟失，從而對優化問題的求解不準確。這種策略引入動態尺度搜索的概念，將前后兩次搜索到的解的變化情況代替尺度收斂的控制條件值c的判斷，可以使尺度靈活縮放，進而改善算法的性能。

2 自由粒子模型求解TSP的基本過程

TSP一直是組合優化問題中經典的NP難問題。TSP的描述是給定若干城市和城市之間的距離，通過算法求解訪問每座城市并返回出發城市的最小回路距離，也稱為最短路徑問題。啟發式優化算法往往更擅長連續空間中的優化問題，在求解組合優化問題時需要先對組合優化問題進行合適的映射，進而按照函數優化求解的思路進行求解。相比連續空間中的函數優化，組合優化問題需要先選擇優化問題解的編碼方式以及可行解的更新方式。TSP解的編碼與給定的城市數目有關，路徑距離基于城市序列計算歐幾里得距離得出，所以城市序列的表示尤為重要，簡單的一種表達方式是對城市進行編號，比如5個城市，分別給每個城市編成1，2，…，5。城市序列123451便表示TSP的一個可行解，表示從城市1出發，分別經過城市2，3，4，5，再回到城市1。在編碼方式的基礎之上，通過隨機選擇兩個位置進行交換產生新解。

TSP作為離散問題求解的代表，一直是組合優化問題中經典問題。如何求解城市序列的最優路徑，使得旅行商經過的總路程最小，也一直以來被用來測試啟發式算法的性能。MFPOA求解TSP的策略主要包含兩個步驟：1）在當前尺度下進行城市序列的迭代并記錄最優解的情況和采樣中心的變化；2）在不同尺度下進行城市序列在已搜索到最優區域的精細搜索。根據多尺度量子自由粒子模型，MFPOA-DS求解TSP的基本工作流程的偽代碼如下所示：

由偽代碼和算法的基本思想，可以將自由粒子模型求解TSP的基本過程概述為以下3個主要步驟：

1）隨機初始化k個城市序列，2，…，k）作為當前尺度（初始尺度）下的搜索中心。

2）在當前尺度δ上的收斂過程，以每次迭代過程中保留的k個較好的采樣位置作為k個局部搜索區域的中心，形成k個區間為δ的均勻分布函數U(xi-δ2，xi+δ2)。分別根據每個中心序列進行隨機位置的交換產生新解，在迭代過程中記錄最優解和最優解的位置，實現單一均勻分布的收斂，多個均勻分布采樣迭代的過程實現當前尺度下的基本收斂。

3）搜索區域的聚焦過程，在上一尺度搜索下達到基本收斂的情況下，根據實際求解到的最優解的變化情況進行當前尺度的變化，實現局部區域的精細搜索，直到搜索區域減小到足夠小或者得到預設定的迭代次數。

當給定n個城市的坐標時，可以通過歐幾里得距離計算城市之間的距離?？尚薪饪臻g表示為n個城市的全排列，一共有n！種方案，算法的求解便在這些方案中進行搜索求解。一一窮舉這些解在城市數目很多的情況下是不可行的，啟發式算法在求解過程中根據當前最優解的分布進行下一步的搜索，有選擇地進行搜索區域的迭代來避免全部空間的窮盡搜索。在求解TSP時，算法中δ的初始值設為城市總數，δmin設置為1，在[1，δ]區間內進行均勻分布隨機位置的產生，當δ小于1時沒有意義。l(x)表示城市序列的路程長度，從一個城市出發，歷經所有城市后回到原始出發的城市的路徑長度。在MFPOA求解中δ隨著迭代過程按照固定倍數進行減小，在MFPOA-DS中δ隨著迭代過程中搜索到的解的情況進行自適應的尺度增減。這一尺度的變化控制著算法求解TSP的收斂情況，當尺度縮小到一定程度(δmin=1)，算法終止。

MFPOA本身模擬量子系統下的自由粒子的特點進行優化問題中的可行域的搜索，搜索策略簡潔，通過均勻分布采樣和最優值替換得到最終整體的最優解。同樣在TSP中通過簡單的城市序列變換在可變的搜索空間中進行最優解的搜索。當迭代足夠充分，算法求解的最優城市序列會逼近最優解。

3 仿真實驗與分析

本章使用組合優化問題中代表性的TSP測試自由粒子模型算法的性能。連續優化問題更容易針對不同的問題在合適的范圍進行采樣，組合優化問題需要先進行合適的映射，將問題中的求解結果用合適的序列表示，同時采樣過程需要根據問題制定。使用自由粒子優化算法求解TSP時，將問題映射成一個在初始解的基礎上進行交換位置的均勻分布采樣產生新解，根據新解比當前解的效果進行采樣尺度的變化（采樣尺度初始為城市序列中的城市個數）。在MFPOA中的尺度變化根據迭代過程中的解不發生變化的次數進行尺度的減半策略。在尺度自適應改變版本MFPOA-DS的算法中，根據搜索中的解的變化進行采樣尺度的放縮變化，相比原始做法，降低了算法的復雜度，使算法更加靈活。

所有仿真實驗均在Intel Xeon CPU E5-1630 v3 3.7 GHz，16 GB內存的計算機上進行，程序采用Matlab 2018a實現。經過反復實驗測試，在本次TSP上測試過程中使用參數設置如下，初始粒子個數k為30，初始的搜索尺度δ為求解TSP時的城市數目，算法停止迭代的條件δmin設置為1，最大迭代次數w0設置為30 000，搜索尺度變化系數λ初始為1，在迭代過程中隨最優解的變化在1.2和0.9之間變化。城市間的距離采用歐幾里得距離進行度量。

3.1 尺度自適應變換策略的有效性

實驗采用不同的對稱性TSP的數據驗證算法求解的可行性。實驗采用了 ulysses16、ulysses22、eli51、burma14、bayg29的實驗數據，通過這些數據測試了MQHOA和使用尺度自適應改進策略的MQHOA（MQHOA-DS）、MFPOA和MFPOA-DS。由于算法的隨機性，采用了10次重復實驗數據進行對比，實驗數據如表1。

表1 10次重復實驗在不同實驗數據上的結果對比Tab.1 Result comparison of 10 repeated experiments on different experimental data

MFPOA在MQHOA的啟發下，考察了量子系統中的不同粒子在微觀系統下的運動形式對算法指導優化求解過程的影響，主要考察了10次獨立重復實驗的最優值（Global Best，GB），10次實驗結果的平均值（Mean Best，MB）和10次實驗結果的平均值與最優值間的差距（MB-GB）。對比MFPOA和MFPOA-DS，MFPOA-DS在5種數據集上性能比MFPOA都有相應的提升。根據迭代過程中得到的解的情況進行尺度的縮減和放大，可以有效提高算法的搜索能力。而且通過觀察重復實驗中最好值和平均值的差距，可以發現使用自適應尺度變化策略后的算法會減弱算法每次運行的隨機性，每次實驗的差距變小，說明實驗的性能更加穩定。對比基于量子諧振子的概率進行尋優的算法，自由粒子模型求解的結果更好，在eli51和bayg29數據上的數據對比尤為明顯。這從側面驗證了自由粒子中的搜索模型更適合求解TSP，均勻分布采樣過程在求解組合優化問題上比高斯采樣具有一定的優勢。在TSP中，由于解空間的離散化很難保證解之間的相互關系，這種相互關系在函數優化問題中體現得較為明顯，一般最優解和次優解緊挨或者距離很近。這也就導致了算法求解過程中依高斯概率在當前次優解周圍采樣后的新解，回路距離減少不代表新解跟最優解的相對距離更近。這種相對距離可以理解成城市序列在相同位置上不一樣的城市的個數。

3.2 拓展實驗

本節是對自由粒子算法框架改進算法和其他啟發式優化算法的對比實驗分析。對比實驗分析了自由粒子模型的算法跟HPSO、SA、GA、ACO性能。HPSO混合了遺傳算法中的交叉變異操作，將標準粒子群中的速度更新和位置更新替換為粒子與個體最優和全體最優進行交叉的過程，算法的初始參數設置如下：粒子數目為1 000，進化次數為1 000；SA的初始參數設置如下初始溫度為1000，終止溫度為0.001，在每個溫度下迭代500次，降溫速率設置為0.9；GA的初始參數設置如下：種群個數為100，最大遺傳代數為1 000，交叉概率為0.9，變異概率為0.05，選擇進行交叉和變異的染色體數目的控制系數（有時被稱為個體間的代溝）為0.9；ACO的初始參數設置如下：螞蟻數量為50，信息素重要程度因子為1，啟發函數重要程度因子為5，信息素揮發因子為0.1，最大迭代次數為1 000。在這些參數下相應的算法均能求解出相對較好的性能解。

為保證算法的穩定性，分別對6種算法進行了10次重復實驗，對比了6種算法的最優值如表2所示，平均花費時間如表3所示，表中MV（Mean Value）表示10次重復實驗中求解結果的平均值，MT（Mean Time）表示10次重復實驗中求解所花時間的平均值。

表2 不同算法求解的最優值結果對比Tab.2 Comparison of theoptimal valuesobtained by different algorithms

表2中加粗數字表示不同算法求出的解中最好的結果。對比一些其他的優化算法，可以看出自由粒子模型算法的求解效果與HPSO、AS、GA、ACO在eli51、bayg29數據集上的求解精度有些差距，在burma14、ulysses16、ulysses22數據集上求解結果更好。在這些數據集中的MFPOA-DS比MFPOA求解的結果都是更優的，自適應尺度的變化有效提高了算法的求解性能。內部機制搜索的有效性保證了算法求解的可行性，自由粒子模型是一種基礎的骨架模型，算法性能在此基礎上通過一些策略機制的加入，可以進一步提高算法的性能。以上是對算法求解精度層面進行的對比，接下來對算法的求解時間上進行對比。

表3 求解TSP實例時不同算法求解時間的對比 單位：sTab.3 Solving timecomparison of different algorithms in solving TSPexamples unit：s

表4表示在不同的測試集上，MFPOA-DS相對HPSO、AS、GA、ACO在求解時間上減少的百分比，分別用IT1、IT2、IT3、IT4表示：

IT2、IT3、IT4同IT1。由表4可見，MFPOA-DS并不比AS算法有時間上的優勢，所以表格中的數據為負值。

從相對節省時間數據可以得出，AS比其他算法求出最優解的求解時間最少，在時間復雜度上求解更優，在時間效率上MFPOA-DS僅次于 AS。相比HPSO、GA和 ACO，MFPOA-DS在5個測試集上將運行速度分別平均提升了80.42%、55.72%、66.45%?？梢娀谧杂闪Ｗ雍唵谓Y構構造的優化算法，具有更小的時間復雜度，能夠快速、相對高效的求解組合優化問題。

值得注意的是，自由粒子優化算法結構簡潔，在求解TSP時卻有著比較好的表現。這種結構更適合對領域定義不是非常明確的問題，MFPOA的均勻分布結構對任何可行解都是同等看待，會使在可行空間上的搜索更加充分。后續工作會繼續驗證。

表4 不同測試集下MFPOA-DS相對其他算法在求解時間上減少的百分比 單位：%Tab.4 MFPOA-DSreducesthepercentageof thesolution timeof other algorithms respectively under different test sets unit：%

4 結語

多尺度量子自由粒子優化算法模擬了量子系統下自由粒子的簡單運動過程，在優化問題中通過粒子在空間中的運動進行可行域中解的采樣，根據采樣信息對搜索空間進行縮放控制，實現在同一搜索尺度下的精細搜索和不同尺度下搜索中心的聚焦。自由粒子優化算法作為基于量子算法模型中的骨架模型，結構簡潔，易于實現，通過對算法進行策略的改進，會使得算法的性能得到很好的提升。在自由粒子模型下的零勢阱和諧振子模型下的高斯勢阱的啟發下，將優化問題中的目標函數進行勢阱等效和泰勒近似實現量子模型下的轉化，通過粒子在空間中的存在形式指導優化算法的尋優過程，實現量子模型的搜索機制。MFPOA和MFPOA-DS求解TSP的結果可以證明自適應尺度改進策略的有效性，MFPOA-DS改進算法與HPSO、AS、ACO和GA在TSP上的求解數據表明，MFPOA-DS算法在不同數據集上的求解結果表現不同，但在保持較好結果精度的情況下，在求解速度上比HPSO，ACO和GA平均提升50%，但是在求解精度方面算法并不總是比其他的優化算法最好，如何使算法的求解精度進一步提升成為下一步研究工作中的重點。

MFPOA的搜索結構簡單，作為優化問題求解的一種基本的搜索框架，能夠衍射出很多基于均勻分布為采樣方式的算法。未來的工作將重點對算法的改進和算法間的聯系進行深入研究，找出算法在不同應用上的優勢。