模拟退火算法（解决多元和tsp问题

模拟退火算法

简介

模拟退火算法得益于材料统计力学的研究成果。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下，粒子的能量较高，可以自由运动和重新排列。在低温条件下，粒子能量较低。如果从高温开始，非常缓慢地降温（这个过程被称为退火），粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于低能状态的晶体。

算法流程（以求最小值为例）

冷却进度表

退火过程由一组初始参数，即冷却进度表控制，它的目的是尽量使系统达到平衡，
以使算法在有限的时间内逼近最优解

参数的问题解决了，那我们怎么根据当前解产生新解？这个是要根据具体的问题去实现的

python代码实现，TSP问题

已知100个目标的经度、纬度如表17.1所示。我方有一个基地，经度和纬度为（70,40）。假设我方飞机的速度为1000公里/小时。我方派一架飞机从基地出发，侦察完所有目标，再返回原来的基地。在每一目标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。

这是一个旅行商问题。给我方基地编号为0，目标依次编号为1，2，…，100，最后我方基地再重复编号为101（这样便于程序中计算）。距离矩阵 ，其中 表示表示 两点的距离， ，这里 为实对称矩阵。则问题是求一个从点0出发，走遍所有中间点，到达点101的一个最短路径。

python代码

from numpy import loadtxt,radians,sin,cos,inf,exp
from numpy import array,r_,c_,arange,savetxt
from numpy.lib.scimath import arccos
from numpy.random import shuffle,randint,rand
from matplotlib.pyplot import plot, show, rc
a=loadtxt("Pdata17_1.txt")
x=a[:,::2]. flatten(); y=a[:,1::2]. flatten()
d1=array([[70,40]]); xy=c_[x,y]
xy=r_[d1,xy,d1]; N=xy.shape[0]
t=radians(xy)  #转化为弧度
d=array([[6370*arccos(cos(t[i,0]-t[j,0])*cos(t[i,1])*cos(t[j,1])+ sin(t[i,1])*sin(t[j,1])) for i in range(N)]
   for j in range(N)]).real
savetxt('Pdata17_2.txt',c_[xy,d]) #前两列是102*2的矩阵代表102个点的位置坐标，后面是102*102的矩阵，代表的是距离矩阵
path=arange(N); L=inf
for j in range(1000):   #蒙特卡洛
	path0=arange(1,N-1); shuffle(path0)
	path0=r_[0,path0,N-1]; L0=d[0,path0[1]]  #初始化
	for i in range(1,N-1):L0+=d[path0[i],path0[i+1]]
	if L0<L: path=path0; L=L0
print(path,'\n',L)  

e=0.1**30; M=20000; at=0.999; T=1   #模拟退火
for k in range(M):
	c=randint(1,101,2); c.sort()
	c1=c[0]; c2=c[1]
	df=d[path[c1-1],path[c2]]+d[path[c1],path[c2+1]]-d[path[c1-1],path[c1]]-d[path[c2],path[c2+1]]
  
	if df<0:
    	path=r_[path[0],path[1:c1],path[c2:c1-1:-1],path[c2+1:102]]; L=L+df
	else:
    	if exp(-df/T)>=rand(1):
        	path=r_[path[0],path[1:c1],path[c2:c1-1:-1],path[c2+1:102]]
        	L=L+df
	T=T*at
	if T<e: break54 
print(path,'\n',L)  #输出巡航路径及路径长度
xx=xy[path,0]; yy=xy[path,1]; rc('font',size=16)
plot(xx,yy,'-*'); show()  #画巡航路径