当前位置: 首页 > news >正文

【读论文】DiffPhyCon 扩散物理系统控制

    • 摘要
    • Introduction
    • 2 Background
      • 2.1 问题设置
      • 2.2 预备知识:扩散模型
    • 3 Method
    • 3.1 扩散模型的生成控制
      • 3.2 先验重加权
    • 4 Experiments
      • 4.1 1D Burgers’ Equation Control
    • 5 Limitation
    • 6 Conclusion
    • A Additional Related work
      • A.1 Physical Systems Simulation
      • A.2 Physical Systems Control
      • A.3 Diffusion Models
    • B 1D Burgers’ Equation Visualization
    • C 附加信息:1D Burgers 方程的控制
      • C.1 数据生成
      • C.2 实验设置
      • C.2.1 部分观测,完全控制(PO-FC)
      • C.2.2 完全观测,部分控制
      • C.2.3 部分观测,部分控制
      • C.3 模型设计
      • C.3.1 DiffPhyCon-lite
      • C.3.2 DiffPhyCon
      • C.4 训练与评估
        • 训练
        • 推理
        • 评估

摘要

控制复杂物理系统的演变是科学和工程领域中的一项基础任务。传统技术通常因适用性有限或计算成本过高而受到限制。另一方面,近年来基于深度学习和强化学习的方法在满足系统动力学约束下优化长期控制序列时往往表现欠佳。在本研究中,我们提出了一种名为“扩散物理系统控制” (DiffPhyCon)的新方法来应对物理系统的控制问题。DiffPhyCon 的优势在于能够同时最小化学习到的生成能量函数和整个轨迹及控制序列中的预设控制目标。因此,它可以进行全局探索并识别接近最优的控制序列。此外,我们通过先验重加权增强了 DiffPhyCon,使其能够发现与训练分布显著偏离的控制序列。我们在一维伯格斯方程和二维水母在流体环境中的运动控制中测试了此方法。我们的方法优于广泛应用的传统方法以及最新的深度学习和强化学习方法。值得注意的是,DiffPhyCon 揭示了一种有趣的“快速闭合-缓慢打开”模式,该模式与流体动力学领域的既有研究结果一致。

Introduction

对复杂物理系统动力学的建模是科学和工程中的重要问题类别。通常,我们不仅关心预测物理系统的行为,还希望注入随时间变化的信号来引导其演变并优化特定目标。这形成了复杂物理控制问题,这是一项具有众多应用的基础任务, 包括受控核聚变 [9]、流体控制 [21]、水下设备 [72] 和航空 [48] 等领域。

尽管其重要性显而易见,但高效地控制复杂物理系统面临着重大的挑战。这类问题继承了模拟复杂物理系统的根本困难,这些系统通常是高维且高度非线性的(详见附录 A.1)。此外,用于优化控制模型的已观测控制信号及相应的系统轨迹通常远未达到特定控制目标的最优解。这个事实带来了一个显著的难题:如何在训练分布之外探索长期控制序列以寻找接近最优的解决方案,同时确保生成的系统轨迹符合物理系统的动力学?

为了解决物理系统控制问题,已有各种技术被提出,但这些方法未能应对上述挑战。对于传统控制方法,比例-积分-微分(PID)控制 [33] 效率较高,但仅适用于有限的问题范围。相比之下,尽管模型预测控制(MPC) [59] 的适用范围更广,但其计算成本高且在全局优化方面面临挑战。近期在监督学习(SL) [21, 23] 和强化学习(RL) [14, 47, 53] 方面的进展,通过系统轨迹和控制信号数据进行训练,在解决物理系统控制问题上展现了令人印象深刻的表现。然而,现有的 SL 和基于模型的 RL 方法往往陷入短视的失败模式 [25],未能实现长期的近最优解,或生成与物理系统动力学相违背的对抗性状态轨迹 [74]。

其主要原因可能在于,这些方法从迭代视角处理动力学的连续演变,缺乏长期视野并在全局优化方面存在困难。更多关于物理系统控制的相关工作见附录 A.2。

这些传统和现代方法大多采用 逐步预测或迭代 的方法,即 每一步都根据当前状态进行预测和调整,缺乏对整个系统演化的全局视野,导致在全局优化上表现不足,无法从整体上优化整个控制序列。

在本研究中,我们引入了“扩散物理系统控制”(DiffPhyCon),一种用于解决物理系统控制问题的新方法。我们 能量优化的角度出发,在整个时间范围内对系统轨迹和控制序列进行优化,以隐式捕获系统动力学中固有的约束。我们通过扩散模型实现这一点,该模型使用系统轨迹数据和控制序列进行训练。在推理阶段,DiffPhyCon 将模拟和控制优化整合为一个统一的能量优化过程。这 防止了生成的系统动力学超出分布,并为长期动力学提供了增强的视角,从而有助于发现优化目标的控制序列。

物理系统控制的一个关键方面在于其生成接近最优控制的能力,即使这些控制可能显著偏离训练分布。我们通过一个重要的见解解决这一挑战,即 学习到的生成能量景观可以分解为两个组成部分:代表控制序列的先验概率分布(即在训练数据中出现的控制序列的统计特性)和在给定控制序列时刻画系统轨迹的条件概率分布(体现系统在特定控制下如何演化)。基于这一见解,我们开发了一种 先验重加权技术,在推理过程中以可调强度从整体生成能量景观中减去控制序列的先验分布效应。

我们通过在一维伯格斯方程二维水母运动控制任务上的广泛实验展示了 DiffPhyCon 的有效性,这些任务还增加了部分观测和部分控制的挑战。特别是,二维水母运动控制问题为水下和空中设备的设计提供了重要的灵感,其挑战来自复杂的涡流行为和流固耦合动力学 [55, 27]。我们基于 CFD(计算流体动力学)软件生成了这样一个数据集,以模拟在控制信号作用下水母的拍动行为。 为促进复杂物理系统控制的研究,我们将该数据集作为基准贡献给学术界。在这两个任务中,我们的方法优于广泛应用的传统控制方法,并且在与最新的监督学习和强强化学习基线的比较中表现出色。值得注意的是,DiffPhyCon 揭示了水母表现出的“快速闭合-缓慢打开”有趣模式,这与流体动力学领域的研究结果一致 [27]。

总结来说,我们的贡献如下:(1) 开发了 DiffPhyCon,这是一种用于控制复杂物理系统的新型生成方法。通过在整个时间范围内联合优化轨迹和控制序列,DiffPhyCon 促进了长期动力学的全局优化,并有助于减少短视失败模式。(2) 引入了先验重加权技术,以生成优于训练集中控制序列的控制。(3) 在一维伯格斯方程和二维水母运动控制任务上展示了我们方法的有效性,尤其是在部分观测/控制场景中显示了优势。(4) 为复杂物理系统控制社区发布了一个新的水母运动控制基准。

2 Background

2.1 问题设置

我们考虑以下复杂物理系统的控制问题:

w ∗ = argmin w J ( u , w ) s.t. C ( u , w ) = 0. (1) \mathbf w^* = \underset{\mathbf w}{\text{argmin}} \, \mathcal{J}( \mathbf u, \mathbf w) \quad \text{s.t.} \quad \mathcal{C}(\mathbf u,\mathbf w) = 0. \tag{1} w=wargminJ(u,w)s.t.C(u,w)=0.(1)

  • u ( t , x ) : [ 0 , T ] × Ω ↦ R d u \mathbf u(t, \mathbf x) : [0, \mathcal T] \times \Omega \mapsto \mathbb{R}^{d_\mathbf u} u(t,x):[0,T]×ΩRdu 表示系统轨迹 { u ( t , ⋅ ) , t ∈ [ 0 , T ] } \{\mathbf u(t, \cdot), t \in [0, \mathcal T]\} {u(t,),t[0,T]},其定义在时间范围 [ 0 , T ] ⊂ R [0,\mathcal T] \subset \mathbb{R} [0,T]R 和空间域 Ω ⊂ R D \Omega \subset \mathbb{R}^D ΩRD 上;
  • w ( t , x ) : [ 0 , T ] × Ω ↦ R d w \mathbf w(t, \mathbf x) : [0, \mathcal T] \times \Omega \mapsto \mathbb{R}^{d_\mathbf w} w(t,x):[0,T]×ΩRdw外部控制信号,维度为 d w d_\mathbf w dw
  • J ( u , w ) \mathcal{J}(u, w) J(u,w) 表示控制目标。例如,可以设计 J \mathcal{J} J 来衡量系统向目标状态 u ∗ \mathbf u^* u 演化的控制性能,具有如下代价约束: J = ∫ ∥ u − u ∗ ∥ 2 d x d t + ∫ ∥ w ∥ 2 d x d t . \mathcal{J} = \int \|\mathbf u - \mathbf u^*\|^2 \, \text{d}\mathbf x \text{d}t + \int \|\mathbf w\|^2 \,\text{d}\mathbf x \text{d}t. J=uu2dxdt+w2dxdt.
  • C ( u , w ) = 0 \mathcal{C}(\mathbf u, \mathbf w) = 0 C(u,w)=0 表示物理约束,可以 通过描述控制信号 w ( t , x ) \mathbf w(t, \mathbf x) w(t,x) 如何在边界和初始条件下驱动系统轨迹 u ( t , x ) \mathbf u(t, \mathbf x) u(t,x) 的偏微分方程显式指定,或者 通过从物理系统观测中收集的控制序列和轨迹数据隐式指定。在后一种情况下,我们可能只能获取部分轨迹数据或进行部分控制。这些情形共同对物理系统控制任务提出了显著挑战。

