zoukankan      html  css  js  c++  java
  • 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()
    
  • 相关阅读:
    一个很好的Delphi博客
    Android开发之adb无法连接
    J2EE--常见面试题总结 -- (二)
    J2EE--常见面试题总结 -- ( 一)
    面向接口编程实现不改代码实现Redis单机/集群之间的切换
    Gradle sync failed 异常
    Dubbo+Zookeeper
    JdbcTemplate的使用
    Spring -- 配置bean
    浅析 @PathVariable 和 @RequestParam(转发,非原创)
  • 原文地址:https://www.cnblogs.com/Martrix-revolution/p/15412470.html
Copyright © 2011-2022 走看看