旅行商问题：寻找一条总距离最短的巡回路线

旅行商问题：寻找一条总距离最短的巡回路线#

本章的内容#

本章要找出这样一条路线：把所有的城市各访问一次，回到出发点，并且总距离最短。

  • 用 add_circuit() 生成一条巡回路线

  • 用 model.minimize() 使距离之和最小

  • 用回调记录目标函数的改进过程

  • 追加约束，求解“必须经过某条线段”之类变形后的 TSP

这四件事覆盖了 TSP 的基本写法与两种常见的变形。

旅行商问题#

给定若干城市，把每个城市各访问一次后回到出发点，求总移动距离最小的访问顺序，这样的问题称为旅行商问题，简称 TSP。

TSP 属于 NP 困难问题，随着城市数 \(N\) 的增加，路径的数目按 \((N-1)!/2\) 的速度增长。

城市数

路径数(约)

10

18 万

15

4 亿

20

1.2 万亿

这种规模用穷举法已经无法应付，只能交给求解器，而在 CP-SAT 中，只要把“每个城市出去 1 条、进来 1 条”的回路约束与距离的最小化组合起来即可。回路约束负责各访问一次，目标函数负责最短。

from matplotlib import pyplot as plt
import numpy as np
from scipy import spatial
from matplotlib.collections import LineCollection
from ortools.sat.python import cp_model

测试数据#

随机生成 100 个点，并算出所有点对之间的距离。由于 CP-SAT 只处理整数，因此要把距离乘以 1000 之后再取整。取整会带来误差，但 1000 倍的精度对本题已经足够。

n = 100
np.random.seed(42)
points = np.random.uniform(0, 100, size=(n, 2))

# 距离矩阵（乘 1000 变成整数）
dist_matrix = (spatial.distance.squareform(spatial.distance.pdist(points)) * 1000).astype(np.int32)

fig, ax = plt.subplots()
ax.scatter(points[:, 0], points[:, 1]);
_images/509d12f031dd3c68589d18f64a0540b01b0b1fb8471dec4c5600553dfd55b585.png

记录改进的过程#

继承 CpSolverSolutionCallback，就可以写出一个在每次找到解时被调用的函数，这里用它记录目标函数即总距离的变化，因此能看到逐次变好的过程。

class ObjectiveTracker(cp_model.CpSolverSolutionCallback):
    """记录搜索过程中找到的目标函数值"""

    def __init__(self):
        super().__init__()
        self.history = []

    def on_solution_callback(self):
        self.history.append(self.objective_value)

求解#

solve_points_tsp() 的流程如下。

  1. 由距离矩阵建立 (i, j) → 距离 的字典 edges，其中 i == j 的情况除外

  2. 为每条边建立布尔变量，为真表示经过这条边

  3. 用 add_circuit() 强制形成一条闭合回路

  4. 用 LinearExpr.weighted_sum() 计算总距离，再用 model.minimize() 使其最小

  5. 从被选中的边取出路径的坐标

第 3 步与第 4 步合起来，就是 TSP 的完整模型：结构由回路固定，优劣由目标函数排序。

只要超过 max_time 搜索就停止，这时得到的解是可行的，但未必最优，从状态是 OPTIMAL 还是 FEASIBLE 即可判断属于哪一种。

def solve_points_tsp(points, dist_matrix, log_history=False, max_time=10):
    """求绕遍点集的最短巡回路线"""
    edges = {}
    for (i, j), d in np.ndenumerate(dist_matrix):
        if i != j:
            edges[i, j] = int(d)

    model = cp_model.CpModel()
    variables = {key: model.new_bool_var(f"e{i}") for i, key in enumerate(edges)}
    model.add_circuit([(s, e, v) for (s, e), v in variables.items()])

    # 让经过的边的距离合计最小
    length = cp_model.LinearExpr.weighted_sum(list(variables.values()), list(edges.values()))
    model.minimize(length)

    solver = cp_model.CpSolver()
    solver.parameters.max_time_in_seconds = max_time
    if log_history:
        callback = ObjectiveTracker()
        status = solver.solve(model, callback)
        print("status:", solver.status_name(status))
        print("距离的变化:", callback.history)
    else:
        status = solver.solve(model)
        print("status:", solver.status_name(status))

    selected = points[np.array([key for key, v in variables.items() if solver.value(v)])]
    return selected

下面用 100 个点求解，并把选中的边连成线，看路线是否闭合成圈。这个规模属于中等，既能看出解的质量，也不至于等太久。

selected_edges = solve_points_tsp(points, dist_matrix, log_history=True)
fig, ax = plt.subplots()
ax.scatter(points[:, 0], points[:, 1])
ax.add_collection(LineCollection(selected_edges));

status: OPTIMAL
距离的变化: [5094854.0, 1132119.0, 779053.0, 777936.0, 775700.0, 774583.0, 772782.0, 755075.0, 754851.0, 746285.0, 745684.0]
_images/2491c9fde3465e6db7b76f64883ea16c6e11a4c39f6278a1a8fe0f2940263c1a.png

线段上的 TSP#

这里考虑的不是巡游“点”，而是巡游“线段”的问题。共有 50 条线段，要追加的约束是从每条线段的某一端进入，再从另一端出来。两个端点之间的移动才算用上了这条线段。

首先要把线段造出来，具体的做法见下面的代码。

n = 50
np.random.seed(42)
start_points = np.random.uniform(0, 100, size=(n, 2))
length = np.random.uniform(5, 10, size=n)
angle = np.random.uniform(0, 2 * np.pi, size=n)
end_points = start_points + np.c_[length * np.cos(angle), length * np.sin(angle)]
segments = np.column_stack([start_points, end_points])
points = segments.reshape(-1, 2)

