python模拟退火求解TSP问题

问题描述

利用模拟退火算法。求解30个城市的TSP问题。两城市之间距离用直角坐标系中的两点距离公式。
{41,94},{37,84},{54,67},{25,62},{7,64},{2,99},{68,58},{71,44},{54,62},{83,69}
{64,60},{18,54},{22,60},{83,46},{91,38},{25,38},{24,42},{58,69},{71,71},{74,78}
{87,76},{18,40},{13,40},{82,7},{62,32},{58,35},{45,21},{41,26},{44,35},{4,50}

思路

模拟退火简介

  • (1)随机挑选一个单元k,并给它一个随机的位移,求出系统因此而产生的能量变化ΔEk。
  • (2)若ΔEk⩽0,该位移可采纳,而变化后的系统状态可作为下次变化的起点;
    若ΔEk>0,位移后的状态可采纳的概率为
    python模拟退火求解TSP问题
    式中T为温度,然后从(0,1)区间均匀分布的随机数中挑选一个数RR,若R<Pk,则将变化后的状态作为下次的起点;否则,将变化前的状态作为下次的起点。
  • (3)转第(1)步继续执行,知道达到平衡状态为止。

TSP问题

模拟退火的关键是在一定概率P下接受比当前差的状态,而概率P又与当前状态有关,当前状态越差,接受概率P越大,当前状态越好,接受概率P越小。
TSP问题,

  1. 初始化,模拟退火中用温度来观察程序进程,这里模拟温度,设置t0=100,从100°开始降温,程序停止条件t0==1,温度下降的速率为alpha=0.9,在每个温度状态下,程序要尝试随机生成路径TIME次,TIME=1000
  2. 创建城市坐标矩阵pointspoint[i][0]表示第i个城市的x坐标,point[i][1]表示第i个城市的y坐标
  3. 先随机生成一条路径way,通过函数calWayDis(way)计算路径总长度waydis
  4. 再随机选择2个城市city1city2,交换way中这两个城市的位置,形成新的路径tempway
tempway[city1],tempway[city2]=tempway[city2],tempway[city1]
  1. calWayDis(tempway)计算tempway的路径长度tempdis
  2. 对比:
    1.如果:
    新的路径长度更小,更新路径,
     并与当前最小路径对比
       如果小于当前最小路径,更新最小路径与最小值
if tempdis<waydis:
	way = tempway.copy()
	waydis = tempdis
if tempdis<bestdis:
	bestway = tempway.copy()
	bestdis = tempdis
	draw(bestway,bestdis)

2. 如果:
新的路径长度不小于现有路径长度,则根据模拟退火的概率,去接受

if tempdis>waydis:
if np.random.rand()<np.exp(-(tempdis-waydis)/(t0)):
way = tempway.copy()
waydis = tempdis
  1. 步骤4、5、6在某个t0状态,迭代TIME次,每次way会交换两个城市生成tempway,然后t0=t0*alph温度状态改变重复步骤4、5、6,直到t0=1,得到目前最优的路径
while t0>t[0]:
	步骤4
	步骤5
	步骤6
	t0 *= alpha
  1. 打印出得到的最优解,路径和路径长度

代码

import numpy as np
import matplotlib.pyplot as plt 
import pdb
import imageio
import shutil
import os 

# 生成距离矩阵
def getDismat(points):
	dismat = np.zeros((N,N))
	for i in range(N):
		for j in range(i,N):
			dismat[i][j] = dismat[j][i] = np.linalg.norm(points[i]-points[j])
	return dismat	

def init():
	alpha = 0.9
	t = (1,100)
	TIME = 1000
	way = np.arange(N)
	waydis = calWayDis(way)
	return alpha,t,TIME,way,waydis

# 计算路径长度
def calWayDis(way0):
	waydis = 0
	for i in range(N-1):
		waydis +=dismat[way0[i]][way0[i+1]] 
	waydis += dismat[way0[N-1]][way0[0]]
	return waydis

