# Python implements a common VRP solution algorithm -- adaptive large neighborhood search algorithm (ALNS)

Keywords: Python Algorithm

Reference notes
https://github.com/PariseC/Algorithms_for_solving_VRP

## 1. Applicable scenarios

• Solve CVRP
• Single vehicle type
• The vehicle capacity shall not be less than the maximum demand of the demand node
• Single vehicle base

# 2. Problem analysis

The solution of CVRP problem is a set of paths of multiple vehicles that meet the needs of demand nodes.

Assuming that there are 10 customer nodes in a physical network, numbered 1 ~ 10, and a vehicle base, numbered 0, a feasible solution to this problem can be expressed as: [0-1-2-0,0-3-4-5-0,0-6-7-8-0,0-9-10-0], that is, four vehicles are required to provide services, and the driving routes of vehicles are 0-1-2-0,0-3-4-5-0,0-6-7-8-0 respectively, 0-9-10-0.

Since the capacity of vehicles is fixed and the base is fixed, the solution of the above problem can be expressed as an ordered sequence of [1-2-3-4-5-6-7-8-9-10], and then the sequence can be cut according to the capacity constraints of vehicles to obtain the driving routes of several vehicles. Therefore, the CVRP problem can be transformed into a TSP problem for solution. After obtaining the optimal solution of the TSP problem, the path cutting is carried out considering the vehicle capacity constraint to obtain the solution of the CVRP problem. Such a treatment may affect the quality of the solution of the CVRP problem, but it simplifies the difficulty of solving the problem.

# 3. Data format

The network data is stored in xlsx file, in which the first line is the title bar and the second line stores the depot data. In the program, vehicle base seq_ The no number is - 1, and the demand node seq_id is numbered from 0. Refer to the relevant documents on the github home page. # 4. Step by step implementation

## (1) Data structure

To facilitate data processing, Sol() class, Node() class and Model() class are defined. Their properties are shown in the following table:

• Sol() class, representing a feasible solution
attributedescribe
nodes_seqDemand node seq_no ordered permutation set, corresponding to the solution of TSP
objOptimization target value
routesVehicle path set, corresponding to the solution of CVRP
• Node() class, representing a network node
attributedescribe
idPhysical node id, optional
namePhysical node name, optional
seq_noPhysical node mapping id, base node is - 1, and demand node is numbered from 0
x_coordPhysical node x coordinate
y_coordPhysical node y coordinate
demandPhysical node requirements
• The Model() class stores algorithm parameters
attributedescribe
best_solGlobal optimal solution, value type is Sol()
node_listPhysical node collection with value type Node()
node_seq_no_listSet of physical node mapping IDS
depotDepot, the value type is Node()
number_of_nodesNumber of demand nodes
opt_typeOptimization target type, 0: minimum number of vehicles, 1: minimum driving distance
vehicle_capVehicle capacity
distanceNetwork arc distance
rand_d_maxUpper limit of random failure degree
rand_d_minLower limit of random failure degree
worst_d_maxUpper limit of worst damage
worst_d_minLower limit of worst failure degree
regret_nNumber of suboptimal positions
r1Operator reward 1
r2Operator reward 2
r3Operator reward 3
rhoOperator weight attenuation coefficient
d_weightDestruction operator weight
d_selectNumber of times destruction operator is selected / round
d_scoreDamage operator is rewarded with score / round
d_history_selectTotal number of times destruction operator history is selected
d_history_scoreTotal awarded score of destruction operator history
r_weightRepair operator weight
r_selectRepair operator selected times / round
r_scoreThe repair operator is rewarded with a score / round
r_history_selectTotal number of times repair operator history is selected
r_history_scoreTotal awarded score of repair operator history

```def readXlsxFile(filepath,model):
# It is recommended that the vehicle depot data be placed in the first line of xlsx file
node_seq_no = -1#the depot node seq_no is -1,and demand node seq_no is 0,1,2,...
for i in range(df.shape):
node=Node()
node.id=node_seq_no
node.seq_no=node_seq_no
node.x_coord= df['x_coord'][i]
node.y_coord= df['y_coord'][i]
node.demand=df['demand'][i]
if df['demand'][i] == 0:
model.depot=node
else:
model.node_list.append(node)
model.node_seq_no_list.append(node_seq_no)
try:
node.name=df['name'][i]
except:
pass
try:
node.id=df['id'][i]
except:
pass
node_seq_no=node_seq_no+1
model.number_of_nodes=len(model.node_list)
```

## (3) Initialization parameters

The network arc distance is calculated when initializing parameters for use by subsequent operators.

