devbox
Monodyee
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

【python】使用DOcplex解决CVRP问题【附源码】

作者：Monodyee | 2024-05-19 14:05:09

赞

踩

docplex

一、CVRP问题简介

CVRP（Capacitated Vehicle Routing Problem）是一种组合优化问题，通常用于优化配送或运输问题。该问题的目标是有效地安排一组车辆，以满足一系列客户的需求，同时考虑车辆的容量限制。

以下是CVRP的一般介绍和要素：

  • 客户需求： 问题涉及到一组客户，每个客户对应一个配送点，有特定的需求量。需求可以表示为每个客户需要的货物数量。

  • 车辆： 问题还涉及到一组车辆，每辆车都有一个特定的容量限制，表示它可以携带的货物总量。

  • 路线： 目标是确定每辆车的路线，使得所有客户的需求都得到满足，并且每辆车的容量限制不被超过。车辆从一个中心点（通常是仓库）出发，访问所有客户后返回中心点。

  • 距离或成本： 优化问题通常还涉及到最小化总行驶距离或成本。这可以包括考虑道路网络、交通状况等因素。

  • 优化目标： 通常，优化的目标是最小化总成本，其中成本可以是行驶距离、时间、或其他相关因素。

解决CVRP一般采用高效的算法，例如启发式算法、遗传算法、模拟退火等。这是一个经典的组合优化问题，在物流和运输领域有很多实际应用，能够帮助提高运输效率、降低成本。

二、DOcplex库介绍

DOcplex是IBM在Python中提供的用于数学优化建模和求解的库。它是IBM Decision Optimization CPLEX Modeling for Python的一部分。这个库专注于提供一种方便的方式来定义和解决线性规划（LP）、整数规划（IP）和混合整数规划（MIP）等数学优化问题。

以下是DOcplex库的一些主要特点和介绍：

  1. 与CPLEX集成： DOcplex库是与IBM CPLEX（Optimization Studio中的求解器之一）紧密集成的。CPLEX是一个高性能的数学优化引擎，专注于解决线性规划、整数规划和混合整数规划问题。

  2. 优雅的建模语言： DOcplex提供了一种优雅的建模语言，使用户能够以自然的方式描述优化问题。这种建模语言允许用户指定目标函数、约束条件和变量，并且在底层使用CPLEX进行求解。

  3. 支持多种问题类型： 除了线性规划和整数规划，DOcplex还支持二次规划（Quadratic Programming）和混合整数二次规划（MIQP）等问题类型。

  4. 集成了约束编程： 除了数学规划，DOcplex还集成了约束编程。这使得它成为一个强大的工具，可以处理多种类型的优化问题。

  5. 与Jupyter Notebook兼容： DOcplex库可以轻松集成到Jupyter Notebook中，支持交互式建模和求解。

  6. 社区和文档： DOcplex有一个积极的社区，提供了详细的文档和示例，以帮助用户更好地理解和使用这个库。

总的来说，cplex是一个商业求解器，以分支定界和割平面理论构建了最优化问题的精确求解算法。DOcplex是IBM公司为了cplex引擎的通用化发展而在cplex Python API的基础上进行二次封装得到的工具库。使用该库构建和解决问题的流程具有普适性，同时求解速度与直接调用cplex API大同小异。IBM公司对cplex API的使用以及DOcplex库的使用都建立了实例文件，详见：Examples of CPLEX - IBM Documentation，然而恕我直言，初学者需要花费大量的精力才能看懂。

三、使用DOcplex解决CVRP问题

为了学习DOcplex的使用，我决定将经典的CVRP问题用这玩意重新写一遍。

1、DOcplex 安装

第一步，需要安装cplex。可以直接去这里：Site，用教育邮箱注册一个账户然后下载安装。（教育免费用，免费万岁）

然后

pip install docplex