def draw(way,dist):
	global N,points ,TIMESIT, PNGFILE, PNGLIST
	plt.cla()
	plt.title('cross=%.4f' % dist)
	xs = [points[i][0] for i in range(N)]
	ys = [points[i][1] for i in range(N)]
	plt.scatter(xs, ys, color='b')
	xs = np.array(xs)
	ys = np.array(ys)
	# plt.plot(xs[[0, solutionpath[0]]], ys[[0, solutionpath[0]]], color='r')
	# 连接路径
	for i in range(N-1):
		plt.plot(xs[[way[i], way[i + 1]]], ys[[way[i], way[i + 1]]], color='r')
	# 将终点与起点连接
	plt.plot(xs[[way[N - 1], 0]], ys[[way[N - 1], 0]], color='r')
	plt.scatter(xs[0], ys[0], color='k', linewidth=10)
	for i, p in enumerate(points):
		plt.text(*p, '%d' % i)
	plt.savefig('%s/%d.png' % (PNGFILE, TIMESIT))
	PNGLIST.append('%s/%d.png' % (PNGFILE, TIMESIT))
	TIMESIT += 1

points = np.array([[41,94],[37,84],[54,67],[25,62],[7,64],[2,99],[68,58],[71,44],[54,62],[83,69],
[64,60],[18,54],[22,60],[83,46],[91,38],[25,38],[24,42],[58,69],[71,71],[74,78],
[87,76],[18,40],[13,40],[82,7],[62,32],[58,35],[45,21],[41,26],[44,35],[4,50]
])
N = points.shape[0]
dismat = getDismat(points)
alpha,t,TIME,way,waydis=init()
t0 = t[1]
K=0.8

# 画图初试化
TIMESIT = 0
PNGFILE = './png/'
PNGLIST = []
if not os.path.exists(PNGFILE):
    os.mkdir(PNGFILE)
else:
    shutil.rmtree(PNGFILE)
    os.mkdir(PNGFILE)
# 记录每次迭代的结果
result = []
tempway = way.copy()
bestway = way.copy()
bestdis = 10000

while t0>t[0]:
	for i in range(TIME):
		if np.random.rand() > 0.5:
		# 两点交换
			while True:
			# 随机生成不同2个点,
				city1 = np.int(np.ceil(np.random.rand()*(N-1)))
				city2 = np.int(np.ceil(np.random.rand()*(N-1)))
				if city1!=city2:
					break
			# 交换
	tempway[city1],tempway[city2]=tempway[city2],tempway[city1]
		else:
			# 3个点
			while True:
				city1 = np.int(np.ceil(np.random.rand()*(N-1)))
				city2 = np.int(np.ceil(np.random.rand()*(N-1))) 
				city3 = np.int(np.ceil(np.random.rand()*(N-1)))
				if((city1 != city2)&(city2 != city3)&(city1 != city3)):
					break
			# 下面的三个判断语句使得city1<city2<city3
			if city1 > city2:
				city1,city2 = city2,city1
			if city2 > city3:
				city2,city3 = city3,city2
			if city1 > city2:
				city1,city2 = city2,city1
			#下面的三行代码将[city1,city2)区间的数据插入到city3之后
			temp = tempway[city1:city2].copy()
			tempway[city1:city3-city2+1+city1] = tempway[city2:city3+1].copy()
			tempway[city3-city2+1+city1:city3+1] = temp.copy()

		tempdis = calWayDis(tempway)
		if tempdis<waydis:
			way = tempway.copy()
			waydis = tempdis
		if tempdis<bestdis:
			bestway = tempway.copy()
			bestdis = tempdis
			draw(bestway,bestdis)
		else:
			if np.random.rand()<np.exp(-(tempdis-waydis)/(t0)):
				way = tempway.copy()
				waydis = tempdis
				# 更新路径
			else: tempway = way.copy()

	t0 *= alpha
	result.append(bestdis)

print("最短路径:%s"%np.array(result[-1]))
for i in bestway:
	print(i,end="-->")
print(bestway[0])

#生成路径gif 
generated_images = []
for png_path in PNGLIST:
    generated_images.append(imageio.imread(png_path))
# shutil.rmtree(PNGFILE)  # 可删掉
generated_images = generated_images + [generated_images[-1]] * 10
imageio.mimsave('Sovle_TSP01.gif', generated_images, 'GIF', duration=0.1)

运行结果

python模拟退火求解TSP问题python模拟退火求解TSP问题