文章目录

• 一，旅行商问题QUBO的两种实现
• 二，方式一：取余操作
• 三、方式二：独立矩阵
• 总结

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

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

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

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

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

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

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

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

二，方式一：取余操作

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

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

旅行问题的QUBO的定义：

%matplotlib inlinefrom pyqubo import Array, Placeholder, Constraintimport matplotlib.pyplot as pltimport networkx as nximport numpy as npimport 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)

下面实现约束部分：

# ‘BINARY’代表目标变量时0或1 x = Array.create('c', (n_city, n_city), 'BINARY')# 约束① : 每个时间步只能访问1个地点time_const = 0.0for 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.0for j in range(n_city):city_const += Constraint((sum(x[i, j] for i in range(n_city)) - 1)**2, label="city{}".format(j))

接下来是目标函数：

# 目标函数的定义distance = 0.0for 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]

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

# 总体的Hamiltonian，这里的A代表约束项的惩罚系数，数值自由指定A = Placeholder("A")H = distance + A * (time_const + city_const)# 求BQMmodel = 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 solutiondecoded_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)

三、方式二：独立矩阵

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

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

import nealimport numpy as npimport networkx as nximport matplotlib.pyplot as pltfrom pyqubo import Array, Placeholderfrom scipy.spatial.distance import cdistdef 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

整体Hamiltonian量：

def create_hamiltonian(G:nx.Graph):""" create QUBO model from graph structure """# generate QUBO variablen_pos = G.number_of_nodes()X = np.array(Array.create("X", shape = (n_pos, n_pos), vartype = "BINARY"))# construct `Y` matrixXtmp = np.concatenate([X, X[:, 0].reshape(-1, 1)], axis = 1)Y = np.delete(Xtmp, 0, 1)# Distance matrixL = np.array(nx.adjacency_matrix(G).todense())# return hamiltonianHbind1 = 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") * Hbind2return H.compile()

求解并输出结果：

def decode(response, Gorigin:nx.Graph):""" return output graph generated from response """# derive circuitsolution = 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 structureGout = 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 Goutdef 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 graphconfig = {"with_labels":True, "node_size":500, "edge_color":"r", "width":1.5,"node_color":"k", "font_color":"w", "font_weight":"bold"}# construct complete graphnum_pos = 8G = construct_graph(num_pos)draw_graph(G, **config)# samplingmodel = create_hamiltonian(G)qubo, offset = model.to_qubo(feed_dict = {"a1":500, "a2":500})response = neal.SimulatedAnnealingSampler().sample_qubo(qubo, num_reads = 1000)# output graphGout = decode(response, G)draw_graph(Gout, **config)

总结

以上就是两种实现方式，大家可以体会一下，怎么实现稍微复杂的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