OR-Tools VRP 问题从入门到升天(二) 从CVRP问题

CVRP问题

有容量限制的车辆路径规划问题(Capacitated Vehicle Routing Problem)是车辆路径规划问题的一类经典变体。在这类问题中，每个节点都有一个需求量，每辆车都有最大的负载，要求分配给每辆车的路径上，所有节点需求量之和都不超过车辆的最大负载。

这里我们用CVRP问题以及一个小变种为例，看一下在Ortools中，如何对目标函数进行自定义。

算例我们仍然采用TSP问题中的gr120.tsp，车辆数设置为4辆，并且随机生成每个节点的需求以及每辆车的起点和终点。

数据准备部分的代码如下：

## 导入库函数
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp
import numpy as np
import tsplib95
import matplotlib.pyplot as plt
%matplotlib inline

## 数据准备
# 给定数据目录并读入数据
fname: str = "gr120.tsp"
input_fpath: str = r"data/" + fname
data = tsplib95.load(input_fpath).as_name_dict()
# 让节点编号从0开始
data["display_data"] = {k-1: v for k, v in data["display_data"].items()}
# 补充每个节点的demand数据
rng = np.random.default_rng(seed = 0)
data["demand"] = rng.integers(1, 11, data["dimension"]) # 随机生成[1, 10]的需求
# 车辆数据
n_vehicles = 4 # 车辆数
capacity = [250] * n_vehicles # 车辆最大负载
vehicle_start_idx = rng.integers(0, data["dimension"], n_vehicles) # 每辆车的起点
# 注意ortools是c++写的，它并不接受 np.int64，因此要转化为python原生的int类型
vehicle_start_idx = [i.item() for i in vehicle_start_idx]
vehicle_end_idx = rng.integers(0, data["dimension"], n_vehicles) # 每辆车的终点
vehicle_end_idx = [i.item() for i in vehicle_end_idx]

在完成数据准备之后，可以通过下面代码对数据进行可视化：

# 节点可视化
x_coor: list = [i[0] for i in data['display_data'].values()]
y_coor: list = [i[1] for i in data['display_data'].values()]
plt.figure(figsize = (8, 6), dpi = 120)
plt.scatter(x_coor, y_coor)
idx = 1
offset = 2.5
for start, end in zip(vehicle_start_idx, vehicle_end_idx):
    plt.plot(data["display_data"][start][0], data["display_data"][start][1], "r*", markersize=16)
    plt.plot(data["display_data"][end][0], data["display_data"][end][1], "b*", markersize=16)
    plt.text(data["display_data"][start][0]+offset, data["display_data"][start][1]+offset, "{}".format(idx),fontsize=12)
    plt.text(data["display_data"][end][0]+offset, data["display_data"][end][1]+offset, "{}".format(idx),fontsize=12)
    idx += 1
plt.xlabel("X coordinate", fontsize = 14)
plt.ylabel("Y coordinate", fontsize = 14)
plt.title("{} Nodes".format(fname), fontsize = 16)

如下图，在途中，红色星形表示车辆起点，蓝色星形表示车辆终点，其他圆形代表一般节点：

gr120.tsp_start_end_vrp.png

模型与求解

## 建立并求解CVRP问题
# 建立RoutingIndexManager
manager = pywrapcp.RoutingIndexManager(data["dimension"], n_vehicles, 
                             vehicle_start_idx, vehicle_end_idx)
# 实例化RoutingModel
routing = pywrapcp.RoutingModel(manager)
# 建立一个callback来记录车辆的总行驶里程
def distance_callback(from_index: int, to_index: int):
    from_node = manager.IndexToNode(from_index)
    to_node = manager.IndexToNode(to_index)
    from_coor = data['display_data'][from_node]
    to_coor = data['display_data'][to_node]
    return (to_coor[0] - from_coor[0])**2 + (to_coor[1] - from_coor[1])**2

