VC-PINN Variable coefficient physics-informed neural network for forward and inverse problems

VC-PINN: Variable coefficient physics-informed neural network for forward and inverse problems of PDEs with variable coefficient

变系数问题

考虑实空间中一类具有时变系数的演化方程，如下：
$$u_t=\mathbf{N}[u]\cdot \mathbf{C}[t{]}^T,\mathbf{x}\in \Omega ,t\in [T_0,T_1]$$
其中 $u=u(x,t)$ 表示方程的实值解，$\Omega$ 是 ${\mathbb{R}}^N$ 的子集，$N$ 维向量 $\mathbf{x}$ 记为 $\mathbf{x}=(x_1,x_2,\dots ,x_N)$ ，因此，此方程实际为一个$N+1$ 维的演化方程。$\mathbf{N}[\cdot ]$ 表示算子向量，即 $\mathbf{N}[u]=(N_1[u],N_2[u],\dots )$ 其中每个分量 $N_i$ 为一个算子。$\mathbf{C}[t]=(c_1(t),c_2(t),\dots )$ 是一个系数向量，其分量 $c_i(t)$ 是时间变量 $t$ 的函数，并且 $\mathbf{C}[t]$ 的维度与 $\mathbf{N}[u]$ 相同。此外，$\mathbf{N}[u]$ 和 $\mathbf{C}[t]$ 相应分量的乘积表示具有时变系数的算子 $(c_i[t]N_i[u])$ 。特别地，这里考虑的可变系数只是时间变量 $t$ 的函数。变量系数也与空间变量 $x$ 相关的情况将在后文中讨论。

对于连续意义上的前向问题，变量系数的表达式是完全已知的并写入方程中。在处理此类问题时，像常系数问题一样使用标准的PINN方法就足够了。但在工程应用中，充分了解变系数表达式的要求十分苛刻。因此，后续基于正向和逆向问题的讨论都是在离散意义上进行的。具体来说，文中以离散意义上的变量系数是否已知作为区分正问题和逆问题的基础。

变系数条件下偏微分方程的正、反问题的形式化定义

鉴于变系数方程与常系数方程在正、反问题描述上的差异，有必要给出变系数版本下偏微分方程的正、反问题的形式化定义。如下：

• 在前向问题中，可变系数 $c_i(t)$ 是已知的。具体地，$[T_0,T_1]$ 上有限个离散点处的系数 $c_i(t)$ 的值是已知信息。这些系数值对应于实际应用中一些可观测的物理量（具有不同程度的噪声）。因此，可变系数的前向问题可被形式化定义为：使用 $u$ 的初始边界值条件（与常数系数一致）和上述离散系数值来求解该区域上的 $u$。
• 在反问题中，待确定的系数不再是几个固定的常数，而是与时间变量 $t$ 有关的一组函数。与常系数问题一致，区域内有限离散点处的 $u$ 值是已知信息。此外，由于反问题的多解性质，需要提供变量系数 $c_i(t)$ 的边界值$(c_i(T_0),c_i(T_1))$。该边界信息对应于实验中观测到的初始值和终值。综上所述，逆问题可被形式化定义为：利用上述两个已知信息，得到时域 $[T_0，T_1]$ 上完整的系数变化，即任意时刻 $c_i(t)$ 的离散值。

网络结构

在变系数问题中，不仅 $u$ 需要用神经网络来表示，而且不同数量的自变量的变系数 $c_i(t)$ 也需要用新的网络来近似。不失一般性，在方法描述中，假设方程仅涉及单变量系数，即 $\mathbf{C}[t]=c_1(t)$ ，对应的算子向量 $\mathbf{N}[u]$ 也简写为 $\mathcal{N}[u]$。该方法同样适用于多变量系数的情况。这种简化只是为了更清晰的描述。下文数值实验中也给出了多变量系数的例子。

由于ResNet在PDE求解中已得到广泛使用，同时，在变系数问题上，ResNet可以更好地统一网络中的线性和非线性，以适应不同的变系数，因此，下文作者选择ResNet来构建整个网络。

