给定约束的列联表填充：混合整数二次规划求解技术问询

API 探路者

2026-8-3

列联表整数求解问题

问题背景

需要在给定边缘分布的前提下，求解列联表的非负整数单元格值$x_{ij}$，列联表如下：

2574	2572	3393	3768	3822	b（行和，近似值）	e（预期加权和）
$x_{11}$	$x_{12}$	$x_{13}$	$x_{14}$	$x_{15}$	187	23846753.74
$x_{21}$	$x_{22}$	$x_{23}$	$x_{24}$	$x_{25}$	3	324024.64
$x_{31}$	$x_{32}$	$x_{33}$	$x_{34}$	$x_{35}$	13755	1489591510.50
$x_{41}$	$x_{42}$	$x_{43}$	$x_{44}$	$x_{45}$	543	76173239.22
$x_{51}$	$x_{52}$	$x_{53}$	$x_{54}$	$x_{55}$	68	8188751.57
$x_{61}$	$x_{62}$	$x_{63}$	$x_{64}$	$x_{65}$	1332	172945247.86
$x_{71}$	$x_{72}$	$x_{73}$	$x_{74}$	$x_{75}$	361	41675606.70

其中列和为精确值，行和b为近似值。额外约束：给定常数因子$F = (13336.41847153, 102412.73466321, 41811.01724119, 78689.83110577, 282353.66682778)^T$，要求满足$X·F ≈ e$。

优化建模

将问题建模为以下优化问题：

目标函数：$\min (| X C – d |^2 + | X F – e |^2)$
约束条件：
1. $A X = b$
2. $x_{ij} ≥ 0$
3. $x_{ij} \in \mathbb{N}$

其中：

$b$：列和的精确值
$A$：全1向量，用于计算每一列的和
$d$：行和的近似值
$C$：全1向量，用于计算每一行的和
$e$：预期加权和向量
$F$：权重常数向量

尝试的求解代码及报错

使用cvxpy尝试求解，代码如下：

import cvxpy as cp
import numpy as np

F = np.array([13336.41847153, 102412.73466321, 41811.01724119, 78689.83110577, 282353.66682778])
e = [23846753.74, 324024.64, 1489591510.50, 76173239.22, 8188751.57, 172945247.86, 41675606.70]
d = [187., 3., 13755., 543., 68., 1332., 361.]
C = np.ones(len(F))
b = np.array([2574, 2572, 3393, 3768, 3822])
A = np.ones(len(d))

x = cp.Variable((len(e), len(b)), integer=True)
cost = cp.sum_squares(x @ C - d) + cp.sum_squares(x @ F - e)
objective = cp.Minimize(cost)
constraint_gt0 = x >= 0
constraint_eq = A @ x == b
    
problem = cp.Problem(objective, [constraint_gt0, constraint_eq])
solution = problem.solve()

运行后报错：

Either candidate conic solvers (['GLPK_MI', 'SCIPY']) do not support the cones output by the problem (SOC, NonNeg, Zero), or there are not enough constraints in the problem.

移除integer=True约束后代码可运行，但无法找到解。需要批量求解多个独立此类表格，因此需要通用代码方案，而非针对此示例的解。

问题

该问题是否规范且可解？
是否如报错提示存在欠约束情况？
为何会出现二阶锥（SOC）？原以为这是带额外约束的最小二乘问题。
如何使用cvxpy或其他Python包求解该问题？
若近似非整数解更易实现，该如何求解？

解答

1. 问题的规范性与可解性

该问题是规范的，但整数版本的可解性取决于约束相容性：

列和总和为$2574+2572+3393+3768+3822=16129$，行和总和为$187+3+13755+543+68+1332+361=16249$，两者存在差异——行和是近似值，因此严格满足列和约束与行和精确匹配不可能，只能通过最小二乘拟合行和近似值，同时满足列和精确约束。
只要存在非负整数矩阵$X$满足列和约束，问题就有可行解。从列和总和来看，这类矩阵存在（比如按比例分配），因此问题可解。

2. 是否欠约束

并非欠约束，核心是约束相容性问题：

列和约束是5个等式约束，变量共$7×5=35$个，看似变量多于约束，但非负整数约束会将可行域限制为离散有界集合（每个$x_{ij}$不超过对应列和）。
报错中的“欠约束”可能是指松弛（非整数）版本中仅列和约束无法唯一确定解，但目标函数的最小二乘项会收敛到唯一解；之前无法找到解是因为代码约束写法错误，导致约束不相容。

3. 为何出现二阶锥（SOC）

cvxpy会自动将平方和目标函数转换为二阶锥约束：

最小二乘的$|y|^{2$可表示为SOC约束（$|y|}2 ≤ t$等价于$(t, y)$属于二阶锥），使用cp.sum_squares时，cvxpy会将其建模为SOC锥优化问题。而默认的MIP求解器（GLPK_MI、SCIPY）不支持SOC锥，因此触发报错。

4. 整数版本求解方案

方案1：更换支持MIP+SOC的求解器

安装商业求解器（如Gurobi、CPLEX，有学术许可证），指定求解器即可：

# 安装Gurobi后修改求解调用
solution = problem.solve(solver=cp.GUROBI)

方案2：将目标函数线性化

将平方和目标替换为绝对值误差的线性和（$|XC-d|_1 + |XF-e|_1$），问题变为混合整数线性规划（MILP），开源求解器GLPK_MI支持：

# 修改目标函数为L1范数
res1 = x @ C - d
res2 = x @ F - e
cost = cp.sum(cp.abs(res1)) + cp.sum(cp.abs(res2))
objective = cp.Minimize(cost)

方案3：使用启发式算法

对于大规模问题，启发式算法（如差分进化）更高效，可使用scipy.optimize.differential_evolution：

from scipy.optimize import differential_evolution
import numpy as np

# 将变量展平为一维数组
def objective_func(x_flat):
    x = x_flat.reshape(7,5)
    res1 = x @ C - d
    res2 = x @ F - e
    return np.sum(res1**2) + np.sum(res2**2)

# 列和约束函数
def constraint_col_sum(x_flat):
    x = x_flat.reshape(7,5)
    col_sums = x.sum(axis=0)
    return col_sums - b

# 变量边界：每个x_ij ∈ [0, 对应列和]
bounds = []
for col_sum in b:
    bounds += [(0, col_sum)]*7

# 差分进化求解，处理整数约束
result = differential_evolution(
    objective_func,
    bounds=bounds,
    constraints={'type': 'eq', 'fun': constraint_col_sum},
    integers=True,
    seed=42
)

# 转换为矩阵形式的解
x_opt = result.x.reshape(7,5)

5. 近似非整数解的求解

之前非整数版本无法找到解是因为约束写法错误：A @ x == b维度不匹配，正确约束应为每一列的和等于b对应元素，修正后的代码如下：

x = cp.Variable((len(e), len(b)))
cost = cp.sum_squares(x @ C - d) + cp.sum_squares(x @ F - e)
objective = cp.Minimize(cost)
constraint_gt0 = x >= 0
# 正确的列和约束：每列求和等于b的对应值
constraint_eq = cp.sum(x, axis=0) == b
    
problem = cp.Problem(objective, [constraint_gt0, constraint_eq])
solution = problem.solve()

# 获取非整数解
x_opt = x.value

若需要将非整数解转为整数，可先四舍五入，再逐列调整单元格值以满足列和约束，或用非整数解作为MILP的初始点优化。

内容的提问来源于stack exchange，提问作者hfs