CVRP問題

有容量限制的車輛路徑規(guī)劃問題(Capacitated Vehicle Routing Problem)是車輛路徑規(guī)劃問題的一類經(jīng)典變體罢绽。在這類問題中峰鄙，每個(gè)節(jié)點(diǎn)都有一個(gè)需求量瓷叫，每輛車都有最大的負(fù)載唤反，要求分配給每輛車的路徑上惫确，所有節(jié)點(diǎn)需求量之和都不超過車輛的最大負(fù)載。

這里我們用CVRP問題以及一個(gè)小變種為例眼耀，看一下在Ortools中英支，如何對目標(biāo)函數(shù)進(jìn)行自定義。

算例我們?nèi)匀徊捎肨SP問題中的gr120.tsp畔塔，車輛數(shù)設(shè)置為4輛潭辈，并且隨機(jī)生成每個(gè)節(jié)點(diǎn)的需求以及每輛車的起點(diǎn)和終點(diǎn)鸯屿。

數(shù)據(jù)準(zhǔn)備部分的代碼如下：

## 導(dǎo)入庫函數(shù)
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

## 數(shù)據(jù)準(zhǔn)備
# 給定數(shù)據(jù)目錄并讀入數(shù)據(jù)
fname: str = "gr120.tsp"
input_fpath: str = r"data/" + fname
data = tsplib95.load(input_fpath).as_name_dict()
# 讓節(jié)點(diǎn)編號從0開始
data["display_data"] = {k-1: v for k, v in data["display_data"].items()}
# 補(bǔ)充每個(gè)節(jié)點(diǎn)的demand數(shù)據(jù)
rng = np.random.default_rng(seed = 0)
data["demand"] = rng.integers(1, 11, data["dimension"]) # 隨機(jī)生成[1, 10]的需求
# 車輛數(shù)據(jù)
n_vehicles = 4 # 車輛數(shù)
capacity = [250] * n_vehicles # 車輛最大負(fù)載
vehicle_start_idx = rng.integers(0, data["dimension"], n_vehicles) # 每輛車的起點(diǎn)
# 注意ortools是c++寫的澈吨，它并不接受 np.int64，因此要轉(zhuǎn)化為python原生的int類型
vehicle_start_idx = [i.item() for i in vehicle_start_idx]
vehicle_end_idx = rng.integers(0, data["dimension"], n_vehicles) # 每輛車的終點(diǎn)
vehicle_end_idx = [i.item() for i in vehicle_end_idx]

在完成數(shù)據(jù)準(zhǔn)備之后寄摆，可以通過下面代碼對數(shù)據(jù)進(jìn)行可視化：

# 節(jié)點(diǎn)可視化
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)

如下圖谅辣，在途中，紅色星形表示車輛起點(diǎn)婶恼，藍(lán)色星形表示車輛終點(diǎn)桑阶，其他圓形代表一般節(jié)點(diǎn)：

gr120.tsp_start_end_vrp.png

模型與求解

## 建立并求解CVRP問題
# 建立RoutingIndexManager
manager = pywrapcp.RoutingIndexManager(data["dimension"], n_vehicles, 
                             vehicle_start_idx, vehicle_end_idx)
# 實(shí)例化RoutingModel
routing = pywrapcp.RoutingModel(manager)
# 建立一個(gè)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)
# 給每條邊設(shè)置cost
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)
# 建立一個(gè)callback來記錄車輛的負(fù)載
def demand_callback(from_index: int):
    # 經(jīng)過一個(gè)節(jié)點(diǎn)之后，累加節(jié)點(diǎn)上的需求量
    from_node = manager.IndexToNode(from_index)
    return data['demand'][from_node]
demand_callback_index = routing.RegisterUnaryTransitCallback(demand_callback)
# 利用這個(gè)callback勾邦，施加capacity約束
routing.AddDimensionWithVehicleCapacity(
        demand_callback_index,
        0,  # 不允許松弛
        capacity,  # 車輛的最大允許負(fù)載
        True,  # 從0開始進(jìn)行累積
        'Capacity')
