39、TSP的几种精确求解方法（DFJ、MTZ、最短边破圈式、行生成式）和启发式方法（插入法、近邻法、2-opt、3-opt、模拟退火）

0、定义

TSP(Traveling Salesman Problem)，旅行商问题，又名旅行推销员问题、货郎担问题

假设一个旅行商要拜访n个城市，每个城市只能经过一次，且最后要回到起点，问如何走路程最短

1、整数模型

设cij为第i个点到第j个点的距离，xij表示是否从走到

1.1、DFJ整数模型

该模型较好理解，前面两个约束保证对每个点入度为1，出度也为1，第三个约束保证每个非空集合到其补集一定有边相连

def TSP_IP_DFJ(Dis, n, time_limit):model = Model("TSP_IP_DFJ")A = [[model.addVar(vtype="B", name="A[%s,%s]" % (i, j)) for j in range(n)] for i in range(n)]# 以总成本最小为目标model.setObjective(quicksum(Dis[i][j] * A[i][j] for i in range(n) for j in range(n)))#一度约束for i in range(n):model.addCons(quicksum(A[i][j] for j in range(n) if i != j) == 1)model.addCons(quicksum(A[j][i] for j in range(n) if i != j) == 1)#无子环约束SS = subset([i for i in range(n)])for S in SS:model.addCons(quicksum(A[i][j] for i in S for j in range(n) if j not in S) >= 1)# 设置求解时间model.setRealParam("limits/time", time_limit)model.optimize()print("\ngap:", model.getGap())A1 = [[round(model.getVal(A[i][j])) for j in range(n)] for i in range(n)]return A1

1.2、MTZ整数模型

两个模型区别在于无子环约束形式的不同，下面证明两者等价

证明

MTZ模型的好处在于其约束相较于DFJ明显变少(2^n到n^2)

def TSP_IP_MTZ(Dis, n, time_limit):model = Model("TSP_IP_MTZ")A = [[model.addVar(vtype="B", name="A[%s,%s]" % (i, j)) for j in range(n)] for i in range(n)]U = [model.addVar(name="U[%s]"%i) for i in range(n)]# 以总成本最小为目标model.setObjective(quicksum(Dis[i][j] * A[i][j] for i in range(n) for j in range(n)))# 一度约束for i in range(n):model.addCons(quicksum(A[i][j] for j in range(n) if i != j) == 1)model.addCons(quicksum(A[j][i] for j in range(n) if i != j) == 1)#MTZ无子环约束for i in range(1, n):for j in range(1, n):if i != j:model.addCons(U[i]-U[j]+(n-1)*A[i][j] <= n-2)# 设置求解时间model.setRealParam("limits/time", time_limit)model.optimize()print("\ngap:", model.getGap())A1 = [[round(model.getVal(A[i][j])) for j in range(n)] for i in range(n)]return A1

1.3、行生成模型

这里先解上面其松弛模型，再判定解是否满足条件，再不断添加不满足的约束，添加约束的方案一般为约束已生成的子环不会发生，下面列举一些方法

让其与补集有边相连，或让其内部连边不够，代码块里通过type区分（精确解）

1.3.1、多次调用模式

• def sub_TSP_row_generate(Dis, n, L_S, time_limit, type=2):model = Model("sub_TSP_row_generate")A = [[model.addVar(vtype="B", name="A[%s,%s]" % (i, j)) for j in range(n)] for i in range(n)]# 以总成本最小为目标model.setObjective(quicksum(Dis[i][j] * A[i][j] for i in range(n) for j in range(n)))# 一度约束for i in range(n):model.addCons(quicksum(A[i][j] for j in range(n) if i != j) == 1)model.addCons(quicksum(A[j][i] for j in range(n) if i != j) == 1)#两个方向限制不生成已有的子环if type == 1:#该集合与补集有边相连for S in L_S:model.addCons(quicksum(A[i][j] for i in S for j in range(n) if j not in S) >= 1)model.addCons(quicksum(A[j][i] for i in S for j in range(n) if j not in S) >= 1)else:#该集合内部边不够for S in L_S:model.addCons(quicksum(A[i][j] for i in S for j in S if i != j) <= len(S)-1)# 设置求解时间model.setRealParam("limits/time", time_limit)model.optimize()print("\ngap:", model.getGap())A1 = [[round(model.getVal(A[i][j])) for j in range(n)] for i in range(n)]return A1def TSP_IP_row_generate(Dis, n, time_limit):L_S = []ring_sub = Truewhile ring_sub:ring_sub = FalseA = sub_TSP_row_generate(Dis, n, L_S, time_limit)D = {i:j for i in range(n) for j in range(n) if A[i][j] == 1}S = {0}a = 0while len(S) < n:b = D[a]if b in S:L_S.append(S)ring_sub = Truebreakelse:S.add(b)a = bif not ring_sub:return A