首先考虑深度为 $D$的前馈神经网络（FNN）。第 $0$ 层和第 $D$ 层分别称为输入层和输出层，自然有 $D-1$ 个隐藏层。在引入快捷连接之前必须考虑的一个特殊要求是要连接的两个向量需要具有相同的维度。常见的做法是使用线性投影来使得尺寸满足上述条件。本着尽可能不引入新的网络参数的原则，可以假设每个隐藏层的节点数为 $N_d$ 。那么，一个具有 $N_B$个残差块且每个残差块包含 $N_h$ 个隐藏层的 ResNet 结构图如下：

$X^{[d]}$ 表示第 $d$ 层节点的状态向量，$R^{[i]}$ 既是第 $i$ 个残差块的输出，也是第 $(i+1)$ 个残差块的输入。第 i 个残差块由网络的第 $[(i-1)N_h+1]$ 层到第 $(i{N}_h)$ 层组成。所以 $R^{[i]}$等于第 $(i{N}_h+1)$ 层节点的状态向量$X^{[i{N}_h+1]}$。为了更清晰地展示网络的结构，也为了统一表达ResNet和普通的FNN，更数学化的表达方式如下：
$$\begin{array}{rl}{X}_L^{[0]}& =W^{[0]}{X}^{[0]}+b^{[0]},\\ X^{[1]}={\mathcal{R}}^{[0]}& =\mathcal{K}{X}_L^{[0]}+(1-\mathcal{K})\mathcal{F}(X_L^{[0]}),\\ X^{[i{N}_h+1]}={\mathcal{R}}^{[i]}& ={\mathcal{L}}_i({\mathcal{R}}^{[i-1]})+\mathcal{K}{\mathcal{R}}^{[i-1]},i=1,2,\dots ,N_B,\\ X^{[D]}& =W^{[D-1]}{X}^{[D-1]}+b^{[D-1]},\\ & =W^{[D-1]}{\mathcal{R}}^{[N_B]}+b^{[D-1]},\end{array}$$
其中 $X_L^{[0]}$ 为中间变量，系数 K ∈ { 0 , 1 } \mathcal K \in \{0, 1\} K∈{0,1} 主要控制是否包含捷径连接。具体来说，当 $\mathcal{K}=1$时为ResNet结构，$\mathcal{K}=0$ 时为普通FNN结构。非线性映射${\mathcal{L}}_i$ 是第 $i$ 个残差块的输入和输出之间的非线性部分，定义如下：
$$\begin{array}{l}{\mathcal{L}}_i\triangleq {\mathcal{T}}_{i{N}_h}\circ {\mathcal{T}}_{i{N}_h-1}\circ \cdots \circ {\mathcal{T}}_{(i-1)N_h+1},i=1,2,\dots ,N_B,\\ {\mathcal{T}}_i(X)\triangleq \mathcal{F}(W^{[i]}X+b^{[i]}),i=1,2,\dots ,N_B{N}_h,\end{array}$$
其中， $W^{[i]}\in {\mathbb{R}}^{N_i\times N_{i+1}}$ 和 $b^{[i]}\in {\mathbb{R}}^{N_{i+1}}$ 分别表示第 $i$ 层与第 $(i+1)$ 层网络之间的权重矩阵和偏置向量，其中 $N_i$ 为第 $i$ 层网络中的节点，且 $N_i=N_d,i=1,2,...,D-1$。

值得注意的是，作者这里使用的是预激活的ResNet。

为了在后续的介绍中区分不同的网络，具体的ResNet被表示为结构参数列表形式 N N { D , N d , N B , N h } K = 1 NN\{D,N_d,N_B,N_h\}_{\mathcal K=1} NN{D,Nd​,NB​,Nh​}K=1​，其中 $\mathcal{K}$ 的定义与先前定义相同。特别地，当 $\mathcal{K}=0$ 时，参数列表 N N { D , N d , N B , N h } K = 0 NN\{D,N_d,N_B,N_h\}_{\mathcal K=0} NN{D,Nd​,NB​,Nh​}K=0​ 表示普通的FNN，则网络结构仅由 $D$ 和 $N_h$ 决定，因此简写为 N N { D , N h } K = 0 NN\{D, Nh\}_{\mathcal K=0} NN{D,Nh}K=0​。

