在Python中，根据转移概率生成TSP的随机路径

Question

对于旅行商问题（TSP），我想用Python从状态转移矩阵M（nxn）生成一个随机旅行（其中n组中的每个节点只访问一次），其中M[I，j]包含最短全局路径将节点I连接到节点j的概率

有人知道怎么用Python做到这一点吗？我在这里问的是实现这一点的一般方法或模块。在

示例：假设M[i，j]=1表示j=（i+1）%n，其他地方为0。将生成的所有路径（总是从0开始）是：（0,1,2，…，n）。如果你稍微改变一下这个矩阵，用1.0代替0.9，把0.1放在M[i，i+2]上，可能的路径是：（0,2,3，…，n，1）。在这个概率为0的特殊例子中，我知道最后一次移动是不可能的（从节点n到节点1），因此可以假设概率总是大于0。在

Answer 1

（如果您提供一点代码，这会是一个更好的问题；特别是因为问题的目标是一个之前有很多步骤的步骤）

这里有一些方法，应该只作为一些想法，因为设置可能有点烦人（cvxpy和一些MIP解算器不错；我在这里使用CBC，请参阅cvxpy的文档以获取替代或设置文档）。也没有经过真正的测试。这个想法在某种程度上是基于概率图形模型中的MAP计算，但是传输可能在数学上是错误的（没有保证）。这也更难，因为在我们的例子中，opt问题是受约束的！

想法

制定一个优化问题，使所用的先验概率最大化=解路径的一部分（平方这个=偏差越大，惩罚越重），同时生成一个有效解（有效路径）。

虽然问题可能是非凸的（我不确定），因此不可能以全局最优的方式解决，但我们在这里使用了一些分析良好的启发式方法来求解凸规划的差分。

备注：这种方法是通过设计搜索全局最优值（并非完全正确，因为我们使用的是非凸优化算法）。这意味着，在这种方法中，使用不同的解决方案进行多个采样并不是很自然（但是使用不同的起点调用DCCP例程将导致高概率的不同解决方案）。

备注2:对于非小型实例（使用非商业解决方案），性能非常差，这使得它成为一种更加理论化的方法。

实施

下面是一些使用python3、numpy、scipy（最短路径）、cvxpy（opt problem formulation）和dccp（cvxpy凸函数优化扩展的差异）的实现

代码

import numpy as np
from scipy.sparse.csgraph import csgraph_from_dense, shortest_path
from scipy.spatial import distance
from cvxpy import *
import dccp
np.random.seed(1)

""" Create random problem """
N = 10
distances = np.random.rand(N, N)                    # assymetric

""" Calculate shortest paths """
csparse_distances = csgraph_from_dense(distances)
shortest = shortest_path(csparse_distances)         # directed

""" Calculate a-prori probabilities based on global shortest paths """
shortest_global = np.copy(shortest)
for row in range(shortest_global.shape[0]):
    # normalize sum to 1
    row_sum = np.sum(shortest_global[row, :])
    shortest_global[row, :] /= row_sum

""" Formulate MIQP problem """
# Based on: https://en.wikipedia.org/wiki/Travelling_salesman_problem
#       and my example here: https://github.com/cvxgrp/cvxpy/blob/master/examples/tsp_mip.py#L16

# Variables
x = Bool(N, N)                                      # edge x,y in PATH
u = Int(N)                                          # aux-var

# Constraints
constraints = []

for j in range(N):                                      # ingoing: exactly 1
    indices = list(range(0, j)) + list(range(j + 1, N))
    constraints.append(sum_entries(x[indices, j]) == 1)
for i in range(N):
    indices = list(range(0, i)) + list(range(i + 1, N)) # outgoing: exactly 1
    constraints.append(sum_entries(x[i, indices]) == 1)

for i in range(1, N):                               # subtour-elimination
    for j in range(1, N):
        if i != j:
            constraints.append(u[i] - u[j] + N*x[i, j] <= N-1)

# Objective
obj = Maximize(sum_entries(square(mul_elemwise(shortest_global, x))))

# Solve
prob = Problem(obj, constraints)
print("problem is DCP: ", prob.is_dcp())
prob.solve(method='dccp', solver=CBC, ccp_times=10)  # do not use default solver!
# Remark: formulation above not accepted by CVX-ruleset
#         -> need "difference of convex function"-extension
#         -> heuristic (which is well-known for good behaviour)!


""" Print solution """
print('Solution path matrix')
print(np.round(x.value).astype('int'))
print('A-priori probability matrix')
print(np.round(shortest_global, 2))

输出

^{pr2}$

编辑：

• 哎哟，不知何故ECOS_-BB似乎仍在使用，这告诉我们有更好的解算器设置的潜力更大。在