用python可视化模拟退火算法

阳光穿透心脏的1/2处 2022-05-31 02:13 212阅读 0赞

关于这个算法的解释的话，直接百度百科吧（摊手）：

[https://baike.baidu.com/item/%E6%A8%A1%E6%8B%9F%E9%80%80%E7%81%AB%E7%AE%97%E6%B3%95/355508?fr=aladdin][https_baike.baidu.com_item_E6_A8_A1_E6_8B_9F_E9_80_80_E7_81_AB_E7_AE_97_E6_B3_95_355508_fr_aladdin]

按我个人的理解的话，是解决组合优化的问题是，使用随机化的方法得到新解，如果新解比旧解要好，那么就接受。如果新解没有旧解好，那么也按一定概率\[exp(-delta\_f/T)\]接受。T是一个温度，内循环就产生新解直到达到平稳，外循环就退火（缓慢的速率温度）。到结束温度时，会收敛到最优解。那么我用的示例是旅行商问题。直接贴代码吧。

from matplotlib import pyplot as plt
    import numpy as np
    
    def coordinate_init(size):
        #产生坐标字典
        coordinate_dict = {}
        coordinate_dict[0] = (0, 0)#起点是（0，0）
        for i in range(1, size + 1):#顺序标号随机坐标
            coordinate_dict[i] = (np.random.uniform(0, size), np.random.uniform(0, size))
        coordinate_dict[size + 1] = (0, 0)#终点是（0,0)
        return coordinate_dict
    
    def distance_matrix(coordinate_dict,size):#生成距离矩阵
        d=np.zeros((size+2,size+2))
        for i in range(size+1):
            for j in range(size+1):
                if(i==j):
                    continue
                if(d[i][j]!=0):
                    continue
                x1 = coordinate_dict[i][0]
                y1 = coordinate_dict[i][1]
                x2 = coordinate_dict[j][0]
                y2 = coordinate_dict[j][1]
                distance=np.sqrt((x1-x2)**2+(y1-y2)**2)
                if(i==0):
                    d[i][j]=d[size+1][j]=d[j][i]=d[j][size+1]=distance
                else:
                    d[i][j]=d[j][i]=distance
        return d
    
    def path_length(d_matrix,path_list,size):#计算路径长度
        length=0
        for i in range(size+1):
            length+=d_matrix[path_list[i]][path_list[i+1]]
        return length
    
    def new_path(path_list,size):
        #二交换法
        change_head = np.random.randint(1,size+1)
        change_tail = np.random.randint(1,size+1)
        if(change_head>change_tail):
            change_head,change_tail=change_tail,change_head
        change_list = path_list[change_head:change_tail + 1]
        change_list.reverse()#change_head与change_tail之间的路径反序
        new_path_list = path_list[:change_head] + change_list + path_list[change_tail + 1:]
        return change_head,change_tail,new_path_list
    
    def diff_old_new(d_matrix,path_list,new_path_list,head,tail):#计算新旧路径的长度之差
        old_length=d_matrix[path_list[head-1]][path_list[head]]+d_matrix[path_list[tail]][path_list[tail+1]]
        new_length=d_matrix[new_path_list[head-1]][new_path_list[head]]+d_matrix[new_path_list[tail]][new_path_list[tail+1]]
        delta_p=new_length-old_length
        return delta_p
    
    
    T_start=2000#起始温度
    T_end=1e-20#结束温度
    a=0.995#降温速率
    Lk=50#内循环次数,马尔科夫链长
    size=20
    coordinate_dict=coordinate_init(size)
    print(coordinate_dict)#打印坐标字典
    path_list=list(range(size+2))#初始化路径
    d=distance_matrix(coordinate_dict,size)#距离矩阵的生成
    best_path=path_length(d,path_list,size)#初始化最好路径长度
    print('初始路径:',path_list)
    print('初始路径长度:',best_path)
    best_path_temp=[]#记录每个温度下最好路径长度
    best_path_list=[]#用于记录历史上最好路径
    balanced_path_list=path_list#记录每个温度下的平衡路径
    balenced_path_temp=[]#记录每个温度下平衡路径(局部最优)的长度
    while T_start>T_end:
        for i in range(Lk):
            head, tail, new_path_list = new_path(path_list, size)
            delta_p = diff_old_new(d, path_list, new_path_list, head, tail)
            if delta_p < 0:#接受状态
                balanced_path_list=path_list = new_path_list
                new_len=path_length(d,path_list,size)
                if(new_len<best_path):
                    best_path=new_len
                    best_path_list=path_list
            elif np.random.random() < np.exp(-delta_p / T_start):#以概率接受状态
                path_list = new_path_list
        path_list=balanced_path_list#继承该温度下的平衡状态（局部最优）
        T_start*=a#退火
        best_path_temp.append(best_path)
        balenced_path_temp.append(path_length(d,balanced_path_list,size))
    print('结束温度的局部最优路径:',balanced_path_list)
    print('结束温度的局部最优路径长度:',path_length(d,balanced_path_list,size))
    print('最好路径:',best_path_list)
    print('最好路径长度:',best_path)
    x=[]
    y=[]
    for point in best_path_list:
        x.append(coordinate_dict[point][0])
        y.append(coordinate_dict[point][1])
    plt.figure(1)
    plt.subplot(311)
    plt.plot(balenced_path_temp)#每个温度下平衡路径长度
    plt.subplot(312)
    plt.plot(best_path_temp)#每个温度下最好路径长度
    plt.subplot(313)
    plt.scatter(x,y)
    plt.plot(x,y)#路径图
    plt.grid()
    plt.show()

说实话，我最讨厌别人的博客只有代码没有输出，所以，输出还是要有的：

{0: (0, 0), 1: (8.0813445997638, 19.378598885734345), 2: (5.9336947559652735, 18.43334356701462), 3: (2.4187354288286844, 11.208318685158517), 4: (11.364680404115896, 17.838011084801682), 5: (14.57899298590396, 10.443803487389676), 6: (14.475687688871977, 6.430692753878264), 7: (15.564547884718248, 1.8489819656041484), 8: (2.762034307781256, 1.3609270597085188), 9: (2.4449326896325507, 13.971590679990305), 10: (6.376590999061471, 13.590205047457003), 11: (2.0008565516450005, 0.25779568392405805), 12: (8.756514935267091, 15.616855386568506), 13: (8.19333746755698, 14.267482040877905), 14: (6.688671060971291, 7.594062408648769), 15: (16.826772155106706, 18.439257279177564), 16: (0.17109735207725185, 9.772548332790812), 17: (10.057710887083482, 4.213811533920735), 18: (14.404568409630825, 16.182396066315434), 19: (6.522071276732497, 4.044136354883117), 20: (3.2651888440785837, 14.12486494048631), 21: (0, 0)}
    初始路径: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]
    初始路径长度: 225.0507842833548
    结束温度的局部最优路径: [0, 11, 8, 19, 14, 17, 7, 6, 5, 18, 15, 4, 1, 2, 12, 13, 10, 20, 9, 3, 16, 21]
    结束温度的局部最优路径长度: 78.12501363830815
    最好路径: [0, 11, 8, 19, 14, 17, 7, 6, 5, 18, 15, 4, 1, 2, 12, 13, 10, 20, 9, 3, 16, 21]
    最好路径长度: 78.12501363830815

