多站点的TSP问题求解-06
https://mit-realm.github.io/CMDTSP/
https://github.com/Brelliothe/CMDTSP
paper：A Hierarchical Framework for Solving the Constrained Multiple Depot Traveling Salesman Problem
最基础的 树形 DP（Tree DP）
什么是树 DP？
树 DP = 在一棵树上做动态规划。
本质：递归 + 记忆化。利用树的层次结构，把子树的答案合并起来，推到父节点。
和线性 DP 不同的是：
线性 DP 一般是「前 i 个元素的最优解」。 0-1背包问题
树 DP 是「以某个节点为根的子树的最优解」。
核心思路
树 DP 的套路一般有两步：
设计状态：
常见设计是 dp[u] = 子树 u 的最优答案
或者 dp[u][0/1] = 子树 u 在某种约束下的最优答案
转移方程：
从孩子结点的结果，推父结点。
因为树是递归结构，这一步天然适合「DFS 遍历 + 合并答案」。
简单的例子：树上最大独立集
树 DP 的「hello world」问题
题目：给你一棵树，每个节点有一个权值，要求选一些点，使得相邻两个点不能同时选，求权值和最大。
状态设计：
dp[u][0] = u 不选 时，子树 u 的最大权值和
dp[u][1] = u 选 时，子树 u 的最大权值和
状态转移：
如果 u 选了，那么它的孩子都不能选：dp[u][1] = w[u] + Σ dp[v][0] （v 是 u 的孩子）
如果 u 不选，那么孩子可以选，也可以不选，取最优：dp[u][0] = Σ max(dp[v][0], dp[v][1])
整棵树的答案就是：max(dp[root][0], dp[root][1])
代码：
from collections import defaultdict
import sys
sys.setrecursionlimit(10**6)
n = 5
w = [0, 1, 2, 3, 4, 5] # 节点权值，1-index
tree = defaultdict(list)
edges = [(1,2),(1,3),(2,4),(2,5)]
for u,v in edges:
tree[u].append(v)
tree[v].append(u)
dp = [[0,0] for _ in range(n+1)]
visited = [False]*(n+1)
def dfs(u):
visited[u] = True
dp[u][1] = w[u] # 选 u
dp[u][0] = 0 # 不选 u
for v in tree[u]:
if not visited[v]:
dfs(v)
dp[u][0] += max(dp[v][0], dp[v][1]) # u 不选
dp[u][1] += dp[v][0] # u 选
dfs(1)
print("最大独立集权值和:", max(dp[1][0], dp[1][1]))
多depots分配cities，TSP求解的过程，核心的PTC分组是如何进行的
# our method
from src.baseline import Baseline
import gurobipy as gp
from gurobipy import GRB
from itertools import islice
import networkx as nx
from networkx.algorithms import approximation
import numpy as np
from python_tsp.heuristics import solve_tsp_local_search, solve_tsp_simulated_annealing
class PTC(Baseline):
def __init__(self, depots, cities, stations, graph, distance, full, limit, method='lkh'):
super().__init__(depots, cities, stations, graph, distance, full, limit, 'Hierarchical Pipeline')
self.k = 10
self.group = []
self.lb, self.ub = 0, 1000
self.shortest_paths = {}
self.method = method
if method == 'neural':
import torch
self.model = load_model('models/gat/pretrained/tsp_50/').cuda()
self.model.set_decode_type("sampling", temp=1)
def partition(self):
"""
allocate the cities to salesman by partition the minimum spanning tree
"""
# find the minimum spanning tree
graph = nx.Graph()
for i in np.concatenate((self.depots, self.cities)):
graph.add_node(i)
for i in graph.nodes:
for j in graph.nodes:
if i == j:
continue
graph.add_edge(i, j, weight=self.distance[i][j])
tree = nx.minimum_spanning_tree(graph)
# turn the spanning tree into a rooted tree
for i in tree.nodes:
tree.nodes[i]['parent'] = -1
def rooted_tree(node):
# con represent the minimum cost in the subtree when the connected depot is its child
# ncon represent the minimum cost in the subtree when the connected depot is not its child
# for the depot node, treat itself as its child, set its ncon as -1 to be invalid
cons, ncons, index, diff = [], [], [], []
for n in tree.neighbors(node):
if n != tree.nodes[node]['parent']:
# indicator whether n will connect to a depot via its parent
tree.nodes[n]['pcon'] = True
tree.nodes[n]['parent'] = node
con, ncon = rooted_tree(n)
tree.nodes[n]['con'] = con
tree.nodes[n]['ncon'] = ncon
# n is a source
if n in self.depots:
cons.append(con + tree[node][n]['weight'])
ncons.append(con)
diff.append(cons[-1] - ncons[-1])
index.append(n)
# n is a target and has source inside the subtree
elif con != -1:
cons.append(con + tree[node][n]['weight'])
ncons.append(min(ncon + tree[node][n]['weight'], con))
# the node will not connect to its parent if not connect source to its parent
diff.append(cons[-1] - ncons[-1])
index.append(n)
if ncons[-1] == con:
tree.nodes[n]['pcon'] = False
# n do not have source within the subtree
else:
ncons.append(ncon + tree[node][n]['weight'])
if node in self.depots:
# it is automatically connected to itself
con = sum(ncons)
ncon = -1
elif len(diff) == 0:
con = -1
ncon = sum(ncons)
else:
id = np.argmin(np.array(diff))
con = min(diff) + sum(ncons)
ncon = sum(ncons)
tree.nodes[node]['child'] = index[id]
return con, ncon
con, ncon = rooted_tree(self.depots[0])
tree.nodes[self.depots[0]]['con'] = con
tree.nodes[self.depots[0]]['ncon'] = ncon
tree.nodes[self.depots[0]]['pcon'] = False
def assign_group(node, value):
if node in self.depots:
# case 1: node is a depot, all children do not connect to both it and a depot inside the subtree
# the node belongs to the group named by itself
tree.nodes[node]['group'] = np.where(self.depots == node)[0].item()
# for all neighbor nodes besides the parent
for n in tree.neighbors(node):
if n != tree.nodes[node]['parent']:
# if node n is also a depot
if n in self.depots:
# it belongs to the group named by itself
tree.nodes[n]['group'] = np.where(self.depots == n)[0].item()
assign_group(n, tree.nodes[n]['con'])
# if node n is a city but connects to 'node'
elif tree.nodes[n]['pcon']:
# assign its group to be node
tree.nodes[n]['group'] = tree.nodes[node]['group']
assign_group(n, tree.nodes[n]['ncon'])
# node n is a city and connects to some depot inside the subtree
else:
tree.nodes[n]['group'] = assign_group(n, tree.nodes[n]['con'])
# case 2: node is a city, and it connects to a depot inside the subtree
elif value == tree.nodes[node]['con']:
# find its child whose subtree contains the depot
n = tree.nodes[node]['child']
# if child n is a depot
if n in self.depots:
index = np.where(self.depots == n)[0].item()
tree.nodes[node]['group'] = index
tree.nodes[n]['group'] = index
assign_group(n, tree.nodes[n]['con'])
# if child n is a city
else:
tree.nodes[node]['group'] = assign_group(n, tree.nodes[n]['con'])
# other neighbors besides child n and parent
for n in tree.neighbors(node):
if n != tree.nodes[node]['parent'] and n != tree.nodes[node]['child']:
# depot should not connect to it
if n in self.depots:
tree.nodes[n]['group'] = np.where(self.depots == n)[0].item()
assign_group(n, tree.nodes[n]['con'])
# city connect to 'node'
elif tree.nodes[n]['pcon']:
tree.nodes[n]['group'] = tree.nodes[node]['group']
assign_group(n, tree.nodes[n]['ncon'])
# city not connect to 'node'
else:
tree.nodes[n]['group'] = assign_group(n, tree.nodes[n]['con'])
else:
# node is a city and the connected depot is outside the subtree
# a city is visited before so node has already been assigned a group
for n in tree.neighbors(node):
if n != tree.nodes[node]['parent']:
if n in self.depots:
tree.nodes[n]['group'] = np.where(self.depots == n)[0].item()
assign_group(n, tree.nodes[n]['con'])
elif tree.nodes[n]['pcon']:
tree.nodes[n]['group'] = tree.nodes[node]['group']
assign_group(n, tree.nodes[n]['ncon'])
else:
tree.nodes[n]['group'] = assign_group(n, tree.nodes[n]['con'])
return tree.nodes[node]['group']
# start from the first node
assign_group(self.depots[0], tree.nodes[self.depots[0]]['con'])
self.group = [[depot] for depot in self.depots]
for node in tree.nodes:
if node not in self.depots:
self.group[tree.nodes[node]['group']].append(node)
def christofides(self, group):
if len(group) == 1:
return [group[0], group[0]]
elif len(group) == 2:
return [group[0], group[1], group[0]]
else:
graph = nx.Graph()
for i in group:
for j in group:
graph.add_edge(i, j, weight=self.distance[i][j])
return approximation.christofides(graph)
def exact(self, group):
if len(group) == 1:
return [group[0], group[0]]
elif len(group) == 2:
return [group[0], group[1], group[0]]
size = len(group)
model = gp.Model('TSP')
model.setParam(GRB.Param.OutputFlag, 0)
edges = model.addVars(size, size, vtype=GRB.BINARY, name='edges')
flows = model.addVars(size, size, vtype=GRB.CONTINUOUS, name='flows')
# no self loop
model.addConstrs(edges[i, i] == 0 for i in range(size))
model.addConstrs(flows[i, i] == 0 for i in range(size))
# one in and one out
model.addConstrs(gp.quicksum([edges[i, j] for i in range(size)]) == 1 for j in range(size))
model.addConstrs(gp.quicksum([edges[i, j] for j in range(size)]) == 1 for i in range(size))
# initial number of items
model.addConstr(gp.quicksum([flows[0, j] for j in range(size)]) == size - 1)
# drop one item at each loc
model.addConstrs(
gp.quicksum([flows[i, j] for i in range(size)]) - gp.quicksum([flows[j, k] for k in range(size)])
== 1 for j in range(1, size))
# set the upper bound of loc
model.addConstrs(flows[i, j] <= edges[i, j] * (size - 1) for i in range(size) for j in range(size))
# minimize the total length
model.setObjective(gp.quicksum([edges[i, j] * self.distance[group[i]][group[j]] for i in range(size)
for j in range(size) if j != i]))
model.optimize()
tour = [0 for _ in group]
for key in edges.keys():
if edges[key].x > 0.5:
tour[key[0]] = key[1]
visit = [0]
for _ in group:
visit.append(tour[visit[-1]])
return [group[i] for i in visit]
def simulated_annealing(self, group):
if len(group) == 1:
return [group[0], group[0]]
elif len(group) == 2:
return [group[0], group[1], group[0]]
distance_matrix = np.array([[self.distance[j][i] for i in group] for j in group])
visit, _ = solve_tsp_simulated_annealing(distance_matrix)
visit, _ = solve_tsp_local_search(distance_matrix, x0=visit, perturbation_scheme="ps3")
return [group[i] for i in visit] + [group[0]]
def neural(self, group):
if len(group) == 1:
return [group[0], group[0]]
elif len(group) == 2:
return [group[0], group[1], group[0]]
import torch
import math
def get_best(sequences, cost, ids=None, batch_size=None):
if ids is None:
idx = cost.argmin()
return sequences[idx:idx + 1, ...], cost[idx:idx + 1, ...]
splits = np.hstack([0, np.where(ids[:-1] != ids[1:])[0] + 1])
mincosts = np.minimum.reduceat(cost, splits)
group_lengths = np.diff(np.hstack([splits, len(ids)]))
all_argmin = np.flatnonzero(np.repeat(mincosts, group_lengths) == cost)
result = np.full(len(group_lengths) if batch_size is None else batch_size, -1, dtype=int)
result[ids[all_argmin[::-1]]] = all_argmin[::-1]
return [sequences[i] if i >= 0 else None for i in result], [cost[i] if i >= 0 else math.inf for i in result]
sequences, costs = self.model.sample_many(torch.tensor([self.graph.nodes[i]['pos'] for i in group],
dtype=torch.float32).unsqueeze(0).cuda(),
batch_rep=1280, iter_rep=1)
ids = torch.arange(1, dtype=torch.int64, device=costs.device)
if sequences is None:
sequences = [None] * 1
else:
sequences, _ = get_best(
sequences.cpu().numpy(), costs.cpu().numpy(),
ids.cpu().numpy() if ids is not None else None,
1
)
return sequences[0].tolist()
def get_shortest_paths(self):
graph = nx.Graph()
for i in self.stations:
graph.add_node(i)
for i in graph.nodes:
for j in graph.nodes:
if i == j:
continue
graph.add_edge(i, j, weight=self.distance[i][j])
for i in self.stations:
for j in self.stations:
if i != j:
self.shortest_paths[i, j] = [self.k_shortest_paths(graph, i, j, self.k), []]
for path in self.shortest_paths[i, j][0]:
self.shortest_paths[i, j][-1].append(self.path_length(graph, path))
else:
self.shortest_paths[i, i] = [[[i] for _ in range(self.k)], [0 for _ in range(self.k)]]
@staticmethod
def k_shortest_paths(graph, i, j, k):
return list(islice(nx.shortest_simple_paths(graph, i, j, weight='weight'), k))
@staticmethod
def path_length(graph, path):
length = 0
for u, v in zip(path[:-1], path[1:]):
length += graph[u][v]['weight']
return length
def heuristic_search(self, visit):
paths, costs, energy = [], [], []
# for each pair of city, find the shortest path
for source, target in zip(visit[:-1], visit[1:]):
paths.append([])
costs.append([])
energy.append([])
if self.distance[source][target] <= self.full:
paths[-1].append([source, target])
costs[-1].append(self.distance[source][target])
energy[-1].append((0, 1e10, self.distance[source][target]))
path_candidates = []
# check all valid paths
for i in self.stations:
for j in self.stations:
for k in range(len(self.shortest_paths[i, j][0])):
path_candidates.append((i, j, k, self.shortest_paths[i, j][-1][k]
+ self.distance[i][source]
+ self.distance[i][target]))
# sort valid paths by the cost
path_candidates.sort(key=lambda t: t[-1])
# take the top num paths, fill empty part by the last valid path
while len(paths[-1]) < self.k:
if len(path_candidates) == 0:
paths[-1].append(paths[-1][-1])
costs[-1].append(costs[-1][-1])
energy[-1].append(energy[-1][-1])
continue
i, j, k, c = path_candidates.pop(0)
path = [source] + [t for t in self.shortest_paths[i, j][0][k]] + [target]
if path not in paths[-1]:
paths[-1].append(path)
costs[-1].append(c)
energy[-1].append((self.distance[source][path[1]], self.full - self.distance[target][path[-2]], -1e10))
return paths, costs, energy
def congestion_control(self, paths, costs, power):
solution, cost = [], 0
a, d, n, s = len(paths), max([len(x) for x in paths]), len(paths[0][0]), len(paths[0][0][0])
model = gp.Model('control')
model.setParam(GRB.Param.OutputFlag, 0)
choices = model.addVars(a, d, n, vtype=GRB.BINARY, name='choices')
energy = model.addVars(a, d + 1, vtype=GRB.CONTINUOUS, name='energy')
model.addConstrs(choices[i, j, k] == 0 for i in range(a) for j in range(len(paths[i]), d) for k in range(n))
model.addConstrs(gp.quicksum([choices[i, j, k] for k in range(n)]) == 1 for i in range(a)
for j in range(len(paths[i])))
model.addConstrs(gp.quicksum([choices[i, j, k] * paths[i][j][k][m] for i in range(a)
for j in range(len(paths[i])) for k in range(n)]) <= self.limit for m in range(s))
model.addConstrs(energy[i, j] >= 0 for i in range(a) for j in range(d + 1))
model.addConstrs(
energy[i, j] >= gp.quicksum([power[i][j][k][0] * choices[i, j, k] for k in range(n)]) for i in range(a)
for j in range(len(paths[i])))
model.addConstrs(
energy[i, j] <= gp.quicksum([power[i][j - 1][k][1] * choices[i, j - 1, k] for k in range(n)])
for i in range(a) for j in range(1, len(paths[i]) + 1))
model.addConstrs(
energy[i, j] - energy[i, j + 1] >= gp.quicksum([power[i][j][k][2] * choices[i, j, k] for k in range(n)])
for i in range(a) for j in range(len(paths[i])))
model.addConstrs(energy[i, j] <= self.full for j in range(d + 1) for i in range(a))
model.addConstrs(energy[i, j] == 0 for i in range(a) for j in range(len(paths[i]) + 1, d + 1))
model.setObjective(
gp.quicksum([choices[i, j, k] * costs[i][j][k] for i in range(a) for j in range(len(paths[i]))
for k in range(n)]))
model.optimize()
if model.Status == GRB.INFEASIBLE:
model.computeIIS()
model.write("model.ilp")
return None, 0
else:
for i in range(a):
solution.append([])
for j in range(d):
for k in range(n):
if choices[i, j, k].x > 0.5:
solution[i].append((j, k))
cost += costs[i][j][k]
return solution, cost
def solve(self):
partitioner = {}
self.partition()
self.get_shortest_paths()
paths, proposal, price, power = [], [], [], []
for group in self.group:
tsp_solver = {'lkh': self.lkh, 'appr': self.christofides, 'exact': self.exact,
'meta': self.simulated_annealing, 'neural': self.neural}[self.method]
visit = tsp_solver(group)
path, cost, energy = self.heuristic_search(visit)
paths.append(path)
proposal.append([[[1 if s in path[i][j] else 0 for s in self.stations] for j in range(self.k)] for i in
range(len(path))])
price.append(cost)
power.append(energy)
solution, cost = self.congestion_control(proposal, price, power)
if solution is not None:
solution = [[loc for index in solution[i] for loc in paths[i][index[0]][index[1]][:-1]]
for i in range(len(solution))]
self.solution = [solution[i] + [self.depots[i]] for i in range(len(solution))]
PTC（基于树的分区）分组
下面把你代码里 PTC 分组的核心思想 + 每个变量/分支的含义 + 算法流程都拆开来讲清楚，rooted_tree → assign_group 的两遍树遍历：自底向上 DP，再自上向下分配）。
把所有点（depots + waypoints）用完全图的边权（距离）建图，取最小生成树（MST）。
在这棵树上做一次自底向上的动态规划（rooted_tree），为每个节点计算两类代价：
- 若子树内有 depot（内部连通）的最小代价 con，
- 若子树不含 depot、必须向外连通（外连）的最小代价 ncon。
再做一次自上向下的遍历（assign_group），按照上一步算出的决策（哪个子树负责“把 depot 带回”）把每个 waypoint 分配给某个 depot（即某架无人机）。
关键变量含义（对应代码）
tree：MST（节点索引是 0..len(all_points)-1，其中前 len(depots) 个是 depot
depot_indices_array：depot 的索引集合（方便判断节点是不是 depot）
对每个节点 node，算法维护：
tree.nodes[node]['parent']：在树中的父节点索引（在 rooted_tree 中设置）。
tree.nodes[node]['pcon']：布尔标志，表示该子节点是否倾向于“与父节点（外部）相连”。
True：该子树更倾向由父/外部 depot 服务（会把子树视为“外连”并继承父组）。 依附外部成一组
False：该子树包含内部 depot 并倾向用自己内部的 depot 服务（会自己被分组）。自成一组
tree.nodes[node]['con']：若子树内存在 depot，把该子树“内部服务”时的最小花费（否则设 -1 作为 sentinel）。-- 自成一组时的最小cost。
tree.nodes[node]['ncon']：若**子树不内部连接（需外连）**时的最小花费（如果子树是 depot，则 ncon = -1）。-- 依附外部成一组最小cost。
tree.nodes[node]['child']：当 con 有意义时，指向 “把 depot 带回父节点/本节点所选的那个子孩子”（该孩子的子树包含被选的 depot）。-- depo作为树根当前node连接的是它的哪个child
rooted_tree 内临时用到的：
ons / ncons / index / diff：分别用于收集每个子分支在“内部连接 vs 外连”下的代价，以及 diff = cons - ncons 用于选择哪个子分支“带回” depot。
rooted_tree(node)（自底向上 DP）详解
对每个节点做后序遍历（先处理所有子节点，再处理自己）。对于当前节点 node：
-
对于每个子邻居 n（排除父节点）递归得到 con_n, ncon_n，并把 con_n / ncon_n 放回 tree.nodes[n] 以供后续使用。
-
根据 n 的类型分三类处理（并向 cons / ncons / diff / index 添加对应条目）：
n 是 depot（子树根就是 depot）：
n 是 waypoint 且其子树中有 depot（con_n != -1）：n 是 waypoint 且其子树中无 depot (con_n == -1)：
这种子必须外连，只有 ncon_n + weight 是有效贡献，加入到 ncons 中（没有 cons/diff/index 条目，因为它不能“带回 depot”）。
直观一句话：
对一个节点，ncon 表示“如果这个子树不在内部有 depot，需要花多少钱（都靠外连）”；
con 表示“如果这个子树内有 depot（或我们选择把某个子分支的 depot 带回来），最小需要的总花费是多少”。
算法在子分支间选择哪一个子分支把 depot 带回来以最小化 con
ssign_group(node, value)（自上向下分配）详解
怎么把组 id 真正打到每个节点？
rooted_tree 给了每个节点 con / ncon / child / pcon 的信息，assign_group 根据这些信息把具体的 depot 索引（即无人机索引）分配到每个节点：
- 如果当前 node 是 depot：
- 如果当前 node 不是 depot 且 value == tree.nodes[node]['con'] 即在这个节点选择了“内部连通”方案
- 如果当前 node 不是 depot 且 value != con（意味着这里按 ncon 情况——外连）
最终 tree.nodes[node]['group'] 会是某个 depot 的索引（depot_indices_array 中的索引），代表该节点所属无人机组。
最终把节点映射成“无人机的 waypoint 列表”
实现里遍历所有 tree.nodes，对非 depot 的节点（即 waypoint nodes）：
读取 group_id = tree.nodes[node]['group']（这是 depot 的序号 0..len(drones)-1），
把 node（注意 node 是相对于 all_points 的索引）加入 drone_groups[group_id]。
并把 self.group_assignments[waypoint_idx] = group_id（其中 waypoint_idx = node - len(self.depots) 用来把 all_points 索引转换回 waypoints 的索引）。
总结
MST 把问题简化到一棵树上，避免全连通图的复杂组合。
con/ncon 的 DP 是典型的“树形二选一”策略：对子树而言，要么内部有 depot（可以在子树内部解决）要么没有（必须靠外边）。
diff = cost_if_selected - cost_if_not_selected 用来衡量 把某个子树的 depot “带到父节点” 的额外开销，从而在父节点做最优选择（选哪条子边来带 depot）。
顶层从某个 depot 开始向下展开，最终每个 waypoint 被分配到一个 depot（组）。
浙公网安备 33010602011771号