transit_callback_index = routing.RegisterTransitCallback(distance_callback)
# 给每条边设置cost
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)
# 建立一个callback来记录车辆的负载
def demand_callback(from_index: int):
    # 经过一个节点之后，累加节点上的需求量
    from_node = manager.IndexToNode(from_index)
    return data['demand'][from_node]
demand_callback_index = routing.RegisterUnaryTransitCallback(demand_callback)
# 利用这个callback，施加capacity约束
routing.AddDimensionWithVehicleCapacity(
        demand_callback_index,
        0,  # 不允许松弛
        capacity,  # 车辆的最大允许负载
        True,  # 从0开始进行累积
        'Capacity')
# 设定求解参数
search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.first_solution_strategy = (
    routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)
search_parameters.local_search_metaheuristic = (
    routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)
search_parameters.time_limit.FromSeconds(100) # 设置最大求解时间100秒

# 求解
solution = routing.SolveWithParameters(search_parameters)

# 输出结果
print(f'Objective: {solution.ObjectiveValue()}')
total_distance = 0
total_load = 0
paths = []
for vehicle_id in range(n_vehicles):
    index = routing.Start(vehicle_id)
    plan_output = '车辆 {} 的路径:\n'.format(vehicle_id)
    route_distance = 0
    route_load = 0
    path: list = []
    while not routing.IsEnd(index):
        node_index = manager.IndexToNode(index)
        path.append(node_index)
        route_load += data['demand'][node_index]
        plan_output += ' {0} 装载({1}) -> '.format(node_index, route_load)
        previous_index = index
        index = solution.Value(routing.NextVar(index))
        route_distance += routing.GetArcCostForVehicle(
            previous_index, index, vehicle_id)
    plan_output += ' {0} 装载({1})\n'.format(manager.IndexToNode(index),
                                             route_load)
    path.append(manager.IndexToNode(index))
    plan_output += '路径总长度: {}m\n'.format(route_distance)
    plan_output += '路径总负载: {}\n'.format(route_load)
    print(plan_output)
    paths.append(path)
    total_distance += route_distance
    total_load += route_load
print('所有路径总长度: {}m'.format(total_distance))
print('所有路径总负载: {}'.format(total_load))

得到解如下：

Objective: 33045
车辆 0 的路径:
 68 装载(9) ->  113 装载(9)
路径总长度: 0m
路径总负载: 9

车辆 1 的路径:
 48 装载(3) ->  12 装载(9) ->  117 装载(12) ->  112 装载(21) ->  64 装载(30) ->  67 装载(38) ->  90 装载(42) ->  57 装载(49) ->  78 装载(55) ->  99 装载(64) ->  32 装载(65) ->  51 装载(69) ->  24 装载(73) ->  116 装载(81) ->  30 装载(90) ->  65 装载(92) ->  84 装载(100) ->  17 装载(106) ->  45 装载(113) ->  97 装载(114) ->  16 装载(121) ->  41 装载(122) ->  40 装载(127) ->  6 装载(128) ->  37 装载(131) ->  33 装载(140) ->  7 装载(141) ->  110 装载(147) ->  95 装载(154) ->  89 装载(158) ->  53 装载(168) ->  54 装载(177) ->  88 装载(185) ->  5 装载(186) ->  83 装载(190) ->  34 装载(191) ->  9 装载(200) ->  98 装载(204) ->  61 装载(211) ->  66 装载(217) ->  82 装载(220) ->  38 装载(225) ->  56 装载(229) ->  4 装载(233) ->  26 装载(239) ->  52 装载(244) ->  81 装载(250) ->  10 装载(250)
路径总长度: 17927m
路径总负载: 250

