PINN：把物理方程塞进神经网络当硬约束

PINN：把物理方程塞进神经网络当硬约束

拿不到完整的初始/边界样本，却清楚地知道背后的偏微分方程——那就把方程本身当作监督信号，塞进损失函数里。

一、背景与痛点

有些工程现场很尴尬：核反应堆包壳的温度场、油藏的压力分布、翼型周围的不可压缩流场，我们往往拿不到完整的初始/边界样本，但描述它们的偏微分方程（PDE）明明白白挂在那里。传统做法是上网格：有限差分、有限元、谱方法。这套体系成熟可靠，可一旦维数升上来，网格点数就指数爆炸（维数灾难）；碰到反问题（已知部分观测、反推物理参数）时，还要额外做一轮繁琐的数据同化。

物理信息神经网络（Physics-Informed Neural Network, PINN）给出一条完全不同的路：用一个标量神经网络 $u_{\theta }(x,t)$ 直接承担整个 PDE 的解算任务，用自动微分把 PDE 的残差当成损失项，和稀疏观测数据一起优化。它的核心诉求只有一句话：让网络在拟合数据的同时，还「知道」自己必须满足物理。这套范式的标志性论文，是 Raissi、Yazdani 与 Karniadakis 于 2019 年发表在《Journal of Computational Physics》上的工作，它第一次把「PDE 残差可微 + 多任务加权 + 反问题协同」这几件事一次性打通。

二、核心原理

自动微分：精确计算残差

传统有限差分离散二阶导时，截断误差是 $\mathcal{O}(\Delta x^2)$；PINN 不需要网格，却仍能拿到任意阶导数，因为它靠自动微分（autograd）在计算图上反推。给定标量输出 $u_{\theta }(x,t)$，其对 $x$ 的一阶、二阶导数可以嵌套求导得到：

$$u_x=\frac{\partial u_{\theta }}{\partial x},u_{xx}=\frac{\partial }{\partial x}\left(\frac{\partial u_{\theta }}{\partial x}\right)$$

对现代框架，只要把输入标记为可微、并在第一次求导时保留计算图，就能顺理成章地拿到高阶导数。这个机制让我们能「零成本」写出任意线性/非线性 PDE 的残差，例如：

$${\mathcal{R}}_{\theta }(x,t):=u_t+\mathcal{N}[u_{\theta };x,t]$$

其中算子 $\mathcal{N}$ 对经典方程各不相同：Burgers 是 $u{u}_x-\nu u_{xx}$，Poisson 是 $u_{xx}+u_{yy}$，Allen–Cahn 是 $u^3-u$。

多任务损失：物理与数据的统一框架

PINN 的训练目标不是逼近一堆样本对，而是把「物理约束」写进损失。Raissi 给出的多任务损失形式是：

$$\mathcal{L}(\theta )={\lambda }_r{\mathcal{L}}_r+{\lambda }_b{\mathcal{L}}_b+{\lambda }_0{\mathcal{L}}_0+{\lambda }_d{\mathcal{L}}_d$$

其中 ${\mathcal{L}}_r$ 是配置点上的 PDE 残差均方，${\mathcal{L}}_b$、${\mathcal{L}}_0$ 分别是边界条件与初始条件的均方误差，${\mathcal{L}}_d$ 是真实测量数据项（仅在反问题或数据驱动场景下出现）：

Lr=1Nr∑i=1Nr∣Rθ(xri,tri)∣2\\mathcal{L}_r = \\frac{1}{N_r}\\sum_{i=1}^{N_r}\\big|\\mathcal{R}_\\theta(x_r^i, t_r^i)\\big|^2Lr​=Nr​1​i=1∑Nr​​​Rθ​(xri​,tri​)​2

权重 $\lambda$ 是 PINN 工程的「暗物质」：太小则物理约束失效，太大则数据约束失声。后续工作用 NTK 或 learning-rate annealing 把 $\lambda$ 从超参数升格为可学习量。因为各项都是同一个 $\theta$ 的可微函数，它们能被同一次反向传播一次性吃掉——这正是自动微分与多任务优化的协同红利。

