Quadratic unconstrained binary optimization，QUBO中文名是二次无约束二元优化，它是在二次规划(QP, Quadratic Programming)的基础上添加了两个限制条件：（1）只有目标函数，没有约束条件，例如等式约束、不等式约束等；（2）决策变量的取值只能是0和1。

下面给出它的标准形式：

图片来源：参考文献1，见文章最后

从上图可知：

（1）x是决策变量，只有两种取值情况，通常用0和1表示。Q是一个常数组成的方阵！

（2）QUBO有两种形式，取决于矩阵Q是一个对称的矩阵还是一个上三角矩阵。两种形式是等价的，对称的形式看起来更加直观。

下面来看一个实际的例子：

大家可以思考：

（1）有几个决策变量？

（2）矩阵Q的大小是多少？

（3）xi的平方和xi有什么关系？

（4）你能试着写出Q矩阵吗？

下面给出答案：

如果大家学过线性代数的话，应该会想到二次型变换：

再次强调：除了决策变量是01约束之外，不能加入其他的约束条件，所有的约束条件都要想办法转换到Q矩阵中。

此外，许多其他看起来与QUBO问题无关的问题可以重新表述为QUBO模型。下面看具体的例子。

Number Partitioning Problem 问题，也称为数字划分问题

给定一个集合S = {3,1,1,2,2,1}

如何将其元素划分给两个集合，使得两个集合中的元素的和相等。

例如集合1：S(1) = {1,1,1,2} 、集合2：S(2) = {2,3}

此时每个集合中的元素相加都为5，这就是问题的一个解。

然而，有时候无论如何划分两个集合的和也不会相等，

例如：S = {8,1,1,2,2,1}，此时我们需要找到两个集合中的元素的和尽可能接近的一组解。

下面我们来看建模过程：

假设S中有m个元素，分别是s1,s2,...,sm，如果第j个元素sj属于第一个集合，那么sj就等于1（否则sj属于第二个集合，此时sj等于0）。

那么，第一个集合中的元素和就是sum1，第二个集合中的元素和就可以用S中所有的元素和减去sum1,得到的结果是sum2。

两个集合的元素和的差异是：

c是S中所有元素的和，因此是一个常数。

我们的目标函数是使得diff尽可能接近于0，因此有两种处理方法：

（1）min abs(diff)  求diff的绝对值，绝对值越小越好。

（2）min (diff)^2   求diff的平方，平方越小越好。

那么，哪个目标函数和QUBO的关系更接近呢？

很明显，第二种方法和QUBO的关系更大！

注意：这里要将平方项展开（具体的展开形式论文中省略了，大家有时间可以自己推导看看，应该有点复杂），然后得到Q矩阵，Q矩阵是一个m阶的对称矩阵！里面所有元素都是已知的。

下面给了一个数值计算的例子：

对于上面这个优化问题，有很多种求解方法：

（1）使用MATLAB求解

由于MATLAB没有求解QUBO的内置函数，因此需要借助第三方求解器。第三方求解器有很多，这里推荐一个免费的工具箱：OPTI Toolbox，这个工具箱中集成了很多第三方的求解器。官网下载地址：

https://www.controlengineering.co.nz/Wikis/OPTI/pmwiki.php/DL/DownloadOPTI

OPTI Toolbox集成了一系列优化求解器，安装好后在Matlab命令窗口输入optiSolver，即可看到自带的求解器。

然而，对于QUBO优化问题，工具箱中自带的求解器很弱，容易求出局部最优解，因此我们我们需要补充安装一个更强大的SCIP求解器。

学术版本的SCIP可以免费使用，只需要提交一个邮箱即可下载scip.mexw64，下载好后需要把文件放到OPTI的Solvers文件夹。

接下来以我们这个问题为例，试试工具箱的计算。

OPTI官网给出了支持的标准的MIQP（混合整数二次规划问题）形式：

https://www.controlengineering.co.nz/Wikis/OPTI/pmwiki.php/Probs/MIQP

我们要求解的QUBO可以看成MIQP问题的特例！

第一步：得到Q矩阵：

第二步：调用求解器

这个警告的意思是，目标函数可能是非凸的，也就是说，它可能有多个局部最优解，而不是一个全局最优解。这样的话，你得到的解可能不是全局最优解，而是一个局部最优解。

从倒数第二段可以看出，使用的是SCIP工具箱，算法是空间分枝定界法。

最优解对应的目标函数是-3444.5（实际上要乘以2等于-6889，工具箱的标准型上将目标函数乘以了1/2）。

