引言

本系列文章是路径优化问题学习过程中一个完整的学习路线。问题从简单的单车场容量约束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)                     # 路径总长度变化

粒子群运行结果：

粒子群路径图：

粒子群路径长度变化：

模拟退火运行结果：

模拟退火群路径图：

模拟退火路径长度变化：

视频演示：

完整代码

完整算法源码+数据：见下方微信公众号：关注后回复：路径优化

# 微信公众号：学长带你飞
# 主要更新方向：1、车辆路径问题求解算法
#              2、学术写作技巧
#              3、读书感悟
# @Author  : Jack Hao

百度云有可能失效，若失效可加微信：diligent-1618，请备注：学校-专业-昵称。