0、前言

TSP问题是VRP问题最早的一个雏形,也就是我们常说的邮递员问题,从邮局出发,选出合适的路线,使一个邮递员所走的总路程最短,一个最简单的约束就是每一个点都要访问到,不难看出,旅行商问题(Traveling Saleman Problem,TSP)是VRP的特例!TSP和VRP问题都为NP-hard!启发式算法是目前常用解决NP问题的算法,而Ortools工具包的求解器就是基于启发式算法开发的!

1、问题描述

从原点出发经过所有需求点并回到原点,使得途经的距离最短。单车辆,多个服务点,约束:每个点服务一次。其数学描述为:设有一个城市集合C= (c1, c2, …, cn) , 其中每对城市之间的距离d (ci, cj) ∈R+, 求一对经过C中每个城市一次的路线 (cπ1, cπ2, …, cπn) ,使得:

节约里程法Java实现 节约里程算法模型_vrp 节约算法 c++

其中π1, π2, …, πn是 (1, 2, …, n) 的一个置换。

2、解决办法

利用python语言借助ortools工具包的RoutingModel模型求解的完整代码如下:

"""Simple travelling salesman problem between cities."""

from __future__ import print_function
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp

def create_data_model():
    """Stores the data for the problem."""

    data = {}

    data['distance_matrix'] = [
        [0, 2451, 713, 1018, 1631, 1374, 2408, 213, 2571, 875, 1420, 2145, 1972],
        [2451, 0, 1745, 1524, 831, 1240, 959, 2596, 403, 1589, 1374, 357, 579],
        [713, 1745, 0, 355, 920, 803, 1737, 851, 1858, 262, 940, 1453, 1260],
        [1018, 1524, 355, 0, 700, 862, 1395, 1123, 1584, 466, 1056, 1280, 987],
        [1631, 831, 920, 700, 0, 663, 1021, 1769, 949, 796, 879, 586, 371],
        [1374, 1240, 803, 862, 663, 0, 1681, 1551, 1765, 547, 225, 887, 999],
        [2408, 959, 1737, 1395, 1021, 1681, 0, 2493, 678, 1724, 1891, 1114, 701],
        [213, 2596, 851, 1123, 1769, 1551, 2493, 0, 2699, 1038, 1605, 2300, 2099],
        [2571, 403, 1858, 1584, 949, 1765, 678, 2699, 0, 1744, 1645, 653, 600],
        [875, 1589, 262, 466, 796, 547, 1724, 1038, 1744, 0, 679, 1272, 1162],
        [1420, 1374, 940, 1056, 879, 225, 1891, 1605, 1645, 679, 0, 1017, 1200],
        [2145, 357, 1453, 1280, 586, 887, 1114, 2300, 653, 1272, 1017, 0, 504],
        [1972, 579, 1260, 987, 371, 999, 701, 2099, 600, 1162, 1200, 504, 0],
    ]  # yapf: disable
    data['num_vehicles'] = 1
    data['depot'] = 0
    return data


def print_solution(manager, routing, solution):
    """Prints solution on console."""
        # Cities
    city_names = ["New York", "Los Angeles", "Chicago", "Minneapolis", "Denver", "Dallas", "Seattle",
            "Boston", "San Francisco", "St. Louis", "Houston", "Phoenix", "Salt Lake City"]
    print('Objective: {} miles'.format(solution.ObjectiveValue()))
    index = routing.Start(0)
    plan_output = 'Route for vehicle 0:\n'
    route_distance = 0
    route_name = ''
    while not routing.IsEnd(index):
        plan_output += ' {} ->'.format(manager.IndexToNode(index))
        route_name += ' {} ->'.format(str(city_names[manager.IndexToNode(index)]))
        previous_index = index
        index = solution.Value(routing.NextVar(index))
        route_distance += routing.GetArcCostForVehicle(previous_index, index, 0)
    plan_output += ' {}\n'.format(manager.IndexToNode(index))

    plan_output += 'Route distance: {}miles\n'.format(route_distance)
    print('route_name:\n%s'%(route_name))
    print(plan_output)


def main():
    """Entry point of the program."""
    # Instantiate the data problem.
    data = create_data_model()

    # Create the routing index manager.
    manager = pywrapcp.RoutingIndexManager(len(data['distance_matrix']),
                                           data['num_vehicles'], data['depot'])

    # Create Routing Model.
    routing = pywrapcp.RoutingModel(manager)


    def distance_callback(from_index, to_index):
        """Returns the distance between the two nodes."""
        # Convert from routing variable Index to distance matrix NodeIndex.
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return data['distance_matrix'][from_node][to_node]

    transit_callback_index = routing.RegisterTransitCallback(distance_callback)

    # Define cost of each arc.
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

    # Setting first solution heuristic.
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)

    # Solve the problem.
    solution = routing.SolveWithParameters(search_parameters)

    # Print solution on console.
    if solution:
        print_solution(manager, routing, solution)