车辆 2 的路径:
 119 装载(1) ->  31 装载(3) ->  28 装载(11) ->  60 装载(20) ->  15 装载(28) ->  0 装载(37) ->  75 装载(46) ->  29 装载(54) ->  58 装载(64) ->  14 装载(74) ->  91 装载(83) ->  27 装载(84) ->  44 装载(85) ->  77 装载(95) ->  85 装载(101) ->  93 装载(104) ->  80 装载(111) ->  21 装载(120) ->  18 装载(126) ->  107 装载(127) ->  42 装载(128) ->  106 装载(129) ->  19 装载(139) ->  49 装载(146) ->  43 装载(148) ->  55 装载(158) ->  71 装载(162) ->  39 装载(167) ->  104 装载(175) ->  73 装载(180) ->  86 装载(186) ->  13 装载(193) ->  74 装载(193)
路径总长度: 6801m
路径总负载: 193

车辆 3 的路径:
 23 装载(1) ->  59 装载(8) ->  115 装载(16) ->  25 装载(25) ->  3 装载(28) ->  118 装载(36) ->  102 装载(40) ->  8 装载(42) ->  22 装载(49) ->  50 装载(57) ->  114 装载(66) ->  1 装载(73) ->  92 装载(76) ->  20 装载(79) ->  108 装载(86) ->  63 装载(90) ->  87 装载(94) ->  96 装载(95) ->  11 装载(105) ->  94 装载(113) ->  76 装载(114) ->  62 装载(122) ->  72 装载(127) ->  36 装载(128) ->  105 装载(137) ->  103 装载(140) ->  35 装载(146) ->  111 装载(148) ->  109 装载(152) ->  100 装载(157) ->  79 装载(161) ->  2 装载(167) ->  101 装载(175) ->  47 装载(182) ->  46 装载(188) ->  70 装载(192) ->  69 装载(192)
路径总长度: 8317m
路径总负载: 192

所有路径总长度: 33045m
所有路径总负载: 644

可以看到，这里我们的路径总长度和目标函数值是相等的。

用如下代码进行路径的可视化：

# 路径可视化
colors = ['r', 'g', 'b', 'k']
# 画出所有节点
x_coor: list = [i[0] for i in data['display_data'].values()]
y_coor: list = [i[1] for i in data['display_data'].values()]
plt.figure(figsize = (8, 6), dpi = 120)
plt.scatter(x_coor, y_coor)
# 画出每辆车的起点和终点
idx = 1
for start, end in zip(vehicle_start_idx, vehicle_end_idx):
    plt.plot(data["display_data"][start][0], data["display_data"][start][1], "r*", markersize=16)
    plt.plot(data["display_data"][end][0], data["display_data"][end][1], "b*", markersize=16)
    plt.text(data["display_data"][start][0]+offset, data["display_data"][start][1]+offset, "{}".format(idx),fontsize=12)
    plt.text(data["display_data"][end][0]+offset, data["display_data"][end][1]+offset, "{}".format(idx),fontsize=12)
    idx += 1
# 画出每辆车的路径
for k in range(n_vehicles):
    path = paths[k]
    for i, j in zip(path, path[1:]):
        start_coor = data['display_data'][i]
        end_coor = data['display_data'][j]
        plt.arrow(start_coor[0], start_coor[1], end_coor[0] - start_coor[0], 
                  end_coor[1] - start_coor[1], fc=colors[k], ec=colors[k])
plt.xlabel("X coordinate", fontsize = 14)
plt.ylabel("Y coordinate", fontsize = 14)
plt.title("CVRP paths", fontsize = 16)
plt.savefig("../pics/vrp_paths_without_load_balance.png")

得到：

vrp_paths_without_load_balance.png

车辆负载的平衡

在上面的求解例子中，可以看到车辆之间的负载是恨不均衡的， 0号车只经过了起点和终点，总共负载为9，而车辆1的负载则达到了250。在实际使用时，我们通常希望我们的车辆有相近的负载，不要让某些员工过量工作而某些员工整天闲逛。因此我们需要对目标函数进行修改，让负载的均衡也成为我们目标的一部分。