```def initParam(model):
for i in range(model.number_of_nodes):
for j in range(i+1,model.number_of_nodes):
d=math.sqrt((model.node_list[i].x_coord-model.node_list[j].x_coord)**2+
(model.node_list[i].y_coord-model.node_list[j].y_coord)**2)
model.distance[i,j]=d
model.distance[j,i]=d
```

## (4) Initial solution

```def genInitialSol(node_seq):
node_seq=copy.deepcopy(node_seq)
random.seed(0)
random.shuffle(node_seq)
return node_seq
```

## (5) Target value calculation

The target value calculation depends on the "splitRoutes" function to segment the TSP feasible solution to obtain the vehicle driving route and the number of vehicles required, and the "calDistance" function calculates the driving distance.

```def splitRoutes(nodes_seq,model):
num_vehicle = 0
vehicle_routes = []
route = []
remained_cap = model.vehicle_cap
for node_no in nodes_seq:
if remained_cap - model.node_list[node_no].demand >= 0:
route.append(node_no)
remained_cap = remained_cap - model.node_list[node_no].demand
else:
vehicle_routes.append(route)
route = [node_no]
num_vehicle = num_vehicle + 1
remained_cap =model.vehicle_cap - model.node_list[node_no].demand
vehicle_routes.append(route)
return num_vehicle,vehicle_routes
def calDistance(route,model):
distance=0
depot=model.depot
for i in range(len(route)-1):
distance+=model.distance[route[i],route[i+1]]
first_node=model.node_list[route]
last_node=model.node_list[route[-1]]
distance+=math.sqrt((depot.x_coord-first_node.x_coord)**2+(depot.y_coord-first_node.y_coord)**2)
distance+=math.sqrt((depot.x_coord-last_node.x_coord)**2+(depot.y_coord - last_node.y_coord)**2)
return distance
def calObj(nodes_seq,model):
num_vehicle, vehicle_routes = splitRoutes(nodes_seq, model)
if model.opt_type==0:
return num_vehicle,vehicle_routes
else:
distance=0
for route in vehicle_routes:
distance+=calDistance(route,model)
return distance,vehicle_routes
```

## (6) Define destroy operator (destroy operator)

Two types of destruction are implemented here:

• random destroy: randomly remove a certain proportion of demand nodes from the current solution
• worst destroy: remove a certain proportion of demand nodes from the current solution that cause a large increase in the objective function
```# define random destory action
def createRandomDestory(model):
d=random.uniform(model.rand_d_min,model.rand_d_max)
reomve_list=random.sample(range(model.number_of_nodes),int(d*model.number_of_nodes))
return reomve_list
# define worse destory action
def createWorseDestory(model,sol):
deta_f=[]
for node_no in sol.nodes_seq:
nodes_seq_=copy.deepcopy(sol.nodes_seq)
nodes_seq_.remove(node_no)
obj,vehicle_routes=calObj(nodes_seq_,model)
deta_f.append(sol.obj-obj)
sorted_id = sorted(range(len(deta_f)), key=lambda k: deta_f[k], reverse=True)
d=random.randint(model.worst_d_min,model.worst_d_max)
remove_list=sorted_id[:d]
return remove_list
```

## (7) Define repair operator (repair operator)

Three types of repair are implemented here:

