自适应大邻域搜索算法

自适应大邻域搜索算法(Adaptive Large Neighborhood Search)是由Ropke与Pisinger在2006年提出的一种启发式方法，其在邻域搜索的基础上增加了的对算子的作用效果的衡量，使算法能够自动选择好的算子对解进行破坏与修复，从而有一定几率得到更好的解。在邻域搜索算法中，有的算法可以只使用一种邻域，如模拟退火算法，因此它仅仅搜索了解空间的一小部分，找到全局最优的概率较小，它的优势之一是可以避免陷入局部最优；而有的算法可以使用多种算子，如变邻域搜索算法，它通过在当前解的多个邻域中寻找更满意的解，能够大大提高算法在解空间的搜索范围，但是它在使用算子时盲目地将每种算子形成的邻域结构都搜索一遍，缺少了一些启发式信息的指导。而自适应大邻域搜索算法就弥补了这种不足，这种算法根据算子的历史表现与使用次数选择下一次迭代使用的算子，通过算子间的相互竞争来生成当前解的邻域结构，而在这种结构中有很大概率能够找到更好的解。

自适应大邻域搜索算法的步骤如下：
1.生成初始解X，最优解X’ = X，初始化算法参数。
2.重复如下步骤直到停止准则：
2.1轮盘赌方法根据算子权重选择破坏与修复算子对当前解进行操作得到新解X’’。
2.2更新算子使用次数
2.3如f(X’’) < f(X)，则X = X’’。
2.4如f(X’’) < f(X’)，则X’ = X’’。
2.5如f(X’’) > f(X)，则以一定的概率接受该解作为当前解。
2.5更新算子的分数。
2.6更新算子的权重。权重更新式子如下：

式中，Wd为算子权重，Sd为算子分数，Ud为算子的使用次数。
3.输出最优解。

对于TSP问题，ALNS往往能在较短时间内找到满意解：

# Adaptive Large Neighborhood Search
import numpy as np
import random as rd
import copy

distmat = np.array([[0,350,290,670,600,500,660,440,720,410,480,970],
                 [350,0,340,360,280,375,555,490,785,760,700,1100],
                 [290,340,0,580,410,630,795,680,1030,695,780,1300],
                 [670,360,580,0,260,380,610,805,870,1100,1000,1100],
                 [600,280,410,260,0,610,780,735,1030,1000,960,1300],
                 [500,375,630,380,610,0,160,645,500,950,815,950],
                 [660,555,795,610,780,160,0,495,345,820,680,830],
                 [440,490,680,805,735,645,495,0,350,435,300,625],
                 [720,785,1030,870,1030,500,345,350,0,475,320,485],
                 [410,760,695,1100,1000,950,820,435,475,0,265,745],
                 [480,700,780,1000,960,815,680,300,320,265,0,585],
                 [970,1100,1300,1100,1300,950,830,625,485,745,585,0]])

def disCal(path):             # calculate distance
    distance = 0
    for i in range(len(path) - 1):
        distance += distmat[path[i]][path[i + 1]]
    distance += distmat[path[-1]][path[0]]
    return distance

def selectAndUseDestroyOperator(destroyWeight,currentSolution):    # select and use destroy operators
    destroyOperator = -1
    sol = copy.deepcopy(currentSolution)
    destroyRoulette = np.array(destroyWeight).cumsum()
    r = rd.uniform(0, max(destroyRoulette))
    for i in range(len(destroyRoulette)):
        if destroyRoulette[i] >= r:
            if i == 0:
                destroyOperator = i
                removedCities = randomDestroy(sol)
                destroyUseTimes[i] += 1
                break
            elif i == 1:
                destroyOperator = i
                removedCities = max3Destroy(sol)
                destroyUseTimes[i] += 1
                break
    return sol,removedCities,destroyOperator

def selectAndUseRepairOperator(repairWeight,destroyedSolution,removeList):    # select and use repair operators
    repairOperator = -1
    repairRoulette = np.array(repairWeight).cumsum()
    r = rd.uniform(0, max(repairRoulette))
    for i in range(len(repairRoulette)):
        if repairRoulette[i] >= r:
            if i == 0:
                repairOperator = i
                randomInsert(destroyedSolution,removeList)
                repairUseTimes[i] += 1
                break
            elif i == 1:
                repairOperator = i
                greedyInsert(destroyedSolution,removeList)
                repairUseTimes[i] += 1
                break
    return destroyedSolution,repairOperator

