本篇通过两个贴近日常生活的例子——整数分割问题和旅行商问题，介绍了如何将经典组合优化问题转化为 QUBO形式，并使用 Python 和模拟量子退火算法进行求解。内容涵盖建模思路、目标函数构建、约束定义及代码实现，帮助读者直观理解量子优化的基本原理与实际应用流程。

在量子退火算法（1）中我们学习了什么是QUBO，那么在实际生活中，哪些问题可以用量子计算来帮助解决？本篇从把几个数分成两组、找出最短旅行路线这些小例子开始，带你了解如何把这些问题转化为量子可处理的 QUBO 形式，并用 Python 实现。别担心，这里没有晦涩难懂的公式推导，有的是贴近生活的例子、清晰的逻辑和可运行的代码，帮你轻松打开量子优化世界的大门。

目录

一、整数分割问题
      1.1 什么是整数分割问题
      1.2 转化为组合优化问题
      1.3 目标函数转化为QUBO
      1.4 PyQUBO实现Ising模型
二、旅行商问题（Traveling Saleman Problem，TSP）
      2.1 什么是旅行商问题
           2.1.1 旅行商问题的定义
           2.1.2 旅行商问题求解的计算量
      2.2 TSP问题的建模
           2.2.1 总体Hamilton量H
           2.2.2 约束条件
           2.2.3 目标函数
      2.3 旅行商问题QUBO的两种实现
      2.4 方式一：取余操作
      2.5 方式二：独立矩阵
三、小结

一、整数分割问题

1.1 什么是整数分割问题

QUBO建模最重要的就是，把建模对象中的变量映射为binary（0/1 或者 -1/+1）的变量。我先从简单的问题开始说明，让大家有些直观感受。整数分割问题就是一个非常简单，并容易理解的例子。此文参考了日本NTT公司的量子计算指南文档[*1]。

整数分割问题定义：
判断能否将一个N 个整数 a₁ , ・・・ a_N 的整数集合分割成两个子集合，并且这两个子集合里的元素之和相等。

例子：

我们可以看到，上面👆的例子，分割后的两个子集合A和B的元素之和都等于6，所以该集合是可以满足整数分割的。

1.2 转化为组合优化问题

之前讲的QUBO的例子里，用的变量都是0或1。其实，还可以是-1或+1，这时候也叫ising模型。求解变量什么时候用0/1，什么用-1/+1，这个之后用例子给大家解释。

这次的问题，是把求解的两个子集合的标签作为-1/+1的变量。如下图所示。

我们的优化目标就成为了，最小化两个子集合的差值的平方和。大家可以思考一下，如果使用0/1作为子集合A和B的标签，我们要怎么定义最小化的目标函数。

目标集合为{ 1, 5, 6 }时的目标函数就是：

我们来枚举一下，所有 x i 从-1或+1取值的组合结果：

我们可以看到，目标函数值为0时，我们得到了正确的整数分割结果。

1.3 目标函数转化为QUBO

因为已经获得了目标函数，那我们先把多项式展开就好了。展开结果如下图所示：

因为【 x_i 从-1或+1取值】这个设定，下面👇的式子是成立的：

所以我们可以得到下面的QUBO结果：

最后献上Python代码。

1.4 PyQUBO实现Ising模型

之前的目标函数都是展开后的二次多项式，大家可以直接计算出QUBO矩阵。这次使用PyQUBO直接定义目标函数，大家就不用手动求解QUBO矩阵了。

import pyqubo

import neal

x = pyqubo.Array.create('x', shape=(3), vartype='SPIN') # 'SPIN' 就表示目标变量是从{-1, 1}取值。目标变量需要从{0, 1}中取值时，就设定为 'BINARY' 



objective_function = (1 * x[0] + 5 * x[1] + 6 * x[2]) ** 2



model = objective_function.compile() 

bqm = model.to_bqm()



print("我们可以将bqm转为ising或qubo输出")

print(bqm.to_ising())



sa = neal.SimulatedAnnealingSampler()

sampleset = sa.sample(bqm, num_reads=10)

samples = model.decode_sampleset(sampleset)

best_sample = min(samples, key=lambda s: s.energy)



print("求解时，pyqubo内部已经将ising模型转换为qubo的0或1，所以输出结果为0或1")

print(best_sample.sample) 

运行结果如下：

我们可以将bqm转为ising或qubo输出

