车辆路径问题:精确和启发式解法

车辆路径问题:解法

了解如何使用Python解决复杂的路径规划问题

Photo by Nik Shuliahin 💛💙 on Unsplash

车辆路径规划问题(VRP)旨在确定一组最佳路线,由一组车辆为给定的一组客户提供服务。由于其多个应用和具有挑战性的组合方面,它是运筹学和数值优化中最研究的问题之一。

容量车辆路径规划问题(CVRP)是其中最常见的变体之一,因为它引入了具有有限载荷能力和可能的持续时间/距离约束的车辆。其他常见的变体还引入了多个配送中心、异构车队、取送货和时间窗口约束。

这些问题的组合方面使得在考虑一个简单的15个点的集合时,有6 × 10¹¹条可能的连接它们的路线(Dantzig&Ramser,1959)。因此,一些现实世界的应用将在计算和算法研究过去几十年的进展之前仍然不实际。分支-剪切-定价算法已经能够证明具有几百个客户的CVRP实例的最优性(Fukasawa等,2006;Pecin等,2017),而与局部搜索技术结合的最先进的元启发式算法可以在几秒钟内为这些实例提供高质量(有时是最优)的解决方案(Vidal等,2012;Vidal,2022)。

在本文中,我们将介绍具有负载(和持续时间)约束的容量车辆路径规划问题,并使用混合整数规划(MIP)和专门的(元)启发式算法来解决它。在第一部分中,我们将使用Python AML Pyomo和HiGHS求解器,而在第二部分中,我们将使用Google OR Tools包。

对于那些更感兴趣于实际应用而不是问题的理论方面的读者,可以快速浏览MIP部分,并更加关注专门的(元)启发式算法和有用的扩展部分。

那些对MIP表述有详细了解但对数值优化尚不熟悉的读者,在继续阅读本文之前,可以参考我之前关于线性规划和分支限界方法的故事。

如往常一样,您可以在此git存储库中找到完整的代码。

现在,让我们开始吧!

混合整数规划

本节中介绍的数学公式将使用Toth&Vigo(2002)在其模型中所提到的“三指标车流公式”中所呈现的相同方程。

考虑一个节点集合V(需求和配送中心)和一个车辆集合K。我们将使用小写字母i和j表示节点索引,小写字母k表示车辆索引。由于该模型适用于非对称情况,让我们假设节点是完全定向图G(V,A)的一部分,其中A是弧。在这个问题中,有一个由0索引的单个配送中心节点,所有车辆都具有相同的容量Q。考虑两组决策变量:

  • x_{i, j, k}:是一个二进制变量,表示由车辆k执行的从节点i到节点j的活动弧。
  • y_{i, k}:是一个二进制变量,表示节点i的需求由车辆k满足。

我们的目标是最小化与活动弧相关联的成本值。总持续时间或距离是常见的示例。假设穿越弧i,j的成本是cᵢⱼ。目标函数可以表示如下。

容量车辆路径规划问题的目标函数。 (作者提供的图像)

我们还需要包括一些约束条件,以确保:

  • 每个客户i只被访问一次,因此具有从它出发的一条活动弧和到达它的一条活动弧。
  • 如果任何由车辆k索引的弧变量进入节点i或从节点i出去,该节点的需求q将分配给车辆k。
  • 分配给车辆的总需求不得超过其容量Q。
  • 正好有|K|个节点从配送中心出发并到达配送中心。
  • 不存在子路径…然而,子路径的数量可能太多,无法从一开始就枚举出来。我们将详细介绍如何做到这一点。
CVRP的约束条件。 (图片由作者提供)

与Python教程一样,让我们从导入本节中使用的库开始:

import timefrom itertools import cycleimport numpy as npfrom scipy.spatial.distance import pdist, squareformimport matplotlib.pyplot as pltimport matplotlib as mplimport networkx as nximport pyomo.environ as pyofrom pyomo.contrib.appsi.solvers.highs import Highs

现在,让我们用N个需求节点实例化一个随机问题。在这个例子中,假设仓库节点是第一个节点(索引为0),所以我们确保它对应的需求也是零。

np.random.seed(42)  # 结果应始终相同N = 10demands = np.random.randint(1, 10, size=N)demands[0] = 0capacity = 15n_vehicles = 4coordinates = np.random.rand(N, 2)distances = squareform(pdist(coordinates, metric="euclidean"))distances = np.round(distances, decimals=4)  # 避免数值误差

可以通过使用装箱问题来计算所需车辆的数量。完整源代码中还包含了如何执行此操作的示例。