关于「软约束」与「硬约束」：把边界点加入训练集算 MSE 是软约束，实现简单，但有两个隐忧——边界外推没有保证，且当 ${\lambda }_b$ 偏小时残差项可能「绕过」边界，学出一个物理上错误却恰好拟合边界的野解。硬约束则把边界显式编码进网络结构（例如用满足 Dirichlet 条件的辅助函数把解空间限制在正确的子流形上），代价是需要为每个方程手写引导项，难以推广。工程上常走软约束，非齐次边界时推荐混合编码。下面代码为了直观看清残差如何反传，走的是以 PDE 残差为核心的软约束路径。

三、代码实战

为了不依赖深度学习框架、也方便逐行看清「残差 → 梯度 → 更新」这一链条，下面用纯 numpy 手写一个 tanh 小网络与 Adam 优化器，用中心差分拿到 $u_x,u_t,u_{xx}$，把 Burgers 方程的残差当损失来训练。目标方程：

$$u_t+u{u}_x=\nu u_{xx},\nu =0.01/\pi ,u(x,0)=-\sin (\pi x)$$

import numpy as np
np.random.seed(0)

nu = 0.01 / np.pi # Burgers 粘性系数
H = 0.05 # 有限差分步长

# ———- 极简 tanh-MLP：2 -> 32 -> 32 -> 1 ———-
layers = [2, 32, 32, 1]
W = [np.random.randn(layers[i], layers[i+1]) * np.sqrt(2/(layers[i]+layers[i+1]))
for i in range(len(layers)–1)]
B = [np.zeros(layers[i+1]) for i in range(len(layers)–1)]

def forward(X):
a, cache = X, []
for i in range(len(W)):
inp = a
z = inp @ W[i] + B[i]
a = np.tanh(z) if i < len(W)–1 else z
cache.append((inp, z, a))
return a, cache

def backward(cache, dout):
dW = [None]*len(W); dB = [None]*len(W)
delta = dout
for i in reversed(range(len(W))):
inp, z, _ = cache[i]
dW[i] = inp.T @ delta
dB[i] = delta.sum(0)
if i > 0:
delta = (delta @ W[i].T) * (1 – np.tanh(cache[i–1][1])**2)
return dW, dB

def burgers(x, t):
"""用中心差分拿 u_x, u_t, u_xx，拼出 PDE 残差 R。"""
x = x.reshape(–1, 1); t = t.reshape(–1, 1)
X0 = np.hstack([x, t])
Xxp = np.hstack([x+H, t]); Xxn = np.hstack([x–H, t])
Xtp = np.hstack([x, t+H]); Xtn = np.hstack([x, t–H])
u0, c0 = forward(X0); uxp, cxp = forward(Xxp); uxn, cxn = forward(Xxn)
utp, ctp = forward(Xtp); utn, ctn = forward(Xtn)
ux = (uxp – uxn) / (2*H)
ut = (utp – utn) / (2*H)
uxx = (uxp – 2*u0 + uxn) / H**2
r = ut + u0*ux – nu*uxx # 这就是要压到 0 的 PDE 残差
return r, u0, ux, c0, cxp, cxn, ctp, ctn

xr = np.random.uniform(–1, 1, 64) # 配点（collocation）
tr = np.random.uniform(0, 1, 64)
x0 = np.random.uniform(–1, 1, 32) # 初始条件点
u0t = –np.sin(np.pi * x0)

mW = [np.zeros_like(w) for w in W]; vW = [np.zeros_like(w) for w in W]
mB = [np.zeros_like(b) for b in B]; vB = [np.zeros_like(b) for b in B]
b1, b2, eps, lr = 0.9, 0.999, 1e-8, 5e-3

