TSP问题——模拟退火法

1.前言

TSP问题，又称旅行商问题，是数学领域一个著名的问题。假设有一个旅行商人要拜访n个城市，每个城市只访问一次，最后要回到原来出发的城市，这样会有很多路径，选择其中路径长度最小的路径。该问题是一个np完全问题，因为该问题实质是在一个带权完全无向图中，找一个权值最小的Hamilton回路，该问题的可行解是所有顶点的全排列，随着顶点数的增加，会产生组合爆炸。早期的研究者使用精确算法求解该问题，常用的方法包括：分枝定界法、线性规划法、动态规划法等。但是，随着问题规模的增大，精确算法将变得无能为力，因此，在后来的研究中，国内外学者重点使用近似算法或启发式算法，主要有遗传算法、模拟退火法、蚁群算法、禁忌搜索算法、贪婪算法和神经网络等。【百度百科】

2.模拟退火法流程

2.1 爬山法

先来看看爬山法。爬山法，又称瞎子爬山法，这种方法只利用到了局部的信息反馈，像瞎子一样只能利用拐棍来探测周围(局部)的高低情况，试探着爬山。以TSP问题为例，爬山法的基本步骤如下：

1. 考虑问题的编码问题，如TSP问题中可以给每个城市一个编号，n个城市的编号从0到n-1。那么所有的解都可以用一个一维数组表示，包含0到n-1的所有数字且不重复。
2. 生成一个可能的解。若达到目标则停止，否则继续。生成一个可能的解。若达到目标(一般对于TSP问题这种np完全问题，我们只要求达到接近最优解的10%即可)则停止，否则继续。
3. 从当前解出发，生成新的解。TSP问题中生成新解的方法有很多，这里简要介绍三种：
   
   如果使用多种生成新解的方法，则可以将这个爬山法称为多领域操作爬山法，与传统的爬山法有很大的差别，稍后介绍一下。
4. 如果上述步骤产生的新解中有最优解，停止。否则找出最接近最优解的最佳解(TSP问题中的最佳解可以是路径长度最小的解)，从该解出发，重新进行步骤2、步骤3、步骤4，直到满意为止。当然，不会一直这样迭代下去，因为个人感到的所谓的满意是不同的，只要在进行爬山法的时候控制步骤2、3、4的迭代次数就行了，还可以在找到满意解或者最佳解长期没有变化后跳出循环。

爬山法有很多问题，由于它只能探测周围的情况来生成新解，所以它很容易陷入局部最优，而且无法处理平顶情况。

2.2 模拟退火法

模拟退火法是克服爬山法缺点的有效方法。这是模仿工业界中为了达到某些特种晶体结构重复将金属不断加热或冷却的过程。模拟退火法能在多项式时间里给出一个近似最优解，而且不会陷入局部最优。模拟退火法的基本思想是在系统朝能量减少的方向减少的时候，偶尔控制系统跳到能量较高的状态，以避免局部极小点，最终稳定在全局极小点。这样讲可能有点抽象，那我们用参数来细化一下吧：

1. 先定义三个常数，作为初温T_initial、末温T_final、温度变化控制参数C。
2. 像爬山法一样产生一个初始解，判断是否满意。
3. 将T_initial不断的乘上C来让T_initial接近T_final。在T_initial变化的每次迭代中可以从当前解产生新解，使用一定次数的循环来探测当前解在解空间中的周围情况，如果有比当前解小的接受它作为当前解，如果比它大的则以一定的概率接受它，以跳出局部最优。接受差解的概率应该随着温度降低而不断减少，越接近最优解，越没必要接受差解。假设当前解与产生的新解的差为diff，这个概率的计算如下：
   p = 1 / (1 + e^-diff/T)，其中T是当前温度。
   最后我们可以定一个条件在得到最够好的解后直接退出迭代，不用进行没用的迭代了。
4. 对了，一般我们都会维护一个全局最优解记录所有迭代所有产生的新解的最佳解，因为局部的最佳解可能会远离最优解。

2.3 比较一下

首先是爬山法，这个方法的缺点已经介绍了，就是局部最优以及平顶问题。而模拟退火法能够解决这个问题，因为模拟退火法能接受差解。当然，多邻域操作爬山法甚至就是简单的单领域爬山法其实也是能从一些局部最优的“陷阱”中出来的，只要它们产生新解的邻域操作有足够的扰动能力(探测周围的能力，跳出这个局部来探测到更好的情况)。

3.程序流程

3.1 环境以及数据

