4. 基于 QUBO 的背包问题建模与求解

问题设定与 QUBO 目标函数

取一个 5 件物品的 0-1 背包实例：重量 w = [2, 3, 4, 5, 9]，价值 v = [3, 4, 5, 8, 10]，容量 C = 12。决策变量 x_i ∈ {0,1} 表示是否装入第 i 件物品。目标是最大化总价值，同时满足总重量不超过容量。

QUBO 的标准形式是 min xᵀQx（x_i ∈ {0,1}）。因此先把最大化取负，再把不等式约束转成惩罚项。引入松弛变量把 Σw_i x_i ≤ C 变成等式：

Σ w_i x_i + s = C,  s ∈ {0,1,...,C}

松弛变量用二进制展开 s = Σ 2^k y_k，需要 ⌈log2(C+1)⌉ = 4 个辅助比特 y_0..y_3（权重 1,2,4,8）。于是 QUBO 目标为：

H(x,y) = -Σ v_i x_i + P · (Σ w_i x_i + Σ 2^k y_k - C)²

其中 P 是惩罚系数。展开平方项后，对角项 Q_ii 来自 a_i² x_i² = a_i² x_i（因 x_i²=x_i），交叉项 Q_ij = 2 a_i a_j（i≠j），常数项 P·C² 不影响最优解位置，可丢弃。下面用代码构造 Q 矩阵并枚举验证。

构造 Q 矩阵并枚举验证

import itertools, numpy as np

w = np.array([2, 3, 4, 5, 9])
v = np.array([3, 4, 5, 8, 10])
C = 12
n = len(w)
# 松弛变量二进制权重：1,2,4,8
s_weights = [1, 2, 4, 8]
m = len(s_weights)
N = n + m  # 总比特数 = 9

def build_Q(P):
    Q = np.zeros((N, N))
    # 价值项：-v_i 放在对角
    for i in range(n):
        Q[i, i] += -v[i]
    # 惩罚项系数向量 a = [w_0..w_4, 1,2,4,8, -C]
    a = np.concatenate([w, s_weights, [-C]])
    # 展开 P*(Σ a_k z_k)^2，z 为 x 与 y 拼接
    for i in range(N):
        Q[i, i] += P * a[i] * a[i]
        for j in range(i+1, N):
            Q[i, j] += 2 * P * a[i] * a[j]
    return Q

def energy(Q, z):
    return float(z @ Q @ z)

def brute_force(Q):
    best, best_z = np.inf, None
    for bits in itertools.product([0,1], repeat=N):
        z = np.array(bits)
        e = energy(Q, z)
        if e < best:
            best, best_z = e, z
    return best, best_z

for P in [1, 2, 5, 10, 20]:
    Q = build_Q(P)
    best, z = brute_force(Q)
    xs = z[:n]
    weight = int(w @ xs)
    value = int(v @ xs)
    feasible = weight <= C
    print(f"P={P:2d}  bestE={best:8.2f}  x={xs}  w={weight:2d}  v={value:2d}  feasible={feasible}")

运行后可以观察到：当 P 较小时，惩罚不足以压制超重解，最优比特串可能违反容量约束；P 增大到某个阈值后，最优解稳定落在可行域内。这个阈值与价值项和约束项的尺度比有关——本例中单件最大价值为 10，而超重 1 单位带来的惩罚增量约为 P·(2C+1) 量级，因此 P 需要明显大于价值量级才能保证可行性。

用退火/QAOA 求解并验证惩罚系数影响

枚举只能处理小规模。实际求解时把 Q 矩阵交给量子退火器（D-Wave 的 dimod + neal 模拟退火）或 QAOA 模拟器。下面用 neal 的模拟退火采样器对每个 P 采样，统计可行解比例，验证惩罚项的作用。

import neal, dimod

def solve_with_sa(Q, num_reads=200, seed=42):
    sampler = neal.SimulatedAnnealingSampler()
    bqm = dimod.BinaryQuadraticModel.from_qubo(
        {(i, j): Q[i, j] for i in range(N) for j in range(i, N) if Q[i, j] != 0},
        offset=0.0
    )
    sampleset = sampler.sample(bqm, num_reads=num_reads, seed=seed)
    return sampleset

for P in [1, 2, 5, 10, 20]:
    Q = build_Q(P)
    ss = solve_with_sa(Q)
    feas = 0
    best_feasible_E = np.inf
    for record in ss.record:
        z = np.array(record.sample, dtype=int)
        xs = z[:n]
        if w @ xs <= C:
            feas += 1
            best_feasible_E = min(best_feasible_E, record.energy)
    ratio = feas / len(ss)
    print(f"P={P:2d}  可行解比例={ratio:.2%}  最优可行能量={best_feasible_E:.2f}")

预期现象：P=1 时可行解比例明显偏低，采样器被超重的高价值组合吸引；P 增大后可行解比例上升并趋于 1，但 P 过大时能量尺度差异被放大，退火器在有限比特精度下可能更难分辨细微的价值差异。因此工程上常取“刚好使约束满足”的 P，而不是一味增大。

结果解码与可行性校验

采样返回的是 9 位比特串，前 5 位是物品选择，后 4 位是松弛变量。解码时只取前 5 位计算重量和价值，并独立校验 Σwᵢxᵢ ≤ C，不要依赖松弛变量是否自洽——退火器可能返回松弛比特与 x 不匹配的串，这类串即使能量低也应判为不可行或需修复。

def decode_and_check(z):
    xs = z[:n]
    weight = int(w @ xs)
    value = int(v @ xs)
    ok = weight <= C
    return xs.tolist(), weight, value, ok

# 取 P=10 的最优样本
Q = build_Q(10)
ss = solve_with_sa(Q, num_reads=500)
best = ss.first.sample
z = np.array([best[i] for i in range(N)], dtype=int)
print(decode_and_check(z))

若最优样本不可行，说明 P 仍偏小或采样不充分；可提高 P、增加 num_reads，或对不可行解做贪心修复（按价值密度降序尝试移除超重物品）后再比较目标值。

惩罚系数选择速查

本例的完整闭环是：建模（含松弛变量）→ 构造 Q → 枚举确认全局最优 → 退火采样 → 统计可行解比例 → 选 P → 解码校验。掌握这条链路后，换任意小规模 0-1 背包实例只需替换 w、v、C 并重算松弛比特数即可复用。