700字范文 > python调用开源求解器SCIP求解带时间窗车辆路径问题（VRPTW）

python调用开源求解器SCIP求解带时间窗车辆路径问题（VRPTW）

时间：2024-04-10 02:37:42

文章目录

1. 问题定义2. 数学模型3. python调用SCIP实现代码4. 结果参考文献

1. 问题定义

带时间窗车辆路径问题（vehicle routing problem with time windows，VRPTW）是在CVRP (Capacitated vehicle routing problem)问题基础上引入了顾客需求点对于车辆到达时间的要求，即每个顾客点都有一个对应的[最早到达时间，最晚到达时间]时间窗限制。

因此在该问题的目标函数中，除了要考虑车辆的行驶成本还包括车辆提前到达客户点等待的时间和顾客需要的服务时间。对于时间窗的处理，常见的有两种方式：

硬时间窗，要求车辆必须在客户要求的服务时间窗口内到达，提前到达需等待，延误到达则拒收；软时间窗，车辆可以不在服务时间窗内到达，但无论提前还是延误都会进行一定的惩罚。

2. 数学模型

建立multi-commodity network ﬂow model，数学模型直接参考教材《Column Generation》第三章

线性化及大M的取值

Solomon算例下载地址：https://www.sintef.no/projectweb/top/vrptw/100-customers/

3. python调用SCIP实现代码

