【建模算法】基于粒子群算法（PSO）求解TSP問題（Python實現）

TSP (traveling salesman problem，旅行商問題)是典型的NP完全問題，即其最壞情況下的時間復雜度隨著問題規模的增大按指數方式增長，到目前為止還未找到一個多項式時間的有效算法。本文探討了基于粒子群算法求解TSP問題的Python實現。

一、問題描述

? 本案例以31個城市為例，假定31個城市的位置坐標如表1所列。尋找出一條最短的遍歷31個城市的路徑。

城市編號X坐標Y坐標城市編號X坐標Y坐標

1	1.304	2.312	17	3.918	2.179
2	3.639	1.315	18	4.061	2.37
3	4.177	2.244	19	3.78	2.212
4	3.712	1.399	20	3.676	2.578
5	3.488	1.535	21	4.029	2.838
6	3.326	1.556	22	4.263	2.931
7	3.238	1.229	23	3.429	1.908
8	4.196	1.044	24	3.507	2.376
9	4.312	0.79	25	3.394	2.643
10	4.386	0.57	26	3.439	3.201
11	3.007	1.97	27	2.935	3.24
12	2.562	1.756	28	3.14	3.55
13	2.788	1.491	29	2.545	2.357
14	2.381	1.676	30	2.778	2.826
15	1.332	0.695	31	2.37	2.975
16	3.715	1.678

二、粒子群算法主要步驟

初始化粒子群。包括粒子的初始位置及速度，慣性因子等參數值

粒子數M一般選取20~40個，不過對于一些特殊的難題需要更多的粒子，粒子數量越多，搜索范圍就越廣，越容易找到全局最優解，但是也會帶來更大的計算消耗。

評價各個粒子的初始適應值。

將初始的適應值作為各個粒子的局部最優解，保存各粒子的最優位置。并找到其中的最優值，作為全局最優解的初值，并記錄其位置

更新粒子速度及位置

計算更新后粒子的適應值，更新每個粒子的局部最優值以及整個粒子群的全局最優值。

重復4~5直至滿足迭代結束條件

流程圖如下：

三、粒子群算法求解TSP問題

1.TSP問題

**旅行商問題（Traveling Salesman Problem，TSP）**是一個典型的NP問題。

問題描述：一個銷售員需要訪問31個城市，恰好訪問每個城市一次，并最終回到出發城市，銷售員希望整個旅行的路程最短（或者成本定義為旅行所需的全部費用是他旅行經過的的各邊費用之和，而銷售員希望使整個旅行費用最低）。從圖論的角度來看，該問題實質是在一個帶權完全無向圖中，找一個權值最小的Hamilton回路。

2.位置更新（position update）

在TSP問題中，我們所求的是一個序列，需要優化的是這個序列的順序，粒子$x_i$ 如何定義是關鍵，我們將其定義為當前粒子 i所走的順序序列，如：

$x_i=(1,3,5,2,4)$代表著該TSP問題一共有五個城市需要訪問，當前粒子i的訪問順序為1→3→5→2→4，

因此全部城市的訪問順序的所有可能序列就構成了問題的可行空間，**粒子位置的更新就意味著粒子$x_i$ 從一種序列$s_i$ 變化成另外一種序列 $s_j$。

3.速度更新（velocity update）

四、求解結果

最優路線與最優值：

最優軌跡圖：

優化過程：

當然，PSO是一個基于概率的隨機自搜索算法，每次的搜索結果都不盡相同，但是若算法的收斂性表現的較好，也可以認為對算法的設計是合理的。

五、Python源代碼

