本文介绍一种远优于暴力枚举的专业方法——基于 SciPy milp 的二元整数线性规划（BILP）建模，可在秒级内判定给定 16 个数能否重排为行、列及两条主对角线之和均为 0 的 4×4 矩阵，并输出可行解（若存在）。

本文介绍一种远优于暴力枚举的专业方法——基于 scipy `milp` 的二元整数线性规划（bilp）建模，可在秒级内判定给定 16 个数能否重排为行、列及两条主对角线之和均为 0 的 4×4 矩阵，并输出可行解（若存在）。

在组合约束求解问题中，面对 16! ≈ 2.1×10¹³ 种全排列的暴力搜索空间，传统穷举法完全不可行。原始代码尝试用 itertools.permutations 截断采样，本质上仍属随机试探，既无法保证解的存在性证明，也无法在合理时间内收敛。真正可扩展、可验证的解法需转向数学规划建模：将“重排”转化为“分配决策”，用二元变量刻画每个数字是否被放置于特定位置，并将所有约束（唯一性、行列和、对角线和）统一表达为线性等式。

核心思想是引入三维二元决策变量 $ x_{k,i,j} \in {0,1} $，其中：

• $ k = 0,\dots,15 $ 表示第 $k$ 个给定数字（共 16 个）；
• $ i,j = 0,\dots,3 $ 表示目标网格的行与列索引。

$ x_{k,i,j} = 1 $ 当且仅当第 $k$ 个数被填入位置 $(i,j)$。由此可构建以下五类线性约束：

✅ 每个数字恰好使用一次：
$$ \sum{i=0}^{3}\sum{j=0}^{3} x_{k,i,j} = 1 \quad \forall\,k $$

✅ 每个格子恰好填一个数字：
$$ \sum{k=0}^{15} x{k,i,j} = 1 \quad \forall\,i,j $$

✅ 每行和为 0：
$$ \sum{k=0}^{15} \sum{j=0}^{3} (\text{nums}[k] \cdot x_{k,i,j}) = 0 \quad \forall\,i $$

✅ 每列和为 0：
$$ \sum{k=0}^{15} \sum{i=0}^{3} (\text{nums}[k] \cdot x_{k,i,j}) = 0 \quad \forall\,j $$

✅ 两条主对角线和为 0：
$$ \sum{k=0}^{15} \sum{i=0}^{3} (\text{nums}[k] \cdot x{k,i,i}) = 0, \quad \sum{k=0}^{15} \sum{i=0}^{3} (\text{nums}[k] \cdot x{k,i,3-i}) = 0 $$

这些约束共同构成一个混合整数线性规划（MILP）问题。由于目标函数仅用于可行性判定（无优化目标），我们设 $c = \mathbf{0}$，调用 scipy.optimize.milp 求解即可。

以下是完整可运行实现（已适配原题数据）：

import numpy as np
from scipy.optimize import milp

# 原始输入（注意：原始数据经验证无解；此处按答案提示修正第二行首元素为 2）
numbers = np.array([5, -5, 3, -2,  # row 0
                    2, -1, 9, -10,  # row 1 ← modified from 1 to 2
                    -7, -4, 7, 4,   # row 2
                    -3, -8, 2, 8]) # row 3

n = 4
N = n * n  # 16

# 构造三维广播模板：(16, 4, 4)，每个 slice 是该数字在全部位置的副本
nums_3d = numbers.reshape(N, 1, 1)
nums_3d = np.tile(nums_3d, (1, n, n))

# 初始化约束矩阵 A 和右端向量 b
A_list, b_list = [], []

# 约束 1：每个数字用且仅用一次 → 16 条等式
for k in range(N):
    A_k = np.zeros((N, n, n))
    A_k[k] = 1
    A_list.append(A_k.flatten())
    b_list.append(1)

# 约束 2：每个格子填且仅填一个数字 → 16 条等式
for i in range(n):
    for j in range(n):
        A_ij = np.zeros((N, n, n))
        A_ij[:, i, j] = 1
        A_list.append(A_ij.flatten())
        b_list.append(1)

# 约束 3：每行和为 0 → 4 条等式
for i in range(n):
    A_row = np.zeros((N, n, n))
    A_row[:, i, :] = nums_3d[:, i, :]
    A_list.append(A_row.flatten())
    b_list.append(0)

# 约束 4：每列和为 0 → 4 条等式
for j in range(n):
    A_col = np.zeros((N, n, n))
    A_col[:, :, j] = nums_3d[:, :, j]
    A_list.append(A_col.flatten())
    b_list.append(0)

# 约束 5：主对角线 & 反对角线和为 0 → 2 条等式
A_diag = np.zeros((N, n, n))
A_diag[:, range(n), range(n)] = nums_3d[:, range(n), range(n)]
A_list.append(A_diag.flatten())
b_list.append(0)

A_anti = np.zeros((N, n, n))
A_anti[:, range(n), range(n-1, -1, -1)] = nums_3d[:, range(n), range(n-1, -1, -1)]
A_list.append(A_anti.flatten())
b_list.append(0)

# 拼接为标准 MILP 形式：A @ x == b
A_eq = np.vstack(A_list)
b_eq = np.array(b_list)

# 求解：变量 x ∈ [0,1] 且为整数
res = milp(
    c=np.zeros(A_eq.shape[1]),           # 无目标，仅求可行解
    constraints=(A_eq, b_eq, b_eq),     # (A, lb, ub) → 等式约束
    bounds=(0, 1),                        # 二元变量
    integrality=1                         # 强制整数解
)

if res.success:
    x_sol = np.round(res.x).astype(int)  # 解码二元解
    # 重构 4x4 网格：对每个 (i,j)，找出满足 x[k,i,j]==1 的 k，取 numbers[k]
    grid_sol = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            k = np.argmax(x_sol.reshape(N, n, n)[:, i, j])
            grid_sol[i, j] = numbers[k]

    print("✅ 找到可行解：")
    print(grid_sol)
    # 验证
    assert np.allclose(grid_sol.sum(axis=0), 0), "列和不为 0"
    assert np.allclose(grid_sol.sum(axis=1), 0), "行和不为 0"
    assert np.allclose(np.diag(grid_sol).sum(), 0), "主对角线和不为 0"
    assert np.allclose(np.diag(np.fliplr(grid_sol)).sum(), 0), "反对角线和不为 0"
    print("✔ 所有约束验证通过。")
else:
    print("❌ 无可行解（原始数据不可行，或修正后仍无解）")

? 关键注意事项：

• 原始数据不可行：对题目所给 [5,-5,3,-2,1,-1,9,-10,-7,-4,7,4,-3,-8,2,8]，该模型会明确返回 success=False，即数学上不存在满足全部约束的重排——这是暴力法永远无法给出的确定性结论。
• 数值稳定性：milp 对浮点精度敏感，建议使用 np.round(res.x) 解码，并用 np.allclose(..., atol=1e-8) 验证。
• 扩展性：本框架天然支持任意 $n\times n$ 网格（只需调整 n 和约束循环），且时间复杂度不随 $n!$ 增长，而是取决于 MILP 求解器性能（对 $n=4$ 通常
• 替代方案：若需更高性能或更大规模，可迁移至 ortools（Google）、PuLP + CBC 或商业求解器（Gurobi/CPLEX）。

综上，面对高维组合约束问题，放弃暴力、拥抱建模，是工程实践中兼顾效率、正确性与可解释性的必由之路。