算法原理

??????旅行商問題是一個經典的NP問題，假設有N個城市，需要確定一個訪問順序，使得每個城市都訪問一面，最后回到起點城市，且保證行走的總距離最短。
?????? 假設隨機生成10個城市坐標，城市之間的距離用歐式距離表示。

?????? 在以往例子中，粒子xi的值是具有意義的值，而這里需要優化序列的順序，因此粒子xi如何表示序列以及序列的改進是粒子群算法的一個關鍵所在。同時也應該認識到，雖然智能算法的原理比較容易理解，但在不同應用里，巧妙的構造解xi的表示形式，以及解xi的進化方式，才是關鍵。

?????? 在求解TSP問題時，可以用一個序列s表示城市的訪問順序，因此xi可以表示成如下形式：
?????? ?????? ?????? ?????? xi=(1,2,3,4,5)
?????? 式中表示有5個城市的TSP問題，其訪問順序為1-3-2-5-4-1.因此全部城市的所有可能序列構成了問題的搜索空間，例子位置的更新意味著粒子xi從所有城市的一種序列si變成另一種序列sj。
?????? 那么序列是如何變化得，這里引入交換子 和交換序列的概念。
??????交換子定義為s=Swap x(i,j)表示交換子s作用在序列x上使位置i的元素和位置j的元素相互交換，如x=(1,3,2,5,4)，對其添加一個交換操作s=Swap x(1,3)，其新的序列為x=(2,3,1,5,4)，在粒子群算法中，交換子可以看作是粒子的速度，能改變x的位置。
??????交換序列 定義為一組有前后順序的交換子的集合，即ss=[Swap1,Swap2…]，表示對序列x進行多次交換子交換操作，如x=(1,3,2,5,4)，交換序列ss=[(3,2),(1,5)]，當序列進行交換序列后，新的序列為x=(4,2,3,5,1)
??????解釋了交換子和交換序列后，重新定義粒子的速度和更新公式：

??????式中Pbesti-xi表示有一個交換序列ss，它使xi經過ss變換得到Pbest，同理Gbesti-xi也是一個交換序列ss，符號?表示兩個交換序列合并一個交換序列，如[(3,2),(1,5)]?[(4,2)]=[(3,2),(1,5),(4,2)]；符號+表示對一個序列按照交換序列ss執行交換操作。因此位置更新公式表示一個粒子經過交換序列Vi(t+1) 作用得到新的序列。r1,r2表示0-1之間的隨機數，表示一定概率進行交換子操作。

粒子群算法原理路見往期：粒子群算法求解無約束優化問題 源碼實現

算例

假設隨機生成10個城市坐標

import numpy as np import matplotlib.pyplot as plt import matplotlib as mpl import seaborn as snssns.set_style("whitegrid") mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默認字體 mpl.rcParams['axes.unicode_minus'] = False # 解決保存圖像是負號'-'顯示為方塊的問題# 固定隨機數種子 np.random.seed(1234)# 一些參數 city_num = 10 # 城市的數量 iter_num = 1000 # 算法最大迭代次數# 隨機生成city_num個城市的坐標，注意是不放回抽樣 X = np.random.choice(list(range(1, 100)), size=city_num, replace=False) Y = np.random.choice(list(range(1, 100)), size=city_num, replace=False)plt.scatter(X, Y, color='r') plt.title('城市坐標圖') plt.show()

步驟1：初始化參數

import numpy as np import matplotlib.pyplot as plt import matplotlib as mpl import seaborn as snssns.set_style("whitegrid") mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默認字體 mpl.rcParams['axes.unicode_minus'] = False # 解決保存圖像是負號'-'顯示為方塊的問題# 固定隨機數種子 np.random.seed(1234)# 一些參數 city_num = 10 # 城市的數量 size = 50 # 種群大小，即粒子的個數 r1 = 0.7 # pbest-xi 的保留概率 r2 = 0.8 # gbest-xi 的保留概率 iter_num = 1000 # 算法最大迭代次數 fitneess_value_list = [] # 每一步迭代的最優解# 隨機生成city_num個城市的坐標，注意是不放回抽樣 X = np.random.choice(list(range(1, 100)), size=city_num, replace=False) Y = np.random.choice(list(range(1, 100)), size=city_num, replace=False)

步驟2：距離函數 和適應度函數

# 計算城市之間的距離 def calculate_distance(X, Y):"""計算城市兩輛之間的歐式距離，結果用numpy矩陣存儲:param X: 城市的X坐標，np.array數組:param Y: 城市的Y坐標，np.array數組"""distance_matrix = np.zeros((city_num, city_num))for i in range(city_num):for j in range(city_num):if i == j:continuedis = np.sqrt((X[i] - X[j]) ** 2 + (Y[i] - Y[j]) ** 2) # 歐式距離計算distance_matrix[i][j] = disreturn distance_matrixdef fitness_func(distance_matrix, xi):"""適應度函數，計算目標函數值.:param distance: 城市的距離矩陣:param xi: PSO的一個解:return: 目標函數值，即總距離"""total_distance = 0for i in range(1, city_num):start = xi[i - 1]end = xi[i]total_distance += distance_matrix[start][end]total_distance += distance_matrix[end][xi[0]] # 從最后一個城市回到出發城市return total_distance