1.3.2、freeTransform模式

• def LL_A_generate(A):LL_A = [[i] for i in range(len(A))]D = {i: i for i in range(len(A))}for i in range(len(A)):for j in range(len(A[i])):if A[i][j] == 1 and D[i] != D[j]:index_min = min(D[i], D[j])index_max = max(D[i], D[j])LL_A[index_min] = LL_A[index_min] + LL_A[index_max]del LL_A[index_max]for k in D:if D[k] > index_max:D[k] -= 1for k in LL_A[index_min]:D[k] = index_minreturn LL_Adef TSP_IP_exact(points, Dis, time_limit=1000, type=2):n = len(Dis)model = Model("TSP_IP_exact")model.hideOutput()A = [[model.addVar(ub=1, name="A[%s,%s]" % (i, j)) for j in range(n)] for i in range(n)]# 以总成本最小为目标model.setObjective(quicksum(Dis[i][j] * A[i][j] for i in range(n) for j in range(n) if i != j), "minimize")for i in range(n):for j in range(n):if i == j:model.addCons(A[i][j] == 0)else:model.addCons(A[i][j] <= 1)model.addCons(A[i][j] >= 0)for i in range(n):model.addCons(quicksum(A[i][j] for j in range(n) if i != j) == 1)model.addCons(quicksum(A[j][i] for j in range(n) if i != j) == 1)# 设置求解时间model.setRealParam("limits/time", time_limit)EPS = 1.e-6loop = 0isMIP = Falsewhile True:loop += 1print("精确求解第", loop, "次迭代")model.optimize()if not isMIP:A1 = [[1 if model.getVal(A[i][j]) > EPS else 0 for j in range(n)] for i in range(n)]else:A1 = [[round(model.getVal(A[i][j])) for j in range(n)] for i in range(n)]LL_A = LL_A_generate(A1)if loop % 10 == 0:A2 = [[A1[i][j] if i < j else 0 for j in range(n)] for i in range(n)]A_plot(points, A2, x_min, y_min, x_max, y_max)if len(LL_A) == 1:if not isMIP:print("转为MIP")model.freeTransform()for i in range(len(A)):for j in range(len(A[i])):model.chgVarType(A[i][j], "B")isMIP = Truecontinueelse:breakmodel.freeTransform()# 已测试过的不能在一起的点，避免再次得到if type == 1:for S in LL_A:model.addCons(quicksum(A[i][j] for i in S for j in range(n) if j not in S) >= 1)else:for S in LL_A:model.addCons(quicksum(A[i][j] for i in S for j in S if i != j) <= len(S) - 1)return A1

1.3.3、freeTransform和networkx模式

• def solve_tsp(V, c, type=1):def addcut(cut_edges):G = networkx.Graph()G.add_edges_from(cut_edges)Components = list(networkx.connected_components(G))if len(Components) == 1:return Falsemodel.freeTransform()if type == 1:# 该集合与补集有边相连for S in Components:model.addCons(quicksum(x[i, j] for i in S for j in V if j not in S) >= 1)print("cut: %s 与其补集相连" % (S))else:# 该集合内部边不够for S in Components:model.addCons(quicksum(x[i, j] for i in S for j in S if j != i) <= len(S) - 1)print("cut: len(%s) <= %s" % (S, len(S) - 1))return Truedef addcut2(cut_edges):G = networkx.Graph()G.add_edges_from(cut_edges)Components = list(networkx.connected_components(G))if len(Components) == 1:return Falsemodel.freeTransform()for S in Components:T = set(V) - set(S)print("S:", S)print("T:", T)model.addCons(quicksum(x[i, j] for i in S for j in T if j > i) +quicksum(x[i, j] for i in T for j in S if j > i) >= 2)print("cut: %s >= 2" % "+".join([("x[%s,%s]" % (i, j)) for i in S for j in T if j > i]))return Truemodel = Model("tsp")x = {}for i in V:for j in V:x[i, j] = model.addVar(ub=1, name="x(%s,%s)" % (i, j))for i in V:model.addCons(quicksum(x[i, j] for j in V if j != i) == 1, "OutDegree(%s)" % i)model.addCons(quicksum(x[j, i] for j in V if j != i) == 1, "InDegree(%s)" % i)model.setObjective(quicksum(c[i, j] * x[i, j] for i in V for j in V if j != i), "minimize")EPS = 1.e-6isMIP = Falsewhile True:model.optimize()edges = []for (i, j) in x:if model.getVal(x[i, j]) > EPS:edges.append((i, j))if addcut(edges) == False:if isMIP:  # integer variables, components connected: solution foundbreakmodel.freeTransform()for (i, j) in x:  # all components connected, switch to integer modelmodel.chgVarType(x[i, j], "B")isMIP = Truereturn model.getObjVal(), edgesdef TSP_IP_exact_official(points, Dis):nums = len(points)V = range(1, nums + 1)c = {(i, j): Dis[i - 1][j - 1] for i in range(1, nums + 1) for j in range(1, nums + 1)}obj, edges = solve_tsp(V, c)A = [[0 for j in range(nums)] for i in range(nums)]for edge in edges:A[edge[0] - 1][edge[1] - 1] = 1return A