在pyomo中,有两种建模问题的方法:抽象模型和具体模型。在第一种方法中,问题的代数表达式在提供一些数据值之前被定义,而在第二种方法中,模型实例在其元素被定义时立即创建。您可以在库文档或Bynum等人的书中了解更多关于这些方法的信息(2021年)。在本文中,我们将采用具体模型的形式。

model = pyo.ConcreteModel()

让我们实例化需求节点V、弧A和车辆K的集合。请注意,仓库节点包括在节点V的集合中,就像在原始的数学公式中一样。

model.V = pyo.Set(initialize=range(len(demands)))model.A = pyo.Set(initialize=[(i, j) for i in model.V for j in model.V if i != j])model.K = pyo.Set(initialize=range(n_vehicles))

现在是我们的容量、需求负载和弧成本的参数。

model.Q = pyo.Param(initialize=capacity)model.q = pyo.Param(model.V, initialize={i: d for (i, d) in enumerate(demands)})model.c = pyo.Param(model.A, initialize={(i, j): distances[i, j] for (i, j) in model.A})

我们还必须包括先前列出的约束条件。首先,让我们使用通常的Pyomo签名来实现它们:function(model, *domain)。

def arcs_in(model, i):    if i == model.V.first():        return sum(model.x[:, i, :]) == len(model.K)    else:        return sum(model.x[:, i, :]) == 1.0def arcs_out(model, i):    if i == model.V.first():        return sum(model.x[i, :, :]) == len(model.K)    else:        return sum(model.x[i, :, :]) == 1.0def vehicle_assignment(model, i, k):    return sum(model.x[:, i, k]) == model.y[i, k]def comp_vehicle_assignment(model, i, k):    return sum(model.x[i, :, k]) == model.y[i, k]def capacity_constraint(model, k):    return sum(model.y[i, k] * model.q[i] for i in model.V) <= model.Q

然后将它们合并到我们的模型中。

model.arcs_in = pyo.Constraint(model.V, rule=arcs_in)model.arcs_out = pyo.Constraint(model.V, rule=arcs_out)model.vehicle_assignment = pyo.Constraint(model.V, model.K, rule=vehicle_assignment)model.comp_vehicle_assignment = pyo.Constraint(model.V, model.K, rule=comp_vehicle_assignment)model.capacity_constraint = pyo.Constraint(model.K, rule=capacity_constraint)

请注意,我还没有包含消除子路径约束。我们应该考虑将N个节点中的每k个节点的所有可能排列,即使对于中等规模的实例来说,这些排列也可能过大以至于无法枚举。作为替代,我们将在找到新解时递归地包含消除子路径约束,如果我们验证该解产生了子路径。在一些商业求解器中,这些被称为“惰性约束”,可以通过回调直接将其纳入求解器。

首先,让我们创建一个函数,给定一个子路径、所有剩余节点、一个来自子路径的节点和一个车辆,返回与之前所述的数学公式相对应的Pyomo表达式。同时,让我们包括一个ConstraintList,我们将在解决方案中不断添加新元素。

def subtour_elimination(model, S, Sout, h, k):    nodes_out = sum(model.x[i, j, k] for i in S for j in Sout)    return model.y[h, k] <= nodes_outmodel.subtour_elimination = pyo.ConstraintList()

我们必须创建一些函数,它们将根据解决方案返回创建的子路径(如果存在)。为此,我们将首先使用find_arcs函数创建一个模型中活动弧的列表。然后,我们将使用Networkx DiGraph类创建一个不完整的有向图。find_subtours函数应返回一个连接组件集的列表。

def find_arcs(model):    arcs = []    for i, j in model.A:        for k in model.K:            if np.isclose(model.x[i, j, k].value, 1):                arcs.append((i, j))    return arcsdef find_subtours(arcs):    G = nx.DiGraph(arcs)    subtours = list(nx.strongly_connected_components(G))    return subtours

我们的目标是消除不包括集散地节点的连接组件组。因此,在下一步中,我们将创建一些函数,迭代集合列表,并在组件集合不包括集散地节点时包括新的约束条件。这将使用ConstraintList类的add方法。

def eliminate_subtours(model, subtours):    proceed = False    for S in subtours:        if 0 not in S:            proceed = True            Sout = {i for i in model.V if i not in S}            for h in S:                for k in model.K:                    model.subtour_elimination.add(                      subtour_elimination(model, S, Sout, h, k)                    )    return proceed

现在,我们已经准备好提出一个解决方案程序,该程序通过迭代求解MIP,验证当前解决方案是否存在子路径,并在存在子路径时包括新的约束条件以消除它们。否则,找到的解决方案是最优的。

