3-opt方法解旅行商问题的三边选取方法

解旅行商问题，最简单有效的算法是k-opt方法。常见的是2-opt和3-opt，相关算法的介绍文章网上很多，这里不加赘述，只对三边选择方法进行说明。
三条边的选取，应当满足两个原则：
原则一：不允许相接；违反这条，无法成功求解。
原则二：如果随机选择，应对服从均匀分布。违反这条，新解的出现频率不同，搜索范围被缩小，不利于找寻全局最优解的目标。
本程序采用随机抽取3个数字[0，n-1]内的数字，分别作为三条边的起点，也就是离开城市。具体抽样方式：
方法一：直接使用随机函数抽三点，如不可行，则重抽。
具体做法是：随机抽[0,n-1]内的3个不重复的整数，如果发现有任意两个数的差小于2（意味着两条边前后相连），则判为不可行的组合，重新再抽一次。
此方法好处是简单，适用于问题规模不太小的情况，推荐用于30城以上。
缺点是如果城市太少，因为符合的概率太低，则程序会陷入不断的重抽，效率低。因为随机抽取3个整数是可行组合（即：互不相连）的概率为$\frac{\frac{n(n-4)(n-5)}{6}}{\left(\begin{array}{c}n\\ 3\end{array}\right)}=\frac{(n-4)(n-5)}{(n-1)(n-2)}$，随城市增加而增加，当达30城时为80%。
如果问题规模大，重抽机会很小，效率高。

方法二：三边依次在所有允许的范围内按均匀分布的概率抽取。
具体做法：
第一条边起点loc1在[0,n-1]范围内随意抽。第二条边起点loc2在除[loc1-1,loc1+1]范围以外的n-3个点内随意抽。第三条边起点loc3的抽取要分两种情况：
如果|loc1-loc2|=2，则允许抽的点还有n-5个；如果|loc1-loc2|>2，则允许抽的点还有n-6个。两者相加一共有$\frac{n(n-4)(n-5)}{6}$种组合。
最后上代码。

import numpy as np
distance_matrix=np.load(r'distance_matrix.npy')
global city_num
city_num=len(distance_matrix)

def distance_total_func(route):
    return distance_matrix[route[np.arange(city_num)-1],route].sum()

def change3(route,a,b,c):#要求a<b<c
    route1=np.repeat(route.reshape((1,city_num)),7,axis=0)
    route1[0]=np.concatenate((route[:a+1],route[c:b:-1],route[a+1:b+1],route[c+1:]))
    route1[1]=np.concatenate((route[:a+1],route[b+1:c+1],route[b:a:-1],route[c+1:]))
    route1[2]=np.concatenate((route[:a+1],route[b:a:-1],route[c:b:-1],route[c+1:]))
    route1[3]=np.concatenate((route[:a+1],route[b+1:c+1],route[a+1:b+1],route[c+1:]))
    #后三种情况其实是2-opt
    route1[4,a+1:b+1]=route[b:a:-1]
    route1[5,b+1:c+1]=route[c:b:-1]
    route1[6,a+1:c+1]=route[c:a:-1]
    return route1

def LOC():#随机选取互不相连的三条边
    if city_num<30:
        loc1=np.random.randint(0,city_num-1)
        loc2=np.random.randint(loc1+2,loc1-2+city_num)%city_num
        if (loc2-loc1)%city_num==2:
            loc3=np.random.randint(loc2+2,loc2-4+city_num-6)%city_num
        elif (loc1-loc2)%city_num==2:
            loc3=np.random.randint(loc1+2,loc1-4+city_num-6)%city_num
        else:
            loc3=np.random.randint(0,city_num-7)
            if loc3<(loc1-loc2-3)%city_num:
                loc3=(loc3+loc2+2)%city_num
            else:
                loc3=(loc3+loc2+5)%city_num
        loc=[loc1,loc2,loc3]
        loc.sort()
        return loc
    else:
        loc=[0,0,0]
        while loc[1]-loc[0]<2 or loc[2]-loc[1]<2:
            loc=np.random.choice(np.arange(city_num-1),size=3,replace=False)
            loc=np.sort(loc)
        return loc.tolist()
    
neighbor=((city_num)*(city_num-4)*(city_num-5))//6#全部可行组合的范围
best_route=np.random.permutation(city_num)
best=distance_total_func(best_route)

i=0
while i<neighbor:#可自行改小，这是上限，
    loc=LOC()
    routes_this=change3(best_route,loc[0],loc[1],loc[2])
    result=np.apply_along_axis(distance_total_func,axis=1,arr=routes_this)
    if min(result)<best:#发现更优解
        best=min(result)
        best_route=routes_this[result.argmin()]
        #print(r"%.2f"%best)#迭代过程展示
        i=0#重新计数
    else:
        i+=1
print(r"best:%.2f"%best,best_route)