求解旅行商问题的三种精确性建模方法，性能差距巨大

旅行商问题介绍

旅行商问题 (Traveling Salesman Problem, TSP) 是一个经典的组合优化问题，目标是找到一条最短路径，该路径必须经过每个城市一次且仅一次，最终返回到起始城市。

问题的输入：

• $N$: 城市的总数。
• $c$: 城市之间的距离，$c_{ij}$ 表示城市 $i$ 和城市 $j$ 之间的距离。

问题的输出：

---

三种模型对比

以下是三种模型的对比：

特性	模型1	模型2	模型3
变量定义	有向弧（$x_{ij}$）	无向弧（$x_{ij}$，仅 $i<j$）	无向弧（$x_{ij}$，$x_{ij}$ 与 $x_{ji}$ 意义相同）
变量数量	$N\times (N-1)$	$N\times (N-1)/2$	$N\times (N-1)$
约束类型	显式约束	隐式无向性 + 惰性约束	显式对称约束 + 惰性约束
子环消除	静态约束	动态惰性约束	动态惰性约束
求解效率	较低	中等	较高
运行时间 (50个城市)	12.0 s	0.2 s	0.1 s

注意：

• 以上测试是基于gurobi 12.0。
• 模型1的求解效率是gurobi 12.0 高于 gurobi 11.0，模型2和3反之。

---

求解模型1

决策变量

• $x_{ij}$: 二进制变量，如果旅行者从城市 $i$ 访问城市 $j$，则 $x_{ij}=1$，否则 $x_{ij}=0$。
• $u_i$: 辅助变量，用于表示城市 $i$ 的访问顺序（顺序编号，整数）。

目标函数

最小化总的旅行距离:
$$\text{minimize}Z=\sum _{i=0}^{N-1}\sum _{j=0,j\ne i}^{N-1}{c}_{ij}\cdot x_{ij}$$

约束条件

1. 每个城市有且仅有一个出度:
   $$\sum _{j=0,j\ne i}^{N-1}{x}_{ij}=1,\forall i=0,1,\dots ,N-1$$

2. 每个城市有且仅有一个入度:
   $$\sum _{i=0,i\ne j}^{N-1}{x}_{ij}=1,\forall j=0,1,\dots ,N-1$$

3. 防止子环的约束:
   $$u_i-u_j+N\cdot x_{ij}\le N-1,\forall i,j=1,\dots ,N-1,i\ne j$$

• 注意此处 $i,j$ 是从 $1$ 开始，不是从 $0$ 开始。这是因为，假设路径为 $0\to 1\to 2\to 0$，那么 $0$ 的路径索引即低于1又高于2，出现矛盾。（~~ 在看公式的时候没注意到此处细节，复现的时候卡在这里了一段时间 ~~）

Python代码