本文需要构建两个网络来分别逼近解 $u$ 和变量系数 $c_1(t)$。两个网络都将采用上述ResNet结构，称为主干网络 N N u { D u , N d u , N B u , N h u } NN_u\{D_u, N^u_d, N^u_B, N^u_h\} NNu​{Du​,Ndu​,NBu​,Nhu​} 和分支网络 N N c { D c , N d c , N B c , N h c } NN_c \{D_c , N^c_d, N^c_B, N^c_h \} NNc​{Dc​,Ndc​,NBc​,Nhc​}

前向问题

考虑涉及单个时变系数的偏微分方程的初始边值问题（狄利克雷边界条件）如下：
$$\begin{array}{rl}& u_t=c_1(t)\mathcal{N}[u],x\in \Omega ,t\in [T_0,T_1],\\ & u(x,T_0)=g_0(x),x\in \Omega ,\\ & u(x,t)=g_{\Gamma }(x,t),x\in \partial \Omega ,t\in [T_0,T_1],\end{array}$$
其中 $\partial \Omega$表示空间域 $\Omega$ 的边界，第一个方程是PDE的特例，后两个方程分别对应初值条件和狄利克雷边界条件。当 $c_1(t)$ 已知（离散意义上已知）时，采用神经网络方法求解初始边值问题，关键是构造优化问题。

主干网络 N N u { D u , N d u , N B u , N h u } NN_u\{D_u, N^u_d, N^u_B, N^u_h\} NNu​{Du​,Ndu​,NBu​,Nhu​} 用于逼近初边界问题的实值解 $u(x,t)$ ，记为 $\tilde{u}(x,t;{\theta }_u)$，用分支网络 N N c { D c , N d c , N B c , N h c } NN_c \{D_c , N^c_d, N^c_B, N^c_h \} NNc​{Dc​,Ndc​,NBc​,Nhc​} 表示的用于逼近变系数函数 $c_1(t)$ ，记为 $\tilde{c}(t;{\theta }_c)$。 ${\theta }_u\in {\Theta }_u$ 和 ${\theta }_c\in {\Theta }_c$ 分别是 $N{N}_u$ 和 $N{N}_c$ 两个网络的参数空间（权重和偏置空间）。为了引入损失函数，将方程在点 $(\tilde{x},\tilde{t})$处的残差定义如下：
$$f(\tilde{\mathbf{x}},\tilde{t};u,c):={\mathcal{N}}_0[u,c]{|}_{\mathbf{x}=\tilde{\mathbf{x}},t=\tilde{t}},{\mathcal{N}}_0[u,c]={\partial }_t[u]-c\mathcal{N}[u].$$
残差由问题的第一个方程导出， ${\mathcal{N}}_0[\cdot ,\cdot ]$ 是由 $\mathcal{N}[\cdot ]$ 和 ${\partial }_t$ 组成的新算子。然而，$f(\sim ,\sim ;u,c)$ 被视为函数的函数，用于将函数空间中的点 $(u,c)$ 映射到域 $\Omega \times [T_0,T_1]$ 上的函数，其中 $u=u(x,t)$ 和 $c=c(t)$ 被解释为函数类型参数。因此，上式定义的残差衡量了在给定函数 $u(x,t)$ 和系数 $c(t)$ 的情况下方程在 $(\tilde{x},\tilde{t})$ 点满足的程度。特别地，若 $u_0$ 是变系数 $c_1(t)$下初始边值问题的解，则显然 $f(\tilde{x},\tilde{t};u_0,c_1)=0,\forall \tilde{x}\in \Omega ,\tilde{t}\in [T_0,T_1]$。参数空间 θ = { θ u , θ c } \theta = \{\theta_u, \theta_c \} θ={θu​,θc​} 中的每组参数 Θ = { Θ u , Θ c } \Theta_ = \{\Theta_u, \Theta_c \} Θ=​{Θu​,Θc​} 定义函数 $ \tilde{u}(x, t; \theta_u)$ 和可变系数 $\tilde{c}(t;{\theta }_c)$。然而，从参数空间 $\Theta$ 中找到合适的参数 ${\theta }^{\star }={\theta }_u^{\star },{\theta }_c^{\star }$ 使得残差 $f(\tilde{x},\tilde{t};{\tilde{u}}^{\star },{\tilde{c}}^{\star })$ 足够接近于零域 $\Omega \times [T_0,T_1]$ 是优化的目标。同时，若 ${\tilde{u}}^{\star }$ 满足初始边界条件那么它足够接近初始边值问题的真实解。

