算法原理

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

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

?????? 在求解TSP問題時，可以用一個序列s表示城市的訪問順序，因此xi可以表示成如下形式：
?????? ?????? ?????? ?????? xi=(1,2,3,4,5)
?????? 式中表示有5個城市的TSP問題，其訪問順序為1-3-2-5-4-1.因此全部城市的所有可能序列構(gòu)成了問題的搜索空間，例子位置的更新意味著粒子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)]，當(dāng)序列進行交換序列后，新的序列為x=(4,2,3,5,1)
??????解釋了交換子和交換序列后，重新定義粒子的速度和更新公式：

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

粒子群算法原理路見往期：粒子群算法求解無約束優(yōu)化問題 源碼實現(xiàn)

算例

假設(shè)隨機生成10個城市坐標(biāo)

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'] # 指定默認(rèn)字體 mpl.rcParams['axes.unicode_minus'] = False # 解決保存圖像是負(fù)號'-'顯示為方塊的問題# 固定隨機數(shù)種子 np.random.seed(1234)# 一些參數(shù) city_num = 10 # 城市的數(shù)量 iter_num = 1000 # 算法最大迭代次數(shù)# 隨機生成city_num個城市的坐標(biāo)，注意是不放回抽樣 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('城市坐標(biāo)圖') plt.show()

步驟1：初始化參數(shù)

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'] # 指定默認(rèn)字體 mpl.rcParams['axes.unicode_minus'] = False # 解決保存圖像是負(fù)號'-'顯示為方塊的問題# 固定隨機數(shù)種子 np.random.seed(1234)# 一些參數(shù) city_num = 10 # 城市的數(shù)量 size = 50 # 種群大小，即粒子的個數(shù) r1 = 0.7 # pbest-xi 的保留概率 r2 = 0.8 # gbest-xi 的保留概率 iter_num = 1000 # 算法最大迭代次數(shù) fitneess_value_list = [] # 每一步迭代的最優(yōu)解# 隨機生成city_num個城市的坐標(biāo)，注意是不放回抽樣 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：距離函數(shù) 和適應(yīng)度函數(shù)

# 計算城市之間的距離 def calculate_distance(X, Y):"""計算城市兩輛之間的歐式距離，結(jié)果用numpy矩陣存儲:param X: 城市的X坐標(biāo)，np.array數(shù)組:param Y: 城市的Y坐標(biāo)，np.array數(shù)組"""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):"""適應(yīng)度函數(shù)，計算目標(biāo)函數(shù)值.:param distance: 城市的距離矩陣:param xi: PSO的一個解:return: 目標(biāo)函數(shù)值，即總距離"""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]] # 從最后一個城市回到出發(fā)城市return total_distance

步驟3：繪制最優(yōu)解的圖像函數(shù)

def plot_tsp(gbest):"""繪制最優(yōu)解的圖形"""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：交換序列函數(shù) 和交換操作

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] #表示去掉重復(fù)計算，不同處只計算一處
# j=1,3,5

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

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

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)

??????第三步：計算每個粒子對應(yīng)適應(yīng)度

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

??????第四步：計算全局適應(yīng)度和對應(yīng)的解gbest

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

??????第五步：迭代尋優(yōu)

# 記錄算法迭代效果 fitneess_value_list.append(gbest_fitness)# 下面開始迭代 for i in range(iter_num):#遍歷迭代次數(shù)for j in range(size): # 對每個粒子迭代pbesti = pbest[j].copy() # 此處要用copy，不然出現(xiàn)淺拷貝問題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# 執(zhí)行交換操作，即 x = x + vxi = do_ss(xi, ss)# 判斷是否更優(yōu)fitness_new = fitness_func(distance_matrix, xi)fitness_old = pbest_fitness[j]if fitness_new < fitness_old:pbest_fitness[j] = fitness_newpbest[j] = xi# 判斷全局是否更優(yōu)# 記錄算法迭代效果# 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)

??????第六步：繪制結(jié)果

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

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

作者：電氣余登武

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。