一，旅行商问题QUBO的两种实现

提示：上篇已经讲过了旅行商问题的QUBO建模，这里直接讲两种编程实现：

看过上篇的读者应该已经注意到，因为旅行商问题需要最终返回到初始点的。所以，下面👇的目标函数里，循环进行到$N$时，最后一个$x_{j,t+1}$应该确定回到初始点的。

针对这个特殊设定，我们可以有两种实现方式：

• 方式一：使用取余操作符%，在$t=N$时，这样的话$(t+1)\%N$=1，也就相当于最后一个时间步回到了初始点。
  

• 方式二：把$x_{i,t}$ 和$x_{j,t+1}$对应的数值矩阵，分为独立的两个矩阵。具体来说，$x_{i,t}$从矩阵X中取值和$x_{j,t+1}$从矩阵Y中取值。矩阵Y的数值就是，$y_{j,t}=x_{j,t+1}$，从式子上看，有点难以理解，大家可以看下面的图示。其实就是把矩阵X的第一列，转移到最后一列。
  
  最后我们的目标函数就就看转换如下：

大家可以品味一下，∑符号怎么转换为矩阵操作。

提示：下面分别献上python实现：

二，方式一：取余操作

这里代码复制于下面的链接，这里只讲解QUBO部分的代码：

https://github.com/recruit-communications/pyqubo/blob/master/notebooks/TSP.ipynb
1

旅行问题的QUBO的定义：

%matplotlib inline
from pyqubo import Array, Placeholder, Constraint
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import neal

# 地点名和坐标 list[("name", (x, y))]
cities = [
    ("a", (0, 0)),
    ("b", (1, 3)),
    ("c", (3, 2)),
    ("d", (2, 1)),
    ("e", (0, 1))
]

n_city = len(cities)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

下面实现约束部分：

# ‘BINARY’代表目标变量时0或1 
x = Array.create('c', (n_city, n_city), 'BINARY')

# 约束① : 每个时间步只能访问1个地点
time_const = 0.0
for i in range(n_city):
    # Constraint(...)函数用来定义约束项
    time_const += Constraint((sum(x[i, j] for j in range(n_city)) - 1)**2, label="time{}".format(i))

# 约束② : 每个地点只能经过1次
city_const = 0.0
for j in range(n_city):
    city_const += Constraint((sum(x[i, j] for i in range(n_city)) - 1)**2, label="city{}".format(j))
1
2
3
4
5
6
7
8
9
10
11
12
13

接下来是目标函数：

# 目标函数的定义

distance = 0.0
for i in range(n_city):
    for j in range(n_city):
        for k in range(n_city):
            d_ij = dist(i, j, cities)
            distance += d_ij * x[k, i] * x[(k+1)%n_city, j]
1
2
3
4
5
6
7
8

最后就是整体的Hamiltonian量定义和输出采样结果。

# 总体的Hamiltonian，这里的A代表约束项的惩罚系数，数值自由指定
A = Placeholder("A")
H = distance + A * (time_const + city_const)

# 求BQM
model = H.compile()
feed_dict = {'A': 4.0}
bqm = model.to_bqm(feed_dict=feed_dict)


sa = neal.SimulatedAnnealingSampler()
# 这里要注意，退火算法是要采样足够的次数，从中取占比最高的结果作为最优解
sampleset = sa.sample(bqm, num_reads=100, num_sweeps=100)

# Decode solution
decoded_samples = model.decode_sampleset(sampleset, feed_dict=feed_dict)
best_sample = min(decoded_samples, key=lambda x: x.energy)


# 如果指定参数only_broken=True，则只会返回损坏的约束。
# 如果best_sample长度为0，就代表没有损坏的约束项，也就满足最优解了。
num_broken = len(best_sample.constraints(only_broken=True))
if num_broken == 0:
print(best_sample.sample)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

三、方式二：独立矩阵

这里的实现复制于以下文章：

https://qiita.com/yufuji25/items/0425567b800443a679f7
1

import neal
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
from pyqubo import Array, Placeholder
from scipy.spatial.distance import cdist


def construct_graph(n_pos:int):
    """ construct complete graph """
    pos = nx.spring_layout(nx.complete_graph(n_pos))
    coordinates = np.array(list(pos.values()))
    G = nx.Graph(cdist(coordinates, coordinates))
    nx.set_node_attributes(G, pos, "pos")
    return G

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

整体Hamiltonian量：

def create_hamiltonian(G:nx.Graph):
    """ create QUBO model from graph structure """

    # generate QUBO variable
    n_pos = G.number_of_nodes()
    X = np.array(Array.create("X", shape = (n_pos, n_pos), vartype = "BINARY"))

    # construct `Y` matrix
    Xtmp = np.concatenate([X, X[:, 0].reshape(-1, 1)], axis = 1)
    Y = np.delete(Xtmp, 0, 1)

    # Distance matrix
    L = np.array(nx.adjacency_matrix(G).todense())

    # return hamiltonian
    Hbind1 = np.sum((1 - X.sum(axis = 0)) ** 2)
    Hbind2 = np.sum((1 - X.sum(axis = 1)) ** 2)
    Hobj = np.sum(X.dot(Y.T) * L)
    H = Hobj + Placeholder("a1") * Hbind1 + Placeholder("a2") * Hbind2

    return H.compile()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

求解并输出结果：

def decode(response, Gorigin:nx.Graph):
    """ return output graph generated from response """

    # derive circuit
    solution = min(response.record, key = lambda x: x[1])[0]
    X = solution.reshape((num_pos, num_pos))
    circuit = np.argmax(X, axis = 0).tolist()
    circuit.append(circuit[0])

    # generate output graph structure
    Gout = nx.Graph()
    Gout.add_edges_from(list(zip(circuit, circuit[1:])))
    positions = nx.get_node_attributes(Gorigin, "pos")
    nx.set_node_attributes(Gout, positions, "pos")

    return Gout


def draw_graph(G:nx.Graph, **config):
    """ drawing graph """
    plt.figure(figsize = (5, 4.5))
    nx.draw(G, nx.get_node_attributes(G, "pos"), **config)
    plt.show()



# configuration for drawing graph
config = {"with_labels":True, "node_size":500, "edge_color":"r", "width":1.5,
          "node_color":"k", "font_color":"w", "font_weight":"bold"}

# construct complete graph
num_pos = 8
G = construct_graph(num_pos)
draw_graph(G, **config)

# sampling
model = create_hamiltonian(G)
qubo, offset = model.to_qubo(feed_dict = {"a1":500, "a2":500})
response = neal.SimulatedAnnealingSampler().sample_qubo(qubo, num_reads = 1000)


# output graph
Gout = decode(response, G)
draw_graph(Gout, **config)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

总结

以上就是两种实现方式，大家可以体会一下，怎么实现稍微复杂的Hamiltonian量，接下来，讲解一下量子退火的物理原理。

参考文献：
[*1] : https://www.nttdata.com/jp/ja/-/media/nttdatajapan/files/news/services_info/2021/012800/012800-01.pdf
[*2] : https://qiita.com/yufuji25/items/0425567b800443a679f7