2、 新建文件，导入必要的包

  1. from docplex.mp.model import Model
  2. import time
  3. import numpy as np
  4. import pandas as pd

 3、创建问题类，导入数据

  1. class Problem:
  2.     def __init__(self, file_path):
  3.  
  4.         # node data file path
  5.         self.file_path = file_path
  6.  
  7.         # some parameter
  8.         # the number of vehicle
  9.         self.m = 20
  10.         # vehicle capacity
  11.         self.Q = 200
  12.         # node number
  13.         self.n = 0
  14.  
  15.         # some matrix
  16.         # node data
  17.         self.data = np.array([])
  18.         # node distance
  19.         self.distance = np.array([])
  20.         # demand matrix
  21.         self.demand = np.array([])
  22.  
  23.     def lord_data(self):
  24.         """
  25.         lord the node data
  26.         get demand and node number
  27.         :return: None
  28.         """
  29.         df = pd.read_excel(self.file_path)
  30.         self.data = df.values
  31.         self.demand = self.data[:, 3].copy().reshape(1, -1)[0]
  32.         self.n = len(self.data)
  33.  
  34.     def calcu_distance_method(self, i, j):
  35.         """
  36.         distance function
  37.         :param i: index of first node
  38.         :param j: index of second node
  39.         :return: distance between i and j
  40.         """
  41.         if i == j:
  42.             return 1000000
  43.         x1 = self.data[i, 1]
  44.         y1 = self.data[i, 2]
  45.         x2 = self.data[j, 1]
  46.         y2 = self.data[j, 2]
  47.         dis = np.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)
  48.         return dis
  49.  
  50.     def calcu_distance_matrix(self):
  51.         """
  52.         calculate distance matrix
  53.         :return: None
  54.         """
  55.         dis = np.zeros((len(self.data), len(self.data)))
  56.         for i in range(len(self.data)):
  57.             for j in range(len(self.data)):
  58.                 dis[i, j] = self.calcu_distance_method(i, j)
  59.         self.distance = dis

问题类主要用于在建模过程中关键参数的提供。

4、构建算法

严格按照线性规划的显式公式进行算法构建，包括构建目标函数，构建约束条件等。对于CVRP问题，目标函数是运输成本（运输距离）最小化，约束包括唯一服务约束、流量平衡约束、子回路消除约束、需求满足约束等。其中最难搞的是子回路消除约束，本文中使用通用的MTZ约束方法来消除子回路。算法对应的CVRP模型公式如下：

（1）目标函数

min \sum_{k=1}^{K}\sum_{i=0}^{N}\sum_{j=0}^{N}c_{ij}x_{ij}^k

（2）每个点只被一辆车服务，每个点均被服务

\sum_{k=1}^{K}\sum_{i=1}^{N}x_{ij}^k = 1, \forall j = 1,...,N;

（3）MTZ消除子回路约束

\mu_{i}^{k}-\mu_{j}^{k} + Mx_{ij}^k \leq M-1,\forall i = 1,...,N; \forall j = 1,...,N; \forall k = 1,...,K;

（4）到达n点的路径数等于从n点出发的路径数

\sum_{i=1}^{N}x_{in}^k =\sum_{i=1}^{N}x_{nj}^k, \forall n = 1,...,N; \forall k =1,...,K;

（5）每条路径的终点必须是原点

\sum_{i=1}^{N}x_{i0}^k = 1, \forall k =1,...,N;

（6）每条路径的起点必须是原点

\sum_{j=1}^{N}x_{0j}^k = 1, \forall k =1,...,N;

（7）每条路径上的需求总量不超过车最大载重

\sum_{j=1}^Nd_jx_{ij}^k \leq Q, \forall i = 1,...,N; \forall k = 1,...,K;

算法构建如下，对应的约束和目标函数都加了注释，虽然是我蹩脚的英文注释。

  1. def create_program(problem):
  2.     """
  3.     create program and solve the problem
  4.     :param problem: problem object
  5.     :return: None
  6.     """
  7.  
  8.     # instantiate
  9.     model = Model("cvrp")
  10.  
  11.     # define variables
  12.     X = model.binary_var_cube(np.arange(problem.m), np.arange(problem.n), np.arange(problem.n), name="X")
  13.     mu = model.integer_var_matrix(np.arange(problem.m), np.arange(problem.n), name="Mu", lb=0, ub=100)
  14.  
  15.     # define objective function
  16.     sum_dis = 0
  17.     for k in range(problem.m):
  18.         for i in range(problem.n):
  19.             for j in range(problem.n):
  20.                 if i == j:
  21.                     sum_dis += 100000
  22.                 sum_dis += problem.distance[i][j] * X[(k, i, j)]
  23.     model.minimize(sum_dis)
  24.  
  25.     # another method
  26.     # model.maximize(model.sum(problem.distance[i][j] * X[(k, i, j)] for k in range(problem.m) for i in range(problem.n)
  27.     #                          for j in range(problem.n)))
  28.  
  29.     # constraint 1: each node can only be accessed once
  30.     for i in range(1, problem.n):
  31.         cons = 0
  32.         for k in range(problem.m):
  33.             for j in range(problem.n):
  34.                 cons += X[(k, i, j)]
  35.         model.add_constraint(cons == 1)
  36.  
  37.     # constraint 2: MTZ constraint
  38.     for k in range(problem.m):
  39.         for i in range(problem.n):
  40.             for j in range(problem.n):
  41.                 if i != j and i != 0 and j != 0:
  42.                     model.add_indicator(X[(k, i, j)], mu[(k, j)] - mu[(k, i)] == 1)
  43.  
  44.     # constraint 3: path constraint
  45.     for k in range(problem.m):
  46.         for n in range(1, problem.n):
  47.             cons1 = 0
  48.             cons2 = 0
  49.             for i in range(problem.n):
  50.                 if i != n:
  51.                     cons1 += X[(k, i, n)]
  52.             for j in range(problem.n):
  53.                 if j != n:
  54.                     cons2 += X[(k, n, j)]
  55.             model.add_constraint(cons1 == cons2)
  56.  
  57.     # constraint 4: origin constraint
  58.     for k in range(problem.m):
  59.         cons = 0
  60.         for j in range(1, problem.n):
  61.             cons += X[(k, 0, j)]
  62.         model.add_constraint(cons == 1)
  63.  
  64.     # constraint 5: destination constraint
  65.     for k in range(problem.m):
  66.         cons = 0
  67.         for i in range(1, problem.n):
  68.             cons += X[(k, i, 0)]
  69.         model.add_constraint(cons == 1)
  70.  
  71.     # constraint 6: capacity constraint
  72.     for k in range(problem.m):
  73.         cons = 0
  74.         for i in range(problem.n):
  75.             for j in range(1, problem.n):
  76.                 cons += X[(k, i, j)] * problem.demand[j]
  77.         model.add_constraint(cons <= problem.Q)
  78.     start_time = time.time()
  79.     sol = model.solve()
  80.     end_time = time.time()
  81.     print(sol)
  82.     print("spend time:", end_time - start_time)