def randomDestroy(sol):    # randomly remove 3 cities
    solNew = copy.deepcopy(sol)
    removed = []
    removeIndex = rd.sample(range(0, distmat.shape[0]), 3)
    for i in removeIndex:
        removed.append(solNew[i])
        sol.remove(solNew[i])
    return removed

def max3Destroy(sol):      # remove city with 3 longest segments
    solNew = copy.deepcopy(sol)
    removed = []
    dis = []
    for i in range(len(sol) - 1):
        dis.append(distmat[sol[i]][sol[i + 1]])
    dis.append(distmat[sol[-1]][sol[0]])
    disSort = copy.deepcopy(dis)
    disSort.sort()
    for i in range(3):
        if dis.index(disSort[len(disSort) - i - 1]) == len(dis) - 1:
            removed.append(solNew[0])
            sol.remove(solNew[0])
        else:
            removed.append(solNew[dis.index(disSort[len(disSort) - i - 1]) + 1])
            sol.remove(solNew[dis.index(disSort[len(disSort) - i - 1]) + 1])
    return removed

def randomInsert(sol,removeList):     # randomly insert 3 cities
    insertIndex = rd.sample(range(0, distmat.shape[0]), 3)
    for i in range(len(insertIndex)):
        sol.insert(insertIndex[i],removeList[i])

def greedyInsert(sol,removeList):     # greedy insertion
    dis = float("inf")
    insertIndex = -1
    for i in range(len(removeList)):
        for j in range(len(sol) + 1):
            solNew = copy.deepcopy(sol)
            solNew.insert(j,removeList[i])
            if disCal(solNew) < dis:
                dis = disCal(solNew)
                insertIndex = j
        sol.insert(insertIndex,removeList[i])
        dis = float("inf")

T = 100
a = 0.97
b = 0.5
wDestroy = [1 for i in range(2)]     # weights of the destroy operators
wRepair = [1 for i in range(2)]      # weights of the repair operators
destroyUseTimes = [0 for i in range(2)]     #The number of times the destroy operator has been used
repairUseTimes = [0 for i in range(2)]      #The number of times the repair operator has been used
destroyScore = [1 for i in range(2)]      # the score of destroy operators
repairScore = [1 for i in range(2)]       # the score of repair operators
solution = [i for i in range(distmat.shape[0])] # initial solution
bestSolution = copy.deepcopy(solution)  # best solution
iterx, iterxMax = 0, 100

if __name__ == '__main__':
    while iterx < iterxMax:       # while stop criteria not met
        while T > 10:
            destroyedSolution,remove,destroyOperatorIndex = selectAndUseDestroyOperator(wDestroy,solution)
            newSolution,repairOperatorIndex = selectAndUseRepairOperator(wRepair,destroyedSolution,remove)

            if disCal(newSolution) <= disCal(solution):
                solution = newSolution
                if disCal(newSolution) <= disCal(bestSolution):
                    bestSolution = newSolution
                    destroyScore[destroyOperatorIndex] += 1.5     # update the score of the operators
                    repairScore[repairOperatorIndex] += 1.5
                else:
                    destroyScore[destroyOperatorIndex] += 1.2
                    repairScore[repairOperatorIndex] += 1.2
            else:
                if rd.random() < np.exp(- disCal(newSolution)/ T):  # the simulated annealing acceptance criteria
                    solution = newSolution
                    destroyScore[destroyOperatorIndex] += 0.8
                    repairScore[repairOperatorIndex] += 0.8
                else:
                    destroyScore[destroyOperatorIndex] += 0.6
                    repairScore[repairOperatorIndex] += 0.6

            wDestroy[destroyOperatorIndex] = wDestroy[destroyOperatorIndex] * b + (1 - b) * \
                                             (destroyScore[destroyOperatorIndex] / destroyUseTimes[destroyOperatorIndex])
            wRepair[repairOperatorIndex] = wRepair[repairOperatorIndex] * b + (1 - b) * \
                                           (repairScore[repairOperatorIndex] / repairUseTimes[repairOperatorIndex])
            # update the weight of the operators

            T = a * T 
        iterx += 1
        T = 100

    print(bestSolution)
    print(disCal(bestSolution))