两个集合包含的元素：

这个答案和参考资料里面给的答案不同：

但是最优结果是一样的，这说明我们这个问题有多个解。

（2）使用Python求解

这里推荐一个很强大的第三方库：PyQUBO。

官网：https://pyqubo.readthedocs.io/en/latest/getting_started.html

它可以将组合优化问题映射为QUBO形式，刚好找到了一篇文章介绍这个库的使用：

https://cloud.tencent.com/developer/article/2111550

如果你将QUBO看成MIQP问题的特例，那么可以求解的第三方包就更多了，例如CVXPY、pyMIQP等。

（3）使用lingo求解

LINGO是Linear Interactive and General Optimizer的缩写，即“交互式的线性和通用优化求解器”，由美国LINDO系统公司（Lindo System Inc.）推出的，可以用于求解非线性规划，也可以用于一些线性和非线性方程组的求解等，功能十分强大，是求解优化模型的最佳选择。其特色在于内置建模语言，提供了许多内部函数，可以允许决策变量是整数（即整数规划，包括 0-1 整数规划），方便灵活，而且执行速度非常快。

关于lingo求解二次整数规划的问题，大家可以自己查阅相关资料。

有约束条件时如何转换？

注意，上面我们的例子中除了要求x是01变量外，是没有其他约束的，然而在实际例子中，往往有其他的约束限制，我们需要对这些约束进行转换，转换的思想很简单，假设我们的目的是计算目标函数的最小值，那么我们可以在目标函数上增加一个非负的惩罚函数（设计为关于决策变量的二次函数形式），该函数满足下面两点要求：

（1）如果当前解满足约束，那么惩罚函数的值等于0；

（2）如果当前解不满足约束，那么惩罚函数是一个很大的正数。

我们来看参考文献的介绍：

P是一个很大的正数，如果违反约束的话，目标函数会受到一个很大的惩罚，以至于我们最终得到的解会倾向于符合约束。

我们还是以之前的数字划分问题为例：

假如我们要求第一个元素25和第四个元素31要分布在不同的集合中，即要求：x1+x4=1

那么我们可以根据上表在原目标函数的基础上要添加惩罚：

不妨假设P=10000，原目标函数中的Q需要进行更改：

那么Q中第一行第一列的元素以及第四行第四列的元素都需要减去10000，第一行第四列的元素以及第四行第一列的元素则需要增加10000。

P*1是一个常数可以暂时忽略，不影响最优解。

最优的方案发生了变化，但是最优的目标函数没有变化：

大家可以思考下：P值的选择有什么要求吗？

（1）P值过小可能有什么问题？

（2）P值过大可能有什么问题？

P值的选取技巧：

英文不好的同学可以看下面的机器翻译：

正如我们所指出的，许多问题转换为QUBO形式需要引入一个惩罚参数P，而且P必须给出一个具体的数值。惩罚参数P值并不是唯一的，可以选择许多不同的P值。对于一个特定的问题，一个可行的P值通常是根据领域知识和需要完成的目标来设定的。通常，我们对所有的约束使用相同的惩罚，但是如果有充分的理由区别对待不同的约束，那么对不同的约束使用不同的惩罚也没有问题。如果一个约束必须绝对满足，即“硬”约束，那么P必须足够大，以防止违反。然而，有些约束是“软”的，意味着满足它们是可取的，但是可以容忍轻微的违反。对于这样的情况，一个更适度的惩罚值就足够了。 

一个太大的惩罚值会阻碍解决过程，因为惩罚项压倒了原始目标函数的信息，使得难以区分一个解和另一个解的质量。另一方面，一个太小的惩罚值会危及寻找可行解的过程。一般来说，有一个相当大的“Goldilocks region”，包含了可以很好工作的惩罚值。对模型进行一些初步思考可以得到原始目标函数值的一个大致估计。将P取为这个估计的一定百分比（75%到150%）通常是一个好的起点。最后，总是可以检查生成的解是否可行，从而根据需要改变惩罚值并进行进一步的解决过程，以找到一个可接受的解。

简单点说就是提前估计一个最优的目标函数，然后将P值设置在这个值的周围，如果最后求出来的解违反了你的约束条件，那么你可以尝试增加P的大小来提高违反约束的惩罚程度。

思考题：将下面的约束问题转换为QUBO问题，式中xi均是变量。

有两个约束，均需要转换为惩罚函数：

转换后：

存在线性约束怎么转换？

一样的加入惩罚项：

例如：

下面给出具体的求解过程：