编程环境：python3、pycharm
数据来源：https://www.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/tsp/ 的d198.tsp，即有198个城市。
(其实直接搜tsplib就可以了)
当然，为了操作方便，我将数据形式修改了一下，如下：

198
0.00000e+00 0.00000e+00
5.51200e+02 9.96400e+02
6.27400e+02 9.96400e+02
7.03600e+02 9.96400e+02
7.03600e+02 1.04720e+03
6.27400e+02 1.04720e+03
5.51200e+02 1.04720e+03
8.81400e+02 1.35200e+03
9.32200e+02 1.35200e+03
9.57600e+02 1.35200e+03
9.83000e+02 1.35200e+03
1.00840e+03 1.35200e+03
1.03380e+03 1.35200e+03
1.31320e+03 1.12340e+03
1.28780e+03 1.09800e+03
1.28780e+03 9.96400e+02
1.31320e+03 9.96400e+02
1.46560e+03 9.96400e+02
1.51640e+03 9.96400e+02
1.59260e+03 9.96400e+02
1.59260e+03 1.09800e+03
1.51640e+03 1.09800e+03
1.46560e+03 1.09800e+03
1.56720e+03 1.12340e+03
1.59260e+03 1.14880e+03
1.56720e+03 1.17420e+03
1.54180e+03 1.17420e+03
1.49100e+03 1.17420e+03
1.44020e+03 1.17420e+03
1.46560e+03 1.19960e+03
1.41480e+03 1.22500e+03
1.44020e+03 1.22500e+03
1.49100e+03 1.22500e+03
1.51640e+03 1.22500e+03
1.59260e+03 1.25040e+03
1.46560e+03 1.25040e+03
1.44020e+03 1.25040e+03
1.38940e+03 1.25040e+03
1.54180e+03 1.27580e+03
1.16080e+03 1.12340e+03
1.16080e+03 1.22500e+03
1.26240e+03 1.30120e+03
1.28780e+03 1.30120e+03
1.33860e+03 1.30120e+03
1.41480e+03 1.30120e+03
1.49100e+03 1.30120e+03
1.54180e+03 1.30120e+03
1.64340e+03 1.30120e+03
1.66880e+03 1.32660e+03
1.61800e+03 1.32660e+03
1.56720e+03 1.32660e+03
1.51640e+03 1.32660e+03
1.46560e+03 1.32660e+03
1.41480e+03 1.32660e+03
1.33860e+03 1.32660e+03
1.31320e+03 1.32660e+03
1.23700e+03 1.32660e+03
1.23700e+03 1.35200e+03
1.31320e+03 1.35200e+03
1.33860e+03 1.35200e+03
1.41480e+03 1.35200e+03
1.46560e+03 1.35200e+03
1.69420e+03 1.35200e+03
1.61800e+03 1.37740e+03
1.51640e+03 1.37740e+03
1.41480e+03 1.37740e+03
1.33860e+03 1.37740e+03
1.31320e+03 1.37740e+03
1.23700e+03 1.37740e+03
1.35770e+03 1.39010e+03
1.23700e+03 1.40280e+03
1.31320e+03 1.40280e+03
1.33860e+03 1.40280e+03
1.41480e+03 1.40280e+03
1.46560e+03 1.40280e+03
1.56720e+03 1.40280e+03
1.59260e+03 1.40280e+03
1.61800e+03 1.40280e+03
1.69420e+03 1.42820e+03
1.66880e+03 1.42820e+03
1.54180e+03 1.42820e+03
1.44020e+03 1.42820e+03
1.41480e+03 1.42820e+03
1.33860e+03 1.42820e+03
1.26240e+03 1.42820e+03
1.23700e+03 1.45360e+03
1.33860e+03 1.45360e+03
1.41480e+03 1.45360e+03
1.46560e+03 1.45360e+03
1.49100e+03 1.45360e+03
1.66880e+03 1.45360e+03
1.69420e+03 1.45360e+03
1.61800e+03 1.47900e+03
1.59260e+03 1.47900e+03
1.56720e+03 1.47900e+03
1.51640e+03 1.47900e+03
1.44020e+03 1.47900e+03
1.41480e+03 1.47900e+03
1.33860e+03 1.47900e+03
1.26240e+03 1.47900e+03
1.59260e+03 1.50440e+03
1.61800e+03 1.50440e+03
1.66880e+03 1.52980e+03
1.69420e+03 1.52980e+03
1.74500e+03 1.52980e+03
1.82120e+03 1.52980e+03
1.84660e+03 1.52980e+03
1.94820e+03 1.52980e+03
1.92280e+03 1.55520e+03
1.89740e+03 1.55520e+03
1.87200e+03 1.55520e+03
1.82120e+03 1.55520e+03
1.79580e+03 1.55520e+03
1.77040e+03 1.55520e+03
1.77040e+03 1.65680e+03
1.79580e+03 1.65680e+03
1.82120e+03 1.65680e+03
1.87200e+03 1.65680e+03
1.89740e+03 1.65680e+03
1.92280e+03 1.65680e+03
1.80850e+03 1.69490e+03
1.75770e+03 1.69490e+03
1.88470e+03 1.73300e+03
1.99900e+03 1.73300e+03
2.07520e+03 1.73300e+03
2.11330e+03 1.73300e+03
2.17680e+03 1.73300e+03
2.23650e+03 1.73300e+03
2.17680e+03 1.78380e+03
2.12600e+03 1.78380e+03
2.10060e+03 1.78380e+03
2.10060e+03 1.80920e+03
2.12600e+03 1.80920e+03
2.10060e+03 1.83460e+03
2.12600e+03 1.83460e+03
2.15140e+03 1.83460e+03
2.23650e+03 1.84730e+03
1.99900e+03 1.84730e+03
1.88470e+03 1.84730e+03
2.10060e+03 1.86000e+03
2.12600e+03 1.86000e+03
2.10060e+03 1.88540e+03
2.12600e+03 1.88540e+03
2.17680e+03 1.88540e+03
2.15140e+03 1.91080e+03
2.12600e+03 1.91080e+03
2.10060e+03 1.91080e+03
2.10060e+03 1.93620e+03
2.12600e+03 1.93620e+03
2.17680e+03 1.93620e+03
2.22760e+03 1.93620e+03
2.12600e+03 1.96160e+03
2.10060e+03 1.96160e+03
1.79580e+03 1.98700e+03
1.82120e+03 1.98700e+03
1.84660e+03 1.98700e+03
1.87200e+03 1.98700e+03
1.89740e+03 1.98700e+03
1.94820e+03 1.98700e+03
2.05620e+03 1.98700e+03
2.10060e+03 1.98700e+03
2.12600e+03 1.98700e+03
2.25300e+03 1.98700e+03
2.30380e+03 1.98700e+03
2.38000e+03 1.98700e+03
2.40540e+03 1.98700e+03
2.02440e+03 1.40280e+03
2.15140e+03 1.40280e+03
2.07520e+03 1.70760e+03
2.17680e+03 1.70760e+03
2.35080e+03 1.73300e+03
2.35080e+03 1.84730e+03
3.65210e+03 1.01030e+03
3.72570e+03 1.01030e+03
3.72570e+03 1.08650e+03
3.65210e+03 1.08650e+03
3.72620e+03 1.14880e+03
3.80240e+03 1.14880e+03
3.85320e+03 1.14880e+03
3.80240e+03 1.17420e+03
3.70080e+03 1.17420e+03
3.60560e+03 1.19960e+03
3.70080e+03 1.19960e+03
3.80240e+03 1.19960e+03
3.85320e+03 1.19960e+03
4.02830e+03 1.31030e+03
3.95210e+03 1.31030e+03
3.72830e+03 1.31030e+03
3.65210e+03 1.31030e+03
3.65210e+03 1.38650e+03
3.72830e+03 1.38650e+03
3.95210e+03 1.38650e+03
4.02830e+03 1.38650e+03
3.85320e+03 1.12340e+03
3.95210e+03 1.08650e+03
4.02830e+03 1.08650e+03
4.02830e+03 1.01030e+03
3.95210e+03 1.01030e+03

