本文介绍如何在 Python 中利用 scipy.optimize.milp 求解一类关键数值优化问题：给定一组浮点数（如原子电荷），寻找一组正整数系数，使加权和在指定精度（如 ±1e−5）内逼近零，同时允许对原始数值进行微小扰动（±0.0005）——但本方案更优地将“误差吸收”直接建模为优化目标，避免手动调参。
本文介绍如何在 python 中利用 `scipy.optimize.milp` 求解一类关键数值优化问题：给定一组浮点数（如原子电荷），寻找一组正整数系数，使加权和在指定精度（如 ±1e−5）内逼近零，同时允许对原始数值进行微小扰动（如 ±0.0005）——但本方案更优地将“误差吸收”直接建模为优化目标，避免手动调参。
在量子化学、材料建模或金融计量等场景中，常需构造整数加权组合，使浮点数据的线性组合严格满足舍入后为零（例如保证分子总电荷四舍五入到 5 位小数后为 0）。用户提供的字典 d 包含各组分的原始值，而标量 n₁, n₂, … 是待求的固定结构但未知取值的正整数系数（即变量维度与 d 长度一致）。核心挑战在于：不能修改系数 nᵢ（它们由物理/业务逻辑决定），也不能任意调整 d[key]（仅允许极小范围扰动），因此需将“残差补偿”转化为可优化的目标。
本方案采用混合整数线性规划（MILP），将问题形式化为：
- 决策变量：n = [n₁, n₂, ..., nₙ] ∈ ℤ⁺ⁿ（正整数） + ε ∈ ℝ（连续误差项）
- 目标函数：最小化 |ε|（即 min ε，通过双不等式约束实现）
-
约束条件：
∑ᵢ nᵢ·dᵢ − small ≤ ε
small − ∑ᵢ nᵢ·dᵢ ≤ ε
等价于 |∑ᵢ nᵢ·dᵢ − small| ≤ ε，确保加权和与目标 small 的偏差不超过 ε
以下为完整可运行实现（要求 scipy ≥ 1.9.0）：
import numpy as np
from scipy.optimize import milp, Bounds, LinearConstraint
def solve(
d: dict[str, float],
small: float = 1e-5,
nmax: int = 20,
) -> tuple[np.ndarray, float]:
"""
求解整数系数 n，使得 |sum(n_i * d_i) - small| 最小化
Parameters:
-----------
d : dict[str, float]
输入字典，键为标识符，值为浮点系数
small : float
目标偏移量（通常为 1e-5 量级）
nmax : int
系数上限（防止搜索空间爆炸，默认 20）
Returns:
--------
n : np.ndarray of int
最优正整数系数向量
error : float
对应的最小绝对误差
"""
keys = list(d.keys())
d_array = np.array([d[k] for k in keys])
N = len(d_array)
# 目标：最小化误差 ε（n 不贡献成本）
c = np.concatenate((np.zeros(N), [1.0])) # [n₁..nₙ, ε]
# 变量类型：前 N 个为整数，最后 1 个为连续变量
integrality = np.ones(N + 1, dtype=np.uint8)
integrality[-1] = 0
# 边界：n_i ∈ [1, nmax]；ε ∈ (-∞, +∞)
bounds = Bounds(
lb=np.concatenate((np.ones(N), [-np.inf])),
ub=np.concatenate((np.full(N, nmax), [np.inf]))
)
# 约束：-ε ≤ sum(n_i * d_i) - small ≤ ε
# 即：sum(n_i * d_i) - small ≤ ε → [-d_array, 1] @ [n; ε] ≥ -small
# small - sum(n_i * d_i) ≤ ε → [d_array, 1] @ [n; ε] ≥ small
A = np.block([
[-d_array, [1.0]], # -∑n_i d_i + ε ≥ -small
[ d_array, [1.0]], # ∑n_i d_i + ε ≥ small
])
constraints = LinearConstraint(A, lb=[-small, small], ub=[np.inf, np.inf])
result = milp(c=c, integrality=integrality, bounds=bounds, constraints=constraints)
if not result.success:
raise RuntimeError(f"Optimization failed: {result.message}")
n_opt, error = np.split(result.x, [-1])
return n_opt.astype(int), float(error)
# 示例调用
if __name__ == "__main__":
d = {
'C1': 0.68759, 'C2': -0.21432,
'H1': 0.49062, 'H2': -0.13267,
'H3': 0.08092, 'O1': -0.86040,
}
n, err = solve(d, small=1e-5, nmax=20)
# 格式化输出（保留 5 位小数，整数系数无小数）
terms = [f"{int(ni)}*{vi:.5f}" for ni, vi in zip(n, d.values())]
print(" + ".join(terms) + f" = {1e-5:.0e} + {err:.1e}")
# 输出示例：1*0.68759 + 3*-0.21432 + 12*0.49062 + ... = 1e-05 + 0.0e+00
✅ 关键优势：
- 完全通用：自动适配任意长度字典 d 和对应数量的整数系数；
- 精度可控：通过 small 参数嵌入目标偏移，无需手动修改 d[key]；
- 强约束保障：nmax 防止组合爆炸，Bounds 保证正整数解；
- 数学严谨：将“舍入后为零”转化为 |∑nᵢdᵢ − small| ≤ ε，ε 最小化即等效于提升舍入鲁棒性。
⚠️ 注意事项：
- 若 nmax 过小导致无可行解，milp 将报错，此时应适度增大 nmax 或检查 small 是否合理；
- 实际应用中，若原始 d 值本身存在测量/计算误差，建议先对 d 做预处理（如统一截断至 6 位小数）以提升数值稳定性；
- 该方法不修改原始 d[key]，所有“调整”均由整数系数 nᵢ 的优化完成，符合问题中“标量 n 不应改变”的本质要求——此处 n 是待求变量，而非输入参数。
总结而言，此方案将看似经验性的数值凑整问题，升维为一个结构清晰、可验证、可扩展的整数优化任务，兼顾科学严谨性与工程实用性。
