SA求解TSP

SA(模拟退火算法)求解TSP问题

步骤

读取文件生产坐标矩阵和距离矩阵

Step1初始化

Step2定义当前状态的状态邻域转移

1 随机产生两点a、b（a<b），交换状态序列中a、b位置的值。

2 随机产生三点a、b、c，将当前状态序列中[a:b]的值移到位置c的后面。

3 随机产生两个点a、b，将当前状态序列[a:b]内的子序列全部反转。

（前两个对于求解15城市的TSP已足够）

Step3计算cost（i），与当前cost比较，并依概率更新当前cost值。

Step4 温度下降、保存当前最优解

求解结果

状态序列运行结果图

状态序列图以及cost损失函数图

算法思想

利用模拟退火算法的思想求解TSP问题，将TSP不同的解（回路）看作不同的状态，其cost值是各个城市总的距离值之和，求解TSP问题其实就是在N！的空间状态下找到cost值最小的那个状态，但是由于N！极为庞大，若进行遍历搜索会导致NP难问题，所以我们采用模拟退火的算法来求解。首先初始化初始状态、温度t、温度下降系数α、马尔科夫链长度、初始cost等，然后在初始状态的邻域里面寻找n个不同的状态（n为马尔科夫链的长度）并分别计算不同状态的cost（i），若cost>cost（i），表明当前状态距离更短，将解更新为当前状态；如果cost<cost（i），则依概率（P=exp⁡(-(cost(i)-cost)/t)）更新解为当前状态；直至（状态数>n）马尔科夫链结束，根据（P=exp⁡(-(cost(i)-cost)/t)）更新下降温度值，直至低于阈值，或长时间cost不变；结束，输出解。

代码：

import numpy as np
import matplotlib.pyplot as plt

# 处理数据
coord = []
with open("p_xy.txt", "r") as lines:
    lines = lines.readlines()
    for line in lines:
        xy = line.split()
        coord.append(xy)

coord = np.array(coord)
w, h = coord.shape
coordinates = np.zeros((w, h), float)
for i in range(w):
    for j in range(h):
        coordinates[i, j] = float(coord[i, j])
# print(coordinates)

# 得到距离矩阵
distance = np.zeros((w, w))
for i in range(w):
    for j in range(w):
        distance[i, j] = distance[j, i] = np.linalg.norm(coordinates[i] - coordinates[j])

# SA

# 初始化
alpha = 0.99
markovlen = 1000
solution_new = np.arange(w)
solution_cur = solution_new.copy()
value_cur = 100000
solution_bst = solution_new.copy()
value_bst = 100000
t_arange = (1, 100)
t = t_arange[1]
result = []

while t > t_arange[0]:
    for i in range(markovlen):
        # 在邻域内随机产生新解
        # 1 位置交换
        if np.random.rand() > 0.5:
            while True:
                loc1 = np.int(np.ceil(np.random.rand() * (w-1)))
                loc2 = np.int(np.ceil(np.random.rand() * (w-1)))
                if loc1 != loc2:
                    break
            solution_new[loc1], solution_new[loc2] = solution_new[loc2], solution_new[loc1]
        # 2 三位置插入 排序 并将[loc1:loc2] 插入到loc3 后面
        else:
            while True:
                loc1 = np.int(np.ceil(np.random.rand() * (w-1)))
                loc2 = np.int(np.ceil(np.random.rand() * (w-1)))
                loc3 = np.int(np.ceil(np.random.rand() * (w-1)))
                if (loc1 != loc2) and (loc2 != loc3) and (loc3 != loc1):
                    break

            if loc1 > loc2:
                loc2, loc1 = loc1, loc2
            if loc2 > loc3:
                loc3, loc2 = loc2, loc3
            if loc1 > loc2:
                loc1, loc2 = loc2, loc1

            temp1 = solution_new[loc1:loc2 + 1].copy()
            temp2 = solution_new[loc2 + 1:loc3 + 1].copy()
            solution_new[loc1:loc1 + (loc3 - loc2)] = temp2.copy()
            solution_new[loc1 + (loc3 - loc2):loc3 + 1] = temp1.copy()
        # # 3 区间内反转
        # if np.random.rand() > 0.5:
        #     while True:
        #         loc1 = np.int(np.ceil(np.random.rand() * w))
        #         loc2 = np.int(np.ceil(np.random.rand() * w))
        #         if loc1 != loc2:
        #             break

        # 计算cost 别忘记回到原点
        value_new = 0
        for k in range(w - 1):
            value_new += distance[solution_new[k], solution_new[k + 1]]
        value_new += distance[solution_new[w - 1], solution_new[0]]

        print("T:{}, markov_node:{} cost:{}---->".format(t, i, value_new))

        if value_new < value_cur:
            # print("---Accept---")
            value_cur = value_new
            solution_cur = solution_new.copy()

            if value_new < value_bst:
                # print("---Best---")
                value_bst = value_new
                solution_bst = solution_new.copy()

        elif np.random.rand() < np.exp(-(value_new - value_cur) / t):
            # print("---概率接受---")
            value_cur = value_new
            solution_cur = solution_new.copy()

        else:
            # print("---拒绝接受---")
            value_new = value_cur
            solution_new = solution_cur.copy()

    # 温度下降
    t = t * alpha
    result.append(value_bst)
    print("---Best_Plan_Cost_Now:{}---".format(value_bst))

for i in range(w):
    print(solution_bst[i] + 1, end="--->")
print(solution_bst[0] + 1)
print("Cost:{}".format(value_bst))

# 画图
cor_cur = np.array((w+1, h))
for i in range(w):
    cor_cur[i, :] = coordinates[solution_bst[i], :].copy()


plt.figure(figsize=(20, 20))

plt.subplot(1, 2, 1)
plt.plot(cor_cur[:, 0], cor_cur[:, 1])

plt.subplot(1, 2, 2)
plt.plot(np.array(result))
plt.ylabel("BestValue")
plt.xlabel("t({}->{})".format(t_arange[1], t_arange[0]))
plt.show()