if __name__ == '__main__':
    main()

结果:

Objective: 7293 miles
route_name:
 New York -> Boston -> Chicago -> Minneapolis -> Denver -> Salt Lake City -> Seattle -> San Francisco -> Los Angeles -> Phoenix -> Houston -> Dallas -> St. Louis ->
Route for vehicle 0:
 0 -> 7 -> 2 -> 3 -> 4 -> 12 -> 6 -> 8 -> 1 -> 11 -> 10 -> 5 -> 9 -> 0
Route distance: 7293miles

代码详解

a、创建数据,在本求解器中,利用的数据类型是字典,其中的距离矩阵是一个数组,它的i, j项是位置i到位置j的距离(以英里为单位),其中数组的下标对应于以下顺序的位置:

0. New York - 1. Los Angeles - 2. Chicago - 3. Minneapolis - 4. Denver - 5. allas - 6. Seattle - 7. Boston - 8. San Francisco - 9. St. Louis - 10. Houston - 11. Phoenix - 12. Salt Lake City
data = create_data_model()

b、创建每个点的索引,其中距离矩阵位置上的数字表示其下标两个点之间的物理距离,用如下函数传输获得

manager = pywrapcp.RoutingIndexManager(len(data['distance_matrix']),
                                       data['num_vehicles'], data['depot'])
routing = pywrapcp.RoutingModel(manager)

c、创建距离计算器,计算任意两点之间的距离,相当于是返回一个两个点的索引对

def distance_callback(from_index, to_index):
    """Returns the distance between the two nodes."""
    # Convert from routing variable Index to distance matrix NodeIndex.
    from_node = manager.IndexToNode(from_index)
    to_node = manager.IndexToNode(to_index)
    return data['distance_matrix'][from_node][to_node]

transit_callback_index = routing.RegisterTransitCallback(distance_callback)

d、把每两个点之间的旅行成本参数传入进去

routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

e、设置了默认的搜索参数和寻找第一个解决方案的启发式方法,后面我把寻找第一个求解方案的算法函数类型罗列出来了,可根据需要自行选择

search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.first_solution_strategy = (
    routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)

f、这步设置,相当于是告诉整个优化器,开始求解的命令

solution = routing.SolveWithParameters(search_parameters)

设置初始可行解的求解算法:

search_parameters.first_solution_strategy = (routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)

初始可行解求解算法的列表:

Function | Descriptions

AUTOMATIC | 让求解器自动选择需要使用的函数
PATH_CHEAPEST_ARC | 从路由“开始”节点开始,将其连接到费用最便宜的路由节点,然后通过迭代一直添加完所有的路由节点
PATH_MOST_CONSTRAINED_ARC | 类似于PATH_CHEAPEST_ARC,但弧是用基于比较的选择器来评估的,它会优先选择约束最大的弧。要为路由模型分配选择器,请使用ArcIsMoreConstrainedThanArc()方法。
EVALUATOR_STRATEGY | 与PATH_CHEAPEST_ARC类似,只是使用传递给SetFirstSolutionEvaluator()的函数来计算ARC成本。
SAVINGS | 储蓄算法 (Clarke & Wright). 参考: Clarke, G. & Wright, J.W.: “Scheduling of Vehicles from a Central Depot to a Number of Delivery Points”, Operations Research, Vol. 12, 1964, pp. 568-581.
SWEEP | 扫描算法 (Wren & Holliday). 参考: Anthony Wren & Alan Holliday: Computer Scheduling of Vehicles from One or More Depots to a Number of Delivery Points Operational Research Quarterly (1970-1977), Vol. 23, No. 3 (Sep., 1972), pp. 333-344.
CHRISTOFIDES | 克里斯托菲德斯算法(实际上是克里斯托菲德斯算法的一种变体,使用最大匹配代替最大匹配,这不能保证度量旅行商的近似系数为3/2)。通过扩展路径直到没有节点可以插入到路径上,对通用车辆路径模型进行处理。参考文献: Nicos Christofides, Worst-case analysis of a new heuristic for the travelling salesman problem, Report 388, Graduate School of Industrial Administration, CMU, 1976.
ALL_UNPERFORMED | 所有节点不活动,仅当节点是可选(是具有有限惩罚成本的析取约束的元素)时才找到解决方案。
BEST_INSERTION | 通过在最便宜的位置插入最便宜的节点来迭代地构建解决方案;插入成本基于路由模型的全局成本函数。截至2012年2月,仅适用于带有可选节点的模型(惩罚成本有限)。
PARALLEL_CHEAPEST_INSERTION | 通过在最便宜的位置插入最便宜的节点来迭代地构建解决方案;插入成本基于arc成本函数。比BEST_INSERTION速度快
LOCAL_CHEAPEST_INSERTION | 通过在最便宜的位置插入最便宜的节点来迭代地构建解决方案;插入成本基于arc成本函数。不同于PARALLEL_CHEAPEST_INSERTION
GLOBAL_CHEAPEST_ARC | 迭代连接产生最便宜的路段的两个节点。
LOCAL_CHEAPEST_ARC | 选择具有未绑定后继的第一个节点,并将其连接到产生最便宜的路由段的节点。
FIRST_UNBOUND_MIN_VALUE | 选择具有未绑定后继节点的第一个节点,并将其连接到第一个可用节点。这相当于选择CHOOSE_FIRST_UNBOUND策略与赋值ASSIGN_MIN_VALUE相结合。