1.3.4、最短的边破除该子环（启发式，效果一般）

• def sub_TSP_row_generate_heuristic(Dis, n, pairs_must, time_limit):model = Model("sub_TSP_row_generate_heuristic")A = [[model.addVar(vtype="B", name="A[%s,%s]" % (i, j)) for j in range(n)] for i in range(n)]# 以总成本最小为目标model.setObjective(quicksum(Dis[i][j] * A[i][j] for i in range(n) for j in range(n)))# 一度约束for i in range(n):model.addCons(quicksum(A[i][j] for j in range(n) if i != j) == 1)model.addCons(quicksum(A[j][i] for j in range(n) if i != j) == 1)#把要连的边连上for pair in pairs_must:i, j = pair[0], pair[1]model.addCons(A[i][j]+A[j][i] == 1)# 设置求解时间model.setRealParam("limits/time", time_limit)model.optimize()print("\ngap:", model.getGap())A1 = [[round(model.getVal(A[i][j])) for j in range(n)] for i in range(n)]return A1def TSP_IP_row_generate_heuristic(points, Dis, n, time_limit):pairs = pair_points_generate(Dis)pairs_must = set()ring_sub = Trueused = {i:0 for i in range(n)}while ring_sub:ring_sub = FalseA = sub_TSP_row_generate_heuristic(Dis, n, pairs_must, time_limit)D = {i:j for i in range(n) for j in range(n) if A[i][j] == 1}S = {0}a = 0while len(S) < n:b = D[a]#若生成子环，则将最短破圈边加入if b in S:ring_sub = Truefor pair in pairs:i, j = pair[0], pair[1]if ((i in S and j not in S) or (i not in S and j in S)) and used[i] <= 1 and used[j] <= 1:pairs_must.add((i, j))used[i] += 1used[j] += 1breakbreakelse:S.add(b)a = bA = deal_overlap(points, A)A = Undirected_to_directed(A)return A

2、近似算法

2.1、插入法

该方法很直观，从一个点开始，不断的把最近的邻居加进来

#近邻插入
def TSP_nearest(points, Dis):n = len(Dis)S = [0 for i in range(len(Dis))]T = [1 for i in range(len(Dis))]# 加入第一个点route = [0]S[0] = 1T[0] = 0# 开始加点while sum(T) > 0:u = route[-1]nexts = [i for i in range(n) if S[i] == 0]v = nexts[0]for j in range(1, len(nexts)):if Dis[u][v] > Dis[u][nexts[j]]:v = nexts[j]route.append(v)S[v] = 1T[v] = 0A = [[0 for j in range(n)] for i in range(n)]for i in range(len(route) - 1):A[route[i]][route[i + 1]] = 1A[route[i + 1]][route[i]] = 1A[route[0]][route[-1]] = 1A[route[-1]][route[0]] = 1A = deal_overlap(points, A)return A

2.2、凸包插入法

先把所有点的凸包算出来，然后把剩余点以损失最小的方式插入路径中

#先着凸包，再插入
def TSP_insert(points, Dis):# 每个坐标对应的索引D = {}for p in points:D[p[0]] = p[1:]ps = list(D.values())b = boundary_find(ps)route = []for p in b:idex = list(D.keys())[list(D.values()).index(p)]if idex not in route:route.append(idex)S = [i for i in range(len(points)) if i not in route]while len(S) > 0:costs = [[Dis[route[-1]][S[i]] + Dis[S[i]][route[0]] - Dis[route[-1]][route[0]]] + [Dis[route[j]][S[i]] + Dis[S[i]][route[j + 1]] - Dis[route[j]][route[j + 1]] for j in range(len(route) - 1)]for i in range(len(S))]new, idex = 0, 0for i in range(len(S)):for j in range(len(route)):if costs[i][j] < costs[new][idex]:new, idex = i, jroute.insert(idex, S[new])del S[new]n = len(points)A = [[0 for j in range(n)] for i in range(n)]for i in range(len(route) - 1):A[route[i]][route[i + 1]] = 1A[route[i + 1]][route[i]] = 1A[route[0]][route[-1]] = 1A[route[-1]][route[0]] = 1A = deal_overlap(points, A)return A