fig, ax = plt.subplots()
ax.scatter(start_points[:, 0], start_points[:, 1])
ax.scatter(end_points[:, 0], end_points[:, 1])
ax.add_collection(LineCollection(segments.reshape(-1, 2, 2)))

<matplotlib.collections.LineCollection at 0x1e9d465ccd0>
_images/683b1b657fe2a87859c236faf2ef7a8c3225dbf292a8cda3bde040f1b31accd6.png

约束偏弱时，路线会绕开线段#

先只做一件事，即把同一条线段两端之间的距离设为 0。这样在线段上移动的代价为 0，但即使如此，也不能保证所有的线段都被走到。原因在于不走代价为 0 的边、改走捷径时，总距离反而更短。代价为 0 只是不额外收费，并不等于必须走。

下图中红线是原来的线段，黑线是被选中的路线，对照即可看出哪些线段没有被走到。

dist_matrix = (spatial.distance.squareform(spatial.distance.pdist(segments.reshape(-1, 2))) * 1000).astype(np.int32)
start_index = np.arange(0, n * 2, 2)
end_index = start_index + 1
# 把线段两端的距离设为 0
dist_matrix[start_index, end_index] = 0
dist_matrix[end_index, start_index] = 0

selected_edges = solve_points_tsp(points, dist_matrix, log_history=True)
fig, ax = plt.subplots()
ax.scatter(start_points[:, 0], start_points[:, 1])
ax.scatter(end_points[:, 0], end_points[:, 1])
ax.add_collection(LineCollection(segments.reshape(-1, 2, 2), linewidth=5, color='red', alpha=0.6))
ax.add_collection(LineCollection(selected_edges, linewidth=2))

status: OPTIMAL
距离的变化: [5215210.0, 5208664.0, 947648.0, 536334.0, 461284.0, 429924.0, 428821.0, 427685.0, 425074.0, 424315.0, 423622.0]

<matplotlib.collections.LineCollection at 0x1e9d447afd0>
_images/8ac679dfacfc330613a33242820298383ce52f1bb27c0396ded275c9eed102ff.png

用约束表达“必须经过”#

对线段 (s, e)，用 add_exactly_one() 约束 s → e 与 e → s 这两条边中必有一条被经过，这样得到的路线就会经过所有的线段。两条边中恰好一条为真，方向由此也被固定。

def solve_segments_tsp(points, dist_matrix, log_history=False, max_time=10):
    """求经过所有线段的最短路线"""
    edges = {}
    for (i, j), d in np.ndenumerate(dist_matrix):
        if i != j:
            edges[i, j] = int(d)

    n = dist_matrix.shape[0]
    start_index = np.arange(0, n, 2)
    end_index = start_index + 1

    model = cp_model.CpModel()
    variables = {key: model.new_bool_var(f"e{i}") for i, key in enumerate(edges)}
    model.add_circuit([(s, e, v) for (s, e), v in variables.items()])

    # 每条线段必须经过其中一个方向
    for s, e in zip(start_index, end_index):
        model.add_exactly_one([variables[s, e], variables[e, s]])

    length = cp_model.LinearExpr.weighted_sum(list(variables.values()), list(edges.values()))
    model.minimize(length)

    solver = cp_model.CpSolver()
    solver.parameters.max_time_in_seconds = max_time
    if log_history:
        callback = ObjectiveTracker()
        status = solver.solve(model, callback)
        print("status:", solver.status_name(status))
        print("距离的变化:", callback.history)
    else:
        status = solver.solve(model)
        print("status:", solver.status_name(status))

    selected = points[np.array([key for key, v in variables.items() if solver.value(v)])]
    return selected

这样一来，就得到了经过全部线段的路线，与上一节相比只多了 add_exactly_one 这一类约束，模型的结构没有变。

selected_edges = solve_segments_tsp(points, dist_matrix, log_history=True)
fig, ax = plt.subplots()
ax.scatter(start_points[:, 0], start_points[:, 1])
ax.scatter(end_points[:, 0], end_points[:, 1])
ax.add_collection(LineCollection(segments.reshape(-1, 2, 2), linewidth=5, color='red', alpha=0.6))
ax.add_collection(LineCollection(selected_edges, linewidth=2))

status: OPTIMAL
距离的变化: [1358864.0, 1129401.0, 823853.0, 711604.0, 630001.0, 556269.0, 544769.0, 540940.0, 538884.0, 535907.0, 532078.0, 524710.0, 498004.0, 493894.0, 489224.0, 456436.0, 455037.0, 438037.0, 437024.0, 436133.0, 434591.0, 429840.0, 429027.0, 425941.0]

<matplotlib.collections.LineCollection at 0x1e9d5755bd0>
_images/8653fe374728115632faf30c3cd86edcb95a3f5d91d7171762ef28fd9b3e16f1.png

小结#

  • TSP 可以写成“add_circuit() 加 minimize(距离之和)”

  • 距离要取整数，例如乘以 1000。CP-SAT 不能处理小数

  • 想观察解的改进过程时，继承 CpSolverSolutionCallback

  • 用 max_time_in_seconds 可以限定搜索时间。超时得到的解是 FEASIBLE，未必最优

  • “必须经过”这类条件，只要在边的变量上追加 add_exactly_one 等约束即可表达

  • 约束偏弱时，得到的路线会与预期不符，线段的例子即是一例，条件要明确写出来