为了更好的衡量真实解 $u_0(x,t)$ 与近似解 ${\tilde{u}}^{\star }$ 之间的差距，损失函数可以设置如下：
$$Loss(\theta )=Los{s}_I(\theta )+Los{s}_b(\theta )+Los{s}_f(\theta )+Los{s}_c(\theta ),$$
其中：
$$\begin{array}{rl}& Los{s}_I(\theta )=\frac{1}{n_I}\sum _{i=1}^{n_I}|\tilde{u}(x_I^i,T_0;{\theta }_u)-g_0({\mathbf{x}}_I^i){|}^2,\\ & Los{s}_b(\theta )=\frac{1}{n_b}\sum _{i=1}^{n_b}|\tilde{u}(x_b^i,t_b^i;{\theta }_u)-g_{\Gamma }(x_b^i,t_b^i){|}^2,\\ & Los{s}_f(\theta )=\frac{1}{n_f}\sum _{i=1}^{n_f}|f(x_f^i,t_f^i;\tilde{u}(x,t;{\theta }_u),\tilde{c}(t;{\theta }_c)){|}^2\\ & Los{s}_c(\theta )=\frac{1}{n_c}\sum _{i=1}^{n_c}|\tilde{c}(t_c^i;{\theta }_c)-c_c^i{|}^2.\end{array}$$
上式中，$Los{s}_I$和$Los{s}_b$分别为初值约束和边界约束，$Los{s}_f$为物理约束，$Los{s}_c$为变系数问题中的唯一系数约束。

反问题

反问题方程的真实解 $u_0$ 在离散意义上是已知的，可变系数 $c_1$ 成为要求解的目标。因此，问题的表述也变为：
$$\begin{array}{rl}& u_t=c_1(t)\mathcal{N}[u],x\in \Omega ,t\in [T_0,T_1],\\ & u(x,t)=u_r(x,t),(x,t)\in D,\\ & c_1(T_0)=C_0,c_1(T_1)=C_1,\end{array}$$
其中第一行是原方程，第二行表示 $u$ 的离散观测点，$u_r$和 $\mathcal{D}$ 分别是实解和观测点的坐标集。第三行表示变系数的两端条件。在简单的问题中，这个条件可以放宽到单个端点，甚至不需要，而在更具挑战性的问题中，需要两个端点处的更高导数的信息，甚至需要内点信息。这里还给出了高阶导数的条件：
$${\frac{{\partial }^k{c}_1}{\partial t^k}|}_{t=T_0}=C_0^{(k)},{\frac{{\partial }^k{c}_1}{\partial t^k}|}_{t=T_1}=C_1^{(k)},k=1,2,\dots ,$$
其中 $C_0^{(k)}$ 和 $C_1^{(k)}$ 是两端对应的高阶导数值。在VC-PINN框架下，处理变系数的逆问题与其正问题具有一定的统一性，变化几乎只发生在损失函数的构成上。具体区别如下：
$$\begin{array}{rl}& Los{s}_I(\theta )+Los{s}_b(\theta )\to Los{s}_s(\theta ),\\ & Los{s}_s(\theta )=\frac{1}{n_s}\sum _{i=1}^{n_s}|\tilde{u}(x_s^i,t_s^i;{\theta }_u)-u_s^i{|}^2,\end{array}$$
其中，真实解 $u_0$ 在全时空区域的采样点称为$s$型点，它们构成集合 { x s i , t s i , u s i } i = 1 n s \{x^i _s, t^i _s, u^i _s \}^{n_s} _{i=1} {xsi​,tsi​,usi​}i=1ns​​ ，$u_s^i$ 为点 $(x_s^i,t_s^i)$ 处的真实解 $u_0$。 $s$ 型点代表已知的解信息（实际中视为可观测量）。于是，损失函数可表达如下：
$$Loss(\theta )=Los{s}_s(\theta )+Los{s}_f(\theta )+Los{s}_c(\theta ).$$

