旅行商问题-模拟退火求解（Python）

旅行商问题-模拟退火求解

使用模拟退火算法解决旅行商问题， 使用Python实现。
旅行商问题是一个典型的NP难问题， 当问题规模很大时，因计算机算力局限， 无法使用如动态规划等全局枚举求解最优解的算法。
这时可以使用启发算法，模拟退火、遗传算法等，在有限时间、有限的解空间，求解最优解（不一定得出全局最优）。

旅行商问题定义

旅行推销员问题（英语：Travelling salesman problem, TSP）是这样一个问题：给定一系列城市和每对城市之间的距离，求解访问每一座城市一次并回到起始城市的最短回路。它是组合优化中的一个NP难问题，在运筹学和理论计算机科学中非常重要。 (百度搜来的)

问题地点代价邻接矩阵

import matplotlib.pyplot as plt
import numpy as np

# 边界条件， 非对称
weights = np.array([[0, 10, 15, 20],
         [5, 0, 9, 10],
          [6, 13,  0,  12],
          [8, 8, 9 ,  0],
         ])
weights

array([[ 0, 10, 15, 20],
       [ 5,  0,  9, 10],
       [ 6, 13,  0, 12],
       [ 8,  8,  9,  0]])

算法流程图

np.random.seed(100)

定义初始解函数与目标函数

dim = weights.shape[0]

def random_path(start=0, dim=4):
    # 随机生成解
    p = list(range(dim))
    p.remove(start)
    d = [start]
    print(f"p: {p}")
    while p:
        index = np.random.randint(0, len(p))
        v = p.pop(index)
        d.append(v)
        
    d.append(start)
    return d

def objective_func_test(path):
    w = 0
    i0 = path[0]
    for i in path[1:]:
#         print(weights[i0, i])
        w += weights[i0, i]
        i0 = i
        
    return w

path = random_path(start=0, dim=dim)
print(path)

p: [1, 2, 3]
[0, 1, 2, 3, 0]

objective_func_test(path)

39

算法实现

class SAA(object):
    
    def __init__(self, path):
        self.path = path
        self._run_data = []
    
    def get_new_path(self, path):
        i = np.random.randint(1, len(path) - 1) # 随机交换两个地点
#         i = 1
        print(f"=={i}")
        if i == len(path) - 2:
            i1 = 1
        else:
            i1 = i + 1

        print(i, i1)
        new = [i for i in path]
        new[i], new[i1] = new[i1], new[i]
    
        return new
    
    def probability(self, delta, T):
        return np.exp(-(delta / T))
    
    def objective_func(self, path):
        w = 0
        i0 = path[0]
        for i in path[1:]:
    #         print(weights[i0, i])
            w += weights[i0, i]
            i0 = i

        return w

    
    def temperature_down(self, T, a = 0.99, num=0):
        #  指数下降温度
        return T * a
    
    def set_run_data(self, n, path, value, T):
        self._run_data.append((n, path, value, T))
    
    def get_run_data(self):
        return self._run_data
    
    def run(self):
        T = weights.sum() * 100
        print(f"初始T：{T}")
        path = self.path
        
        all_best = None
        all_best_path = None
        trange_best = None
        
        n = 0
        while True:
            old_w = self.objective_func(path)
            new_path = self.get_new_path(path)
            new_w = self.objective_func(new_path)
            
            print(f"=={n}: {old_w} -> {new_w}")
            
            value = old_w
            p = 0
            delta =  new_w - old_w
            if delta < 0:
                path = new_path
                trange_best = new_w
                value = new_w
                all_best = new_w
                all_best_path = path
                
            else:
                p = self.probability(delta, T)
                pr = np.random.random()
                if pr < p:
                    path = new_path
                    trange_best = new_w
                    value = new_w
                    print(f"==={n}概率接受解p: {p} T: {T} δ：{delta} ({pr},{p})")
            
            self.set_run_data(n, path, value, T)
            T = self.temperature_down(T, num=n)
            if T < 1 or p >= 1 or n > 100_00:
                print(f"==={n}退出：T：{T}， P：{p}")
                break
            
            n += 1
        
        self.path = all_best_path
        return path, trange_best, all_best, all_best_path
    
    

saa = SAA(path)
saa.run()

初始T：12500
==1
1 2
==0: 39 -> 46
===0概率接受解p: 0.9994401567707347 T: 12500 δ：7 (0.14860484007536512,0.9994401567707347)
==1
1 2
==1: 46 -> 39
==3
3 1
==2: 39 -> 47
===2概率接受解p: 0.9993472185621816 T: 12251.25 δ：8 (0.18646718316338273,0.9993472185621816)
==3
3 1
==3: 47 -> 39
==3
3 1
==4: 39 -> 47
===4概率接受解p: 0.9993339688696272 T: 12007.450125 δ：8 (0.45273990285923016,0.9993339688696272)
==1
1 2
==5: 47 -> 40
==1
1 2
==6: 40 -> 47
===6概率接受解p: 0.9994053687832825 T: 11768.5018675125 δ：7 (0.06368104287476795,0.9994053687832825)
2 3
==915: 35 -> 39
===915概率接受解p: 0.042654647195248654 T: 1.2679819454842354 δ：4 (0.029660839697256525,0.042654647195248654)

3 1
==937: 35 -> 40
==1
1 2
==938: 35 -> 43
===938退出：T：0.9962256975051745， P：0.00035266014513277063





([0, 1, 3, 2, 0], 35, 35, [0, 1, 3, 2, 0])

绘制求解过程曲线

data = saa.get_run_data()
Yt = [i[-1] for i in data]
Yv = [i[2] for i in data]
X = [i[0] for i in data]

# X = np.linspace(0, 10, 100)
# Y = np.exp(-X)
plt.figure(figsize=(15, 5))
plt.plot(X, Yt)
plt.show()

plt.figure(figsize=(15, 5))
plt.plot(X, Yv, ".b")
plt.show()