返回值有以下情况:

Value | Descriptions
0 | ROUTING_NOT_SOLVED: Problem not solved yet.
1 | ROUTING_SUCCESS: Problem solved successfully.
2 | ROUTING_FAIL: No solution found to the problem.
3 | ROUTING_FAIL_TIMEOUT: Time limit reached before finding a solution.
4 | ROUTING_INVALID: Model, model parameters, or flags are not valid.

3、实例

为了更好的使用,看一个官网实例:自动钻孔机在电路板上钻孔。问题是要找到最短的路线,钻取板上所需的所有孔。

"""Simple travelling salesman problem on a circuit board."""

from __future__ import print_function
import math
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp


def create_data_model():
    """Stores the data for the problem."""
    data = {}
    # Locations in block units
    data['locations'] = [
        (288, 149), (288, 129), (270, 133), (256, 141), (256, 157), (246, 157),
        (236, 169), (228, 169), (228, 161), (220, 169), (212, 169), (204, 169),
        (196, 169), (188, 169), (196, 161), (188, 145), (172, 145), (164, 145),
        (156, 145), (148, 145), (140, 145), (148, 169), (164, 169), (172, 169),
        (156, 169), (140, 169), (132, 169), (124, 169), (116, 161), (104, 153),
        (104, 161), (104, 169), (90, 165), (80, 157), (64, 157), (64, 165),
        (56, 169), (56, 161), (56, 153), (56, 145), (56, 137), (56, 129),
        (56, 121), (40, 121), (40, 129), (40, 137), (40, 145), (40, 153),
        (40, 161), (40, 169), (32, 169), (32, 161), (32, 153), (32, 145),
        (32, 137), (32, 129), (32, 121), (32, 113), (40, 113), (56, 113),
        (56, 105), (48, 99), (40, 99), (32, 97), (32, 89), (24, 89),
        (16, 97), (16, 109), (8, 109), (8, 97), (8, 89), (8, 81),
        (8, 73), (8, 65), (8, 57), (16, 57), (8, 49), (8, 41),
        (24, 45), (32, 41), (32, 49), (32, 57), (32, 65), (32, 73),
        (32, 81), (40, 83), (40, 73), (40, 63), (40, 51), (44, 43),
        (44, 35), (44, 27), (32, 25), (24, 25), (16, 25), (16, 17),
        (24, 17), (32, 17), (44, 11), (56, 9), (56, 17), (56, 25),
        (56, 33), (56, 41), (64, 41), (72, 41), (72, 49), (56, 49),
        (48, 51), (56, 57), (56, 65), (48, 63), (48, 73), (56, 73),
        (56, 81), (48, 83), (56, 89), (56, 97), (104, 97), (104, 105),
        (104, 113), (104, 121), (104, 129), (104, 137), (104, 145), (116, 145),
        (124, 145), (132, 145), (132, 137), (140, 137), (148, 137), (156, 137),
        (164, 137), (172, 125), (172, 117), (172, 109), (172, 101), (172, 93),
        (172, 85), (180, 85), (180, 77), (180, 69), (180, 61), (180, 53),
        (172, 53), (172, 61), (172, 69), (172, 77), (164, 81), (148, 85),
        (124, 85), (124, 93), (124, 109), (124, 125), (124, 117), (124, 101),
        (104, 89), (104, 81), (104, 73), (104, 65), (104, 49), (104, 41),
        (104, 33), (104, 25), (104, 17), (92, 9), (80, 9), (72, 9),
        (64, 21), (72, 25), (80, 25), (80, 25), (80, 41), (88, 49),
        (104, 57), (124, 69), (124, 77), (132, 81), (140, 65), (132, 61),
        (124, 61), (124, 53), (124, 45), (124, 37), (124, 29), (132, 21),
        (124, 21), (120, 9), (128, 9), (136, 9), (148, 9), (162, 9),
        (156, 25), (172, 21), (180, 21), (180, 29), (172, 29), (172, 37),
        (172, 45), (180, 45), (180, 37), (188, 41), (196, 49), (204, 57),
        (212, 65), (220, 73), (228, 69), (228, 77), (236, 77), (236, 69),
        (236, 61), (228, 61), (228, 53), (236, 53), (236, 45), (228, 45),
        (228, 37), (236, 37), (236, 29), (228, 29), (228, 21), (236, 21),
        (252, 21), (260, 29), (260, 37), (260, 45), (260, 53), (260, 61),
        (260, 69), (260, 77), (276, 77), (276, 69), (276, 61), (276, 53),
        (284, 53), (284, 61), (284, 69), (284, 77), (284, 85), (284, 93),
        (284, 101), (288, 109), (280, 109), (276, 101), (276, 93), (276, 85),
        (268, 97), (260, 109), (252, 101), (260, 93), (260, 85), (236, 85),
        (228, 85), (228, 93), (236, 93), (236, 101), (228, 101), (228, 109),
        (228, 117), (228, 125), (220, 125), (212, 117), (204, 109), (196, 101),
        (188, 93), (180, 93), (180, 101), (180, 109), (180, 117), (180, 125),
        (196, 145), (204, 145), (212, 145), (220, 145), (228, 145), (236, 145),
        (246, 141), (252, 125), (260, 129), (280, 133)
    ]  # yapf: disable
    data['num_vehicles'] = 1
    data['depot'] = 0
    return data

