一.问题描述
考虑带容量限制的VRP问题，即CVRP。
1.1 假设
配送中心：只有一个配送中心；
配送方式：任何车辆配送货物结束后需要返回配送中心；
车型：只考虑一种车型，且容量约束相同；
需求不可拆分：客户需求只能由一辆车满足，即任一客户需求小于车辆容量；
车辆充足：不限制车辆数量，即配送车辆需求均能满足；
1.2 要求
优化目标：最小化车辆启动成本和车辆行驶成本之和；
约束条件：客户访问约束，容积约束；
已知信息：配送中心位置、客户点位置、客户点需求、车辆最大容积、车辆启动成本、车辆单位距离行驶成本；（数据如有需要请后台私信我）
二.模拟退火算法（SA)
2.1 SA算法原理
模拟退火算法包含两个部分即Metropolis算法和退火过程,，分别对应内循环和外循环。外循环就是退火过程，将固体达到较高的温度（初始温度T(0))，然后按照降温系数alpha使温度按照一定的比例下降，当达到终止温度T(f)时，冷却结束，即算法过程结束。
Metropolis算法是内循环，即在每次温度下，迭代L次，寻找在该温度下能量的最小值（即最优解）。在一次温度下，整个迭代过程中温度不发生变化，能量发生变化，如果当前一个状态x(n)的能量大于后一个状态x(n+1)的能量时，即状态x(n)的解没有状态x(n+1)的解好，所以接受状态x(n+1)。但是如果下一状态的能量比前一个状态的能量高时，该不该接受下一状态呢？在这里设置一个接受概率P，计算公式如下：
其中能量的变化量和T进行决定概率P的大小，所以这个值是动态的。然后进行概率操作，在区间【0,1】产生一个均匀分布的随机数ϵ ，如果ϵ<P，则此种转移接受，否则拒绝转移，进入下一步，往复循环。
2.2 SA算法流程
- 初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起点)，每个T值的迭代次数L
- 对k=1, …, L做第(3)至第6步：
- 产生新解S′，通常选择由当前解经过简单变换产生新解S′，例如对构成解的全部或部分元素进行置换、互换等。产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。
- 计算增量ΔT=C(S′)-C(S)，其中C(S)为目标函数,C(S)相当于能量
- 若ΔT<0则接受S′作为新的当前解，否则以概率exp(-ΔT/T)接受S′作为新的当前解.
- 如果满足终止条件则输出当前解作为最优解，结束程序。
- T逐渐减少，且T->0，然后转第2步。
三.Python算法求解
step.0 导入模块
import pandas as pd
import math
import random
import numpy as np
import copy
import xlsxwriter
import matplotlib.pyplot as pltstep.1 定义解、网络节点、全局参数的数据结构
# 数据结构：解
class Sol():
def __init__(self):
self.nodes_seq=None # 解的编码
self.obj=None # 目标函数
self.routes=None # 解的解码
# 数据结构：网络节点
class Node():
def __init__(self):
self.id=0 # 节点id
self.name='' # 节点名称，可选
self.seq_no=0 # 节点映射id
self.x_coord=0 # 节点平面横坐标
self.y_coord=0 # 节点平面纵坐标
self.demand=0 # 节点需求
# 数据结构：全局参数
class Model():
def __init__(self):
self.best_sol=None # 全局最优解
self.node_list=[] # 需求节点集合
self.node_seq_no_list=[] #需求节点映射id集合
self.depot=None # 车场节点
self.number_of_nodes=0 # 需求节点数量
self.opt_type=0 # 优化目标类型
self.vehicle_cap=0 # 车辆最大容量step.2处理表格数据
def readXlsxFile(filepath,model):
# 建议在xlsx文件中，第一行为表头，其中: x_coord,y_coord,demand是必须项；车辆基地数据放在表头下首行
node_seq_no = -1#车辆基地的seq_no值为-1,剩余需求节点的seq_no 依次编号为 0,1,2,...
df = pd.read_excel(filepath)
for i in range(df.shape[0]):
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)step.4构造初始解（基于Split思想）
def genInitialSol(node_seq):
node_seq=copy.deepcopy(node_seq)
random.seed(0)
random.shuffle(node_seq)
return node_seq
# 采用Split思想对TSP序列进行切割，得到可行车辆路径
def splitRoutes(nodes_seq,model):
"""
采用简单的分割方法：按顺序依次检查路径的容量约束，在超出车辆容量限制的位置插入车场。
例如某TSP解为：[1,2,3,4,5,6,7,8,9,10],累计需求为：[10,20,30,40,50,60,70,80,90,10]，车辆容量为：30，则应在3,6,9节点后插入车场，
即得到：[0,1,2,3,0,4,5,6,0,7,8,9,0,10,0]
"""
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_routesstep.5计算行驶距离以及目标函数
# 计算路径总行驶距离
def calDistance(route,model):
distance=0
depot=model.depot
for i in range(len(route)-1):
from_node=model.node_list[route[i]]
to_node=model.node_list[route[i+1]]
distance+=math.sqrt((from_node.x_coord-to_node.x_coord)**2+(from_node.y_coord-to_node.y_coord)**2)
first_node=model.node_list[route[0]]
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_routesstep.6构造新解
# 定义邻域算子
def createActions(n):
action_list=[]
nswap=n//2
# 第一种算子（Swap）：前半段与后半段对应位置一对一交换
for i in range(nswap):
action_list.append([1, i, i + nswap])
# 第二种算子（DSwap）：前半段与后半段对应位置二对二交换
for i in range(0, nswap, 2):
action_list.append([2, i, i + nswap])
# 第三种算子（Reverse）：指定长度的序列反序
for i in range(0, n, 4):
action_list.append([3, i, i + 3])
return action_list
# 执行邻域搜索
def doACtion(nodes_seq,action):
nodes_seq=copy.deepcopy(nodes_seq)
if action[0]==1:
index_1=action[1]
index_2=action[2]
temporary=nodes_seq[index_1]
nodes_seq[index_1]=nodes_seq[index_2]
nodes_seq[index_2]=temporary
return nodes_seq
elif action[0]==2:
index_1 = action[1]
index_2 = action[2]
temporary=[nodes_seq[index_1],nodes_seq[index_1+1]]
nodes_seq[index_1]=nodes_seq[index_2]
nodes_seq[index_1+1]=nodes_seq[index_2+1]
nodes_seq[index_2]=temporary[0]
nodes_seq[index_2+1]=temporary[1]
return nodes_seq
elif action[0]==3:
index_1=action[1]
index_2=action[2]
nodes_seq[index_1:index_2+1]=list(reversed(nodes_seq[index_1:index_2+1]))
return nodes_seqstep.7收敛曲线与车辆路径
# 绘制目标函数收敛曲线
def plotObj(obj_list):
plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
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()
# 绘制优化车辆路径
def plotRoutes(model):
for route in model.best_sol.routes:
x_coord = [model.depot.x_coord]
y_coord = [model.depot.y_coord]
for node_no in route:
x_coord.append(model.node_list[node_no].x_coord)
y_coord.append(model.node_list[node_no].y_coord)
x_coord.append(model.depot.x_coord)
y_coord.append(model.depot.y_coord)
plt.plot(x_coord, y_coord, marker='s', color='b', linewidth=0.5)
plt.show()step.8定义参数、输出结果
# 输出优化结果
def outPut(model):
work = xlsxwriter.Workbook('result.xlsx')
worksheet = work.add_worksheet()
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):
route.insert(0, model.depot.id)
route.append(model.depot.id)
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()
def run(filepath,T0,Tf,deltaT,v_cap,opt_type):
"""
:param filepath: xlsx文件路径
:param T0: 初始温度
:param Tf: 终止温度
:param deltaT: 温度下降步长或下降比例
:param v_cap: 车辆容量
:param opt_type: 优化类型:0:最小化车辆数,1:最小化行驶距离
:return:
"""
model=Model()
model.vehicle_cap=v_cap
model.opt_type=opt_type
readXlsxFile(filepath,model)
action_list=createActions(model.number_of_nodes)
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)
Tk=T0
nTk=len(action_list)
while Tk>=Tf:
for i in range(nTk):
new_sol = Sol()
new_sol.nodes_seq = doACtion(sol.nodes_seq, action_list[i])
new_sol.obj, new_sol.routes = calObj(new_sol.nodes_seq, model)
detla_f=new_sol.obj-sol.obj
if detla_f<0 or math.exp(-detla_f/Tk)>random.random():
sol=copy.deepcopy(new_sol)
if sol.obj<model.best_sol.obj:
model.best_sol=copy.deepcopy(sol)
if deltaT<1:
Tk=Tk*deltaT
else:
Tk = Tk - deltaT
history_best_obj.append(model.best_sol.obj)
print("当前温度：%s，local obj:%s best obj: %s" % (Tk,sol.obj,model.best_sol.obj))
plotObj(history_best_obj)
plotRoutes(model)
outPut(model)
if __name__=='__main__':
file=r'filepath'
run(filepath=file,T0=6000,Tf=0.001,deltaT=0.95,v_cap=80,opt_type=1)最终结果如下：
扫一扫
1901
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