def solve_step(model, solver):    sol = solver.solve(model)    arcs = find_arcs(model)    subtours = find_subtours(arcs)    time.sleep(0.1)    proceed = eliminate_subtours(model, subtours)    return sol, proceed def solve(model, solver):    proceed = True    while proceed:        sol, proceed = solve_step(model, solver)    return sol

现在,让我们实例化求解器并解决我们的模型。如果已安装highspy软件包,则Pyomo中可用Highs求解器(请检查导入)。因此,请确保运行pip install highspy

solver = Highs()solver.highs_options = {    "log_file": "Highs.log",    "mip_heuristic_effort": 0.2,    "mip_detect_symmetry": True,    "mip_rel_gap": 1e-6,}solution = solve(model, solver)

再编写一个函数来查找创建的旅行路线,然后我们就可以准备绘制结果了。

def find_tours(model):    tours = []    for k in model.K:        node = 0        tours.append([0])        while True:            for j in model.V:                if (node, j) in model.A:                    if np.isclose(model.x[node, j, k].value, 1):                        node = j                        tours[-1].append(node)                        break            if node == 0:                break    return tours
Tours produced in CVRP using MIP. (Image by the author).

对于包含10个节点的小实例,结果非常出色。然而,即使对于这样一个小实例,求解器也花费了将近半分钟来获得解决方案,并且随着需求点的增加,复杂度显著增加。幸运的是,现有专门的算法可以在短时间内为更大的实例找到高质量的解决方案。让我们在下一节中来看一下它们。

专用(元)启发式算法

多年来,已经提出了几种专用的(元)启发式算法来解决VRP的变体问题。其中大部分算法都强烈依赖于局部搜索算法,因此在给定一个解决方案后,尝试不同的扰动以依次改进其成本,直到在给定邻域中无法进一步改进为止。使用Google OR工具时,我们也将使用与构建算法相关的局部搜索方法。

在本节中,将使用Rochat和Taillard(1995)的第150d实例。其数据来自CVRPLIB。该实例有150个客户和一个仓库节点,所以我们肯定无法使用之前介绍的MIP策略来解决它。

让我们再次从导入所使用的库开始。

from itertools import cycleimport numpy as npimport pandas as pdfrom scipy.spatial.distance import pdist, squareformimport matplotlib.pyplot as pltimport matplotlib as mplfrom ortools.constraint_solver import routing_enums_pb2from ortools.constraint_solver import pywrapcp

然后,让我们从包含每个节点坐标和需求的文件中加载问题数据。

dataset = pd.read_csv("./data/tai150d.csv", index_col=0)coordinates = dataset.loc[:, ["x", "y"]]demands = dataset.d.valuescapacity = 1874n_vehicles = 15N = coordinates.shape[0]distances = squareform(pdist(coordinates, metric="euclidean"))distances = np.round(distances, decimals=4)

使用OR Tools VRP求解器的第一步是实例化路由管理器和模型。

# 创建路由索引管理器:节点数量,车辆数量,仓库节点manager = pywrapcp.RoutingIndexManager(    N, n_vehicles, 0)# 创建路由模型routing = pywrapcp.RoutingModel(manager)

接下来,我们将包括回调函数来量化与弧/边和节点相关的维度。我们的路由实例的RegisterTransitCallback方法可以用于量化与弧/边相关的任何维度,而RegisterUnaryTransitCallback方法可以量化与节点相关的值。

# 与弧/边相关的任何回调函数都可以使用相同的有效方法def distance_callback(from_index, to_index):    from_node = manager.IndexToNode(from_index)    to_node = manager.IndexToNode(to_index)    return distances[from_node, to_node]transit_callback_index = routing.RegisterTransitCallback(distance_callback)# 与节点相关的任何回调函数都可以使用相同的有效方法def demand_callback(from_index):    from_node = manager.IndexToNode(from_index)    return demands[from_node]demand_callback_index = routing.RegisterUnaryTransitCallback(demand_callback)

现在,我们将使用之前定义的demand_callback_index包括一个容量约束。注意,可以使用相同的语法定义持续时间约束,只需将RegisterTransitCallback的实例作为第一个参数传递即可。此外,重要的是要注意,路由模型处理异构车队,因此我们必须在第三个参数中传递一个值的列表。

# 与车辆相关的任何约束都可以使用相同的参数routing.AddDimensionWithVehicleCapacity(    demand_callback_index,    0,  # 空容量    [capacity,] * n_vehicles,  # 每个车辆的最大容量(每个车辆的列表)    True,  # 开始累积为零    'Capacity')

类似地,目标的定义也需要一个回调函数作为主要参数。在这个例子中,让我们最小化在transit_callback_index中定义的距离。

routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

最后,我们必须定义求解器参数并解决我们的模型。

