1. 项目概述从数学建模到量子计算前沿的实战跨越去年带队参加MathorCup选了A题最后拿了个一等奖。赛题的核心是把一个经典的组合优化问题——大规模网络中的资源分配与路径规划转化成了量子计算领域里一个非常热门的模型QUBO。这听起来有点跨界数学建模怎么就和量子计算扯上关系了其实这正是当前交叉学科研究的一个缩影。我们队伍里既有搞传统运筹优化的也有对量子计算感兴趣的大家一合计觉得这个方向既有挑战性未来潜力也大就决定啃下这块硬骨头。简单来说QUBO模型全称是“二次无约束二值优化”它之所以在量子计算圈特别是量子退火和相干伊辛机这类专用硬件上备受青睐是因为它的形式极其简洁所有变量只能是0或1目标函数是一个二次型。很多NP难问题比如旅行商问题、图着色、设备布局等都能被“编码”成QUBO形式。一旦编码成功理论上就可以丢给量子退火机去求解利用量子隧穿、量子纠缠等效应在特定问题上寻找比经典算法更优的解。我们这次做的就是把一个实际的、带有多重约束的数学建模问题“翻译”成QUBO这门“量子硬件能听懂的语言”。整个过程充满了探索和试错。网上关于QUBO建模的教程大多停留在理论和小型玩具案例真正把一个中等规模、带有复杂现实约束的MathorCup赛题完整走通的分享并不多。我们花了大量时间在“建模转化”和“求解策略”上如何精确地将约束条件“惩罚”进目标函数如何调整惩罚权重避免无效解以及在没有真实量子硬件的情况下如何用经典模拟器和启发式算法来验证和求解我们构建的QUBO模型。这篇文章我就把这套从问题分析、QUBO建模、到经典求解与结果分析的完整思路和核心代码分享出来希望能给未来想涉足量子优化或者参加类似交叉学科竞赛的朋友们一些实实在在的参考。2. 赛题核心与QUBO建模思路拆解2.1 问题重述与经典建模视角我们遇到的A题本质上是一个多商品流网络优化问题。场景可以抽象为一个有向图节点代表不同的资源点或任务点边代表连接路径每条边有容量限制。有多种不同的“商品”或任务需要在网络中从各自的起点运输到终点每种商品有确定的流量需求。目标是在满足所有边容量约束的前提下找到一种路由方案使得总运输成本或距离、时间等最小。如果完全用经典方法比如混合整数线性规划来建模我们会定义一系列决策变量x_{ijk}表示商品k是否经过边(i, j)流量是多少。然后列出流量守恒约束、容量约束、变量类型约束构建一个大规模的MILP模型丢给CPLEX、Gurobi这类求解器去算。对于规模不大的问题这很有效。但一旦节点和商品数量上去问题规模会指数级膨胀求解时间可能变得不可接受。这时量子计算特别是量子退火提供了一种不同的思路。它不擅长直接处理复杂的线性约束但它非常擅长在巨大的解空间中进行“能量”最小化的搜索。QUBO模型正是描述这种“能量景观”的完美数学形式。我们的核心任务就是将带有约束的经典优化问题转化为一个无约束的二次二值优化问题。转化的关键在于把所有的硬约束必须满足的条件都转换成“惩罚项”加到目标函数里。如果某个解违反了约束它的“能量”即目标函数值就会因为惩罚项而急剧增高从而在最小化过程中被淘汰。2.2 QUBO模型的核心形式与转化哲学QUBO模型的标准形式如下Minimize: y x^T Q x其中x是一个由二值变量0或1组成的列向量Q是一个实对称矩阵或上三角矩阵。我们的所有智慧都体现在如何构造这个Q矩阵上。对于我们的网络流问题直接使用0/1变量表示“是否选择某条边”是自然的。但难点在于流量是连续值而QUBO变量是离散二值的。这里就需要引入二进制展开的技巧。例如如果一条边的容量上限是7我们可以用3个量子比特对应3个二值变量的二进制组合来表示这条边上分配的流量0到7。假设变量x1, x2, x3分别代表2^0, 2^1, 2^2位那么流量f x1 2*x2 4*x3。这样我们就用一组二值变量描述了一个有限范围内的整数量。接下来是重头戏约束的惩罚项构造。以边容量约束为例对于边e其上的总流量f_e不能超过容量C_e。我们可以将其转化为惩罚项P * (f_e - C_e)^2但这里f_e - C_e可能为负未超容平方后仍为正这会错误地惩罚未超容的好解。因此更标准的做法是引入松弛变量。定义s_e为一个非负整数同样用二进制展开表示将约束改写为f_e s_e C_e。那么惩罚项就是P * (f_e s_e - C_e)^2。当约束被严格满足时此项为0否则为正起到惩罚作用。惩罚系数P是一个需要精心调整的超参数它必须足够大以确保任何违反约束的解的目标值都差于任何可行解但又不能太大以免数值溢出或导致能量地形过于陡峭难以优化。流量守恒约束对于每个商品在每个中间节点流入等于流出是另一个核心约束处理方式类似也需要为每个节点和每个商品引入类似的等式约束惩罚项。最终我们的目标函数将由两部分组成原始的最小化总成本项通常是流量的线性函数通过二进制展开后转化为二值变量的二次型和所有约束对应的惩罚项。将它们全部加起来就得到了最终的QUBO形式x^T Q x。构建Q矩阵的过程就是将这些二次项和一次项可以转化为二次项如x_i可视为x_i * x_i的系数整理到一个矩阵里的过程。注意惩罚系数P的选取是艺术也是科学。一个实用的技巧是让其大小与原始目标函数值的量级相匹配。可以先忽略惩罚项求解一个松弛问题估算目标值范围然后设定P比这个范围大一个数量级例如10倍再进行调试。我们当时采用了迭代调整策略从一个小值开始逐步增加直到求得的解100%满足所有约束。3. 模型构建的详细步骤与关键实现3.1 决策变量定义与二进制编码这是整个建模的基石变量定义的好坏直接影响到模型规模和求解难度。我们面临一个网络G(V, E)和商品集合K。第一步为每条边上的每个商品流量创建二进制变量组。这不是一个容易的决定。最直接的想法是为每条边(i,j)和每个商品k分配一组变量来表示商品k在边(i,j)上的流量。但这样变量数将是|E| * |K| * (二进制位数)规模极易爆炸。我们观察赛题数据发现很多商品的路由路径是稀疏的即一种商品不太可能经过所有边。因此我们采用了基于路径预选的策略。预选潜在路径对于每个商品k使用经典的K最短路径算法如Yens Algorithm找出其起点s_k到终点t_k的Top M条最短路径M根据网络规模选取如10-20条。我们假设最优解很可能出现在这些较短的路径组合中。这一步极大地缩减了变量规模。变量定义对于商品k我们为其每一条预选路径p创建一个二值决策变量x_{kp}。x_{kp} 1表示商品k选择路径p来运输其全部需求量d_kx_{kp} 0则表示不选。同时我们约束每个商品必须且只能选择一条路径∑_p x_{kp} 1。这样我们巧妙地将连续流量分配问题转化为了离散路径选择问题。变量总数降至|K| * M个二值变量规模可控。第二步处理边容量约束。边e上的总流量是其上所有经过该边的商品流量之和。商品k若选择路径p且p包含边e则其对边e的流量贡献为d_k。因此边e的总流量f_e ∑_k ∑_{p: e∈p} d_k * x_{kp}。 容量约束为f_e ≤ C_e。我们引入非负松弛变量s_e同样用二进制展开表示假设容量C_e已知s_e的位数足以表示0到C_e的值将其转化为等式f_e s_e C_e。 对应的惩罚项为P_capacity * (f_e s_e - C_e)^2。 将f_e的表达式代入展开后我们会得到关于x_{kp}和s_e的二进制变量的二次型。3.2 Q矩阵的系统化构建方法有了变量和惩罚项表达式构建Q矩阵就是一个系统性的“会计”工作。我们编写了一个函数来自动化这个过程。假设我们有N个二值变量Q就是一个N x N的矩阵初始化为全零。目标函数项最小化总成本假设边(i,j)的单位流量成本为c_{ij}那么商品k选择路径p所产生的成本为d_k * (∑_{e∈p} c_e)。总成本就是∑_k ∑_p [ d_k * (∑_{e∈p} c_e) ] * x_{kp}。这是一个线性项。在QUBO中线性项c_i * x_i可以表示为c_i * x_i * x_i因为x_i^2 x_i对于二值变量成立。因此对于变量x_{kp}我们在Q矩阵的对角线位置Q[idx(k,p), idx(k,p)]上加上其成本系数。单路径选择约束惩罚项对于每个商品k约束∑_p x_{kp} 1。惩罚项为P_path * (∑_p x_{kp} - 1)^2。 展开P_path * [ (∑_p x_{kp})^2 - 2*(∑_p x_{kp}) 1 ]。(∑_p x_{kp})^2 ∑_p x_{kp}^2 ∑_{p≠q} 2 * x_{kp} * x_{kq}。由于x_{kp}^2 x_{kp}这部分会产生线性项和二次交叉项。-2*(∑_p x_{kp})产生线性项。1是常数不影响优化可以忽略。 我们需要将展开后每一项的系数累加到Q矩阵的对应位置。例如对于线性项(P_path * 1 - 2*P_path) * x_{kp}将其加到Q[idx(k,p), idx(k,p)]。对于二次项2*P_path * x_{kp} * x_{kq}将其加到Q[idx(k,p), idx(k,q)]和Q[idx(k,q), idx(k,p)]因为Q通常视为上三角只需加一次但需注意对称性。边容量约束惩罚项惩罚项P_capacity * (f_e s_e - C_e)^2。f_e是线性组合∑_k ∑_{p: e∈p} d_k * x_{kp}。 展开后将包含关于x_{kp}的二次项来自f_e^2。关于x_{kp}和s_e的二次交叉项来自2 * f_e * s_e。关于s_e的二次项和线性项。关于x_{kp}的线性项来自-2 * C_e * f_e和2 * f_e * s_e的一部分。常数项。每一项都需要根据其系数精确地累加到Q矩阵中对应变量索引的位置上。这个过程非常繁琐必须仔细处理下标和系数。实操心得我们最初手动推导公式很容易出错。后来我们改用符号计算工具如Python的sympy库来辅助。先定义符号变量然后让sympy展开惩罚项表达式再编程解析展开后的多项式结构自动提取每个单项式如x_i,x_i*x_j的系数最后映射到Q矩阵。这大大减少了人为错误并且当问题规模或约束形式变化时代码更容易调整。3.3 代码实现骨架与关键片段以下是构建QUBO模型核心部分的Python伪代码/骨架使用了dimod库D-Wave开源工具包的一部分来帮助管理QUBO模型但求解部分我们先使用经典模拟退火。import numpy as np import dimod from collections import defaultdict import itertools class NetworkFlowToQUBO: def __init__(self, graph, commodities, paths_for_commodity, edge_capacities, edge_costs, P_path, P_capacity): graph: 网络图节点和边列表 commodities: 商品列表每个商品是字典 {id: k, demand: d_k, origin: s_k, dest: t_k} paths_for_commodity: 字典key为商品idvalue为该商品预选的路径列表每条路径是边的序列 edge_capacities: 字典key为边value为容量C_e edge_costs: 字典key为边value为单位成本c_e P_path, P_capacity: 惩罚系数 self.graph graph self.commodities commodities self.paths paths_for_commodity self.C edge_capacities self.c edge_costs self.P_path P_path self.P_capacity P_capacity # 变量映射为每个决策变量分配一个唯一的索引 self.var_index {} index 0 # 1. 路径选择变量 x_{kp} for k in commodities: for p_idx, path in enumerate(self.paths[k[id]]): self.var_index[(‘x’, k[‘id’], p_idx)] index index 1 # 2. 边容量松弛变量 s_e (二进制展开) self.slack_vars {} for e in edge_capacities: # 计算表示0到C_e所需的二进制位数 bits_needed int(np.ceil(np.log2(edge_capacities[e] 1))) for b in range(bits_needed): self.var_index[(‘s’, e, b)] index self.slack_vars.setdefault(e, []).append(index) index 1 self.num_vars index self.Q np.zeros((self.num_vars, self.num_vars)) def add_objective(self): 添加总成本目标项 for k in self.commodities: demand k[‘demand’] for p_idx, path in enumerate(self.paths[k[‘id’]]): # 计算路径p的总成本 path_cost sum(self.c[e] for e in path) cost_coeff demand * path_cost var_idx self.var_index[(‘x’, k[‘id’], p_idx)] # 线性项 c*x 体现在对角线上 self.Q[var_idx, var_idx] cost_coeff def add_single_path_constraints(self): 添加每个商品必须选且仅选一条路径的约束惩罚项 for k in self.commodities: k_id k[‘id’] path_vars [self.var_index[(‘x’, k_id, p_idx)] for p_idx in range(len(self.paths[k_id]))] n_paths len(path_vars) # 惩罚项 P * (sum(x) - 1)^2 # 展开后的线性部分P * (1 - 2) * x_i -P * x_i for i in path_vars: self.Q[i, i] -1.0 * self.P_path # 展开后的二次部分2P * x_i * x_j for i_idx, i in enumerate(path_vars): for j_idx, j in enumerate(path_vars): if i_idx j_idx: # 上三角部分 self.Q[i, j] 2.0 * self.P_path def add_capacity_constraints(self): 添加边容量约束惩罚项 for e in self.C: capacity self.C[e] # 收集所有经过边e的路径变量及其贡献的流量 contributing_vars [] # 存储 (变量索引, 流量贡献d_k) for k in self.commodities: demand k[‘demand’] for p_idx, path in enumerate(self.paths[k[‘id’]]): if e in path: var_idx self.var_index[(‘x’, k[‘id’], p_idx)] contributing_vars.append((var_idx, demand)) # 获取该边松弛变量的索引 slack_indices self.slack_vars.get(e, []) # 计算松弛变量代表的数值s sum(2^b * s_b) # 我们需要构建表达式 f_e s - C_e然后计算其平方的系数 # 这里简化处理采用循环计算各项系数对于大规模问题可用sympy展开 # 构建变量列表和系数列表用于表示 (f_e s - C_e) all_vars [] coeffs [] # f_e 部分 for var_idx, demand in contributing_vars: all_vars.append(var_idx) coeffs.append(demand) # s 部分 for bit_pos, s_idx in enumerate(slack_indices): all_vars.append(s_idx) coeffs.append(2**bit_pos) # 二进制权重 # 常数部分 -C_e const_coeff -capacity # 计算 (sum(coeff_i * x_i) const)^2 的系数并累加到Q # 遍历所有变量对包括自身 n_vars_in_expr len(all_vars) for i in range(n_vars_in_expr): idx_i all_vars[i] coeff_i coeffs[i] # 二次项系数: P * coeff_i * coeff_i self.Q[idx_i, idx_i] self.P_capacity * coeff_i * coeff_i # 交叉项系数: 2 * P * coeff_i * coeff_j for j in range(i1, n_vars_in_expr): idx_j all_vars[j] coeff_j coeffs[j] self.Q[idx_i, idx_j] 2.0 * self.P_capacity * coeff_i * coeff_j # 与常数项的交叉产生线性项: 2 * P * coeff_i * const_coeff self.Q[idx_i, idx_i] 2.0 * self.P_capacity * coeff_i * const_coeff # 常数项平方不影响优化忽略 def get_qubo_matrix(self): 返回完整的QUBO矩阵上三角形式 # 确保矩阵是对称的dimod等库通常接受上三角 Q_upper np.triu(self.Q, k0) # 下三角部分可以通过对称性得到但很多求解器只需要上三角 return Q_upper def solve_with_simulated_annealing(self, num_reads1000): 使用经典模拟退火求解QUBO替代量子退火模拟 import neal # 获取QUBO矩阵 qubo self.get_qubo_matrix() # 将矩阵转换为dimod标准的字典格式 {(i,j): coeff} qubo_dict {} n self.num_vars for i in range(n): for j in range(i, n): coeff qubo[i, j] if abs(coeff) 1e-10: # 忽略接近零的系数 qubo_dict[(i, j)] coeff bqm dimod.BinaryQuadraticModel.from_qubo(qubo_dict) sampler neal.SimulatedAnnealingSampler() sampleset sampler.sample(bqm, num_readsnum_reads) # 取能量最低的解 best_sample sampleset.first.sample best_energy sampleset.first.energy return best_sample, best_energy, sampleset4. 求解策略、调参与结果分析4.1 经典求解器作为验证基准在真正尝试量子求解或模拟之前我们先用经典的数学规划求解器如Gurobi求解了原始的MILP模型。这有两个目的一是验证我们QUBO转化逻辑的正确性。如果我们构建的QUBO模型其最优解在惩罚系数足够大的情况下经过解码后与Gurobi求出的最优解一致或接近就说明我们的建模是准确的。二是为我们调整QUBO惩罚系数和评估求解质量提供一个“黄金标准”。对于中小规模算例Gurobi可以在可接受时间内求得精确最优解。我们将其作为基准对比QUBO模拟退火得到的结果。4.2 模拟退火求解与参数调优由于比赛期间无法使用真实的量子退火机我们使用经典模拟退火SA来模拟量子退火的过程。我们选择了D-Wave的neal库。模拟退火的效果高度依赖于退火计划初始温度、终止温度、降温速率、迭代次数等。关键调参经验退火计划我们采用了指数降温策略。初始温度T0需要设置得足够高以使算法在初期有足够概率接受差解跳出局部最优。一个经验法则是T0应该大于目标函数值典型差异的量级。我们通过几次试跑观察能量下降曲线来调整T0和num_sweeps总迭代步数。采样次数num_reads模拟退火是随机算法单次运行可能陷入不同的局部最优。我们设置num_reads1000甚至更高从大量独立退火过程中选取能量最低的解。这模仿了量子退火机多次读取multiple reads的操作。惩罚系数P的精细调整这是最耗时的部分。我们编写了一个自动化脚本固定P_path为一个较大的值确保路径选择约束必须满足。让P_capacity在一个范围内变化例如从1, 10, 100, 1000, ...。对每个P_capacity运行模拟退火获取若干最优解。解码这些解将二值变量x_{kp}还原为路径选择将松弛变量s_e还原为剩余容量。检查约束满足情况是否所有商品都只选了一条路径是否所有边流量都不超容s_e 0记录可行解的比例和平均目标成本。 我们发现当P_capacity太小时得到的解几乎总是违反容量约束。随着P_capacity增大可行解比例上升但目标成本也可能变差因为惩罚项主导了能量函数迫使求解器过于关注满足约束而忽略了原始成本。我们需要选择一个能稳定产生高比例可行解且成本接近经典最优解的P_capacity值。4.3 结果解码与方案验证从模拟退火得到的最佳样本一个二值变量字典后我们需要将其“解码”成可读的调度方案。def decode_solution(self, sample): 将求解器返回的样本解码为路径分配方案和边负载情况 solution {} flow_on_edge defaultdict(float) # 1. 解码路径选择 for k in self.commodities: k_id k[‘id’] selected_path_idx None for p_idx in range(len(self.paths[k_id])): var_key (‘x’, k_id, p_idx) if var_key in self.var_index and sample.get(self.var_index[var_key], 0) 0.5: # 视为1 if selected_path_idx is not None: print(f警告商品 {k_id} 选择了多条路径) selected_path_idx p_idx if selected_path_idx is not None: path self.paths[k_id][selected_path_idx] solution[k_id] path # 累加边流量 demand k[‘demand’] for e in path: flow_on_edge[e] demand else: print(f警告商品 {k_id} 未选择任何路径) solution[k_id] None # 2. 解码松弛变量验证容量 capacity_violation {} for e in self.C: slack_value 0 if e in self.slack_vars: for bit_pos, s_idx in enumerate(self.slack_vars[e]): if sample.get(s_idx, 0) 0.5: slack_value 2**bit_pos used_capacity flow_on_edge.get(e, 0) remaining_capacity self.C[e] - slack_value # 理论上应有used_capacity slack_value self.C[e] if abs(used_capacity slack_value - self.C[e]) 1e-5: print(f边 {e} 容量约束惩罚项可能未完全满足使用量{used_capacity}, 松弛量{slack_value}, 容量{self.C[e]}) capacity_violation[e] max(0, used_capacity - self.C[e]) # 实际超容量 return solution, flow_on_edge, capacity_violation解码后我们对比了QUBOSA方案与Gurobi精确解。在调优后的参数下对于中等规模问题我们方案的总成本与最优解的差距可以控制在5%以内且所有约束都得到满足。对于更大规模的问题Gurobi可能无法在时限内求得最优解而我们的启发式方法仍能在较短时间内给出一个质量不错的可行解这体现了其应用潜力。4.4 向真实量子硬件的延伸思考虽然比赛用的是经典模拟但我们完整走通了QUBO建模流程这为未来使用真实量子退火机如D-Wave奠定了基础。要上真机还需考虑嵌入问题QUBO模型的变量需要映射到量子退火机的物理量子比特上而量子比特之间存在连接限制。D-Wave的Chimera或Pegasus图结构并非全连接因此复杂的QUBO交互项可能需要通过“链”将多个物理量子比特耦合起来代表一个逻辑变量。这需要额外的步骤minor-embedding可能会引入错误并消耗更多量子比特资源。问题规模当前量子退火机的量子比特数数千个仍然有限能直接求解的QUBO问题规模变量数受限于其图结构。对于大规模问题需要结合经典算法进行分解例如使用量子退火作为子问题求解器。参数调优真实量子退火有退火时间、退火路径等更多参数需要调节且存在噪声。这需要更复杂的调参和误差缓解技术。5. 常见问题、挑战与应对策略在整个项目过程中我们遇到了不少坑这里总结一下希望能帮你避雷。1. 问题规模爆炸与变量定义策略挑战最直接的变量定义每条边*每个商品会导致变量数急剧膨胀使QUBO矩阵大到无法处理。应对采用基于路径的变量定义是降维的关键。通过预选K短路径将连续流量分配转化为离散路径选择极大减少了变量数。关键在于预选路径的数量M要取得合理太少可能丢失最优解太多则失去降维意义。需要根据网络直径和稀疏性进行试验。2. 惩罚系数P的设定难题挑战P太小约束无效P太大数值问题突出矩阵元素过大且可能使能量地形过于陡峭让求解器无论是经典还是量子难以搜索。应对分层调整对不同约束使用不同的惩罚系数。例如单路径选择约束P_path通常需要设定为绝对优先满足的硬约束可以设一个很大的固定值。容量约束的P_capacity则需要精细调节。自适应方法实现一个简单的循环从一个较小的P开始求解检查解是否可行如果不可行增大P再次求解直到得到可行解。或者可以参考一些文献中的公式根据约束违反的期望程度来估算P的下界。归一化在构建QUBO前对原始目标函数的成本系数进行归一化例如缩放到[0,1]区间这样惩罚系数的选取范围会更稳定。3. 经典求解器与QUBO求解结果对比时的“差距”挑战即使约束都满足了QUBOSA解的成本可能仍高于Gurobi求出的最优解。应对首先要明白模拟退火是启发式算法不保证最优。其次QUBO模型本身是原问题的近似例如路径预选可能排除了真正的最优路径。差距在可接受范围内如5%-10%是合理的。评估时应更关注QUBO方法在大规模问题或实时求解场景下的可用性而不是在小规模问题上与精确解拼小数点。4. 松弛变量的二进制展开位数选择挑战边容量C_e可能很大用二进制展开需要很多位增加变量。应对分析数据如果很多边容量远大于实际可能的总流量可以适当减少位数。或者可以采用对数编码等更高效的整数编码方式但这会增加QUBO的复杂度。一个稳妥的做法是位数b满足2^b C_e的最小整数。5. 模拟退火耗时与解的质量权衡挑战增加num_reads和num_sweeps可以提高找到好解的概率但也会增加计算时间。应对不要盲目追求最大迭代次数。可以先做参数扫描固定num_reads观察解的能量随num_sweeps的变化曲线在曲线进入平台区的位置选择一个性价比高的值。同样可以调整num_reads观察多次独立运行最佳能量的分布当分布稳定时即可。6. 代码实现中的精度与效率问题挑战构建大型QUBO矩阵时双重循环容易导致效率低下且浮点数累加可能产生精度误差。应对使用稀疏矩阵格式如scipy.sparse存储Q因为大多数变量之间没有交互Q矩阵非常稀疏。如前所述利用符号计算sympy自动生成系数映射避免手动推导和编码错误。对于系数累加使用np.add.at或类似方法进行无缓冲区的累加避免精度损失。这次MathorCup A题的实战让我深刻体会到将现实问题转化为QUBO模型是一套系统性的“翻译”功夫需要同时对原问题领域和QUBO建模技巧有深入理解。它不是一个黑箱工具而是一种需要精心设计的建模语言。虽然过程充满挑战但看到经典优化问题能够通过这种形式与前沿的量子计算硬件对话并为未来可能带来的加速充满期待这一切的付出都是值得的。对于想入门量子优化的同学我的建议是从小规模的、经典的NP-Hard问题如最大割、旅行商问题开始练习QUBO建模使用dimod和模拟退火进行验证然后再挑战像网络流这样带有复杂约束的实际问题。这条路很长但每一步都算数。 1. 项目概述从数学建模到量子计算前沿的实战跨越去年带队参加MathorCup选了A题最后拿了个一等奖。赛题的核心是把一个经典的组合优化问题——大规模网络中的资源分配与路径规划转化成了量子计算领域里一个非常热门的模型QUBO。这听起来有点跨界数学建模怎么就和量子计算扯上关系了其实这正是当前交叉学科研究的一个缩影。我们队伍里既有搞传统运筹优化的也有对量子计算感兴趣的大家一合计觉得这个方向既有挑战性未来潜力也大就决定啃下这块硬骨头。简单来说QUBO模型全称是“二次无约束二值优化”它之所以在量子计算圈特别是量子退火和相干伊辛机这类专用硬件上备受青睐是因为它的形式极其简洁所有变量只能是0或1目标函数是一个二次型。很多NP难问题比如旅行商问题、图着色、设备布局等都能被“编码”成QUBO形式。一旦编码成功理论上就可以丢给量子退火机去求解利用量子隧穿、量子纠缠等效应在特定问题上寻找比经典算法更优的解。我们这次做的就是把一个实际的、带有多重约束的数学建模问题“翻译”成QUBO这门“量子硬件能听懂的语言”。整个过程充满了探索和试错。网上关于QUBO建模的教程大多停留在理论和小型玩具案例真正把一个中等规模、带有复杂现实约束的MathorCup赛题完整走通的分享并不多。我们花了大量时间在“建模转化”和“求解策略”上如何精确地将约束条件“惩罚”进目标函数如何调整惩罚权重避免无效解以及在没有真实量子硬件的情况下如何用经典模拟器和启发式算法来验证和求解我们构建的QUBO模型。这篇文章我就把这套从问题分析、QUBO建模、到经典求解与结果分析的完整思路和核心代码分享出来希望能给未来想涉足量子优化或者参加类似交叉学科竞赛的朋友们一些实实在在的参考。2. 赛题核心与QUBO建模思路拆解2.1 问题重述与经典建模视角我们遇到的A题本质上是一个多商品流网络优化问题。场景可以抽象为一个有向图节点代表不同的资源点或任务点边代表连接路径每条边有容量限制。有多种不同的“商品”或任务需要在网络中从各自的起点运输到终点每种商品有确定的流量需求。目标是在满足所有边容量约束的前提下找到一种路由方案使得总运输成本或距离、时间等最小。如果完全用经典方法比如混合整数线性规划来建模我们会定义一系列决策变量x_{ijk}表示商品k是否经过边(i, j)流量是多少。然后列出流量守恒约束、容量约束、变量类型约束构建一个大规模的MILP模型丢给CPLEX、Gurobi这类求解器去算。对于规模不大的问题这很有效。但一旦节点和商品数量上去问题规模会指数级膨胀求解时间可能变得不可接受。这时量子计算特别是量子退火提供了一种不同的思路。它不擅长直接处理复杂的线性约束但它非常擅长在巨大的解空间中进行“能量”最小化的搜索。QUBO模型正是描述这种“能量景观”的完美数学形式。我们的核心任务就是将带有约束的经典优化问题转化为一个无约束的二次二值优化问题。转化的关键在于把所有的硬约束必须满足的条件都转换成“惩罚项”加到目标函数里。如果某个解违反了约束它的“能量”即目标函数值就会因为惩罚项而急剧增高从而在最小化过程中被淘汰。2.2 QUBO模型的核心形式与转化哲学QUBO模型的标准形式如下Minimize: y x^T Q x其中x是一个由二值变量0或1组成的列向量Q是一个实对称矩阵或上三角矩阵。我们的所有智慧都体现在如何构造这个Q矩阵上。对于我们的网络流问题直接使用0/1变量表示“是否选择某条边”是自然的。但难点在于流量是连续值而QUBO变量是离散二值的。这里就需要引入二进制展开的技巧。例如如果一条边的容量上限是7我们可以用3个量子比特对应3个二值变量的二进制组合来表示这条边上分配的流量0到7。假设变量x1, x2, x3分别代表2^0, 2^1, 2^2位那么流量f x1 2*x2 4*x3。这样我们就用一组二值变量描述了一个有限范围内的整数量。接下来是重头戏约束的惩罚项构造。以边容量约束为例对于边e其上的总流量f_e不能超过容量C_e。我们可以将其转化为惩罚项P * (f_e - C_e)^2但这里f_e - C_e可能为负未超容平方后仍为正这会错误地惩罚未超容的好解。因此更标准的做法是引入松弛变量。定义s_e为一个非负整数同样用二进制展开表示将约束改写为f_e s_e C_e。那么惩罚项就是P * (f_e s_e - C_e)^2。当约束被严格满足时此项为0否则为正起到惩罚作用。惩罚系数P是一个需要精心调整的超参数它必须足够大以确保任何违反约束的解的目标值都差于任何可行解但又不能太大以免数值溢出或导致能量地形过于陡峭难以优化。流量守恒约束对于每个商品在每个中间节点流入等于流出是另一个核心约束处理方式类似也需要为每个节点和每个商品引入类似的等式约束惩罚项。最终我们的目标函数将由两部分组成原始的最小化总成本项通常是流量的线性函数通过二进制展开后转化为二值变量的二次型和所有约束对应的惩罚项。将它们全部加起来就得到了最终的QUBO形式x^T Q x。构建Q矩阵的过程就是将这些二次项和一次项可以转化为二次项如x_i可视为x_i * x_i的系数整理到一个矩阵里的过程。注意惩罚系数P的选取是艺术也是科学。一个实用的技巧是让其大小与原始目标函数值的量级相匹配。可以先忽略惩罚项求解一个松弛问题估算目标值范围然后设定P比这个范围大一个数量级例如10倍再进行调试。我们当时采用了迭代调整策略从一个小值开始逐步增加直到求得的解100%满足所有约束。3. 模型构建的详细步骤与关键实现3.1 决策变量定义与二进制编码这是整个建模的基石变量定义的好坏直接影响到模型规模和求解难度。我们面临一个网络G(V, E)和商品集合K。第一步为每条边上的每个商品流量创建二进制变量组。这不是一个容易的决定。最直接的想法是为每条边(i,j)和每个商品k分配一组变量来表示商品k在边(i,j)上的流量。但这样变量数将是|E| * |K| * (二进制位数)规模极易爆炸。我们观察赛题数据发现很多商品的路由路径是稀疏的即一种商品不太可能经过所有边。因此我们采用了基于路径预选的策略。预选潜在路径对于每个商品k使用经典的K最短路径算法如Yens Algorithm找出其起点s_k到终点t_k的Top M条最短路径M根据网络规模选取如10-20条。我们假设最优解很可能出现在这些较短的路径组合中。这一步极大地缩减了变量规模。变量定义对于商品k我们为其每一条预选路径p创建一个二值决策变量x_{kp}。x_{kp} 1表示商品k选择路径p来运输其全部需求量d_kx_{kp} 0则表示不选。同时我们约束每个商品必须且只能选择一条路径∑_p x_{kp} 1。这样我们巧妙地将连续流量分配问题转化为了离散路径选择问题。变量总数降至|K| * M个二值变量规模可控。第二步处理边容量约束。边e上的总流量是其上所有经过该边的商品流量之和。商品k若选择路径p且p包含边e则其对边e的流量贡献为d_k。因此边e的总流量f_e ∑_k ∑_{p: e∈p} d_k * x_{kp}。 容量约束为f_e ≤ C_e。我们引入非负松弛变量s_e同样用二进制展开表示假设容量C_e已知s_e的位数足以表示0到C_e的值将其转化为等式f_e s_e C_e。 对应的惩罚项为P_capacity * (f_e s_e - C_e)^2。 将f_e的表达式代入展开后我们会得到关于x_{kp}和s_e的二进制变量的二次型。3.2 Q矩阵的系统化构建方法有了变量和惩罚项表达式构建Q矩阵就是一个系统性的“会计”工作。我们编写了一个函数来自动化这个过程。假设我们有N个二值变量Q就是一个N x N的矩阵初始化为全零。目标函数项最小化总成本假设边(i,j)的单位流量成本为c_{ij}那么商品k选择路径p所产生的成本为d_k * (∑_{e∈p} c_e)。总成本就是∑_k ∑_p [ d_k * (∑_{e∈p} c_e) ] * x_{kp}。这是一个线性项。在QUBO中线性项c_i * x_i可以表示为c_i * x_i * x_i因为x_i^2 x_i对于二值变量成立。因此对于变量x_{kp}我们在Q矩阵的对角线位置Q[idx(k,p), idx(k,p)]上加上其成本系数。单路径选择约束惩罚项对于每个商品k约束∑_p x_{kp} 1。惩罚项为P_path * (∑_p x_{kp} - 1)^2。 展开P_path * [ (∑_p x_{kp})^2 - 2*(∑_p x_{kp}) 1 ]。(∑_p x_{kp})^2 ∑_p x_{kp}^2 ∑_{p≠q} 2 * x_{kp} * x_{kq}。由于x_{kp}^2 x_{kp}这部分会产生线性项和二次交叉项。-2*(∑_p x_{kp})产生线性项。1是常数不影响优化可以忽略。 我们需要将展开后每一项的系数累加到Q矩阵的对应位置。例如对于线性项(P_path * 1 - 2*P_path) * x_{kp}将其加到Q[idx(k,p), idx(k,p)]。对于二次项2*P_path * x_{kp} * x_{kq}将其加到Q[idx(k,p), idx(k,q)]和Q[idx(k,q), idx(k,p)]因为Q通常视为上三角只需加一次但需注意对称性。边容量约束惩罚项惩罚项P_capacity * (f_e s_e - C_e)^2。f_e是线性组合∑_k ∑_{p: e∈p} d_k * x_{kp}。 展开后将包含关于x_{kp}的二次项来自f_e^2。关于x_{kp}和s_e的二次交叉项来自2 * f_e * s_e。关于s_e的二次项和线性项。关于x_{kp}的线性项来自-2 * C_e * f_e和2 * f_e * s_e的一部分。常数项。每一项都需要根据其系数精确地累加到Q矩阵中对应变量索引的位置上。这个过程非常繁琐必须仔细处理下标和系数。实操心得我们最初手动推导公式很容易出错。后来我们改用符号计算工具如Python的sympy库来辅助。先定义符号变量然后让sympy展开惩罚项表达式再编程解析展开后的多项式结构自动提取每个单项式如x_i,x_i*x_j的系数最后映射到Q矩阵。这大大减少了人为错误并且当问题规模或约束形式变化时代码更容易调整。3.3 代码实现骨架与关键片段以下是构建QUBO模型核心部分的Python伪代码/骨架使用了dimod库D-Wave开源工具包的一部分来帮助管理QUBO模型但求解部分我们先使用经典模拟退火。import numpy as np import dimod from collections import defaultdict import itertools class NetworkFlowToQUBO: def __init__(self, graph, commodities, paths_for_commodity, edge_capacities, edge_costs, P_path, P_capacity): graph: 网络图节点和边列表 commodities: 商品列表每个商品是字典 {id: k, demand: d_k, origin: s_k, dest: t_k} paths_for_commodity: 字典key为商品idvalue为该商品预选的路径列表每条路径是边的序列 edge_capacities: 字典key为边value为容量C_e edge_costs: 字典key为边value为单位成本c_e P_path, P_capacity: 惩罚系数 self.graph graph self.commodities commodities self.paths paths_for_commodity self.C edge_capacities self.c edge_costs self.P_path P_path self.P_capacity P_capacity # 变量映射为每个决策变量分配一个唯一的索引 self.var_index {} index 0 # 1. 路径选择变量 x_{kp} for k in commodities: for p_idx, path in enumerate(self.paths[k[id]]): self.var_index[(‘x’, k[‘id’], p_idx)] index index 1 # 2. 边容量松弛变量 s_e (二进制展开) self.slack_vars {} for e in edge_capacities: # 计算表示0到C_e所需的二进制位数 bits_needed int(np.ceil(np.log2(edge_capacities[e] 1))) for b in range(bits_needed): self.var_index[(‘s’, e, b)] index self.slack_vars.setdefault(e, []).append(index) index 1 self.num_vars index self.Q np.zeros((self.num_vars, self.num_vars)) def add_objective(self): 添加总成本目标项 for k in self.commodities: demand k[‘demand’] for p_idx, path in enumerate(self.paths[k[‘id’]]): # 计算路径p的总成本 path_cost sum(self.c[e] for e in path) cost_coeff demand * path_cost var_idx self.var_index[(‘x’, k[‘id’], p_idx)] # 线性项 c*x 体现在对角线上 self.Q[var_idx, var_idx] cost_coeff def add_single_path_constraints(self): 添加每个商品必须选且仅选一条路径的约束惩罚项 for k in self.commodities: k_id k[‘id’] path_vars [self.var_index[(‘x’, k_id, p_idx)] for p_idx in range(len(self.paths[k_id]))] n_paths len(path_vars) # 惩罚项 P * (sum(x) - 1)^2 # 展开后的线性部分P * (1 - 2) * x_i -P * x_i for i in path_vars: self.Q[i, i] -1.0 * self.P_path # 展开后的二次部分2P * x_i * x_j for i_idx, i in enumerate(path_vars): for j_idx, j in enumerate(path_vars): if i_idx j_idx: # 上三角部分 self.Q[i, j] 2.0 * self.P_path def add_capacity_constraints(self): 添加边容量约束惩罚项 for e in self.C: capacity self.C[e] # 收集所有经过边e的路径变量及其贡献的流量 contributing_vars [] # 存储 (变量索引, 流量贡献d_k) for k in self.commodities: demand k[‘demand’] for p_idx, path in enumerate(self.paths[k[‘id’]]): if e in path: var_idx self.var_index[(‘x’, k[‘id’], p_idx)] contributing_vars.append((var_idx, demand)) # 获取该边松弛变量的索引 slack_indices self.slack_vars.get(e, []) # 计算松弛变量代表的数值s sum(2^b * s_b) # 我们需要构建表达式 f_e s - C_e然后计算其平方的系数 # 这里简化处理采用循环计算各项系数对于大规模问题可用sympy展开 # 构建变量列表和系数列表用于表示 (f_e s - C_e) all_vars [] coeffs [] # f_e 部分 for var_idx, demand in contributing_vars: all_vars.append(var_idx) coeffs.append(demand) # s 部分 for bit_pos, s_idx in enumerate(slack_indices): all_vars.append(s_idx) coeffs.append(2**bit_pos) # 二进制权重 # 常数部分 -C_e const_coeff -capacity # 计算 (sum(coeff_i * x_i) const)^2 的系数并累加到Q # 遍历所有变量对包括自身 n_vars_in_expr len(all_vars) for i in range(n_vars_in_expr): idx_i all_vars[i] coeff_i coeffs[i] # 二次项系数: P * coeff_i * coeff_i self.Q[idx_i, idx_i] self.P_capacity * coeff_i * coeff_i # 交叉项系数: 2 * P * coeff_i * coeff_j for j in range(i1, n_vars_in_expr): idx_j all_vars[j] coeff_j coeffs[j] self.Q[idx_i, idx_j] 2.0 * self.P_capacity * coeff_i * coeff_j # 与常数项的交叉产生线性项: 2 * P * coeff_i * const_coeff self.Q[idx_i, idx_i] 2.0 * self.P_capacity * coeff_i * const_coeff # 常数项平方不影响优化忽略 def get_qubo_matrix(self): 返回完整的QUBO矩阵上三角形式 # 确保矩阵是对称的dimod等库通常接受上三角 Q_upper np.triu(self.Q, k0) # 下三角部分可以通过对称性得到但很多求解器只需要上三角 return Q_upper def solve_with_simulated_annealing(self, num_reads1000): 使用经典模拟退火求解QUBO替代量子退火模拟 import neal # 获取QUBO矩阵 qubo self.get_qubo_matrix() # 将矩阵转换为dimod标准的字典格式 {(i,j): coeff} qubo_dict {} n self.num_vars for i in range(n): for j in range(i, n): coeff qubo[i, j] if abs(coeff) 1e-10: # 忽略接近零的系数 qubo_dict[(i, j)] coeff bqm dimod.BinaryQuadraticModel.from_qubo(qubo_dict) sampler neal.SimulatedAnnealingSampler() sampleset sampler.sample(bqm, num_readsnum_reads) # 取能量最低的解 best_sample sampleset.first.sample best_energy sampleset.first.energy return best_sample, best_energy, sampleset4. 求解策略、调参与结果分析4.1 经典求解器作为验证基准在真正尝试量子求解或模拟之前我们先用经典的数学规划求解器如Gurobi求解了原始的MILP模型。这有两个目的一是验证我们QUBO转化逻辑的正确性。如果我们构建的QUBO模型其最优解在惩罚系数足够大的情况下经过解码后与Gurobi求出的最优解一致或接近就说明我们的建模是准确的。二是为我们调整QUBO惩罚系数和评估求解质量提供一个“黄金标准”。对于中小规模算例Gurobi可以在可接受时间内求得精确最优解。我们将其作为基准对比QUBO模拟退火得到的结果。4.2 模拟退火求解与参数调优由于比赛期间无法使用真实的量子退火机我们使用经典模拟退火SA来模拟量子退火的过程。我们选择了D-Wave的neal库。模拟退火的效果高度依赖于退火计划初始温度、终止温度、降温速率、迭代次数等。关键调参经验退火计划我们采用了指数降温策略。初始温度T0需要设置得足够高以使算法在初期有足够概率接受差解跳出局部最优。一个经验法则是T0应该大于目标函数值典型差异的量级。我们通过几次试跑观察能量下降曲线来调整T0和num_sweeps总迭代步数。采样次数num_reads模拟退火是随机算法单次运行可能陷入不同的局部最优。我们设置num_reads1000甚至更高从大量独立退火过程中选取能量最低的解。这模仿了量子退火机多次读取multiple reads的操作。惩罚系数P的精细调整这是最耗时的部分。我们编写了一个自动化脚本固定P_path为一个较大的值确保路径选择约束必须满足。让P_capacity在一个范围内变化例如从1, 10, 100, 1000, ...。对每个P_capacity运行模拟退火获取若干最优解。解码这些解将二值变量x_{kp}还原为路径选择将松弛变量s_e还原为剩余容量。检查约束满足情况是否所有商品都只选了一条路径是否所有边流量都不超容s_e 0记录可行解的比例和平均目标成本。 我们发现当P_capacity太小时得到的解几乎总是违反容量约束。随着P_capacity增大可行解比例上升但目标成本也可能变差因为惩罚项主导了能量函数迫使求解器过于关注满足约束而忽略了原始成本。我们需要选择一个能稳定产生高比例可行解且成本接近经典最优解的P_capacity值。4.3 结果解码与方案验证从模拟退火得到的最佳样本一个二值变量字典后我们需要将其“解码”成可读的调度方案。def decode_solution(self, sample): 将求解器返回的样本解码为路径分配方案和边负载情况 solution {} flow_on_edge defaultdict(float) # 1. 解码路径选择 for k in self.commodities: k_id k[‘id’] selected_path_idx None for p_idx in range(len(self.paths[k_id])): var_key (‘x’, k_id, p_idx) if var_key in self.var_index and sample.get(self.var_index[var_key], 0) 0.5: # 视为1 if selected_path_idx is not None: print(f警告商品 {k_id} 选择了多条路径) selected_path_idx p_idx if selected_path_idx is not None: path self.paths[k_id][selected_path_idx] solution[k_id] path # 累加边流量 demand k[‘demand’] for e in path: flow_on_edge[e] demand else: print(f警告商品 {k_id} 未选择任何路径) solution[k_id] None # 2. 解码松弛变量验证容量 capacity_violation {} for e in self.C: slack_value 0 if e in self.slack_vars: for bit_pos, s_idx in enumerate(self.slack_vars[e]): if sample.get(s_idx, 0) 0.5: slack_value 2**bit_pos used_capacity flow_on_edge.get(e, 0) remaining_capacity self.C[e] - slack_value # 理论上应有used_capacity slack_value self.C[e] if abs(used_capacity slack_value - self.C[e]) 1e-5: print(f边 {e} 容量约束惩罚项可能未完全满足使用量{used_capacity}, 松弛量{slack_value}, 容量{self.C[e]}) capacity_violation[e] max(0, used_capacity - self.C[e]) # 实际超容量 return solution, flow_on_edge, capacity_violation解码后我们对比了QUBOSA方案与Gurobi精确解。在调优后的参数下对于中等规模问题我们方案的总成本与最优解的差距可以控制在5%以内且所有约束都得到满足。对于更大规模的问题Gurobi可能无法在时限内求得最优解而我们的启发式方法仍能在较短时间内给出一个质量不错的可行解这体现了其应用潜力。4.4 向真实量子硬件的延伸思考虽然比赛用的是经典模拟但我们完整走通了QUBO建模流程这为未来使用真实量子退火机如D-Wave奠定了基础。要上真机还需考虑嵌入问题QUBO模型的变量需要映射到量子退火机的物理量子比特上而量子比特之间存在连接限制。D-Wave的Chimera或Pegasus图结构并非全连接因此复杂的QUBO交互项可能需要通过“链”将多个物理量子比特耦合起来代表一个逻辑变量。这需要额外的步骤minor-embedding可能会引入错误并消耗更多量子比特资源。问题规模当前量子退火机的量子比特数数千个仍然有限能直接求解的QUBO问题规模变量数受限于其图结构。对于大规模问题需要结合经典算法进行分解例如使用量子退火作为子问题求解器。参数调优真实量子退火有退火时间、退火路径等更多参数需要调节且存在噪声。这需要更复杂的调参和误差缓解技术。5. 常见问题、挑战与应对策略在整个项目过程中我们遇到了不少坑这里总结一下希望能帮你避雷。1. 问题规模爆炸与变量定义策略挑战最直接的变量定义每条边*每个商品会导致变量数急剧膨胀使QUBO矩阵大到无法处理。应对采用基于路径的变量定义是降维的关键。通过预选K短路径将连续流量分配转化为离散路径选择极大减少了变量数。关键在于预选路径的数量M要取得合理太少可能丢失最优解太多则失去降维意义。需要根据网络直径和稀疏性进行试验。2. 惩罚系数P的设定难题挑战P太小约束无效P太大数值问题突出矩阵元素过大且可能使能量地形过于陡峭让求解器无论是经典还是量子难以搜索。应对分层调整对不同约束使用不同的惩罚系数。例如单路径选择约束P_path通常需要设定为绝对优先满足的硬约束可以设一个很大的固定值。容量约束的P_capacity则需要精细调节。自适应方法实现一个简单的循环从一个较小的P开始求解检查解是否可行如果不可行增大P再次求解直到得到可行解。或者可以参考一些文献中的公式根据约束违反的期望程度来估算P的下界。归一化在构建QUBO前对原始目标函数的成本系数进行归一化例如缩放到[0,1]区间这样惩罚系数的选取范围会更稳定。3. 经典求解器与QUBO求解结果对比时的“差距”挑战即使约束都满足了QUBOSA解的成本可能仍高于Gurobi求出的最优解。应对首先要明白模拟退火是启发式算法不保证最优。其次QUBO模型本身是原问题的近似例如路径预选可能排除了真正的最优路径。差距在可接受范围内如5%-10%是合理的。评估时应更关注QUBO方法在大规模问题或实时求解场景下的可用性而不是在小规模问题上与精确解拼小数点。4. 松弛变量的二进制展开位数选择挑战边容量C_e可能很大用二进制展开需要很多位增加变量。应对分析数据如果很多边容量远大于实际可能的总流量可以适当减少位数。或者可以采用对数编码等更高效的整数编码方式但这会增加QUBO的复杂度。一个稳妥的做法是位数b满足2^b C_e的最小整数。5. 模拟退火耗时与解的质量权衡挑战增加num_reads和num_sweeps可以提高找到好解的概率但也会增加计算时间。应对不要盲目追求最大迭代次数。可以先做参数扫描固定num_reads观察解的能量随num_sweeps的变化曲线在曲线进入平台区的位置选择一个性价比高的值。同样可以调整num_reads观察多次独立运行最佳能量的分布当分布稳定时即可。6. 代码实现中的精度与效率问题挑战构建大型QUBO矩阵时双重循环容易导致效率低下且浮点数累加可能产生精度误差。应对使用稀疏矩阵格式如scipy.sparse存储Q因为大多数变量之间没有交互Q矩阵非常稀疏。如前所述利用符号计算sympy自动生成系数映射避免手动推导和编码错误。对于系数累加使用np.add.at或类似方法进行无缓冲区的累加避免精度损失。这次MathorCup A题的实战让我深刻体会到将现实问题转化为QUBO模型是一套系统性的“翻译”功夫需要同时对原问题领域和QUBO建模技巧有深入理解。它不是一个黑箱工具而是一种需要精心设计的建模语言。虽然过程充满挑战但看到经典优化问题能够通过这种形式与前沿的量子计算硬件对话并为未来可能带来的加速充满期待这一切的付出都是值得的。对于想入门量子优化的同学我的建议是从小规模的、经典的NP-Hard问题如最大割、旅行商问题开始练习QUBO建模使用dimod和模拟退火进行验证然后再挑战像网络流这样带有复杂约束的实际问题。这条路很长但每一步都算数。 