5、求解

  1. if __name__ == "__main__":
  2.     # define the problem
  3.     p = Problem("cvrp.xlsx")
  4.     p.lord_data()
  5.     p.calcu_distance_matrix()
  6.  
  7.     # create the solving program
  8.     create_program(p)

6、数据

自己新建一个excel文件，取名“cvrp.xlsx”，把这一堆复制进去，文件放在与代码同一个目录下。（够详细了吧）

idx_coordy_coorddemand
140500
2258520
3227530
4228510
5208040
6208520
7187520
8157520
9158010
10103520
11104030
1284040
1384520
1453510
1554510
1624020
1704020
1804520
1944520
20421040
21421510
2240510
23401540
2438530
25381510
2635520
27953030
28953520
29923010
30903510
31883010
32883520
33873010
34852510
35853530
36678520
37658540
38658210
39628030
40608010
41608530
42587520
43558010
44558520
45558210
46208210
47188010
4824510
4942510
50421210
51723530
52552019
5325303
5420505
55556016
56306016
57503519
58302523
59151020
60102019
61156017
6245659
6365353
6465206
65453017
66354016
67413716
6864429
69406021
70315227
71356923
72655514
7363658
742605
7520208
765516
77601231
782337
7985627
8066830
81474713
82495810
8327439
84373114
85572918
8663232
87212428
88122413
89245819
9067525
9137476
92494213
93534314
9461523
95574823
9656376
97555426
9841835
9926529
100263515
10131673

7、结果

solution for: cvrp
objective: 2.02002e+08
status: OPTIMAL_SOLUTION(2)

略
spend time: 186.94465684890747

8、分析

首先，用该代码求出的是精确解，在论文中比用启发式算法认可度高了不知道多少倍。其次，100个顾客点的CVRP问题才跑了187秒，速度还是可以的。要知道很多论文的算法算例特别小，就是因为他们用的精确算法，而本身模型复杂度贼高，导致算例大了就跑不出来（可以忍受的时间之内）。总结，cplex牛批，省去我手撸单纯形法写分支定界了。

四、展望

下期预告，过去的作业，手撸单纯形法。

