ALNS自适应大邻域搜索(Adaptive Large Neighborhood Search)是由LNS(Large Neighborhood Search)发展而来的。LNS又是属于NS(Neighborhood Search)下的VLNS(Very Large-Scale Neighborhood Search)的一种。
- NS
简而言之是将求解得出的bit值,依次对其次位进行取反操作。
- VLNS
随着算例规模的增大,邻域搜索算法的邻域规模呈指数增长或者当邻域太大而不能在实际中明确搜索时,我们把这类邻域搜索算法归类于VLNS。与NS并无太大区别。
- LNS
在LNS中,邻域是由destroy和repair方法隐式定义的。destroy方法会破坏当前解的一部分,而后repair方法会对被破坏的解进行重建。所以,LNS是首先通过利用destroy方法破坏解x,然后利用repair方法重建解x,从而得到的一系列解的集合。LNS不会搜索解的整个邻域,而只是对该邻域进行采样搜索。
- ALNS
ALNS在LSN的基础上,允许在同一个搜索中使用多个destroy和repair方法来获得当前解的邻域。相较于LNS,ALNS在初始化时会对destory和repair方法进行加权。求邻域时会对算子权重进行比较,权重越大,选中该算子的概率就越大。除此之外,ALNS会对destroy和repair方法的权重进行动态调整。具体是:完成一次迭代搜索后,会使用函数对每组结果的好坏进行评估。
ω d = { ω d , u d = 0 ( 1 − ρ ) ω d + ρ s d u d , u d > 0 \omega_d = \begin{cases} \omega_d,\,\,u_d=0\\ (1-\rho)\omega_d+\rho\frac{s_d}{u_d},\,\,u_d>0\\ \end{cases} ωd={ωd,ud=0(1−ρ)ωd+ρudsd,ud>0
其中:
ω d \omega_d ωd为算子权重, s d s_d sd为算子分数, u d u_d ud为算子的使用次数, ρ \rho ρ为权重更新系数(控制权重变化的速度)。自适应大领域搜索中通常采用模拟退火算法Metropolis准则,在一定概率下接受劣解。
轮盘赌选择算子
简而言之:得到新的权重后,算法基于轮盘赌的思想对算子进行选择,使算子被选中的概率与其权重表现成正比。
关于destory和repair
在LNS和ALNS中,destory和repair是需要成对对一个解使用的。每个解经过一次destory和repair后,实则是通过该变换获得了一个新的落在邻域中的解,其过程属于CVRP(Capacitated Vehicle Routing Problem - 仅考虑车辆容量约束的 [VRP - 车辆路径规划] 问题)。
比如在一个CVRP问题中,将destroy方法定义为:移除当前解x中15%的customers。假如当前的解x中有100名customers,那么就有C(100,15)= 100!/(15!×85!) =2.5×1017种移除的方式。并且,根据每一种移除方式,又可以有很多种修复的方法。这样下来,一对destroy和repair方法能生成非常多的邻居解,而这些邻居解的集合,就是邻域了。
destory 破坏算子
破坏解的方法主要有随机移除、最差移除、相似移除等。其中,随机移除删除当前解决方案中的任意节点;最差移除删除当前解决方案中距离较长的路段。
先破坏再修复,这一步得到破坏后的解和移除的节点列表。如,当前解(TSP城市访问顺序)为【1, 4, 2, 3】,随机删除节点4后,得到破坏后的解【1, 2, 3】
repair 修复算子
修复解的方法主要有随机插入、贪婪插入等。其中,随机插入将移除的节点逐个插入到破坏后解的任意位置;贪婪插入将移除的节点插到距离成本最小,即插入后总路径最短的位置中。
修复操作的对象是破坏后的解和移除的节点列表,将移除的节点重新插入到破坏后的解中,得到完整的、新的一组解。如,将2.1中的移除节点4,随机插入到破坏后的解【1, 2, 3】中,可得【1, 2, 3, 4】
需要选择至少两组破坏和修复算子加入算法,更多的算子介绍详见参考文献,也可以根据问题特性自行设计其他的破坏或修复算子。
算法介绍
是由Ropke与Pisinger在2006年提出的一种启发式方法,其在邻域搜索的基础上增加了对算子的作用效果的衡量,使算法能够自动选择好的算子对解进行破坏与修复,得到一个成本(cost)最低的一个解。由于算法的自动反馈更新权重机制,根据算子的历史表现与使用次数选择下一次迭代使用的算子,通过算子间的相互竞争来生成当前解的邻域结构,而在这种结构中有很大概率能够找到更好的解,同时可以尽可能地避免出现局部最优解的情况。
算法设计[4]
- 将随机生成的船舶出港次序作为一个可行的初始化解,用 S0 表示,令最优解 S*=S0;
- 设计破坏算子和修复算子。本文采用的2种破坏算子分别为基于最差船舶的破坏算子和随机破坏算子,基于最差船舶的破坏算子是指将疏散时间最长的船舶作为“最差”船舶,而随机破坏算子是指从当前解中随机选取 η \eta η艘船舶移除。本文采用的2种修复算子分别为深度贪婪插入算子和随机插入算子,深度贪婪插入算子是将被移除的船舶按照一定顺序依次插入到它们的最佳位置,而此处的船舶顺序根据船舶的疏散时间值从小到大排序,随机插入算子是从被移除了 η \eta η艘船舶后剩余的解结构中随机选取 η \eta η个位置插入 η \eta η艘船舶;
- 初始化模拟退火温度T,以及初始化所有破坏算子di和修复算子ri的权重w(di),w(ri) 以及得分 s(di),s(ri),并计算各个算子对应的被选择概率p(di),p(ri);
- 通过轮盘赌的方式选择一个破坏算子和修复算子并更新算子被选择次数tm(di)和tm(ri),分别通过破坏算子和修复算子对当前解S进行操作得到新的解S’;
- 通过模拟退火保留机制决定是否接受新的解S’以及是否更新最优解S*,如果Ob(S’)<Ob(S),则接受该新的解;否则以概率e接受该新的解,然后更新破坏算子和修复算子对应的得分 s(di),s(ri) 以及模拟退火温度T;
- 达到一定代数pu之后,更新每个算子在下一个pu代的权重;
- 判断是否达到终止条件,如果是,则输出最优解;否则,返回步骤4。
Python代码实现部分
此处为跨市路径规划示例
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):
"""
定义目标函数计算距离
:param path: 已知访问顺序
:return: 求路径长度
"""
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 randomDestroy(sol):
"""
定义第一个摧毁算子-随机移除:随机移除3个城市
:param sol: 当前解
:return: 需要移除的城市
"""
solNew = copy.deepcopy(sol)
removed = []
removeIndex = rd.sample(range(0, distmat.shape[0]), 3) # sample: 随机选取3个城市存入,distmat.shape[0]为距离矩阵的维度
for i in removedIndex:
removed.append(solNew[i]) # 需要移除的城市添加到removed
sol.remove(solNew[i]) # 更新剩下的城市
return removed
def max3Destroy(sol):
"""
定义第二个摧毁算子-最大3段距离移除
:param sol: 当前解
:return: 需要移除的城市
"""
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])
'''
len(disSort)-i-1得到排在最后面即距离最长的路段。len-1,len-2,...
disSort[]得到最大值,dis.index得到最大值在距离列表中的索引
solNew得到最大值路段所在起点的索引,+1删除路段终点
如solNew = 【1,3,2】,dis = 【d13,d32,d21】=【9,7,8】
disSort = 【7,8,9】,先移除最后面的9,再移除8,...
9对应的dis索引为0,solnew对应索引0上的城市为1,移除1-3段,去掉3
'''
sol.remove(solNew[dis.index(disSort[len(disSort)-i-1])+1])) # 更新剩下的城市
return removed
def randomInsert(sol, removeList):
"""
定义第一个修复算子-随机插入
:param sol: 移除后的列表sol
:param removed: 移除列表
"""
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):
"""
定义第二个修复算子-贪婪插入(插入生成距离成本最小的解)
:param sol: 移除后的列表sol
:param removed: 移除列表
"""
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: # solNew作为temp,仅作为计算是否为更优解并记录索引位置
dis = disCal(solNew) # 更新距离,将原解换成新解
insertIndex = j
sol.insert(insertIndex, removeList[i]) # 更新列表
dis = float('inf') # 完成插入后需要初始化距离
def selectAndUseDestroyOperator(destroyWeight, currentSolution):
"""
定义破坏算子选择轮盘赌,其中需要调用摧毁算子randomDestroy() / max3Destroy()
:param destroyWeight : 摧毁算子权重
:param currentSolution: 当前解
:return: 需要移除的城市, 移除城市列表, 选择的摧毁算子序号
"""
des=troyOperator = -1 # 算子值初始化
sol = copy.deepcopy(currentSolution)
destroyRoulette = np.array(destroyWeight).cumsum() # 轮盘赌,cumsum()为依次累加权重值生成的新列表
r = rd.uniform(0, max(destroyRoulette)) # 随机生成权重树间任意的浮点值
for i in range(len(destroyWeight)): # 判断命中的是哪个算子
if destroyRoulette[i] >= r: # 判断是否在某个算子范围内,防止报错
if i == 0: # 序列为0,选择随机移除
destroyOperator = i
removedCities = randomDestroy(sol) # 得到随机移除的列表
destroyUseTimes[i] += 1 # 记录使用次数
break
elif i == 1: # 序列为1,选择最大3段距离移除
destroyOperator = i
removedCities = max3Destroy(sol)
destroyUseTimes[i] += 1 # 记录使用次数
break
return sol, removedCities, destroyOperator
def selectAndUseRepairOperator(repairWeight, destroyedSolution, removeList):
"""
定义修复算子选择轮盘赌,其中需要调用修复算子randomInsert()和greedyInsert()
:param repairWeight : 修复算子权重
:param destroyedSolution : 摧毁后的列表
:param removeList: 移除列表
:return: 修复解, 修复算子序号
"""
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: # 序列为0,选择随机插入
repairOperator = i
randomInsert(destroyedSolution, removeList)
repairUseTimes[i] += 1 # 记录使用次数
break
elif i == 1: # 序列为1,选择贪婪插入
repairOperator= i
greedyInsert(destroyedSolution, removeList)
repairUseTimes[i] += 1 # 记录使用次数
break
return destroyedSolution, repairOperator
'''
常数初始化
此处为枚举出各常数的初始值
'''
T = 100 # 模拟退火温度
a = 0.97 # 降火速度,越小降温越快。如同淬火,过快会导致算法来不及寻得最优
b = 0.5 # 权重更新系数,控制权重变化速度
wDestroy = [1 for i in range(2)] # 摧毁算子初始权重[1, 1]
wRepair = [1 for i in range(2)] # 修复算子初始权重[1, 1]
destroyUseTime = [0 for i range(2)] # 摧毁算子使用次数初始化
repairUseTime = [0 for i range(2)] # 修复算子使用次数初始化
destroyScore = [1 for i in range(2)] # 摧毁算子初始得分
repairScore = [1 for i in range(2)] # 修复算子初始得分
solution = [i for i in range(distmat.shape[0])] # 初始解
bestSolution = copy.deepcopy(solution) # 最优解
iterx = 0 # 初始迭代次数
iterMax = 100 # 最大迭代次数
'''
主程序
设置终止条件,内嵌选择轮盘赌及模拟退火算法接受准则
'''
if __name__ == '__main__':
while iterx < iterMax: # 终止条件,不满足则根据降火速度继续搜索
while T > 10: # 终止温度
destroyedSolution, remove, destroyOperatorIndex = selectAndUseDestroyOperator(wDestroy, solution)
repairOperatorIndex = selectAndUseRepairOperator(wRepair, destroyedSolution, remove)
if disCal(newSolution,) <= disCal(solution): # 比较新旧解的距离
solution = newSolution
if disCal(solution) <= disCal(bestSolution): # 新解<最优解
bestSolution = solution
destroyScore[destroyOperatorIndex] += 1.5
repairScore[destroyOperatorIndex] += 1.5
else:
destroyScore[destroyOperatorIndex] += 1.2
repairScore[destroyOperatorIndex] += 1.2
else:
if rd.random() < np.exp(disCal(newSolution)-disCal(solution)/T): # 使用模拟退火算法接受准则在一定标准下接受劣解
solution = newSolution
destroyScore[destroyOperatorIndex] += 0.8
repairScore[destroyOperatorIndex] += 0.8
else:
destroyScore[destroyOperatorIndex] += 0.6
repairScore[destroyOperatorIndex] += 0.6
wDestroy[destroyOperatorIndex] = (1-b)*(destroyScore[destroyOperatorIndex]/destroyUseTimes[destroyOperatorIndex]) + wDestroy[destroyOperatorIndex]*b # 更新权重
wRepair[repairOperatorIndex] = (1-b)*(repairScore[repairOperatorIndex]/repairUseTimes[repairOperatorIndex]) + wRepair[repairOperatorIndex]*b # 更新权重
T = a*T # 降温
iterx += 1 # 小于终止温度后, 增加一次迭代
T = 100 # 完成迭代后温度初始化
print(bestSolution) # [8, 11, 7, 10, 9, 0, 2, 1, 4, 3, 5, 6]
print(disCal(bestSolution)) # 4140
GO代码实现部分
结构体
/**
* 算法参数: 算法执行过程中一些初始化参数
*/
type Parameters struct {
MaxExecTime int json:"max_exec_time" // 最长执行时间(单位s)
TimeSegmentsIt int json:"time_segments_iteration" // 重新计算权重迭代次数
ReloadFrequency int json:"reload_frequency" // 重置当前解的迭代次数
Sigma1 int json:"sigma1" // 算子提升最优解 增加分数
Sigma2 int json:"sigma2" // 算子提升当前解 增加分数
Sigma3 int json:"sigma3" // 算子更新当前解 增加分数
Rho float64 json:"rho" // 重新计算算子权重的系数
MinimumWeight float64 json:"min_weight" // 最小权重值
MaximumWeight float64 json:"max_weight" // 最大权重值
MaxTemperature float64 json:"max_temperature" // 研究对象-最大温度
MinTemperature float64 json:"min_temperature" // 研究对象-最小温度
TemperatureAlpha float64 json:"temperature_alpha" // 研究对象-降温系数
MaxNoImproveRatio float64 json:"max_no_improve_ratio" // 最大没有提升解的迭代次数占比(超过停止)
}
/**
* 状态管理器: 管理计数的状态变量
*/
type Status struct {
IterationId int // 迭代次数
NIterationAvaliable int // 迭代次数中可行解的迭代次数
NIterationWithoutImprovement int // 距离上一次改善最优解的迭代次数
NIterationWithoutImprovementSinceLastReload int // 距离上一次重置当前解后改善最优解的迭代次数
NIterationWithoutImprovementCurrent int // 没有改善当前解的迭代次数
NIterationWithoutTransition int // 没有更新当前解的迭代次数
NIterationReloadCurrent int // 重置当前解的迭代次数
NIterationUpdateBest int // 更新最优解的迭代次数
NIterationRecomputeWeights int // 更新算子权重的迭代次数
NewBestSolution int // 当前解是否是最优解
AcceptedAsCurrentSolution int // 是否更新了当前解
ImproveCurrentSolution int // 是否提升了当前解
}
算子管理器
包含如下接口:
- 添加摧毁算子
- 添加修复算子
- 选择摧毁算子
- 选择修复算子
- 更新算子分数
- 更新算子权重
- 更新算子调用次数
/**
* 根据算子权重选择destory/repair
*/
func (m *OperatorManager) selectOperator(vecOp []IOperator, sumW float64) IOperator {
randomVal := rand.Float64()
randomWeightPos := randomVal * sumW
cumulSum := 0.0
for (i:=0; i<len(vecOp); i++ {
cumulSum += vecOp[i].GetWeight()
if cumulSum >= randomWeightPos {
if m.noise {
vecOp[i].SetNoise()
} else {
vecOp[i].UnsetNoise()
}
vecOp[i].IncreaseNumberOfCalls() // 更新算子使用次数
return vecOp[i]
}
}
return vecOp[len(vecOp)-1]
}
/**
* 获取邻域解,先destory再repair
*/
func (s *AlnsProcess) generateNeighborSolution(repair IRepairOperator, destory IDestoryOperator) *alg.Solution {
neighbor := s.currentSolution.CopySolution() // 从局部最优中生成新解
removeJobs := destory.DestorySolution(neighbor) // destory
repair.RepairSolution(neighbor, removeJobs) // repair
return neighbor
}
/**
* 更新当前解
* 新解 < 当前解,一定接受
* 新解 > 当前解,根据温度与成本(cost)变化值随机接受
*/
func (a *Acceptance) TransitionAccepted(currentSolution, newSolution *alg.Solution) bool {
a.temperature *= a.parameters.TemperatureAlpha // 降温
if newSolution.SolutionCost < currentSolution.SolutionCost {
return true
} else {
difference := newSolution.SolutionCost - currentSolution.SolutionCost
random := rand.Float64()
return math.Exp(-1.0*difference/a.temperature) >= random // 花费更少的新解算得的difference是负数
}
}
/**
* 更新最优解
* 新解 < 最优解,一定接受
*/
func (s *AlnsProcess) updateBestSolution(newSol *alg.Solution) bool {
bestSolution := s.bestSolManager.GetBestSolution()
accept := false
if newSol.SolutionCost < bestSolution.SolutionCost {
accept = true
}
if accept {
s.bestSolManager.AddNewBestSolution(newSol)
s.status.NewBestSolution = TRUE
s.status.NIterationWithoutImprovement = s.nIterationWithoutImprovement
s.status.NIterationWithoutImprovementSinceLastReload = 0
s.status.NIterationUpdateBest += 1
return true
} else {
s.status.NewBestSolution = FALSE
s.status.NIterationWithoutImprovement = s.nIterationWithoutImprovement
s.status.NIterationWithoutImprovementSinceLastReload++
return false
}
}
/**
* 更新算子的分数
*/
func (m *OperatorManager) UpateScores(des IDestoryOperator, rep IRepairOperator, status *Status) {
// 新解是最优解
if status.NewBestSolution == TRUE {
rep.SetScore(rep.GetScore() + float64(m.parameters.Sigma1))
des.SetScore(des.GetScore() + float64(m.parameters.Sigma1))
m.perfRepairOperatorsWithNoise += 1
m.perfRepairOperatorsWithoutNoise += 1
}
// 新解改善了当前解
if status.ImproveCurrentSolution == TRUE {
rep.SetScore(rep.GetScore() + float64(m.parameters.Sigma2))
des.SetScore(des.GetScore() + float64(m.parameters.Sigma2))
m.perfRepairOperatorsWithNoise += 1
m.perfRepairOperatorsWithoutNoise += 1
}
// 新解更新了当前解
if status.ImproveCurrentSolution == FALSE &&
status.AcceptedAsCurrentSolution == TRUE &&
status.AlreadyKnownSolution == FALSE {
rep.SetScore(rep.GetScore() + float64(m.parameters.Sigma3))
des.SetScore(des.GetScore() + float64(m.parameters.Sigma3))
m.perfRepairOperatorsWithNoise += 1
m.perfRepairOperatorsWithoutNoise += 1
}
if m.parameters.Noise {
performanceRepairOperatorsGlobal := 0.0
performanceRepairOperatorsGlobal += m.perfRepairOperatorsWithNoise
performanceRepairOperatorsGlobal += m.perfRepairOperatorsWithoutNoise
randomVal := rand.Float64()
randomWeightPos := randomVal * performanceRepairOperatorsGlobal
m.noise = (randomWeightPos < performanceRepairOperatorsGlobal)
}
}
/**
* 更新算子的权重
* 每迭代TimeSegmentsIt次,更新所有算子的权重,新的权重和算子分数、算子调用次数等有关
*/
func (m *OperatorManager) recomputeWeight(op IOperator, sumW *float64) {
prevWeight := op.GetWeight()
*sumW -= prevWeight
currentScore := op.GetScore()
nbCalls := op.GetNumberOfCallsSinceLastEvaluation()
newWeight := (1-m.parameters.Rho)*prevWeight + m.parameters.Rho*(float64(nbCalls)/float64(m.parameters.TimeSegmentsIt))*currentScore
// We ensure that the weight is within the bounds.
if newWeight > m.parameters.MaximumWeight {
newWeight = m.parameters.MaximumWeight
}
if newWeight < m.parameters.MinimumWeight {
newWeight = m.parameters.MinimumWeight
}
*sumW += newWeight
op.SetWeight(newWeight)
op.ResetScore()
op.ResetNumberOfCalls()
}
/**
* 重置当前解(防止陷入局部最优解)
* 每迭代 ReloadFrequency 次并且没有改善最优解,就重置当前解
*/
func (s *AlnsProcess) reloadCurrentSolution() *alg.Solution {
if s.status.NIterationWithoutImprovementSinceLastReload > 0 &&
s.status.NIterationWithoutImprovementSinceLastReload%s.params.ReloadFrequency == 0 {
s.status.NIterationWithoutImprovementSinceLastReload = 0
s.status.NIterationReloadCurrent += 1
return s.bestSolManager.GetBestSolution().CopySolution()
}
return s.currentSolution
}
参考文献及算法应用:
[1]王新. 车辆和无人机联合配送路径问题研究[D].大连海事大学,2020.
[2]李婷玉. 多商户多车程同城物流配送车辆调度问题研究[D].大连理工大学,2018.
[3]张梦颖. 不确定因素下路径规划问题研究[D].中国科学技术大学,2016.
[4]蔡文学. 台风天气下广州港出海航道船舶疏散调度研究[D].华南理工大学,2021.