本篇围绕两类经典组合优化问题——车辆路径问题（VRP）和护士调度问题（NSP），探讨了其在量子计算中的 QUBO 建模与 Python 实现过程，并使用量子退火算法进行求解与验证。同时在实验的基础上，指出量子退火算法的缺陷。

在量子退火算法（2）中，我们学习了整数分割问题和旅行商问题的QUBO 建模与 Python 实现过程。本篇将在系列前文的基础上，探讨车辆路径问题和护士调度问题，具体内容如目录所示。

目录

一、车辆路径问题 (VRP)
      1.1 什么是车辆路径问题
      1.2 车辆路径问题（VRP）的QUBO建模
           1.2.1 目标变量
           1.2.2 目标函数与约束条件
      1.3 Python实现VRP的QUBO
           1.3.1 引入库
           1.3.2 设置参数和距离矩阵
           1.3.2 设置参数和距离矩阵
           1.3.3 QUBO实现
           1.3.4 OpenJij求解QUBO目标变量
           1.3.5 输出求解结果路径
      1.3 小结
二、护士调度问题(NSP)
      2.1 什么是护士调度问题
      2.2 护士调度问题（NSP）的QUBO建模
           2.2.1 目标变量
           2.2.2 约束条件定义
           2.2.3【约束a】的补充说明
           2.2.4【约束c】的补充说明
      2.3 Python实现NSP的QUBO
           2.3.1 参数设定表
           2.3.2 QUBO实现
           2.3.3 实验结果检验和可视化
      2.4 实验结果
三、总结

一、车辆路径问题 (VRP)

1.1 什么是车辆路径问题

QUBO为解决组合优化问题的利器，车辆路径问题是最经常被提起的现实应用。车辆路径问题 (VRP) 是优化多辆送货车辆的送货顺序的组合优化问题，是旅行商问题 (TSP) 的推广。这是一个可以应用于各种现实世界问题的问题，例如物流中的货物配送、工厂中的货物搬运，以及按需运输服务中的车辆分配计划。VRP 是一个典型的组合优化问题，也是经典的NP-Hard问题，普通计算机很难在多项式时间内找到精确解。量子退火有望利用量子隧穿效应实时解决 NP-Hard 组合优化问题。

1.2 车辆路径问题（VRP）的QUBO建模

作为一个基本的VRP，我们考虑了当多辆运送车辆访问所有目标基地时，搜索总旅行成本最小的路线的组合优化问题。我们假设所有的送货车辆都从一个叫做 Depot 的基地出发，并且总是在最后返回到 Depot。移动成本使用基地之间的距离计算。 VRP的QUBO的定义为：把地点1作为Depot，然后假设有2辆车用来配送4个地点。2辆车的配送示意图如下：

其实大家可以可以把VRP问题理解为，有V辆车，V辆车各自负责不重叠的地点子集，所有车辆负责的地点子集的路径总和最小。

所以这里我们有三个参数需要设定：

• V：车辆个数


• P：地点个数


• S：时间步（TSP建模中，我们不需要设置时间步，因为S=P+1）

小写字母代表变量：

1.2.1 目标变量

我们期待的QUBO矩阵如下图所示，

• 红色矩阵代表车辆1，从地点1==>地点3==>地点1；


• 蓝色矩阵代表车辆2，从地点1==>地点2==>地点4==>地点1。

1.2.2 目标函数与约束条件

目标函数：其实就是TSP的目标函数，又增加了一个车辆的维度

约束条件：

条件a：第1个时间步的地点必须是Depot地点
条件b：第S个时间步的地点必须是Depot地点
条件c：除了Depot地点，剩余地点必须被访问过一次
条件d：一个时间步，车辆v不能同时出现在多个地点p

各个约束条件对应的表达式如下图：

最后转换后的哈密顿量为：

1.3 Python实现VRP的QUBO

本章节主要来自以下代码库：

https://github.com/yskmt2018/quantum2020/blob/master/vehicle_routing_problem.ipynb

1.3.1 引入库

import math

import pyqubo

import numpy as np

import pandas as pd

from openjij import SQASampler

1.3.2 设置参数和距离矩阵

# 车辆个数

V = 3

# 地点个数

L = len(distance[0])

# 最大时间步

S = 10



# 距离矩阵

