案例解析-旅行商问题

什么是旅行商问题

旅行商问题(TravelingSalesmanProblem，TSP)是一个经典的组合优化问题。经典的TSP可以描述为：一个商品推销员要去若干个城市推销商品，该推销员从一个城市出发，需要经过所有城市后，回到出发地。应如何选择行进路线，以使总的行程最短。

从图论的角度来看，该问题实质是在一个带权完全无向图中，找一个权值最小的Hamilton回路。由于该问题的可行解是所有顶点的全排列，随着顶点数的增加，会产生组合爆炸，它是一个NP完全问题。

由于其在交通运输、电路板线路设计以及物流配送等领域内有着广泛的应用，国内外学者对其进行了大量的研究。早期的研究者使用精确算法求解该问题，常用的方法包括:分支定界法、线性规划法、动态规划法等。

但是，随着问题规模的增大，精确算法将变得无能为力，因此，在后来的研究中，国内外学者重点使用近似算法或启发式算法，主要有遗传算法、模拟退火法、蚁群算法、禁忌搜索算法、贪婪算法和神经网络等。如今，量子计算的出现为该问题的求解提供了全新的方式。

旅行商问题分析

1.  问题本质

旅行商问题有一点像“最短路径问题”，然后我们就会自然地想到用Dikstra算法去求得“从某一个城市出发，到其他所有剩余城市的最短路径”，再或者如果是个真实地图，我们可以用启发式的“A星算法”快速搜索出“从某-个城市到另一个指定城市间的最短路径"。但仔细想，这个问题并非单纯这么简单，它还要求去寻找“从某个城市开始，分别经过其它城市一次且仅一次，最后再回到这个出发城市的最短巡回路径”。

2.  深入分析