import time
import numpy as np
from gurobipy import *
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdistdef generate_random_city_coordinates(num_cities):"""生成随机的城市坐标及距离矩阵"""np.random.seed(3)  # 锁定随机种子以确保可重复性city_coordinates = np.random.rand(num_cities, 2)    # 生成随机城市坐标（0到1之间的浮点数）c = cdist(city_coordinates, city_coordinates, metric='euclidean') # 计算城市之间的欧几里得距离return c,city_coordinatesdef plot_route(city_coordinates, solution):"""可视化城市和路径"""# 画出路径plt.plot(city_coordinates[solution][:, 0], city_coordinates[solution][:, 1], color='black', marker='o')plt.plot([city_coordinates[solution[0], 0], city_coordinates[solution[-1], 0]],[city_coordinates[solution[0], 1], city_coordinates[solution[-1], 1]], color='black', marker='o')  # 回到起点# 去掉坐标轴黑框ax = plt.gca()ax.spines['top'].set_color('none')ax.spines['right'].set_color('none')ax.spines['left'].set_color('none')ax.spines['bottom'].set_color('none')# 隐藏坐标轴刻度ax.xaxis.set_ticks_position('none')ax.yaxis.set_ticks_position('none')# 隐藏坐标轴刻度标签ax.set_xticks([]) ax.set_yticks([])plt.show()def solve_tsp(num_cities):"""解决旅行商问题 (TSP)"""# 生成距离矩阵c,city_coordinates = generate_random_city_coordinates(num_cities)# 创建模型TSP = Model("Traveling Salesman Problem")# 定义决策变量x = TSP.addVars(num_cities, num_cities, vtype=GRB.BINARY, name='visit')  # 边的访问变量u = TSP.addVars(num_cities, vtype=GRB.INTEGER, lb=0, name='aux')  # 辅助变量用于限制子环# 设置目标函数：最小化总的旅行距离TSP.setObjective(quicksum(x[i, j] * c[i, j] for i in range(num_cities) for j in range(num_cities)), GRB.MINIMIZE)# 设置约束条件# 1. 每个城市有且仅有一个出度TSP.addConstrs(quicksum(x[i, j] for j in range(num_cities) if i != j) == 1 for i in range(num_cities))# 2. 每个城市有且仅有一个入度TSP.addConstrs(quicksum(x[j, i] for j in range(num_cities) if i != j) == 1 for i in range(num_cities))# 3. 防止子环的约束TSP.addConstrs(u[i] - u[j] + x[i, j] * num_cities <= num_cities - 1 for i in range(1, num_cities) for j in range(1, num_cities) if i != j)# 求解模型TSP.optimize()if TSP.status == GRB.OPTIMAL:print("找到最优解。")# 输出选定的路径route = []for i in range(num_cities):for j in range(num_cities):if x[i, j].x > 0.5:  # 判断是否选择了这条边route.append((i, j))# 寻找完整路径current_city = 0solution = [current_city]while True:current_city = next((j for i,j in route if i == current_city ),None)solution.append(current_city)if current_city == 0:breakprint("最优路径为路径为：","->".join(map(str,solution)))print(f"总旅行距离: {TSP.ObjVal:.2f}")plot_route(city_coordinates, solution)return TSP# 主程序入口
if __name__ == "__main__":start_time = time.time()            # 标记开始时间number_of_city_coordinates = 50               # 城市数量solve_tsp(number_of_city_coordinates)         # 调用解决TSP的函数runtime = time.time() - start_time  # 计算运行时间print(f"程序运行时间: {runtime:.2f}秒")

---

求解模型2

决策变量

• $x_{ij}$: 二进制变量，如果旅行者经过城市 $i$ 和城市 $j$ 之间的弧，则 $x_{ij}=1$，否则 $x_{ij}=0$。 （注意：仅定义 $i<j$ 的边，即 $x_{ij}$ 表示无向边。）

目标函数

最小化总的旅行距离:
$$\text{minimize}Z=\sum _{i=0}^{N-1}\sum _{j=i+1}^{N-1}{c}_{ij}\cdot x_{ij}$$

约束条件

1. 每个城市有且仅有两个邻接边（入度和出度之和为2）:
   $$\sum _{j=0,j<i}^{N-1}{x}_{ji}+\sum _{j=0,j>i}^{N-1}{x}_{ij}=2,\forall i=0,1,\dots ,N-1$$

2. 防止子环的约束（作为惰性约束，通过回调函数动态添加）:
   $$\sum _{i,j\in S}{x}_{ij}\le |S|-1.$$
   （其中 $S$ 是任意子集，表示子环中的城市集合。）

Python代码

• 注意代码中的 tsp_model.update() 在Gurobi 12.0 中必须加，在 Gurobi 11.0 中可加可不加。（~~ 此处又是卡时间的一个小坑 ~~）