({'x[2]': 0.0, 'x[0]': 0.0, 'x[1]': 0.0}, {('x[0]', 'x[2]'): 12.0, ('x[1]', 'x[2]'): 60.0, ('x[1]', 'x[0]'): 10.0}, 62.0)

求解时，pyqubo内部已经将ising模型转换为qubo的0或1，所以输出结果为0或1

{'x[2]': 1, 'x[0]': 0, 'x[1]': 0}

以上就是一个简单建模的例子。下面讲旅行商问题的建模。

二、旅行商问题（Traveling Saleman Problem，TSP）

2.1 什么是旅行商问题

2.1.1 旅行商问题的定义

旅行商问题，是一个经典的组合优化问题，而且是著名NP问题之一。如下图所示，可以想象，有A，B，C，D，E 五个地点，我们想找到一条路径，从地点A出发，经过剩余四个地点，然后回到地点A，从所有可能路径中找到距离最短的一条路径。本章借用了文献[1]的图表。

2.1.2 旅行商问题求解的计算量

最简单的求解方式就是，如下图所示把所有的求解路径全部计算一遍，然后算出每条路径的长度，求出最短路径。

如下图所示，所有的枚举路径总共有24条，我们可以很快找到最短路径。

如果下面A～Z的情况，这个计算量，日本的第一超级计算机富岳，每秒的计算速度约为44.2京次（京是10的16次方，即万兆）。一年的秒数是3600×24×365=3153.6万秒。有兴趣的可以计算一下要算多少年。

2.2 TSP问题的建模

2.2.1 总体Hamilton量H

该问题输入有两个，这里借用了文章[*2]的图表：

 · 地点数目：N 

 · 地点之间的距离：l _{i , j} ​(i=1,・・・,N)

约束条件：

 · 每个时间步只能访问一个地点。

 · 每个地点都访问过一次。

整体的Hamilton量H如下：

目标变量x _{i , j} 的两个下标的意思如下图👇所示，绿色的圆圈代表在某个时间步访问了某个第地点，所以我们的目标变量就可以用0或1表示了，0代表未访问，1代表访问。

2.2.2 约束条件

约束条件比较简单，先从约束条件解释，这里有2个约束可以解释如下：

（1）每个时间步只能访问一个地点。

         => 上图矩阵里的每列元素之和必须为1。也就是每列中只有一个元素为1。

（2）每个地点都访问过一次。

         => 上图矩阵里的每行元素之和必须为1。也就是每行中只有一个元素为1。

具体表达式如下：

2.2.3 目标函数

解析：

 · x _{i , t} x _{j , t + 1} ：这里的目标函数，最难理解的是x _{i , t} x _{j , t + 1} 。可以理解为【t时间步访问地点 i ，t+1 时间步访问地点 j 时，x _{i , t} x _{j , t + 1} =1；其他的情况，x _{i , t} x _{j , t + 1} =0】。

 · ∑ i= 1^ N ∑ j = 1^N ∑ t = 1^N 

该表达式代表了，【t时间步访问地点 i，t+1 时间步访问地点 j 时，地点 i 和 j 之间的距离 L _{i , j} 之和】。所以，这个目标函数就代表了，从初始地点，经过所有地点后，回到初始地点的距离总和。

2.3 旅行商问题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实现。

2.4 方式一：取余操作

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

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

旅行问题的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)

下面实现约束部分：

# ‘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))

接下来是目标函数：

# 目标函数的定义



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]

最后就是整体的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)

2.5 方式二：独立矩阵

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

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

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

整体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()

求解并输出结果：

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)

以上就是两种实现方式，大家可以体会一下，怎么实现稍微复杂的Hamiltonian量。

三、小结

本篇基于整数分割问题和旅行商问题，讲解了量子退火算法的QUBO设计和Python实现。下篇我们将在此基础上，介绍如何用这种方法解决车辆路径问题和护士调度问题：量子退火算法（3）。

在阅读参考文献时，经常会发现资料里的一些小错误，大家以后阅读资料时也要小心啊。

参考文献：

