物理信息神经网络PINN这两年几乎成了科学计算领域最热的词之一尤其在偏微分方程求解、流动建模和反问题方向。以前我做流动计算遇到强非线性方程第一反应是加密网格、换格式、调时间推进方案后来认真把PINN跑通才发现深度学习和物理定律结合之后能做的事比想象中多很多。这篇文章我会分享一套基于MATLAB的PINN求解方案用Burger方程作为第一个算例再延伸到其他复杂流动模型。如果你手上有MATLAB但又不想专门去搭一套Python深度学习环境这篇文章的路子可以直接照着走。先说清楚这个内容能解决什么问题。传统数值方法在网格生成、边界处理和参数扫描上往往要花大量人力PDE换个初边值条件很可能就要重新生成一套流程。PINN的思路刚好相反把偏微分方程的残差直接写进损失函数让神经网络在拟合数据的同时遵守物理规律。它尤其适合参数化研究、稀疏数据反演这类场景。文章后面会用Burger方程跑通完整流程再给出向Navier-Stokes等复杂流动模型扩展的具体思路。1. 先搞清楚PINN到底在解什么方程1.1 把物理方程写进损失函数而不是写进网络结构很多新手第一次接触PINN总以为它是在网络结构上做了什么手脚让网络天生满足物理定律。其实不是。PINN的核心思路非常朴素你有一个待求函数u(x,t)用一个神经网络u_θ(x,t)去逼近它然后在损失函数里显式地惩罚它偏离物理方程的程度。以Burger方程为例u_t u * u_x ν * u_xx如果u_θ是精确解那残差项 r u_t u * u_x - ν * u_xx 应该处处等于零。PINN的做法就是对整个时空域采样一堆点算这个残差然后让残差的平方均值尽可能小。再加上初始条件的拟合误差和边界条件的拟合误差一起作为总的损失函数。比较关键的一点是网络训练过程中需要的导数并不是用有限差分近似出来的而是通过自动微分Automatic Differentiation精确计算的。MATLAB里的dlarray和dlgradient就是干这个用的。自动微分沿着网络的计算图链式求导所以u_t、u_x、u_xx都是机器算出来的不是我们手推的。这个思路用生活化的说法解释网络可以自由变形但每次变形都要交“物理罚金”偏离方程越远罚金越重。训练的本质就是在有限的数据约束和物理罚金之间找一个平衡点。1.2 为什么我用MATLAB而不是Python框架提到PINN很多人第一反应是Python里的各种开源框架生态确实丰富。但我最后选择用MATLAB有几个很现实的原因。第一很多课题组和工程师手里已经有大量MATLAB代码数据预处理、后处理、可视化全在MATLAB里。如果单独为PINN搭一套Python环境数据要导来导去光格式转换就能折腾半天。MATLAB深度学习工具箱从R2019b开始支持dlnetwork后续版本又完善了自定义训练循环写一个完整的PINN并不困难。第二MATLAB调试体验好。PINN训练不收敛的时候经常需要检查中间某个张量的数值、看某一层的输出分布Python里你要么print要么上debuggerMATLAB工作区里点开变量一目了然这对从传统科学计算转过来的研究者特别友好。第三PINN这类问题网络规模通常不大不需要上千层的Transformer也不需要超大规模分布式训练。一个小型全连接网络在MATLAB里用CPU跑几千步迭代耗时也就几分钟到几十分钟。GPU有更好没有也不会是瓶颈。当然MATLAB不是万能的。如果你要做非常复杂的算子库、大规模GPU并行或者想直接用社区里现成的PINN组件Python生态确实更丰富。我的建议是项目数据和团队习惯在MATLAB就用MATLAB不要为了“时髦”额外引入一套技术栈。1.3 PINN擅长和不擅长的问题边界我在实际项目里踩过不少坑发现PINN不是万能药。先说说它擅长什么。高维问题是一个典型场景。传统数值方法在高维空间会遇到网格爆炸比如五维六维的PDE网格法基本没法做。PINN用点采样代替网格剖分只要采样点覆盖得住高维问题也能尝试。还有就是参数化代理模型这个概念后面会细说——把多个工况打包成一个网络推理时秒出解。再一个是反问题已知稀疏的场观测数据反推方程系数或源项这是PINN相比传统数值方法的突出优势。但PINN也有明显的短板。强激波、低粘性这类问题解里存在近间断的陡峭剖面网络很难用一组光滑基函数精确表达容易在激波附近产生抹平或振荡。超高雷诺数湍流直接硬解目前也不太现实。另外对于简单的标准PDE传统有限差分或有限元无论在精度还是计算效率上往往都更胜一筹。所以我的建议是不要一上来就想着用PINN替代传统CFD而是把它当成第二套工具在参数扫描、反演、数据融合这些传统方法比较吃力的场景里发挥作用。2. 用MATLAB搭建Burger方程的完整PINN求解器2.1 方程设定与参考解准备Burger方程是PINN入门最经典的测试用例非线性对流项加粘性扩散项既能体现出流动问题的难度又不像完整Navier-Stokes方程那样复杂得让人劝退。我采用的设定是x ∈ [-1, 1], t ∈ [0, 1] u_t u * u_x ν * u_xx u(0, x) -sin(πx) u(t, -1) u(t, 1) 0粘性系数ν取0.01。这个参数下有较明显的激波形成过程对网络有一定挑战性但对展示PINN能力又足够。如果你第一次跑我建议先用ν 0.1试一遍解会更平滑容易把整个流程跑通等确认代码没有问题再切回ν 0.01。做PINN训练之前先准备一套参考解用于事后对比。Burger方程理论上可以通过Cole-Hopf变换求解析解但实际实现还是要做数值积分。更省事的方式是直接用高分辨率守恒型差分格式算一个基准解。我用的是迎风加中心差分的混合格式空间网格取2000个点时间步数取4000步得到的解足够精细。示意图大致是这样% 高分辨率差分求解Burger方程作为参考解示意代码 % 注意这里为了展示思路做了简化正式使用请采用能达到网格无关的格式 nu 0.01; Nx 2000; Nt 4000; x linspace(-1, 1, Nx); dx x(2) - x(1); dt 1 / Nt; u -sin(pi * x); for n 1:Nt % 对流项采用迎风思路粘性项用中心差分 flux 0.5 * u.^2; % 一阶迎风/中心混合实际需要CFL判断 u(2:end-1) u(2:end-1) - dt/dx * (flux(2:end-1) - flux(1:end-2)) ... nu * dt/dx^2 * (u(3:end) - 2*u(2:end-1) u(1:end-2)); % 边界条件 u(1) 0; u(end) 0; end这个参考解不要求达到论文级精度只要网格够细、格式稳定能给出可靠的对比基准就够了。2.2 网络定义与自动微分残差在MATLAB里定义网络使用dlnetwork。我常用的结构是4层隐藏层、每层50个神经元激活函数用tanh。layers [ featureInputLayer(2, Normalization, none) fullyConnectedLayer(50) tanhLayer fullyConnectedLayer(50) tanhLayer fullyConnectedLayer(50) tanhLayer fullyConnectedLayer(1) ]; net dlnetwork(layers);为什么要用tanh而不用ReLU原因很简单PINN的损失里要求u对x的二阶导ReLU是分段线性函数二阶导恒等于零网络就失去了表达粘性项的能力。tanh这类光滑激活函数是更合理的选择。损失函数是PINN的核心。下面的示意代码展示了如何用dlgradient计算一阶导和二阶导从而形成PDE残差function [loss, grads] modelLoss(net, X, T, X0, T0, U0, nu) % 网络输入是 [X; T]输出是u U forward(net, [X; T]); % 计算一阶和二阶导数 Ux dlgradient(U, X); Uxx dlgradient(Ux, X); Ut dlgradient(U, T); % PDE残差 f Ut U .* Ux - nu * Uxx; pdeLoss mean(f(:).^2); % 初始条件损失 U0_pred forward(net, [X0; T0]); icLoss mean((U0_pred(:) - U0(:)).^2); % 边界条件损失x-1 和 x1 处近似为0 % 实际实现要传入左右边界点这里示意 U_left forward(net, [X_left; T_left]); U_right forward(net, [X_right; T_right]); bcLoss mean(U_left(:).^2) mean(U_right(:).^2); % 总损失 loss pdeLoss icLoss bcLoss; grads dlgradient(loss, net.Learnables); end注意这里X、T、X0、T0都必须是dlarray而且格式要保持列对齐。比如残差点采样时X和T是同一条样本的横坐标和时间戳拼接成2×N的输入矩阵网络输出的列才和样本一一对应。2.3 采样策略与训练流程PINN不需要网格但需要“点”。点怎么采直接影响训练效果。残差点在x∈[-1,1]、t∈[0,1]的整个定义域内随机采样。我用lhsdesign拉丁超立方或者简单rand都行采样量Nf取10000。这些点的作用是让网络在全局范围满足PDE。初始点在t0时刻x方向上采样200个点u值是已知的由初始条件给出。边界点在x-1和x1的边界上各采样100~200个时间点u值等于0。采样代码大致如下rng(0); Nf 10000; X dlarray(-1 2 * rand(1, Nf), CB); T dlarray(rand(1, Nf), CB); Nic 200; X0 dlarray(-1 2 * rand(1, Nic), CB); T0 dlarray(zeros(1, Nic), CB); U0 dlarray(-sin(pi * X0), CB);训练的主循环用adamupdate。这里有个很实用的建议不要迷信单阶段长训练我推荐两阶段策略。先用一个相对大的学习率1e-2跑15000到20000步让loss快速下降再降到1e-3精调。两阶段搭配起来往往比始终用一个学习率收敛得更快、更稳。3. 训练实操与收敛细节3.1 拆开损失函数学会看训练曲线很多新手看到总loss开始下降就觉得万事大吉这是很危险的。总loss只是三个损失项的加权和同一个数值可能对应完全不同的内部状态。我建议在训练循环里把pdeLoss、icLoss、bcLoss分别打印出来。numIter 20000; learnRate 1e-2; avgGrad []; avgSqGrad []; for iter 1:numIter [loss, grads] dlfeval(modelLoss, net, X, T, X0, T0, U0, nu); [net, avgGrad, avgSqGrad] adamupdate(net, grads, avgGrad, avgSqGrad, iter, learnRate); if mod(iter, 1000) 0 fprintf(Iter %6d, loss %.3e, pde %.3e, ic %.3e, bc %.3e\n, ... iter, extractdata(loss), extractdata(pdeLoss), ... extractdata(icLoss), extractdata(bcLoss)); end end如果你看到pdeLoss一直压在1e-2量级下不去大概率是残差点不够密或者网络宽度不够如果icLoss总是比pdeLoss高一个量级说明初始点附近的拟合出了问题需要增加初始点密度或者对初始条件做硬约束。3.2 两阶段训练与学习率调整两阶段训练是PINN实战中很有效的手段。第一轮用1e-2的学习率跑20000步观察loss曲线。当loss下降变得非常缓慢时把学习率降到1e-3再跑10000到20000步。这个过程中不需要改其他参数。需要注意PINN的loss曲线不像传统深度学习那样稳定平滑经常会出现小幅震荡尤其是PDE残差点重新采样之后。遇到这种情况先看总体趋势不要因为单次跳变就急着调参。如果训练到一半loss突然升高最常见的元凶是学习率太大导致参数越过最优区域。这时候不要盲目重启先把学习率降一个量级继续从当前模型接着训练很多时候就能救回来。3.3 结果对比与误差判断训练完成后需要在稠密网格上推理出整个解曲面。方法很简单x_eval linspace(-1, 1, 400); t_eval linspace(0, 1, 200); [Xm, Tm] meshgrid(x_eval, t_eval); Xt dlarray([Xm(:); Tm(:)], CB); U_pred predict(net, Xt); U_pred reshape(extractdata(U_pred), size(Xm));接下来和参考解做三件事第一固定几个时间层画剖面线比如t0.25、0.5、0.75、1.0看解的激波位置和厚度是否和参考解基本吻合第二画整个时空域的预测热图直观感受解曲面第三算相对L2误差作为量化指标。从实际测试来看ν0.01的Burger方程PINN在激波附近通常会有轻微的抹平效应这是正常现象不必过度焦虑。如果发现激波附近误差特别大可以在激波轨迹附近增加采样密度或者把隐藏层宽度从50加到80。4. 从Burger方程走向复杂流动模型4.1 不可压缩Navier-Stokes方程的PINN形式从Burger方程迁移到Navier-Stokes方程核心思路完全一样只是网络输入和输出的维数变了。以二维不可压缩NS方程为例u_x v_y 0 u_t u*u_x v*u_y -p_x / ρ ν*(u_xx u_yy) v_t u*v_x v*v_y -p_y / ρ ν*(v_xx v_yy)网络输入从(x,t)变成(x,y,t)输出从一个标量u变成三个分量[u, v, p]。损失函数里除了动量方程的两个残差项还要加入连续性方程残差项。这里有一个非常值得注意的坑压力p没有独立的进化方程它整体可以加减一个常数而不影响速度场。如果训练过程中不固定压力参考点p的均值会漂移训练时经常会看到一个不断变化的常数偏移。解决方法是固定某一点的p值或者在网络输出端让p减去当前batch的平均值。初始条件和边界条件的处理方式与Burger方程一致都是构造拟合损失。如果几何比较复杂只需要在边界上采样点不需要生成贴体网格这是PINN相对于传统CFD方法的一大优势。4.2 参数化代理模型与多工况覆盖当网络结构确定之后一个很自然的扩展就是把方程里的参数也当成输入。比如想让网络覆盖多个粘性系数ν的工况可以把网络输入从(x,t)扩展到(x,t,ν)训练时每个batch随机抽取若干个ν值。这个做法对工程场景特别有意义。传统CFD要覆盖多个雷诺数需要逐个工况仿真每个工况都要重新生成网格、迭代求解。PINN参数化代理模型训练一次推理时输入不同的ν就能秒出对应流场。我第一次跑通这个模式时还挺惊喜的一个网络顶替掉了原先一堆重复算例。不过要注意参数跨度大时网络输入尺度差异会非常明显。建议先把ν这类参数做归一化到0到1否则网络会因为输入尺度过大而难以收敛。当然也可以把这类参数做log缩放后归一化效果往往更好。4.3 反问题与数据融合PINN更大的价值在反问题。很多时候我们拿不到完整的初始条件和边界条件但手里有部分稀疏的观测数据希望反推方程的未知参数这时候传统方法很棘手PINN却能自然地处理。做法是把方程中的未知系数比如粘性ν、某个源项强度定义成网络的可训练变量损失函数里除了PDE残差和初边值条件再加上观测数据的拟合误差。训练过程中网络权重和这个未知参数同时被Adam更新。训练结束后未知参数的取值就是反演结果。这里有个经验之谈观测数据项的权重通常要设得大一些比如10或者50否则物理残差会把数据拟合给“带偏”。如果数据本身有噪声权重还要根据信噪比适当调整。这种把实测数据和物理方程结合在一起的思路在工业测量、医疗影像、场重建等场景里非常实用。5. 常见问题与排查技巧5.1 训练不收敛时的排查路径我见过太多人第一步就栽在“loss怎么都降不下去”上。这里给出一个排查路径基本可以覆盖大多数情况。第一步先看网络输出有没有NaN。如果有八成是学习率太高或者输入没有归一化导致梯度爆炸。把学习率降到1e-3以下试试。第二步看总loss的内部构成。用前面说的拆解打印定位是pdeLoss、icLoss还是bcLoss的问题。三个都高大概率是采样点不够或者网络容量不足只有某一个高就针对性地处理。第三步看残差点的采样范围。确定有没有把时间和空间范围设错比如把t采到了[0,1]之外网络在定义域外学习了一堆没用的信息。第四步检查激活函数。如果用了ReLU别犹豫换成tanh。ReLU做PINN基本等于自断一臂。5.2 边界条件失效与激波附近误差过大边界条件不满足是很常见的问题。软约束方式下网络有可能为了降低PDE残差而牺牲边界拟合。如果你发现边界附近误差始终压不下来最简单的办法是加大边界损失的权重比如把bcLoss乘上一个5到10的系数。更彻底的办法是硬约束。比如Burger方程的边界条件是u(t,-1)u(t,1)0可以在网络输出后面做一个仿射变换u_out (1 - x)/2 * u_left(t) (1 x)/2 * u_right(t) (1 - x^2) * u_net这样在x-1处u_out自动等于u_left(t)在x1处自动等于u_right(t)边界条件被严格满足连边界损失项都可以去掉。这个技巧在很多PINN工程里非常实用尤其适合边界条件比较简单的情况。至于激波附近误差大我的习惯是先增大采样点数量并在激波轨迹附近加密采样。如果还不行再把网络宽度往上加一档。同时可以把PDE残差中粘性项的权重适当调高因为激波位置的二阶导数变化非常剧烈网络需要更强的监督才能捕捉。5.3 几个值得养成的调参习惯最后分享几个我自己总结的调参习惯。第一每次训练固定随机种子。PINN对采样点和初始化很敏感不固定随机种子两次训练曲线完全不同你会分不清是参数改好了还是运气问题。第二把各个损失分量分开打印。总loss不能说明任何问题拆开看才定位得了瓶颈。第三先用简单case验证流程再上难度。比如先跑ν0.1的平滑解确认整套代码没问题再挑战ν0.01的激波问题。不要一上来就挑战最难的配置否则代码有bug你都分不清是bug还是网络能力不足。第四训练过程中定期保存模型。我经常遇到一个模型跑了很久突然学习率没有调好loss跳上去的情况。如果中途保存了较好的一次结果至少不用全部重来。第五不要迷信网络越深越好。PINN的问题难点往往不在网络容量而在采样策略、损失函数权重和学习率。加宽一档如果没用优先检查损失函数的设计。写到这里说点我自己的实际体会。第一次在MATLAB里把PINN完整跑通时我并没有立刻替换掉原来的CFD流程而是把它当成“第二套工具”来用传统数值方法继续负责高精度基准解PINN负责快速预测、参数扫描和反问题。这种混合工作模式比单纯争论某种方法的输赢要实用得多。如果你也准备在自己的项目里试PINN我建议先用Burger方程把这条链路彻底走通再决定要不要往Navier-Stokes方程上迁移。思路一通后面的路会顺很多。