【进阶六】Python实现SDVRPTW常见求解算法—

基于python语言，采用经典遗传算法（GA）对 带硬时间窗的需求拆分车辆路径规划问题（SDVRP） 进行求解。

往期优质资源

经过一年多的创作，目前已经成熟的代码列举如下，如有需求可私信联系，表明需要的 **问题与算法**，原创不宜，有偿获取。

VRP问题	GA	ACO	ALNS	DE	DPSO	QDPSO	TS	SA
CVRP	√	√	√	√	√	√	√	√
VRPTW	√	√	√	√	√	√	√	√
MDVRP	√	√	√	√	√	√	√	√
MDHVRP	√	√	√	√	√	√	√	√
MDHVRPTW	√	√	√	√	√	√	√	√
SDVRP	√	√	√	√	√	√	√	√
SDVRPTW	√

1. 适用场景

求解SDVRPTW
车辆类型单一
车辆容量小于部分需求节点需求
单一车辆基地
带硬时间窗

2. 代码调整

2.1 需求拆分

与SDVRP问题相比，SDVRPTW问题不仅允许客户需求大于车辆载重，而且考虑了客户节点的时间窗约束。为了使得每个客户的需求得到满足，必须派遣一辆或多辆车辆在规定时间窗内对客户进行服务。对于需求节点的拆分，这里依然采取先验拆分策略，本文采用文献[1]提出的先验分割策略，表述如下：

（1）20/10/5/1拆分规则

m₂₀ =max{ $m\in Z^+ \cup \{0\} | 0.20Qm <= D_i$ }
m₁₀ =max{ $m\in Z^+ \cup \{0\} | 0.10Qm <= D_i-0.20Qm_{20}~$ }
m₅ =max{ $m\in Z^+ \cup \{0\} | 0.05Qm <= D_i-0.20Qm_{20}-0.10Qm_{10}$ }
m₁ =max{ $m\in Z^+ \cup \{0\} | 0.01Qm <= D_i-0.20Qm_{20}-0.10Qm_{10}-0.05Qm_{5}$ }

（2）25/10/5/1拆分规则

m₂₅ =max{ $m\in Z^+ \cup \{0\} | 0.25Qm <= D_i$ }
m₁₀ =max{ $m\in Z^+ \cup \{0\} | 0.10Qm <= D_i-0.25Qm_{25}~$ }
m₅ =max{ $m\in Z^+ \cup \{0\} | 0.05Qm <= D_i-0.25Qm_{25}-0.10Qm_{10}$ }
m₁ =max{ $m\in Z^+ \cup \{0\} | 0.01Qm <= D_i-0.25Qm_{25}-0.10Qm_{10}-0.05Qm_{5}$ }

在实现过程中，对于需求超过车辆容量的客户必须进行需求拆分，而对于未超过车辆容量的客户可以拆分也可以不拆分，这里设置了参数比例进行限制。

2.2 需求拆分后的服务时长取值问题

节点的服务时长会影响车辆的行进时间，进而会影响与节点时间窗的匹配问题。一般来说，节点的服务时长与需求量成正比关系，在进行节点需求拆分后，新节点的需求量降低，其服务时长理应也降低。但从标准数据集来看，各需求节点的服务时长均采用同一数值。因此本文在代码实现过程中也采用固定值，不考虑新节点服务时长的变化。当然，如有需要，也可以设置单位货物的服务时长，根据拆分后节点的具体需求量设置相应的服务时长。

3. 求解结果

（1）收敛曲线

在这里插入图片描述

（2）车辆路径
在这里插入图片描述
（3）输出内容

在这里插入图片描述

4. 代码片段

（1）数据结构

import math
import random
import numpy as np
import copy
import xlsxwriter
import matplotlib.pyplot as plt
import csv
import time
# 数据结构：解
class Sol():def __init__(self):self.obj=None # 目标函数值self.fitness = 0self.node_no_seq=[] # 解的编码self.route_list=[] # 解的解码self.timetable_list=[] # 车辆访问各点的时间self.route_distance_list = None
# 数据结构：需求节点
class Node():def __init__(self):self.id=0 # 节点idself.x_coord=0 # 节点平面横坐标self.y_coord=0  # 节点平面纵坐标self.demand=0 # 节点需求self.start_time=0 # 节点开始服务时间self.end_time=1440 # 节点结束服务时间self.service_time=0 # 单次服务时长self.vehicle_speed = 0 # 行驶速度
# 数据结构：车场节点
class Depot():def __init__(self):self.id=0 # 节点idself.x_coord=0 # 节点平面横坐标self.y_cooord=0  # 节点平面纵坐标self.start_time=0 # 节点开始服务时间self.end_time=1440 # 节点结束服务时间self.v_speed = 0 # 行驶速度self.v_cap = 80 # 车辆容量
# 数据结构：全局参数
class Model():def __init__(self):self.best_sol=None # 全局最优解self.sol_list=[] # 解的集合self.demand_dict = {}  # 需求节点集合self.depot = None  # 车场节点集合self.demand_id_list = [] # 需求节点id集合self.distance_matrix = {}  # 距离矩阵self.time_matrix = {}  # 时间矩阵self.number_of_demands = 0 # 需求点数量self.popsize=100 # 种群规模self.pc = 0.5  # 交叉概率self.pm = 0.2  # 变异概率self.n_select = 80  # 种群选择数量self.demand_id_list_ = []  # 经先验需求分割后的节点集合self.demand_dict_ = {}  # 需求分割后的节点需求集合self.distance_matrix_ = {}  # 原始节点id间的距离矩阵self.time_matrix_ = {}  # 原始节点id间的时间矩阵self.mapping = {}  # 需求分割前后的节点对应关系self.split_rate = 0.5 # 控制需求分割的比例（需求超出车辆容量的除外）