旅行商问题要从图G的所有周游路线中求取最小成本的周游路线，而从初始点出发的周游路线一共有(n-1)!条，即等于除初始结点外的n-1个结点的排列数，因此旅行商问题是一个排列问题。排列问题比子集合的选择问题通常要难于求解得多，这是因为n个物体有n!种排列，只有 个子集合(n!>O()）。通过枚举(n-1)!条周游路线，从中找出一条具有最小成本的周游路线的算法，其计算时间显然为O(n!)。

所以该怎么求解呢，我们很容易想到一种类似于穷举的思路：

现在假设我们要拜访11个城市，从城市1出发，最后回到城市1。显然，从城市1出来后，我们随即可以选择剩余的10个城市之一进行拜访(这里所有城市都是连通的，总是可达的，而不连通的情况属于个人特殊业务的装饰处理，不是本案例考虑范畴)，那么很显然这里就有10种选择，以此类推，下一次就有9种选择.总的可选路线数就是：10!。也就是说需要用for循环迭代10!次，才能找出所有的路线，进而筛选出最短的那条路线。如果只拜访10个城市的话(需要迭代3628800次)或许还好，那要拜访100个城市(需要迭代9.3326215443944*10^157)简直就是经典计算机的噩梦！更多个城市的话，计算的时间开销可想而知！

更一般地，如果要拜访n+1个城市，总的可选路线数就是n!，进而时间复杂度就是O(n!)，从这里我们同理也可以看出，这个算法的时间复杂度是非多项式的，它的开销大是显而易见的。所以问题的关键不在于寻找两两城市间的最短路径，而在于去寻找一那条最短的巡回路径，换句话说，就是寻找一组拜访城市的先后次序序列。

旅行商问题求解方案

传统的旅行商求解方案主要有:

●  TSP旅行商问题-蛮力法( 深度遍历优先算法DFS)

●  TSP_旅行商问题-动态规划

●  TSP旅行商问题-模拟退火算法

●  TSP 旅行商问题-遗传算法

●  TSP 旅行商问题-粒子群算法

●  TSP 旅行商问题-神经网络

旅行商问题是个NP完全问题，穷举算法的效率又不高，那我们该如何通过一个多项式时间复杂度的算法快速求出这个先后次序呢？目前比较主流的方法是采用一些随机的、启发式的搜索算法，比如遗传算法、蚁群算法、模拟退货算法、粒子群算法等。但这些算法都有一个缺点，就是不一定能求出最优解，只能收敛于(近似逼近)最优解，得到一个次优解，因为他们本质都是随机算法，大多都会以类似“一定概率接受或舍去"的思路去筛选解。各算法的实现思路都有不同，但也或多或少有互相借鉴的地方，有的与随机因子有关、有的与初始状态有关、有的与随机函数有关、有的与选择策略有关。

综合上述分析，TSP问题的求解大概是由以下两步构成：

1.  计算两两城市间的最短路径：利用类似Diikstra、Flord、A星的算法求出最短路线。

2.  计算最短巡回路径：利用类似遗传算法、蚁群算法的搜索算法求巡回拜访的次序。

关于1中需要说明一点，就是现实生活中我们的地图往往不是一个完全图，而是一个非完全图，甚至有些节点仅仅是道路的分岔口，而不是城市节点。完全图和非完全图的区别如下所示。

完全图和非完全图

完全图：两两城市间都有直达的路线，这条路线不需要经过中间其他节点；

非完全图：偶尔有两个城市间的路线需要经过其他中间节点。

应用量子计算求解旅行商问题

1.  变量定义

给定一张有权重的连接图，V是节点的集合，E是边的集合。

我们使用距离矩阵来表示图 ，当两点间不存在边时，其权重值为 。

令 ，N表示节点的数量。

令代表点在规划路径中第个点被访问，否则。 其中, 。

2.  约束条件处理

首先，我们需要对于每一个，只对应一个顺序，则有：

其次，我们需要对于每一个顺序，只对应一个点，则有：

将以上两个约束写成QUBO形式，有：

对于，我们希望、不连续出现在位置、上，则有：

3.  QUBO模型构建

如此，由以上的约束条件，我们可以构建一个遍历所有的点的路径的模型（又称哈密尔顿环）：

然而对于旅行商问题，在此基础上，我们想要获得上述路径中最短的一个。令为边的权重，则有：

将两个QUBO表达式相加，即可得到最终的优化模型：

优化目标是求 H 的最小值。

4.  使用开物SDK进行Ising矩阵生成

加载依赖包

# Import numpy and kaiwu
import numpy as np
import kaiwu as kw

加载图数据，并做预处理

# Import distance matrix
w = np.array([[ 0, 13, 11, 16,  8],
              [13,  0,  7, 14,  9],
              [11,  7,  0, 10,  9],
              [16, 14, 10,  0, 12],
              [ 8,  9,  9, 12,  0]])

# Get the number of nodes
n = w.shape[0]

# Create qubo variable matrix
x = kw.qubo.ndarray((n, n), "x", kw.qubo.Binary)

# Get sets of edge and non-edge pairs
edges = [(u, v) for u in range(n) for v in range(n) if w[u, v] != 0]
no_edges = [(u, v) for u in range(n) for v in range(n) if w[u, v] == 0]

加载约束项，构建目标函数

# Node constraint: Each node must belong to exactly one position
sequence_cons = kw.qubo.quicksum([(1 - kw.qubo.quicksum([x[v, j] for j in range(n)])) ** 2 for v in range(n)])

# Position constraint: Each position can have only one node
node_cons = kw.qubo.quicksum([(1 - kw.qubo.quicksum([x[v, j] for v in range(n)])) ** 2 for j in range(n)])

# Edge constraint: Pairs without edges cannot appear in the path
connect_cons = kw.qubo.quicksum([kw.qubo.quicksum([x[u, j] * x[v, j + 1] for j in range(n - 1)]) + x[u, n - 1] * x[v, 0] for u, v in no_edges])

# Hamiltonian cycle constraint: Sum of the above three constraints
ham_cycle = sequence_cons + node_cons + connect_cons

# TSP path cost
path_cost = kw.qubo.quicksum([w[u, v] * (kw.qubo.quicksum([x[u, j] * x[v, j + 1] for j in range(n - 1)]) + x[u, n - 1] * x[v, 0]) for u, v in edges])

# Final objective function with penalty factor 100 for the Hamiltonian constraint
obj = 100 * ham_cycle + path_cost

模型编译，生成Ising矩阵

# Parse QUBO
obj = kw.qubo.make(obj)

# Convert to Ising model
obj_ising = kw.qubo.qubo_model_to_ising_model(obj)

# Extract the Ising matrix
matrix = obj_ising.get_ising()["ising"]

5.  使用光量子计算模拟器进行计算

该问题有三个硬约束（点约束、位置约束和边约束）。硬约束的取值必须严格为才能形成合法的路径。而使用光量子计算模拟器的结果有时硬约束并不为，因此我们通过在光量子计算模拟器的结果上运行一遍禁忌搜索在该结果的邻域中找到满足硬约束的解。

# Perform calculation using CIM simulator
worker = kw.cim.SimulatedCIMOptimizer(
    pump=1.3,
    noise=0.2,
    laps=5000,
    delta_time=0.05,
    normalization=0.3,
    iterations=50
)
output = worker.solve(matrix)

# Sort the results
opt = kw.sampler.optimal_sampler(matrix, output, bias=0, negtail_ff=False)

# Select the best solution
cim_best = opt[0][0]

# If the linear term variable is -1, perform a flip
cim_best = cim_best * cim_best[-1]

print(cim_best)

通过模拟器求解，得到Ising解向量

[ 1 -1 -1 -1 -1 -1  1 -1 -1 -1 -1 -1  1 -1 -1 -1 -1 -1  1 -1 -1 -1 -1 -1  1  1]

6.  还原Ising解向量的spin结果为QUBO变量值

以上我们得到了spin结果，需要将其转化为原始变量的值，才能进一步给出求得的路径。我们使用get_val模块进行转化。

# Get the list of variable names
vars = obj_ising.get_variables()

# Substitute the spin vector and obtain the result dictionary
sol_dict = kw.qubo.get_sol_dict(cim_best, vars)

# Check the hard constraints for validity and path length
seq_val = kw.qubo.get_val(sequence_cons, sol_dict)
node_val = kw.qubo.get_val(node_cons, sol_dict)
ham_val = kw.qubo.get_val(ham_cycle, sol_dict)
print('position cons: {}'.format(seq_val))
print('node_cons cons: {}'.format(node_val))
print('ham_cycle: {}'.format(ham_val))

# Calculate the path length using path_cost
path_val = kw.qubo.get_val(path_cost, sol_dict)
print('path_cost: {}'.format(path_val))

position cons: 0.0
node_cons cons: 0.0
ham_cycle: 0.0
path_cost: 50.0

以上我们获得了路径的合法性和长度，下面我们恢复变量的值并获取路径。

if ham_val == 0:
    print('valid path')

    # Get the numerical value matrix of x
    x_val = kw.qubo.get_array_val(x, sol_dict)

    # Find the indices of non-zero items
    nonzero_index = np.array(np.nonzero(x_val)).T

    # Sort the indices of non-zero items by the second column
    nonzero_index = nonzero_index[nonzero_index[:, 1].argsort()]

    # Print the path order
    print(nonzero_index[:,0].flatten())
else:
    print('invalid path')

valid path
[0 1 2 3 4]

——end——