我们可以这样处理负载均衡的要求：

1. 通过累加所有节点上的需求量，得到每辆车平均需要承担的负载（这里更精细一点的处理应该去除掉depot的需求量，因为通常定义上我们不需要服务车库）；
2. 统计每辆车路径上的负载，当车辆负载超过平均负载时，施加一个惩罚函数

具体实现如下，demand_dimension上累加了车辆在路径上的负载，我们在车辆到达最后一个节点时，检查车辆负载是不是超过我们设定的平均值，如果出现超载，则通过penalty_coef施加一个惩罚项：

# 通过软约束来平衡车辆的负载
demand_dimension = routing.GetDimensionOrDie('Capacity')
penalty_coef: int = 1000
average_demand: int = sum(data['demand']).item()//n_vehicles
for vehicle_id in range(manager.GetNumberOfVehicles()):
    index = routing.End(vehicle_id)
    demand_dimension.SetCumulVarSoftUpperBound(index, average_demand, penalty_coef)

得到的结果如下：

Objective: 37059
车辆 0 的路径:
 68 装载(9) ->  42 装载(10) ->  107 装载(11) ->  18 装载(17) ->  84 装载(25) ->  65 装载(27) ->  30 装载(36) ->  116 装载(44) ->  24 装载(48) ->  51 装载(52) ->  32 装载(53) ->  99 装载(62) ->  78 装载(68) ->  57 装载(75) ->  90 装载(79) ->  67 装载(87) ->  64 装载(96) ->  112 装载(105) ->  12 装载(111) ->  117 装载(114) ->  16 装载(121) ->  97 装载(122) ->  55 装载(132) ->  86 装载(138) ->  73 装载(143) ->  39 装载(148) ->  25 装载(157) ->  46 装载(163) ->  109 装载(167) ->  113 装载(167)
路径总长度: 14260m
路径总负载: 167

车辆 1 的路径:
 48 装载(3) ->  41 装载(4) ->  40 装载(9) ->  6 装载(10) ->  22 装载(17) ->  50 装载(25) ->  114 装载(34) ->  92 装载(37) ->  63 装载(41) ->  52 装载(46) ->  4 装载(50) ->  56 装载(54) ->  38 装载(59) ->  82 装载(62) ->  66 装载(68) ->  36 装载(69) ->  61 装载(76) ->  98 装载(80) ->  103 装载(83) ->  9 装载(92) ->  34 装载(93) ->  83 装载(97) ->  35 装载(103) ->  5 装载(104) ->  88 装载(112) ->  54 装载(121) ->  89 装载(125) ->  95 装载(132) ->  53 装载(142) ->  70 装载(146) ->  47 装载(153) ->  101 装载(161) ->  81 装载(167) ->  10 装载(167)
路径总长度: 10055m
路径总负载: 167

车辆 2 的路径:
 119 装载(1) ->  28 装载(9) ->  60 装载(18) ->  15 装载(26) ->  0 装载(35) ->  75 装载(44) ->  58 装载(54) ->  14 装载(64) ->  29 装载(72) ->  31 装载(74) ->  91 装载(83) ->  27 装载(84) ->  44 装载(85) ->  77 装载(95) ->  85 装载(101) ->  93 装载(104) ->  80 装载(111) ->  21 装载(120) ->  17 装载(126) ->  106 装载(127) ->  19 装载(137) ->  45 装载(144) ->  49 装载(151) ->  43 装载(153) ->  13 装载(160) ->  74 装载(160)
路径总长度: 4011m
路径总负载: 160