2.3、2-opt

随机初始化，不断的交换两个点的位置，更短则更新

#2-opt
def TSP_2_opt(points, Dis):route = [i for i in range(len(points))]n = len(points)flag = 1while flag == 1:flag = 0for i in range(n - 1):for j in range(i + 1, n):route1 = copy.deepcopy(route)route1[i], route1[j] = route1[j], route1[i]if sub_k_opt(route, Dis) > sub_k_opt(route1, Dis):route = route1flag = 1breakif flag == 1:breakA = [[0 for j in range(n)] for i in range(n)]for i in range(len(route) - 1):A[route[i]][route[i + 1]] = 1A[route[i + 1]][route[i]] = 1A[route[0]][route[-1]] = 1A[route[-1]][route[0]] = 1A = deal_overlap(points, A)return A

2.4、3-opt

在2-opt的基础上，多加入一个换的点

def TSP_3_opt(points, Dis):route = [i for i in range(len(points))]n = len(points)flag = 1while flag == 1:flag = 0for i in range(n - 2):for j in range(i + 1, n - 1):for k in range(j + 1, n):route1, route2, route3, route4, route5 = copy.deepcopy(route), copy.deepcopy(route), copy.deepcopy(route), copy.deepcopy(route), copy.deepcopy(route)route1[i], route1[j], route1[k] = route1[i], route1[k], route1[j]route2[i], route2[j], route2[k] = route2[j], route2[i], route2[k]route3[i], route3[j], route3[k] = route3[j], route3[k], route3[i]route4[i], route4[j], route4[k] = route4[k], route4[i], route4[j]route5[i], route5[j], route5[k] = route5[k], route5[j], route5[i]if sub_k_opt(route, Dis) > sub_k_opt(route1, Dis):route = route1flag = 1if sub_k_opt(route, Dis) > sub_k_opt(route2, Dis):route = route2flag = 1if sub_k_opt(route, Dis) > sub_k_opt(route3, Dis):route = route3flag = 1if sub_k_opt(route, Dis) > sub_k_opt(route4, Dis):route = route4flag = 1if sub_k_opt(route, Dis) > sub_k_opt(route5, Dis):route = route5flag = 1if flag == 1:breakif flag == 1:breakif flag == 1:breakA = [[0 for j in range(n)] for i in range(n)]for i in range(len(route) - 1):A[route[i]][route[i + 1]] = 1A[route[i + 1]][route[i]] = 1A[route[0]][route[-1]] = 1A[route[-1]][route[0]] = 1A = deal_overlap(points, A)return A

2.5、SA

引入模拟退火思想，一定程度接受质量不高的解

def TSP_SA(points, Dis):route = [i for i in range(len(points))]n = len(points)T = 1000 * sub_k_opt(route, Dis)T0 = 1alpha = 0.99flag = 0limit = n * math.log(n)loop = 0while flag < limit or T <= T0:loop += 1if loop % 10 == 0:print("模拟退火进行至:", loop)i = int(rd.random() * (n - 1))j = int(rd.random() * (n - i - 1)) + i + 1route1 = copy.deepcopy(route)route1[i], route1[j] = route1[j], route1[i]if sub_k_opt(route, Dis) > sub_k_opt(route1, Dis):route = route1T = T * alphaflag = 0else:dc = sub_k_opt(route1, Dis) - sub_k_opt(route, Dis)if math.exp(-dc / T) > rd.random():route = route1T = T * alphaflag = 0else:flag += 1A = [[0 for j in range(n)] for i in range(n)]for i in range(len(route) - 1):A[route[i]][route[i + 1]] = 1A[route[i + 1]][route[i]] = 1A[route[0]][route[-1]] = 1A[route[-1]][route[0]] = 1A = deal_overlap(points, A)return A

3、测试结果

3.1、各方法适用范围

3.2、测试结果

这里采用公共数据集中有最优解的case测试，总体来说凸包插入法在求解效率(秒级别出结果)和求解精度(gap4.43%)上表现较好

对比图