使用gurobi解决旅行商问题（TSP）
TSP问题
旅行商问题（Traveling Salesman Problem）是最著名的组合优化问题之一。这个问题很容易解释，但解决起来非常复杂——即使是在少数城市的情况下也是如此。
问题描述
TSP 可以定义如下：对于给定的城市列表和每对城市之间的距离，我们希望找到到达每个城市一次并返回原点城市的最短路线。
有一类旅行商问题是这样的：假设从城市$i$到城市$j$的距离和城市$j$到城市$i$的距离是相等的，这种类型的旅行商问题也称为对称旅行商问题。在接下来的例子中我们使用欧几里得距离，模型是有效的，与确定距离的方式无关。
解决方法
数学编程是一种声明性方法，建模者制定数学优化模型，该模型捕获复杂决策问题的关键方面。
数学优化模型有五个组成部分，即：
- 集合和索引。
- 参数。
- 决策变量。
- 目标函数。
- 约束。
TSP模型制定
集合和索引
$i,j\in Capitals$:索引以及美国首府城市的集合。
Pairings=$\{(i,j)\in Capitals\times Capitals\}$:可采取的组合集合。
$S\subset Capticals$:一个美国首府城市集合的子集。
$G=(Capitals,Pairings)$:一幅由集合Capticals定义的点、集合Pairings定义的边构成的图。
参数
$d_{i,j}\in \mathbb{R}^+$: 从城市$i$到城市$j$的距离，其中$(i,j)\in Pairings$。
注意：从城市$i$到城市$j$的距离和城市$j$到城市$i$的距离是相等的，也就是$d_{i,j}=d_{j,i}$,因此，这个TSP也被称为对称旅行推销员问题。
决策变量
$x_{i,j}\in\{0,1\}$：如果我们决定连接城市$i$与城市$j$，则此变量等于1，否则，决策变量等于零。
目标函数
最短路线，最小化路线的总距离。路线是一系列首都城市，销售人员只访问每个城市一次，然后返回起始的首都城市。
$$ \begin{align} \mathrm{Min}\quad Z=\sum\limits_{(i,j)\in\mathrm{Pairings}}d_{i,j}\cdot x_{i,j} \end{align} $$
约束
- 对称约束。对于每条边$(i,j)$，确保首府$i$和$j$是相连接的，如果前者在访问后者之前或之后，需立即被访问（ps:这句话翻译的不好）。
$$ \begin{align} x_{i,j}=x_{j,i}\forall(i,j)\in Pairings \end{align} $$
进出首都。对于每个首都$i$，确保该城市与其他两个城市相连.
$$ \begin{align} \sum\limits_{(i,j)\in\mathrm{Pairings}}x_{i,j}=2\forall i\in Capitals \end{align} $$
- 子回路的消除。这些约束确保了对于任何属于集合$Cpitals$的城市子集$S$，是没有循环的。也就是说，不存在访问子集中所有城市并返回原始城市的路线。
$$ \begin{align} \sum\limits_{(i\neq j)\in S}\boldsymbol{x}_{i,j}\leq|S|-1\forall\quad S\subset Capitals \end{align} $$
Python 实现
此建模示例需要导入以下库：
- math 访问数学函数。
- itertools 实现了许多迭代器构建块。
- folium 创建地图。
gurobipy 调用Gurobi算法来求解MIP模型。
读取和输入数据
从json文件中读取大写名称和坐标。
import json # Read capital names and coordinates from json file capitals_json = json.load(open('capitals.json')) capitals = [] coordinates = {} for state in capitals_json: if state not in ['AK', 'HI']: capital = capitals_json[state]['capital'] capitals.append(capital) coordinates[capital] = (float(capitals_json[state]['lat']), float(capitals_json[state]['long']))预处理
以下函数计算每对州府组合的距离。
import math from itertools import combinations,product # Compute pairwise distance matrix def distance(city1, city2): c1 = coordinates[city1] c2 = coordinates[city2] diff = (c1[0]-c2[0], c1[1]-c2[1]) return math.sqrt(diff[0]*diff[0]+diff[1]*diff[1]) dist = {(c1, c2): distance(c1, c2) for c1, c2 in product(capitals, capitals) if c1 != c2}回调定义
以下函数确定延迟约束以消除子回路。
# Callback - use lazy constraints to eliminate sub-tours def subtourelim(model, where): if where == GRB.Callback.MIPSOL: # make a list of edges selected in the solution vals = model.cbGetSolution(model._vars) selected = gp.tuplelist((i, j) for i, j in model._vars.keys() if vals[i, j] > 0.5) # find the shortest cycle in the selected edge list tour = subtour(selected) if len(tour) < len(capitals): # add subtour elimination constr. for every pair of cities in subtour model.cbLazy(gp.quicksum(model._vars[i, j] for i, j in combinations(tour, 2)) <= len(tour)-1)查找最短路线
以下函数确定从给定的一组边开始的最短路径。
# Given a tuplelist of edges, find the shortest subtour def subtour(edges): unvisited = capitals[:] cycle = capitals[:] # Dummy - guaranteed to be replaced while unvisited: # true if list is non-empty thiscycle = [] neighbors = unvisited while neighbors: current = neighbors[0] thiscycle.append(current) unvisited.remove(current) neighbors = [j for i, j in edges.select(current, '*') if j in unvisited] if len(thiscycle) <= len(cycle): cycle = thiscycle # New shortest subtour return cycle模型部署
我们现在通过定义决策变量、约束和目标函数来确定TSP的模型。接下来，我们开始优化，Gurobi找到TSP的最佳路径。
import gurobipy as gp from gurobipy import GRB # tested with Python 3.7 & Gurobi 9.0.0 m = gp.Model() # Variables: is city 'i' adjacent to city 'j' on the tour? vars = m.addVars(dist.keys(), obj=dist, vtype=GRB.BINARY, name='e') for i, j in vars.keys(): m.addConstr(vars[j, i] == vars[i, j]) # edge in opposite direction # Constraints: two edges incident to each city m.addConstrs(vars.sum(c, '*') == 2 for c in capitals) # Optimize the model m._vars = vars m.Params.lazyConstraints = 1 m.optimize(subtourelim)分析
我们检索TSP的最优解，并验证最优路线（或旅游）到达所有城市并返回到原始城市。
# Retrieve solution vals = m.getAttr('x', vars) selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5) tour = subtour(selected) assert len(tour) == len(capitals)最佳路线显示在下图中。
# Map the solution import folium map = folium.Map(location=[40,-95], zoom_start = 4) points = [] for city in tour: points.append(coordinates[city]) points.append(points[0]) folium.PolyLine(points).add_to(map) map
总结
旅行商问题（TSP）是最流行的组合优化问题。这个问题很容易解释，尽管解决起来很复杂。解决的最大TSP问题有85900个城市。TSP是解决复杂组合优化问题的新方法的发现来源，并导致了许多应用。
在这个建模示例中，我们展示了如何将对称旅行推销员问题公式化为MIP问题。我们还展示了如何通过使用惰性约束来动态消除子任务。参考文献
[1] D. L. Applegate, R. E. Bixby, V. Chvatal and W. J. Cook , The Traveling Salesman Problem: A Computational Study, Princeton University Press, Princeton, 2006.
[2] http://www.math.uwaterloo.ca/tsp/index.html