（2）距离矩阵

# 初始化参数：计算距离矩阵时间矩阵
def calDistanceTimeMatrix(model):for i in range(len(model.demand_id_list)):from_node_id = model.demand_id_list[i]for j in range(len(model.demand_id_list)):to_node_id = model.demand_id_list[j]dist = math.sqrt((model.demand_dict[from_node_id].x_coord - model.demand_dict[to_node_id].x_coord) ** 2+ (model.demand_dict[from_node_id].y_coord - model.demand_dict[to_node_id].y_coord) ** 2)model.distance_matrix[from_node_id, to_node_id] = distmodel.time_matrix[from_node_id,to_node_id] = math.ceil(dist/model.depot.v_speed)dist = math.sqrt((model.demand_dict[from_node_id].x_coord - model.depot.x_coord) ** 2 +(model.demand_dict[from_node_id].y_coord - model.depot.y_coord) ** 2)model.distance_matrix[from_node_id, model.depot.id] = distmodel.distance_matrix[model.depot.id, from_node_id] = distmodel.time_matrix[from_node_id,model.depot.id] = math.ceil(dist/model.depot.v_speed)model.time_matrix[model.depot.id,from_node_id] = math.ceil(dist/model.depot.v_speed)

（3）邻域搜索

# 计算适应度
def calFit(model):#calculate fit value：fit=Objmax-objmax_obj=-float('inf')best_sol=Sol()#record the local best solutionbest_sol.obj=float('inf')for sol in model.sol_list:node_no_seq=sol.node_no_seqsol.route_list= splitRoutes(node_no_seq, model)sol.timetable_list,sol.obj, sol.route_distance = calTravelCost(sol.route_list, model)if sol.obj > max_obj:max_obj = sol.objif sol.obj < best_sol.obj:best_sol = copy.deepcopy(sol)#calculate fit valuefor sol in model.sol_list:sol.fitness = max_obj-sol.obj#update the global best solutionif best_sol.obj<model.best_sol.obj:model.best_sol=best_sol
# 二元锦标赛
def selectSol(model):sol_list=copy.deepcopy(model.sol_list)model.sol_list=[]for _ in range(model.n_select):f1_index=random.randint(0,len(sol_list)-1)f2_index=random.randint(0,len(sol_list)-1)f1_fit=sol_list[f1_index].fitnessf2_fit=sol_list[f2_index].fitnessif f1_fit<f2_fit:model.sol_list.append(sol_list[f2_index])else:model.sol_list.append(sol_list[f1_index])
# OX交叉
def crossSol(model):sol_list=copy.deepcopy(model.sol_list)model.sol_list=[]while True:[f1_index,f2_index] = random.sample(range(len(sol_list)),2)f1 = copy.deepcopy(sol_list[f1_index])f2 = copy.deepcopy(sol_list[f2_index])if random.random() <= model.pc:cro1_index = random.randint(0,model.number_of_demands-1)cro2_index = random.randint(cro1_index,model.number_of_demands-1)new_c1_f = []new_c1_m=f1.node_no_seq[cro1_index:cro2_index+1]new_c1_b = []new_c2_f = []new_c2_m=f2.node_no_seq[cro1_index:cro2_index+1]new_c2_b = []for index in range(model.number_of_demands):if len(new_c1_f)<cro1_index:if f2.node_no_seq[index] not in new_c1_m:new_c1_f.append(f2.node_no_seq[index])else:if f2.node_no_seq[index] not in new_c1_m:new_c1_b.append(f2.node_no_seq[index])for index in range(model.number_of_demands):if len(new_c2_f)<cro1_index:if f1.node_no_seq[index] not in new_c2_m:new_c2_f.append(f1.node_no_seq[index])else:if f1.node_no_seq[index] not in new_c2_m:new_c2_b.append(f1.node_no_seq[index])new_c1=copy.deepcopy(new_c1_f)new_c1.extend(new_c1_m)new_c1.extend(new_c1_b)f1.node_no_seq=new_c1new_c2=copy.deepcopy(new_c2_f)new_c2.extend(new_c2_m)new_c2.extend(new_c2_b)f2.node_no_seq=new_c2model.sol_list.append(copy.deepcopy(f1))model.sol_list.append(copy.deepcopy(f2))else:model.sol_list.append(copy.deepcopy(f1))model.sol_list.append(copy.deepcopy(f2))if len(model.sol_list)>model.popsize:break
# 变异
def muSol(model):sol_list=copy.deepcopy(model.sol_list)model.sol_list=[]while True:f1_index = random.randint(0, len(sol_list) - 1)f1 = copy.deepcopy(sol_list[f1_index])m1_index=random.randint(0,model.number_of_demands-1)m2_index=random.randint(0,model.number_of_demands-1)if m1_index!=m2_index:if random.random() <= model.pm:node1=f1.node_no_seq[m1_index]f1.node_no_seq[m1_index]=f1.node_no_seq[m2_index]f1.node_no_seq[m2_index]=node1model.sol_list.append(copy.deepcopy(f1))if len(model.sol_list)>model.popsize:break