引言
本系列文章是路径优化问题学习过程中一个完整的学习路线。问题从简单的单车场容量约束CVRP问题到多车场容量约束MDCVRP问题,再到多车场容量+时间窗口复杂约束MDCVRPTW问题,复杂度是逐渐提升的。
如果大家想学习某一个算法,建议从最简单的TSP问题开始学习,一个阶梯一个阶梯走。
如果不是奔着学习算法源码的思路,只想求解出个结果,请看历史文章,有ortools、geatpy、scikit-opt等求解器相关文章,点击→路径优化历史文章,可直接跳转。
本系列文章汇总:
- 1、路径优化历史文章
- 2、路径优化丨带时间窗和载重约束的CVRPTW问题-改进遗传算法:算例RC108
- 3、路径优化丨带时间窗和载重约束的CVRPTW问题-改进和声搜索算法:算例RC108
- 4、路径优化丨复现论文-网约拼车出行的乘客车辆匹配及路径优化
- 5、多车场路径优化丨遗传算法求解MDCVRP问题
- 6、论文复现详解丨多车辆路径优化问题:粒子群+模拟退火算法求解
问题描述
多车场车辆路径问题 (MDVRP)是基本车辆路径问题(VRP)的扩展,指的是有数个车场同时对多个用户进行服务,各用户有一定的货物需求,每个车场都可提供货物,并有车队负责执行运输任务,要求对各车场的车辆和行驶路线进行适当的安排,在保证满足各用户的需求的前提下,使总的运输成本最低。
问题满足假设:
- 各车场派出车辆的数目不能超过该车场所拥有的车辆数;
- 各车辆都是从各自的车场出发,并回到原车场;
- 每个用户只能被一辆车服务一次;
- 车辆满足载重约束。
数学模型
数学模型:
运行环境
windows系统,python3.6.0,第三方库及版本号如下:
numpy==1.18.5
matplotlib==3.2.1
xlsxwriter==1.4.3
xlrd==1.2.0
第三方库需要在安装完python之后,额外安装,以前文章有讲述过安装第三方库的解决办法。
算法数据处理
1、数据:
本文是数据是3个车场和100个需求量,各个车场的车辆数是20,具体车场和需求点的坐标及需求点的需求量见两个depot和demand文件,车辆载重暂时设为50。模型是一般的路径优化模型,目标函数设置为3个车场派出的所有车辆的路径总长度。
以3个车场为原点,分别和100个需求点组合,计算组合后的任意两点的距离,形成3个101×101的距离矩阵。距离矩阵计算比较简单,计算矩阵右上角的值即可,左下角的值和右上角的值对称。
2、距离矩阵
def distance(dep_demand): # 3个点和需求点距离矩阵计算
dis = np.zeros((len(dep_demand),len(dep_demand)))
for i in range(len(dep_demand)-1):
for j in range(i+1):
point1 = dep_demand[i]
point2 = dep_demand[j]
di = np.sqrt((point1[0]-point2[0])**2+(point1[1]-point2[1])**2)
dis[i,j], dis[j,i] = di, di # 距离关于对角线对称
return dis
编码及解码
1、编码
采用实数编码的方式,用1到100共100个数表示100个需求点,随机打乱即是一个可行的编码:
road = random.sample(range(1,101),100)
一个可行编码:
road = [95, 34, 71, 23, 18, 17, 80, 91, 96, 60, 50, 28, 72, 37, 30, 41, 33, 2, 15, 66, 36, 73, 38, 82, 49, 84, 31, 16, 3, 7, 43, 6, 4, 20, 13, 5, 93, 52, 70, 53, 29, 89, 67, 98, 26, 45, 100, 79, 87, 64, 59, 81, 44, 61, 63, 83, 8, 76, 92, 12, 32, 1, 9, 77, 75, 46, 65, 25, 85, 24, 69, 86, 58, 19, 22, 74, 88, 90, 94, 21, 55, 39, 42, 51, 27, 10, 48, 35, 97, 57, 11, 56, 40, 54, 62, 99, 68, 78, 14, 47]
2、解码
核心是遍历到某个需求点时,寻找最优车场的最优车辆。一个车场的最优车辆寻找如下:
step1:遍历该车场的所有车辆,如果至少有一辆车的剩余载重满足该需求点,计算对应车辆的最后一个点到该需求点的距离,再加上该需求点到对应回归车场的距离,转入步骤3,如果没有车辆的剩余载重满足,转入step2
step2:新增一辆车,
计算车场到该需求点的距离,再加上需求点到车场的距离,车场的最优车辆只有一个,即最优
step3:比较各个车辆对应的计算距离,最优个体是距离最短对应的车辆
3个车场的最优个体做比较,距离最短是整个系统的最优个体。这种车辆选择方式采用一种贪婪的思想,寻找使路径最短的车辆,一定程度上会加快最优路径的寻找。
核心代码如下:
def fitness(self,road):
p1, p2, p3 = [[]], [[]], [[]] # 存车场到需求点的路径
dis_p =[0, 0, 0]
d1, d2, d3 = [[]], [[]], [[]]
for w in range(len(road)):
point = road[w] # 新的点
p = [p1, p2, p3]
dis = [self.dis1, self.dis2, self.dis3] # 3个车场到需求点的距离矩阵
s = []
d = [d1, d2, d3]
location = []
for i in range(3):
p_now = p[i] # 依次从第1车场到第3车场
dis_now = dis[i]
d_now = d[i]
signal = 0
if len(p_now[-1]) > 0: # 如果最后一个车的装载情况非空
dx = []
loc = []
for j in range(len(p_now)): # 依次遍历对应车场的各个车
now_load = sum(d_now[j]) # 计算每个车的装载情况
if now_load + self.demand.iloc[point-1,-1] <= self.load: # 如果可以装
dw = dis_now[p_now[j][-1], point] + dis_now[point,0] # 计算其与上一个点的距离和回到车场的距离
dx.append(dw)
loc.append(j)
if len(dx)>0:
di = min(dx) # 挑选最小距离对应的车
idx = dx.index(min(dx))
best_idx = loc[idx]
if len(dx)==0: # 如果各个车都不能装,进入下一步
signal=1
if len(p_now[-1]) == 0 or signal==1: # 如果车场最后一个车是空车,或者上一步没装成功
now_car = len(p_now) # 现在的车场对应的车辆使用情况
if now_car + 1 <= self.car_number: # 如果使用车辆加1小于等于20辆,增加新车
if i ==0 :
p1.append([])
d1.append([])
if i ==1 :
p2.append([])
d2.append([])
if i ==2:
p3.append([])
d3.append([])
else: # 如果使用车辆加1大于20辆,距离赋予一个很大的值,对应车不能装,也没有增加新车
di = 100000
best_idx = 111111
if len(p_now[-1]) == 0: # 如果有新车
di = dis_now[0, point] + dis_now[point, 0] # 计算车场到需求点及回去的距离
best_idx = -1
s.append(di)
location.append(best_idx)
idx = s.index(min(s)) # 挑选最小距离对应的车场
idd = location[idx] # 对应车场对应的车辆选择
if idx ==0 : # 车场1对应车辆添加需求点,下面相同
p1[idd].append(point)
d1[idd].append(self.demand.iloc[point-1,-1])
if idx ==1 :
p2[idd].append(point)
d2[idd].append(self.demand.iloc[point-1,-1])
if idx ==2:
p3[idd].append(point)
d3[idd].append(self.demand.iloc[point-1,-1])
dis_p[idx] += min(s) # 各个车场对应的距离累计
return p1, p2, p3, dis_p, d1, d2, d3
粒子群算法
PSO初始化为一群随机粒子(随机解)。然后通过迭代找到最优解。在每一次的迭代中,粒子通过跟踪两个“极值”(pbest,gbest)来更新自己,涉及的参数有:粒子速度v、粒子位置(解)x 、学习因子c1和c2通常为2、惯性因子w等:
粒子更新如下
v = self.w*v + self.c1*(p_best - x) + self.c2*(g_best - x) # 速度更新
x = x +v # 解更新
参考博客:http://t.csdn.cn/bmK17
1、编码转化
粒子群算法主要求解连续问题,求解离散的路径问题需要转化编码,设计一种转化编码方式如下:
初始100个-5到5的数:
x = -np.random.rand(100)*5
[-4.222859939408526, -0.7378758325871354, -2.531496290674307, -3.081918679218121, -0.9450792327097823, -4.569715731364068, -2.7400011815357543, -2.99650451798587, -0.41349659498989333, -2.2288969851054077, -0.36786270261017273, -2.4468354315767797, -3.910812142856349, -3.0637623402366927, -2.332315091979022, -2.7177602325404404, -4.541267756871481, -4.69697312686115, -3.3059924360977395, -3.4841095762643457, -2.877548844377139, -4.032728456165723, -3.4404545972375136, -2.8036370303251474, -4.951591253304101, -4.910429728258995, -0.6112195838203671, -2.6354585990909767, -2.0455702528915625, -0.17405514344649486, -1.0852636011496548, -2.9215978086986176, -0.7458159804853137, -3.382498681066565, -4.3725684498094886, -4.258597260063095, -4.984384929086641, -0.5502586595988224, -3.279919057173333, -2.0696624557821193, -4.882834003166256, -1.1836213215106584, -0.5368030576588034, -0.8599180541783724, -4.481553817796012, -3.1365242704040797, -1.101455323388079, -1.2610501451145457, -0.7999963652901365, -4.1597085217213055, -0.31215414224349314, -0.48409957805062853, -1.6225821453698264, -4.071147896543535, -0.7783607060733533, -3.4928005090201375, -2.5489785717902302, -4.202306093344149, -4.019341001826704, -0.5439515321829191, -1.9713569149098138, -3.59038876561098, -3.0362106346135787, -0.5740948051844869, -3.7563446101821634, -4.47351967056031, -4.727089965125739, -3.2135077433468706, -2.918742796048436, -1.8208682778518437, -4.4087235517825345, -1.9292196302268976, -3.2576439754676185, -3.225210925945394, -4.635485881541316, -0.6367643314114257, -1.13145866646567, -4.426101187578101, -3.4110716571536894, -0.7656637896100027, -3.6487082371118253, -0.9115766194806096, -0.6561367289343345, -2.2503671433024715, -0.42554179259105596, -2.6159236296092225, -1.7939737830191467, -0.5509601511604301, -4.173784388393538, -3.269503035137562, -1.1500767686495, -2.86664525312616, -4.351846319047835, -1.4664609078101716, -1.6230513700840459, -1.1637671361175472, -1.116896850654382, -3.8808744039222383, -3.894695781089468, -0.7318115397339747]
每一个数对应从小到大的索引
idx = np.argsort(x) # 从小到大位置索引
[36, 24, 25, 40, 66, 17, 74, 5, 16, 44, 65, 77, 70, 34, 92, 35, 0, 57, 88, 49, 53, 21, 58, 12, 98, 97, 64, 80, 61, 55, 19, 22, 78, 33, 18, 38, 89, 72, 73, 67, 45, 3, 13, 62, 7, 31, 68, 20, 91, 23, 6, 15, 27, 85, 56, 2, 11, 14, 83, 9, 39, 28, 60, 71, 69, 86, 94, 52, 93, 47, 41, 95, 90, 76, 96, 46, 30, 4, 81, 43, 48, 54, 79, 32, 1, 99, 82, 75, 26, 63, 87, 37, 59, 42, 51, 84, 8, 10, 50, 29]
初始road是1到100按顺序的排列,转化一个可行的编码:
road_init = np.arange(1,101)
road = road_init[idx] # 编码转化
[37, 25, 26, 41, 67, 18, 75, 6, 17, 45, 66, 78, 71, 35, 93, 36, 1, 58, 89, 50, 54, 22, 59, 13, 99, 98, 65, 81, 62, 56, 20, 23, 79, 34, 19, 39, 90, 73, 74, 68, 46, 4, 14, 63, 8, 32, 69, 21, 92, 24, 7, 16, 28, 86, 57, 3, 12, 15, 84, 10, 40, 29, 61, 72, 70, 87, 95, 53, 94, 48, 42, 96, 91, 77, 97, 47, 31, 5, 82, 44, 49, 55, 80, 33, 2, 100, 83, 76, 27, 64, 88, 38, 60, 43, 52, 85, 9, 11, 51, 30]
2、局部最优更新和全局最优更新
设置初始的路径为局部最优解,以后的每一次迭代:如果新的路径优于原路径,更新局部最优,多个局部最优的最优值,就是全局最优。
初始v是-2到2之间的100个数。
核心代码:
for i in range(self.popsize):
x = Xx[i]
p_best = P_best[i]
v = Vx[i]
v = self.w*v + self.c1*(p_best - x) + self.c2*(g_best - x) # 速度更新
x = x +v # 解更新
V[i] = v
X[i] = x
idx = np.argsort(x)
road = road_init[idx] # 新路径
_, _, _, dis_p,_, _, _ = self.de.fitness(road) # 路径长度计算
if sum(dis_p) < answer1[i]:
P_best[i] = x # 局部最优更新
answer1[i] = sum(dis_p)
idx = answer1.index(min(answer1))
g_best = P_best[idx] # 全局最优
模拟退火算法
模拟退火算法(SA)来源于固体退火原理,是一种基于概率的算法。将固体加温至充分高的温度,再让其徐徐冷却,加温时,固体内部粒子随温升变为无序状,内能增大,分子和原子越不稳定。而徐徐冷却时粒子渐趋有序,能量减少,原子越稳定。在冷却(降温)过程中,固体在每个温度都达到平衡态,最后在常温时达到基态,内能减为最小。
简单来说:对一组解进行更新产生新解,如果新的解劣于原解,在一定概率下选择新解,温度高时概率大,退火后概率小。
Metropolis算法下的概率计算公式如下图:
算法步骤:
(1) 初始化:初始温度T(充分大),初始解状态S(是算法迭代的起点),每个T值的迭代次数L
(2) 对k=1, …, L做第(3)至第6步:
(3) 产生新解S′
(4) 计算增量ΔT=C(S′)-C(S),其中C(S)为目标函数,C(S)相当于能量
(5) 若ΔT<0则接受S′作为新的当前解,否则以概率exp(-ΔT/T)接受S′作为新的当前解.
(6) 如果满足终止条件则输出当前解作为最优解,结束程序。
(7) T逐渐减少,且T->0,然后转第2步。
参考博客:http://t.csdn.cn/4mDLx
1、Metropolis算法
比较简单,核心代码:
def Metrospolis(self, f, f_new): #Metropolis准则
if f_new <= f:
return 1
else:
p = math.exp((f - f_new) / self.T)
if random.random() < p:
return 1
else:
return 0
2、个体更新(邻域路径寻找)
对于可行的路径编码,邻域路径的寻找:1、随机生成两个位置,逆序两个位置间的编码;2、随机生成两个位置。交换两个位置的编码。
核心代码:
def serch(self, road):
index = random.sample(range(len(road)), 2) # 随机生成两个位置
idx = list(set(index)) # 位置按先后顺序调整
d = road[idx[0]:idx[1] + 1]
d.reverse() # 逆序
road[idx[0]:idx[1] + 1] = d
index = random.sample(range(len(road)), 2) # 随机生成两个位置
road[index[0]], road[index[1]] = road[index[1]], road[index[0]]
return road
3、最优个体保留设计
因为当新解劣于原解时,仍然有一定可能接受新解,设计最优个体保留规则:最优个体对应的新解不优于原解时,不更新。
核心代码:
just =self. Metrospolis(answer1[i], sum(dis_p))
if just == 1:
if best_idx == i and sum(dis_p) > answer1[i]: #当种群的最优个体的新解劣于原解时,即使Metrospolis结束,也不更新。
continue
else:
Road1[i] = road
answer1[i] = sum(dis_p)
结果
主函数:
from decode import decoding
from pso import pso
from sa import sa
parm_car = [50,20] # 载重和车场车辆数
de = decoding(parm_car) # 解码模块
information = ['粒子群', '模拟退火']
choose = [1, 1] # 非零对应算法的运行,1,0运行粒子群,0,1运行模拟退火,1,1两个都运行
if choose[0] : # 选择粒子群
print('\n',information[0])
parm_pso = [15, 10, 2, 2, .5] # 迭代次数,种群规模,学习因子c1,c2和惯性因子w
ps = pso(parm_pso,de)
road2, result = ps.total()
de.draw(road2) # 每次迭代的路径动画展示
de.draw_change(result) # 路径总长度变化
if choose[1]:
print('\n',information[1])
parm_sa = [100,90,0.99,10] # 其实温度,终止温度,温度下降系数、每个温度迭代规模
sa = sa(parm_sa,de)
road2, result = sa.total()
de.draw(road2) # 每次迭代的路径动画展示
de.draw_change(result) # 路径总长度变化
粒子群运行结果:
粒子群路径图:
粒子群路径长度变化:
模拟退火运行结果:
模拟退火群路径图:
模拟退火路径长度变化:
视频演示:
多车场车辆路径优化丨MDCVRP的粒子群+模拟退火算法求解
完整代码
完整算法源码+数据:见下方微信公众号:关注后回复:路径优化
# 微信公众号:学长带你飞
# 主要更新方向:1、车辆路径问题求解算法
# 2、学术写作技巧
# 3、读书感悟
# @Author : Jack Hao
百度云有可能失效,若失效可加微信:diligent-1618,请备注:学校-专业-昵称。