[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

 

---

本文转载自CSDN博主：gang_unerry

版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。

原文链接：

https://blog.csdn.net/gangshen1993/article/details/127594967

https://blog.csdn.net/gangshen1993/article/details/127602524

https://blog.csdn.net/gangshen1993/article/details/127638370

本篇通过两个贴近日常生活的例子——整数分割问题和旅行商问题，介绍了如何将经典组合优化问题转化为 QUBO形式，并使用 Python 和模拟量子退火算法进行求解。内容涵盖建模思路、目标函数构建、约束定义及代码实现，帮助读者直观理解量子优化的基本原理与实际应用流程。

在量子退火算法（1）中我们学习了什么是QUBO，那么在实际生活中，哪些问题可以用量子计算来帮助解决？本篇从把几个数分成两组、找出最短旅行路线这些小例子开始，带你了解如何把这些问题转化为量子可处理的 QUBO 形式，并用 Python 实现。别担心，这里没有晦涩难懂的公式推导，有的是贴近生活的例子、清晰的逻辑和可运行的代码，帮你轻松打开量子优化世界的大门。

目录

一、整数分割问题
      1.1 什么是整数分割问题
      1.2 转化为组合优化问题
      1.3 目标函数转化为QUBO
      1.4 PyQUBO实现Ising模型
二、旅行商问题（Traveling Saleman Problem，TSP）
      2.1 什么是旅行商问题
           2.1.1 旅行商问题的定义
           2.1.2 旅行商问题求解的计算量
      2.2 TSP问题的建模
           2.2.1 总体Hamilton量H
           2.2.2 约束条件
           2.2.3 目标函数
      2.3 旅行商问题QUBO的两种实现
      2.4 方式一：取余操作
      2.5 方式二：独立矩阵
三、小结

一、整数分割问题

1.1 什么是整数分割问题

QUBO建模最重要的就是，把建模对象中的变量映射为binary（0/1 或者 -1/+1）的变量。我先从简单的问题开始说明，让大家有些直观感受。整数分割问题就是一个非常简单，并容易理解的例子。此文参考了日本NTT公司的量子计算指南文档[*1]。

整数分割问题定义：
判断能否将一个N 个整数 a₁ , ・・・ a_N 的整数集合分割成两个子集合，并且这两个子集合里的元素之和相等。

例子：

我们可以看到，上面👆的例子，分割后的两个子集合A和B的元素之和都等于6，所以该集合是可以满足整数分割的。

1.2 转化为组合优化问题

之前讲的QUBO的例子里，用的变量都是0或1。其实，还可以是-1或+1，这时候也叫ising模型。求解变量什么时候用0/1，什么用-1/+1，这个之后用例子给大家解释。

这次的问题，是把求解的两个子集合的标签作为-1/+1的变量。如下图所示。

我们的优化目标就成为了，最小化两个子集合的差值的平方和。大家可以思考一下，如果使用0/1作为子集合A和B的标签，我们要怎么定义最小化的目标函数。

目标集合为{ 1, 5, 6 }时的目标函数就是：

我们来枚举一下，所有 x i 从-1或+1取值的组合结果：

我们可以看到，目标函数值为0时，我们得到了正确的整数分割结果。

1.3 目标函数转化为QUBO

因为已经获得了目标函数，那我们先把多项式展开就好了。展开结果如下图所示：

因为【 x_i 从-1或+1取值】这个设定，下面👇的式子是成立的：

所以我们可以得到下面的QUBO结果：

最后献上Python代码。

1.4 PyQUBO实现Ising模型

之前的目标函数都是展开后的二次多项式，大家可以直接计算出QUBO矩阵。这次使用PyQUBO直接定义目标函数，大家就不用手动求解QUBO矩阵了。

import pyqubo

import neal

x = pyqubo.Array.create('x', shape=(3), vartype='SPIN') # 'SPIN' 就表示目标变量是从{-1, 1}取值。目标变量需要从{0, 1}中取值时，就设定为 'BINARY' 



objective_function = (1 * x[0] + 5 * x[1] + 6 * x[2]) ** 2



model = objective_function.compile() 

bqm = model.to_bqm()



print("我们可以将bqm转为ising或qubo输出")

print(bqm.to_ising())



sa = neal.SimulatedAnnealingSampler()

sampleset = sa.sample(bqm, num_reads=10)

samples = model.decode_sampleset(sampleset)

best_sample = min(samples, key=lambda s: s.energy)



print("求解时，pyqubo内部已经将ising模型转换为qubo的0或1，所以输出结果为0或1")

print(best_sample.sample) 

运行结果如下：

我们可以将bqm转为ising或qubo输出

({'x[2]': 0.0, 'x[0]': 0.0, 'x[1]': 0.0}, {('x[0]', 'x[2]'): 12.0, ('x[1]', 'x[2]'): 60.0, ('x[1]', 'x[0]'): 10.0}, 62.0)

求解时，pyqubo内部已经将ising模型转换为qubo的0或1，所以输出结果为0或1

{'x[2]': 1, 'x[0]': 0, 'x[1]': 0}

以上就是一个简单建模的例子。下面讲旅行商问题的建模。

二、旅行商问题（Traveling Saleman Problem，TSP）

2.1 什么是旅行商问题

2.1.1 旅行商问题的定义

旅行商问题，是一个经典的组合优化问题，而且是著名NP问题之一。如下图所示，可以想象，有A，B，C，D，E 五个地点，我们想找到一条路径，从地点A出发，经过剩余四个地点，然后回到地点A，从所有可能路径中找到距离最短的一条路径。本章借用了文献[1]的图表。

2.1.2 旅行商问题求解的计算量

最简单的求解方式就是，如下图所示把所有的求解路径全部计算一遍，然后算出每条路径的长度，求出最短路径。

如下图所示，所有的枚举路径总共有24条，我们可以很快找到最短路径。

如果下面A～Z的情况，这个计算量，日本的第一超级计算机富岳，每秒的计算速度约为44.2京次（京是10的16次方，即万兆）。一年的秒数是3600×24×365=3153.6万秒。有兴趣的可以计算一下要算多少年。

2.2 TSP问题的建模

2.2.1 总体Hamilton量H

该问题输入有两个，这里借用了文章[*2]的图表：

 · 地点数目：N 

 · 地点之间的距离：l _{i , j} ​(i=1,・・・,N)

约束条件：

 · 每个时间步只能访问一个地点。

 · 每个地点都访问过一次。

整体的Hamilton量H如下：

目标变量x _{i , j} 的两个下标的意思如下图👇所示，绿色的圆圈代表在某个时间步访问了某个第地点，所以我们的目标变量就可以用0或1表示了，0代表未访问，1代表访问。

2.2.2 约束条件

约束条件比较简单，先从约束条件解释，这里有2个约束可以解释如下：

（1）每个时间步只能访问一个地点。

         => 上图矩阵里的每列元素之和必须为1。也就是每列中只有一个元素为1。

（2）每个地点都访问过一次。

         => 上图矩阵里的每行元素之和必须为1。也就是每行中只有一个元素为1。

具体表达式如下：

2.2.3 目标函数

解析：

 · x _{i , t} x _{j , t + 1} ：这里的目标函数，最难理解的是x _{i , t} x _{j , t + 1} 。可以理解为【t时间步访问地点 i ，t+1 时间步访问地点 j 时，x _{i , t} x _{j , t + 1} =1；其他的情况，x _{i , t} x _{j , t + 1} =0】。

 · ∑ i= 1^ N ∑ j = 1^N ∑ t = 1^N 

该表达式代表了，【t时间步访问地点 i，t+1 时间步访问地点 j 时，地点 i 和 j 之间的距离 L _{i , j} 之和】。所以，这个目标函数就代表了，从初始地点，经过所有地点后，回到初始地点的距离总和。

2.3 旅行商问题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实现。

2.4 方式一：取余操作

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

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

旅行问题的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)

下面实现约束部分：

# ‘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))

接下来是目标函数：

# 目标函数的定义



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]

最后就是整体的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)

2.5 方式二：独立矩阵

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

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

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

整体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()

求解并输出结果：

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)

以上就是两种实现方式，大家可以体会一下，怎么实现稍微复杂的Hamiltonian量。

三、小结

本篇基于整数分割问题和旅行商问题，讲解了量子退火算法的QUBO设计和Python实现。下篇我们将在此基础上，介绍如何用这种方法解决车辆路径问题和护士调度问题：量子退火算法（3）。

在阅读参考文献时，经常会发现资料里的一些小错误，大家以后阅读资料时也要小心啊。

参考文献：

[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

 

---

本文转载自CSDN博主：gang_unerry

版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。

原文链接：

https://blog.csdn.net/gangshen1993/article/details/127594967

https://blog.csdn.net/gangshen1993/article/details/127602524

https://blog.csdn.net/gangshen1993/article/details/127638370