distance = np.array([

  [0, 548, 776, 696, 582, 274, 502, 194, 308, 194, 536, 502, 452, 1020, 393, 208, 1196, 983, 427, 482, 479, 517],

  [548, 0, 684, 308, 194, 502, 730, 354, 696, 742, 1084, 594, 1102, 298, 741, 582, 712, 1085, 1024, 202, 220, 1050],

  [776, 684, 0, 992, 878, 502, 274, 810, 468, 742, 400, 1278, 1152, 819, 510, 514, 399, 537, 579, 509, 285, 555],

  [696, 308, 992, 0, 114, 650, 878, 502, 844, 890, 1232, 514, 300, 263, 593, 153, 801, 1128, 729, 1230, 317, 946],

  [582, 194, 878, 114, 0, 536, 764, 388, 730, 776, 1118, 400, 511, 886, 168, 514, 396, 594, 886, 190, 729, 1042],

  [274, 502, 502, 650, 536, 0, 228, 308, 194, 240, 582, 776, 270, 343, 716, 535, 1269, 410, 1177, 494, 1133, 651],

  [502, 730, 274, 878, 764, 228, 0, 536, 194, 468, 354, 1004, 144, 353, 408, 220, 636, 514, 1180, 1146, 130, 648],

  [194, 354, 810, 502, 388, 308, 536, 0, 342, 388, 730, 468, 476, 1296, 338, 1100, 1125, 986, 710, 320, 279, 460],

  [308, 696, 468, 844, 730, 194, 194, 342, 0, 274, 388, 810, 263, 138, 145, 740, 1251, 1262, 402, 569, 189, 877],

  [194, 742, 742, 890, 776, 240, 468, 388, 274, 0, 342, 536, 1190, 757, 293, 580, 1195, 124, 268, 505, 309, 602],

  [536, 1084, 400, 1232, 1118, 582, 354, 730, 388, 342, 0, 878, 369, 161, 306, 1074, 1110, 1130, 347, 163, 817, 637],

  [502, 594, 1278, 514, 400, 776, 1004, 468, 810, 536, 878, 0, 381, 1041, 1004, 1251, 406, 578, 1231, 584, 727, 902],

  [452, 1102, 1152, 300, 511, 270, 144, 476, 263, 1190, 369, 381, 0, 720, 1011, 991, 957, 233, 758, 934, 1054, 154],

  [1020, 298, 819, 263, 886, 343, 353, 1296, 138, 757, 161, 1041, 720, 0, 156, 958, 1202, 919, 401, 685, 716, 994],

  [393, 741, 510, 593, 168, 716, 408, 338, 145, 293, 306, 1004, 1011, 156, 0, 798, 270, 298, 618, 246, 766, 957],

  [208, 582, 514, 153, 514, 535, 220, 1100, 740, 580, 1074, 1251, 991, 958, 798, 0, 1244, 516, 710, 1086, 946, 1098],

  [1196, 712, 399, 801, 396, 1269, 636, 1125, 1251, 1195, 1110, 406, 957, 1202, 270, 1244, 0, 1112, 1031, 1106, 785, 913],

  [983, 1085, 537, 1128, 594, 410, 514, 986, 1262, 124, 1130, 578, 233, 919, 298, 516, 1112, 0, 593, 109, 649, 1150],

  [427, 1024, 579, 729, 886, 1177, 1180, 710, 402, 268, 347, 1231, 758, 401, 618, 710, 1031, 593, 0, 676, 785, 602],

  [482, 202, 509, 1230, 190, 494, 1146, 320, 569, 505, 163, 584, 934, 685, 246, 1086, 1106, 109, 676, 0, 394, 263],

  [479, 220, 285, 317, 729, 1133, 130, 279, 189, 309, 817, 727, 1054, 716, 766, 946, 785, 649, 785, 394, 0, 826],

  [517, 1050, 555, 946, 1042, 651, 648, 460, 877, 602, 637, 902, 154, 994, 957, 1098, 913, 1150, 602, 263, 826, 0],

])

1.3.3 QUBO实现

def build_objective(x: pyqubo.array.Array) -> pyqubo.core.express.AddList:

  H = sum(distance[l_from][l_to] * x[v][l_from] * x[v][s+1][l_to]

          for v in range(V)

          for s in range(S-1)

          for l_from in range(L)

          for l_to in range(L)

          )

  return H

def build_depot_start_rule(x: pyqubo.array.Array) -> pyqubo.core.express.Constraint:

  H = pyqubo.Constraint(

      sum((x[v][0][0] - 1)**2 for v in range(V)),

      'w_0')

  return H

def build_depot_end_rule(x: pyqubo.array.Array) -> pyqubo.core.express.Constraint:

  H = pyqubo.Constraint(

      sum((x[v][S-1][0] - 1)**2 for v in range(V)),

      'w_1')

  return H

def build_vehicle_stay_rule(x: pyqubo.array.Array) -> pyqubo.core.express.Constraint:

  H = pyqubo.Constraint(

      sum(

          (sum(x[v][l] for l in range(L)) - 1)**2

          for v in range(V)

          for s in range(S)),

          'w_2')

  return H

def build_location_visit_rule(x: pyqubo.array.Array) -> pyqubo.core.express.Constraint:

  H = pyqubo.Constraint(

      sum(

          (sum(x[v][l] for v in range(V) for s in range(S)) - 1)**2

          for l in range(1, L)),

          'w_3')

  return H

# Vehicle Routing

x = pyqubo.Array.create('x', shape=(V,S,L), vartype='BINARY')

%%time

H = build_objective(x) + \

  pyqubo.Placeholder('w_1') * build_depot_start_rule(x) + \

  pyqubo.Placeholder('w_2') * build_depot_end_rule(x) + \

  pyqubo.Placeholder('w_3') * build_vehicle_stay_rule(x) + \

  pyqubo.Placeholder('w_4') * build_location_visit_rule(x)



feed_dict = {'w_1': 1000,

             'w_2': 1000,

             'w_3': 1000,

             'w_4': 1000}

model = H.compile()

qubo, constant = model.to_qubo(feed_dict=feed_dict)

 1.3.4 OpenJij求解QUBO目标变量

OpenJij（GitHub链接）是使用普通计算机提供SQA（Simulated Quantum Annealing）API来求解，经常用来搭配OpenJij使用。

sampler = SQASampler(trotter=20, num_reads=10)

response = sampler.sample_qubo(qubo)

1.3.5 输出求解结果路径

