最近做结构力学小实验时碰到一个有意思的问题一根悬臂梁一端固定在墙上另一端悬空在均布载荷作用下弯曲怎么用神经网络把它的挠度曲线求出来用Python写一个PINN物理信息神经网络求解一维悬臂梁挠度问题听起来挺唬人其实核心思路就是把欧拉-伯努利梁方程作为惩罚项塞进损失函数剩下的交给自动微分和优化器。这篇博文会给出完整可运行的Python代码以及我在调参过程中实打实踩过的坑。这个算例的妙处在于它不需要任何真实标签数据只需要控制方程和边界条件网络就能自己“长”出挠曲线。特别适合两类人一是刚接触PINN的工科生想找一个物理意义清晰、随时能跑通的最小Demo二是结构分析方向的同学想看看无网格、无标签数据、纯物理约束的求解方式到底怎么落地。读完你应该能掌握一套从控制方程到损失函数、再到训练验证的完整工作流以后换成梁、板、流体问题套路都是相通的。1. 什么是PINN把物理方程写进损失函数1.1 为什么悬臂梁挠度问题也能用神经网络来解传统算悬臂梁挠度无非两条路查公式或者上有限元。查公式方便但只适用于等截面、规则载荷的简单情况有限元通用但需要画网格、定单元类型、组集刚度矩阵前处理一套流程下来至少小半天过去了。PINN走的是第三条路把物理规律当成一种“软约束”让神经网络自己逼近微分方程的解。具体到悬臂梁问题物理规律就是四阶常微分方程EI·w(x) q(x)。PINN用一个神经网络来表示挠度函数w(x)输入是坐标x输出是挠度预测值。然后对输出求各阶导数代回控制方程计算残差。残差越小说明网络输出越接近真实解。整个过程不需要网格生成也不需要单元积分本质上是一个无网格方法。有人可能会问神经网络不是需要大量标签数据吗这里恰恰相反PINN一点真实数据都不要它只靠两个东西控制方程残差和边界条件。方程和边界就是它的“教师”训练就是让网络在所有采样点上都尽量满足这两个物理约束。用个生活化的类比神经网络像个来应聘的人物理方程是面试官的评分表答得越符合物理规律分数越高。1.2 自动微分神经网络的“求导超能力”PINN能成立靠的是深度学习框架里的自动微分技术。传统数值解法求导要么手推解析式要么用有限差分近似前者费脑后者有截断误差。而神经网络本质上是一个多层复合函数每层都可导框架通过链式法则能自动算出任意阶导数精度能达到机器精度。在PyTorch里这一句torch.autograd.grad就能让网络输出对输入求导。想求四阶导数连续调用四次grad就行。这在十年前是不可想象的——那时候求一个复杂函数的四阶导数得自己推导半天现在框架帮你全干了。记住这句话PINN的“物理信息”就是靠自动微分把微分方程嵌入到网络训练中的。理解这点看所有PINN代码都不会懵。2. 悬臂梁挠度的数学建模控制方程与边界条件2.1 控制方程从材料力学到四阶微分方程悬臂梁挠度问题材料力学里对应的是欧拉-伯努利梁理论。梁的挠曲线w(x)满足四阶常微分方程EI * w(x) q(x)其中E是弹性模量I是截面惯性矩二者乘积EI称为抗弯刚度q(x)是单位长度上的分布载荷方向与挠度w的正方向一致。本文按这个无符号差异的列式来写不同教材里可能会因为坐标方向约定差一个负号但只要你把解析解和PINN控制方程保持同一约定最后结果就是一致的。这个四阶方程不是天上掉下来的它的物理链条很清晰弯矩M与挠度曲率的关系是M EI·w剪力V是弯矩沿梁长的变化率V M EI·w而分布载荷q又是剪力的变化率q V EI·w。从载荷到剪力、从剪力到弯矩、从弯矩到转角、从转角到挠度一共经历了四次求导或积分所以最终方程是四阶的。这就是为什么悬臂梁挠度问题需要“四阶”导数——物理层级决定的不是数学上故意刁难。2.2 四个边界条件的物理含义四阶微分方程需要四个边界条件才能定解。悬臂梁两端各有各的约束固定端x0处挠度为零w(0)0转角为零w(0)0。自由端xL处弯矩为零即w(L)0剪力为零即w(L)0。这四个条件对应着悬臂梁的真实物理状态固定端被死死卡住既不能下沉也不能转动自由端没有支撑也没有外力作用所以弯矩和剪力都为零。在PINN里边界条件不像有限元那样“强制施加”而是变成损失函数里的惩罚项让网络主动去满足。这也是PINN软约束思想的体现——每个边界条件都是一条平方误差网络在训练中会自动权衡满足程度。需要注意一个细节固定端和自由端的边界条件类型不同固定端给的是“函数值和一阶导为零”自由端给的是“二阶导和三阶导为零”。在代码实现时必须把它们正确映射到对应阶数的导数上。这个映射错了结果一定不对。3. 环境准备与无量纲化处理3.1 依赖库与运行环境这次实现只依赖三个库都非常常见torch深度学习框架提供自动微分和优化器numpy数值计算主要用于生成采样点matplotlib可视化画挠度曲线和误差曲线。版本上没有特殊要求PyTorch 1.13以上、2.x都行。最友好的一点是整个训练全程CPU即可完成不需要GPU。因为网络结构很浅三层隐藏层、配点数也不多单次迭代毫秒级总共训练两万步也就一两分钟。这算是我喜欢用这个算例入门PINN的原因之一——门槛低跑得快随时能试。3.2 无量纲化为什么我把均布载荷设成8而不是1理论上悬臂梁参数可以任意取但如果你直接把真实工程数值扔进网络大概率会碰到收敛问题。原因很简单神经网络对“量级”非常敏感。比如实际中E取200GPa、截面惯性矩I取个微小值挠度量级可能只有毫米甚至微米网络输出初始量级却可能是0.1或1这种水平二者差距过大损失函数就像在“跷跷板”上失衡训练会非常挣扎。所以我在这个例子里做了无量纲化直接取一组让解析解很“整齐”的数值参数取值含义L1.0梁长度无量纲EI1.0抗弯刚度无量纲q8.0均布载荷无量纲这里q取8不是随手写的有讲究。悬臂梁在均布载荷下的自由端最大挠度为qL⁴/(8EI)代入L1、EI1、q8最大挠度正好等于1。这样网络输出的目标范围就是0到1之间和初始权重初始化出的输出量级基本对齐训练难度大幅降低。如果q取1自由端最大挠度只有0.125输出量级太小网络需要很长时间才能把预测值压到这么小的范围还会和边界条件损失产生竞争收敛又慢又容易卡在错误解上。这是一个非常实用的“物理量级对齐”技巧把算例的期望输出量级调整到1附近比调网络结构、调学习率都管用。你在处理自己问题的时候不妨先算一下目标解的典型幅值再决定怎么归一无量纲化。4. 基于PyTorch的PINN完整实现4.1 网络模型与初始化网络结构我用了一个非常朴素的三层MLP每层64个神经元输入输出都只有一维import torch import torch.nn as nn import numpy as np import matplotlib.pyplot as plt # 随机种子 torch.manual_seed(42) np.random.seed(42) # 无量纲化后的物理参数 L 1.0 # 梁长度 EI 1.0 # 抗弯刚度 q 8.0 # 均布载荷集度class PINN(nn.Module): def __init__(self): super().__init__() self.net nn.Sequential( nn.Linear(1, 64), nn.Tanh(), nn.Linear(64, 64), nn.Tanh(), nn.Linear(64, 64), nn.Tanh(), nn.Linear(64, 1) ) def forward(self, x): return self.net(x) model PINN()激活函数我选了Tanh这里先留个悬念后面第5章我会重点讲为什么ReLU在这种高阶导数问题里完全不行。4.2 四阶导数计算autograd.grad的写法与坑PINN中最关键、也最容易写错的一步就是用自动微分求网络输出对输入的四阶导数。直接看代码def derivs(model, x): x x.clone().requires_grad_(True) w model(x) ones torch.ones_like(w) dw torch.autograd.grad(w, x, grad_outputsones, create_graphTrue, retain_graphTrue)[0] d2w torch.autograd.grad(dw, x, grad_outputstorch.ones_like(dw), create_graphTrue, retain_graphTrue)[0] d3w torch.autograd.grad(d2w, x, grad_outputstorch.ones_like(d2w), create_graphTrue, retain_graphTrue)[0] d4w torch.autograd.grad(d3w, x, grad_outputstorch.ones_like(d3w), create_graphTrue, retain_graphTrue)[0] return w, dw, d2w, d3w, d4w这里面有两个坑必须强调。第一个坑create_graphTrue绝对不能不写。因为我们要对四阶导数继续做反向传播也就是计算损失对网络参数的梯度这要求每次求导都保留计算图。如果你漏了create_graphTrue求二阶导数时图就断了后面的反向传播直接报错。第二个坑grad_outputs的维度要和被求导张量一致一般用torch.ones_like(w)就行。这其实是雅可比向量积的标准写法理解成“给每个输出分量一个权重再按链式法则把梯度传回去”即可不用深挖。4.3 损失函数方程残差与边界约束的加权损失函数分两部分控制方程残差和边界条件残差。控制方程残差在所有内部配点上计算边界条件残差在梁两端各两个点上计算。# 固定端 x0自由端 xL x_left torch.tensor([[0.0]]) x_right torch.tensor([[L]]) # 配点均匀分布 微小扰动让采样更自然 x_collocation torch.linspace(0, L, 200).reshape(-1, 1) x_collocation torch.clamp( x_collocation 0.01 * torch.randn_like(x_collocation), 0, L ) def loss_fn(model, x_col): w, dw, d2w, d3w, d4w derivs(model, x_col) # 控制方程残差EI*w - q 0 loss_eq torch.mean((EI * d4w - q) ** 2) # 固定端边界w0, w0 wL, dwL, _, _, _ derivs(model, x_left) loss_fixed wL[0] ** 2 dwL[0] ** 2 # 自由端边界w0, w0 _, _, d2wR, d3wR, _ derivs(model, x_right) loss_free d2wR[0] ** 2 d3wR[0] ** 2 loss_bc loss_fixed loss_free total loss_eq 10.0 * loss_bc return total, loss_eq, loss_bc关于权重10.0这是PINN调参里非常常见的手感。边界条件只靠两个点约束在总损失中的天然占比很小如果不加重网络会先“顾着”满足内点方程边界收敛得很慢。把边界权重加到10左右能让网络在训练前期就先把边界条件拉满然后再慢慢调内点。很多PINN实现会进一步用自适应权重但在这里10倍边界惩罚已经够用。4.4 训练循环Adam加学习率衰减最后L-BFGS精修优化器我用了两段式先Adam快速把损失拉下来再用L-BFGS做精细优化。这是PINN社区公认的稳定套路。optimizer torch.optim.Adam(model.parameters(), lr1e-3) scheduler torch.optim.lr_scheduler.StepLR(optimizer, step_size5000, gamma0.5) for it in range(20001): optimizer.zero_grad() total, loss_eq, loss_bc loss_fn(model, x_collocation) total.backward() optimizer.step() scheduler.step() if it % 2000 0: print(fiter {it:5d} total {total.item():.3e} feq {loss_eq.item():.3e} bc {loss_bc.item():.3e})# Adam结束后用L-BFGS精修几十步 optimizer_lbfgs torch.optim.LBFGS(model.parameters(), lr0.1, max_iter100) def closure(): optimizer_lbfgs.zero_grad() total, _, _ loss_fn(model, x_collocation) total.backward() return total optimizer_lbfgs.step(closure)为什么最后要切L-BFGSAdam用的是随机梯度下降的惯性思路到后期损失接近极小值时它会在极小点附近来回震荡很难继续压低损失。L-BFGS是一种拟牛顿法利用损失函数的曲率信息收敛路径更“直”特别适合这种小规模的平滑优化问题。实测中Adam跑到一定精度后再让L-BFGS跑几十步误差往往能再降一个数量级。4.5 与解析解对比的验证脚本悬臂梁均布载荷下的解析解是经典公式直接拿来当“标准答案”w(x) q * x^2 * (x^2 - 4*L*x 6*L^2) / (24*EI)验证代码很简单def analytic_w(x): return q * x**2 * (x**2 - 4*L*x 6*L**2) / (24.0 * EI) x_test torch.linspace(0, L, 200).reshape(-1, 1) x_test_req x_test.clone().requires_grad_(True) w_pred, _, _, _, _ derivs(model, x_test_req) w_pred w_pred.detach().numpy() x_plot x_test.numpy() w_exact analytic_w(x_plot) err_max np.max(np.abs(w_pred - w_exact)) rel_err np.linalg.norm(w_pred - w_exact) / np.linalg.norm(w_exact) print(fmax error: {err_max:.3e}) print(frelative L2 error: {rel_err:.3e})fig, axes plt.subplots(1, 2, figsize(10, 4)) axes[0].plot(x_plot, w_exact, k-, labelAnalytic) axes[0].plot(x_plot, w_pred, r--, labelPINN) axes[0].set_xlabel(x) axes[0].set_ylabel(w) axes[0].legend() axes[1].plot(x_plot, np.abs(w_pred - w_exact), b-) axes[1].set_xlabel(x) axes[1].set_ylabel(absolute error) plt.tight_layout() plt.show()按我的经验这套配置跑完L2相对误差大体能落到1e-4量级左右自由端附近的最大绝对误差可能稍大但视觉上曲线和解析解几乎完全重合。对一个不用任何标签数据的纯物理求解来说这个精度已经相当能打。5. 调参与收敛性分析我实际跑出来的经验5.1 边界条件权重为什么10倍惩罚最顺手我一开始做这个算例时边界权重用的是1结果训练几千步后发现固定端挠度不为零自由端也有明显偏差。打印各分项损失才知道边界损失只占总量极少一部分网络只顾着压内点方程残差边界被忽略了。把边界权重提到10后情况立刻好转再提到50边界是满足了但内点残差又变大训练出现“顾此失彼”。最终我用10作为默认值平衡度最好。这也引出一个PINN调参的通用方法训练时一定分开打印各项损失不要只看总损失。如果哪项降不下去就给哪项加权重。5.2 激活函数与网络宽度为什么ReLU在高阶导数问题里直接废掉激活函数的选择对PINN是致命的。我在这个算例里试过ReLU结果完全没法用——ReLU的二阶导数恒为0三阶四阶更是0网络在数学上根本表达不出“弯曲”的信息控制方程残差永远压不下去。你想想ReLU的导数只在0处跳变二阶导数就是一堆脉冲四阶导数是脉冲的脉冲这对光滑挠度曲线来说毫无意义。Tanh的好处是无穷阶可导且导数在训练初期不容易饱和。Sigmoid也可导但收敛速度明显慢于Tanh。Sin激活函数在周期性问题里很强但悬臂梁这类单调光滑曲线用Sin有点杀鸡用牛刀而且前期更容易振荡。综合下来Tanh是求解梁类问题最稳的默认选择。网络宽度方面我试过从32到128神经元三到四层隐藏层结论是64×3就够用了。更宽的网络不会带来明显精度提升反而增加训练时间和高阶导数的数值波动。5.3 配点数、学习率与典型收敛趋势配点数量影响也不大从100试到2000最终精度几乎一样只是2000个点训练会慢一些。我最终取200个配点已经能让控制方程残差在整根梁上平均分布。关键是配点不能只在固定点采样加一点随机扰动会让采样更自然还能稍微提升边界附近的稳定性。学习率上初始1e-3配StepLR每5000步衰减0.5是稳妥的。如果你发现训练前期损失不降可以试试1e-4如果发现损失震荡厉害说明学习率偏大。典型的收敛过程是前500步总损失从几十快速降到个位数2000步左右降到1e-110000步后降到1e-3量级最后L-BFGS精修后总损失可以到1e-5甚至更低。6. 常见问题与排查技巧实录6.1 网络输出恒为0、损失卡住不动这是新手最容易碰到的“假死”现象。表面上看损失在缓慢下降但画出来的挠度曲线永远贴着0轴。这个时候先检查两件事一是输出量级是否和目标挠度严重不匹配建议按第3章的方式调整无量纲参数让最大目标挠度接近1二是边界损失权重是否太低导致网络只满足方程却不管边界把所有内点都“压”到0附近去降低方程残差。还有一种可能网络初始化后输出本来就接近0而方程残差项又在初始阶段主导了梯度方向。解决方法是把q适当调大让方程残差的量级和边界损失的量级尽量可比。6.2 四阶导数震荡、梯度爆炸、出现NaN高阶导数的数值刚性是PINN的老大难。网络输出稍微波动一下四阶导数就会被放大导致损失很大、梯度爆炸。常见解决三板斧对输入x做归一化让特征范围落在0到1之间这是最基础的一步网络深度不要太深层数越多高阶级联越大越容易炸三层足够了给优化器加梯度裁剪或者把学习率调小到1e-4量级。如果仍然震荡最有效的兜底方案是切换激活函数——把Tanh训练好的权重拿到Sin网络里继续微调这个迁移技巧偶尔会救场但我用得不多因为工程上更推荐下面的降阶方案。6.3 高阶导数求解的降阶替代方案直接硬解四阶方程虽然能跑通但四阶导数的数值刚性始终是隐患。一个工程上更稳的做法是“降阶”把四阶问题拆成两个二阶问题。引入辅助变量M表示弯矩原方程变成EI * w -M M -q这样网络需要输出两个值挠度w和弯矩M。损失函数包含两组方程残差和四个边界条件。二阶导数的数值稳定性远好于四阶导数训练起来明显更顺滑。代价是网络输出变成了二维边界条件的写法也要跟着调整但代码逻辑和多输出网络很常见改动不大。如果你的问题涉及变截面梁、复杂载荷或者想进一步扩展到逆问题我建议直接上手降阶版本省心。6.4 常见问题速查表现象可能原因排查与修复思路损失不下降学习率不合理、边界权重过低打印分项损失调lr加大边界权重输出始终为0目标挠度量级太小、方程与边界失衡无量纲化到最大挠度≈1调整q四阶导数震荡高阶导数数值刚性、网络过深输入归一化减少层数加梯度裁剪边界满足但内点误差大内点方程权重低、配点稀疏提高loss_eq权重增加配点数训练后期loss抖动Adam步长无法继续收敛切L-BFGS精修或继续衰减学习率NaN梯度爆炸降学习率、降低网络宽度、换Tanh结尾这个算例还能怎么用从悬臂梁这个最小算例入手我发现PINN真正的价值不在于和成熟有限元软件比精度、比效率——那没有意义。它最打动我的地方是能把“物理规律”变成一行可微分的损失函数让网络在没有任何实验数据的前提下自己学会力学行为。以后遇到变截面梁就把EI从常数改成x的函数方程残差那行代码改几个字母而已遇到集中力载荷就把q改成脉冲函数想反推载荷分布把均布载荷q设为可训练参数再加几组实测挠度数据项这就是一个标准的参数识别算例。我在实际使用中最深的体会是遇到四阶问题不要慌也别急着堆网络复杂度先检查量级对齐和损失权重这两个因素对结果的影响远大于网络结构本身。希望这份代码和踩坑记录能帮你少走点弯路。如果你也跑出了不错的收敛曲线或者碰到了我没提到的奇怪现象欢迎按这套流程继续折腾调参的乐趣往往就藏在这些看似不起眼的细节里。