• random repair: randomly insert the removed demand node into the allocated node sequence;
• Green repair: according to the increment size of the objective function of each possible position in the sequence of Allocated Nodes inserted by the removed demand node, select the demand node with the smallest increment of the objective function and the insertion position in turn until all the removed demand nodes are re inserted (which can be simply understood as, Select the demand node that minimizes the increment of the objective function and its optimal insertion position in turn); Among them, P P P is the unassigned demand node set, I P IP IP is a collection of locations that may be inserted, f ( s p , i ) f(s_{p,i}) f(sp,i) is the demand node p p p insert s s s i i Objective function value at i positions.
• regret repair: calculate the sum of the difference between the objective function value and the objective function value of the optimal position when the removed node is inserted back to n suboptimal positions in the allocated node sequence, and then select the demand node with the largest difference and its optimal position( It can be simply understood that n suboptimal locations and demand nodes far from the optimal location and their optimal locations are preferred.); Among them, P P P is the unassigned demand node set, f ( s i ( p ) ) − f ( s 1 ( p ) ) f(s_{i}(p))-f(s_{1}(p)) f(si (p)) − f(s1 (p)) is the demand node p p p insert s s The optimal position in s and i i i difference of objective function at suboptimal position.
```def createRandomRepair(remove_list,model,sol):
unassigned_nodes_seq=[]
assigned_nodes_seq = []
# remove node from current solution
for i in range(model.number_of_nodes):
if i in remove_list:
unassigned_nodes_seq.append(sol.nodes_seq[i])
else:
assigned_nodes_seq.append(sol.nodes_seq[i])
# insert
for node_no in unassigned_nodes_seq:
index=random.randint(0,len(assigned_nodes_seq)-1)
assigned_nodes_seq.insert(index,node_no)
new_sol=Sol()
new_sol.nodes_seq=copy.deepcopy(assigned_nodes_seq)
new_sol.obj,new_sol.routes=calObj(assigned_nodes_seq,model)
return new_sol

def createGreedyRepair(remove_list,model,sol):
unassigned_nodes_seq = []
assigned_nodes_seq = []
# remove node from current solution
for i in range(model.number_of_nodes):
if i in remove_list:
unassigned_nodes_seq.append(sol.nodes_seq[i])
else:
assigned_nodes_seq.append(sol.nodes_seq[i])
#insert
while len(unassigned_nodes_seq)>0:
insert_node_no,insert_index=findGreedyInsert(unassigned_nodes_seq,assigned_nodes_seq,model)
assigned_nodes_seq.insert(insert_index,insert_node_no)
unassigned_nodes_seq.remove(insert_node_no)
new_sol=Sol()
new_sol.nodes_seq=copy.deepcopy(assigned_nodes_seq)
new_sol.obj,new_sol.routes=calObj(assigned_nodes_seq,model)
return new_sol

def findGreedyInsert(unassigned_nodes_seq,assigned_nodes_seq,model):
best_insert_node_no=None
best_insert_index = None
best_insert_cost = float('inf')
assigned_nodes_seq_obj,_=calObj(assigned_nodes_seq,model)
for node_no in unassigned_nodes_seq:
for i in range(len(assigned_nodes_seq)):
assigned_nodes_seq_ = copy.deepcopy(assigned_nodes_seq)
assigned_nodes_seq_.insert(i, node_no)
obj_, _ = calObj(assigned_nodes_seq_, model)
deta_f = obj_ - assigned_nodes_seq_obj
if deta_f<best_insert_cost:
best_insert_index=i
best_insert_node_no=node_no
best_insert_cost=deta_f
return best_insert_node_no,best_insert_index

def createRegretRepair(remove_list,model,sol):
unassigned_nodes_seq = []
assigned_nodes_seq = []
# remove node from current solution
for i in range(model.number_of_nodes):
if i in remove_list:
unassigned_nodes_seq.append(sol.nodes_seq[i])
else:
assigned_nodes_seq.append(sol.nodes_seq[i])
# insert
while len(unassigned_nodes_seq)>0:
insert_node_no,insert_index=findRegretInsert(unassigned_nodes_seq,assigned_nodes_seq,model)
assigned_nodes_seq.insert(insert_index,insert_node_no)
unassigned_nodes_seq.remove(insert_node_no)
new_sol = Sol()
new_sol.nodes_seq = copy.deepcopy(assigned_nodes_seq)
new_sol.obj, new_sol.routes = calObj(assigned_nodes_seq, model)
return new_sol

def findRegretInsert(unassigned_nodes_seq,assigned_nodes_seq,model):
opt_insert_node_no = None
opt_insert_index = None
opt_insert_cost = -float('inf')
for node_no in unassigned_nodes_seq:
n_insert_cost=np.zeros((len(assigned_nodes_seq),3))
for i in range(len(assigned_nodes_seq)):
assigned_nodes_seq_=copy.deepcopy(assigned_nodes_seq)
assigned_nodes_seq_.insert(i,node_no)
obj_,_=calObj(assigned_nodes_seq_,model)
n_insert_cost[i,0]=node_no
n_insert_cost[i,1]=i
n_insert_cost[i,2]=obj_
n_insert_cost= n_insert_cost[n_insert_cost[:, 2].argsort()]
deta_f=0
for i in range(1,model.regret_n):
deta_f=deta_f+n_insert_cost[i,2]-n_insert_cost[0,2]
if deta_f>opt_insert_cost:
opt_insert_node_no = int(n_insert_cost[0, 0])
opt_insert_index=int(n_insert_cost[0,1])
opt_insert_cost=deta_f
return opt_insert_node_no,opt_insert_index
```

## (8) Random selection operator

The roulette method is used to select destruction and repair operators.

