Python | 聚焦MDVRPTW系列求解算法

摘要

     七大优化算法源码投递【蚁群算法、自适应大邻域算法、差分进化算法、粒子群算法、遗传算法、模拟退火算法、禁忌搜索算法】

致力于打造交通领域优质代码分享平台

祝大家五一快乐，外出遵守防疫政策，做好个人防护。

今天为大家分享MDVRPTW系列求解算法源码。欢迎：点赞+收藏+关注。后续将推送更多优质内容，为避免错过一个亿，建议星标。

前两期应用元启发式算法求解CVRP和MDCVRP时采用了简单的路径分割规则，其优点是简单易实现，但求解质量一般。若在此基础上进一步考虑时间窗约束，将进一步降低寻优效果。为了能够在时间窗约束下寻求更好的求解效果，今天在CVRP和MDCVRP代码基础上复现文献中的算法对Split函数进行改进，使其能够处理多（单）车场条件下带有时间窗约束的车辆路径规划问题。

01

注意事项

      · 求解MDVRPTW或VRPTW

· 车辆类型单一

· 车辆容量不小于需求节点最大需求

· （多）车辆基地

· 各车场车辆总数满足实际需求

· 节点时间窗至少满足车辆路径只包含一个需求节点的情况

· 车场车辆总数满足实际需求

· 服务必须在要求时间窗内开始和结束

02

实现思路

Christian Prins[2]给出的CVRP的Split方法过程如下：

在此基础上Lin Chengs[3]给出了Split在VRPTW问题上的改进：

我们将在以上算法流程上略加改进，使其具备处理MDVRPTW和VRPTW的能力。

03

数据结构

继续以类的形式对算法参数、编码、解码、目标函数等进行封装。相比MDCVRP系列内容，有以下调整：

 · Node数据结构类中增加‘start_time’属性表示开始服务时间，‘end_time’属性表示结束服务时间，‘service_time’表示服务时长；

· Sol数据结构中增加‘timetable_list’属性记录车辆到达每个节点的时刻，‘cost_of_distance’属性表示路径的距离成本，‘cost_of_time’表示路径的时间成本

· Model数据结构类中增加‘time_matrix’属性记录网络弧旅行时间

· splitRoutes函数参考经典方法实现

       基本数据结构定义如下：

      · Sol()类，表示一个可行解

属性     描述

obj 优化目标值

node_id_list       需求节点id有序排列集合

cost_of_distance       距离成本

cost_of_time      时间成本

route_list     车辆路径集合，对应MDVRPTW的解

timetable_list      车辆节点访问时间集合，对应MDVRPTW的解

       · Node()类，表示一个网络节点

属性     描述

id   物理节点id，需唯一

x_coord       物理节点x坐标

y_coord       物理节点y坐标

demand       物理节点需求

depot_capacity

车辆基地车队规模

start_time    最早开始服务（被服务）时间

end_time     最晚结束服务（被服务）时间

service_time       需求节点服务时间

        蚁群算法（ACO）存储算法参数部分的数据结构定义如下：

属性     描述

best_sol       全局最优解，值类型为Sol()

sol_list  蚁群集合，值类型为Sol()

demand_dict      需求节点集合（字典），值类型为Node()

depot_dict   车场节点集合（字典），值类型为Node()

depot_id_list      车场节点id集合

demand_id_list  需求节点id集合

opt_type      优化目标类型，0：最小旅行距离，1：最小时间成本

vehicle_cap 车辆容量

vehicle_speed     车辆行驶速度，用于计算旅行时间

distance_matrix  网络弧距离

time_matrix 节点旅行时间矩阵

popsize 蚁群规模

alpha     信息启发式因子

beta      期望启发式因子

Q   信息素总量

rho 信息素挥发系数

tau 网络弧信息素

tau0      路径初始信息素

       自适应大邻域搜索算法(ALNS)存储算法参数部分的数据结构定义如下（仅显示特有的数据结构定义部分）：

属性     描述

rand_d_max

随机破坏程度上限

rand_d_min

随机破坏程度下限

worst_d_max

最坏破坏程度上限

