Python — 人工蜂群演算法(Artificial Bee Colony, ABC)求解旅行推銷員問題(Traveling Salesman Problem, TSP)
人工蜂群演算法（Artificial Bee Colony, ABC）是由 Karaboga於2005年提出的全局優化演算法，主要來源自蜂群採蜜的行為，是屬於群體智能算法的其中一種。
何謂人工蜂群演算法(Artificial Bee Colony, ABC)?
ABC是一種全局優化演算法，來源自蜂群採蜜的行為，蜜蜂根據各自所擔任的角色分別進行不同的活動，訊息會在蜂群間共享和交流，進而找到問題的最佳解。優點為穩定性高、控制參數較少、較易理解，缺點則為存在著過早收斂的疑慮。
筆者曾對於ABC進行過詳細的講解，並用以求解最佳化問題，若是有興趣深入了解，可以參考筆者以下這篇文章。
何謂旅行推銷員問題(Traveling Salesman Problem, TSP)?
此部分主要來自搜尋維基百科，在此只會針對問題定義上大致做一個講解，若想詳細了解可以點擊以下連結，但相信大部分會點擊此篇文章的讀者，應都已具備相對應的觀念。
TSP定義為，給定一系列城市和每對城市之間的距離，求解存取每一座城市一次並回到起始城市的最短路徑。它是組合最佳化中的一個NP困難問題，在作業研究和理論電腦科學中非常重要。
在本文你將會學會以下技巧：
1.培養分析問題的能力。
2.以Python實現簡單的演算法。
3.了解旅行推銷員問題。
4.使用蜂群演算法解決旅行推銷員問題。求解TSP問題簡述
ABC屬於全局優化演算法，所以以下求解的流程採用方法為，優先生成一距離矩陣（distance matrix），此矩陣是求解最佳路徑中，所需點中彼此間的距離，並以此矩陣為基礎，在每個點不重複的情況下，求出一總距離最短的路徑。
嚴格來說，ABC演算法並不是一個求解TSP問題最佳的選擇。但是在實務中，解決任何問題往往需要採取大量的方法才能求得最佳解，所以在以下本文中會告訴大家筆者所使用的方法外，也會在最後總結其適合使用的情況和條件。
以下正文開始～
決定城市數量並生成距離矩陣（Distance Matrix）
首先，決定本次要求解的城市點數，以下設定為10，意指在10個城市間選擇出距離最短的路徑解法。可以看到以下是10x10的距離矩陣，分別代表第 i 個城市到第 j 個城市的距離，以及相對應的第 j 個城市到第 i 個城市的距離。
基本參數設定
- D、 N：維度以及組數，藉由此生成 Food Source。
- ub、lb：區域上下界。
- Ub、 Lb：全域上下界 。
- limit：試驗上限，用以決定是否放棄該蜜源。
- exchange_num : 生成 Xnew 時，Food Source中一次要轉換的個數。
- max_iter：迭代次數。
如何求解？
和上一篇介紹的使用ABC求解最佳化問題所採用的方法不同，在這次要解決的問題是TSP，必須要將 Food Source 的概念轉換成路徑，並且計算每次路徑的距離，才能在迭代過程中求得最短路徑，進而得出整體最佳解。
老實說我一開始也沒想到這個解法，網路上也幾乎沒什麼人針對這一方面下去講解其關連性和解法，直到筆者一位很厲害的朋友告訴我後才恍然大悟，以下介紹使用此方法的條件以及限制。
實數解取代路徑手法
此法實際上就是針對所需的路徑城市點數量，生成一相對應的點數的亂數數列，並且按照數列中數字的大小進行排列，就是為所求的路徑。
條件及限制
- 實際轉換後的路徑必須是0開頭0結尾的結構：TSP求解大多在最後會回到原出發點，所以在實際解中，如果是一10個城市的TSP，在此並不是考慮10個點而是一次考量11個點下去生成串列，進而產生對應頭尾相同的路徑。
生成的數列，頭尾數字一定須相同，代表從相同城市出發再回到相同的城市。
[ 0 7 5 6 9 2 4 1 3 8 0 ]- 實數串列開頭數字最小，結尾數字最大：由於實際路徑是按照此實數串列中的大小進行排序的，所以在生成一實數串列時，務必要控制好開頭及結尾的大小，才能確保此路徑在之後任何階段變動都頭尾都不受到影響，進而確保其開頭結尾城市相同的條件。
生成的數列，將頭尾數字控制在(-10, 10)，使其在迭代過程中不會被影響實際順序，其餘數字設定成0~1的均勻隨機亂數，以確保不會有重複值的產生。
[ -10 0.767 0.453 0.532 0.890 0.111 0.325 0.022 0.231 0.800 ]
[ 0 7 5 6 9 2 4 1 3 8 0 ]初始化蜜源
本次實作範例所採用的參數設定如下：
- N：10。
- D：9。上述雖提及有11個點在求解，但頭尾需固定數字，所以一組解中實際有在變動的只有9個點。
- ub、lb：1、 -1。計算Xnew時，藉由區域上下界控制其數字變動大小。另一種則是Xnew計算時也許會超過全域上下界的值，這部分的變數控制在下文中會提及。
- Ub、 Lb：10、 -10。生成實數串列時，頭尾的值。
- limit：D*N。偵查蜂階段蜜源的替換上限。
- exchange_num : 2。
- max_iter：10000。
以下為初始化的蜜源，以實數串列型態表現：
將每組串列按照大小排序，得到各組的路徑，並進而求出其距離和fitness：
Minimum: 482 (Distance)
雇佣蜂階段 (Employed Bee Phase)
其算法流程大致上和在求最佳化問題相同，除了實數數列的置入外，其他詳細流程可參照上一篇，有較詳細的說明，以下只會針對計算實數解的步驟詳細說明。
雇佣蜂階段主要目的為，當前蜜蜂A透過隨機向附近其他隻蜜蜂探查，假設選擇的蜜蜂為蜜蜂B，若是蜜蜂B的蜜源較好，則蜜蜂A會在蜜蜂B附近生成一新蜜源進行採集。
而新蜜源位置的產生，是雇佣蜂根據以下公式進行更新的。算法上為將蜜蜂A中其中一參數選擇出來和蜜蜂B相對應的參數進行Xnew的計算，生成出的Xnew再取代原蜜蜂A中參數位置的值，重新排序藉以求出新的路徑。
其中 ϕ 需特別注意，為了確保其實數數字能在串列中改變位置，所以設置為1~-1之間的亂數。其他組合並沒有不好，但經過不斷的試驗，發現 ϕ 在(1~-1)之間生成的亂數在最後結果中收斂速度較隨機。
原是設定成(0~-0.5)，原因是收斂速度最快，但是後來發現如果按照此設定，數列中的亂數值會被控制在0.4多左右，只要迭代次數一高，就會發生路徑上城市點重複經過的問題。
Xnew上下界值的設定
這裡有一重要問題是如何控制生成後Xnew解的大小，如果超過蜜源中均勻隨機亂數設定值0~1的界限怎麼辦？如果將超過界限的值設定成定值呢？
當然不行，迭代次數少說是介於100～100000，過程中絕對會出現一樣超過界限值的數字很多次。所以在打code時，筆者將其設定成隨著迭代次數的增加，其超過界限的值會根據其迭代的次數進行改變。
設定如下，將超過上界值控制在0.9999999，並且隨著迭代次數(it)和Xnew第幾個運算次數(i)進的微調，可100%確保數字不會重複。
Xnew = num[i]+(random.uniform(lb,ub))*(num[i]-challenge_num_list[i])
if Xnew>=1:
Xnew = np.array(0.9999999-(0.000000001*it)+(i*0.0000000001))
elif Xnew<=0:
Xnew = np.array(0.0000001+(0.000000001*it)+(i*0.0000000001))雇佣蜂過程中一次須交替多少參數？
沒錯，一次要替換多少參數也是必須要考量的問題之一。其替換的多寡將會決定以下兩件事：
- 一次只交換 1 個參數：這是最保險的做法，優點是可確保新生成的路徑變動不會太大，其distance也會穩定收斂。缺點是收斂速度過慢，容易落在局部最佳解，且由於ABC本質上是全局優化演算法，只交換 1 個效率會過低，如果電腦性能較差，迭代次數一大就會吃不消。今天示範範例是10個城市點，實務絕對是上百上千個，如果一次只交換 1 個城市點的參數，那真的別鬧了，你絕對是在瞧不起蜂群演算法。
- 一次交換多個參數：變動太大，雖說ABC的優點是收斂快速，前期是可快速收斂到定值，但到後期就會開始變動，主要是取決交換參數的數量，經過測試，如若是 3 個以上，完全不會收斂，控制試驗次數上限(limit)也沒用，除非迭代次數(iteration)設小才有可能隨機找到最佳解，注意是“隨機”，代表這是一件非常吃運氣的事情。
經過筆者的試驗，對10個城市點來說，一次交換 2 個參數可大大達到上述兩點所提及的效益，收斂快速又可以很快找到最佳路徑，到後期也不大影響收斂結果。
如要找出城市點和交換參數數量的關係，需要進行大量的測試才可找出一最佳組合關係，因本文範例是10個城市點，故只會採用 2 個參數進行實作。順帶一提，在測試100個點時，交換參數值一超過 2 就不會收斂了。
經過上述所有限制的介紹，現在又回到雇佣蜂階段的部分繼續講解。其新生成的路徑，若是distance的fitness較大，則將其路徑取代；反之若是新生成distance的fitness較小，則維持原路徑，並且試驗次數(Trial)加 1。
關於蜜源豐富程度(fitness)大小和試驗次數(Trial)的規範
TSP問題絕對為“求解最小化問題”的一種，故在寫code時需讓演算法永遠保留最小的值，所以針對fitness公式來說：
f(x)越小，代表fitness越大；f(x)越大，代表其fitness就越小
故為了只保留最小值，若生成的f(x)沒有比較小，則保留原值，trial必須加1，反之若是生成的f(x)有比較小，則將其替換且trial保持不變。觀察蜂階段(On-looker Bee Phase)
此部分主要是蜜蜂依據其蜜源豐富程度的大小，採用輪盤法(Roulette Wheel Selection)決定是否要進行蜜源的更換。簡單來說就是蜜源如果採集次數過多就會枯竭，此時就必須要決定是否需要更換至其他蜜源進行採集，此更換的流程和雇佣蜂過程相同，在此就不多以贅述。
建立Probability list
透過上述第一階段結束的fitness求出各別的機率p值，並建立一機率列表。
針對每一組數依序生成一0~1之間的亂數 r，透過和其組數對應的機率p比大小，若是 r 比 p 小，則進行到觀察蜂階段，若沒有則保持不變，且Trial增加1；由上述可知，若是fitness的機率p越大，則可能需替換的機率就會越大。
觀察蜂階段
其主階段內容和雇佣蜂相同，以下直接附上代碼，值得一提的是在設定Xnew超過上下界所控制的值，避免和雇佣蜂階段產生的Xnew重疊而改成如下：
基本上數值只是進行微調而已，目的是為了不重複而已，但又不能控制得太誇張。
Xnew = num[i]+(random.uniform(lb,ub))*(num[I]-challenge_num_list[i])
if Xnew>=1:
Xnew = np.array(0.9999998-(0.0000000011*it)+(i*0.0000000001))
elif Xnew<=0:
Xnew = np.array(0.0000001+(0.0000000011*it)+(i*0.0000000001))以下為判斷是否需要進入觀察蜂階段，會應用到上述提及的probability_list和Onlooker_Bee_Phase函式，通常在此階段就會儲存每次迭代的最佳解。
for i in range(N):
if random.uniform(0,0.1) < Probability_list[i]:
onlooker_out = Onlooker_Bee_Phase(Parameters)
foodsource_realvalue2 = onlooker_out[]
trial_afe = onlooker_out[]
else:
foodsource_realvalue2 = foodsource_realvalue
trial_afe[i] += 1偵察蜂階段(Scout Bee Phase)
偵察蜂判斷每個蜜源是否到達試驗上限，若達到則放棄該蜜源，並生成一蜜源代替。
這部分需要針對各組的Trial值(試驗次數)判斷是否有超過limit上限，若是有 則隨機生成全域範圍內的一組全新的Food Source，因這部分只需做簡單的判斷即可，所以就不以gists呈現。
#Scout bee phase
for i in range(N):
if trial_afe[i] > limit:
foodsource_realvalue2[i] = np.array(foodsource(D)[0])
trial_afe[i] = 0迭代結束，輸出最佳蜜源位置。
為判斷是否能求出最佳解，採用窮舉法(Method of exhaustion)，將頭尾都是 0 的路徑中全部都嘗試一遍。
得出最佳路徑和最短距離為：
實際丟進去ABC迭代10000次得出的最佳解如下：
將結果視覺化如下，可看出在本次試驗，於迭代800左右就收斂到最佳解了，且也沒有因為一次交換兩個參數導致收斂線動盪的問題出現。
可看出這次試驗成功，最佳結果有迭代到窮舉法的最佳解。
完整代碼分享
結語
用ABC解決旅行推銷員問題並不是一個好的做法，需要考量的點非常多，但是拿來訓練思維能力還是非常好用的。實務上如果要應用就會較為吃力，要校正的數值也會比較多，像筆者在雇佣蜂階段和觀察蜂階段就花費非常多時間校正“Xnew上下界值”和“一次交換多少值”，避免發生經過相同城市點的問題出現，直到目前才調整出可以求得最佳解且不會發生經過一樣城市點的參數組。
如要利用ABC解決更複雜的旅行推銷員問題，最好的方法建議是增加 N ，也就是初始路徑的組數，但是相對應的要從數量這麼多的路徑迭中代出最佳解就必須使用大量的迭代次數才能解決。
筆者其實正在解決一實務問題，主要是100城市點的旅行推銷員問題，因點數實在過多，採用100組初始路徑並迭代10000次試圖找出最佳解，跑了一天後結果出來卻發現距離還沒收斂，透過此可以發現初始組數初始組數 N 越大就會花更多時間找最佳解，所以後來筆者將迭代次數調整成30000讓演算法下去運算，而目前還正在跑。因為實在是太久了，所以我就抽空寫了這篇文章，結果寫完還在跑。
等最佳結果真的跑出來後，我在抓時間寫一篇看看吧，謝謝大家。
