本文介绍如何在Python中利用scipy.optimize.milp求解一类关键数值优化问题：给定浮点系数字典和固定小目标值（如1e-5），寻找一组正整数权重，使加权和经5位舍入后严格等于零，且各系数扰动不超过±0.0005。
本文介绍如何在python中利用`scipy.optimize.milp`求解一类关键数值优化问题：给定浮点系数字典和固定小目标值（如1e-5），寻找一组正整数权重，使加权和经5位舍入后严格等于零，且各系数扰动不超过±0.0005。
在量子化学、计量校准或金融建模等场景中，常需满足“加权和舍入后精确为零”这一强约束条件。例如，原子电荷分配需满足总电荷守恒（∑nᵢ·qᵢ ≈ 0），但原始浮点电荷值经标准舍入（如保留5位小数）后，因累积误差可能导致∑round(nᵢ·qᵢ, 5) ≠ 0。此时，若权重 nᵢ（如原子个数、重复次数）必须为正整数且不可修改，而允许对原始系数 qᵢ 进行微调（如 ±0.0005 内浮动），则问题转化为一个混合整数线性规划（MILP）问题。
核心思想是：将权重 nᵢ 设为整数决策变量，引入连续变量 error 表征实际加权和与目标值 small 的偏差，并通过线性约束强制 |∑nᵢ·qᵢ − small| ≤ error。目标是最小化 error（即逼近精度），同时保证 nᵢ ∈ ℤ⁺ 且有合理上界（避免过大权重导致物理意义失效）。
以下为完整可运行实现：
import numpy as np
from scipy.optimize import milp, Bounds, LinearConstraint
def solve_rounding_balance(
d: dict[str, float],
small: float = 1e-5,
max_weight: int = 20,
max_perturb: float = 0.0005,
) -> tuple[np.ndarray, float]:
"""
求解加权舍入平衡问题：找到正整数权重 n_i，使得 sum(n_i * q_i) ≈ small，
且当所有 q_i 被舍入至5位小数时，加权和严格为0（等价于 |sum(n_i * q_i) - small| 极小）。
Parameters
----------
d : dict[str, float]
系数字典，键为标识符，值为原始浮点系数
small : float
目标偏移量（通常为1e-5量级）
max_weight : int
权重上限，防止解无意义地增大（默认20）
max_perturb : float
单个系数允许的最大调整量（用于后续舍入可行性分析）
Returns
-------
n : np.ndarray
最优正整数权重数组（按d.values()顺序）
error : float
实际残差 |sum(n_i * q_i) - small|
"""
values = np.array(list(d.values()))
N = len(values)
# 目标：最小化绝对误差（error变量）
c = np.concatenate((np.zeros(N), [1.0])) # [n₁,...,nₙ, error]
# 变量类型：前N个为整数，最后一个（error）为连续变量
integrality = np.ones(N + 1, dtype=np.uint8)
integrality[-1] = 0
# 边界：n_i ∈ [1, max_weight]；error ≥ 0（由约束隐含，此处设为无下界更稳健）
bounds = Bounds(
lb=np.concatenate((np.ones(N), [-np.inf])),
ub=np.concatenate((np.full(N, max_weight), [+np.inf]))
)
# 约束：|∑ n_i * q_i - small| ≤ error
# 等价于两个线性不等式：
# ∑ n_i * q_i - small ≤ error → -∑ n_i * q_i + error ≥ -small
# -∑ n_i * q_i + small ≤ error → ∑ n_i * q_i + error ≥ small
A = np.block([
[-values, [1.0]], # -q₁n₁ -... -qₙnₙ + 1·error ≥ -small
[ values, [1.0]], # q₁n₁ +... +qₙnₙ + 1·error ≥ 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_opt,) = np.split(result.x, [-1])
return n_opt.astype(int), float(error_opt)
# 示例调用
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_rounding_balance(d, small=1e-5, max_weight=20)
# 验证：计算原始加权和、舍入后加权和
weighted_sum = np.sum(n * np.array(list(d.values())))
rounded_sum = np.sum(n * np.round(list(d.values()), 5))
print("=== 优化结果 ===")
print(f"权重 n = {list(n)}")
print(f"原始加权和 ∑nᵢ·qᵢ = {weighted_sum:.8f}")
print(f"目标偏移 small = {1e-5}")
print(f"残差 |∑nᵢ·qᵢ − small| = {abs(weighted_sum - 1e-5):.2e}")
print(f"5位舍入后加权和 = {rounded_sum:.0f} （应为0）")
# 展示各分项（便于人工核对）
terms = [f"{ni}*{qi:.5f}" for ni, qi in zip(n, d.values())]
print("\n分项展开：", " + ".join(terms))
关键注意事项与实践建议：
- ✅ 舍入可行性保障：本解法确保 |∑nᵢ·qᵢ − small| 极小（常达机器精度），故对 qᵢ 做 ±0.0005 内微调即可满足 round(∑nᵢ·qᵢ, 5) == 0。例如，若最优解残差为 2.3e-6，只需将任一 qᵢ 增加 2.3e-6 / nᵢ 即可完全消除舍入误差。
- ⚠️ 权重上界设定：max_weight 是算法收敛性与物理合理性的平衡点。过小（如 5）可能导致无可行解；过大（如 100）易得数学最优但化学/工程上不可接受的解（如某原子重复100次）。建议从 10–20 起步，结合领域知识调整。
- ? 多解与鲁棒性：milp 返回首个最优解。若需全部最优解或带权重的帕累托前沿（如同时最小化 error 和 ∑nᵢ），需改用分支定界框架（如 pyomo + glpk）或添加正则项。
- ? 扩展性说明：该函数天然支持任意长度字典 d，只要 len(d) 与权重维度一致。对于超大规模问题（>100变量），建议启用 method='highs'（SciPy 1.9+ 默认）并预处理系数缩放以提升数值稳定性。
综上，此方案将领域约束精准建模为标准MILP问题，兼顾严谨性与实用性，是解决“舍入致零”类优化任务的通用范式。