worst_d_min

最坏破坏程度下限

regret_n

次优位置个数

r1

算子奖励1

r2

算子奖励2

r3

算子奖励3

rho

算子权重衰减系数

d_weight

破坏算子权重

d_select

破坏算子被选中次数/每轮

d_score

破坏算子被奖励得分/每轮

d_history_select

破坏算子历史共计被选中次数

d_history_score

破坏算子历史共计被奖励得分

r_weight

修复算子权重

r_select

修复算子被选中次数/每轮

r_score

修复算子被奖励得分/每轮

r_history_select

修复算子历史共计被选中次数

r_history_score

修复算子历史共计被奖励得分

       差分进化算法(DE)存储算法参数部分的数据结构定义如下（仅显示特有的数据结构定义部分）：

属性     描述

Cr

交叉概率

F

缩放因子

popsize

种群规模

       离散粒子群算法(DPSO)存储算法参数部分的数据结构定义如下（仅显示特有的数据结构定义部分）：

属性     描述

pl

可行解历史最优位置

pg

全局最优解历史最优位置

v

可行解更新速度

Vmax

最大速度

w

惯性权重

c1

学习因子

c2

学习因子

       遗传算法(GA)存储算法参数部分的数据结构定义如下（仅显示特有的数据结构定义部分）：

属性     描述

pc

交叉概率

pm

突变概率

n_select

优良个体选择数量

popsize

种群规模

       模拟退火算法(SA)存储算法参数部分的数据结构定义均在前文出现过，不再赘述。

       禁忌搜索算法(TS)存储算法参数部分的数据结构定义如下（仅显示特有的数据结构定义部分）：

属性     描述

tabu_list

禁忌表

TL

算子禁忌长度

04

主程序

蚁群算法（ACO）主程序如下：

def run(demand_file,depot_file,Q,tau0,alpha,beta,rho,epochs,v_cap,opt_type,popsize):

    """

    :param demand_file: demand file path

    :param depot_file: depot file path

    :param Q:信息素总量

    :param tau0: 路径信息素初始值

    :param alpha:信息启发式因子

    :param beta:期望启发式因子

    :param rho:信息挥发因子

    :param epochs:迭代次数

    :param v_cap:车辆容量

    :param opt_type:优化类型:0:最小化车辆数,1:最小化行驶距离

    :param popsize:蚁群规模

    :return:

    """

    model=Model()

    model.vehicle_cap=v_cap

    model.opt_type=opt_type

    model.alpha=alpha

    model.beta=beta

    model.Q=Q

    model.tau0=tau0

    model.rho=rho

    model.popsize=popsize

    sol=Sol()

    sol.obj=float('inf')

    model.best_sol=sol

    history_best_obj = []

    readCSVFile(demand_file,depot_file,model)

    calDistanceTimeMatrix(model)

    for ep in range(epochs):

        movePosition(model)

        upateTau(model)

        history_best_obj.append(model.best_sol.obj)

        print("%s/%s， best obj: %s" % (ep,epochs, model.best_sol.obj))

    plotObj(history_best_obj)

    plotRoutes(model)

    outPut(model)

if __name__=='__main__':

    demand_file = r'D:\code\py_projects\2021-05-24-algorithms-demos\datasets\MDVRPTW\demand.csv'

    depot_file = r'D:\code\py_projects\2021-05-24-algorithms-demos\datasets\MDVRPTW\depot.csv'

    run(demand_file,depot_file,Q=10,tau0=10,alpha=1,beta=5,rho=0.1,epochs=100,v_cap=80,opt_type=0,popsize=100)

自适应大邻域算法（ALNS）主程序如下：

