模拟退火算法解决TSP问题

计算智能第二次大作业

用模拟退火算法(SA)解决旅行商问题(TSP)。第二次用python写的作业,挺有收获的。

路径调整有坑,调了好久……

问题描述

TSP问题(旅行商问题)是指旅行家要旅行n个城市,要求各个城市经历且仅经历一次然后回到出发城市,并要求所走的路程最短。

转换为图论问题即:给定一个 $n$ 个节点的加权完全图 $d$ ,求一权重最小的哈密尔顿回路 $p$

动态规划求解的经典问题。但是不用动态规划了,用模拟退火来实现。

要点

参数

  1. 符号说明
    • 总共 $n$ 个节点,节点用数字编号 $0,1,\cdots,n$ 表示
    • 任意两个节点距离用邻接矩阵 $ d[i,j] $ ,其中 $i,j\in {0,1,\cdots,n} $
    • 访问路径$ path[0:n-1] $,第 $i$ 次访问的城市标号记录在 $path[i]$ 中
  2. 冷却表
    • 控制参数 $t$ 初始值设置为280
    • Mapkob链长即每次搜索新解数为$100n$
    • 衰减函数 $\alpha$ 设置为 0.92
    • 停止条件为连续两次Mapkob链路径无变化

细节思路

  • 目标函数为解的路径长度length。

  • 每次生成M链新解采用随机交换当前解路径上的两个节点Generate()。

  • 新旧两个状态的目标函数的差df=f(new)-f(now)的时候不需要遍历一遍路径,路径值变动只和 当前/新 路径上节点左右两边的边有关,仅需加上新连接的边长,减去旧的边长即可。注意分两个交换的节点是否相邻/非相邻,相邻不是非相邻的特例,两者要分开写,相邻节点的交换仅影响4条边,非相邻影响8条边。如图。并且还要注意B节点有可能是最后一个节点,使用B+1节点需mod node

    image-20181025165242420

image-20181025165359465

  • 若满足转移条件,记得在解路径path上交换两个节点。并更新路径长度length +=df

python代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Oct 26 00:03:44 2018

@author: jiang
"""

import numpy as np
import pylab as pl
import matplotlib as plt

def initialize(node,L,seed = None):
#随机初始化生成完全图 d 路径path
#节点数 node 初始温度t0 Mapkob链长L
#length 当前路径长度
# city from 0 to node-1
np.random.seed(seed)
d = np.random.randint(1,1000, (node, node))
d = np.triu(d,1) + np.triu(d,1).T #生成完全图
#d = np.array([[0,3,1,5,8],[3,0,6,7,9],[1,6,0,4,2],[5,7,4,0,3],[8,9,2,3,0]])
path = [i for i in range(0,node)] #路径
length = np.sum(d.diagonal(offset = 1)) + d[node-1][0] #路径长度
t = 0.28 * np.sum(d)//(node*node-node)
print(np.sum(d))
return d,t,path,length

def Generate(node):
#产生交换的两个节点
a,b = 0,0
np.random.seed()
while a == b:
a = np.random.randint(1,node)
b = np.random.randint(1,node)
if a>b: #顺序
a,b = b,a
return a,b

if __name__=='__main__':

node = 5 #节点数
L = node *100 #冷却表参数
alpha = 0.92 #冷却表参数
d,t,path,length = initialize(node,L,seed = 1) #初始化
print('d=\n',d,'\n','t0=',t)
tim = 0 #计数器
s = 0
log = []
bChange = 1
while bChange == 1 :
bChange = 0
#print(t,length)
for i in range(1,L):
a,b = Generate(node) #生成随机转移位置a,b

#生成目标函数差 df = f(j)-f(i)
if b-a == 1:
df = d[path[b]][path[a-1]] + d[path[a]][path[(b+1)%(node)]] \
- d[path[a]][path[a-1]] - d[path[b]][path[(b+1)%(node)]]
else:
df = d[path[a]][path[(b+1)%(node)]] + d[path[b]][path[(a+1)]] + d[path[a]][path[(b-1)]] + d[path[b]][path[(a-1)]] - d[path[b]][path[(b-1)]] - d[path[b]][path[(b+1)%node]] - d[path[a]][path[(a-1)]] - d[path[a]][path[(a+1)]]

#判断转移条件
if df<0 or np.exp(-df/t)>np.random.rand():
length += df
path[a],path[b] = path[b],path[a]
bChange = 1
log.append(length)
tim = tim+1
t = t * alpha
if bChange == 0:
s = s + 1
else:
s = 0
print('Path:',path)
print('Length:',length)
print('Times:',tim)
print('Final temp:',t)
1
2
3
4
5
6
7
8
9
10
11
12
13
#输出图像函数
%matplotlib inline
plt.rcParams['figure.figsize'] = (5.0,5.0)
plt.rcParams['savefig.dpi'] = 100
plt.rcParams['figure.dpi'] = 100

pl.grid()
pl.xlabel('Time')
pl.ylabel('Length')
pl.plot(log[:],'-')
# ax = pl.gca()
# ax.xaxis.get_major_formatter().set_powerlimits((0,1))
pl.show()

测试: 500节点 边权重为rand(1,1000) 随机种子seed = 1

结果:image-20181027141840754

image-20181027142919939