如果有线性不等式约束应该怎么处理？

这里面涉及到运筹学里面的概念：松弛变量。

它的作用是把不等式约束转换成等式约束，当约束条件为“≤”（或者“≥”）类型时的不等式约束时，可在不等式左边加上（或者减去）一个非负的新变量s，即可将这个不等式约束转换成等式约束。这个新增的非负变量称为松弛变量（或剩余变量），也可统称为松弛变量。在目标函数中一般认为新增的松弛变量的系数为0。

下面我们来看参考资料中给的例子：

显然，第一个约束和第三个约束都是不等式约束，我们需要进行转换：

对于第一个约束条件，可以在不等式左边加上一个非负的新变量s1，即可化为等式：2x1+2x2+4x3+3x4+2x5+s1=7，其中s1≥0。

对于第三个约束条件，可以在不等式左边减去一个非负的新变量s3，即可化为等式：3x1+3x2+2x3+4x4+4x5-s3=5，其中s3≥0。

由于我们要求解QUBO问题，因此这里的s1和s3不能直接写成这种形式，需要将其进行转换：转换成由01变量构成。

为此我们需要先估计出s1和s2的上界：

为啥s1的上界是3呢？我们可以将2x1+2x2+4x3+3x4+2x5+s1=7减去第二个等式约束x1+2x2+2x3+x4+2x5=4，得到：x1+2x3+2x4+s1 = 3，即s1 =3-(x1+2x3+2x4)

又因为x1+2x3+2x4理论上的最小值是0，因此可以估计出s1的上界是3。

为啥s3的上界是6呢？我们将3x1+3x2+2x3+4x4+4x5-s3=5减去两倍的第二个等式约束，得到：x1-x2-2x3+2x4-s3 = -3，那么s3 =x1-x2-2x3+2x4+3，当x1和x4为1，x2和x3为0时，可以估计出s3的上界是6。

（注意：这里估计的上界可以更大，比如s3的上界我们还可以这样估计：

由3x1+3x2+2x3+4x4+4x5-s3=5得s3 = 3x1+3x2+2x3+4x4+4x5-5，当xi全取1时s3为11，那么为什么我们不将s3的上界估计成11呢？这是因为s3的上界我们给的越大，后续求解时你的解空间中也可能遇到更多的情况，会降低搜索的效率，因此这里松弛变量的上界我们尽可能去找一个小的满足原来约束的，这里需要很多的技巧！）

接下来对每个松弛变量进行了变形，让他们变成01变量的和：

注意：为啥s1和s3要这样拆分呢？能不能直接让s1取为3x6呢？注意，如果你直接让s1=3x6，那么就限制了s1只能为0或者3，但实际上s1是可能取为1和2的，因此让s1=3x6不可取；那能否让s1=x6+x7+x8呢？答案是肯定的，但这样多引入了三个变量，而直接让s1=x6+2x7既可以保证s1取到0 1 2 3，又可以只引入两个变量，可以节省计算的开销。

类似的，s3 = x8+2x9+4x10也能保证s3取到0 1 2 3 4 5 6，又是引入新的变量最少的一种方案。

后续的计算步骤就很简单了，引入惩罚项转换为标准的QUBO形式：

二次分配（指派）问题

很明显这里面存在n的平方个01决策变量，因此当n很大时问题的求解规模很大。下面给出n=3的例子：

这里将双下标转换为单下标的形式，然后进行建模：

添加惩罚项：

以上介绍来自参考论文：

GLOVER, FRED, KOCHENBERGER, GARY, DU, YU. Quantum Bridge Analytics I: a tutorial on formulating and using QUBO models[J]. 4OR: Quarterly Journal of the Belgian, French and Italian Operations Research Societies,2019,17(4):335-371. DOI:10.1007/s10288-019-00424-y.

上面这篇文章的内容对于初学者而言是一个很好的资料，下面这篇博客是一个学者整理的QUBO在经典优化问题上的应用，有100多个经典优化问题，包括最大割、图着色、旅行商问题、背包问题等：

https://blog.xa0.de/post/List-of-QUBO-formulations/

里面比较有名的问题有：

最大割问题（Max-Cut Problem）

图着色问题（Graph Coloring Problem）

旅行商问题（Travelling Salesman Problem）

背包问题（Knapsack Problem）

作业调度问题（Job Shop Scheduling Problem）

车间调度问题（Flow Shop Scheduling Problem）

二次指派问题（Quadratic Assignment Problem）

集合覆盖问题（Set Covering Problem）

博客后面给出了这些问题的参考文献：