具有状态依赖事件的微分代数方程的可微参数优化

# 带有状态依赖事件的 DAE 的可微参数优化

来源: https://arxiv.org/html/2605.05395

Ion Matei  
Fujitsu Research of America  
[email protected]

& Maksym Zhenirovskyy  
Fujitsu Research of America  
[email protected]

& Anthony Wong  
Fujitsu Research of America  
[email protected]

###### 摘要

带有状态依赖事件的微分代数方程（DAEs）出现在其连续动态受代数方程约束并被模式切换、切换逻辑、冲击或状态重置所中断的系统中。对于此类系统进行基于梯度的参数学习具有挑战性，因为代数变量是隐式定义的，事件时间依赖于参数，且重置映射引入了不连续性。本文研究了带有事件的半显式 DAE 的可微参数优化。我们将学习问题表述为一个带有 DAE 动态、代数约束、守卫方程和重置映射的约束最小二乘问题。然后，我们开发了两种互补的梯度计算策略。第一种是通过模拟的自动微分方法，它在向量场内部求解代数变量，使用隐函数定理对代数求解过程进行微分，并通过分段可微积分处理事件。第二种是显式离散伴随方法，它将前向模拟表示为一个事件分割的残差系统，并通过求解光滑段和事件残差的拉格朗日乘子来计算梯度。该 formulation 明确了伴随方法中的残差项是等式约束，而非启发式惩罚项。我们从梯度解释、事件时间处理、实现复杂性和局部有效性方面比较了这两种方法。两种方法都为前向模拟所选的事件路径提供梯度，并在固定事件顺序和横截守卫穿越的条件下有效。

## 1 引言

微分代数方程（DAEs）是物理和工程系统常见的建模形式，其动态结合了微分演化与代数约束。它们出现在电路、受约束的机械系统、化学工艺、电网以及基于方程的建模环境中。关于数值 DAE 积分和灵敏度分析的经典文献为这一设置奠定了基础 (Ascher and Petzold,1998; Hairer and Wanner,1996)。在这些应用中，许多情况的连续动态会被离散事件所中断：开关打开或关闭、控制器改变模式、约束变为活跃、发生冲击，或者当阈值被跨越时状态变量被重置。由此产生的模型是混合 DAE：它们在事件之间连续演化，但其状态和代数变量必须在事件触发的不连续处保持一致。这种观点与 Modelica DAE 表示法一致，其中方程可以定义带有不连续性、可变结构和离散事件控制的 DAEs (Modelica Association,2025)。

混合 DAE 的参数估计和设计优化具有挑战性，因为模拟器不仅仅是一个从参数到输出的平滑映射。代数变量由约束隐式定义，必须在积分期间以及重置后重新计算。事件时间也依赖于参数；改变参数可以移位守卫条件达到的瞬间，从而改变连续动态积分的时间间隔。最后，事件通过重置映射引入不连续性。即使事件序列是固定的，梯度也必须考虑连续动态、事件时间敏感性、重置映射以及事件后的代数一致性。

最近的微分 ODE 求解器使得基于梯度的学习对于连续时间模型变得切实可行。如 torchdiffeq (Chen,2018) 和 diffrax (Kidger,2021) 等库支持通过 ODE 求解器的反向传播，并通过 odeint_event 包含可微的事件处理，标准示例包括弹跳球模型 (Chen et al.,2021, 2018)。相关的 Neural ODE 工作也考虑了事件序列建模，其中事件作为潜在连续时间模型的观察值或输入，而非硬性物理约束 (Kuleshov et al.,2024)。这些发展表明，带有事件的可微模拟是有用的且越来越易于获取。

然而，大多数现有工作是为 ODE 或潜在 ODE 架构 formulated 的。它们没有直接解决 DAE 的参数学习问题，其中代数约束必须在事件驱动的模拟中一致地求解和微分。

本文研究了带有状态依赖事件的半显式 DAE 的可微参数学习。我们考虑如下形式的系统：

$$
\begin{aligned}
\dot{x}(t) &= f(t, x(t), z(t), p), \\
0 &= g(t, x(t), z(t), p), \\
y(t) &= h(t, x(t), z(t), p),
\end{aligned}
$$