```#select destory action and repair action
def selectDestoryRepair(model):
d_weight=model.d_weight
d_cumsumprob = (d_weight / sum(d_weight)).cumsum()
d_cumsumprob -= np.random.rand()
destory_id= list(d_cumsumprob > 0).index(True)

r_weight=model.r_weight
r_cumsumprob = (r_weight / sum(r_weight)).cumsum()
r_cumsumprob -= np.random.rand()
repair_id = list(r_cumsumprob > 0).index(True)
return destory_id,repair_id
```

## (9) Execute destroy operator

The index sequence of nodes to be removed is returned according to the selected destruction operator.

```# do destory action
def doDestory(destory_id,model,sol):
if destory_id==0:
reomve_list=createRandomDestory(model)
else:
reomve_list=createWorseDestory(model,sol)
return reomve_list
```

## (10) Execute repair operator

Repair the current connection according to the selected repair operator.

```def doRepair(repair_id,reomve_list,model,sol):
if repair_id==0:
new_sol=createRandomRepair(reomve_list,model,sol)
elif repair_id==1:
new_sol=createGreedyRepair(reomve_list,model,sol)
else:
new_sol=createRegretRepair(reomve_list,model,sol)
return new_sol
```

## (11) Reset operator score

Every pu time destruction and repair are executed, the score of the reset operator and the number of times it is selected.

```# reset action score
def resetScore(model):

model.d_select = np.zeros(2)
model.d_score = np.zeros(2)

model.r_select = np.zeros(3)
model.r_score = np.zeros(3)
```

## (12) Update operator weight

There are two strategies for updating operator weights. One is to update the weights successively every pu time destruction and repair are executed, and the other is to update the weights every time destruction and repair are executed. The former can ensure that the weight is updated in time, but it needs more calculation time; In the latter, by reasonably setting pu parameters, the calculation time is saved, and the weight update does not lag too much. The latter update strategy is adopted here. Among them, w ( h ) w(h) w(h) represents the operator weight, ρ ρ ρ Is the attenuation coefficient, s ( h ) s(h) s(h) is the total reward score of the operator, u ( h ) u(h) u(h) is the number of times the operator is selected.

```# update action weight
def updateWeight(model):
for i in range(model.d_weight.shape):
if model.d_select[i]>0:
model.d_weight[i]=model.d_weight[i]*(1-model.rho)+model.rho*model.d_score[i]/model.d_select[i]
else:
model.d_weight[i] = model.d_weight[i] * (1 - model.rho)
for i in range(model.r_weight.shape):
if model.r_select[i]>0:
model.r_weight[i]=model.r_weight[i]*(1-model.rho)+model.rho*model.r_score[i]/model.r_select[i]
else:
model.r_weight[i] = model.r_weight[i] * (1 - model.rho)
model.d_history_select = model.d_history_select + model.d_select
model.d_history_score = model.d_history_score + model.d_score
model.r_history_select = model.r_history_select + model.r_select
model.r_history_score = model.r_history_score + model.r_score
```

## (13) Draw convergence curve

```def plotObj(obj_list):
plt.rcParams['font.sans-serif'] = ['SimHei']  #show chinese
plt.rcParams['axes.unicode_minus'] = False   # Show minus sign
plt.plot(np.arange(1,len(obj_list)+1),obj_list)
plt.xlabel('Iterations')
plt.ylabel('Obj Value')
plt.grid()
plt.xlim(1,len(obj_list)+1)
plt.show()
```

## (14) Output results

```def outPut(model):
work=xlsxwriter.Workbook('result.xlsx')
worksheet.write(0,0,'opt_type')
worksheet.write(1,0,'obj')
if model.opt_type==0:
worksheet.write(0,1,'number of vehicles')
else:
worksheet.write(0, 1, 'drive distance of vehicles')
worksheet.write(1,1,model.best_sol.obj)
for row,route in enumerate(model.best_sol.routes):
worksheet.write(row+2,0,'v'+str(row+1))
r=[str(i)for i in route]
worksheet.write(row+2,1, '-'.join(r))
work.close()
```

(15) Main function
Two layers of loops are set in the main function. The outer loop is controlled by the epochs parameter and the inner loop is controlled by the pu parameter. The operator weight is updated every time the inner loop is executed. There are many types of acceptance criteria for new neighborhood solutions, such as RW, GRE, SA, Ta, oba and GDA. Here, the TA criterion is adopted, i.e f ( n e w x ) − f ( x ) < T f(newx)-f(x)<T Accept the new solution when f(newx) − f (x) < t, T T T according to coefficient φ φ φ Attenuation, initial T T T is the current solution * 0.2. In addition, when defining the destruction operator, there is no specific damage degree value specified here, but an interval is used to randomly select a damage degree each time to destroy the current solution. For more algorithm details, please refer to reference 1 at the end of the article. For acceleration techniques, please refer to reference 4 and 5 at the end of the article.