2.2 预备知识:扩散模型

扩散模型 [18] 是一类从数据中学习数据分布的生成模型。扩散模型由两个相反的过程组成:

  • 正向过程 q ( x k + 1 ∣ x k ) = N ( x k + 1 ; α k x k , ( 1 − α k ) I ) q(x_{k+1} | x_k) = \mathcal{N}(x_{k+1}; \sqrt{\alpha_k} x_k, (1 - \alpha_k) \mathbf I) q(xk+1xk)=N(xk+1;αk xk,(1αk)I) 将干净数据 x 0 \mathbf x_0 x0 扰乱成高斯噪声 x K ∼ N ( 0 , I ) \mathbf x_K \sim \mathcal{N}(0, \mathbf I) xKN(0,I)

  • 以及逆向参数化过程 p θ ( x k − 1 ∣ x k ) = N ( x k − 1 ; μ θ ( x k , k ) , σ k I ) p_\theta(x_{k-1} | x_k) = \mathcal{N}(x_{k-1}; \mu_\theta(x_k, k), \sigma_k \mathbf I) pθ(xk1xk)=N(xk1;μθ(xk,k),σkI) 从标准高斯噪声 x K ∼ N ( 0 , I ) \mathbf x_K \sim \mathcal{N}(0, \mathbf I) xKN(0,I) 去噪。

  • 其中 { α k } k = 1 K \{\alpha_k\}_{k=1}^K {αk}k=1K 是方差调度。

  • 为了训练扩散模型,[18] 提出了 DDPM 方法,通过 最小化去噪网络 ϵ θ \epsilon_\theta ϵθ 的损失函数,该损失是数据对数似然的证据下界(ELBO)的简化版本
    L = E k ∼ U ( 1 , K ) , x 0 ∼ p ( x ) , ϵ ∼ N ( 0 , I ) [ ∥ ϵ − ϵ θ ( α ˉ k x 0 + 1 − α ˉ k ϵ , k ) ∥ 2 2 ] , (2) \mathcal L = \mathbb{E}_{k \sim U(1, K), \mathbf x_0 \sim p(x), \epsilon \sim \mathcal{N}(0, \mathbf I)} \left[ \| \epsilon - \epsilon_\theta(\sqrt{\bar{\alpha}_k} \mathbf x_0 + \sqrt{1 - \bar{\alpha}_k} \epsilon, k) \|_2^2 \right], \tag{2} L=EkU(1,K),x0p(x),ϵN(0,I)[ϵϵθ(αˉk x0+1αˉk ϵ,k)22],(2)

    其中 α ˉ k : = ∏ i = 1 k α i \bar{\alpha}_k := \prod_{i=1}^k \alpha_i αˉk:=i=1kαi ϵ θ \epsilon_\theta ϵθ 估计需要移除的噪声以恢复数据 x 0 \mathbf x_0 x0

  • 在推理过程中,通过从高斯噪声开始迭代应用 ϵ θ \epsilon_\theta ϵθ,可以生成一个近似遵循数据分布 p ( x ) p(\mathbf x) p(x) 的新样本 x 0 \mathbf x_0 x0。更多关于扩散模型的相关研究见附录 A.3。

符号说明:用 v [ n , m ] = [ v n , … , v m ] \mathbf{v}_{[n,m]} = [\mathbf v_n, \ldots,\mathbf v_m] v[n,m]=[vn,,vm] **表示变量序列。用 z [ 0 , T − 1 ] , k \mathbf z_{[0, T-1], k} z[0,T1],k 表示扩散步骤 k k k 中的隐变量 z [ 0 , T − 1 ] \mathbf z_{[0, T-1]} z[0,T1]。为简便起见,将 w [ 0 , T − 1 ] , u [ 1 , T ] \mathbf w_{[0, T-1]}, \mathbf u_{[1, T]} w[0,T1],u[1,T] 简写为 w , u \mathbf w, \mathbf u w,u。变量的拼接用例如 [ u , w ] [\mathbf u, \mathbf w] [u,w] 表示。

3 Method

在本节中,我们详细介绍了我们的方法 DiffPhyCon。在第 3.1 节中,我们介绍了我们的方法,包括其训练和推理。在第 3.2 节中,我们进一步提出了一种先前的重新加权技术来改进 DiffPhyCon。DiffPhyCon 的概览如图 1 所示。

在这里插入图片描述

图 1:DiffPhyCon 概述。图中展示了 DiffPhyCon 的训练过程(顶部)、推理过程(左下角)和评估过程(右下角)。橙色和蓝色分别代表学习联合分布 p θ ( u , w ) p_\theta(\mathbf u,\mathbf w) pθ(u,w) 和先验分布 p ϕ ( w ) p_\phi(\mathbf w) pϕ(w) 的模型。通过先验重加权和引导,DiffPhyCon 能够生成更优的控制序列。

3.1 扩散模型的生成控制

DiffPhyCon 采用能量优化的角度来解决方程 (1) 中的问题,其中 偏微分方程(PDE)约束可建模为参数化的基于能量的模型 (EBM) E θ ( u , w , c ) E_{\theta}(\mathbf u,\mathbf w, \mathbf c) Eθ(u,w,c),它通过关系 p ( u , w ∣ c ) ∝ exp ⁡ ( − E θ ( u , w , c ) ) p(\mathbf u, \mathbf w|\mathbf c) \propto \exp(-E_{\theta}(\mathbf u,\mathbf w,\mathbf c)) p(u,wc)exp(Eθ(u,w,c)) 来描述条件 c \mathbf c c u \mathbf u u w \mathbf w w 的分布。较低的 E θ ( u , w , c ) E_{\theta}(\mathbf u,\mathbf w,\mathbf c) Eθ(u,w,c) (或等价地,较高的 p ( u , w ∣ c ) p(\mathbf u, \mathbf w|\mathbf c) p(u,wc) )表示对 PDE 约束的满足程度更好。

因此,方程 (1) 可转换为:

u ∗ , w ∗ = arg ⁡ min ⁡ u , w [ E θ ( u , w , c ) + λ ⋅ J ( u , w ) ] (3) \mathbf u^*, \mathbf w^* = \arg\min_{\mathbf u, \mathbf w} [E_{\theta}(\mathbf u,\mathbf w,\mathbf c) + \lambda \cdot \mathcal{J}(\mathbf u, \mathbf w)] \tag{3} u,w=argu,wmin[Eθ(u,w,c)+λJ(u,w)](3)