import time
import math
import numpy as np
import gurobipy as gp
import matplotlib.pyplot as plt
from itertools import combinations,permutations
from gurobipy import *  # 使用gurobipy库
from scipy.spatial.distance import cdistdef generate_random_cities(num_cities):"""生成随机的城市坐标及距离矩阵"""np.random.seed(3)  # 锁定随机种子以确保可重复性city_coordinates = np.random.rand(num_cities, 2)    # 生成随机城市坐标（0到1之间的浮点数）c = cdist(city_coordinates, city_coordinates, metric='euclidean') # 计算城市之间的欧几里得距离return c,city_coordinates# 计算两个城市之间的距离
def distance(city_index1, city_index2, distance_matrix):"""计算城市 city_index1 和 city_index2 之间的距离"""return distance_matrix[city_index1, city_index2]def plot_route(city_coordinates, solution):"""可视化城市和路径"""# 画出路径plt.plot(city_coordinates[solution][:, 0], city_coordinates[solution][:, 1], color='black', marker='o')plt.plot([city_coordinates[solution[0], 0], city_coordinates[solution[-1], 0]],[city_coordinates[solution[0], 1], city_coordinates[solution[-1], 1]], color='black', marker='o')  # 回到起点# 去掉坐标轴黑框ax = plt.gca()ax.spines['top'].set_color('none')ax.spines['right'].set_color('none')ax.spines['left'].set_color('none')ax.spines['bottom'].set_color('none')# 隐藏坐标轴刻度ax.xaxis.set_ticks_position('none')ax.yaxis.set_ticks_position('none')# 隐藏坐标轴刻度标签ax.set_xticks([]) ax.set_yticks([])plt.show()# 创建 Gurobi 模型
def create_model(num_cities, distance_matrix):"""创建旅行商问题的 Gurobi 模型"""model = Model("Traveling Salesman Problem")# 定义变量：只使用单向变量 (i < j)city_pairs = list(combinations(range(num_cities), 2))vars = model.addVars(city_pairs, vtype = GRB.BINARY, name='x')# 每个城市的边数为 2model.addConstrs(vars.sum(c, '*') + vars.sum('*', c) == 2 for c in range(num_cities))# 设置目标函数：最小化总的旅行距离model.setObjective(quicksum(vars[i, j] * distance_matrix[i, j] for i, j in city_pairs), GRB.MINIMIZE)model.update()return model, vars# 回调函数 - 用于消除子巡环
def subtourelim(model, where):"""回调函数，用于切断子巡环"""if where == GRB.Callback.MIPSOL:vals = model.cbGetSolution(model._vars)  # 获取当前解中选择的边selected_edges = gp.tuplelist((i, j) for (i, j), val in vals.items() if val > 0.5)tour = find_shortest_subtour(selected_edges)  # 寻找短子巡环if len(tour) < len(capitals):pairs_tour = [(tour[i], tour[i+1]) if tour[i] < tour[i+1] else (tour[i+1], tour[i]) for i in range(len(tour) - 1)]if tour[-1] < tour[0]:pairs_tour.append((tour[-1],tour[0]))else:pairs_tour.append((tour[0],tour[-1]))# 对于子巡环中的每对城市，添加子巡环消除约束model.cbLazy(gp.quicksum(model._vars[i, j] for i, j in pairs_tour) <= len(pairs_tour) - 1)# 寻找给定边的最短子巡环
def find_shortest_subtour(edges):unvisited = capitals[:]shortest_cycle = capitals[:]  # 初始占位，后续会替换while unvisited:this_cycle = []neighbors = unvisitedwhile neighbors:current = neighbors[0]this_cycle.append(current)unvisited.remove(current)neighbors = [j for i, j in edges.select(current, '*') if j in unvisited] + [i for i, j in edges.select('*', current) if i in unvisited]if len(this_cycle) <= len(shortest_cycle):shortest_cycle = this_cycle  # 更新为新的最短子巡环return shortest_cycle# 主程序
if __name__ == "__main__":start_time = time.time()  # 开始时间记录number_of_cities = 50  # 城市数量capitals = list(range(number_of_cities))# 生成随机城市坐标和距离矩阵distance_matrix,city_coordinates = generate_random_cities(number_of_cities)# 创建模型并优化tsp_model, vars = create_model(number_of_cities, distance_matrix)tsp_model._vars = varstsp_model.Params.lazyConstraints = 1  # 启用懒约束tsp_model.optimize(subtourelim)  # 优化模型# 检查模型状态，输出结果if tsp_model.status == GRB.OPTIMAL:print("找到最优解！")selected_edges = [(i, j) for i, j in vars.keys() if vars[i, j].x > 0.5]shortest_cycle = find_shortest_subtour(gp.tuplelist(selected_edges))total_distance = tsp_model.ObjValprint('最优路径', shortest_cycle)print('最优长度', total_distance)plot_route(city_coordinates, shortest_cycle)else:print("未找到可行解或模型求解失败。")elapsed_time = time.time() - start_time  # 计算运行时间print(f"程序运行时间: {elapsed_time:.2f}秒")