#粒子群算法求解31座城市TSP問題完整代碼： import numpy as np import matplotlib.pyplot as plt import seaborn as sns import time#計算距離矩陣 def clac_distance(X, Y):"""計算兩個城市之間的歐氏距離，二范數:param X: 城市X的坐標.np.array數組:param Y: 城市Y的坐標.np.array數組:return:"""distance_matrix = np.zeros((city_num, city_num))for i in range(city_num):for j in range(city_num):if i == j:continuedistance = np.sqrt((X[i] - X[j]) ** 2 + (Y[i] - Y[j]) ** 2)distance_matrix[i][j] = distancereturn distance_matrix#定義總距離(路程即適應度值) def fitness_func(distance_matrix, x_i):"""適應度函數:param distance_matrix: 城市距離矩陣:param x_i: PSO的一個解（路徑序列）:return:"""total_distance = 0for i in range(1, city_num):start_city = x_i[i - 1]end_city = x_i[i]total_distance += distance_matrix[start_city][end_city]total_distance += distance_matrix[x_i[-1]][x_i[0]] # 從最后的城市返回出發的城市return total_distance#定義速度更新函數 def get_ss(x_best, x_i, r):"""計算交換序列，即x2結果交換序列ss得到x1，對應PSO速度更新公式中的 r1(pbest-xi) 和 r2(gbest-xi):param x_best: pbest or gbest:param x_i: 粒子當前的解:param r: 隨機因子:return:"""velocity_ss = []for i in range(len(x_i)):if x_i[i] != x_best[i]:j = np.where(x_i == x_best[i])[0][0]so = (i, j, r) # 得到交換子velocity_ss.append(so)x_i[i], x_i[j] = x_i[j], x_i[i] # 執行交換操作return velocity_ss# 定義位置更新函數 def do_ss(x_i, ss):"""執行交換操作:param x_i::param ss: 由交換子組成的交換序列:return:"""for i, j, r in ss:rand = np.random.random()if rand <= r:x_i[i], x_i[j] = x_i[j], x_i[i]return x_idef draw(best):result_x = [0 for col in range(city_num+1)]result_y = [0 for col in range(city_num+1)]for i in range(city_num):result_x[i] = city_x[best[i]]result_y[i] = city_y[best[i]]result_x[city_num] = result_x[0]result_y[city_num] = result_y[0]plt.xlim(0, 5) # 限定橫軸的范圍plt.ylim(0, 4) # 限定縱軸的范圍plt.plot(result_x, result_y, marker='>', mec='r', mfc='w',label=u'路線')plt.legend() # 讓圖例生效plt.margins(0)plt.subplots_adjust(bottom=0.15)for i in range(len(best)):plt.text(result_x[i] + 0.05, result_y[i] + 0.05, str(best[i]+1), color='red')plt.xlabel('橫坐標')plt.ylabel('縱坐標')plt.title('軌跡圖')plt.show()def print_route(route):result_cur_best=[]for i in route:result_cur_best+=[i]for i in range(len(result_cur_best)):result_cur_best[i] += 1result_path = result_cur_bestresult_path.append(result_path[0])return result_pathif __name__=="__main__":#讀取31座城市坐標coord = []with open("data.txt", "r") as lines:lines = lines.readlines()for line in lines:xy = line.split()coord.append(xy)coord = np.array(coord)w, h = coord.shapecoordinates = np.zeros((w, h), float)for i in range(w):for j in range(h):coordinates[i, j] = float(coord[i, j])city_x=coordinates[:,0]city_y=coordinates[:,1]#城市數量city_num = coordinates.shape[0]# 參數設置size = 500 #粒子數量r1 = 0.7 #個體學習因子r2 = 0.8 #社會學習因子iter_max_num = 1000 #迭代次數fitness_value_lst = []distance_matrix = clac_distance(city_x, city_y)# 初始化種群各個粒子的位置，作為個體的歷史最優pbestpbest_init = np.zeros((size, city_num), dtype=np.int64)for i in range(size):pbest_init[i] = np.random.choice(list(range(city_num)), size=city_num, replace=False)# 計算每個粒子對應的適應度pbest = pbest_initpbest_fitness = np.zeros((size, 1))for i in range(size):pbest_fitness[i] = fitness_func(distance_matrix, x_i=pbest_init[i])# 計算全局適應度和對應的gbestgbest = pbest_init[pbest_fitness.argmin()]gbest_fitness = pbest_fitness.min()# 記錄算法迭代效果fitness_value_lst.append(gbest_fitness)# 迭代過程for i in range(iter_max_num):# 控制迭代次數for j in range(size):# 遍歷每個粒子pbest_i = pbest[j].copy()x_i = pbest_init[j].copy()# 計算交換序列，即 v = r1(pbest-xi) + r2(gbest-xi)ss1 = get_ss(pbest_i, x_i, r1)ss2 = get_ss(gbest, x_i, r2)ss = ss1 + ss2x_i = do_ss(x_i, ss)fitness_new = fitness_func(distance_matrix, x_i)fitness_old = pbest_fitness[j]if fitness_new < fitness_old:pbest_fitness[j] = fitness_newpbest[j] = x_igbest_fitness_new = pbest_fitness.min()gbest_new = pbest[pbest_fitness.argmin()]if gbest_fitness_new < gbest_fitness:gbest_fitness = gbest_fitness_newgbest = gbest_new fitness_value_lst.append(gbest_fitness)# 輸出迭代結果print("最優路線：", print_route(gbest))print("最優值：", gbest_fitness)#繪圖 sns.set_style('whitegrid')plt.rcParams['font.sans-serif'] = 'SimHei' # 設置中文顯示plt.rcParams['axes.unicode_minus'] = Falsedraw(gbest)plt.figure(2)plt.plot(fitness_value_lst)plt.title('優化過程')plt.ylabel('最優值')plt.xlabel('迭代次數({}->{})'.format(0, iter_max_num))plt.show()

總結

以上是生活随笔為你收集整理的【建模算法】基于粒子群算法求解TSP问题（Python实现）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。