本文介绍一种远优于暴力枚举的专业方法——基于 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）。
综上，面对高维组合约束问题，放弃暴力、拥抱建模，是工程实践中兼顾效率、正确性与可解释性的必由之路。