车辆 3 的路径:
 23 装载(1) ->  110 装载(7) ->  7 装载(8) ->  59 装载(15) ->  104 装载(23) ->  71 装载(27) ->  37 装载(30) ->  102 装载(34) ->  8 装载(36) ->  1 装载(43) ->  20 装载(46) ->  108 装载(53) ->  87 装载(57) ->  96 装载(58) ->  11 装载(68) ->  94 装载(76) ->  76 装载(77) ->  62 装载(85) ->  72 装载(90) ->  105 装载(99) ->  111 装载(101) ->  100 装载(106) ->  79 装载(110) ->  26 装载(116) ->  2 装载(122) ->  118 装载(130) ->  3 装载(133) ->  33 装载(142) ->  115 装载(150) ->  69 装载(150)
路径总长度: 8733m
路径总负载: 150

所有路径总长度: 37059m
所有路径总负载: 644

可以看到，我们的平均负载是 670/4=167（输出的644是减去了4个车库之后的量，和我们代码中的处理稍有不同）。如果我们施加一个较小的平均负载值，我们就可以看到路径长度和目标函数之间的差异，差异的部分也就是惩罚项。

对路径进行可视化如下：

vrp_paths_with_load_balance.png

改变目标函数的方式

从上面这个负载均衡问题引申出来的问题就是：我们怎么样去在不同的需求下，对目标函数进行定制。

Ortools的VRP求解器中，并没有明确的定义目标函数。而是在通用的意义上，将最小化所有车辆在经过的边(arc)上的代价(cost)作为目标函数。

在实际应用中，比较常见的改变目标函数的方式有下面三种：

• 通过软约束在目标函数中添加惩罚项

  • RoutingDimension.SetCumulVarSoftUpperBound(index: int64, upper_bound: int64, coefficient: int64) - 为一个累积量设置一个上界，如果累积量的值超过上界，则在目标函数中施加一个和coefficient成比例的惩罚项，也就是说cumulVar <= upper_bound -> cost = 0，cumulVar > upper_bound -> cost = coefficient * (cumulVar - upper_bound)；这个方法还有对应的设置下界的版本。
  • RoutingModel.AddDisjunction(indices: List[int64], penalty = 0: int64, max_cardinality: int64 = 1) - 添加一个disjuction constraint，要求在传入的indices代表的节点中有max_cardinality个被访问到（注意这里在indices中，我们不能传入车辆的起点和终点），如果没有足够数量的节点被访问，那么我们的惩罚为 p*penalty，其中p=max_cardinality - sum(active[i])这里要求传入的penalty为正数，否则将会强制这些节点不会被访问；它的一个变种应用是当我们有一些节点可选时，我们将每个节点的index传入，这样当我们drop掉这个节点时，我们会受到惩罚，但是不再强制每个节点都会被访问。

• 通过系数设置在目标函数中添加自定义项

  • RoutingDimension.SetGlobalSpanCostCoefficient(coefficient: int64) - 这里所谓的Global，指的是全局所有车辆中该累积量在终点的最大值和在起点的最小值之差，也就是global_span_cost = coefficient * (Max(dimension end value) - Min(dimension start value))；在限定车辆服务的节点数时，我们常用这种方法。它还有两个类似的方法SetSpanCostCoefficientForAllVehicles和SetSpanCostCoefficientForVehicle。

• 通过设定边上/车辆上的代价改变目标函数

  • RoutingModel.SetArcCostEvaluatorOfAllVehicles(evaluator_index: int) - 为所有车辆设置访问arc时的代价。注意这里传入的参数是evaluator_index而非evaluate function，因此在准备好评价函数之后，还需要先进行一步register call back，获取index才能传入这个方法。
  • RoutingModel.SetArcCostEvaluatorOfVehicle(evaluator_index: int, vehicle: int) - 为单独的某辆车设置访问arc时的代价，这个方法在有driver-dependent cost时非常有用。
  • RoutingModel.SetFixedCostOfAllVehicles(cost: int64) - 很好理解的方法，为每辆车添加相同的固定成本。它还有个变种RoutingModel.SetFixedCostOfVehicle(cost: int64, vehicle: int) ，可以单独为每辆车设定固定成本。