3.2 具体程序实现

首先给出要运行各个算法的主程序以及导入的包：

import copy
import time
import numpy as np
# import matplotlib
# matplotlib.use("Agg")
import matplotlib.pyplot as plt
# import threading
from multiprocessing import Process
from matplotlib import animation

if __name__ == "__main__":
    continue_flag = True
    while continue_flag:
        print("这里有四个选择可以用于解TSP问题，你可以选择其中一个用于解答TSP问题，你也可以输入quit离开这里：\n"
              "    A. 修改圈算法\n"
              "    B. 遗传算法\n"
              "    C. 模拟退火法1\n"
              "    D. 模拟退火法2\n"
              "请选择其中一个编号，并按照\"编号 样例个数\"的格式输入：\n")
        command = input()
        if command == "quit":
            break
        args = command.split()
        print()
        if args[0] == "A":
            # Solution(args[1], modify_flag=True, num=int(args[2]))
            Solution(modify_flag=True, num=int(args[1]))
        if args[0] == "B":
            # Solution(args[1], ga_flag=True, num=int(args[2]))
            Solution(ga_flag=True, num=int(args[1]))
        if args[0] == "C":
            # Solution(args[1], sa_flag=True, num=int(args[2]))
            Solution(sa_flag=True, num=int(args[1]))
        if args[0] == "D":
            # Solution(args[1], sa_flag2=True, num=int(args[2]))
            Solution(sa_flag2=True, num=int(args[1]))
        print()