---

求解模型3

• 该模型来自于此处：https://github.com/Gurobi/modeling-examples/blob/master/traveling_salesman/tsp.ipynb。
• 因为感觉该模型有些冗余，所以将该模型改写为模型2，然而还是原模型的求解效率高。

决策变量

• $x_{ij}$: 二进制变量，如果旅行者经过城市 $i$ 和城市 $j$ 之间的弧，则 $x_{ij}=1$，否则 $x_{ij}=0$。
  （注意：$x_{ij}$ 和 $x_{ji}$ 是两个独立的变量，表示相同的无向弧。）

目标函数

最小化总的旅行距离:
$$\text{minimize}Z=\sum _{i=0}^{N-1}\sum _{j=0,j\ne i}^{N-1}{c}_{ij}\cdot x_{ij}$$

约束条件

1. 每个城市有且仅有两个邻接边（入度和出度之和为2）:
   $$\sum _{j=0,j\ne i}^{N-1}{x}_{ij}+\sum _{j=0,j\ne i}^{N-1}{x}_{ji}=2,\forall i=0,1,\dots ,N-1$$

2. 双向路径对称性约束:
   $$x_{ij}=x_{ji},\forall i,j=0,1,\dots ,N-1,i\ne j$$

3. 防止子环的约束（作为惰性约束，通过回调函数动态添加）:
   $$\sum _{i,j\in S}{x}_{ij}\le |S|-1.$$
   （其中 $S$ 是任意子集，表示子环中的城市集合。）

Python代码