距离矩阵的话，太大没必要输出了，还有个图：

![SouthEast][]

还是说下思路吧：为了不用什么文本文件输入什么的，我打算先是随机生成一个坐标字典，然后通过这个坐标字典，算出对应的距离矩阵。然后按顺序初始化路径。好了，然后就是退火了。初始温度我定了2000，其实定大一点比较好吧，虽说这样算法就比较慢，退火的速率我定了0.995，这样算是比较慢的退火，这样比较能找到全局最优解吧而不是局部最优解。生成新的路径我是用了二交换法，也就是随机选两个点，然后这两点之间的路径反序。比较跟别人的博客不一样的是，我做了可视化处理，记录了每个温度下的最好路径长度和平衡路径长度。然后最优路径图也plot出来了。还有一点要注意的点是，一定要弄距离矩阵，不然几万次循环都开方的话，速度太慢了，而且开方多了，误差叠加。。。真的会很大。以上。

[https_baike.baidu.com_item_E6_A8_A1_E6_8B_9F_E9_80_80_E7_81_AB_E7_AE_97_E6_B3_95_355508_fr_aladdin]: https://baike.baidu.com/item/%E6%A8%A1%E6%8B%9F%E9%80%80%E7%81%AB%E7%AE%97%E6%B3%95/355508?fr=aladdin
[SouthEast]: /images/20220531/62b04ab36b04459da028fb494194a372.png