物理约束MCMC与化学信息高斯过程协同用于反应网络发现

# 协同物理约束MCMC与化学信息高斯过程实现反应网络发现
来源: https://arxiv.org/html/2606.23757
###### 摘要

从稀疏、含噪声的化学时间序列数据中提取可解释的控制方程仍然困难，因为离散反应拓扑与连续动力学参数紧密耦合。我们提出PC-MCMC-CIGP，一个可复现的灰箱工作流，它结合了尖峰-板层拓扑采样、严格的守恒和热力学筛选，以及一个用于参数校准和实验设计的化学信息高斯过程(CIGP)残差模型。此方法的贡献并非孤立地提出新的MCMC或GP家族，而是将这些组件整合到一个具有明确不确定性感知采集选择的、受物理约束的工作流中。在H2+Br2基准测试中，约束采样器在我们的实验中能够区分基本自由基路径与具有欺骗性的现象学拟合。在苯乙烯环氧化反应中，CIGP优化循环相比于已报告的GP-BO基线，最终产率提升了12.5%。一项新的10种子采集研究表明，EI、GWU、PC-EI、不确定性采样、差异搜索和随机搜索具有不同的权衡关系：PC-EI大幅减少了低产率的BO建议，而EI类准则在最终产率性能上表现最强。

反应网络发现, 贝叶斯推断, 高斯过程, 实验设计

## 1 引言