def extract_samples(response) -> (list, list):

  solutions = []

  energies = []

  

  for record in response.record:

    sol, num_occ = record[0], record[2]

    solution, broken, energy = model.decode_solution(dict(zip(response.variables, sol)), vartype='BINARY', feed_dict=feed_dict)

    if len(broken) == 0:

      solutions += [solution] * num_occ

      energies += [energy] * num_occ

  

  return solutions, energies



solutions, energies = extract_samples(response)



def vehicle_movement(solution: list) -> list:

  vehicles = []



  for v in range(V):

    solution_sorted = sorted(solution['x'][v].items(), key=lambda x:x[0])

    s = [i[1] for i in solution_sorted]

    df = pd.DataFrame(s).astype(int)

    df.columns = [chr(l) for l in range(65, 65+L)]

    vehicles.append(df)

  

  return vehicles

  

best_solution = solutions[energies.index(min(energies))]

best_vehicles = vehicle_movement(best_solution)



def highlight_positive(val: int) -> str:

  if val == 1:

    return 'color: {}; font-weight: bold'.format('red')

  else:

    return 'color: {}'.format('gray')



for v in range(V):

  display_html('<b>vehicle-{}</b>'.format(v+1) + best_vehicles[v].style.applymap(highlight_positive)._repr_html_() + '<br>', raw=True)

 最后的输出结果如下：

 

1.3 小结

本章讲解以最短距离为优化目标的VRP的求解，约束条件也比较简单，真正的应用里，也有以最短时间为优化目标的，也有很多考虑更多约束条件的：①优先送货服务②车辆载货量有限制③指定时间配送……至于建模以上类似的约束条件，下面两篇文章可以给大家一些启示：

[1]  A Hybrid Solution Method for the Capacitated Vehicle Routing Problem Using a Quantum Annealer   https://arxiv.org/abs/1811.07403

[2] Quantum Annealing of Vehicle Routing Problem with Time, State and Capacity   https://arxiv.org/abs/1903.06322

二、护士调度问题(NSP)

2.1 什么是护士调度问题

调度就是人和机器的调度和排序。由于许多人的日程安排相互交织，还要考虑机器的运行状态等因素，实际业务中出现的日程安排极其复杂。

手动排程时，熟悉现场知识的老手通常要花几个小时来规划一个小规模的项目，甚至可能需要几天时间。但是，如果我们可以把这些问题建模成 组合优化问题，这样便能在几分钟到几十分钟内自动设置复杂的调度。实际的工业界，有工人的工作时间调度和机器的使用时间调度。

下面的论文里可知，NSP 从 1969 年之前就开始被研究，被认为是 NP-hard问题。作为现实生活中常见调度问题，我们会有很多约束，最主要的就是轮班约束和护士约束。本文将介绍护士调度问题中，一些常见约束如何建模为QUBO。

Solos, Ioannis; Tassopoulos, Ioannis; Beligiannis, Grigorios (21 May 2013). “A Generic Two-Phase Stochastic Variable Neighborhood Approach for Effectively Solving the Nurse Rostering Problem”. Algorithms. 6 (2): 278–308.doi:10.3390/a6020278.

2.2 护士调度问题（NSP）的QUBO建模

提示：本节参考了https://qard.is.tohoku.ac.jp/T-Wave/?p=1756，是日本东北大学专门研究量子退火研究室的学生写的量子退火的各种论文解析，这次解析的论文如下：

Title : Application of Quantum Annealing to Nurse Scheduling Problem Author : Kazuki Ikeda, Yuma Nakamura & Travis S.Humble Reference : Scientific Reports 9, Article number: 12837 (2019) https://doi.org/10.1038/s41598-019-49172-3

下面是一个简单的，4个护士，3个工作日的排班示例，我们称下面的表格为排班矩阵：

我们发现右边的排班，会导致护士3没有被安排工作，这是现实生活中不允许的，接下来讲目标变量和约束条件的建模。

2.2.1 目标变量

2.2.2 约束条件定义

因为护士调度问题不存在直接最小化的目标函数，只需要排班矩阵满足以下约束即可。

在 NSP 中可以想象各种约束示例。例如，可以完成的工作量（劳动力）取决于护士的经验。因此，在工作量大的日子里，不可能只有新人在一起工作。因此，我们需要【保证每天所需的劳动力】这一约束条件。可以考虑许多其他约束，但本文作为示例，定义了以下三个约束：

约束a. 确保没有护士连续工作超过 2 天（硬约束）

约束b. 确保每天都有所需的劳动力（硬约束）

约束c. 可以根据护士的忙碌程度调整上班的次数（软约束）

这里把约束分为硬约束和软约束。硬约束是必须满足的约束。另一方面，软约束是不必满足但需要满足的约束。

为了描述这些约束，需要的变量如下表所示：

★变量的注解：

每个约束对应的QUBO项如下：

上面的约束全部满足的话，总哈密顿量的值为0。

2.2.3【约束a】的补充说明

其实我们也可以不使用矩阵J来表示该约束，可以用下面的式子替代：

2.2.4【约束c】的补充说明

约束c的,,通过下面的实例来帮助大家理解。

上面的具体值代入后：

通过该项，我们可以通过每个护士的繁忙程度来优化每个人的出勤时间。

2.3 Python实现NSP的QUBO

本节参考了https://qard.is.tohoku.ac.jp/T-Wave/?p=2764，日本人经常会把实验结果可视化出来，值得借鉴。

2.3.1 参数设定表

因为NSP问题需要先设定一些参数的初始值，所以整理称表格如下：

from itertools import product

import numpy as np

from itertools import product

