建模实战|第二期:快速掌握ortools求解VRPTW数学模型(附python代码)

经过上一篇无形忍者-禁忌搜索算法求解带时间窗的车辆路径规划问题详解(附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;

⑤车辆需要在时间窗[a_i,b_i]  内服务客户i,并且需要在[a_0,b_0] 时间窗内返回仓库点

因此,VRPTW需要考虑更多约束条件和限制,使得问题更具挑战性。

2.VRPTW问题建模

2.1问题背景

VRPTW 问题是不仅考虑CVRP的所有约束,还需要考虑时间窗约束,也就是每个顾客对应一个时间窗[a_i,b_i] ,其中 和a_ib_i 分别代表该点的最早到达时间和最晚到达时间。通过下面介绍的参数说明和决策变量,结合参考文献(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         90

5.运行结果

-------------------- 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:无形忍者的个人空间-无形忍者个人主页-哔哩哔哩视频

Python是一种功能强大的编程语言,它可以用于解决各种问题,包括VRPTW(Vehicle Routing Problem with Time Windows)。VRPTW是指在考虑供应商、司机和用户的时间窗口约束下,通过合理调度车辆来完成物流配送的问题。 Python有许多优秀的库和算法可以用于解决VRPTW。以下是使用Python求解VRPTW的一般步骤: 1. 数据准备:首先,需要收集供应商、司机和用户的相关信息,例如位置坐标、货物数量、时间窗口等。这些数据可以保存在Excel、CSV或其他格式的文件中。 2. 数据读取:使用Python的pandas库或其他文件读取库,将准备好的数据文件导入到Python中,并存储为适当的数据结构,例如DataFrame或列表。 3. 算法选择:根据问题的规模和复杂性,选择适合的算法。常用的VRPTW算法包括遗传算法、模拟退火算法和粒子群算法Python中有许多优秀的开源库,例如DEAP、PyGMO和Particle Swarm Optimization,可以用于实现这些算法。 4. 编码实现:根据选择的算法,使用Python编写相应的代码实现。这可能包括定义适应度函数、编写遗传算子或迭代过程等。 5. 优化求解:运行算法,通过不断迭代和调整参数,寻找最佳的物流路径和调度方案。这个过程可能会消耗一定的计算资源和时间,但Python的高效性和可扩展性使得它成为求解VRPTW的理想选择。 6. 结果分析:在求解完成后,使用Python的可视化库如matplotlib或seaborn,将求解结果呈现出来。这样可以更好地理解和分析最佳路径和调度方案的有效性。 总之,Python求解VRPTW问题的理想工具之一。它提供了丰富的库和算法,使得这一问题的求解和分析变得更加高效和便捷。使用Python,我们可以快速而准确地找到最佳的物流路径和调度方案,帮助企业提高运输效率、降低成本。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值