从含噪声的时间序列数据中发现控制方程仍然是科学机器学习中的一个基本挑战(Brunton等，2016 (https://arxiv.org/html/2606.23757#bib.bib1); Schmidt和Lipson，2009 (https://arxiv.org/html/2606.23757#bib.bib3); Champion等，2019 (https://arxiv.org/html/2606.23757#bib.bib4))。尽管深度学习范式能够有效近似连续动态，但可解释性常常为黑箱插值而牺牲(Cranmer，2023 (https://arxiv.org/html/2606.23757#bib.bib5))。在化学动力学中，主要目标是*结构辨识*，定义为重建控制系统演化的离散反应拓扑，以恢复一个物理上可解释的因果图(Burnham等，2008 (https://arxiv.org/html/2606.23757#bib.bib11); Jiang等，2022 (https://arxiv.org/html/2606.23757#bib.bib15))。

机制发现被数学形式化为一个不适定逆问题，其特点是离散结构与连续动力学常数之间的强耦合(Walter和Pronzato，1997 (https://arxiv.org/html/2606.23757#bib.bib19); Davidescu和Jørgensen，2008 (https://arxiv.org/html/2606.23757#bib.bib16); Audoly等，2002 (https://arxiv.org/html/2606.23757#bib.bib18))。传统的贝叶斯推断常常受到这种相互依赖性的阻碍。标准参数估计技术缺乏明确的物理编码(Williams和Rasmussen，2006 (https://arxiv.org/html/2606.23757#bib.bib31))，而结构采样则常常受到巨大的组合空间和不足的物理约束的阻碍，导致数学上有效但物理上禁止的解(Enciso等，2021 (https://arxiv.org/html/2606.23757#bib.bib51))。

为了解决这些局限性，我们提出了一个统一的、受物理约束的贝叶斯框架，旨在协同物理流形上的离散采样与化学信息连续建模(Brubaker等，2012 (https://arxiv.org/html/2606.23757#bib.bib44))。该问题被视为结构推断和参数估计的协作求解。我们引入了一个多物理约束的尖峰-板层马尔可夫链蒙特卡洛(MCMC)采样器用于结构辨识。质量守恒、电子守恒和热力学势环等约束被编码为边界条件，以修剪后验分布并消除虚假机制。

对于参数估计，我们利用*化学信息高斯过程*(CIGP)构建一个混合灰箱模型。基于机制常微分方程被嵌入为先验均值函数，以捕获主导动力学趋势，而结构差异则通过非参数核进行建模(Raissi等，2017 (https://arxiv.org/html/2606.23757#bib.bib33); Chang和Zeng，2023 (https://arxiv.org/html/2606.23757#bib.bib38); Ma等，2020 (https://arxiv.org/html/2606.23757#bib.bib34))。物理参数和核超参数被联合优化以确保稳健估计(Cross等，2024 (https://arxiv.org/html/2606.23757#bib.bib37); Dalton等，2024 (https://arxiv.org/html/2606.23757#bib.bib53))。此外，我们还开发了物理感知的采集函数以促进自动实验设计。图1总结了两个阶段的工作流：PC-MCMC首先识别一个物理上允许的反应拓扑，然后CIGP在校准和主动学习期间将识别出的ODE模型用作GP均值函数。

本工作的主要贡献总结如下：

1. 1. 构建了一个固定维度的尖峰-板层MCMC工作流用于候选反应选择，并将质量、电荷和细致平衡约束作为硬性允许性检查。
2. 2. 使用化学信息高斯过程校准阶段作为灰箱残差模型，保留机理ODE作为先验均值，同时对系统差异进行建模。
3. 3. 提供了一个可复现的主动学习基准，包括10种子比较PC-EI、EI、GWU、差异搜索、不确定性采样和随机搜索。

图1：所提出的PC-MCMC-CIGP框架概述。上层面板展示了受物理约束的结构发现阶段，其中候选反应网络使用受质量守恒和热力学约束的尖峰-板层MCMC方案进行采样。下层面板展示了化学信息高斯过程(CIGP)，它将发现到的机理ODE模型嵌入为GP均值函数，同时非参数化地对结构差异进行建模。统一的框架使得同时进行机制识别、参数推断和物理感知的主动学习成为可能。
## 2 相关工作

### 2.1 现象学发现与机制发现

数据驱动的动力学定律发现通常通过符号回归或稀疏识别技术来实现(Brunton等，2016 (https://arxiv.org/html/2606.23757#bib.bib1); Rudy等，2017 (https://arxiv.org/html/2606.23757#bib.bib2); Schmidt和Lipson，2009 (https://arxiv.org/html/2606.23757#bib.bib3); Burlacu等，2020 (https://arxiv.org/html/2606.23757#bib.bib6))。诸如遗传编程和稀疏回归等算法可以从观测数据中恢复准确的宏观表达式，包括自由基链反应中出现的分数阶动力学(Cranmer，2023 (https://arxiv.org/html/2606.23757#bib.bib5); Otte，2014 (https://arxiv.org/html/2606.23757#bib.bib52))。它们在本工作中的局限性在于，一个准确的表达式并不一定等于一个反应机理：一个半阶速率定律可能拟合了HBr的形成，但却没有识别出基本的自由基拓扑和隐藏的中间体。

### 2.2 反应网络与神经ODE模型

反应网络重建也被形式化为使用MINLP、通量平衡类约束或超图搜索的确定性超结构选择(Langary和Nikoloski，2019 (https://arxiv.org/html/2606.23757#bib.bib13); Searson等，2007 (https://arxiv.org/html/2606.23757#bib.bib14); Willis和von Stosch，2016 (https://arxiv.org/html/2606.23757#bib.bib12); Pal等，2025 (https://arxiv.org/html/2606.23757#bib.bib29); Bonvin和Rippin，1990 (https://arxiv.org/html/2606.23757#bib.bib20))。这些方法可以显式地编码化学计量，但通常对替代拓扑的后验不确定性提供有限。最近的神经ODE和化学反应神经网络模型改进了动力学系统的柔性动力学学习(Kim等，2021 (https://arxiv.org/html/2606.23757#bib.bib54); Owoyele和Pal，2021 (https://arxiv.org/html/2606.23757#bib.bib55); Chang和Zeng，2023 (https://arxiv.org/html/2606.23757#bib.bib38))。它们对我们的目标是互补的：PC-MCMC-CIGP不是学习一个无约束的神经向量场，而是假设一个候选基本反应集，并估计哪些反应是激活的，因此输出仍然是一个可解释的反应步骤子集。

### 2.3 概率建模与实验设计

高斯过程(GP)被广泛用作校准和贝叶斯优化的概率替代模型(Williams和Rasmussen，2006 (https://arxiv.org/html/2606.23757#bib.bib31); Kocijan，2016 (https://arxiv.org/html/2606.23757#bib.bib35); Snoek等，2012 (https://arxiv.org/html/2606.23757#bib.bib49); Shields等，2021 (https://arxiv.org/html/2606.23757#bib.bib30))。标准GP-BO对于产率最大化是有效的，但它通常将实验视为一个黑箱响应面。物理信息与混合GP模型将机制结构纳入均值或核中(Kennedy和O'Hagan，2001 (https://arxiv.org/html/2606.23757#bib.bib32); Raissi等，2017 (https://arxiv.org/html/2606.23757#bib.bib33); Ma等，2020 (https://arxiv.org/html/2606.23757#bib.bib34); Cross等，2024 (https://arxiv.org/html/2606.23757#bib.bib37))。我们的CIGP组件遵循这一灰箱传统。因此，这里声称的新颖性是经过刻意限定的：我们提供了一个统一的、受物理约束的反应发现与主动学习工作流，而不是一个全新的GP推断算法。

## 3 PC-MCMC-CIGP框架

### 3.1 问题形式化

我们考虑一个具有Ns种物种和预先枚举的Nr个基本反应候选集R={Rj}j=1Nr的化学体系。此候选集被认为是已知的；未知的是二元活动向量γ∈{0,1}Nr、活动反应图G(γ)以及动力学参数k。令c(t)∈R≥0Ns表示浓度，令S(G)∈ZNs×Nr表示化学计量矩阵，其第j列是反应Rj的净化学计量向量，并在模拟过程中由γj掩蔽。质量作用动力学为

dc/dt = S(G) r(c(t); k, γ), (1)
其中 rj(c; k, γ) = γj kj ∏i ci^{νij^{-}}，νij^{-}为反应物化学计量阶数。反应速率遵循阿伦尼乌斯参数化 k = A exp(-Ea/RT) 或以下描述的细致平衡参数化(Horn和Jackson，1972 (https://arxiv.org/html/2606.23757#bib.bib21); Vlad和Ross，2009 (https://arxiv.org/html/2606.23757#bib.bib25))。

机制发现数据是稀疏的含噪声轨迹 D = {(tm, ym)}m=1M，其中ym包含c(tm)在高斯观测模型下的观测分量。主动学习阶段使用不同的符号：u∈U表示可控实验条件，如初始浓度、温度和停留时间，y(u)表示标量目标，如产物产率。这分离了状态变量c(t)和设计变量u。

### 3.2 受物理约束的拓扑搜索

为了缓解反应网络发现中固有的组合爆炸问题，我们在一个物理允许流形上形式化了一个贝叶斯变量选择问题(Mitchell和Beauchamp，1988 (https://arxiv.org/html/2606.23757#bib.bib45); George和McCulloch，1993 (https://arxiv.org/html/2606.23757#bib.bib46))。第j个候选反应的存在由二元潜变量γj∈{0,1}控制。

有效速率形式化。为了确保与主导动力学一致的结构稀疏性，我们将有效速率常数定义为一个掩蔽变量：
kj^{eff} = γj · exp(θj), (2)
其中θj是辅助对数动力学参数。这种乘法形式确保了当反应被停用时(γj=0)，无论θj的值如何，它对微分方程(式1)的贡献恰好为零(kj^{eff}≡0)。

先验指定。为了引入严格稀疏性，在辅助参数θj上施加一个尖峰-板层先验以正则化搜索空间：
p(θj | γj) = (1-γj) N(θj; 0, ε²) + γj N(θj; μθ, σθ²), (3)
其中尖峰分量(由狭窄高斯近似，ε→0)约束非活动路径的辅助参数以实现可辨识性，而板层分量捕获活动反应的参数不确定性(Ishwaran和Rao，2005 (https://arxiv.org/html/2606.23757#bib.bib47); Piironen和Vehtari，2017 (https://arxiv.org/html/2606.23757#bib.bib10))。对γ和θ的后验推断使用Metropolis–Hastings马尔可夫链蒙特卡洛(MCMC)采样器进行(Metropolis等，1953 (https://arxiv.org/html/2606.23757#bib.bib41); Peskun，1973 (https://arxiv.org/html/2606.23757#bib.bib42); Girolami和Calderhead，2011 (https://arxiv.org/html/2606.23757#bib.bib43))。

#### 物理约束强制执行。

为了消除数学上允许但物理上禁止的解，通过拒绝采样施加三类硬约束。

质量和电子守恒。令A∈R^{Na×Ns}为原子组成矩阵。一个有效的反应j必须满足
A·S·j = 0, (4)
确保原子物种和电荷的严格守恒。

热力学一致性。对于每对检测到的正反反应，细致平衡通过潜无量纲化学势μ约束速率比：
ln(k_fwd / k_rev) = - ∑_{i=1}^{Ns} ν_{ij} μ_i, (5)
其中ν_{ij}是正向反应中物种i的净化学计量系数。势{μ_i}与动力学参数在宽盒先验下联合推断，这避免了对瞬态自由基需要表格化热力学数据。违反质量、电荷、速率界限或细致平衡约束的候选移动在似然评估之前被拒绝。这种基于拒绝的实现很简单，但随着可逆对数量的增加可能会降低接受率；我们将其报告为一个局限性，而不是假设可扩展性是自动的。

### 3.3 化学信息高斯过程

为了将参数估计与数值分离