#下面的函数计算数据中任意两点之间的欧式距离,并将其存储在一个数组中。因为路由求解器处理整数,所以函数将计算的距离四舍五入为整数。在本例中,舍入不会影响解决方案,但在其他情况下可能会影响。
def compute_euclidean_distance_matrix(locations):
    """Creates callback to return distance between points."""
    distances = {}
    for from_counter, from_node in enumerate(locations):
        distances[from_counter] = {}
        for to_counter, to_node in enumerate(locations):
            if from_counter == to_counter:
                distances[from_counter][to_counter] = 0
            else:
                # Euclidean distance
                distances[from_counter][to_counter] = (int(
                    math.hypot((from_node[0] - to_node[0]),
                               (from_node[1] - to_node[1]))))
    return distances


def print_solution(manager, routing, solution):
    """Prints solution on console."""
    print('Objective: {}'.format(solution.ObjectiveValue()))
    index = routing.Start(0)
    plan_output = 'Route:\n'
    route_distance = 0
    while not routing.IsEnd(index):
        plan_output += ' {} ->'.format(manager.IndexToNode(index))
        previous_index = index
        index = solution.Value(routing.NextVar(index))
        route_distance += routing.GetArcCostForVehicle(previous_index, index, 0)
    plan_output += ' {}\n'.format(manager.IndexToNode(index))
    print(plan_output)
    plan_output += 'Objective: {}m\n'.format(route_distance)


def main():
    """Entry point of the program."""
    # Instantiate the data problem.
    data = create_data_model()

    # Create the routing index manager.
    manager = pywrapcp.RoutingIndexManager(len(data['locations']),
                                           data['num_vehicles'], data['depot'])

    # Create Routing Model.
    routing = pywrapcp.RoutingModel(manager)

    distance_matrix = compute_euclidean_distance_matrix(data['locations'])

    def distance_callback(from_index, to_index):
        """Returns the distance between the two nodes."""
        # Convert from routing variable Index to distance matrix NodeIndex.
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return distance_matrix[from_node][to_node]

    transit_callback_index = routing.RegisterTransitCallback(distance_callback)

    # Define cost of each arc.
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

    # Setting first solution heuristic.
    search_parameters.local_search_metaheuristic = (
    routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)
    search_parameters.time_limit.seconds = 30
    search_parameters.log_search = True

    # Solve the problem.
    solution = routing.SolveWithParameters(search_parameters)

    # Print solution on console.
    if solution:
        print_solution(manager, routing, solution)