import osimport reimport mathimport pandas as pdimport numpy as npimport pyscipopt as optimport matplotlib.pyplot as pltclass VrpTW(object):def __int__(self):self.customerNum = 0 # 顾客点数量self.vehicleNum = 0 # 车辆数self.capacity = 0 # 车容量self.cor_X = [] # x轴坐标self.cor_Y = [] # y轴坐标self.demand= [] # 顾客点的需求self.readyTime = [] # 顾客点左时间窗self.dueTime= [] # 顾客点右时间窗self.serviceTime = [] # 车辆到达节点需服务的时间self.nodeNum= 0 # 顶点数量self.disMatrix = [[]] # 路径成本def read_data(self):# 通用算例下载地址：https://www.sintef.no/projectweb/top/vrptw/100-customers/with open('./solomon/solomon_25/R101.txt') as file:lines = file.readlines()# 第5行 NUMBER CAPACITYself.vehicleNum, self.capacity = list(map(int, re.split(r" +", lines[4][:-1].strip())))# 第10行到最后一行分别存储. CUST NO, XCOORD, YCOORD, DEMAND, READY TIME, DUE DATE, SERVICE TIMEcust_data = np.array([list(map(int, re.split(r" +", i[:-1].strip()))) for i in lines[9:]])self.cor_X = cust_data[:, 1]self.cor_Y = cust_data[:, 2]self.demand = cust_data[:, 3]self.readyTime = cust_data[:, 4]self.dueTime = cust_data[:, 5]self.serviceTime = cust_data[:, 6]def prepare_date(self):"""根据不同的算例需要手动配置的车辆数及顾客点"""self.vehicleNum = 8 # 手动设置车辆数self.customerNum = 25 # 手动设置取出顾客点的数量"""# 处理取出顾客点及仓库点的数据，并为了方便处理将仓库点复制插入至最后# 假设有n个顾客点 C = {1,2, ...N}，并将仓库点（0）和仓库点的复制作为第n+1一个虚拟点 (这样nodeNum=customerNum + 2)# 举例顾客点集合 C = {1,2, ...N}, 顶点集合 N = C U {0, n+1}"""row_num = self.customerNum + 1 # 设置截取数据集的row_num行，其中第1行为仓库点，所以行数=customerNum+1self.cor_X = np.append(self.cor_X[:row_num], self.cor_X[0])self.cor_Y = np.append(self.cor_Y[:row_num], self.cor_Y[0])self.demand = np.append(self.demand[:row_num], self.demand[0])self.readyTime = np.append(self.readyTime[:row_num], self.readyTime[0])self.dueTime = np.append(self.dueTime[:row_num], self.dueTime[0])self.serviceTime = np.append(self.serviceTime[:row_num], self.serviceTime[0])self.nodeNum = self.customerNum + 2# 计算欧式距离矩阵self.disMatrix = np.full((self.nodeNum, self.nodeNum), 0)for i in range(self.nodeNum):for j in range(self.nodeNum):self.disMatrix[i][j] = math.sqrt((self.cor_X[i] - self.cor_X[j]) ** 2 + (self.cor_Y[i] - self.cor_Y[j]) ** 2)# 大M的取值, 拉紧约束, 提高求解效率M = np.full((self.nodeNum, self.nodeNum), 0)for i in range(self.nodeNum):for j in range(self.nodeNum):M[i][j] = self.dueTime[i] + self.disMatrix[i, j] - self.readyTime[i]self.bigM = np.max(M)print('customerNum={}, nodeNum={}, vehicleNum={}, bigM={}'.format(self.customerNum, self.nodeNum, self.vehicleNum, self.bigM))def plt_route_pic(self, vrp_route):plt.figure(figsize=(10, 8))# 绘制坐标的散点图for i in range(0, self.nodeNum-1):if i == 0:plt.scatter(self.cor_X[i], self.cor_Y[i], marker='^', c='r')plt.annotate('depot', (self.cor_X[i], self.cor_Y[i]))else:plt.scatter(self.cor_X[i], self.cor_Y[i], marker='o', c='black')plt.annotate(str(i), (self.cor_X[i], self.cor_Y[i]))# 绘制箭头for idx in vrp_route:i,j,k = idxplt.arrow(self.cor_X[i], self.cor_Y[i], self.cor_X[j] - self.cor_X[i], self.cor_Y[j] - self.cor_Y[i],length_includes_head = True, head_width =0.5, head_length=0.8, fc='red', ec='red',linewidth=0.3)plt.xlabel('x')plt.ylabel('y')plt.title('VRPTW')plt.savefig('./vrptw.png')def solve_vrptw(self):model = opt.Model()# ==========定义变量==========x = {}s = {}# x_i_j_k: 0-1变量，表示车辆k经过弧（i，j）, i != jfor i in range(self.nodeNum):for j in range(self.nodeNum):if i != j:for k in range(self.vehicleNum):x[i, j, k] = model.addVar(vtype='B', name='x_' + str(i) + '_' + str(j) + '_' + str(k))# s_i_k: 连续变量，表示车辆k开始服务顾客点i的开始时间for i in range(self.nodeNum):for k in range(self.vehicleNum):s[i, k] = model.addVar(vtype='C', name='s_' + str(i) + '_' + str(k))# ==========定义约束==========# 约束1: 对于每个车辆: 驶出仓库for k in range(self.vehicleNum):model.addCons(opt.quicksum(x[0, j, k] for j in range(self.nodeNum) if j != 0) == 1, name='vehicle_depart_' + str(k))# 约束2: 对于每个车辆: 驶回仓库for k in range(self.vehicleNum):model.addCons(opt.quicksum(x[i, self.nodeNum-1, k] for i in range(self.nodeNum) if i != self.nodeNum-1) == 1, name='vehicle_return_' + str(k))# 约束3：对于每个车辆: 车辆的容量约束 (3.3)for k in range(self.vehicleNum):model.addCons(opt.quicksum(self.demand[i] * x[i, j, k] for i in range(1, self.nodeNum-1) for j in range(self.nodeNum) if i!=j) <= self.capacity,name='capacity_vehicle' + str(i))# 约束4：对于顾客点: 进的车数量=1，出的车的数量=1, 即保证每个客户点都被服务 for i in range(1, self.nodeNum-1):model.addCons(opt.quicksum(x[i, j, k] for j in range(self.nodeNum) for k in range(self.vehicleNum) if i!=j) == 1,name='customer_depart_' + str(i))model.addCons(opt.quicksum(x[j, i, k] for j in range(self.nodeNum) for k in range(self.vehicleNum) if i!=j) == 1,name='customer_enter_' + str(i))# 约束5: 对于顾客点: 进出是同一辆车，流平衡，服务之后需离开 for k in range(self.vehicleNum):for i in range(1, self.nodeNum-1):model.addCons(opt.quicksum(x[i, j, k] for j in range(self.nodeNum) if i!=j) ==opt.quicksum(x[j, i, k] for j in range(self.nodeNum) if i!=j), name='flow_' + str(k) + '_' + str(i))"""约束6: 对于顾客点和仓库的时间窗口约束对于顾客点的时间窗口: 只有车辆k服务了顾客点i此约束,才会有时间窗口的约束；服务时间需在窗口范围内"""ready_time, due_time = {}, {} # 临时变量for i in range(1, self.nodeNum-1):for k in range(self.vehicleNum):ready_time[i, k] = model.addVar(vtype='C', name='ready_time_' + str(i) + '_' + str(k))due_time[i, k] = model.addVar(vtype='C', name='due_time_' + str(i) + '_' + str(k))model.addCons(ready_time[i, k] == opt.quicksum(self.readyTime[i] * x[i, j, k] for j in range(self.nodeNum) if i!=j))model.addCons(due_time[i, k] == opt.quicksum(self.dueTime[i] * x[i, j, k] for j in range(self.nodeNum) if i!=j))model.addCons(ready_time[i, k] <= s[i, k], name='cons_ready_time_' + str(i) + '_' + str(k))model.addCons(due_time[i, k] >= s[i, k], name='cons_due_time_' + str(i) + '_' + str(k))# 对于仓库点的时间窗口约束: 在时间点0在仓库，最终有回来时间时间要求for k in range(self.vehicleNum):model.addCons(s[0, k] == 0) # 所有车出发时间为0model.addCons(s[self.nodeNum-1, k] >= self.readyTime[self.nodeNum-1])model.addCons(s[self.nodeNum-1, k] <= self.dueTime[self.nodeNum-1])"""约束7: 保证被服务的相邻节点间访问时间顺序(消除除子回路)x_i_j_k * (s_i_k + t_i_j - s_j_k) <= 0 线性化为 (s_i_k + t_i_j - s_j_k) - （1 - x_i_j_k）* M <= 0"""for k in range(self.vehicleNum):for i in range(self.nodeNum):for j in range(self.nodeNum):if(i != j):model.addCons(s[i,k] + self.disMatrix[i][j] + self.serviceTime[i] - s[j,k] - self.bigM + self.bigM * x[i,j,k] <= 0,name= 'time_windows_' + str(i) + '_' + str(j) + '_' + str(k))# ==========定义目标==========model.setObjective(opt.quicksum(x[i, j, k] * self.disMatrix[i, j] for (i, j, k) in x))model.setMinimize()model.optimize()# ==========输出结果==========print('model_status = ', model.getStatus())print('model_gap =', model.getGap())print('model_obj =', model.getObjVal())# model.writeProblem('vrpTW.lp')vrptw_route = []for k in range(self.vehicleNum):for i in range(self.nodeNum):for j in range(self.nodeNum):if(i != j and model.getVal(x[i,j,k]) > 0):vrptw_route.append((i,j,k))print('vrptw_route:\n', vrptw_route)# 绘制结果路径self.plt_route_pic(vrptw_route)if __name__ == '__main__':vrp = VrpTW()# 读取数据vrp.read_data()# 预处理数据vrp.prepare_date()# 求解模型并绘制结果vrp.solve_vrptw()