实验结果

作者在正向问题中对变系数Sine-Gordon方程和广义变系数Kadomtsev-Petviashvili方程进行了测试，在反问题中对广义变系数Kadomtsev-Petviashvili方程以及变量系数 Korteweg-de Vries 方程和变量系数 Sawada-Kotera 方程进行了测试。

正问题：变系数Sine-Gordon方程

给出具有可变系数的 SG 方程 (vSG)：
$$u_{xt}+h(t)\sin (u)=0$$
其中 $h(t)$ 是一个解析函数，表示方程的系数如何随时间变化。其单孤子解为：
$$u(x,t)=4\text{arctan}\left(e^{k_{1}x-{\omega }_1(t)}\right)$$
其中，${\omega }_1(t)=\int \frac{h(t)}{k_1}dt$ ，因此选择不同的 $h(t)$ 会产生许多具有丰富动态行为的解，这在常系数问题中是找不到的。接下来，作者讨论了三种形式（线性系数、二次多项式、三角函数）的系数函数，并使用所提出的 VC-PINN 方法获得相应初始边值问题的数据驱动解。

• 线性系数：假设 $h(t)=t$，并取不定积分中积分常数为0，则该问题的精确解如下
  $$u_1^{(vSG)}=4\arctan \left(e^{-\frac{t^2}{2k_1}+k_{1}x}\right)$$
  在自由参数 $k_1=\pm 1$ 的情况下，VC-PINN 方法找到的两个数据驱动的解决方案以及伴随的误差结果如下图所示：

  

• 二次多项式：假设系数函数为 $h(t)=t^2$，则方程的精确解为
  $$u_2^{(vSG)}=4\arctan \left(e^{-\frac{t^3}{3k_1}+k_{1}x}\right)$$
  

• 三角函数：当系数函数为具有周期性质的余弦函数，即 $h(t)=3\cos (2t)$ 时，对应的精确解为
  $$u_3^{(vSG)}=4\arctan \left(e^{-\frac{3\sin (2t)}{2k_1}+k_{1}x}\right)$$
  

反问题：广义变系数Kadomtsev-Petviashvili方程

变系数Korteweg-de Vries方程（vKdV）具体形式如下：
$$u_t+f(t)u{u}_x+g(t)u_{xxx}=0$$
其中 $f(t)$ 和 $g(t)$ 是任意解析函数。假设方程中的可变系数满足如下约束条件
$$g(t)=cf(t)$$
其中 $c$ 是任意常数。那么在这个约束下，其精确解具有以下形式：
$$u(x,t)=3c{\alpha }^2{\mathrm{sech}}^2\left[\frac{1}{2}\alpha (x-c{\alpha }^2\int f(t)dt)\right]$$
其中 $\alpha$ 是自由参数。显然，一旦确定了变系数 $f(t)$、参数 $c$ 和 $\alpha$，解析解就完全确定了。令 $c=\alpha =1$，在此参数设置下，作者讨论了3种不同形式的 $f(t)$ 来测试所提出的方法在反问题上的性能。