from pyqubo import Array, Constraint, Placeholder

from neal import SimulatedAnnealingSampler



N = 3  # 护士人数

D = 14  # 待排班日数

a = 7 / 2  # 连续工作2天以上的惩罚值

F = [4, 5, 5]  # 每个护士希望的出勤日数

2.3.2 QUBO实现

惩罚项矩阵J:

J = np.zeros((N, D, N, D))

for n1, d1, n2, d2 in product(range(N), range(D), range(N), range(D)):

    if n1 == n2 and d1 + 1 == d2:

        J[n1, d1, n2, d2] = a

总哈密顿量H

# binary变量

q = Array.create('q', shape=(N, D), vartype='BINARY')



# 连续工作2天以上的惩罚项矩阵

H1 = np.sum([J[n1, d1, n2, d2] * q[n1, d1] * q[n2, d2]

             for n1, n2, d1, d2 in product(range(N), range(N), range(D), range(D))])



# 确保每天d都有一个人的劳动力

H2 = np.sum([(np.sum([q[n,d] for n in range(N)]) - 1)**2 for d in range(D)])



# 确保全员出勤次数相等

H3 = np.sum([(np.sum([q[n,d] for d in range(D)]) - F[n])**2 for n in range(N)])



# 最小化的目标哈密顿量H

H = Placeholder('alpha') * Constraint(H1, 'H1') + Placeholder('lam') * Constraint(H2, 'H2') + Placeholder('gamma') * H3

model = H.compile()

退火算法采样

#这次采用neal的模拟退火

from neal import SimulatedAnnealingSampler

feed_dict = {'alpha': 1.0, 'lam': 1.3, 'gamma': 0.3} # 制約項の係数

bqm = model.to_bqm(feed_dict=feed_dict)

sampler = SimulatedAnnealingSampler()



# D-Wave的量子退火机可以这么使用

#from dwave.system import DWaveSampler, EmbeddingComposite

#sampler_config = {'solver': 'DW_2000Q_6', 'token': 'YOUR_TOKEN'}

#sampler = EmbeddingComposite(DWaveSampler(**sampler_config))



num_reads = 1000

sampleset = sampler.sample(bqm, num_reads=num_reads)

sampleset.record[:10]

输出结果如下：

sampleset.record有以下几部分：

– 解的状态

– 能力值

– 解的个数

rec.array([([1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0], 4.76063633e-13, 1),

           ([0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1], 4.76063633e-13, 1),

           ([1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1], 4.68958206e-13, 1),

           ([0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0], 4.68958206e-13, 1),

           ([1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1], 4.68958206e-13, 1),

           ([0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0], 4.76063633e-13, 1),

           ([0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0], 4.68958206e-13, 1),

           ([0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1], 4.76063633e-13, 1),

           ([0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0], 4.76063633e-13, 1),

           ([0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0], 4.68958206e-13, 1)],

          dtype=[('sample', 'i1', (42,)), ('energy', '<f8'), ('num_occurrences', '<i8')])

2.3.3 实验结果检验和可视化

检查解是否满足约束：

def extract_feasible_samples(samples, print_broken=False):

    feasible_samples = []

    for sample in decoded_samples:

        constraints = sample.constraints(only_broken=True)

        if len(constraints) == 0:

            feasible_samples.append(sample)

        elif print_broken:

            print(constraints)

    return feasible_samples



decoded_samples = model.decode_sampleset(sampleset.aggregate(), feed_dict)

feasible_samples = extract_feasible_samples(decoded_samples)

print('满足约束的解的个数：', len(feasible_samples))

输出结果：

满足约束的解的个数： 894

import matplotlib.pyplot as plt



lowest_sample = feasible_samples[0].sample



schedules = np.zeros((N, D))

for n in range(N):

    for d in range(D):

        if lowest_sample[f'q[{n}][{d}]'] == 1:

            schedules[n, d] = 1



plt.imshow(schedules, cmap="Greys")

plt.xlabel('D')

plt.ylabel('N')

plt.show()

打印出的结果如下：

print('实际的工作日数：', np.sum(schedules, axis=1))

print('希望的工作日数：', F)

实际的工作日数：[4. 5. 5.]

希望的工作日数：[4, 5, 5]

2.4 实验结果

重复实验的这位学生的最后总结如下：

令人失望的是，我无法为多达 56 个变量的问题找到一个好的基态解决方案。通过执行比本文更详细的系数搜索，精度有提高的余地，但随着排班日数的增加，它被认为几乎不可能解决。另一种提高准确性的方法是自己设置QUBO嵌入，但我不知道以我目前的知识做出什么样的嵌入。另外，这样的结果，我觉得很难投入实际使用，因为最多只能安排一个星期的班次。一旦我可以制定至少一个月的时间表，我想进行另一个实验。

其实，量子退火算法在QUBO建模上存在困难的同时，由于量子退火机和伊辛机还没有达到大规模实际问题的全结合bit数，所以，大家也要怀着谨慎的态度看待量子退火的未来。

三、总结

尽管量子退火在小规模问题上展现了较好的潜力，但在建模复杂度、规模扩展性和结果稳定性等方面仍需进一步突破。在量子硬件性能尚未成熟前，仍需要探索可落地的优化路径。

 

---

本文改编转载自CSDN，作者：gang_unerry

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

原文链接[1]：https://blog.csdn.net/gangshen1993/article/details/128164636

原文链接[2]：https://blog.csdn.net/weixin_49155137/article/details/131686724