# 設(shè)定求解參數(shù)
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) # 設(shè)置最大求解時(shí)間100秒

# 求解
solution = routing.SolveWithParameters(search_parameters)

# 輸出結(jié)果
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 += '路徑總負(fù)載: {}\n'.format(route_load)
    print(plan_output)
    paths.append(path)
    total_distance += route_distance
    total_load += route_load
print('所有路徑總長度: {}m'.format(total_distance))
print('所有路徑總負(fù)載: {}'.format(total_load))

得到解如下：

Objective: 33045
車輛 0 的路徑:
 68 裝載(9) ->  113 裝載(9)
路徑總長度: 0m
路徑總負(fù)載: 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
路徑總負(fù)載: 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
路徑總負(fù)載: 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
路徑總負(fù)載: 192

所有路徑總長度: 33045m
所有路徑總負(fù)載: 644

可以看到蚣录，這里我們的路徑總長度和目標(biāo)函數(shù)值是相等的。

用如下代碼進(jìn)行路徑的可視化：

# 路徑可視化
colors = ['r', 'g', 'b', 'k']
# 畫出所有節(jié)點(diǎn)
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)
# 畫出每輛車的起點(diǎn)和終點(diǎn)
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

車輛負(fù)載的平衡

在上面的求解例子中眷篇，可以看到車輛之間的負(fù)載是恨不均衡的萎河， 0號車只經(jīng)過了起點(diǎn)和終點(diǎn)，總共負(fù)載為9，而車輛1的負(fù)載則達(dá)到了250虐杯。在實(shí)際使用時(shí)玛歌，我們通常希望我們的車輛有相近的負(fù)載，不要讓某些員工過量工作而某些員工整天閑逛擎椰。因此我們需要對目標(biāo)函數(shù)進(jìn)行修改支子，讓負(fù)載的均衡也成為我們目標(biāo)的一部分。

我們可以這樣處理負(fù)載均衡的要求：

1. 通過累加所有節(jié)點(diǎn)上的需求量达舒，得到每輛車平均需要承擔(dān)的負(fù)載（這里更精細(xì)一點(diǎn)的處理應(yīng)該去除掉depot的需求量值朋，因?yàn)橥ǔ６x上我們不需要服務(wù)車庫）；
2. 統(tǒng)計(jì)每輛車路徑上的負(fù)載巩搏，當(dāng)車輛負(fù)載超過平均負(fù)載時(shí)吞歼，施加一個(gè)懲罰函數(shù)

具體實(shí)現(xiàn)如下，demand_dimension上累加了車輛在路徑上的負(fù)載塔猾，我們在車輛到達(dá)最后一個(gè)節(jié)點(diǎn)時(shí)篙骡，檢查車輛負(fù)載是不是超過我們設(shè)定的平均值，如果出現(xiàn)超載丈甸，則通過penalty_coef施加一個(gè)懲罰項(xiàng)：

# 通過軟約束來平衡車輛的負(fù)載
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)

得到的結(jié)果如下：

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
路徑總負(fù)載: 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
路徑總負(fù)載: 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
路徑總負(fù)載: 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
路徑總負(fù)載: 150

所有路徑總長度: 37059m
所有路徑總負(fù)載: 644

可以看到糯俗，我們的平均負(fù)載是 670/4=167（輸出的644是減去了4個(gè)車庫之后的量，和我們代碼中的處理稍有不同）睦擂。如果我們施加一個(gè)較小的平均負(fù)載值得湘，我們就可以看到路徑長度和目標(biāo)函數(shù)之間的差異，差異的部分也就是懲罰項(xiàng)顿仇。

對路徑進(jìn)行可視化如下：

vrp_paths_with_load_balance.png

改變目標(biāo)函數(shù)的方式

從上面這個(gè)負(fù)載均衡問題引申出來的問題就是：我們怎么樣去在不同的需求下淘正，對目標(biāo)函數(shù)進(jìn)行定制。