def run(demand_file,depot_file,rand_d_max,rand_d_min,worst_d_min,worst_d_max,regret_n,r1,r2,r3,rho,phi,epochs,pu,v_cap,

        v_speed,opt_type):

    """

    :param demand_file: demand file path

    :param depot_file: depot 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 v_speed Vehicle free speed

    :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

    model.vehicle_speed=v_speed

    readCSVFile(demand_file,depot_file, model)

    calDistanceTimeMatrix(model)

    history_best_obj = []

    sol = Sol()

    sol.node_id_list = genInitialSol(model.demand_id_list)

    calObj(sol, 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=copy.deepcopy(new_sol)

                if new_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

                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)

    plotRoutes(model)

    outPut(model)

    print("random destory weight is {:.3f}\tselect is {}\tscore is {:.3f}".format(model.d_weight[0],

                                                                        model.d_history_select[0],

                                                                        model.d_history_score[0]))

    print("worse destory weight is {:.3f}\tselect is {}\tscore is {:.3f} ".format(model.d_weight[1],

                                                                        model.d_history_select[1],

                                                                        model.d_history_score[1]))

    print("random repair weight is {:.3f}\tselect is {}\tscore is {:.3f}".format(model.r_weight[0],

                                                                       model.r_history_select[0],

                                                                       model.r_history_score[0]))

    print("greedy repair weight is {:.3f}\tselect is {}\tscore is {:.3f}".format(model.r_weight[1],

                                                                       model.r_history_select[1],

                                                                       model.r_history_score[1]))

    print("regret repair weight is {:.3f}\tselect is {}\tscore is {:.3f}".format(model.r_weight[2],

                                                                       model.r_history_select[2],

                                                                   model.r_history_score[2]))

if __name__=='__main__':

    demand_file = r'D:\code\py_projects\2021-05-24-algorithms-demos\datasets\MDVRPTW\demand.csv'

    depot_file = r'D:\code\py_projects\2021-05-24-algorithms-demos\datasets\MDVRPTW\depot.csv'

    run(demand_file=demand_file,depot_file=depot_file,rand_d_max=0.4,rand_d_min=0.1,

        worst_d_min=5,worst_d_max=20,regret_n=5,r1=30,r2=20,r3=10,rho=0.4,

        phi=0.9,epochs=5,pu=5,v_cap=80,v_speed=1,opt_type=0)

差分进化算法(DE)主程序如下：

def run(demand_file,depot_file,epochs,Cr,F,popsize,v_cap,opt_type):

    """

    :param demand_file: 需求节点文件路径

    :param depot_file: 车场节点文件路径

    :param epochs:迭代次数

    :param Cr:差分交叉概率

    :param F:缩放因子

    :param popsize:种群规模

    :param v_cap:车辆容量

    :param opt_type:优化类型:

    :return:

    """

    model=Model()

    model.vehicle_cap=v_cap

    model.Cr=Cr

    model.F=F

    model.popsize=popsize

    model.opt_type=opt_type

    readCSVFile(demand_file,depot_file,model)

    calDistanceTimeMatrix(model)

    best_sol = Sol()

    best_sol.obj = float('inf')

    model.best_sol = best_sol

    genInitialSol(model)

    history_best_obj = []

    for ep in range(epochs):

        for i in range(popsize):

            v1=random.randint(0,len(model.demand_id_list)-1)

            sol=model.sol_list[v1]

            mu_x=muSol(model,v1)

            u=crossSol(model,sol.node_id_list,mu_x)

            new_sol=Sol()

            new_sol.node_id_list=u

            calObj(new_sol,model)

            if new_sol.obj<=sol.obj:

                sol=copy.deepcopy(new_sol)

                if sol.obj

                    model.best_sol=copy.deepcopy(sol)

            history_best_obj.append(model.best_sol.obj)

        print("%s/%s， best obj: %s" % (ep, epochs, model.best_sol.obj))

    plotObj(history_best_obj)

    plotRoutes(model)

    outPut(model)

if __name__ == '__main__':

    demand_file=r'D:\code\py_projects\2021-05-24-algorithms-demos\datasets\MDVRPTW\demand.csv'

    depot_file=r'D:\code\py_projects\2021-05-24-algorithms-demos\datasets\MDVRPTW\depot.csv'

    run(demand_file,depot_file, epochs=100, Cr=0.6,F=0.5, popsize=100, v_cap=80, opt_type=0)

粒子群算法(DPSO)主程序如下：

def run(demand_file,depot_file,epochs,popsize,Vmax,v_cap,v_speed,opt_type,w,c1,c2):

    """

    :param demand_file: demand file path

    :param depot_file: depot file path

    :param epochs: 迭代次数

    :param popsize: 邻域规模

    :param v_cap: 车辆容量

    :param Vmax :最大速度

    :param opt_type: 优化类型:0:最小化车辆数,1:最小化行驶距离

    :param w: 惯性权重

    :param c1:学习因子

    :param c2:学习因子

    :return:

    """

    model=Model()

    model.vehicle_cap=v_cap

    model.Vmax=Vmax

    model.vehicle_speed=v_speed

    model.opt_type=opt_type

    model.w=w

    model.c1=c1

    model.c2=c2

    readCSVFile(demand_file,depot_file,model)

    calDistanceTimeMatrix(model)

    history_best_obj=[]

    genInitialSol(model,popsize)

    history_best_obj.append(model.best_sol.obj)

    for ep in range(epochs):

        updatePosition(model)

        history_best_obj.append(model.best_sol.obj)

        print("%s/%s: best obj: %s"%(ep,epochs,model.best_sol.obj))

    plotObj(history_best_obj)

    plotRoutes(model)

    outPut(model)

if __name__ == '__main__':

    demand_file=r'D:\code\py_projects\2021-05-24-algorithms-demos\datasets\MDVRPTW\demand.csv'

    depot_file=r'D:\code\py_projects\2021-05-24-algorithms-demos\datasets\MDVRPTW\depot.csv'

    run(demand_file=demand_file,depot_file=depot_file,epochs=100,popsize=100,Vmax=2,

        v_cap=80,v_speed=1,opt_type=0,w=0.9,c1=1,c2=5)

遗传算法主程序(GA)主程序如下：

def run(demand_file,depot_file,epochs,pc,pm,popsize,n_select,v_cap,v_speed,opt_type):

    """

    :param demand_file: demand file path

    :param depot_file: depot file path

    :param epochs: Iterations

    :param pc: Crossover probability

    :param pm: Mutation probability

    :param popsize: Population size

    :param n_select: Number of excellent individuals selected

    :param v_cap: Vehicle capacity

    :param v_speed: Vehicle free speed

    :param opt_type: Optimization type:0:Minimize the cost of travel distance;1:Minimize the cost of travel time

    :return:

    """

    model=Model()

    model.vehicle_cap=v_cap

    model.pc=pc

    model.pm=pm

    model.popsize=popsize

    model.n_select=n_select

    model.opt_type=opt_type

    readCSVFile(demand_file,depot_file,model)

    calDistanceTimeMatrix(model)

    generateInitialSol(model)

    history_best_obj = []

    best_sol=Sol()

    best_sol.obj=float('inf')

    model.best_sol=best_sol

    for ep in range(epochs):

        calFitness(model)

        selectSol(model)

        crossSol(model)

        muSol(model)

        history_best_obj.append(model.best_sol.obj)

        print("%s/%s， best obj: %s" % (ep,epochs,model.best_sol.obj))

    plotObj(history_best_obj)

    plotRoutes(model)

    outPut(model)

if __name__=='__main__':

    demand_file = r'D:\code\py_projects\2021-05-24-algorithms-demos\datasets\MDVRPTW\demand.csv'

    depot_file = r'D:\code\py_projects\2021-05-24-algorithms-demos\datasets\MDVRPTW\depot.csv'

    run(demand_file=demand_file,depot_file=depot_file,epochs=100,pc=0.8,pm=0.1,popsize=100,

        n_select=80,v_cap=80,v_speed=1,opt_type=0)

模拟退火算法(SA)主程序如下：

def run(demand_file,depot_file,T0,Tf,deltaT,v_cap,opt_type):

    """

    :param demand_file: 需求节点数据文件

    :param depot_file: 车场节点数据文件

    :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

    readCSVFile(demand_file,depot_file,model)

    calDistanceTimeMatrix(model)

    action_list=createActions(len(model.demand_id_list))

    history_best_obj=[]

    sol=Sol()

    sol.node_id_list=genInitialSol(model.demand_id_list)

    calObj(sol,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.node_id_list = doAction(sol.node_id_list, action_list[i])

            calObj(new_sol, 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=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__':

    demand_file = r'D:\code\py_projects\2021-05-24-algorithms-demos\datasets\MDVRPTW\demand.csv'

    depot_file = r'D:\code\py_projects\2021-05-24-algorithms-demos\datasets\MDVRPTW\depot.csv'

    run(demand_file=demand_file,depot_file=depot_file,T0=6000,Tf=0.001,deltaT=0.9,v_cap=80,opt_type=0)

禁忌搜索算法(TS)主程序如下：

def run(demand_file,depot_file,epochs,v_cap,opt_type):

    """

    :param demand_file: demand file path

    :param depot_file: depot file path

    :param epochs: Iterations

    :param v_cap: Vehicle capacity

    :param opt_type: Optimization type:0:Minimize the number of vehicles,1:Minimize travel distance

    :return: 无

    """

    model=Model()

    model.vehicle_cap=v_cap

    model.opt_type=opt_type

    readCSVFile(demand_file,depot_file,model)

    calDistanceTimeMatrix(model)

    action_list=createActions(len(model.demand_id_list))

    model.tabu_list=np.zeros(len(action_list))

    history_best_obj=[]

    sol=Sol()

    sol.node_id_list=generateInitialSol(model.demand_id_list)

    calObj(sol,model)

    model.best_sol=copy.deepcopy(sol)

    history_best_obj.append(sol.obj)

    for ep in range(epochs):

        local_new_sol=Sol()

        local_new_sol.obj=float('inf')

        for i in range(len(action_list)):

            if model.tabu_list[i]==0:

                new_sol=Sol()

                new_sol.node_id_list=doAction(sol.node_id_list,action_list[i])

                calObj(new_sol,model)

                new_sol.action_id=i

                if new_sol.obj

                    local_new_sol=copy.deepcopy(new_sol)

        sol=local_new_sol

        for i in range(len(action_list)):

            if i==sol.action_id:

                model.tabu_list[sol.action_id]=model.TL

            else:

                model.tabu_list[i]=max(model.tabu_list[i]-1,0)

        if sol.obj

            model.best_sol=copy.deepcopy(sol)

        history_best_obj.append(model.best_sol.obj)

        print("%s/%s: best obj: %s"%(ep,epochs,model.best_sol.obj))

    plotObj(history_best_obj)

    plotRoutes(model)

    outPut(model)

if __name__=='__main__':

    demand_file = r'D:\code\py_projects\2021-05-24-algorithms-demos\datasets\MDVRPTW\demand.csv'

    depot_file = r'D:\code\py_projects\2021-05-24-algorithms-demos\datasets\MDVRPTW\depot.csv'

    run(demand_file=demand_file,depot_file=depot_file,epochs=100,opt_type=0,v_cap=80)

05

结果对比

这里采用solomon开源数据对算法进行测试，数据文件结构如下图所示，共100个需求点，3个车场，以.csv文件储存网络数据。

demand.csv文件内容如下图，其中第一行为标题栏，节点id从0开始递增编号。

depot.csv文件内容如下图，其中第一行为标题栏，节点id建议采用图中命名方式，需要注意的是’ capacity’字段指的是车场配备的车辆数量。

七大算法算法收敛图与优化路线图分别如下：

06

源码获取

   本文的算法程序和数据集可通过以下百度网盘链接获取：

   链接：https://pan.baidu.com/s/1DyK77M-qeeDIZn3fUvD1Fw（或点击【阅读原文】）

   提取码：xafr

08

参考

[1].汪定伟. 智能优化方法[M]. 高等教育出版社, 2007.

[2]Prins C. A simple and effective evolutionary algorithm for the vehicle routing problem[J]. Computers & Operations Research, 2004, 31(12): 1985-2002.

[3]. Cheng L . A genetic algorithmfor the vehicle routing problem with time windows[J]. university of northcarolina wilmington, 2009.

转自：交通科研Lab 2022-05-12 18:00

文章来源于Python助力交通 ，作者Python助力交通