本篇围绕两类经典组合优化问题——车辆路径问题（VRP）和护士调度问题（NSP），探讨了其在量子计算中的 QUBO 建模与 Python 实现过程，并使用量子退火算法进行求解与验证。同时在实验的基础上，指出量子退火算法的缺陷。

在量子退火算法（2）中，我们学习了整数分割问题和旅行商问题的QUBO 建模与 Python 实现过程。本篇将在系列前文的基础上，探讨车辆路径问题和护士调度问题，具体内容如目录所示。

目录

一、车辆路径问题 (VRP)
      1.1 什么是车辆路径问题
      1.2 车辆路径问题（VRP）的QUBO建模
           1.2.1 目标变量
           1.2.2 目标函数与约束条件
      1.3 Python实现VRP的QUBO
           1.3.1 引入库
           1.3.2 设置参数和距离矩阵
           1.3.2 设置参数和距离矩阵
           1.3.3 QUBO实现
           1.3.4 OpenJij求解QUBO目标变量
           1.3.5 输出求解结果路径
      1.3 小结
二、护士调度问题(NSP)
      2.1 什么是护士调度问题
      2.2 护士调度问题（NSP）的QUBO建模
           2.2.1 目标变量
           2.2.2 约束条件定义
           2.2.3【约束a】的补充说明
           2.2.4【约束c】的补充说明
      2.3 Python实现NSP的QUBO
           2.3.1 参数设定表
           2.3.2 QUBO实现
           2.3.3 实验结果检验和可视化
      2.4 实验结果
三、总结

一、车辆路径问题 (VRP)

1.1 什么是车辆路径问题

QUBO为解决组合优化问题的利器，车辆路径问题是最经常被提起的现实应用。车辆路径问题 (VRP) 是优化多辆送货车辆的送货顺序的组合优化问题，是旅行商问题 (TSP) 的推广。这是一个可以应用于各种现实世界问题的问题，例如物流中的货物配送、工厂中的货物搬运，以及按需运输服务中的车辆分配计划。VRP 是一个典型的组合优化问题，也是经典的NP-Hard问题，普通计算机很难在多项式时间内找到精确解。量子退火有望利用量子隧穿效应实时解决 NP-Hard 组合优化问题。

1.2 车辆路径问题（VRP）的QUBO建模

作为一个基本的VRP，我们考虑了当多辆运送车辆访问所有目标基地时，搜索总旅行成本最小的路线的组合优化问题。我们假设所有的送货车辆都从一个叫做 Depot 的基地出发，并且总是在最后返回到 Depot。移动成本使用基地之间的距离计算。 VRP的QUBO的定义为：把地点1作为Depot，然后假设有2辆车用来配送4个地点。2辆车的配送示意图如下：

其实大家可以可以把VRP问题理解为，有V辆车，V辆车各自负责不重叠的地点子集，所有车辆负责的地点子集的路径总和最小。

所以这里我们有三个参数需要设定：

• V：车辆个数


• P：地点个数


• S：时间步（TSP建模中，我们不需要设置时间步，因为S=P+1）

小写字母代表变量：

1.2.1 目标变量

我们期待的QUBO矩阵如下图所示，

• 红色矩阵代表车辆1，从地点1==>地点3==>地点1；


• 蓝色矩阵代表车辆2，从地点1==>地点2==>地点4==>地点1。

1.2.2 目标函数与约束条件

目标函数：其实就是TSP的目标函数，又增加了一个车辆的维度

约束条件：

条件a：第1个时间步的地点必须是Depot地点
条件b：第S个时间步的地点必须是Depot地点
条件c：除了Depot地点，剩余地点必须被访问过一次
条件d：一个时间步，车辆v不能同时出现在多个地点p

各个约束条件对应的表达式如下图：

最后转换后的哈密顿量为：

1.3 Python实现VRP的QUBO

本章节主要来自以下代码库：

https://github.com/yskmt2018/quantum2020/blob/master/vehicle_routing_problem.ipynb

1.3.1 引入库

import math

import pyqubo

import numpy as np

import pandas as pd

from openjij import SQASampler

1.3.2 设置参数和距离矩阵

# 车辆个数

V = 3

# 地点个数

L = len(distance[0])

# 最大时间步

S = 10



# 距离矩阵