• 线性系数：当系数 $f(t)$ 是线性时（即 $f(t)=t$ ），精确解如下
  $$u_1^{(vKdV)}=3{\mathrm{sech}}^2\left[\frac{1}{4}(t^2-2x)\right]$$
  

• 三次多项式：假设变系数$f(t)$是三次多项式，即$f(t)=t^3$，则方程的精确解
  $$u_2^{(vKdV)}=3{\mathrm{sech}}^2\left[\frac{1}{8}(t^4-4x)\right]$$
  

• 三角函数：当系数函数为余弦函数，即$f(t)=cos(t)$时，对应vKdV方程的精确解为
  $$u_3^{(vKdV)}=3{\mathrm{sech}}^2\left[\frac{1}{2}(x-\sin (t))\right]$$
  

ResNet 对比

上图为vKdV方程线性系数和非线性系数下的实验结果， (a) 和 (b) 为不使用ResNet下不同网络层数带来的最终误差， © 和 (d) 为使用ResNet下不同网络层数带来的最终误差

可以看到，使用ResNet带来了更好的表现。

上表为对线性与非线性的对比，可以看出，ResNet在二者上都有所提升。

变系数凸性与学习的关系

作者在实验中发现特定情况下可变系数也会使神经网络的学习变得极其困难。

为了分析神经网络学习失败的原因，作者考虑观察训练过程中变量系数的预测结果的变化。下图显示了vKdV方程在五个不同变量系数下的学习过程（示例中为成功和失败）。

上图中五个例子的学习过程有一个共同的规律，即神经网络首先学习变系数的两个端点，然后逐渐学习中间区域。添加到损失函数中的可变系数边界的软约束可以合理地解释这种学习行为。此外，二次和四次多项式下的学习失败，并且神经网络在预学习（学习损失函数中系数的边界约束）后停滞不前。大梯度似乎并不是学习失败的合理解释，因为使用三次多项式系数进行学习是完全成功的。这似乎意味着神经网络的学习在面对强凸目标时受到阻碍。

为了进一步探究神经网络学习失败的原因，作者在前面的二次和四次多项式系数的例子中做出了不同的调整策略。具体如下：（1）减小一侧的时间间隔（从[−4, 4]到[−4, 2]）； (2) 减小两侧的时间间隔（从[−4, 4]到[−2, 2]）； (3)和(5)Lossc添加中间点信息((0, 0))；(4)和(6)添加系数边界的一阶导数信息。

上图中的(1)和(2)都没有改变曲线，只是缩短了时间范围，曲率最大的 $t=0$ 仍然包含在区间内。但这样的调整使得系数能够成功学习，这说明影响学习的是整个区间的凸性，而不是某一点的凸性。整个区间上凸性的积累使得预学习后的系数曲线与真实的系数曲线相差甚远，可能导致学习失败。

上图中的(3)-(6)涉及的两种调整策略提供了可变系数信息，使得预学习曲线能够更接近真实系数曲线，从而避免凸性积累带来的问题。两种策略提供的信息量很小。例如，中间点的策略只添加一个信息点就可以从失败中学习到成功。这表明预学习曲线和真实系数曲线之间的差距很重要，似乎对学习的成败起着决定性的作用。所有调整策略的二次和四次多项式系数学习都是成功的。 （虽然上图(4)的学习还是有点缺陷）。

总结

本文针对变系数问题提出了VC-PINN结构，并给出了离散形式下的正反问题定义。通过丰富的实验证实了VC-PINN在可变系数问题上的出色表现。同时对ResNet结构与传统FNN进行对比，展示了其在解决梯度消失问题以及统一线性与非线性上的必要性。文章中作者还针对实验中出现的训练失败现象进行了更加深入的研究，并根据实验结果推测系数函数在区间上的累计凸性会影响训练结果的准确性。

这篇文章的实验部分也太过详实了，令人叹为观止。

相关链接：

• 原文：VC-PINN: Variable coefficient physics-informed neural network for forward and inverse problems of PDEs with variable coefficient - ScienceDirect