然后看看Solution类里面是如何通过__init__方法来读取文件并存储响应的数据的。

	# 文件、是否使用修改圈算法、是否使用遗传算法、是否使用模拟退火法1,、是否使用模拟退火法2、每个算法的运行案例个数
    def __init__(self, filename="d198.tsp", modify_flag=False, ga_flag=False, sa_flag=False, sa_flag2=False, num=1):
        file = open(filename, "r")
        lines = file.read().splitlines()
        file.close()
        # 获取城市数量，并准备存储数据的结构，城市编号默认从0~city_num - 1
        city_num = np.int(lines.pop(0))
        self.cities = np.zeros([city_num, 2], dtype=np.float)
        self.city_states = np.zeros([city_num, city_num], dtype=np.float)
        # 城市分布坐标
        for i in range(city_num):
            self.cities[i] = np.array(lines.pop(0).split(" ")).astype(np.float)
        # 城市邻接矩阵
        for i in range(city_num):
            x1, y1 = self.cities[i, 0], self.cities[i, 1]
            for j in range(city_num):
                x2, y2 = self.cities[j, 0], self.cities[j, 1]
                self.city_states[i, j] = np.sqrt(np.power(x1 - x2, 2) + np.power(y1 - y2, 2))
        # 各个算法的所有案例的结果
        results = {"ga_result": list(), "modify_result": list(), "sa_result": list(), "sa_result2": list()}
        for i in range(num):
            start = time.time()
            # 改良圈算法
            if modify_flag:
                results["modify_result"].append(self.modify_circle(i))
            # 遗传算法
            if ga_flag:
                results["ga_result"].append(self.ga_answer(50, 50, 0.2, 1, 10000, i))
            # 模拟退火法
            if sa_flag:
                results["sa_result"].append(self.sa_answer(i))
            if sa_flag2:
                results["sa_result2"].append(self.sa_answer2(i))
            print(time.time() - start)
        # print(results)
        if modify_flag:
            _results = results["modify_result"]
            for i in range(len(_results)):
                print("修改圈算法结果{0}".format(i), _results[i])
            print("均值", np.average(_results))
        if ga_flag:
            _results = results["ga_result"]
            for i in range(len(_results)):
                print("遗传算法结果{0}".format(i), _results[i])
            print("均值", np.average(_results))
        if sa_flag:
            _results = results["sa_result"]
            for i in range(len(_results)):
                print("模拟退火算法结果{0}".format(i), _results[i])
            print("均值", np.average(_results))
        if sa_flag2:
            _results = results["sa_result2"]
            for i in range(len(_results)):
                print("模拟退火算法结果{0}".format(i), _results[i])
            print("均值", np.average(_results))