distance = np.array([

  [0, 548, 776, 696, 582, 274, 502, 194, 308, 194, 536, 502, 452, 1020, 393, 208, 1196, 983, 427, 482, 479, 517],

  [548, 0, 684, 308, 194, 502, 730, 354, 696, 742, 1084, 594, 1102, 298, 741, 582, 712, 1085, 1024, 202, 220, 1050],

  [776, 684, 0, 992, 878, 502, 274, 810, 468, 742, 400, 1278, 1152, 819, 510, 514, 399, 537, 579, 509, 285, 555],

  [696, 308, 992, 0, 114, 650, 878, 502, 844, 890, 1232, 514, 300, 263, 593, 153, 801, 1128, 729, 1230, 317, 946],

  [582, 194, 878, 114, 0, 536, 764, 388, 730, 776, 1118, 400, 511, 886, 168, 514, 396, 594, 886, 190, 729, 1042],

  [274, 502, 502, 650, 536, 0, 228, 308, 194, 240, 582, 776, 270, 343, 716, 535, 1269, 410, 1177, 494, 1133, 651],

  [502, 730, 274, 878, 764, 228, 0, 536, 194, 468, 354, 1004, 144, 353, 408, 220, 636, 514, 1180, 1146, 130, 648],

  [194, 354, 810, 502, 388, 308, 536, 0, 342, 388, 730, 468, 476, 1296, 338, 1100, 1125, 986, 710, 320, 279, 460],

  [308, 696, 468, 844, 730, 194, 194, 342, 0, 274, 388, 810, 263, 138, 145, 740, 1251, 1262, 402, 569, 189, 877],

  [194, 742, 742, 890, 776, 240, 468, 388, 274, 0, 342, 536, 1190, 757, 293, 580, 1195, 124, 268, 505, 309, 602],

  [536, 1084, 400, 1232, 1118, 582, 354, 730, 388, 342, 0, 878, 369, 161, 306, 1074, 1110, 1130, 347, 163, 817, 637],

  [502, 594, 1278, 514, 400, 776, 1004, 468, 810, 536, 878, 0, 381, 1041, 1004, 1251, 406, 578, 1231, 584, 727, 902],

  [452, 1102, 1152, 300, 511, 270, 144, 476, 263, 1190, 369, 381, 0, 720, 1011, 991, 957, 233, 758, 934, 1054, 154],

  [1020, 298, 819, 263, 886, 343, 353, 1296, 138, 757, 161, 1041, 720, 0, 156, 958, 1202, 919, 401, 685, 716, 994],

  [393, 741, 510, 593, 168, 716, 408, 338, 145, 293, 306, 1004, 1011, 156, 0, 798, 270, 298, 618, 246, 766, 957],

  [208, 582, 514, 153, 514, 535, 220, 1100, 740, 580, 1074, 1251, 991, 958, 798, 0, 1244, 516, 710, 1086, 946, 1098],

  [1196, 712, 399, 801, 396, 1269, 636, 1125, 1251, 1195, 1110, 406, 957, 1202, 270, 1244, 0, 1112, 1031, 1106, 785, 913],

  [983, 1085, 537, 1128, 594, 410, 514, 986, 1262, 124, 1130, 578, 233, 919, 298, 516, 1112, 0, 593, 109, 649, 1150],

  [427, 1024, 579, 729, 886, 1177, 1180, 710, 402, 268, 347, 1231, 758, 401, 618, 710, 1031, 593, 0, 676, 785, 602],

  [482, 202, 509, 1230, 190, 494, 1146, 320, 569, 505, 163, 584, 934, 685, 246, 1086, 1106, 109, 676, 0, 394, 263],

  [479, 220, 285, 317, 729, 1133, 130, 279, 189, 309, 817, 727, 1054, 716, 766, 946, 785, 649, 785, 394, 0, 826],

  [517, 1050, 555, 946, 1042, 651, 648, 460, 877, 602, 637, 902, 154, 994, 957, 1098, 913, 1150, 602, 263, 826, 0],

])

1.3.3 QUBO实现

def build_objective(x: pyqubo.array.Array) -> pyqubo.core.express.AddList:

  H = sum(distance[l_from][l_to] * x[v][l_from] * x[v][s+1][l_to]

          for v in range(V)

          for s in range(S-1)

          for l_from in range(L)

          for l_to in range(L)

          )

  return H

def build_depot_start_rule(x: pyqubo.array.Array) -> pyqubo.core.express.Constraint:

  H = pyqubo.Constraint(

      sum((x[v][0][0] - 1)**2 for v in range(V)),

      'w_0')

  return H

def build_depot_end_rule(x: pyqubo.array.Array) -> pyqubo.core.express.Constraint:

  H = pyqubo.Constraint(

      sum((x[v][S-1][0] - 1)**2 for v in range(V)),

      'w_1')

  return H

def build_vehicle_stay_rule(x: pyqubo.array.Array) -> pyqubo.core.express.Constraint:

  H = pyqubo.Constraint(

      sum(

          (sum(x[v][l] for l in range(L)) - 1)**2

          for v in range(V)

          for s in range(S)),

          'w_2')

  return H

def build_location_visit_rule(x: pyqubo.array.Array) -> pyqubo.core.express.Constraint:

  H = pyqubo.Constraint(

      sum(

          (sum(x[v][l] for v in range(V) for s in range(S)) - 1)**2

          for l in range(1, L)),

          'w_3')

  return H

# Vehicle Routing

x = pyqubo.Array.create('x', shape=(V,S,L), vartype='BINARY')

%%time

H = build_objective(x) + \

  pyqubo.Placeholder('w_1') * build_depot_start_rule(x) + \

  pyqubo.Placeholder('w_2') * build_depot_end_rule(x) + \

  pyqubo.Placeholder('w_3') * build_vehicle_stay_rule(x) + \

  pyqubo.Placeholder('w_4') * build_location_visit_rule(x)



feed_dict = {'w_1': 1000,

             'w_2': 1000,

             'w_3': 1000,

             'w_4': 1000}

model = H.compile()

qubo, constant = model.to_qubo(feed_dict=feed_dict)

 1.3.4 OpenJij求解QUBO目标变量

OpenJij（GitHub链接）是使用普通计算机提供SQA（Simulated Quantum Annealing）API来求解，经常用来搭配OpenJij使用。

sampler = SQASampler(trotter=20, num_reads=10)

response = sampler.sample_qubo(qubo)

1.3.5 输出求解结果路径