# 设置启发式策略search_parameters = pywrapcp.DefaultRoutingSearchParameters()search_parameters.first_solution_strategy = (    routing_enums_pb2.FirstSolutionStrategy.CHRISTOFIDES)search_parameters.local_search_metaheuristic = (    routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)search_parameters.time_limit.FromSeconds(300)# 解决问题solution = routing.SolveWithParameters(search_parameters)

以下代码片段可用于提取我们解决方案中使用的路线。

tours = []for vehicle_id in range(n_vehicles):    index = routing.Start(vehicle_id)    tours.append([])    while not routing.IsEnd(index):        node_index = manager.IndexToNode(index)        previous_index = index        index = solution.Value(routing.NextVar(index))        tours[-1].append(node_index)    else:        node_index = manager.IndexToNode(index)        tours[-1].append(node_index)

可以通过运行以下简单代码行solution.ObjectiveValue()来访问目标值。使用提供的配置,我得到了一个目标值为2679,这与已知的最优值2645非常接近(1.2%的差距)。下面是所得到的路径的表示。

在tai150d实例中使用ortools获得的路径。(图片由作者提供)

完整的代码(包括图表)可以在此git存储库中获得。

有用的扩展

OR Tools库非常适用于解决各种变体的VRP问题,如时间窗口、异构车队和多个配送点。然而,适用于经典CVRP的算法可能表现得更好。值得一看的是HGS-CVRP软件包(Vidal, 2022),它将最先进的遗传算法与多个局部搜索策略相结合。该算法能够在不到20秒内找到tai150d实例的最优解。

在一些实际应用中,连接位置的距离(或持续时间)可能更适合使用道路距离而不是欧几里得距离。Google提供了一个很好的付费接口,您可以在本教程中了解相关信息。然而,如果您正在寻找开源替代方案,值得检查OpenStreetMap API。以下是一些有用的请求:

  • https://router.project-osrm.org/table/v1/driving/<LOCATIONS>
  • http://router.project-osrm.org/route/v1/car/<LOCATIONS>?overview=false&steps=true

其中<LOCATIONS>应为由逗号分隔的经度和纬度对列表,并且不同对之间用分号分隔。您还可以在表格请求中指定源和目标,这在完整表格太大以至于无法处理的情况下非常有用。

除了进行精确的路径计算之外,可视化也是一个重要的工具。Python库folium在这方面非常有用,值得一看。

进一步阅读

在本文的前面部分,我们实现了一个用于CVRP的精确MIP模型,但不适用于中等规模的实例。然而,将列生成与分支定界相结合的算法已经成功地解决了具有数百个客户的实例。值得阅读Fukasawa等人(2006)和Pecin等人(2017)的研究论文。

对于对列生成感兴趣的人,可以在我之前的VoAGI文章中找到相关介绍。

关于元启发式算法,Vidal等人(2012)和Vidal(2022)的论文非常出色。两者也可以作为技术报告获得,HGS-CVRP存储库中提供了链接。

结论

本文介绍了两种解决容量车辆路径问题(CVRP)的方法:混合整数规划和(元)启发式算法。第一种方法用于解决一个小规模实例,取得了成功,但无法处理中等规模或大规模实例。第二种方法用于解决文献中具有150个客户的具有挑战性的问题,求解器在300秒内找到了一个与已知最优解相差1.2%的高质量解。

参考文献

Bynum, M. L.等人,2021。Pyomo-optimization modeling in python. Springer。

Dantzig, G. B.和Ramser, J. H., 1959。The truck dispatching problem. Management science, 6(1), 80–91。

Fukasawa, R.,Longo, H.,Lysgaard, J.,Aragão, M. P. D.,Reis, M.,Uchoa, E.和Werneck, R. F., 2006。Robust branch-and-cut-and-price for the capacitated vehicle routing problem. Mathematical programming, 106, 491–511。

Pecin, D., Pessoa, A., Poggi, M., & Uchoa, E., 2017. 提高容量车辆路径规划的分支切割定价算法. 数学规划计算, 9, 61–100.

Rochat, Y., & Taillard, É. D., 1995. 车辆路径规划中的概率多样性和集中化的局部搜索. 启发式学报, 1, 147–167.

Toth, P., & Vigo, D., 2002. 车辆路径规划问题概述. 车辆路径规划问题, 1–26.

Vidal, T., 2022. CVRP的混合遗传搜索: 开源实现和SWAP*邻域. 计算与运筹学, 140, 105643.

Vidal, T., Crainic, T. G., Gendreau, M., Lahrichi, N., & Rei, W., 2012. 一种用于多仓库和周期性车辆路径规划问题的混合遗传算法. 运筹学, 60(3), 611–624.