3. 对流-扩散方程的对流速度与扩散系数联合辨识

任务目标

本案例以含源项的一维稳态对流-扩散方程为物理模型，给定浓度场的离散观测数据，通过PINN同时反演未知的对流速度 v 与扩散系数 D。要求反演结果与预设真值相对误差小于5%，验证多参数联合辨识的可行性。

物理模型与正问题

控制方程为：

D * u''(x) - v * u'(x) + f(x) = 0,  x ∈ [0, 1]

边界条件：u(0) = 0，u(1) = 1。取源项 f(x) = 1，真值设为 v_true = 1.0，D_true = 0.2。该方程存在解析解，用于生成训练数据并验证反演精度。

方法要点

采用硬约束边界条件：令网络输出 u_net(x)，构造 u(x) = x + x(1-x)·u_net(x)，自动满足两端边界条件。损失函数由三部分构成：

  • 控制方程残差（在配点上计算）；
  • 观测数据拟合误差（在观测点上计算）；
  • 参数正则项（可选，用于增强稳定性）。

反演参数 v、D 作为可训练变量，与网络权重一同通过Adam优化器更新。

关键完整代码

import numpy as np
import torch
import torch.nn as nn

# ---------- 1. 生成合成观测数据（解析解） ----------
def exact_solution(x, v, D):
    # 解析解: u(x) = C1 + C2*exp(v*x/D) + x/v
    C2 = (1 - 1/v) / (np.exp(v/D) - 1)
    C1 = -C2
    return C1 + C2 * np.exp(v * x / D) + x / v

v_true, D_true = 1.0, 0.2
x_obs = np.linspace(0.02, 0.98, 40)          # 避开边界，模拟内部观测
u_obs = exact_solution(x_obs, v_true, D_true)
u_obs += 0.002 * np.random.randn(len(x_obs))  # 加入轻微噪声

# 转换为张量
x_obs_t = torch.tensor(x_obs, dtype=torch.float32).reshape(-1, 1)
u_obs_t = torch.tensor(u_obs, dtype=torch.float32).reshape(-1, 1)

# ---------- 2. 定义PINN网络 ----------
class PINN(nn.Module):
    def __init__(self):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(1, 40), nn.Tanh(),
            nn.Linear(40, 40), nn.Tanh(),
            nn.Linear(40, 40), nn.Tanh(),
            nn.Linear(40, 1)
        )
        # 待反演参数（初始猜测远离真值）
        self.v = nn.Parameter(torch.tensor(0.5, dtype=torch.float32))
        self.D = nn.Parameter(torch.tensor(0.5, dtype=torch.float32))

    def forward(self, x):
        return self.net(x)

model = PINN()
optimizer = torch.optim.Adam(model.parameters(), lr=0.005)

# ---------- 3. 损失函数 ----------
def loss_fn(model, x_obs_t, u_obs_t):
    # 硬约束边界: u(x) = x + x(1-x)*N(x)
    x_colloc = torch.linspace(0, 1, 200, requires_grad=True).reshape(-1, 1)
    N = model(x_colloc)
    u = x_colloc + x_colloc * (1 - x_colloc) * N

    # 控制方程残差
    u_x = torch.autograd.grad(u, x_colloc, grad_outputs=torch.ones_like(u), create_graph=True)[0]
    u_xx = torch.autograd.grad(u_x, x_colloc, grad_outputs=torch.ones_like(u_x), create_graph=True)[0]
    f = torch.ones_like(x_colloc)
    pde_res = model.D * u_xx - model.v * u_x + f
    loss_pde = torch.mean(pde_res**2)

    # 观测数据拟合
    N_obs = model(x_obs_t)
    u_pred = x_obs_t + x_obs_t * (1 - x_obs_t) * N_obs
    loss_data = torch.mean((u_pred - u_obs_t)**2)

    return loss_pde + loss_data

# ---------- 4. 训练 ----------
for epoch in range(3000):
    optimizer.zero_grad()
    loss = loss_fn(model, x_obs_t, u_obs_t)
    loss.backward()
    optimizer.step()
    if epoch % 500 == 0:
        print(f"Epoch {epoch:4d} | Loss {loss.item():.2e} | v={model.v.item():.4f} | D={model.D.item():.4f}")

# ---------- 5. 输出反演结果 ----------
v_inv = model.v.item()
D_inv = model.D.item()
print(f"\n反演结果: v = {v_inv:.4f} (真值 {v_true}), D = {D_inv:.4f} (真值 {D_true})")
print(f"相对误差: v: {abs(v_inv-v_true)/v_true*100:.2f}%, D: {abs(D_inv-D_true)/D_true*100:.2f}%")

运行方式与依赖

需要 numpy 与 torch（CPU版本即可）。将上述代码保存为 pinn_advection_diffusion.py，在终端执行：

python pinn_advection_diffusion.py

验证方法与判定标准

训练约3000轮后观察终端输出：

  • 损失值应持续下降并趋于稳定（约 1e-5 量级或更低）；
  • 反演得到的 v 与 D 应分别收敛至接近1.0和0.2；
  • 判定通过条件：v 相对误差 < 5%，D 相对误差 < 5%。

若 D 收敛偏慢，可适当增加训练轮数至5000或调整学习率为0.003。由于观测含噪声，最终误差在1%～3%范围内属正常现象。

常见问题

  • 参数不收敛：检查是否已对 v、D 设置 requires_grad=True（本代码通过 nn.Parameter 自动实现）；尝试增大配点数量或降低学习率。
  • 损失震荡：可对观测数据做归一化处理，或减小Adam的 lr 至0.001。
  • 辨识结果偏移：观测噪声过大时，可增加数据点数量或添加L2正则项约束参数范围。
← 上一章 📖 返回目录 下一章 →