Ortools的VRP求解器中臼闻，并沒有明確的定義目標(biāo)函數(shù)鸿吆。而是在通用的意義上，將最小化所有車輛在經(jīng)過的邊(arc)上的代價(jià)(cost)作為目標(biāo)函數(shù)述呐。

在實(shí)際應(yīng)用中惩淳，比較常見的改變目標(biāo)函數(shù)的方式有下面三種：

• 通過軟約束在目標(biāo)函數(shù)中添加懲罰項(xiàng)

  • RoutingDimension.SetCumulVarSoftUpperBound(index: int64, upper_bound: int64, coefficient: int64) - 為一個(gè)累積量設(shè)置一個(gè)上界，如果累積量的值超過上界乓搬，則在目標(biāo)函數(shù)中施加一個(gè)和coefficient成比例的懲罰項(xiàng)思犁，也就是說cumulVar <= upper_bound -> cost = 0，cumulVar > upper_bound -> cost = coefficient * (cumulVar - upper_bound)进肯；這個(gè)方法還有對應(yīng)的設(shè)置下界的版本激蹲。
  • RoutingModel.AddDisjunction(indices: List[int64], penalty = 0: int64, max_cardinality: int64 = 1) - 添加一個(gè)disjuction constraint，要求在傳入的indices代表的節(jié)點(diǎn)中有max_cardinality個(gè)被訪問到（注意這里在indices中江掩，我們不能傳入車輛的起點(diǎn)和終點(diǎn)）学辱，如果沒有足夠數(shù)量的節(jié)點(diǎn)被訪問含蓉，那么我們的懲罰為 p*penalty，其中p=max_cardinality - sum(active[i])這里要求傳入的penalty為正數(shù)项郊，否則將會(huì)強(qiáng)制這些節(jié)點(diǎn)不會(huì)被訪問馅扣；它的一個(gè)變種應(yīng)用是當(dāng)我們有一些節(jié)點(diǎn)可選時(shí)，我們將每個(gè)節(jié)點(diǎn)的index傳入着降，這樣當(dāng)我們drop掉這個(gè)節(jié)點(diǎn)時(shí)差油，我們會(huì)受到懲罰，但是不再強(qiáng)制每個(gè)節(jié)點(diǎn)都會(huì)被訪問任洞。

• 通過系數(shù)設(shè)置在目標(biāo)函數(shù)中添加自定義項(xiàng)

  • RoutingDimension.SetGlobalSpanCostCoefficient(coefficient: int64) - 這里所謂的Global蓄喇，指的是全局所有車輛中該累積量在終點(diǎn)的最大值和在起點(diǎn)的最小值之差，也就是global_span_cost = coefficient * (Max(dimension end value) - Min(dimension start value))交掏；在限定車輛服務(wù)的節(jié)點(diǎn)數(shù)時(shí)妆偏，我們常用這種方法。它還有兩個(gè)類似的方法SetSpanCostCoefficientForAllVehicles和SetSpanCostCoefficientForVehicle盅弛。

• 通過設(shè)定邊上/車輛上的代價(jià)改變目標(biāo)函數(shù)

  • RoutingModel.SetArcCostEvaluatorOfAllVehicles(evaluator_index: int) - 為所有車輛設(shè)置訪問arc時(shí)的代價(jià)钱骂。注意這里傳入的參數(shù)是evaluator_index而非evaluate function，因此在準(zhǔn)備好評價(jià)函數(shù)之后挪鹏，還需要先進(jìn)行一步register call back见秽，獲取index才能傳入這個(gè)方法。
  • RoutingModel.SetArcCostEvaluatorOfVehicle(evaluator_index: int, vehicle: int) - 為單獨(dú)的某輛車設(shè)置訪問arc時(shí)的代價(jià)讨盒，這個(gè)方法在有driver-dependent cost時(shí)非常有用解取。
  • RoutingModel.SetFixedCostOfAllVehicles(cost: int64) - 很好理解的方法，為每輛車添加相同的固定成本返顺。它還有個(gè)變種RoutingModel.SetFixedCostOfVehicle(cost: int64, vehicle: int) 禀苦，可以單獨(dú)為每輛車設(shè)定固定成本。