import time
import numpy as np
import gurobipy as gp
import matplotlib.pyplot as plt
from itertools import combinations,permutations
from gurobipy import *  # 使用gurobipy库
from scipy.spatial.distance import cdistdef generate_random_cities(num_cities):"""生成随机的城市坐标及距离矩阵"""np.random.seed(3)  # 锁定随机种子以确保可重复性city_coordinates = np.random.rand(num_cities, 2)    # 生成随机城市坐标（0到1之间的浮点数）c = cdist(city_coordinates, city_coordinates, metric='euclidean') # 计算城市之间的欧几里得距离return c,city_coordinates# 计算两个城市之间的距离
def distance(city_index1, city_index2, distance_matrix):"""计算城市 city_index1 和 city_index2 之间的距离"""return distance_matrix[city_index1, city_index2]def plot_route(city_coordinates, solution):"""可视化城市和路径"""# 画出路径plt.plot(city_coordinates[solution][:, 0], city_coordinates[solution][:, 1], color='black', marker='o')plt.plot([city_coordinates[solution[0], 0], city_coordinates[solution[-1], 0]],[city_coordinates[solution[0], 1], city_coordinates[solution[-1], 1]], color='black', marker='o')  # 回到起点# 去掉坐标轴黑框ax = plt.gca()ax.spines['top'].set_color('none')ax.spines['right'].set_color('none')ax.spines['left'].set_color('none')ax.spines['bottom'].set_color('none')# 隐藏坐标轴刻度ax.xaxis.set_ticks_position('none')ax.yaxis.set_ticks_position('none')# 隐藏坐标轴刻度标签ax.set_xticks([]) ax.set_yticks([])plt.show()# 创建 Gurobi 模型
def create_model(num_cities, distance_matrix):"""创建旅行商问题的Gurobi模型"""# 创建模型tsp_model = Model("Traveling Salesman Problem")# 单向城市对Pairings = combinations(range(num_cities), 2)# 双向城市对city_pairs = list(permutations(range(num_cities), 2))# 添加变量：城市 i 和城市 j 是否相邻vars = tsp_model.addVars(city_pairs, vtype=GRB.BINARY, name='x')# 每个城市的边数为 2tsp_model.addConstrs(vars.sum(c, '*') == 2 for c in range(num_cities))# 无向边tsp_model.addConstrs(vars[i,j]==vars[j,i] for i,j in Pairings)# 设置目标函数：最小化2倍的总旅行距离tsp_model.setObjective(quicksum(vars[i, j] * distance(i, j, distance_matrix) for i,j in city_pairs), GRB.MINIMIZE)tsp_model.update()return tsp_model, vars# 回调函数 - 用于消除子巡环
def subtourelim(model, where):"""回调函数，用于切断子巡环"""if where == GRB.Callback.MIPSOL:vals = model.cbGetSolution(model._vars)  # 获取当前解中选择的边selected_edges = gp.tuplelist((i, j) for i, j in model._vars.keys() if vals[i, j] > 0.5)tour = find_shortest_subtour(selected_edges)  # 寻找短子巡环if len(tour) < len(capitals):# 对于子巡环中的每对城市，添加子巡环消除约束model.cbLazy(gp.quicksum(model._vars[i, j] for i, j in combinations(tour, 2)) <= len(tour) - 1)# 寻找给定边的最短子巡环
def find_shortest_subtour(edges):unvisited = capitals[:]shortest_cycle = capitals[:]  # 初始占位，后续会替换while unvisited:this_cycle = []neighbors = unvisitedwhile neighbors:current = neighbors[0]this_cycle.append(current)unvisited.remove(current)neighbors = [j for i, j in edges.select(current, '*') if j in unvisited]if len(this_cycle) <= len(shortest_cycle):shortest_cycle = this_cycle  # 更新为新的最短子巡环return shortest_cycle# 主程序
if __name__ == "__main__":start_time = time.time()  # 开始时间记录number_of_cities = 50  # 城市数量capitals = list(range(number_of_cities))# 生成随机城市坐标和距离矩阵distance_matrix,city_coordinates = generate_random_cities(number_of_cities)# 创建模型并优化tsp_model, vars = create_model(number_of_cities, distance_matrix)tsp_model._vars = varstsp_model.Params.lazyConstraints = 1  # 启用懒约束tsp_model.optimize(subtourelim)  # 优化模型# 检查模型状态，输出结果if tsp_model.status == GRB.OPTIMAL:print("找到最优解！")selected_edges = [(i, j) for i, j in vars.keys() if vars[i, j].x > 0.5]shortest_cycle = find_shortest_subtour(gp.tuplelist(selected_edges))total_distance = tsp_model.ObjVal/2print('最优路径', shortest_cycle)print('最优长度', total_distance)plot_route(city_coordinates, shortest_cycle)else:print("未找到可行解或模型求解失败。")elapsed_time = time.time() - start_time  # 计算运行时间print(f"程序运行时间: {elapsed_time:.2f}秒")

---

三个模型的优势与不足

模型1的优势与不足

• 优势：
  • 模型简单直观，易于实现。
  • 使用静态约束消除子环，适合初学者理解。
• 不足：
  • 静态约束可能导致松弛解质量较差，求解效率低。
  • 变量和约束数量较多。

模型2的优势与不足

• 优势：
  • 变量数量较少，约束较少。
• 不足：
  • 回调函数需要额外处理单向变量的双向性，效率略低。
  • 隐式无向性可能导致求解器无法充分利用对称性优化。

模型3的优势与不足

• 优势：
  • 显式对称约束帮助求解器更快识别问题结构，求解效率高。
  • 回调函数直接处理双向变量，效率更高。
  • 适合大规模问题。
• 不足：
  • 模型不简洁，变量数量较多（$N\times (N-1)$），但通过显式约束弥补了效率损失。