步驟3：繪制最優解的圖像函數

def plot_tsp(gbest):"""繪制最優解的圖形"""plt.scatter(X, Y, color='r')for i in range(1, city_num):start_x, start_y = X[gbest[i - 1]], Y[gbest[i - 1]]end_x, end_y = X[gbest[i]], Y[gbest[i]]plt.plot([start_x, end_x], [start_y, end_y], color='b', alpha=0.8)start_x, start_y = X[gbest[0]], Y[gbest[0]]plt.plot([start_x, end_x], [start_y, end_y], color='b', alpha=0.8)

步驟4：交換序列函數 和交換操作

import numpy as np
xi=np.array([1,2,3,4,5,6,7,8,9,10])
xbest=np.array([2,1,4,3,6,5,7,8,9,10])
velocity_ss = []
for i in range(len(xi)):
??????if xi[i] != xbest[i]:
??????j = np.where(xi == xbest[i])[0][0]
?????? print(j)
??????xi[i], xi[j] = xi[j], xi[i] #表示去掉重復計算，不同處只計算一處
# j=1,3,5

#交換序列函數 def get_ss(xbest, xi, r):#r=r1或r2"""計算交換序列，即 x2 經過交換序列ss得到x1，對應PSO速度更新公式的：r1(pbest-xi) 和 r2(gbest-xi):param xbest: pbest 或者 gbest:param xi: 例子當前解:return:xbest#最優解"""velocity_ss = []for i in range(len(xi)):if xi[i] != xbest[i]:#如果當前索引粒子元素不等于最優粒子元素j = np.where(xi == xbest[i])[0][0]#找到當前粒子和最優粒子不同的元素所在索引so = (i, j, r) # 得到交換子,表示以r的概率對i,j進行操作velocity_ss.append(so)xi[i], xi[j] = xi[j], xi[i] # 執行交換操作return velocity_ss#返回交換子序列#交換操作，返回序列xi解 def do_ss(xi, ss):"""執行交換操作:param xi::param ss: 由交換子組成的交換序列:return:"""for i, j, r in ss:rand = np.random.random()if rand <= r:xi[i], xi[j] = xi[j], xi[i]return xi

步驟5：PSO
??????第一步：計算城市之間的距離矩陣

# 計算城市之間的距離矩陣 distance_matrix = calculate_distance(X, Y)

??????第二步：初始化種群的各個粒子的位置，作為個體的歷史最優解pbest。用一個 50*10 的矩陣表示種群，每行表示一個粒子, 每一行是0-9不重復隨機數，表示城市訪問的順序

XX = np.zeros((size, city_num), dtype=np.int) for i in range(size):XX[i] = np.random.choice(list(range(city_num)), size=city_num, replace=False)

??????第三步：計算每個粒子對應適應度

pbest = XX#初始化粒子最優解 pbest_fitness = np.zeros((size, 1)) for i in range(size):pbest_fitness[i] = fitness_func(distance_matrix, xi=XX[i])

??????第四步：計算全局適應度和對應的解gbest

gbest = XX[pbest_fitness.argmin()] gbest_fitness = pbest_fitness.min()

??????第五步：迭代尋優

# 記錄算法迭代效果 fitneess_value_list.append(gbest_fitness)# 下面開始迭代 for i in range(iter_num):#遍歷迭代次數for j in range(size): # 對每個粒子迭代pbesti = pbest[j].copy() # 此處要用copy，不然出現淺拷貝問題xi = XX[j].copy()# 計算交換序列，即 v = r1(pbest-xi) + r2(gbest-xi)ss1 = get_ss(pbesti, xi, r1)ss2 = get_ss(gbest, xi, r2)ss = ss1 + ss2# 執行交換操作，即 x = x + vxi = do_ss(xi, ss)# 判斷是否更優fitness_new = fitness_func(distance_matrix, xi)fitness_old = pbest_fitness[j]if fitness_new < fitness_old:pbest_fitness[j] = fitness_newpbest[j] = xi# 判斷全局是否更優# 記錄算法迭代效果# fitneess_value_list.append(gbest_fitness)gbest_fitness_new = pbest_fitness.min()gbest_new = pbest[pbest_fitness.argmin()]if gbest_fitness_new < gbest_fitness:gbest_fitness = gbest_fitness_newgbest = gbest_new# 加入到列表fitneess_value_list.append(gbest_fitness)

??????第六步：繪制結果

# 打印迭代的結果 print('迭代最優結果是：', gbest_fitness) print('迭代最優變量是：', gbest) # 迭代最優結果是： 230.344 # 迭代最優變量是： [5 8 2 3 6 1 7 0 4 9]# 繪制TSP路徑圖 plot_tsp(gbest) plt.title('TSP路徑規劃結果') plt.show()

往期TSP回顧
基于組合遺傳粒子群算法的旅行商問題求解
遺傳算法求最短路徑（旅行商問題）python實現

作者：電氣余登武

總結

以上是生活随笔為你收集整理的粒子群算法求解旅行商问题的全部內容，希望文章能夠幫你解決所遇到的問題。

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