其中 $x$ 是微分状态，$z$ 是代数变量，$p$ 是要学习的参数。事件由守卫函数 $\phi_e(t, x, z, p) = 0$ 指定，以及将事件前变量映射到事件后变量的重置规则。学习问题被表述为在 DAE 动态、代数方程、事件守卫和重置映射下的输出匹配。

我们开发并比较了两种梯度计算策略。第一种是通过模拟的自动微分（AD）方法。它将模拟器视为从参数到预测输出的缩减映射，并直接对复合损失进行微分。对于 DAEs，这需要微分向量场内使用的代数求解过程。我们使用 JAX 和 diffrax (Kidger,2021) 实现了这一点，其中代数求解的自定义 JVP 基于隐函数定理，包括分段事件处理、复合守卫检测以及右连续目标时间评估。

第二种策略是显式离散伴随方法。它不是将通过模拟器的黑盒计算图进行微分，而是构建一个事件分割的离散残差系统。光滑段受 DAE 积分残差的约束，而事件块则执行守卫条件、重置方程、未改变变量的连续性以及事件后的代数一致性。残差不是作为启发式惩罚项添加的。它们是等式约束，或者等效地，在可行轨迹上消失的拉格朗日项。它们的乘子定义了伴随变量，产生的后向扫描计算参数梯度，而无需显式形成完整的轨迹和事件时间敏感性矩阵。

本文的贡献是对带有状态依赖事件的 DAE 的可微参数学习进行了统一处理。具体而言，我们：

1. 将混合 DAE 的参数学习表述为一个带有代数约束、事件守卫和重置映射的约束最小二乘问题；
2. 提出了一种通过模拟的 AD 实现，其中代数变量通过基于隐函数定理的自定义 JVP 进行微分，事件通过分段可微积分处理；
3. 推导了一种显式离散伴随 formulation，其中残差是等式约束，相应的乘子通过事件分割的后向扫描求解；
4. 在固定事件序列下，从梯度解释、事件时间处理、实现复杂性和局部有效性方面比较了这两种方法。

两种方法计算的梯度都局部于前向模拟所选的事件路径。它们假设横截事件穿越，并且在足够小的参数扰动下，事件的数量、身份或顺序不会发生变化。这些假设对于微分分段光滑混合轨迹是标准的，并阐明了基于梯度的优化在掠过事件、同时事件或事件序列变化附近的局限性。这种局部路径解释与混合系统中的灵敏度更新密切相关，通常通过事件时间的跳跃条件或盐化矩阵来表达 (Kong et al.,2024; Saccon et al.,2014)。这些条件描述了一阶扰动如何跨越状态触发的不连续性传播，并在混合系统灵敏度分析与本文研究的梯度 formulation 中出现的事件时间项之间提供了有用的概念联系。

我们将这两种方法应用于两个例子：一个带有事件的电路和一个平面弹跳球例子，其中事件由球的撞击触发。生成本文总结的实验结果的代码可在 https://github.com/ionmatei/diff-dae-events/ 找到。

论文结构：第 2 节将带有状态依赖事件的 DAE 的参数学习表述为一个约束最小二乘问题。第 3 节介绍了通过模拟的 AD 方法，包括代数求解的隐式微分和分段事件处理。第 4 节开发了基于事件分割残差约束的显式离散伴随 formulation。第 6 节将所提出的方法与相关的微分 ODE 和事件学习方法进行了比较。第 7 节总结了论文并概述了未来方向。附录包括额外的数学和实现细节。

## 2 问题表述

令 $0 = \tau_0 < \tau_1 < \dots < \tau_M < \tau_{M+1} = T$ 表示由混合 DAE 生成的事件时间。在每个区间 $(\tau_m, \tau_{m+1})$ 上，系统满足：

$$
\begin{aligned}
\dot{x}_m(t) &= f(t, x_m(t), z_m(t), p), \\
0 &= g(t, x_m(t), z_m(t), p).
\end{aligned}
$$

在事件时间 $\tau_m$，活跃的守卫满足 $\phi_{e_m}(\tau_m, x_{m-}, z_{m-}, p) = 0$，并且重置映射或重置方程关联事件前和事件后变量：

