TSP

1. 介绍

问题描述：

给定一个有向完全图 G = (V, E)

有向完全图：一个节点可以直接到达其他节点

V：图中节点集合 |V| = N

E：图中边集合

当边 eij 与边 eji 对称时，称为对称 TSP 问题

目标：

找到一条从起点（任意点）出发依次不重复地经过所有其他节点，最终返回起点的最小权和（最短）的路径

输入： 一系列坐标点的集合

输出： 一个访问序列。每个节点被且只被访问一次；每个节点被且只被离开一次。

2. 建模

变量

模型框架

公式解释

• (1) 最小化解的总权重，在最优解中所有边的权重相加
• (2) 访问约束。被且只被访问一次。对于任意 j 节点，最优解中流入该节点的边总数等于 1
• (3) 离开约束。被且只被离开一次。对于任意 i 节点，最优接种流出该节点的边的总数等于 1
• (4) 决策变量的 0-1

不足

求解时不会出现一个环路，而是出现若干个子圈，子圈的权重之和与大圈相比更小

改进：破子圈约束

方法：MTZ 消除环路约束

加入辅助变量。破坏所有环路的存在。对每个节点引入一个决策变量 $\mu_i$ ，$\forall i \in V, \mu_i \ge 0$。
$$\mu_i - \mu_j + N\times x_{ij} \le N-1, \quad \forall i,j\in V, i,j\ne 0,i\ne j $$
其中，$\mu_i$ 是节点 i 被访问的次序，正常情况下，先访问 i 后访问 j，$\mu_i-\mu_j=-1$，若 j 被访问过了，j 的序号比 i 小，那么 $\mu_i-\mu_j>0$

最终模型

目标

• $$min \sum_{i\in V} \sum_{j \in V} c_{ij}x_{ij}$$

约束

• $$\sum_{i \in V} x_{ij} =1, \quad\quad \forall j \in V, i\ne j$$
• $$\sum_{j\in V}x_{ij} = 1, \quad\quad \forall i \in V, i\ne j$$
• $$\mu_i - \mu_j + N\times x_{ij} \le N-1, \quad\quad \forall i,j\in V, i,j\ne 0,i\ne j$$
• $$x_{ij} \in \{0, 1\}, \quad\quad \forall i,j \in V, i\ne j$$
• $$\mu_i \ge 0, \quad \mu_i \in R, \quad \forall i\in V$$

实验

随机生成30个点，计算两两节点 i 和 j 之间的距离，若 i=j 则将距离设为无限大。使用MTZ方法破子圈。

代码展示

import gurobipy as gp
import numpy as np
import matplotlib.pyplot as plt

# ------------------------------------------ #
# 数据集
# ------------------------------------------ #

coords = []  # 保存所有的节点坐标
n = 0  # 记录产生节点的个数
while n<30:
    coord = [np.random.randint(1,10), np.random.randint(1,10)]
    if coord not in coords:
        coords.append(coord)
        n += 1
grids = gp.tuplelist(coords)  # 转换为gurobi列表类型
N = len(grids)  # 节点个数
# 创建距离矩阵
dis_matrix = np.zeros(shape=[N,N])
for i in range(N):
    for j in range(N):
        if i==j:  # 如果节点自身到自身
            dis_matrix[i][j] = 1000000  # 无限远
        else:
            dis_matrix[i][j] = np.sqrt((grids[i][0]-grids[j][0])**2 \
                                + (grids[i][1]-grids[j][1])**2)

# ------------------------------------------ #
# 模型构建
# ------------------------------------------ #
m = gp.Model()

# 变量
x = m.addMVar((N,N), lb=0, ub=1, vtype=gp.GRB.BINARY)  # 代表i和j节点之间是否可达
u = m.addVars((N), vtype=gp.GRB.CONTINUOUS)  # 给每个节点一个编号 len=N
m.update()

# 目标
# m.setObjective(gp.quicksum(dis_matrix[i][j]*x[i][j] for j in range(N) for i in range(N)), gp.GRB.MINIMIZE)
obj = 0
for i in range(N):
    for j in range(N):
        obj += x[i][j] * dis_matrix[i][j]
m.setObjective(obj, gp.GRB.MINIMIZE)

# 约束
# m.addConstr(gp.quicksum(x[:,j] for j in range(N))==1, name='con_in')  # 流入约束
# m.addConstr(gp.quicksum(x[i,:] for i in range(N))==1, name='con_out')  # 流出约束
# 流入约束
for j in range(N):
    m.addConstr(gp.quicksum(x[:,j])==1)
# 流出约束
for i in range(N):
    m.addConstr(gp.quicksum(x[i,:])==1)


# 破子圈约束
for i in range(1,N):
    for j in range(1,N):
        if i != j:
            m.addConstr(u[i]-u[j]+N*x[i][j]<=N-1)

# 求解
m.optimize()

# 结果展示
var_lst = m.getVars()  # 变量值
for i in range(0, len(var_lst)):
    if var_lst[i].X != 0:
        print(var_lst[i].VarName, var_lst[i].X)
print("目标函数值：", m.ObjVal)
best=(x.X).astype(int)
print(f"最优分配为：\n {best}")

# 记录按顺序遍历的节点编号
ans = np.where(best==1)  # 获取被选节点的索引
ans_i = ans[0]
ans_j = ans[1]
n0 = ans_i[0]  # 起始节点
anslist = [n0]  # 记录编号
coordlist = [coords[n0]]  # 记录坐标
for _ in range(N):
    n1 = ans_j[n0]  # 下一个节点
    anslist.append(n1) 
    coordlist.append(coords[n1])
    n0 = n1  # 更新当前节点，当前走到n1

# 绘图
plt.plot(np.array(coords)[:,0], np.array(coords)[:,1], 'o')  # 散点[x,y]
plt.plot(np.array(coordlist)[:,0], np.array(coordlist)[:,1], '--')  # 折线
plt.show()