def extract_samples(response) -> (list, list):

  solutions = []

  energies = []

  

  for record in response.record:

    sol, num_occ = record[0], record[2]

    solution, broken, energy = model.decode_solution(dict(zip(response.variables, sol)), vartype='BINARY', feed_dict=feed_dict)

    if len(broken) == 0:

      solutions += [solution] * num_occ

      energies += [energy] * num_occ

  

  return solutions, energies



solutions, energies = extract_samples(response)



def vehicle_movement(solution: list) -> list:

  vehicles = []



  for v in range(V):

    solution_sorted = sorted(solution['x'][v].items(), key=lambda x:x[0])

    s = [i[1] for i in solution_sorted]

    df = pd.DataFrame(s).astype(int)

    df.columns = [chr(l) for l in range(65, 65+L)]

    vehicles.append(df)

  

  return vehicles

  

best_solution = solutions[energies.index(min(energies))]

best_vehicles = vehicle_movement(best_solution)



def highlight_positive(val: int) -> str:

  if val == 1:

    return 'color: {}; font-weight: bold'.format('red')

  else:

    return 'color: {}'.format('gray')



for v in range(V):

  display_html('<b>vehicle-{}</b>'.format(v+1) + best_vehicles[v].style.applymap(highlight_positive)._repr_html_() + '<br>', raw=True)

 最后的输出结果如下：

 

1.3 小结

本章讲解以最短距离为优化目标的VRP的求解，约束条件也比较简单，真正的应用里，也有以最短时间为优化目标的，也有很多考虑更多约束条件的：①优先送货服务②车辆载货量有限制③指定时间配送……至于建模以上类似的约束条件，下面两篇文章可以给大家一些启示：

[1]  A Hybrid Solution Method for the Capacitated Vehicle Routing Problem Using a Quantum Annealer   https://arxiv.org/abs/1811.07403

[2] Quantum Annealing of Vehicle Routing Problem with Time, State and Capacity   https://arxiv.org/abs/1903.06322

二、护士调度问题(NSP)

2.1 什么是护士调度问题

调度就是人和机器的调度和排序。由于许多人的日程安排相互交织，还要考虑机器的运行状态等因素，实际业务中出现的日程安排极其复杂。

手动排程时，熟悉现场知识的老手通常要花几个小时来规划一个小规模的项目，甚至可能需要几天时间。但是，如果我们可以把这些问题建模成 组合优化问题，这样便能在几分钟到几十分钟内自动设置复杂的调度。实际的工业界，有工人的工作时间调度和机器的使用时间调度。

下面的论文里可知，NSP 从 1969 年之前就开始被研究，被认为是 NP-hard问题。作为现实生活中常见调度问题，我们会有很多约束，最主要的就是轮班约束和护士约束。本文将介绍护士调度问题中，一些常见约束如何建模为QUBO。

Solos, Ioannis; Tassopoulos, Ioannis; Beligiannis, Grigorios (21 May 2013). “A Generic Two-Phase Stochastic Variable Neighborhood Approach for Effectively Solving the Nurse Rostering Problem”. Algorithms. 6 (2): 278–308.doi:10.3390/a6020278.

2.2 护士调度问题（NSP）的QUBO建模

提示：本节参考了https://qard.is.tohoku.ac.jp/T-Wave/?p=1756，是日本东北大学专门研究量子退火研究室的学生写的量子退火的各种论文解析，这次解析的论文如下：

Title : Application of Quantum Annealing to Nurse Scheduling Problem Author : Kazuki Ikeda, Yuma Nakamura & Travis S.Humble Reference : Scientific Reports 9, Article number: 12837 (2019) https://doi.org/10.1038/s41598-019-49172-3

下面是一个简单的，4个护士，3个工作日的排班示例，我们称下面的表格为排班矩阵：

我们发现右边的排班，会导致护士3没有被安排工作，这是现实生活中不允许的，接下来讲目标变量和约束条件的建模。

2.2.1 目标变量

2.2.2 约束条件定义

因为护士调度问题不存在直接最小化的目标函数，只需要排班矩阵满足以下约束即可。

在 NSP 中可以想象各种约束示例。例如，可以完成的工作量（劳动力）取决于护士的经验。因此，在工作量大的日子里，不可能只有新人在一起工作。因此，我们需要【保证每天所需的劳动力】这一约束条件。可以考虑许多其他约束，但本文作为示例，定义了以下三个约束：

约束a. 确保没有护士连续工作超过 2 天（硬约束）

约束b. 确保每天都有所需的劳动力（硬约束）

约束c. 可以根据护士的忙碌程度调整上班的次数（软约束）

这里把约束分为硬约束和软约束。硬约束是必须满足的约束。另一方面，软约束是不必满足但需要满足的约束。

为了描述这些约束，需要的变量如下表所示：

★变量的注解：

每个约束对应的QUBO项如下：

上面的约束全部满足的话，总哈密顿量的值为0。

2.2.3【约束a】的补充说明

其实我们也可以不使用矩阵J来表示该约束，可以用下面的式子替代：

2.2.4【约束c】的补充说明

约束c的,,通过下面的实例来帮助大家理解。

上面的具体值代入后：

通过该项，我们可以通过每个护士的繁忙程度来优化每个人的出勤时间。

2.3 Python实现NSP的QUBO

本节参考了https://qard.is.tohoku.ac.jp/T-Wave/?p=2764，日本人经常会把实验结果可视化出来，值得借鉴。

2.3.1 参数设定表

因为NSP问题需要先设定一些参数的初始值，所以整理称表格如下：

from itertools import product

import numpy as np

from itertools import product