$$
(x_{m+}, z_{m+}) = \Psi_{e_m}(\tau_m, x_{m-}, z_{m-}, p),
$$

其中理解为在微分状态重置后，代数变量可能会被重新计算。

给定观察时间 $t_i$ 处的输出数据 $y_i^{\text{data}}$，参数估计问题可以写为：

$$
\begin{aligned}
&\min_{p, \{x_m, z_m\}, \{\tau_m\}} \frac{1}{2} \sum_{i=1}^{N_{\text{obs}}} \left\| h(t_i, x(t_i), z(t_i), p) - y_i^{\text{data}} \right\|_2^2 \\
\text{s.t.} \quad & \dot{x}_m(t) - f(t, x_m(t), z_m(t), p) = 0, \\
& g(t, x_m(t), z_m(t), p) = 0, \\
& E_m(\tau_m, w_{m-}, w_{m+}, p) = 0, \quad m=1, \dots, M,
\end{aligned}
$$

其中 $w = (x, z)$，且 $E_m$ 收集了守卫条件、重置方程、未被事件修改的变量的连续性条件以及事件后的代数一致性。

这种约束 formulation 是本文研究的两种方法的共同起点。通过模拟的 AD 方法通过执行模拟器并微分从参数到输出的缩减映射来消除轨迹变量。离散伴随方法在离散化后保持残差约束明确，并求解相应的拉格朗日乘子。

## 3 带有事件的 DAE 的通过模拟的自动微分

通过模拟的 AD 方法将事件驱动的 DAE 模拟器视为一个可微程序。对于选定的参数向量 $p_{\text{opt}}$，模拟器产生预测输出：

$$
\widehat{Y}(p_{\text{opt}}) = \begin{bmatrix} \hat{y}(t_1; p_{\text{opt}}) \\ \vdots \\ \hat{y}(t_N; p_{\text{opt}}) \end{bmatrix}.
$$

缩减后的学习问题为：

$$
\min_{p_{\text{opt}}} \widehat{\mathcal{J}}(p_{\text{opt}}) = \mathcal{J}(\widehat{Y}(p_{\text{opt}})),
$$

梯度通过自动微分计算：

$$
\nabla_{p_{\text{opt}}} \widehat{\mathcal{J}} = \text{AD}[\mathcal{J}(\widehat{Y}(p_{\text{opt}}))].
$$

主要的 DAE 特定问题是代数变量。在积分期间，$z$ 由 $g(t, x, z, p) = 0$ 隐式定义。我们将解写为 $z = \zeta(t, x, p)$，使得 $g(t, x, \zeta(t, x, p), p) = 0$。微分状态使用缩减的向量场推进：

$$
\dot{x} = f(t, x, \zeta(t, x, p), p).
$$

在实现中，代数求解通过弦牛顿迭代执行，并使用基于隐函数定理的自定义 JVP 进行微分。因此，梯度通过代数约束传播，而无需对每次牛顿迭代进行微分。

事件通过分段积分处理。在每个段的开始，评估所有守卫函数，并屏蔽已经活跃的守卫以避免立即重新触发。其余活跃守卫组合成一个标量事件条件：

$$
\Phi(t, x, p) = \min_{e \in \mathcal{A}} \phi_e(t, x, \zeta(t, x, p), p),
$$

其中 $\mathcal{A}$ 是活跃守卫的集合。根查找器定位该复合守卫的第一个零交叉点。选择相应的事件索引，应用重置映射，重新计算代数变量，下一段从事件后状态开始。

为了使此过程兼容 JIT 编译，模拟器被编写为对规定的最大数量事件边界段进行固定长度的扫描。物理段被记录直到达到终端时间；剩余的迭代是填充段。目标时间评估通过对每个记录的段重新积分完成，请求的目标时间被裁剪到该段。每个目标时间的输出从包含该目标的最新真实段的闭区间中选择，给出右连续惯例 $\hat{y}(\tau_m) = h(\tau_m, x(\tau_m^+), z(\tau_m^+), p)$。

产生的损失和梯度通过 JIT 编译的 `value_and_grad` 调用进行评估。这些梯度通过代数求解、ODE 积分、事件时间计算、重置映射、目标时间选择和输出重建传递。它们局部于前向模拟所选的事件序列。