(如果您提供一些代码,这将是一个更好的问题;尤其是因为该问题针对的是之前有很多步骤的步骤)
这里有一些方法,应该只作为一些想法,因为设置可能有点烦人(cvxpy 和一些不错的 MIP 求解器;我在这里使用 CBC,请参阅 cvxpy 的文档以获取替代方案或设置-文档)。它也没有真正经过测试。这个想法在某种程度上基于概率图形模型中的 MAP 计算,但转移可能在数学上是错误的(不能保证)。这也更难,因为我们案例中的 opt-problem 必然会受到约束!
想法
制定一个优化问题,最大化使用的先验概率 = 解决方案路径的一部分(并且平方 this = 更大的偏差受到更严重的惩罚),同时生成一个有效的解决方案(有效路径)。
虽然问题可能是非凸的(我不确定),因此无法以全局最优的方式解决),但我们正在使用一些经过充分分析的启发式算法来这里凸编程的差异时间>。
备注:这种方法是通过设计搜索全局最优的(不完全正确,因为我们使用的是非凸优化算法)。这意味着,在这种方法中,具有不同解决方案的多次采样并不是很自然(但调用具有不同起点的 DCCP 例程将很可能导致不同的解决方案)。
备注 2: 对于非小型实例(使用非商业求解器)而言,性能非常糟糕,这使其成为一种更具理论性的方法。
实施
这里是一些使用 Python 3、numpy、scipy(最短路径)、cvxpy(选择问题公式)和 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 似乎仍在使用,这告诉我们,更好的求解器设置还有更多潜力。