for step in range(1, 4001):
r, u0, ux, c0, cxp, cxn, ctp, ctn = burgers(xr, tr)
N = r.shape[0]
g = 2 * r / N
dW = [np.zeros_like(w) for w in W]; dB = [np.zeros_like(b) for b in B]
pairs = [(c0, g*(ux + 2*nu/H**2)),
(cxp, g*(u0/(2*H) – nu/H**2)),
(cxn, g*(–u0/(2*H) – nu/H**2)),
(ctp, g*(1/(2*H))),
(ctn, g*(–1/(2*H)))]
for cache, dout in pairs:
dWi, dBi = backward(cache, dout)
for i in range(len(W)):
dW[i] += dWi[i]; dB[i] += dBi[i]
# 初始条件 u(x,0) = -sin(pi x)
X0 = np.hstack([x0.reshape(–1, 1), np.zeros((len(x0), 1))])
uic, cic = forward(X0)
if step % 1000 == 0:
print(f"step {step:4d} loss_pde={np.mean(r**2):.3e} "
f"loss_ic={np.mean((uic.reshape(–1)–u0t)**2):.3e}")
dic = 2 * (uic – u0t.reshape(–1, 1)) / len(x0)
dWi, dBi = backward(cic, dic)
for i in range(len(W)):
dW[i] += dWi[i]; dB[i] += dBi[i]
for i in range(len(W)):
mW[i] = b1*mW[i] + (1–b1)*dW[i]; vW[i] = b2*vW[i] + (1–b2)*dW[i]**2
W[i] -= lr * (mW[i]/(1–b1**step)) / (np.sqrt(vW[i]/(1–b2**step)) + eps)
mB[i] = b1*mB[i] + (1–b1)*dB[i]; vB[i] = b2*vB[i] + (1–b2)*dB[i]**2
B[i] -= lr * (mB[i]/(1–b1**step)) / (np.sqrt(vB[i]/(1–b2**step)) + eps)

# 验证：训练后 t=0 处的预测应逼近 -sin(pi x)
xt = np.linspace(–1, 1, 9)
Xv = np.hstack([xt.reshape(–1, 1), np.zeros((len(xt), 1))])
pred, _ = forward(Xv)
print("x :", np.round(xt, 2))
print("预测:", np.round(pred.reshape(–1), 3))
print("真值:", np.round(–np.sin(np.pi*xt), 3))

运行结果：

step 1000 loss_pde=7.442e-03 loss_ic=1.002e-02
step 2000 loss_pde=6.862e-03 loss_ic=2.476e-03
step 3000 loss_pde=4.951e-04 loss_ic=5.280e-04
step 4000 loss_pde=1.750e-03 loss_ic=3.173e-04
x : [-1. -0.75 -0.5 -0.25 0. 0.25 0.5 0.75 1. ]
预测: [ 0.109 0.727 0.989 0.71 -0.026 -0.752 -0.985 -0.723 -0.033]
真值: [ 0. 0.707 1. 0.707 -0. -0.707 -1. -0.707 -0. ]

可以看到：训练前网络输出只是一团噪声；4000 步之后，$t=0$ 处的网络预测已经非常贴近真值 $-\sin (\pi x)$，残差损失从 $10^{-2}$ 一路被压到 $10^{-3}$ 乃至更低。这说明「把方程塞进损失」确实能把网络往物理解上拉。注意这里只用了 64 个配点和 32 个初始点，没有网格、没有手工解析引导项——这正是 PINN 最迷人的地方：物理规律本身成了最强的监督信号。

四、关键经验/避坑

五、完整系列推荐

📚 本文选自《可微分编程与科学计算》100 期系统教程（第 041 期：PINN 原理），每期配可运行 Python 代码。

完整系列（100 期正文 + 3 篇番外，每期文章+代码）已在 ima 知识号【Kruptos】持续更新：

• 🗂 70+ 技术知识库：操作系统、图神经网络、强化学习、数据库系统、推荐系统、编译原理……几乎覆盖全部软硬件技术栈
• 🧠 8 款 AI 技能：系列生产、知识库管理、CMMI 受管开发、自进化 Agent 等，已在 ima 技能广场上架，即装即用
• ✅ 全部免费订阅，后续更新自动推送

🔍 订阅方式：打开 ima（腾讯智能工作台）→ 搜索「Kruptos」→ 一键订阅；或在 ima 内直接搜索《可微分编程与科学计算》。

作者：Kruptos（西电毕业，13 年无线通信/DSP/嵌入式科研，现深耕 AI 与云原生）
原创内容，转载注明出处。