在Python中，根据转移概率为TSP生成随机路径
In Python, generate a random path for the TSP according to transition probabilities
对于旅行商问题 (TSP)，我想在 Python 中从状态转换矩阵 M ( n x n)，其中 M[i,j] 包含最短全局路径将节点 i 连接到节点 j 的概率。
有人知道如何在 Python 中做到这一点吗？我在这里询问执行此操作的一般方法或模块。
示例：假设 M[i,j]=1 for 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。
（如果您提供一些代码，这将是一个更好的问题；特别是因为该问题针对的是一个之前有很多步骤的步骤）
这里有一些方法，应该只作为一些想法，因为设置可能有点烦人（cvxpy 和一些 MIP-solver 这还不错；我在这里使用 CBC，请参阅 cvxpy 的文档备选方案或 setup-docs）。它也没有真正经过测试。这个想法有点基于概率图形模型中的 MAP-calculation，但转移可能在数学上是错误的（不保证）。这也更难，因为在我们的例子中 opt-problem 受到约束！
想法
制定一个 optimization-problem 最大化 a-priori-probabilities 使用的 = solution-path 的一部分（并且平方这个 = 更大的偏差受到更严重的惩罚），同时生成有效的解决方案（有效路径） .
虽然问题可能是 non-convex（我不确定），因此无法以 globally-optimal 方式解决），但我们正在使用一些 well-analyzed 启发式 此处凸规划的区别。
备注： 这种方法是设计用来搜索 global-optimum（不完全正确，因为我们使用的是 nonconvex-optimization 算法）。这意味着，在这种方法中，使用不同解决方案的多次采样并不是很自然（但是使用不同的 starting-points 调用 DCCP 例程将很可能导致不同的解决方案）。
备注 2： non-small 实例（使用 non-commercial 求解器）的性能非常差，这使得它更 theoretical-approach。
实施
这里是使用Python3、numpy、scipy（最短路径）、cvxpy（opt-problem公式）和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))
输出
...
...
iteration= 1 cost value = 0.34508891154470694 tau = 0.005
iteration= 2 cost value = 0.5092119781611304 tau = 0.006
iteration= 3 cost value = 0.5092119781611304 tau = 0.0072
Solution path matrix
[[0 0 0 0 0 0 0 0 1 0]
[1 0 0 0 0 0 0 0 0 0]
[0 0 0 1 0 0 0 0 0 0]
[0 0 0 0 1 0 0 0 0 0]
[0 0 0 0 0 1 0 0 0 0]
[0 0 0 0 0 0 0 0 0 1]
[0 0 1 0 0 0 0 0 0 0]
[0 0 0 0 0 0 1 0 0 0]
[0 0 0 0 0 0 0 1 0 0]
[0 1 0 0 0 0 0 0 0 0]]
A-priori probability matrix
[[ 0. 0.14 0. 0.24 0.11 0.07 0.07 0.03 0.13 0.21]
[ 0.12 0. 0.09 0.22 0.01 0.17 0.13 0.11 0.06 0.07]
[ 0.11 0.1 0. 0.27 0.08 0.12 0.05 0.02 0.1 0.15]
[ 0.06 0.17 0.06 0. 0.15 0.12 0.11 0.09 0.01 0.23]
[ 0.09 0.16 0.09 0.18 0. 0.13 0.13 0.11 0.05 0.05]
[ 0.01 0.15 0.02 0.21 0.12 0. 0.08 0.05 0.15 0.22]
[ 0.06 0.17 0.06 0.25 0.03 0.12 0. 0.09 0.11 0.11]
[ 0.09 0.07 0.07 0.21 0.08 0.08 0.11 0. 0.14 0.15]
[ 0.11 0.16 0.11 0.09 0.07 0.14 0.11 0.12 0. 0.1 ]
[ 0.07 0.17 0.07 0.21 0.15 0.12 0.12 0.09 0. 0. ]]
编辑：
- 哎呀，不知何故 ECOS_BB 似乎仍在使用，这告诉我们更好的 solver-setup 有更多的潜力。
对于旅行商问题 (TSP)，我想在 Python 中从状态转换矩阵 M ( n x n)，其中 M[i,j] 包含最短全局路径将节点 i 连接到节点 j 的概率。
有人知道如何在 Python 中做到这一点吗？我在这里询问执行此操作的一般方法或模块。
示例：假设 M[i,j]=1 for 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。
（如果您提供一些代码，这将是一个更好的问题；特别是因为该问题针对的是一个之前有很多步骤的步骤）
这里有一些方法，应该只作为一些想法，因为设置可能有点烦人（cvxpy 和一些 MIP-solver 这还不错；我在这里使用 CBC，请参阅 cvxpy 的文档备选方案或 setup-docs）。它也没有真正经过测试。这个想法有点基于概率图形模型中的 MAP-calculation，但转移可能在数学上是错误的（不保证）。这也更难，因为在我们的例子中 opt-problem 受到约束！
想法
制定一个 optimization-problem 最大化 a-priori-probabilities 使用的 = solution-path 的一部分（并且平方这个 = 更大的偏差受到更严重的惩罚），同时生成有效的解决方案（有效路径） .
虽然问题可能是 non-convex（我不确定），因此无法以 globally-optimal 方式解决），但我们正在使用一些 well-analyzed 启发式 此处凸规划的区别。
备注： 这种方法是设计用来搜索 global-optimum（不完全正确，因为我们使用的是 nonconvex-optimization 算法）。这意味着，在这种方法中，使用不同解决方案的多次采样并不是很自然（但是使用不同的 starting-points 调用 DCCP 例程将很可能导致不同的解决方案）。
备注 2： non-small 实例（使用 non-commercial 求解器）的性能非常差，这使得它更 theoretical-approach。
实施
这里是使用Python3、numpy、scipy（最短路径）、cvxpy（opt-problem公式）和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))
输出
...
...
iteration= 1 cost value = 0.34508891154470694 tau = 0.005
iteration= 2 cost value = 0.5092119781611304 tau = 0.006
iteration= 3 cost value = 0.5092119781611304 tau = 0.0072
Solution path matrix
[[0 0 0 0 0 0 0 0 1 0]
[1 0 0 0 0 0 0 0 0 0]
[0 0 0 1 0 0 0 0 0 0]
[0 0 0 0 1 0 0 0 0 0]
[0 0 0 0 0 1 0 0 0 0]
[0 0 0 0 0 0 0 0 0 1]
[0 0 1 0 0 0 0 0 0 0]
[0 0 0 0 0 0 1 0 0 0]
[0 0 0 0 0 0 0 1 0 0]
[0 1 0 0 0 0 0 0 0 0]]
A-priori probability matrix
[[ 0. 0.14 0. 0.24 0.11 0.07 0.07 0.03 0.13 0.21]
[ 0.12 0. 0.09 0.22 0.01 0.17 0.13 0.11 0.06 0.07]
[ 0.11 0.1 0. 0.27 0.08 0.12 0.05 0.02 0.1 0.15]
[ 0.06 0.17 0.06 0. 0.15 0.12 0.11 0.09 0.01 0.23]
[ 0.09 0.16 0.09 0.18 0. 0.13 0.13 0.11 0.05 0.05]
[ 0.01 0.15 0.02 0.21 0.12 0. 0.08 0.05 0.15 0.22]
[ 0.06 0.17 0.06 0.25 0.03 0.12 0. 0.09 0.11 0.11]
[ 0.09 0.07 0.07 0.21 0.08 0.08 0.11 0. 0.14 0.15]
[ 0.11 0.16 0.11 0.09 0.07 0.14 0.11 0.12 0. 0.1 ]
[ 0.07 0.17 0.07 0.21 0.15 0.12 0.12 0.09 0. 0. ]]
编辑：
- 哎呀，不知何故 ECOS_BB 似乎仍在使用，这告诉我们更好的 solver-setup 有更多的潜力。