现在来看看最为关键的sa_answer这个方法吧：

    # 模拟退火法1
    def sa_answer(self, tag):
    	 # 随机的解：我选择从0~197再回到0这样的路径
        city_num = self.city_states.shape[0]
        greed_path = [0]
        rest = list(range(1, city_num))
        greed_path.extend(rest)
        greed_path.append(0)
        all_result = [greed_path]
        cost = np.sum([self.city_states[greed_path[i], greed_path[i + 1]] for i in range(city_num)])
        best_cost = cost
        # print(cost)
        # 初温
        current_temperature = 100
        while current_temperature > 0.1:	# 末温
            for _ in range(1000):
                index1, index2 = np.random.randint(1, city_num, size=2).tolist()
                while index1 == index2:
                    index2 = np.random.randint(1, city_num)
                if index2 < index1:
                    index1, index2 = index2, index1
                path1, path2, path3 = copy.deepcopy(greed_path), copy.deepcopy(greed_path), copy.deepcopy(greed_path)
                # 交换index1与index2位置的城市
                path1[index1], path1[index2] = path1[index2], path1[index1]
                cost1 = np.sum([self.city_states[path1[i], path1[i + 1]] for i in range(city_num)])
                result_path = path1
                result_cost = cost1
                # 将index2的城市插入到index1前
                # path2.insert(index1, path2.pop(index2))
                temp = path2[index2]
                for i in range(index2 - index1):
                    path2[index2 - i] = greed_path[index2 - i - 1]
                path2[index1] = temp
                cost2 = np.sum([self.city_states[path2[i], path2[i + 1]] for i in range(city_num)])
                if result_cost > cost2:
                    result_path = path2
                    result_cost = cost2
                # 将index1与index2间的城市逆转
                # path3[index1: index2] = path3[index2 - 1: index1 - 1: -1] 这个好像有问题，虽然可以达到一样的变换效果，但是在退火过程中最佳解死活退不下去
                for i in range(index2 - index1 + 1):
                    path3[index1 + i] = greed_path[index2 - i]
                cost3 = np.sum([self.city_states[path3[i], path3[i + 1]] for i in range(city_num)])
                if result_cost > cost3:
                    result_path = path3
                    result_cost = cost3
                # 选出最优解
                if result_cost < cost or np.exp(-1 * np.abs(cost - result_cost) / current_temperature) >= np.random.rand():
                    greed_path = result_path
                    cost = result_cost
                    if cost < best_cost:
                        best_path = copy.deepcopy(greed_path)
                        best_cost = cost
                        all_result.append(best_path)
            # 温度变化
            current_temperature *= 0.95
            # print("当前局部最优", cost, "当前最优解", best_cost, "当前温度", current_temperature)
        print(best_cost)
        # print("results长度", len(all_result))
        show_thread = Process(target=self.show, args=(all_result, tag, ))
        show_thread.start()
        return best_cost

其实这里的三种方法可以按照一定的比例来分配次数，看看他们的那个效果更好。

最后贴一贴画图/保存视频的函数：

    # 画动图/保存视频
    def show(self, all_result, tag):
        # time.sleep(10)
        len1 = len(all_result)
        len2 = min(300, len1)
        sampling = np.floor(np.linspace(0, len1 - 1, len2, endpoint=True)).astype(np.int)
        figure, ax = plt.subplots()
        ax.scatter(self.cities[:, 0], self.cities[:, 1])
        _line = np.array([self.cities[i] for i in all_result[0]])
        line, = ax.plot(_line[:, 0], _line[:, 1], color="r")

        def init():
            return line

        def update(frame):
            frame = all_result[frame]
            _line2 = np.array([self.cities[i] for i in frame])
            line.set_ydata(_line2[:, 1])
            line.set_xdata(_line2[:, 0])
            return line

        anim = animation.FuncAnimation(fig=figure, func=update, init_func=init, interval=50, frames=sampling,
                                       repeat=False)
        # writer = animation.writers['ffmpeg'](fps=15, metadata=dict(artist='Me'), bitrate=1800)
        # anim.save("video{0}.mp4".format(tag), writer=writer)

        # plt.ion()
        # plt.pause(10)
        # plt.close("all")

        plt.title("figure{0}".format(tag))
        plt.show()

3.3 结果

现在可以看看输出结果，参考最优解是15780：

显然效果很好，不仅超过了10%，甚至还可以达到2%这样的好结果，平均基本是3%这样的结果了，要知道198个城市的路线可是可以组合出198的阶乘(?)个结果的，如果用深度搜索或者广度搜索这样的方法，耗费的时间难以估量，而这里只要120来秒，还是我每次领域操作后都重新计算新解，没有在原来解的基础上计算的原因，而且我用的是比C/C++要满30到100倍的python语言，更重要的是我开了很多的应用，这些理所当然的会影响到pycharm中代码的执行效率(平时用pycharm运行这个也就耗费70到90秒而已)，可见其实这样的效果已经不错了。
其中最优解的路线图如下：

4.后续感悟

其实模拟退火法中关键的地方有初温、末温、内层循环遍历次数以及温度变化控制参数，当然，那些产生新解的领域操作也是重点。其中初温、末温、内层循环遍历次数、温度变化参控制数应该与要解答的问题的规模有关，比如TSP中城市的规模达到400，那么这些参数都得重新调整(具体怎么调整我也不清楚，惭愧，或者可以看看一些相关博客或者论文？)。
至于领域操作关系到产生新解的可能，多个邻域操作的扰动性会更大，当然也要调整好使用它们的比例(我这里就不继续调整了，因为结果好像已经很不错)。