建模实战|第二期：快速掌握ortools求解VRPTW数学模型（附python代码)
经过上一篇无形忍者-禁忌搜索算法求解带时间窗的车辆路径规划问题详解（附python代码）的介绍，我们对带时间窗的车辆路径问题（VRPTW问题）已经有了一定的了解。另外经过这一篇无形忍者-OR-Tools-VRPTW带时间窗的车辆路径问题的介绍我们对ortools求解器也有了一定的了解，这篇文章是使用ortools的自带库来解决VRPTW问题的。但是实际工作中遇到的VRPTW问题更复杂，比如约束更
·
经过上一篇无形忍者-禁忌搜索算法求解带时间窗的车辆路径规划问题详解（附python代码）的介绍，我们对带时间窗的车辆路径问题（VRPTW问题）已经有了一定的了解。另外经过这一篇无形忍者-OR-Tools-VRPTW带时间窗的车辆路径问题 的介绍我们对ortools求解器也有了一定的了解，这篇文章是使用ortools的自带库来解决VRPTW问题的。
但是实际工作中遇到的VRPTW问题更复杂，比如约束更多，或者是多目标的。这时候使用ortools中解决车辆路线的自带库会有一定的局限性。所以，本次我们尝试用ortools中的CP-SAT求解器来自行建模求解。接下来，我们先复习一下之前的介绍。
1.什么是VRPTW？
VRPTW（Vehicle Routing Problem with Time Windows，带时间窗的车辆路径问题）是一种经典的组合优化问题，是基于VRP（Vehicle Routing Problem，车辆路径问题）的扩展。在VRPTW中，除了考虑车辆如何有效地访问一组客户点以最小化总行驶成本外，还必须满足每个客户点的时间窗约束。
在VRPTW中，问题的目标是找到一组车辆的最佳路径，以便在满足以下条件的情况下最小化总成本：
①每辆车必须从仓库出发，最终返回仓库；
②每个客户只被1辆车服务；
③每辆车的载重不超过总容量Q；
④每辆车的总时间不超过Dk；
⑤车辆需要在时间窗 内服务客户i,并且需要在
时间窗内返回仓库点
因此，VRPTW需要考虑更多约束条件和限制，使得问题更具挑战性。
2.VRPTW问题建模
2.1问题背景
VRPTW 问题是不仅考虑CVRP的所有约束，还需要考虑时间窗约束，也就是每个顾客对应一个时间窗 ,其中 和
和
分别代表该点的最早到达时间和最晚到达时间。通过下面介绍的参数说明和决策变量，结合参考文献（Desaulniers et al.，2006），我们将给出VRPTW 问题的标准模型
2.2参数说明
2.3决策变量
2.4目标函数
2.5约束条件
-
约束1~4分别约束如下：每个客户只被1辆车服务；每辆车必须从仓库出发；经过一个点就离开那个点；车辆最终返回仓库，其中点n+1是仓库点的一个备份，是为了方便实现。
-
约束5为车辆的容量约束
-
约束6~7是时间窗约束，分别约束如下：车辆到达每个顾客点的时间均在时间窗内；
3.部分python代码
# -*- coding: utf-8 -*-#
# 作者: 无形忍者
# 时间: 2024/3/18 13:14
# 描述: Python调用ortools求解VRPTW数学模型
import math
import re
import time
import matplotlib.pyplot as plt
from ortools.sat.python import cp_model
class InputData:
customer_num = 0
node_num = 0
vehicle_num = 0
capacity = 0
cor_x = []
cor_y = []
demand = []
service_time = []
ready_time = []
due_time = []
dis_matrix = [[]]
def read_data(path, customer_num):
input_data = InputData()
input_data.customer_num = customer_num
if customer_num is not None:
input_data.node_num = customer_num + 2
with open(path, 'r') as f:
lines = f.readlines()
count = 0
for line in lines:
count += 1
if count == 5:
line = line[:-1]
s = re.split(r" +", line)
input_data.vehicle_num = int(s[1])
input_data.capacity = int(s[2])
elif count >= 10 and (customer_num is None or count <= 10 + customer_num):
line = line[:-1]
s = re.split(r" +", line)
input_data.cor_x.append(int(s[2]))
input_data.cor_y.append(int(s[3]))
input_data.demand.append(int(s[4]))
input_data.ready_time.append(int(s[5]))
input_data.due_time.append(int(s[6]))
input_data.service_time.append(int(s[7]))
input_data.node_num = len(input_data.cor_x) + 1
input_data.customer_num = input_data.node_num - 2
# 回路
input_data.cor_x.append(input_data.cor_x[0])
input_data.cor_y.append(input_data.cor_y[0])
input_data.demand.append(input_data.demand[0])
input_data.ready_time.append(input_data.ready_time[0])
input_data.due_time.append(input_data.due_time[0])
input_data.service_time.append(input_data.service_time[0])
# 计算距离矩阵
input_data.dis_matrix = [([0] * input_data.node_num) for _ in range(input_data.node_num)]
for i in range(input_data.node_num):
for j in range(i + 1, input_data.node_num):
distance = int(math.sqrt(
(input_data.cor_x[i] - input_data.cor_x[j]) ** 2 + (input_data.cor_y[i] - input_data.cor_y[j]) ** 2))
input_data.dis_matrix[i][j] = input_data.dis_matrix[j][i] = distance
return input_data
class Answer:
obj_val = 0
X = [[]]
S = [[]]
routes = [[]]
route_num = 0
def __init__(self, input_data, solver):
self.obj_val = solver.ObjectiveValue()
# X_ijk
self.X = [[([0] * input_data.vehicle_num) for _ in range(input_data.node_num)] for _ in
range(input_data.node_num)]
# S_ik
self.S = [([0] * input_data.vehicle_num) for _ in range(input_data.node_num)]
# routes
self.routes = []
def get_answer(input_data, solver, X, S):
answer = Answer(input_data, solver) # 假设 Answer 类需要传入 input_data
for i in range(input_data.node_num):
for k in range(input_data.vehicle_num):
answer.S[i][k] = solver.Value(S[i][k])
for j in range(input_data.node_num):
answer.X[i][j][k] = solver.Value(X[i][j][k])
for k in range(input_data.vehicle_num):
i = 0
subRoute = []
subRoute.append(i)
finish = False
while not finish:
for j in range(1, input_data.node_num):
if answer.X[i][j][k] > 0.5:
subRoute.append(j)
i = j
if j == input_data.node_num - 1:
finish = True
if len(subRoute) >= 3:
subRoute[-1] = 0
answer.routes.append(subRoute)
answer.route_num += 1
return answer
def plot_answer(answer, customer_num):
plt.xlabel("x")
plt.ylabel("y")
plt.title(f"{data_type} : {customer_num} Customers")
plt.scatter(input_data.cor_x[0], input_data.cor_y[0], c='blue', alpha=1, marker=',', linewidths=3,
label='depot') # 起点
plt.scatter(input_data.cor_x[1:-1], input_data.cor_y[1:-1], c='black', alpha=1, marker='o', linewidths=3,
label='customer') # 普通站点
for k in range(answer.route_num):
for i in range(len(answer.routes[k]) - 1):
a = answer.routes[k][i]
b = answer.routes[k][i + 1]
x = [input_data.cor_x[a], input_data.cor_x[b]]
y = [input_data.cor_y[a], input_data.cor_y[b]]
plt.plot(x, y, 'k', linewidth=1)
plt.grid(False)
plt.legend(loc='best')
plt.show()
def print_answer(answer, input_data):
for index, subRoute in enumerate(answer.routes):
distance = 0
load = 0
for i in range(len(subRoute) - 1):
distance += input_data.dis_matrix[subRoute[i]][subRoute[i + 1]]
load += input_data.demand[subRoute[i]]
print(f"Route-{index + 1} : {subRoute} , distance: {distance} , load: {load}")
def solve(input_data):
# 声明模型
model = cp_model.CpModel()
# 定义变量
X = [[[[] for _ in range(input_data.vehicle_num)] for _ in range(input_data.node_num)] for _ in
range(input_data.node_num)]
S = [[[] for _ in range(input_data.vehicle_num)] for _ in range(input_data.node_num)]
for i in range(input_data.node_num):
for k in range(input_data.vehicle_num):
S[i][k] = model.NewIntVar(input_data.ready_time[i], input_data.due_time[i], name=f'S_{i}_{k}')
for j in range(input_data.node_num): # 创建连续变量
X[i][j][k] = model.NewBoolVar(name=f"X_{i}_{j}_{k}") # 创建0-1变量
# 目标函数
objective_terms = []
for i in range(input_data.node_num):
for j in range(input_data.node_num):
if i != j:
for k in range(input_data.vehicle_num):
objective_terms.append(X[i][j][k] * input_data.dis_matrix[i][j])
model.Minimize(sum(objective_terms))
# 约束1：每个客户只被1辆车服务
for i in range(1, input_data.node_num - 1): # 客户编号1 - 50
cus_expr = []
for j in range(input_data.node_num): # 客户和仓库编号0 - 51
if i != j:
for k in range(input_data.vehicle_num):
if i != 0 and i != input_data.node_num - 1:
cus_expr.append(X[i][j][k])
model.Add(sum(cus_expr) == 1)
# 约束2：车辆必须从仓库出发
for k in range(input_data.vehicle_num):
veh_expr1 = []
for j in range(1, input_data.node_num):
veh_expr1.append(X[0][j][k])
model.Add(sum(veh_expr1) == 1)
# 记录求解开始时间
start_time = time.time()
# 求解
solver = cp_model.CpSolver()
status = solver.Solve(model)
# 输出目标函数值
if status == cp_model.OPTIMAL:
# 输出求解总用时
print(f"求解总用时: {time.time() - start_time} s")
print(f"路径总距离: {solver.ObjectiveValue()}")
answer = get_answer(input_data, solver, X, S)
plot_answer(answer, input_data.customer_num)
print_answer(answer, input_data)
else:
print('求解未成功')
if __name__ == '__main__':
# 哪个数据集
data_type = "c101"
# 数据集路径
data_path = f'./data/{data_type}.txt'
# 顾客个数设置（从上往下读取完 customer_num 个顾客为止，例如c101文件中有100个顾客点，
# 但是跑100个顾客点太耗时了，设置这个数是为了只选取一部分顾客点进行计算，用来快速测试算法）
# 如果想用完整的顾客点进行计算，设置为None即可
customer_num = 50
# 一个很大的正数
M = 10000000
# 读取数据
input_data = read_data(data_path, customer_num)
# 输出相关数据
print("-" * 20, "Problem Information", '-' * 20)
print(f'Data Type: {data_type}')
print(f'Node Num: {input_data.node_num}')
print(f'Customer Num: {input_data.customer_num}')
print(f'Vehicle Num: {input_data.vehicle_num}')
print(f'Vehicle Capacity: {input_data.capacity}')
# 建模求解
solve(input_data)4.算例文件
上面代码中使用到的算例文件如下：
C101
VEHICLE
NUMBER CAPACITY
25 200
CUSTOMER
CUST NO. XCOORD. YCOORD. DEMAND READY TIME DUE DATE SERVICE TIME
0 40 50 0 0 1236 0
1 45 68 10 912 967 90
2 45 70 30 825 870 90
3 42 66 10 65 146 90
4 42 68 10 727 782 90
5 42 65 10 15 67 90
6 40 69 20 621 702 90
7 40 66 20 170 225 90
8 38 68 20 255 324 90
9 38 70 10 534 605 90
10 35 66 10 357 410 90
11 35 69 10 448 505 90
12 25 85 20 652 721 90
13 22 75 30 30 92 90
14 22 85 10 567 620 90
15 20 80 40 384 429 90
16 20 85 40 475 528 90
17 18 75 20 99 148 90
18 15 75 20 179 254 90
19 15 80 10 278 345 90
20 30 50 10 10 73 90
21 30 52 20 914 965 90
22 28 52 20 812 883 90
23 28 55 10 732 777 90
24 25 50 10 65 144 90
25 25 52 40 169 224 90
26 25 55 10 622 701 90
27 23 52 10 261 316 90
28 23 55 20 546 593 90
29 20 50 10 358 405 90
30 20 55 10 449 504 90
31 10 35 20 200 237 90
32 10 40 30 31 100 90
33 8 40 40 87 158 90
34 8 45 20 751 816 90
35 5 35 10 283 344 90
36 5 45 10 665 716 90
37 2 40 20 383 434 90
38 0 40 30 479 522 90
39 0 45 20 567 624 90
40 35 30 10 264 321 90
41 35 32 10 166 235 90
42 33 32 20 68 149 90
43 33 35 10 16 80 90
44 32 30 10 359 412 90
45 30 30 10 541 600 90
46 30 32 30 448 509 90
47 30 35 10 1054 1127 90
48 28 30 10 632 693 90
49 28 35 10 1001 1066 90
50 26 32 10 815 880 905.运行结果
-------------------- Problem Information --------------------
Data Type: c101
Node Num: 52
Customer Num: 50
Vehicle Num: 25
Vehicle Capacity: 200
求解总用时: 121.85331320762634 s
路径总距离: 353.0
Route-1 : [0, 5, 3, 7, 8, 10, 11, 9, 6, 4, 2, 1, 0] , distance: 56 , load: 160
Route-2 : [0, 13, 17, 18, 19, 15, 16, 14, 12, 0] , distance: 95 , load: 190
Route-3 : [0, 32, 33, 31, 35, 37, 38, 39, 36, 34, 0] , distance: 95 , load: 200
Route-4 : [0, 43, 42, 41, 40, 44, 46, 45, 48, 50, 49, 47, 0] , distance: 57 , load: 140
Route-5 : [0, 20, 24, 25, 27, 29, 30, 28, 26, 23, 22, 21, 0] , distance: 50 , load: 170
Process finished with exit code 0
6.相关阅读
【运筹优化】带时间窗约束的车辆路径规划问题（VRPTW）详解 + Python 调用 Gurobi 建模求解
干货|十分钟快速掌握CPLEX求解VRPTW数学模型（附JAVA代码及CPLEX安装流程）
相关视频会在后续发布，欢迎关注我的bilibili：无形忍者的个人空间-无形忍者个人主页-哔哩哔哩视频
腾讯云面向开发者汇聚海量精品云计算使用和开发经验，营造开放的云计算技术生态圈。
更多推荐
26
0- 0
已为社区贡献2条内容
所有评论(0)