Python 数据与优化耦合：结果诊断与可信度验证

Python 侧耦合的典型失败不是“算错了”，而是“每次算得略微不同”：求解器容差、线程数、哈希随机化都会让同一个 FMU 在两次调用中给出末位不同的输出。一旦目标函数带噪，有限差分梯度就被噪声吞掉，优化器会朝错误方向走。诊断集中在三件事：调用契约是否干净、梯度是否被噪声主导、参数组合是否可辨识。

调用契约：谁持有状态

pyfmi.load_fmu 返回的模型对象是有状态的：reset() 清空、initialize() 重新初始化、simulate() 推进。优化循环里忘记 reset()，第二次评估会从上次末态继续，目标函数值随迭代历史变化。这类错误表现为“同一参数两次评估结果不同”，最容易定位也最容易复现。

import numpy as np
from scipy.optimize import minimize
from pyfmi import load_fmu

y_ref = np.loadtxt("exp_T_out.txt")
fmu = load_fmu("plant.fmu")                   # Co-Simulation FMU

def J(theta):
    fmu.reset()                               # 必须复位，否则状态跨迭代泄漏
    fmu.set("k", float(theta[0]))
    fmu.set("c", float(theta[1]))
    fmu.initialize()
    res = fmu.simulate(start_time=0.0, final_time=20.0,
                       options={"ncp": 2000, "result_handling": "memory"})
    r = res["T_out"] - y_ref
    return float(np.sqrt(np.mean(r**2)))

fmpy.simulate_fmu 走无状态路径，每次调用重建实例，天然避免状态泄漏，但启动开销更大。短时程小模型用它更省心，长时程批量扫描则复用 pyfmi 实例更划算。

噪声目标函数的梯度步长

设每次评估的目标值带标准差 $\sigma$ 的随机扰动，中心差分梯度的误差为

$$\varepsilon_{\nabla}(h)\approx \frac{\sigma}{\sqrt{2}\,h}+\frac{h^{2}}{6}\max\lvert\partial^{3}J\rvert,\qquad h^{*}=\left(\frac{3\sigma}{\lvert\partial^{3}J\rvert}\right)^{1/3}$$

实测目标 $J=\mathrm{RMS}(T_{out}-T_{ref})=1.234\,\mathrm{K}$，重复评估标准差 $\sigma=2.0\times10^{-6}\,\mathrm{K}$（来自 FMU 求解器容差 $10^{-8}$ 与输出插值），$\max\lvert\partial^3 J\rvert=1.0\times10^{-3}\,\mathrm{K}$。按上式 $h^{*}=(3\times2.0\times10^{-6}/10^{-3})^{1/3}=0.182$，代入得 $\varepsilon_{\nabla}=2.0\times10^{-6}/(1.414\times0.182)+0.182^2/6\times10^{-3}=7.8\times10^{-6}+5.5\times10^{-6}=1.3\times10^{-5}$。

对照组用 scipy 默认的 eps=1e-8：噪声项为 $2.0\times10^{-6}/(1.414\times10^{-8})=141$，而真实梯度量级只有 $0.8\,\mathrm{K}$ 每单位参数，梯度被噪声完全淹没。这就是“L-BFGS-B 在前 5 次迭代就报收敛、参数却离真值很远”的直接原因。

参数可辨识性的定量检查

即使梯度可靠，若灵敏度矩阵接近奇异，参数组合也无法分离。用下式判断：

$$\mathrm{Cov}(\hat{\theta})\approx \sigma^{2}\left(\mathbf{S}^{\mathsf{T}}\mathbf{S}\right)^{-1},\qquad \kappa=\frac{\sigma_{\max}(\mathbf{S})}{\sigma_{\min}(\mathbf{S})}$$

某三参数问题中 $\mathbf{S}$ 的奇异值为 $4.2$、$0.91$、$0.014$，$\kappa=300$。第三个参数的方差被放大 $300^2=9\times10^4$ 倍，其置信区间跨过整个物理可行域。此时增加迭代次数无用，只能改实验设计：增加对该参数敏感的激励段，或把它固定为标称值后再辨识其余两个。

症状、根因与判定试验

现象	根因	判定试验
同一参数两次评估目标值差 $3\times10^{-6}$	未调用 reset()，状态跨迭代泄漏	连续调用 J 两次，差值应小于 $10^{-12}$
优化器 5 次迭代即报收敛	默认差分步长被噪声淹没	把 eps 从 $10^{-8}$ 改到 $0.18$ 后重跑
单机可复现、集群不可复现	线程数或哈希随机化不同	固定 OMP_NUM_THREADS=1、PYTHONHASHSEED=0 后重跑
某参数置信区间跨全域	灵敏度矩阵条件数过大	计算 $\kappa$，大于 100 时改激励设计
优化中途被杀且无中间结果	无超时与检查点机制	记录每次评估的输入哈希与耗时，复现被杀时刻的参数
目标函数随迭代单调下降但物理量违反守恒	惩罚项权重过大，把物理残差压成了常数	单独输出接口功残差曲线，与 $J$ 对照

进程隔离、超时与重试

单次评估耗时 $1.8\,\mathrm{s}$，400 次评估的预算为 $720\,\mathrm{s}$；若每次用 subprocess 隔离，进程启动另加 $0.35\,\mathrm{s}$，合计 $140\,\mathrm{s}$，占总预算 19%。隔离的收益是失败可恢复：FMU 因数值问题崩溃时，主进程仍能拿到退出码并重试。约定为单次评估超时 $10\,\mathrm{s}$、最多重试 2 次、失败比例超过 5% 即终止优化并保留检查点。

export PYTHONHASHSEED=0 OMP_NUM_THREADS=1 MKL_NUM_THREADS=1
python -c "import numpy,scipy,pyfmi,fmpy;print(numpy.__version__,scipy.__version__)"
sha256sum plant.fmu exp_T_out.txt

参考文献

1. Nocedal J., Wright S. J., Numerical Optimization, 2nd ed., Springer, 2006, §8.1.
2. Moré J. J., Wild S. M., “Estimating Derivatives of Noisy Simulations”, ACM Transactions on Mathematical Software, 38(3), 2012, Article 20.
3. Modelon AB, PyFMI User Guide, 2023.
4. Zuluaga C., fmpy — Simulate Functional Mock-up Units in Python, 2019.
5. Virtanen P., Gommers R., Oliphant T. E., et al., “SciPy 1.0: fundamental algorithms for scientific computing in Python”, Nature Methods, 17, 2020, pp. 261-272.
6. Andersson J. A. E., Gillis J., Horn G., Rawlings J. B., Diehl M., “CasADi: a software framework for nonlinear optimization and optimal control”, Mathematical Programming Computation, 11, 2019, pp. 1-36.
7. Beck J. V., Arnold K. J., Parameter Estimation in Engineering and Science, Wiley, 1977.