if __name__ == '__main__':
    main()

注意:路由求解器并不总是返回TSP的最优解,因为路由问题在计算上是难以处理的。为了找到更好的解决方案,您可以使用一种更高级的搜索策略,称为引导局部搜索,它使求解器能够避开局部最小值,即比所有邻近路径都短但不是全局最小值的解决方案,然后从局部最小值移开后,求解器继续进行搜索,所以设置初始解部分,有不一样的变动。

结果:

Objective: 2672
Route:
 0 -> 1 -> 279 -> 2 -> 278 -> 277 -> 247 -> 248 -> 249 -> 246 -> 244 -> 243 -> 242 -> 241 -> 240 -> 239 -> 238 -> 245 -> 250 -> 229 -> 228 -> 231 -> 230 -> 237 -> 236 -> 235 -> 234 -> 233 -> 232 -> 227 -> 226 -> 225 -> 224 -> 223 -> 222 -> 221 -> 220 -> 219 -> 218 -> 217 -> 216 -> 215 -> 214 -> 213 -> 212 -> 211 -> 210 -> 209 -> 208 -> 251 -> 254 -> 255 -> 257 -> 256 -> 253 -> 252 -> 207 -> 206 -> 205 -> 204 -> 203 -> 202 -> 143 -> 199 -> 201 -> 200 -> 195 -> 194 -> 193 -> 191 -> 190 -> 189 -> 188 -> 187 -> 163 -> 164 -> 165 -> 166 -> 167 -> 168 -> 169 -> 171 -> 170 -> 172 -> 105 -> 106 -> 104 -> 103 -> 107 -> 109 -> 110 -> 113 -> 114 -> 116 -> 117 -> 61 -> 62 -> 63 -> 65 -> 64 -> 84 -> 85 -> 115 -> 112 -> 86 -> 83 -> 82 -> 87 -> 111 -> 108 -> 89 -> 90 -> 91 -> 102 -> 101 -> 100 -> 99 -> 98 -> 97 -> 96 -> 95 -> 94 -> 93 -> 92 -> 79 -> 88 -> 81 -> 80 -> 78 -> 77 -> 76 -> 74 -> 75 -> 73 -> 72 -> 71 -> 70 -> 69 -> 66 -> 68 -> 67 -> 57 -> 56 -> 55 -> 54 -> 53 -> 52 -> 51 -> 50 -> 49 -> 48 -> 47 -> 46 -> 45 -> 44 -> 43 -> 58 -> 60 -> 59 -> 42 -> 41 -> 40 -> 39 -> 38 -> 37 -> 36 -> 35 -> 34 -> 33 -> 32 -> 31 -> 30 -> 29 -> 124 -> 123 -> 122 -> 121 -> 120 -> 119 -> 118 -> 156 -> 157 -> 158 -> 159 -> 174 -> 173 -> 162 -> 161 -> 160 -> 180 -> 175 -> 176 -> 177 -> 150 -> 151 -> 155 -> 152 -> 154 -> 153 -> 128 -> 129 -> 130 -> 19 -> 20 -> 127 -> 126 -> 125 -> 28 -> 27 -> 26 -> 25 -> 21 -> 24 -> 22 -> 23 -> 13 -> 12 -> 14 -> 11 -> 10 -> 9 -> 7 -> 6 -> 8 -> 274 -> 273 -> 272 -> 271 -> 270 -> 15 -> 16 -> 17 -> 18 -> 131 -> 132 -> 133 -> 269 -> 268 -> 134 -> 135 -> 267 -> 266 -> 136 -> 137 -> 138 -> 148 -> 149 -> 178 -> 179 -> 181 -> 182 -> 183 -> 184 -> 186 -> 185 -> 192 -> 196 -> 197 -> 198 -> 144 -> 145 -> 142 -> 141 -> 146 -> 147 -> 140 -> 139 -> 265 -> 264 -> 263 -> 262 -> 261 -> 260 -> 258 -> 259 -> 275 -> 276 -> 3 -> 5 -> 4 -> 0

4、总结

从上面的代码中就可以看到,TSP其实就是一个最最简化的VRP问题了!其中最重要学习的点是:如何设置参数,来进行优化!尤其注意数据结构,使用的是字典!涉及字典中的取值、赋值以及增改查!

其中的主函数(main函数:也就是我详解的部分),在后面问题的解决中是基本不变的!让我知道