求解旅行商问题的三种精确性建模方法，性能差距巨大
文章目录
- 旅行商问题介绍
- 三种模型对比
- 求解模型1
- 决策变量
- 目标函数
- 约束条件
- Python代码
- 求解模型2
- 决策变量
- 目标函数
- 约束条件
- Python代码
- 求解模型3
- 决策变量
- 目标函数
- 约束条件
- Python代码
- 三个模型的优势与不足
旅行商问题介绍
旅行商问题 (Traveling Salesman Problem, TSP) 是一个经典的组合优化问题，目标是找到一条最短路径，该路径必须经过每个城市一次且仅一次，最终返回到起始城市。
问题的输入：
- N N N: 城市的总数。
- c c c: 城市之间的距离， c i j c_{ij} cij 表示城市 i i i 和城市 j j j 之间的距离。
问题的输出：
三种模型对比
以下是三种模型的对比：
| 特性 | 模型1 | 模型2 | 模型3 |
|---|---|---|---|
| 变量定义 | 有向弧（ x i j x_{ij} xij） | 无向弧（ x i j x_{ij} xij，仅 i < j i < j i<j） | 无向弧（ x i j x_{ij} xij， x i j x_{ij} xij 与 x j i x_{ji} xji 意义相同） |
| 变量数量 | N × ( N − 1 ) N \times (N-1) N×(N−1) | N × ( N − 1 ) / 2 N \times (N-1)/2 N×(N−1)/2 | N × ( N − 1 ) N \times (N-1) N×(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 i j x_{ij} xij: 二进制变量，如果旅行者从城市 i i i 访问城市 j j j，则 x i j = 1 x_{ij} = 1 xij=1，否则 x i j = 0 x_{ij} = 0 xij=0。
- u i u_i ui: 辅助变量，用于表示城市 i i i 的访问顺序（顺序编号，整数）。
目标函数
最小化总的旅行距离:
minimize
Z
=
∑
i
=
0
N
−
1
∑
j
=
0
,
j
≠
i
N
−
1
c
i
j
⋅
x
i
j
\text{minimize} \quad Z = \sum_{i=0}^{N-1}\sum_{j=0,j \neq i}^{N-1} c_{ij} \cdot x_{ij}
minimizeZ=i=0∑N−1j=0,j=i∑N−1cij⋅xij
约束条件
-
每个城市有且仅有一个出度:
∑ j = 0 , j ≠ i N − 1 x i j = 1 , ∀ i = 0 , 1 , … , N − 1 \sum_{j=0, j \neq i}^{N-1} x_{ij} = 1, \quad \forall i = 0, 1, \ldots, N-1 j=0,j=i∑N−1xij=1,∀i=0,1,…,N−1 -
每个城市有且仅有一个入度:
∑ i = 0 , i ≠ j N − 1 x i j = 1 , ∀ j = 0 , 1 , … , N − 1 \sum_{i=0, i \neq j}^{N-1} x_{ij} = 1, \quad \forall j = 0, 1, \ldots, N-1 i=0,i=j∑N−1xij=1,∀j=0,1,…,N−1 -
防止子环的约束:
u i − u j + N ⋅ x i j ≤ N − 1 , ∀ i , j = 1 , … , N − 1 , i ≠ j u_i - u_j + N \cdot x_{ij} \leq N - 1, \quad \forall i, j = 1, \ldots, N-1, \; i \neq j ui−uj+N⋅xij≤N−1,∀i,j=1,…,N−1,i=j
- 注意此处 i , j i,j i,j 是从 1 1 1 开始，不是从 0 0 0 开始。这是因为，假设路径为 0 → 1 → 2 → 0 0 \rightarrow 1 \rightarrow 2 \rightarrow 0 0→1→2→0，那么 0 0 0 的路径索引即低于1又高于2，出现矛盾。（~~ 在看公式的时候没注意到此处细节，复现的时候卡在这里了一段时间 ~~）
Python代码
import time
import numpy as np
from gurobipy import *
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdist
def 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_coordinates
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()
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 = 0
solution = [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:
break
print("最优路径为路径为：","->".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 i j x_{ij} xij: 二进制变量，如果旅行者经过城市 i i i 和城市 j j j 之间的弧，则 x i j = 1 x_{ij} = 1 xij=1，否则 x i j = 0 x_{ij} = 0 xij=0。 （注意：仅定义 i < j i < j i<j 的边，即 x i j x_{ij} xij 表示无向边。）
目标函数
最小化总的旅行距离:
minimize
Z
=
∑
i
=
0
N
−
1
∑
j
=
i
+
1
N
−
1
c
i
j
⋅
x
i
j
\text{minimize} \quad Z = \sum_{i=0}^{N-1}\sum_{j=i+1}^{N-1} c_{ij} \cdot x_{ij}
minimizeZ=i=0∑N−1j=i+1∑N−1cij⋅xij
约束条件
-
每个城市有且仅有两个邻接边（入度和出度之和为2）:
∑ j = 0 , j < i N − 1 x j i + ∑ j = 0 , j > i N − 1 x i j = 2 , ∀ i = 0 , 1 , … , N − 1 \sum_{j=0, j < i}^{N-1} x_{ji} + \sum_{j=0, j > i}^{N-1} x_{ij} = 2, \quad \forall i = 0, 1, \ldots, N-1 j=0,j<i∑N−1xji+j=0,j>i∑N−1xij=2,∀i=0,1,…,N−1 -
防止子环的约束（作为惰性约束，通过回调函数动态添加）:
∑ i , j ∈ S x i j ≤ ∣ S ∣ − 1. \sum_{i, j \in S} x_{ij} \leq |S| - 1. i,j∈S∑xij≤∣S∣−1.
（其中 S S 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 cdist
def 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')
# 每个城市的边数为 2
model.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 = unvisited
while 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 = vars
tsp_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
print('最优路径', 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
i
j
x_{ij}
xij: 二进制变量，如果旅行者经过城市
i
i
i 和城市
j
j
j 之间的弧，则
x
i
j
=
1
x_{ij} = 1
xij=1，否则
x
i
j
=
0
x_{ij} = 0
xij=0。
（注意： x i j x_{ij} xij 和 x j i x_{ji} xji 是两个独立的变量，表示相同的无向弧。）
目标函数
最小化总的旅行距离:
minimize
Z
=
∑
i
=
0
N
−
1
∑
j
=
0
,
j
≠
i
N
−
1
c
i
j
⋅
x
i
j
\text{minimize} \quad Z = \sum_{i=0}^{N-1}\sum_{j=0, j \neq i}^{N-1} c_{ij} \cdot x_{ij}
minimizeZ=i=0∑N−1j=0,j=i∑N−1cij⋅xij
约束条件
-
每个城市有且仅有两个邻接边（入度和出度之和为2）:
∑ j = 0 , j ≠ i N − 1 x i j + ∑ j = 0 , j ≠ i N − 1 x j i = 2 , ∀ i = 0 , 1 , … , N − 1 \sum_{j=0, j \neq i}^{N-1} x_{ij} + \sum_{j=0, j \neq i}^{N-1} x_{ji} = 2, \quad \forall i = 0, 1, \ldots, N-1 j=0,j=i∑N−1xij+j=0,j=i∑N−1xji=2,∀i=0,1,…,N−1 -
双向路径对称性约束:
x i j = x j i , ∀ i , j = 0 , 1 , … , N − 1 , i ≠ j x_{ij} = x_{ji}, \quad \forall i, j = 0, 1, \ldots, N-1, \; i \neq j xij=xji,∀i,j=0,1,…,N−1,i=j -
防止子环的约束（作为惰性约束，通过回调函数动态添加）:
∑ i , j ∈ S x i j ≤ ∣ S ∣ − 1. \sum_{i, j \in S} x_{ij} \leq |S| - 1. i,j∈S∑xij≤∣S∣−1.
（其中 S S 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 cdist
def 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')
# 每个城市的边数为 2
tsp_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 = unvisited
while 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 = vars
tsp_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/2
print('最优路径', 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 × ( N − 1 ) N \times (N-1) N×(N−1)），但通过显式约束弥补了效率损失。