五、附完整代码

  1. from docplex.mp.model import Model
  2. import time
  3. import numpy as np
  4. import pandas as pd
  5.  
  6.  
  7. # create by Grasshopper9527
  8. class Problem:
  9.     def __init__(self, file_path):
  10.  
  11.         # node data file path
  12.         self.file_path = file_path
  13.  
  14.         # some parameter
  15.         # the maximum number of vehicle
  16.         self.m = 20
  17.         # vehicle capacity
  18.         self.Q = 200
  19.         # node number
  20.         self.n = 0
  21.  
  22.         # some matrix
  23.         # node data
  24.         self.data = np.array([])
  25.         # node distance
  26.         self.distance = np.array([])
  27.         # demand matrix
  28.         self.demand = np.array([])
  29.  
  30.     def lord_data(self):
  31.         """
  32.         lord the node data
  33.         get demand and node number
  34.         :return: None
  35.         """
  36.         df = pd.read_excel(self.file_path)
  37.         self.data = df.values
  38.         self.demand = self.data[:, 3].copy().reshape(1, -1)[0]
  39.         self.n = len(self.data)
  40.  
  41.     def calcu_distance_method(self, i, j):
  42.         """
  43.         distance function
  44.         :param i: index of first node
  45.         :param j: index of second node
  46.         :return: distance between i and j
  47.         """
  48.         if i == j:
  49.             return 1000000
  50.         x1 = self.data[i, 1]
  51.         y1 = self.data[i, 2]
  52.         x2 = self.data[j, 1]
  53.         y2 = self.data[j, 2]
  54.         dis = np.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)
  55.         return dis
  56.  
  57.     def calcu_distance_matrix(self):
  58.         """
  59.         calculate distance matrix
  60.         :return: None
  61.         """
  62.         dis = np.zeros((len(self.data), len(self.data)))
  63.         for i in range(len(self.data)):
  64.             for j in range(len(self.data)):
  65.                 dis[i, j] = self.calcu_distance_method(i, j)
  66.         self.distance = dis
  67.  
  68.  
  69. def create_program(problem):
  70.     """
  71.     create program and solve the problem
  72.     :param problem: problem object
  73.     :return: None
  74.     """
  75.  
  76.     # instantiate
  77.     model = Model("cvrp")
  78.  
  79.     # define variables
  80.     X = model.binary_var_cube(np.arange(problem.m), np.arange(problem.n), np.arange(problem.n), name="X")
  81.     mu = model.integer_var_matrix(np.arange(problem.m), np.arange(problem.n), name="Mu", lb=0, ub=100)
  82.  
  83.     # define objective function
  84.     sum_dis = 0
  85.     for k in range(problem.m):
  86.         for i in range(problem.n):
  87.             for j in range(problem.n):
  88.                 if i == j:
  89.                     sum_dis += 100000
  90.                 sum_dis += problem.distance[i][j] * X[(k, i, j)]
  91.     model.minimize(sum_dis)
  92.  
  93.     # another method
  94.     # model.maximize(model.sum(problem.distance[i][j] * X[(k, i, j)] for k in range(problem.m) for i in range(problem.n)
  95.     #                          for j in range(problem.n)))
  96.  
  97.     # each node can only be accessed once
  98.     for i in range(1, problem.n):
  99.         cons = 0
  100.         for k in range(problem.m):
  101.             for j in range(problem.n):
  102.                 cons += X[(k, i, j)]
  103.         model.add_constraint(cons == 1)
  104.  
  105.     # i != j
  106.     """for k in range(problem.m):
  107.         for i in range(problem.n):
  108.             for j in range(problem.n):
  109.                 if i == j:
  110.                     model.add_constraint(X[(k, i, j)] == 0)"""
  111.     # MTZ constraint
  112.     for k in range(problem.m):
  113.         for i in range(problem.n):
  114.             for j in range(problem.n):
  115.                 if i != j and i != 0 and j != 0:
  116.                     model.add_indicator(X[(k, i, j)], mu[(k, j)] - mu[(k, i)] == 1)
  117.     """for k in range(problem.m):
  118.         for i in range(problem.n):
  119.             for j in range(problem.n):
  120.                 if i != j and i != 0 and j != 0:
  121.                     model.add_indicator(X[(k, i, j)], mu[(k, i)] - mu[(k, j)] + problem.n * X[(k, i, j)]
  122.                                         <= problem.n - 1)"""
  123.  
  124.     # path constraint
  125.     for k in range(problem.m):
  126.         for n in range(1, problem.n):
  127.             cons1 = 0
  128.             cons2 = 0
  129.             for i in range(problem.n):
  130.                 if i != n:
  131.                     cons1 += X[(k, i, n)]
  132.             for j in range(problem.n):
  133.                 if j != n:
  134.                     cons2 += X[(k, n, j)]
  135.             model.add_constraint(cons1 == cons2)
  136.  
  137.     # origin constraint
  138.     for k in range(problem.m):
  139.         cons = 0
  140.         for j in range(1, problem.n):
  141.             cons += X[(k, 0, j)]
  142.         model.add_constraint(cons == 1)
  143.  
  144.     # destination constraint
  145.     for k in range(problem.m):
  146.         cons = 0
  147.         for i in range(1, problem.n):
  148.             cons += X[(k, i, 0)]
  149.         model.add_constraint(cons == 1)
  150.  
  151.     # capacity constraint
  152.     for k in range(problem.m):
  153.         cons = 0
  154.         for i in range(problem.n):
  155.             for j in range(1, problem.n):
  156.                 cons += X[(k, i, j)] * problem.demand[j]
  157.         model.add_constraint(cons <= problem.Q)
  158.     start_time = time.time()
  159.     sol = model.solve()
  160.     end_time = time.time()
  161.     print(sol)
  162.     print("spend time:", end_time - start_time)
  163.  
  164.  
  165. if __name__ == "__main__":
  166.     # define the problem
  167.     p = Problem("cvrp.xlsx")
  168.     p.lord_data()
  169.     p.calcu_distance_matrix()
  170.  
  171.     # create the solving program
  172.     create_program(p)

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/Monodyee/article/detail/593205
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号