其中, λ \lambda λ 是一个超参数。该公式同时优化所有物理时间步的 u \mathbf u u w \mathbf w w

  • 第一项促使生成的 w \mathbf w w u \mathbf u u 满足 PDE 约束 C ( u , w ) = 0 \mathcal C(\mathbf u, \mathbf w) = 0 C(u,w)=0
  • 第二项引导优化朝着最优目标前进。(表示当前轨迹或控制序列下的目标函数值

该优化框架的优势在于,它 对所有时间步的 u \mathbf u u w \mathbf w w 进行全局优化,从而获得更符合物理系统动力学的解。有关超参数 λ \lambda λ 影响的详细信息,请参见附录 J。

训练过程:为了训练 E θ E_{\theta} Eθ,我们利用扩散模型来估计其梯度。为简化表示,我们引入新变量 z \mathbf z z 表示 u \mathbf u u w \mathbf w w 的合并,即 z = [ u , w ] \mathbf z = [\mathbf u, \mathbf w] z=[u,w]。我们使用去噪网络 ϵ θ \epsilon_{\theta} ϵθ [18],它近似 ∇ E θ ( z , c ) \nabla E_{\theta}(\mathbf z, \mathbf c) Eθ(z,c) [12],以学习每个扩散步骤 k = 1 , … , K k = 1, \ldots, K k=1,,K 中应去除的噪声。去噪网络 ϵ θ \epsilon_{\theta} ϵθ 的训练损失函数为:

L = E k ∼ U ( 1 , K ) , ( z , c ) ∼ p ( z , c ) , ϵ ∼ N ( 0 , I ) [ ∥ ϵ − ϵ θ ( α ˉ k z + 1 − α ˉ k ϵ , c , k ) ∥ 2 2 ] , (4) \mathcal{L} = \mathbb{E}_{k \sim {U}(1, K), (\mathbf z,\mathbf c) \sim p(\mathbf z,\mathbf c), \epsilon \sim \mathcal{N}(0, \mathbf I)} [\| \epsilon - \epsilon_{\theta} (\sqrt{\bar{\alpha}_k} \mathbf z + \sqrt{1 - \bar{\alpha}_k} \epsilon, \mathbf c, k) \|_2^2], \tag{4} L=EkU(1,K),(z,c)p(z,c),ϵN(0,I)[ϵϵθ(αˉk z+1αˉk ϵ,c,k)22],(4)

其中 ϵ θ \epsilon_{\theta} ϵθ 依赖于 c \mathbf c c。对于训练数据集,如果存在明确的 PDE 形式,可以通过设计条件 c \mathbf c c 和控制序列并进行仿真来生成数据集。否则,应从观测到的控制序列和轨迹对中收集数据。

控制优化:在训练完成去噪网络后,可以使用朗之万采样过程来优化方程 (3)。我们从初始样本 z K ∼ N ( 0 , I ) \mathbf z_K \sim \mathcal{N}(0, \mathbf I) zKN(0,I) 开始,并迭代运行以下过程:

z k − 1 = z k − η ∇ z ( E θ ( z k , c ) + λ J ( z ^ k ) ) + ξ , ξ ∼ N ( 0 , σ k 2 I ) (5) \mathbf z_{k-1} = \mathbf z_k - \eta \nabla_\mathbf z (E_{\theta}(\mathbf z_k, \mathbf c) + \lambda \mathcal{J}(\hat{\mathbf z}_k) ) + \xi, \quad \xi \sim \mathcal{N}(0, \sigma_k^2 \mathbf I) \tag{5} zk1=zkηz(Eθ(zk,c)+λJ(z^k))+ξ,ξN(0,σk2I)(5)

其中 σ k 2 \sigma_k^2 σk2 η \eta η 分别对应于扩散过程中的噪声调度和缩放因子。这里, z ^ k \hat{\mathbf z}_k z^k 是通过以下公式从 z k \mathbf z_k zk 估计出的近似无噪声的 z 0 \mathbf z_0 z0

z ^ k = ( z k − 1 − α ˉ k ϵ θ ( z k , c , k ) ) / α ˉ k . (6) \hat{\mathbf z}_k = (\mathbf z_k - \sqrt{1 - \bar \alpha_k} \epsilon_{\theta}(\mathbf z_k,\mathbf c, k)) / \sqrt{\bar \alpha_k}. \tag{6} z^k=(zk1αˉk ϵθ(zk,c,k))/αˉk .(6)

我们基于 z ^ k \hat{\mathbf z}_k z^k 而非直接使用 z k \mathbf z_k zk 计算 J \mathcal{J} J,因为 z k \mathbf z_k zk 中的噪声可能会引入误差。然后,用训练好的去噪网络 ϵ θ \epsilon_{\theta} ϵθ 替换 ∇ E θ \nabla E_{\theta} Eθ 如下:

z k − 1 = z k − η ( ϵ θ ( z k , c , k ) + λ ∇ z J ( z ^ k ) ) + ξ , ξ ∼ N ( 0 , σ k 2 I ) . (7) \mathbf z_{k-1} = \mathbf z_k - \eta (\epsilon_{\theta} (\mathbf z_k, \mathbf c, k) + \lambda \nabla_\mathbf z \mathcal{J}(\hat{\mathbf z}_k)) + \xi, \quad \xi \sim \mathcal{N}(0, \sigma_k^2 \mathbf I). \tag{7} zk1=zkη(ϵθ(zk,c,k)+λzJ(z^k))+ξ,ξN(0,σk2I).(7)

该去噪过程的迭代从 k = K , K − 1 , … , 1 k = K, K - 1, \ldots, 1 k=K,K1,,1 进行,最终得到问题 (3) 的解 z 0 = { u [ 1 , T ] , 0 , w [ 0 , T − 1 ] , 0 } \mathbf z_0 = \{\mathbf u_{[1, T], 0},\mathbf w_{[0, T-1], 0}\} z0={u[1,T],0,w[0,T1],0}

引导条件:除了上述引入的显式引导外,条件也被广泛用于引导扩散模型中的采样 [19, 60]。当控制目标可以自然地以条件形式表达时,例如要求生成的轨迹 u \mathbf u u 与期望的目标轨迹 u ∗ \mathbf u^* u 一致,我们可以 u = u ∗ \mathbf u = \mathbf u^* u=u 作为条件包含在 c \mathbf c c 中,使得从扩散模型中采样的轨迹 u \mathbf u u 自动满足 u = u ∗ \mathbf u = \mathbf u^* u=u。总体而言,在我们提出的 DiffPhyCon 框架中,控制目标 J \mathcal J J 可以通过显式引导 ∇ J \nabla \mathcal J J 或引导条件来优化,这取决于具体的控制目标。

3.2 先验重加权

动机:在物理系统控制中,一个关键的挑战是获得优于训练数据集中的控制序列,这些控制序列通常偏离了实现最优控制目标的路径。尽管在我们的扩散模型中已经融入了控制目标的引导生成的控制序列仍然受到训练数据集中控制序列先验分布的强烈影响。这启发我们探索减轻该先验影响的策略,旨在生成近优的控制序列。

我们通过以下关键洞察来解决这一挑战:基于能量的模型 E ( u , w , c ) E(\mathbf u, \mathbf w,\mathbf c) E(u,w,c) 可以分解为两个组成部分:一个是由控制序列的先验分布 p ( w ∣ c ) p(\mathbf w|\mathbf c) p(wc) 导出的 E ( p ) ( w , c ) E^{(p)}(\mathbf w,\mathbf c) E(p)(w,c),另一个是表示给定控制序列条件下轨迹的条件概率分布 p ( u ∣ w , c ) p(\mathbf u|\mathbf w,\mathbf c) p(uw,c) E ( c ) ( u , w , c ) E^{(c)}(\mathbf u,\mathbf w,\mathbf c) E(c)(u,w,c)

通过概率分布的乘法性质: p ( u , w ∣ c ) = p ( w ∣ c ) p ( u ∣ w , c ) p(\mathbf u, \mathbf w|\mathbf c) = p(\mathbf w|\mathbf c) p(\mathbf u|\mathbf w, \mathbf c) p(u,wc)=p(wc)p(uw,c) ,结合能量函数和概率分布的关系 p ( x ) ∝ exp ⁡ ( − E ( x ) ) p(\mathbf x) \propto \exp(-E(\mathbf x)) p(x)exp(E(x)),因此能量的分解为:

E ( u , w , c ) = E ( p ) ( w , c ) + E ( c ) ( u , w , c ) , (8) E(\mathbf u,\mathbf w, \mathbf c) = E^{(p)}(\mathbf w, \mathbf c) + E^{(c)}(\mathbf u,\mathbf w, \mathbf c), \tag{8} E(u,w,c)=E(p)(w,c)+E(c)(u,w,c),(8)

我们提出了一种先验重加权技术,引入一个可调的超参数 γ > 0 \gamma > 0 γ>0,作为 p ( w ∣ c ) p(\mathbf w|\mathbf c) p(wc) 的指数,用以调整先验分布的影响。然后我们得到一个重加权版本的 p ( u , w ∣ c ) p(\mathbf u, \mathbf w|\mathbf c) p(u,wc),即

p γ ( u , w ∣ c ) = p ( w ∣ c ) γ p ( u ∣ w , c ) / Z p_\gamma(\mathbf u, \mathbf w|\mathbf c) = p(\mathbf w|\mathbf c)^\gamma p(\mathbf u|\mathbf w,\mathbf c) / Z pγ(u,wc)=p(wc)γp(uw,c)/Z

,这仍然是一个概率分布,并根据 p ( u ∣ w , c ) = p ( u , w ∣ c ) p ( w ∣ c ) p(\mathbf u|\mathbf w,\mathbf c) = \frac{p(\mathbf u, \mathbf w|\mathbf c)}{p(\mathbf w|\mathbf c)} p(uw,c)=p(wc)p(u,wc) 进一步转化为

p γ ( u , w ∣ c ) = p ( w ∣ c ) γ − 1 p ( u , w ∣ c ) / Z , (9) p_\gamma(\mathbf u, \mathbf w|\mathbf c) = p(\mathbf w|\mathbf c)^{\gamma-1} p(\mathbf u,\mathbf w|\mathbf c) / Z, \tag{9} pγ(u,wc)=p(wc)γ1p(u,wc)/Z,(9)

其中 Z Z Z 是一个归一化常数。特别地,当 0 < γ < 1 0 < \gamma < 1 0<γ<1 时,这种方法具有优势,因为它压平了分布 p ( u , w ∣ c ) p(\mathbf u, \mathbf w|\mathbf c) p(u,wc),从而增加了从 p ( u , w ∣ c ) p(\mathbf u, \mathbf w|\mathbf c) p(u,wc) 低概率区域采样的可能性,而该区域可能包含问题 Eq. (3) 的最优解。将 Eq. (9) 代入 Eq. (8),我们得到

E ( γ ) ( u , w , c ) = ( γ − 1 ) E ( p ) ( w , c ) + E θ ( u , w , c ) − log ⁡ Z , (10) E^{(\gamma)}(\mathbf u,\mathbf w,\mathbf c) = (\gamma - 1) E^{(p)}(\mathbf w,\mathbf c) + E_\theta(\mathbf u,\mathbf w,\mathbf c) - \log Z, \tag{10} E(γ)(u,w,c)=(γ1)E(p)(w,c)+Eθ(u,w,c)logZ,(10)

其中 E ( γ ) ( u , w , c ) = − log ⁡ ( p γ ( u , w ∣ c ) ) + c o n s t E^{(\gamma)}(\mathbf u,\mathbf w,\mathbf c) = -\log(p_\gamma(\mathbf u, \mathbf w|\mathbf c)) + const E(γ)(u,w,c)=log(pγ(u,wc))+const 是与 Eq. (3) 中的 E θ ( u , w , c ) E_\theta(\mathbf u,\mathbf w,\mathbf c) Eθ(u,w,c) 相关的重加权能量基础模型,依赖于超参数 γ \gamma γ

然后,优化问题 Eq. (3) 可以转化为

u ∗ , w ∗ = arg min ⁡ u , w [ E ( γ ) ( u , w , c ) + λ ⋅ J ( u , w ) ] . (11) \mathbf u^*, \mathbf w^* = \argmin_{\mathbf u,\mathbf w} \left[ E^{(\gamma)}( \mathbf u, \mathbf w,\mathbf c) + \lambda \cdot \mathcal J(\mathbf u, \mathbf w) \right]. \tag{11} u,w=u,wargmin[E(γ)(u,w,c)+λJ(u,w)].(11)

该问题的优化鼓励从 p ( w ∣ c ) p(\mathbf w| \mathbf c) p(wc) 的低概率区域进行采样,同时最小化控制目标,这具备生成比原优化问题中 γ = 1 \gamma = 1 γ=1 的退化版本更接近最优的控制序列的能力。先验重加权的直觉如图 2 所示。

在这里插入图片描述
图 2:先验重加权的直观理解。顶部表面展示了 J ( u , w ) \mathcal J(\mathbf u, \mathbf w) J(u,w) 的景观,其中高维变量 u \mathbf u u w \mathbf w w 用一维表示。中间和平面下方展示了重加权分布 p γ ( w ) p ( u ∣ w ) / Z p^\gamma(\mathbf w)p(\mathbf u|\mathbf w)/Z pγ(w)p(uw)/Z 的概率热图。通过调整 γ \gamma γ γ = 1 \gamma = 1 γ=1(中间平面)到 0 < γ < 1 0 < \gamma < 1 0<γ<1(下方平面),更好的最小值 J \mathcal J J(下方平面的红点)获得了被采样的机会。这与中间平面中受到先验 p ( w ) p(\mathbf w) p(w) 强烈影响的次优红点形成对比。

训练:为了学习重加权能量 E ( γ ) ( u , w , c ) E^{(\gamma)}(\mathbf u,\mathbf w,\mathbf c) E(γ)(u,w,c),我们通过对 Eq. (10) 两边求梯度,将其梯度参数化为两部分的和:

∇ E ϕ , θ ( γ ) ( u , w , c ) = ( γ − 1 ) ∇ E ϕ ( p ) ( w , c ) + ∇ E θ ( u , w , c ) , (12) \nabla E^{(\gamma)}_{\phi, \theta}(\mathbf u,\mathbf w,\mathbf c) = (\gamma - 1) \nabla E^{(p)}_\phi(\mathbf w, \mathbf c) + \nabla E_\theta(\mathbf u, \mathbf w,\mathbf c), \tag{12} Eϕ,θ(γ)(u,w,c)=(γ1)Eϕ(p)(w,c)+Eθ(u,w,c),(12)

其中 E ϕ ( p ) ( w , c ) E^{(p)}_\phi(\mathbf w, \mathbf c) Eϕ(p)(w,c) **参数化了对应于 p ( w ∣ c ) p(\mathbf w|\mathbf c) p(wc) 的能量基础模型 E ( p ) ( w , c ) E^{(p)}(\mathbf w,\mathbf c) E(p)(w,c)。注意, ∇ log ⁡ Z \nabla \log Z logZ 在此处为常数并消失。需要注意的是, ** ∇ E θ ( u , w , c ) \nabla E_\theta(\mathbf u,\mathbf w,\mathbf c) Eθ(u,w,c) 已通过 Eq. (4) 进行训练。 ∇ E ϕ ( p ) ( w , c ) \nabla E^{(p)}_\phi(\mathbf w,\mathbf c) Eϕ(p)(w,c) 可以通过以下损失函数类似地进行训练:

L = E k ∼ U ( 1 , K ) , ( w , c ) ∼ p ( w , c ) , ϵ ∼ N ( 0 , I ) [ ∥ ϵ − ϵ ϕ ( α ˉ t w + 1 − α ˉ k ϵ , c , k ) ∥ 2 2 ] , (13) \mathcal L = \mathbb{E}_{k \sim U(1, K), (\mathbf w,\mathbf c) \sim p(\mathbf w,\mathbf c), \epsilon \sim \mathcal N(0, \mathbf I)} \left[ \|\epsilon - \epsilon_\phi(\sqrt{\bar \alpha_t} \mathbf w + \sqrt{1 - \bar \alpha_k} \epsilon, \mathbf c, k)\|_2^2 \right], \tag{13} L=EkU(1,K),(w,c)p(w,c),ϵN(0,I)[ϵϵϕ(αˉt w+1αˉk ϵ,c,k)22],(13)

其中 ϵ ϕ \epsilon_\phi ϵϕ条件去噪网络,用于逼近 ∇ E ϕ ( p ) ( w , c ) \nabla E^{(p)}_\phi(\mathbf w, \mathbf c) Eϕ(p)(w,c)

控制优化:当 ϵ θ \epsilon_\theta ϵθ ϵ ϕ \epsilon_\phi ϵϕ 都经过训练后,可以通过以下迭代更新 Eq. (11):

z k − 1 = z k − η ( ϵ θ ( z k , c , k ) + λ ∇ z J ( z ^ k ) ) + ξ 1 , ξ 1 ∼ N ( 0 , σ k 2 I ) ( 14 ) w k − 1 = w k − 1 − η ( γ − 1 ) ϵ ϕ ( w k , c , k ) + ξ 2 , ξ 2 ∼ N ( 0 , σ k 2 I ) ( 15 ) \mathbf z_{k-1} =\mathbf z_k - \eta (\epsilon_\theta(\mathbf z_k,\mathbf c, k) + \lambda \nabla_\mathbf z \mathcal J(\hat{\mathbf z}_k)) + \xi_1, \quad \xi_1 \sim \mathcal N(0, \sigma_k^2 \mathbf I) \quad (14) \\ \mathbf w_{k-1} = \mathbf w_{k-1} - \eta (\gamma - 1) \epsilon_\phi(\mathbf w_k,\mathbf c, k) + \xi_2, \quad \xi_2 \sim \mathcal N(0, \sigma_k^2 \mathbf I) \quad (15) zk1=zkη(ϵθ(zk,c,k)+λzJ(z^k))+ξ1,ξ1N(0,σk2I)(14)wk1=wk1η(γ1)ϵϕ(wk,c,k)+ξ2,ξ2N(0,σk2I)(15)

其中 z k = [ u k , w k ] \mathbf z_k = [\mathbf u_k,\mathbf w_k] zk=[uk,wk]。与 Eq. (7) 的迭代方案不同,这里使用额外的步骤来更新 w k \mathbf w_k wk,基于 ϵ ϕ \epsilon_\phi ϵϕ 预测的噪声。这引导 z k = [ u k , w k ] \mathbf z_k = [\mathbf u_k,\mathbf w_k] zk=[uk,wk] 朝着重加权分布 p γ ( u , w ∣ c ) p_\gamma(\mathbf u, \mathbf w|\mathbf c) pγ(u,wc) 移动,同时在迭代中调整 z k \mathbf z_k zk 的方向以最小化目标函数。

完整的算法见 Algorithm 1。关于如何设置超参数 γ \gamma γ 的详细讨论和结果,请参见附录 J。

在这里插入图片描述

DiffPhyCon-lite:DiffPhyCon 中引入的先验重加权技术涉及训练和评估两个模型,因此带来了额外的计算开销。值得欣慰的是,我们可以通过调整参数 γ \gamma γ 来平衡 DiffPhyCon 的控制性能和计算开销。 γ = 1 \gamma = 1 γ=1 时,不再需要模型 ϵ ϕ \epsilon_\phi ϵϕ,我们将这种简化版本的 DiffPhyCon 称为 DiffPhyCon-lite

4 Experiments

在本节中,我们旨在回答以下问题:(1)DiffPhyCon 是否能在物理系统控制上优于传统的监督学习和强化学习方法?(2)所提出的先验重加权技术是否有助于实现更好的控制目标?(3)对于更具挑战性的部分观测或部分控制场景,(1) 和 (2) 的答案是否可以推广?为了解答这些问题,我们在 1D Burgers 方程和 2D 水母运动控制问题上进行了实验,这两个问题都是至关重要且具有挑战性的。

以下的最先进控制方法被选作基准。对于 1D Burgers 方程,我们使用:(1)经典且广泛使用的控制算法比例-积分-微分(PID)[33] 与我们训练的代理模型交互;(2)监督学习方法(SL)[23];包括:(3)软演员评论家(SAC)[15] 的离线版本和代理求解器版本;(4)行为克隆(BC)[50];以及(5)行为近端策略优化(BPPO)[75]。具体来说,SAC 的代理求解器版本与我们训练的代理模型交互,而离线版本只使用给定数据。BC 和 BPPO 也为离线版本。对于 2D 水母运动控制,基准方法包括 SL、SAC(离线)、SAC(代理求解器)、BC、BPPO,以及另一个经典的多输入多输出算法模型预测控制(MPC)[59]。PID 不适用于此数据驱动任务。基准方法的详细描述见附录 F 和附录 G。我们提供了代码,以便重现我们的方法和基准方法的结果。

4.1 1D Burgers’ Equation Control

实验设置:Burgers 方程是许多物理系统中出现的控制方程。我们考虑的是具有狄利克雷边界条件外部力 w ( t , x ) \mathbf w(t, x) w(t,x)1D Burgers 方程,如 [23, 43] 中所研究的:

{ ∂ u ∂ t = − u ⋅ ∂ u ∂ x + ν ∂ 2 u ∂ x 2 + w ( t , x ) in  [ 0 , T ] × Ω u ( t , x ) = 0 in  [ 0 , T ] × ∂ Ω u ( 0 , x ) = u 0 ( x ) in  { t = 0 } × Ω (16) \begin{cases} \frac{\partial u}{\partial t} = -u \cdot \frac{\partial u}{\partial x} + \nu \frac{\partial^2 u}{\partial x^2} + w(t, x) \quad & \text{in}\ [0, T] \times \Omega \\ u(t, x) = 0 \quad & \text{in}\ [0, T] \times \partial\Omega \\ u(0, x) = u_0(x) \quad & \text{in} \ \{t = 0\} \times \Omega\end{cases} \tag{16} tu=uxu+νx22u+w(t,x)u(t,x)=0u(0,x)=u0(x)in [0,T]×Ωin [0,T]×Ωin {t=0}×Ω(16)

这里, ν \nu ν 是粘性参数, u 0 ( x ) u_0(\mathbf x) u0(x) 是初始条件。根据方程 (16),给定目标状态 u d ( x ) u_d(x) ud(x),控制目标是最小化最终时刻 u T u_T uT 与目标状态 u d u_d ud 之间的控制误差 J actual \mathcal J_{\text{actual}} Jactual,同时约束控制序列 w ( t , x ) \mathbf w(t, x) w(t,x) 的能量成本 J energy \mathcal J_{\text{energy}} Jenergy

J actual : = ∫ Ω ∣ u ( T , x ) − u d ( x ) ∣ 2 d x , J energy : = ∫ [ 0 , T ] × Ω ∣ w ( t , x ) ∣ 2 d t d x (17) \mathcal J_{\text{actual}} := \int_{\Omega} |u(T, x) - u_d(x)|^2 \, dx, \quad \mathcal J_{\text{energy}} := \int_{[0,T] \times \Omega} |w(t, x)|^2 \, dt dx \tag{17} Jactual:=Ωu(T,x)ud(x)2dx,Jenergy:=[0,T]×Ωw(t,x)2dtdx(17)

为了使评估更具挑战性,我们选择了三种不同的实验设置,分别对应不同的现实场景:

  • 部分观测,完全控制 (PO-FC);
  • 完全观测,部分控制 (FO-PC);
  • 部分观测,部分控制 (PO-PC),

这些设置在附录 C.2 中进行了详细说明,并在附录 B 中进行了说明。这些设置对传统控制方法(如 PID)构成挑战,因为它们需要捕捉系统动力学中的长程依赖性。请注意,不同设置下报告的指标不能直接比较。

在此实验中,DiffPhyCon 使用了指导条件(第 3.1 节)来优化 J actual \mathcal J_{\text{actual}} Jactual,并使用显式指导来优化能量成本,具体内容见附录 C。

在这里插入图片描述
表格 1:1D Burgers 方程控制中实现的最佳 J a c t u a l \mathcal J_{actual} Jactual。加粗字体表示最佳模型,下划线表示第二最佳模型。

结果:在表 1 中,我们报告了不同方法的控制误差 J actual \mathcal J_{\text{actual}} Jactual 结果。可以观察到,DiffPhyCon 在所有基准方法中提供了最佳的结果。具体来说,在 PO-FC、FO-PC 和 PO-PC 设置下,DiffPhyCon 分别将最佳基准方法的 J actual \mathcal J_{\text{actual}} Jactual 减少了 30.1%、52.6% 和 44.6%。从表 1 可以看出,DiffPhyCon 和 DiffPhyCon-lite 之间的性能差距很小。这是因为,我们的 DiffPhyCon-lite 学习到的 w \mathbf w w 的先验分布是以 u 0 u_0 u0 u T u_T uT 为条件的,这完全决定了最优的 w \mathbf w w。因此, p ( w ∣ u 0 , u T ) p(w|u_0, u_T) p(wu0,uT) 本质上就是最优的分布,因此 DiffPhyCon-lite 已经能够提供令人满意的性能。

为了比较不同方法在约束能量成本 J energy \mathcal J_{\text{energy}} Jenergy 下优化 J actual \mathcal J_{\text{actual}} Jactual 的能力,我们在图 3 中比较了不同方法的帕累托前沿。我们通过调整超参数 λ \lambda λ 来控制 J actual \mathcal J_{\text{actual}} Jactual 和能量成本之间的权衡,因为大多数基准方法都具有这个超参数。如图 3 所示,DiffPhyCon 的帕累托前沿始终处于最佳水平,在大多数能量预算设置下实现了最低的 J actual \mathcal J_{\text{actual}} Jactual。尽管在完全观测设置 (b) 下,SL 表现良好,因为系统动态可以更容易预测,但它在部分观测场景 (a)© 中遇到了困难。结果证明,DiffPhyCon 相比基准方法具有生成近似最优控制序列的能力。更多可视化结果请参见附录 B。有关效率评估,请参考附录 I 中的表 19。

在这里插入图片描述
图 3:1D Burgers 方程中不同方法的 J e n e r g y \mathcal J_{energy} Jenergy J a c t u a l \mathcal J_{actual} Jactual 的 Pareto frontier。

5 Limitation

我们的方法 DiffPhyCon 仍然存在一些局限性,这为未来的工作提供了激动人心的机会。首先,我们的方法是数据驱动的,因此无法保证实现最优解,并且缺乏对生成的控制序列与最优序列之间误差差距的估计。其次,DiffPhyCon 的训练目前是在离线方式下进行的,未与真实解算器进行交互。将解算器纳入训练框架可以促进实时反馈,使得模型能够动态适应环境,并发现新颖的策略和解决方案。此外,我们提出的 DiffPhyCon 目前以开环方式运行,因为它没有考虑来自解算器的实时反馈。整合这种反馈将使算法能够根据环境的不断变化来调整其控制决策。

6 Conclusion

在这项工作中,我们提出了 DiffPhyCon,一种用于控制复杂物理系统的新方法。它通过联合优化生成的能量和控制目标来生成控制序列和状态轨迹。我们进一步引入了先验重加权技术,以便发现那些与训练数据显著不同的控制序列。通过全面的实验,我们展示了我们的方法在挑战性的物理系统控制任务中,相较于传统方法、深度学习和强化学习基准的优越性能。我们在附录 L 中讨论了未来的工作,并在附录 M 中阐述了社会影响。

A Additional Related work

A.1 Physical Systems Simulation

复杂物理系统的模拟是系统控制的基础。尽管经典的物理系统模拟数值技术以其准确性著称,但通常伴随着巨大的计算开销[42, 30]。最近,基于神经网络的求解器在加速模拟方面表现出相对于经典求解器的显著优势。这些方法大致可以分为三类:数据驱动方法[34, 58, 49, 5, 69, 4, 29]、物理信息神经网络(PINNs)[54, 7, 67]和求解器内循环方法[63, 65]。大多数方法采用迭代的横向预测框架。与此不同,我们将系统轨迹视为一个整体变量,并使用扩散模型来学习一个基于控制序列的显式模拟器。值得注意的工作是[6],该工作引入了扩散模型进行时间预测。虽然我们和[6]的工作都采用了扩散模型,但我们处理的是物理系统控制这一不同任务。此外,我们将控制目标纳入推理,并引入先验重加权来调整先验的影响

A.2 Physical Systems Control

对于由偏微分方程(PDE)描述的物理系统,伴随方法[36, 39, 52]在过去几十年中一直是系统控制的最广泛应用方法。它虽然准确,但计算开销大。基于深度学习的方法已成为建模物理系统动态的强大工具。监督学习(SL)[21, 23]训练参数化模型,利用时间反向传播直接优化控制,覆盖整个轨迹。例如,[21]提出了一个层次化的预测-校正方案,用于在长时间帧上控制复杂的非线性物理系统。一项最近的工作[23]设计了两个阶段,分别学习解算子并搜索最优控制。与这些方法不同,我们没有使用代理模型,而是以集成的方式学习状态轨迹和控制序列。强化学习(RL)[14, 47, 53]将控制信号视为动作,并学习策略做出序列决策。特别是在流体动力学领域[31],强化学习已应用于许多特定问题,包括阻力减少[53, 13]、共轭热传递[3, 16]和游泳[46, 64]。但它们隐式考虑了物理信息,并顺序地做出决策。相比之下,我们概括了整个轨迹,从而实现了全局优化,考虑了模型学习的物理信息。最近,PINNs 也被应用于 PDE 控制[43],但它们需要 PDE 动态的显式形式,而我们的方法是数据驱动的,能够处理不依赖显式 PDE 动态的更广泛的复杂物理系统控制问题。

A.3 Diffusion Models

扩散模型[18]在图像和文本生成[10, 45]、逆向设计[70, 66]、逆向问题[22]、物理模拟[6, 51]和决策制定[26, 1, 17]等应用中取得了显著进展。特别是,最近在机器人控制中的进展表明,扩散模型在行动规划上显著优于现有的强化学习方法[8, 71]。生成多样而一致的样本是一个挑战。为了实现多样性,结合各种模型的得分估计的方法[38, 2, 73, 12]已被证明是有效的。为了保持一致性,指导扩散技术[10, 19]被用来生成特定条件的样本。我们的方法通过平坦化联合分布来实现更好的控制,稍微扩展超出先验分布范围,从而优化控制

B 1D Burgers’ Equation Visualization

我们在图 6、图 7 和图 8 中展示了在三种设置下(FOPC、POFC 和 POPC)我们方法和基准方法的更多可视化结果。在每种设置下,我们展示了从测试数据集中随机选择的五个样本的结果。控制目标是使最终状态 u T u_T uT T = 10 T = 10 T=10 )尽可能接近目标状态。

在这里插入图片描述
图 6:在 FO-PC(完全观察,部分控制)设置下,1D Burgers 方程控制的可视化结果。我们的方法(DiffPhyCon)和基准方法在控制下的每个时间步 t = 0 , ⋯ , 10 t = 0, \cdots, 10 t=0,,10 的系统状态 u t u_t ut 曲线被绘制出来。x 轴是空间坐标,y 轴是系统状态的值。

C 附加信息:1D Burgers 方程的控制

C.1 数据生成

我们采用有限差分法(以下简称为求解器或真实求解器)生成 1D Burgers 方程的训练数据。具体而言,初值 u 0 ( x ) \mathbf u_0(x) u0(x) 和控制序列 w ( t , x ) \mathbf w(t,x) w(t,x) 均为随机生成,然后使用求解器数值计算得到状态 u ( t , x ) \mathbf u(t, x) u(t,x)

在数值模拟(使用真实求解器)中,模拟了 x = [ 0 , 1 ] x=[0,1] x=[0,1] t = [ 0 , 1 ] t=[0,1] t=[0,1] 范围内的区域。空间被离散为 128 个网格,时间被离散为 10,000 个时间步。然而,在数据集中,仅保存了 10 个时间点的数据。对于控制序列 w \mathbf w w,其刷新频率为 0.1 − 1 0.1^{-1} 0.11,即在 t ∈ [ 0.1 k , 0.1 ( k + 1 ) ] , k ∈ { 0 , . . . , 9 } t \in [0.1k, 0.1(k + 1)], k \in \{0, ..., 9\} t[0.1k,0.1(k+1)],k{0,...,9} 范围内, w ( t , x ) \mathbf w(t, x) w(t,x) 不随时间 t t t 改变。因此,每条轨迹的数据大小为 [ 11 , 128 ] [11, 128] [11,128](状态 u \mathbf u u)和 [ 10 , 128 ] [10, 128] [10,128](控制 w \mathbf w w)。

在所有设置中,初始值 u ( 0 , x ) \mathbf u(0,x) u(0,x) 是两个高斯函数的叠加 u ( 0 , x ) = ∑ i = 1 2 a i e − ( x − b i ) 2 2 σ i 2 , \mathbf u(0, x) = \sum_{i=1}^2 a_i e^{-\frac{(x-b_i)^2}{2\sigma_i^2}}, u(0,x)=i=12aie2σi2(xbi)2,其中 a i a_i ai b i b_i bi σ i \sigma_i σi 均从以下均匀分布中随机采样: a 1 ∼ U ( 0 , 2 ) , a_1 \sim U(0, 2), a1U(0,2), a 2 ∼ U ( − 2 , 0 ) , a_2 \sim U(-2, 0), a2U(2,0), b 1 ∼ U ( 0.2 , 0.4 ) , b_1 \sim U(0.2, 0.4), b1U(0.2,0.4), b 2 ∼ U ( 0.6 , 0.8 ) , b_2 \sim U(0.6, 0.8), b2U(0.6,0.8), σ 1 ∼ U ( 0.05 , 0.15 ) , \sigma_1 \sim U(0.05, 0.15), σ1U(0.05,0.15), σ 2 ∼ U ( 0.05 , 0.15 ) \sigma_2 \sim U(0.05, 0.15) σ2U(0.05,0.15)

类似地,控制序列 w ( t , x ) \mathbf w(t, x) w(t,x) 也是 8 个高斯函数的叠加

w ( t , x ) = ∑ i = 1 8 a i e − ( x − b 1 , i ) 2 2 σ 1 , i 2 e − ( t − b 2 , i ) 2 2 σ 2 , i 2 , (20) \mathbf w(t, x) = \sum_{i=1}^8 a_i e^{-\frac{(x-b_{1,i})^2}{2\sigma_{1,i}^2}} e^{-\frac{(t-b_{2,i})^2}{2\sigma_{2,i}^2}}, \tag{20} w(t,x)=i=18aie2σ1,i2(xb1,i)2e2σ2,i2(tb2,i)2,(20)

其中每个参数独立生成: b 1 , i ∼ U ( 0 , 1 ) , b_{1,i} \sim U(0, 1), b1,iU(0,1), b 2 , i ∼ U ( 0 , 1 ) , b_{2,i} \sim U(0, 1), b2,iU(0,1), σ 1 , i ∼ U ( 0.05 , 0.2 ) , \sigma_{1,i} \sim U(0.05, 0.2), σ1,iU(0.05,0.2), σ 2 , i ∼ U ( 0.05 , 0.2 ) , \sigma_{2,i} \sim U(0.05, 0.2), σ2,iU(0.05,0.2), a 1 ∼ U ( − 1.5 , 1.5 ) a_1 \sim U(-1.5, 1.5) a1U(1.5,1.5),对于 i ≥ 2 i \geq 2 i2 a i ∼ U ( − 1.5 , 1.5 ) a_i \sim U(-1.5, 1.5) aiU(1.5,1.5) 或 0(以相等概率生成)。

在基于方程 (16) 的数值求解中,给定 u ( 0 , x ) \mathbf u(0, x) u(0,x) w ( t , x ) \mathbf w(t, x) w(t,x) 后计算得到 u ( t , x ) , t ≠ 0 \mathbf u(t, x), t \neq 0 u(t,x),t=0。数据集生成的设置基于之前的研究 [23]。

我们为训练集生成了 90,000 条轨迹,为测试集生成了 50 条轨迹。每条轨迹占用 32KB 空间,数据集总大小为 2GB。

C.2 实验设置

在推理阶段,与控制序列 w ( t , x ) \mathbf w(t, x) w(t,x) 一起,我们的扩散模型生成状态 μ ( t , x ) \mu(t, x) μ(t,x)。某些模型在将控制 w ( t , x ) \mathbf w(t, x) w(t,x) 输入到相应的代理模型时,会生成代理状态 μ ( t , x ) \mu(t, x) μ(t,x)。但我们报告的评估指标 J actual \mathcal J_{\text{actual}} Jactual 始终通过将控制 w ( t , x ) \mathbf w(t, x) w(t,x) 输入真实数值求解器获得 u g.t. ( t , x ) \mathbf u_{\text{g.t.}}(t, x) ug.t.(t,x) 并按照公式 (17) 计算得到。

以下是三种不同实验设置的描述。

C.2.1 部分观测,完全控制(PO-FC)

在现实场景中,系统通常无法完全观测。一般来说,在系统的每个位置放置传感器是不可行的,因此模型从不完整数据中学习的能力至关重要。为了评估这一点,在此设置中我们隐藏了部分 u \mathbf u u 的数据,并测量模型控制的 J actual \mathcal J_{\text{actual}} Jactual

具体来说,训练数据集中 u ( t , x ) , x ∈ [ 1 4 , 3 4 ] \mathbf u(t, x), x \in [\frac{1}{4}, \frac{3}{4}] u(t,x),x[41,43] 被设置为零,测试时 u 0 ( x ) , x ∈ [ 1 4 , 3 4 ] \mathbf u_0(x), x \in [\frac{1}{4}, \frac{3}{4}] u0(x),x[41,43] 也被设置为零。在此部分观测设置中, Ω = [ 1 , 1 4 ] ∪ [ 3 4 , 1 ] \Omega = [1, \frac{1}{4}] \cup [\frac{3}{4}, 1] Ω=[1,41][43,1]。由于中间一半空间的数据从未被观测到,模型无法得知这些未观测状态对控制结果的影响。因此,控制未观测状态并非合理任务,它们被排除在评估指标之外。

这种设置特别具有挑战性,因为未观测状态带来的不确定性,以及中间区域的控制生成会隐式影响受控区域 x ∈ Ω x\in \Omega xΩ 内的 u \mathbf u u

C.2.2 完全观测,部分控制

在这一实际相关的场景中,系统的控制仅限于部分区域。在控制序列中,中央区域 x ∈ [ 1 4 , 3 4 ] x \in \left[\frac{1}{4}, \frac{3}{4}\right] x[41,43] 的值被强制设为零。然而,整个空间 Ω = [ 0 , 1 ] \Omega = [0, 1] Ω=[0,1] 的观测状态 J \mathcal J J 都会被用于评估。

在数据生成过程中,首先按照之前的方式生成控制序列 w \mathbf w w,然后将中央 1 2 \frac{1}{2} 21 区域的 w \mathbf w w 设置为零。为了补偿控制强度的下降,从而使状态 u \mathbf u u 的幅值大致与完全控制的设置相当,控制序列 w \mathbf w w 的幅值被加倍。在评估阶段,输出的控制序列同样会被处理,使中央区域 x ∈ [ 1 4 , 3 4 ] x \in \left[\frac{1}{4}, \frac{3}{4}\right] x[41,43] 的值为零。

需要注意的是,即使在控制能量没有任何限制的情况下,这种设置仍然很难实现完美的控制。因为模型需要学会如何通过间接方式影响中央区域的状态,从而达到全局控制的目标。

C.2.3 部分观测,部分控制

这一设置是前两种场景的结合。系统仅在 Ω = [ 0 , 1 4 ] ∪ [ 3 4 , 1 ] \Omega = \left[0, \frac{1}{4}\right] \cup \left[\frac{3}{4}, 1\right] Ω=[0,41][43,1] 区域内被观测、控制和评估。

需要注意的是,某些模型在生成输出时需要访问当前状态。如果模型与真实解算器交互,而不是使用替代模型(如仿真模型),则可能会导致结果异常优越,因为未观测区域的信息可能通过交互被泄露。

C.3 模型设计

由于模型 ϵ ϕ ≈ ∇ w log ⁡ p ( w ) \epsilon_\phi \approx \nabla_\mathbf w \log p(\mathbf w) ϵϕwlogp(w) ϵ θ ≈ ∇ u , w log ⁡ p ( u , w ) \epsilon_\theta \approx \nabla_{\mathbf u, \mathbf w} \log p(\mathbf u,\mathbf w) ϵθu,wlogp(u,w) 的训练方式基本相同,而后者正是 DiffPhyCon-lite 的核心,我们将首先介绍 DiffPhyCon-lite。

C.3.1 DiffPhyCon-lite

总体而言,DiffPhyCon-lite 遵循文献 [18] 的公式,并在正文中有所描述。状态 u \mathbf u u 和控制 w \mathbf w w 的数据被分别处理为大小为 ( N t , N x ) (N_t, N_x) (Nt,Nx) 的图像,其中 N t N_t Nt 表示时间步数, N x N_x Nx 表示空间网格数(128)。由于 u \mathbf u u w \mathbf w w 的时间步数不一致(分别为 11 和 10),数据会被填充为大小为 16 的统一格式。然后, u \mathbf u u w \mathbf w w 被堆叠为两个通道输入到 2D DDPM 模型中。

模型使用一个 2D UNet ϵ θ \epsilon_\theta ϵθ 预测噪声 ϵ \epsilon ϵ。其结构包括三部分:降采样编码器中心模块升采样解码器。降采样编码器包含四层,每层包括两个 ResNet 块、一个线性注意力块和一个降采样卷积块。中心模块由两个 ResNet 块和一个线性注意力块组成。升采样解码器的每层结构与降采样编码器类似,只是降采样块被替换为升采样卷积块。

实验发现,当学习条件概率分布 p ( w [ 0 , T − 1 ] , u [ 1 , T − 1 ] ∣ u 0 , u T ) p(\mathbf w_{[0,T-1]}, \mathbf u_{[1,T-1]} | \mathbf u_0, \mathbf u_T) p(w[0,T1],u[1,T1]u0,uT) 时,控制效果最佳。总结来说, ϵ θ \epsilon_\theta ϵθ 接收当前轨迹 u \mathbf u u、控制序列 w \mathbf w w、时间步数 k k k、初始状态 u 0 \mathbf u_0 u0 和终态 u T \mathbf u_T uT 作为输入,预测 u \mathbf u u w \mathbf w w 的噪声。需要注意的是, ϵ θ \epsilon_\theta ϵθ 并不直接预测 u 0 \mathbf u_0 u0 u T \mathbf u_T uT,它们只是作为条件输入,但模型输出仍会在相应位置保持一致的数据形状。不同设置的超参数详见表 4。

C.3.2 DiffPhyCon

DiffPhyCon 的实现仅是在推断阶段将 ϵ ϕ ( w ) \epsilon_\phi(\mathbf w) ϵϕ(w) 添加到 DiffPhyCon-lite 的输出 ϵ θ ( u , w ) \epsilon_\theta(\mathbf u, \mathbf w) ϵθ(u,w),如正文 3.2 节所述。其中, ϵ θ \epsilon_\theta ϵθ 是 DiffPhyCon-lite 的去噪网络输出,而 ϵ ϕ \epsilon_\phi ϵϕ 是一个新的去噪网络,用于根据数据集分布生成控制序列 w \mathbf w w。因此,这里只描述 ϵ ϕ \epsilon_\phi ϵϕ 的模型设计。

ϵ ϕ \epsilon_\phi ϵϕ 的输入包括控制序列 w \mathbf w w、时间步数 k k k 以及引导条件 u 0 \mathbf u_0 u0 u T \mathbf u_T uT ϵ ϕ ( w ) \epsilon_\phi(\mathbf w) ϵϕ(w) 的输出形状与 w \mathbf w w 一致,因此它可以被看作一个学习从 p ( w ) = ∫ p ( u , w ) d u p(\mathbf w) = \int p(\mathbf u, \mathbf w) \text d\mathbf u p(w)=p(u,w)du 中采样的网络。在 u \mathbf u u 所在的位置, ϵ ϕ ( w ) \epsilon_\phi(\mathbf w) ϵϕ(w) 的输出填充为零。模型的超参数同样列于表 4。

C.4 训练与评估

训练

在训练阶段,无噪声的 u 0 \mathbf u_0 u0 u T \mathbf u_T uT 作为输入被传递给模型,并且在这些位置的模型输出不会被纳入损失计算。在部分观测的场景中,如附录 C.2 所述,模型在训练和测试中都无法访问未观测区域的数据。这些位置的输入被填充为零,并从训练损失中排除。因此,模型仅学习观测状态与控制序列之间的关联。

在部分控制的设置下,我们使用 x ∈ [ 1 4 , 3 4 ] x \in \left[\frac{1}{4}, \frac{3}{4}\right] x[41,43] 的控制值为零的数据集对 DiffPhyCon 进行训练。这种方式使得模型自然学会在“不可控”区域输出零。

去噪 UNet 的训练采用均方误差(MSE)损失,其他训练超参数详见表 4。

在这里插入图片描述

推理

在推理阶段, u 0 \mathbf u_0 u0 u T \mathbf u_T uT 被设为目标状态 u 0 \mathbf u_0 u0 u T \mathbf u_T uT,从而使 DDPM 生成的样本既满足物理约束,又以目标 u T \mathbf u_T uT 或条件 u 0 \mathbf u_0 u0 为条件。在部分观测场景中,从测试集中抽取的 u 0 \mathbf u_0 u0 u T \mathbf u_T uT 会在未观测的 x ∈ [ 1 4 , 3 4 ] x \in \left[\frac{1}{4}, \frac{3}{4}\right] x[41,43] 区域被填充为零,与训练 UNet 使用的数据一致。

在部分控制场景下,推理过程中我们将去噪网络的输出 ϵ θ ( u , w ) \epsilon_\theta(\mathbf u,\mathbf w) ϵθ(u,w) 替换为 ϵ θ ( u , w ) + ( γ − 1 ) ϵ ϕ ( w ) \epsilon_\theta(\mathbf u,\mathbf w) + (\gamma - 1)\epsilon_\phi(\mathbf w) ϵθ(u,w)+(γ1)ϵϕ(w)。需要注意的是, ϵ θ ( u , w ) \epsilon_\theta(\mathbf u, \mathbf w) ϵθ(u,w) 同时去噪 u \mathbf u u w \mathbf w w,而 ϵ ϕ ( w ) \epsilon_\phi(\mathbf w) ϵϕ(w) 仅去噪 w \mathbf w w。实验表明,为 w \mathbf w w 网络的输出添加一个调度机制是有益的。表 1 中的结果采用了基于逆 Sigmoid 的调度,公式如下:

ϵ k = ϵ θ ( u k , w k , k , u 0 , u T ) + ( γ − 1 ) β K − k ϵ ϕ ( w k , k , u 0 , u T ) , (21) \epsilon_k = \epsilon_\theta(\mathbf u_k, \mathbf w_k, k, \mathbf u_0, \mathbf u_T) + (\gamma - 1)\beta_{K-k}\epsilon_\phi(\mathbf w_k, k, \mathbf u_0, \mathbf u_T), \tag{21} ϵk=ϵθ(uk,wk,k,u0,uT)+(γ1)βKkϵϕ(wk,k,u0,uT),(21)

其中, β \beta β 是噪声调度参数,参见文献 [18]。DiffPhyCon-lite 的推理过程仅需设置 γ = 1 \gamma = 1 γ=1,忽略模型 ϵ ϕ ( w ) \epsilon_\phi(\mathbf w) ϵϕ(w) 的影响。

当需要对控制能量进行调节时,学习条件模型并非一种自然的选择。因此,我们引入外部引导项 J = ∫ w ( x , t ) d x d t \mathcal J = \int \mathbf w(x, t) \,\text dx \text dt J=w(x,t)dxdt 来生成能量受限的控制序列,其结果展示在图 3 中。引导项的梯度被计算后加到公式 (21) 中的 ϵ k \epsilon_k ϵk。需要注意的是,梯度计算基于去噪后的预测样本 w ^ k \hat{\mathbf w}_k w^k u ^ k \hat{\mathbf u}_k u^k,因为它们的噪声较少且引导偏差较小。 u k \mathbf u_k uk w k \mathbf w_k wk 按如下公式计算:

x 0 ≈ x ^ k = 1 α ˉ k x k − 1 − α ˉ k α ˉ k ϵ k , x_0 \approx \hat {\mathbf x}_k = \frac{1}{\sqrt{\bar \alpha_k}} \mathbf{x}_k - \sqrt{\frac{1 - \bar \alpha_k}{\bar \alpha_k}} \epsilon_k, x0x^k=αˉk 1xkαˉk1αˉk ϵk,

其中 x \mathbf x x 表示 u \mathbf u u w \mathbf w w,详见文献 [18]。实验中我们对外部引导项采用余弦调度,因此最终的预测噪声为 ϵ k + λ α k ∇ J \epsilon_k + \lambda\alpha_k\nabla \mathcal J ϵk+λαkJ,其中 β k \beta_k βk 是余弦调度的噪声参数,参见文献 [18]。

评估

生成的状态轨迹 u \mathbf u u 和控制序列 w \mathbf w w 会被输入到真实解算器中,以模拟基于生成的 w \mathbf w w 和测试数据集中初始条件 u 0 \mathbf u_0 u0 的最终状态 u \mathbf u u。解算器与附录 C.1 数据生成过程中使用的相同。最终,通过公式 (17) 计算实际的 J actual \mathcal J_{\text{actual}} Jactual

在部分观测场景下,仅在观测区域计算均方误差(MSE)。在部分控制场景中,生成的控制序列 w \mathbf w w 在输入解算器前会首先将不可控区域设置为零。

相关文章:

  • 【Django】性能优化-普通版
  • C++ 网络编程(11)服务器逻辑层设计和消息完善
  • 7.7 Extracting and saving responses
  • 【医学目标检测】LN-DETR:一种基于多尺度特征融合的肺结节检测高效Transformer架构
  • Ceph分布式存储方案
  • Barcode解码 一维码、二维码识别 物流单号识别
  • ss928v100模型的导出、量化和转换
  • Kotlin的MutableList和ArrayList区别
  • Kotlin基础语法三
  • IntelliJ IDEA代码提示忽略大小写设置详解
  • 容器化包允许应用程序使用 Linux 容器
  • PyTorch:让深度学习像搭积木一样简单!!!
  • 中兴B860AV1.1_MSO9280_降级后开ADB-免刷机破解教程(非刷机)
  • IntelliJ IDEA 豆沙绿护眼色设置
  • ADB(Android Debug Bridge)Android官方调试工具及常用命令
  • 广告推荐系统中模型训练中模型的结构信息、Dense数据、Sparse数据
  • SQL ConcurrencyControl(并发控制)
  • 【机器学习-线性回归-7】中心极限定理在机器学习线性回归中的重要性
  • 从认识AI开始-----生成对抗网络(GAN):通过博弈机制,引导生成
  • 机器学习与深度学习19-线性代数02
  • 大朗仿做网站/谷歌浏览器2021最新版
  • 北京网站建设要多少钱/好的搜索引擎推荐
  • 网站测试域名301怎么做/企业所得税优惠政策
  • c mvc制作网站开发/口碑营销的产品有哪些
  • flash网站建设技术.../新手怎样做网络推广
  • 做网站赚不到钱了/蓝牙耳机网络营销推广方案