```def run(filepath,rand_d_max,rand_d_min,worst_d_min,worst_d_max,regret_n,r1,r2,r3,rho,phi,epochs,pu,v_cap,opt_type):
"""
:param filepath: Xlsx file path
:param rand_d_max: max degree of random destruction
:param rand_d_min: min degree of random destruction
:param worst_d_max: max degree of worst destruction
:param worst_d_min: min degree of worst destruction
:param regret_n:  n next cheapest insertions
:param r1: score if the new solution is the best one found so far.
:param r2: score if the new solution improves the current solution.
:param r3: score if the new solution does not improve the current solution, but is accepted.
:param rho: reaction factor of action weight
:param phi: the reduction factor of threshold
:param epochs: Iterations
:param pu: the frequency of weight adjustment
:param v_cap: Vehicle capacity
:param opt_type: Optimization type:0:Minimize the number of vehicles,1:Minimize travel distance
:return:
"""
model=Model()
model.rand_d_max=rand_d_max
model.rand_d_min=rand_d_min
model.worst_d_min=worst_d_min
model.worst_d_max=worst_d_max
model.regret_n=regret_n
model.r1=r1
model.r2=r2
model.r3=r3
model.rho=rho
model.vehicle_cap=v_cap
model.opt_type=opt_type
initParam(model)
history_best_obj = []
sol = Sol()
sol.nodes_seq = genInitialSol(model.node_seq_no_list)
sol.obj, sol.routes = calObj(sol.nodes_seq, model)
model.best_sol = copy.deepcopy(sol)
history_best_obj.append(sol.obj)
for ep in range(epochs):
T=sol.obj*0.2
resetScore(model)
for k in range(pu):
destory_id,repair_id=selectDestoryRepair(model)
model.d_select[destory_id]+=1
model.r_select[repair_id]+=1
reomve_list=doDestory(destory_id,model,sol)
new_sol=doRepair(repair_id,reomve_list,model,sol)
if new_sol.obj<sol.obj:
sol=copy.deepcopy(new_sol)
if new_sol.obj<model.best_sol.obj:
model.best_sol=copy.deepcopy(new_sol)
model.d_score[destory_id]+=model.r1
model.r_score[repair_id]+=model.r1
else:
model.d_score[destory_id]+=model.r2
model.r_score[repair_id]+=model.r2
elif new_sol.obj-sol.obj<T:
sol=copy.deepcopy(new_sol)
model.d_score[destory_id] += model.r3
model.r_score[repair_id] += model.r3
T=T*phi
print("%s/%s:%s/%s， best obj: %s" % (ep,epochs,k,pu, model.best_sol.obj))
history_best_obj.append(model.best_sol.obj)
updateWeight(model)

plotObj(history_best_obj)
outPut(model)
print("random destory weight is {:.3f}\tselect is {}\tscore is {:.3f}".format(model.d_weight,
model.d_history_select,
model.d_history_score))
print("worse destory weight is {:.3f}\tselect is {}\tscore is {:.3f} ".format(model.d_weight,
model.d_history_select,
model.d_history_score))
print("random repair weight is {:.3f}\tselect is {}\tscore is {:.3f}".format(model.r_weight,
model.r_history_select,
model.r_history_score))
print("greedy repair weight is {:.3f}\tselect is {}\tscore is {:.3f}".format(model.r_weight,
model.r_history_select,
model.r_history_score))
print("regret repair weight is {:.3f}\tselect is {}\tscore is {:.3f}".format(model.r_weight,
model.r_history_select,
model.r_history_score))
```

reference resources

• https://oparu.uniulm.de/xmlui/bitstream/handle/123456789/3264/vts_9568_14474.pdf?sequence=1
• https://homes.di.unimi.it/righini/Didattica/AlgoritmiEuristici/MaterialeAE/15%20-%20ALNS.pdf
• https://mp.weixin.qq.com/s/fpeNCKo62d_7gqBnWd2tTw
• https://mp.weixin.qq.com/s/l7fTHRye4pPBtWrlPrZ6mA
• https://mp.weixin.qq.com/s/K6EgtMdarCgr_vK4qZfGlw https://github.com/PariseC/Algorithms_for_solving_VRP

Posted by hwttdz on Thu, 02 Sep 2021 11:21:44 -0700