from pyqubo import Array, Constraint, Placeholder

from neal import SimulatedAnnealingSampler



N = 3  # 护士人数

D = 14  # 待排班日数

a = 7 / 2  # 连续工作2天以上的惩罚值

F = [4, 5, 5]  # 每个护士希望的出勤日数

2.3.2 QUBO实现

惩罚项矩阵J:

J = np.zeros((N, D, N, D))

for n1, d1, n2, d2 in product(range(N), range(D), range(N), range(D)):

    if n1 == n2 and d1 + 1 == d2:

        J[n1, d1, n2, d2] = a

总哈密顿量H

# binary变量

q = Array.create('q', shape=(N, D), vartype='BINARY')



# 连续工作2天以上的惩罚项矩阵

H1 = np.sum([J[n1, d1, n2, d2] * q[n1, d1] * q[n2, d2]

             for n1, n2, d1, d2 in product(range(N), range(N), range(D), range(D))])



# 确保每天d都有一个人的劳动力

H2 = np.sum([(np.sum([q[n,d] for n in range(N)]) - 1)**2 for d in range(D)])



# 确保全员出勤次数相等

H3 = np.sum([(np.sum([q[n,d] for d in range(D)]) - F[n])**2 for n in range(N)])



# 最小化的目标哈密顿量H

H = Placeholder('alpha') * Constraint(H1, 'H1') + Placeholder('lam') * Constraint(H2, 'H2') + Placeholder('gamma') * H3

model = H.compile()

退火算法采样

#这次采用neal的模拟退火

from neal import SimulatedAnnealingSampler

feed_dict = {'alpha': 1.0, 'lam': 1.3, 'gamma': 0.3} # 制約項の係数

bqm = model.to_bqm(feed_dict=feed_dict)

sampler = SimulatedAnnealingSampler()



# D-Wave的量子退火机可以这么使用

#from dwave.system import DWaveSampler, EmbeddingComposite

#sampler_config = {'solver': 'DW_2000Q_6', 'token': 'YOUR_TOKEN'}

#sampler = EmbeddingComposite(DWaveSampler(**sampler_config))



num_reads = 1000

sampleset = sampler.sample(bqm, num_reads=num_reads)

sampleset.record[:10]

输出结果如下：

sampleset.record有以下几部分：

– 解的状态

– 能力值

– 解的个数

rec.array([([1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0], 4.76063633e-13, 1),

           ([0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1], 4.76063633e-13, 1),

           ([1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1], 4.68958206e-13, 1),

           ([0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0], 4.68958206e-13, 1),

           ([1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1], 4.68958206e-13, 1),

           ([0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0], 4.76063633e-13, 1),

           ([0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0], 4.68958206e-13, 1),

           ([0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1], 4.76063633e-13, 1),

           ([0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0], 4.76063633e-13, 1),

           ([0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0], 4.68958206e-13, 1)],

          dtype=[('sample', 'i1', (42,)), ('energy', '<f8'), ('num_occurrences', '<i8')])

2.3.3 实验结果检验和可视化

检查解是否满足约束：

def extract_feasible_samples(samples, print_broken=False):

    feasible_samples = []

    for sample in decoded_samples:

        constraints = sample.constraints(only_broken=True)

        if len(constraints) == 0:

            feasible_samples.append(sample)

        elif print_broken:

            print(constraints)

    return feasible_samples



decoded_samples = model.decode_sampleset(sampleset.aggregate(), feed_dict)

feasible_samples = extract_feasible_samples(decoded_samples)

print('满足约束的解的个数：', len(feasible_samples))

输出结果：

满足约束的解的个数： 894

import matplotlib.pyplot as plt



lowest_sample = feasible_samples[0].sample



schedules = np.zeros((N, D))

for n in range(N):

    for d in range(D):

        if lowest_sample[f'q[{n}][{d}]'] == 1:

            schedules[n, d] = 1



plt.imshow(schedules, cmap="Greys")

plt.xlabel('D')

plt.ylabel('N')

plt.show()

打印出的结果如下：

print('实际的工作日数：', np.sum(schedules, axis=1))

print('希望的工作日数：', F)

实际的工作日数：[4. 5. 5.]

希望的工作日数：[4, 5, 5]

2.4 实验结果

重复实验的这位学生的最后总结如下：

令人失望的是，我无法为多达 56 个变量的问题找到一个好的基态解决方案。通过执行比本文更详细的系数搜索，精度有提高的余地，但随着排班日数的增加，它被认为几乎不可能解决。另一种提高准确性的方法是自己设置QUBO嵌入，但我不知道以我目前的知识做出什么样的嵌入。另外，这样的结果，我觉得很难投入实际使用，因为最多只能安排一个星期的班次。一旦我可以制定至少一个月的时间表，我想进行另一个实验。

其实，量子退火算法在QUBO建模上存在困难的同时，由于量子退火机和伊辛机还没有达到大规模实际问题的全结合bit数，所以，大家也要怀着谨慎的态度看待量子退火的未来。

三、总结

尽管量子退火在小规模问题上展现了较好的潜力，但在建模复杂度、规模扩展性和结果稳定性等方面仍需进一步突破。在量子硬件性能尚未成熟前，仍需要探索可落地的优化路径。

 

---

本文改编转载自CSDN，作者：gang_unerry

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

原文链接[1]：https://blog.csdn.net/gangshen1993/article/details/128164636

原文链接[2]：https://blog.csdn.net/weixin_49155137/article/details/131686724