物理信息神经网络求解悬臂梁挠度：原理、实现与调参实战

做结构分析的人对“悬臂梁挠度”这几个字应该都不陌生。以前算这类问题我第一反应是翻材料力学教材里的解析公式或者直接上有限元软件。但后来在一个非均布荷载、变截面梁的预研项目里传统方法绕来绕去都不够顺手我才认真把物理信息神经网络PINN用在悬臂梁挠度计算上。这篇内容就是把我用PINN求解一维悬臂梁挠度问题的完整思路、Python实现和踩坑过程整理出来适合正在学PINN但不想只跑教科书例子的朋友也适合做结构仿真想换一种求解思路的工程师。我尽量把代码写得能直接跑把损失函数、高阶导数、边界条件这些容易卡住的细节讲透最后附上验证结果和调参教训。说实话PINN在简单梁问题上不会比有限元快但它把“物理规律”直接塞进神经网络这个思路在东想西想的工程场景里真的很有价值。1. 为什么选PINN算悬臂梁挠度传统数值方法的问题与PINN的思路先说一下我遇到的场景。某结构仿真预研项目里需要频繁计算一根悬臂梁在复杂荷载下的挠度荷载不是简单的均布或集中力而是沿着梁长变化的分布力截面惯性矩也在变化。用解析公式做每换一次荷载分布就要重新推导积分非常容易出错用有限元做网格剖分、前处理、结果导出这套流程虽然成熟但每次参数变化都要重新建模迭代速度跟不上。这时候我想到PINN的核心思想把控制方程这里是欧拉-伯努利梁方程和边界条件写成损失函数的一部分让神经网络在拟合数据的同时必须满足物理规律。换句话说神经网络的输出就是挠度函数 w(x)我要求它在整个梁长范围内满足 EI·w(x) q(x) 这个控制方程在固定端满足位移和转角为零在自由端满足弯矩和剪力为零。求解过程变成优化问题不需要网格不需要解析推导荷载分布怎么变只需要改损失函数里的 q(x) 表达式就行。这个方法最有意思的地方在于它把“求解微分方程”变成了“训练神经网络”。传统数值方法是在离散点上逼近微分方程的解PINN则是在连续空间里找一个函数由神经网络参数化让残差尽量小。对悬臂梁这种一维问题来说PINN的优势不是精度更高而是灵活性换个荷载、改个截面参数重新训练一遍就行代码几乎不用动。我实际对比过同样的参数下PINN用CPU训练几千步单根梁的挠度结果和有限元解析解能对上误差在百分之一左右。对于需要快速试算的预研阶段这个精度完全够用。如果你只是要一根等截面梁在均布荷载下的精确挠度直接套公式更快但如果你想试试“不加网格也能解结构问题”的路子或者后续要把物理规律嵌入更复杂的智能计算流程PINN值得好好研究。2. 悬臂梁挠度问题的数学建模与PINN损失函数设计2.1 欧拉-伯努利梁方程四阶导数的物理含义先明确我们要求解的问题。一根长度为 L 的悬臂梁左端固定在墙体x0 处右端自由xL 处梁上作用沿长度分布的荷载 q(x)。在小变形线弹性假设下挠度 w(x) 由欧拉-伯努利梁方程控制EI·(d⁴w/dx⁴) q(x)这里 EI 是抗弯刚度对均质等截面梁来说是常数。四阶导数的物理含义是挠度的二阶导对应曲率弯矩三阶导对应剪力四阶导对应荷载集度。材料的力学行为被浓缩成这一个方程听起来有点抽象但细想很优雅——一根梁的变形完全由“荷载”和“刚度”决定两者通过四阶导数建立联系。对于悬臂梁边界条件是四类固定端 x0挠度为零 w(0)0转角为零 dw/dx(0)0自由端 xL弯矩为零 d²w/dx²(L)0剪力为零 d³w/dx³(L)0四个边界条件对应四阶微分方程需要的四个约束数学上刚好封闭。PINN的做法就是让神经网络同时满足方程和这四个边界条件。2.2 损失函数的四项组成如何让网络“既懂数据又懂物理”PINN的损失函数由四个部分组成。第一项是控制方程残差在梁内部采样若干点计算 EI·w(x) - q(x) 的平方和这一项强制网络输出满足平衡方程第二项是固定端位移残差w(0) 的平方第三项是固定端转角残差dw/dx(0) 的平方第四项是自由端边界残差d²w/dx²(L) 和 d³w/dx³(L) 的平方和。写成公式就是L L_pde λ_bc · (L_w0 L_θ0 L_M L_V)其中 λ_bc 是边界条件的权重系数通常要设得比方程残差权重更大。为什么因为我在实践中发现如果不给边界条件足够的权重神经网络很容易“聪明地”找到一个在内部满足方程、但在边界上完全跑偏的函数。边界条件是微分方程问题的“锚”锚不钉牢船就会漂。具体到我的参数取值用的是接近工程实际的量级参数符号取值梁长L3.0 m均布荷载集度q(x)5000 N/m暂取常数以便解析验证抗弯刚度EI1.0e6 N·m²理论最大挠度w(L)q·L⁴/(8EI) 0.0506 m这个算例最大的好处是有解析解可以对照w(x) q·x²·(6L² - 4Lx x²)/(24EI)。先用它验证PINN程序的正确性再换复杂荷载才有底气。2.3 为什么要用高阶导数从一阶导到四阶导的链式推导PINN实现里最需要小心的是高阶导数的计算。神经网络输出 w 是输入坐标 x 的函数PyTorch 的 autograd 机制可以自动求导但默认只求一阶。要得到四阶导数需要连续使用 autograd.grad每次求导时设置 create_graphTrue这样求导结果仍然可以继续求导形成高阶导数的计算图。实际写代码时为了防止手动写四层嵌套太乱我会写一个循环函数输入 u、x 和阶数 n返回对应的 n 阶导数。理论上输入输出都是一维时可以用张量操作但为了保险我用的是逐个样本的批量计算方式这样在采样点数量变化时不容易出错。这里有个容易踩的坑高阶导数会导致计算图非常深内存开销不小而且梯度值会变得很小。这也是为什么 PINN 训练比普通神经网络慢——它不仅要前向传播还要在计算图上反复求导。对一维问题这个开销可以接受但到了二维板壳问题就非常吃显存了。3. Python实现从网络结构到高阶导数求导的完整代码3.1 网络设计为什么选tanh而不是ReLU网络结构我用的很简单输入层 1 个神经元坐标 x三到四个隐藏层每层 50 个神经元输出层 1 个神经元挠度 w。激活函数选的是 tanh不是因为时髦而是试过之后发现它最适合这个任务。ReLU 系列的激活函数有个致命问题它在正区间是线性函数二阶导数直接为 0无法表达弯矩曲率。梁的方程需要四阶导数ReLU 类的网络根本算不出非零的四阶导。Swish、GELU 这类函数有高阶导数但数值稳定性不如 tanh。Sigmoid 也可以但它输出范围是 (0,1)对挠度这种可能跨越正负的量不太方便。tanh 输出范围 (-1,1)本身原点对称且导数可以递归地用自身表示求高阶导时数值很稳。我实际对比过 tanh 和 Swish同样的训练轮数tanh 的收敛更平滑最终误差也更小。代码结构如下import torch import torch.nn as nn import numpy as np class PINN_Beam(nn.Module): 用于悬臂梁挠度求解的物理信息神经网络 def __init__(self, n_hidden4, n_neurons50): super().__init__() layers [nn.Linear(1, n_neurons)] for _ in range(n_hidden - 1): layers.append(nn.Linear(n_neurons, n_neurons)) layers.append(nn.Linear(n_neurons, 1)) self.net nn.Sequential(*layers) self.act nn.Tanh() def forward(self, x): for i, layer in enumerate(self.net): if i len(self.net) - 1: x self.act(layer(x)) else: x layer(x) return x3.2 高阶导数与损失函数autograd.grad循环的完整写法接下来是最关键的部分——计算挠度对坐标的高阶导数。我封装了一个函数输入网络输出 u、坐标 x、目标阶数 n输出 u 对 x 的 n 阶导数def ddx(u, x, n): 计算 u 对 x 的 n 阶导数create_graphTrue 保证可以继续求导 for _ in range(n): u torch.autograd.grad( u, x, grad_outputstorch.ones_like(u), create_graphTrue, retain_graphTrue )[0] return u注意这里的两个参数create_graphTrue 是必需的否则求导结果不再参与后续梯度计算retain_graphTrue 是因为在同一个计算图上我们可能要多次求导不及时释放缓存会内存耗尽。损失函数的写法我按四个边界条件分别计算方便单独观察哪一项没有收敛def compute_loss(model, x_pde, x_fixed, x_free, L_const, EI_const, q_const): # 控制方程残差 w model(x_pde) d4w ddx(w, x_pde, 4) loss_pde torch.mean((EI_const * d4w - q_const) ** 2) # 固定端 w(0) 0, dw/dx(0) 0 w_fixed model(x_fixed) dw_fixed ddx(w_fixed, x_fixed, 1) loss_fixed torch.mean(w_fixed ** 2 dw_fixed ** 2) # 自由端 d2w/dx2(L) 0, d3w/dx3(L) 0 w_free model(x_free) d2w_free ddx(w_free, x_free, 2) d3w_free ddx(w_free, x_free, 3) loss_free torch.mean(d2w_free ** 2 d3w_free ** 2) # 组合边界权重系数设置 lambda_bc 10.0 loss loss_pde lambda_bc * (loss_fixed loss_free) return loss, loss_pde, loss_fixed, loss_free我把边界权重系数设成 10.0这个数不是随手写的。权重如果跟 PDE 残差同一量级比如都是 1.0边界条件需要 8000 轮迭代才能勉强满足设成 10 之后大约 2000 轮就能看到固定端位移降到 1e-4 以下。但也不能设成 100 以上否则网络会把所有精力放在边界上导致梁内部的荷载-挠度关系严重失真。这个平衡是 PINN 调参最核心的手感后面专门展开讲。3.3 训练主循环采样策略与Adam优化器训练过程就是反复进行“采样-计算损失-反向传播”。每次迭代从梁内部随机采样一批点我常用 256 个计算 PDE 残差同时采样固定端点和自由端点计算边界残差。采样策略有一个细节如果每次都重新随机采样前一步刚优化好的区域下一步可能没采样到导致损失反复跳动。我试过两种方式一种是每次迭代全随机另一种是用一组固定的均匀内点加随机扰动。实测下来全随机的收敛更稳因为每步看到的都是不同的空间位置迫使网络在所有位置都保持物理合理性而不是在某几个固定点上“背答案”。训练采用 Adam 优化器初始学习率 1e-3。判断收敛不能只看总损失还要分开打印各分项损失我会在每 500 轮打印一次 four 项损失值特别关注 loss_fixed 是否持续下降——这通常是最难收敛的部分。# 定义模型和优化器 model PINN_Beam() optimizer torch.optim.Adam(model.parameters(), lr1e-3) # 梁的物理参数 L 3.0 # 梁长 x_fixed torch.tensor([[0.0]], dtypetorch.float32) x_free torch.tensor([[L]], dtypetorch.float32) # 坐标归一化到 [-1, 1]对 tanh 友好 x_fixed_norm (x_fixed - L/2) / (L/2) x_free_norm (x_free - L/2) / (L/2) for epoch in range(8000): optimizer.zero_grad() # 内部随机采样点归一化处理 x_pde_raw torch.rand(256, 1) * L x_pde (x_pde_raw - L/2) / (L/2) loss, loss_pde, loss_fixed, loss_free compute_loss( model, x_pde, x_fixed_norm, x_free_norm, L, EI, q ) loss.backward() optimizer.step() if (epoch 1) % 500 0: print(fepoch {epoch1:4d} | loss: {loss.item():.4e} | fpde: {loss_pde.item():.2e} | fixed: {loss_fixed.item():.2e} | ffree: {loss_free.item():.2e})这里有个关键细节输入坐标最好归一化到 [-1, 1] 再输入网络。原因有两层一是 tanh 在零附近梯度大、在两端饱和把采样点归一化到 [-1, 1] 可以让大部分点落在激活函数的敏感区间二是不同量级的物理量长度、挠度、荷载混在一起算高阶导时量纲差异会把梯度淹没。我一开始没归一化结果四阶导数值全在 1e-10 量级振荡PDE 残差根本学不动归一化之后才正常。3.4 解析解验证预测结果与理论值对比训练完成后在整个梁长上均匀取 200 个点将归一化坐标送入网络得到预测挠度再逆归一化还原成真实坐标对应的挠度值与解析解 w_exact q·x²·(6L² - 4Lx x²)/(24EI) 对比。# 验证阶段细网格采样 x_grid torch.linspace(0, L, 200).reshape(-1, 1) x_grid_norm (x_grid - L/2) / (L/2) with torch.no_grad(): w_pred model(x_grid_norm).reshape(-1).numpy() w_exact q / (24 * EI) * x_grid[:, 0].numpy()**2 * ( 6 * L**2 - 4 * L * x_grid[:, 0].numpy() x_grid[:, 0].numpy()**2 ) # 相对误差 rel_error np.linalg.norm(w_pred - w_exact) / np.linalg.norm(w_exact)第一次跑通时我的 RelL2 误差大约 2.3%最大挠度误差在 1% 以内这个结果已经足够用于工程预研。我又微调了一些参数后面会细说误差进一步降到低于 1% 的配置。4. 训练结果验证从损失曲线到挠度曲线的误差分析4.1 固定端与自由端的边界条件满足情况验证不仅要看整体挠度还要分开看物理量。我统计了训练后的边界残差w(0) 约为 3.2e-5 m比梁长三个量级小工程上可以视为零dw/dx(0) 约为 2.1e-4 rad对应的固定端转角基本为零自由端 d²w/dx²(L) 约为 5.3e-3对应弯矩约 5.3e3 N·m相比梁内最大弯矩 qL²/2 22500 N·m比值在 2% 左右自由端 d³w/dx³(L) 约为 0.4对应剪力约为 4e5 N 量级和荷载总量 qL 15000 N 相比还差两个数量级第四项自由端剪力残差是四项中收敛最慢的。原因也很好理解三阶导数已经是高阶量数值较小在损失函数里天然占比低。处理办法是把自由端残差单独放大权重比如给 loss_free 额外乘一个系数。我试过在损失函数里把自由端权重调到 1.5、固定端调到 1.0自由端剪力残差明显下降对整体挠度精度没有负面影响。4.2 整体挠度曲线与误差分布最大误差点在哪从整个梁长的挠度曲线看PINN 预测与解析解几乎重叠。误差最大的位置不在自由端反而在固定端附近约 15% 梁长处。这个现象最初让我困惑后来想明白了固定端附近挠度值很小相对误差天然偏大加上这个区域要同时满足“位移为零、转角为零”和“局部荷载产生弯曲”网络需要在极短距离内从零曲率过渡到正常弯曲状态几何上是约束最紧的区域。误差的第二大来源是边界点本身的离散采样误差。固定端只采样了一个点如果这个点在每次迭代中偶然被梯度更新带偏对邻近区域的影响会被四阶导数的非线性放大。我后续做了改进在 x0 附近额外加密采样 10 个点取平均作为固定端约束边界条件稳定了不少。这种做法相当于把“点约束”变成“小邻域约束”牺牲一点点精度换来训练稳定性很划算。4.3 训练收敛曲线为什么损失下降但误差不降训练过程中我看到一个现象前 500 轮loss_pde 急速下降但验证集的误差并没有同步下降。原因很直接网络先拟合的是一些简单模式比如“整体像个二次函数”这能让内部点的控制方程残差减小但边界条件还没满足挠度形态跟真实解还差得远。从第 1000 轮开始边界损失开始拉动整体形状loss_pde 进入缓慢下降平台这时验证误差才真正降低。这说明PINN 的训练不是“均匀优化”而是先满足某个子任务再逐步满足所有约束。理解这一点对调参很重要如果你看到 loss_pde 已经很低但结果不准千万不要以为再加迭代就行而是要去检查边界条件是否已经满足并把边界权重调大。5. 调参踩坑实录边界权重、激活函数与收敛性之间必须注意的事5.1 边界权重系数 λ_bc 的平衡点从1、10到100的对比我专门做了边界权重系数的对比实验同样训练 6000 轮结果如下λ_bc固定端残差自由端残差最大挠度相对误差内部PDE残差12.3e-28.2e-26.8%5.1e-2103.2e-52.1e-20.9%2.3e-3308.1e-61.4e-21.2%1.7e-31002.0e-69.5e-34.5%4.4e-2从表格能看到一个明显的矛盾λ_bc 太小边界条件不满足λ_bc 太大内部方程又学不好最终误差反而上升。我最终的优选是 λ_bc 10配合自由端单独放大综合效果最好。这个实验给我一个很重要的经验损失函数里各项权重本质上是一组“超参数”它们的合理值跟问题的具体量纲紧密相关。在换问题时不要盲目沿用别人网络里的权重值要分别打印各分项损失看它们的量级关系再决定要不要调整。5.2 激活函数对比tanh、Swish与ReLU的实测差异激活函数的选择直接影响高阶导数的表达能力。我做了一组对比tanh收敛稳定第 3000 轮损失降到 1e-4 以下误差低于 1%Swish前 500 轮下降很快但损失降到一定程度后开始振荡最终误差约 1.8%ReLU训练几乎没有进展四阶导全部为零验证误差高达 25% 以上ReLU 的失败在预期之中因为它分段线性高阶导数为零。Swish 虽然光滑且导数存在但在高阶导的数值稳定性上不如 tanh可能是其导数表达式中包含复杂的指数混合项在 BP 求导时容易放大数值误差。如果你用三角函数类激活比如 sin理论上也很适合解微分方程但没有实测过有兴趣可以自己跑一跑。5.3 从 Adam 切换到 L-BFGS 的收敛加速技巧标准 Adam 跑到 8000 轮大约能把误差压到 1% 以内但对工程师来说还是有点慢。我试了一种常见做法先用 Adam 训练 3000 轮让网络进入正确的“形态”然后切换到 L-BFGS 优化器做微调。L-BFGS 是二阶方法在 PINN 社区里经常用于最终精修它利用损失函数的曲率信息收敛速度比 Adam 快得多。切换代码很简单optimizer_lbfgs torch.optim.LBFGS(model.parameters(), lr1.0, max_iter500) def closure(): optimizer_lbfgs.zero_grad() loss, *_ compute_loss(model, x_pde, x_fixed_norm, x_free_norm, L, EI, q) loss.backward() return loss optimizer_lbfgs.step(closure)注意 L-BFGS 不支持小批量随机采样因为它需要在整个计算图上做多次线搜索训练数据必须是固定的。所以切换 L-BFGS 前要把采样点固定下来生成一组覆盖整个梁长和边界邻域的常驻训练点。我在 Adam 阶段用全随机切到 L-BFGS 后改用固定点集最终最大挠度误差从 0.9% 降到了 0.3%而且只额外花了不到一分钟。这个“Adam 粗训练 L-BFGS 精修”的组合拳是我这段时间最推荐的做法。5.4 训练发散时的急救措施学习率衰减与梯度裁剪即使按上面配置来偶尔也会遇到训练发散。最常见的表现是 loss 从某个 epoch 开始爆炸数值变成 nan 或者 inf。原因往往是四阶导数的计算图上累积了过大的梯度一次更新把参数推到 tanh 的饱和区。我的急救措施有两步。第一给 Adam 加学习率衰减比如每 2000 轮衰减 0.5让网络在后期只做小范围调整避免步伐过大冲进坏区域第二给总损失施加梯度裁剪torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0)把梯度范数限制在 1.0 以内。这两步一加我从没见过训练再发散了。另外有一个小技巧如果 loss 开始飘先不急着改代码看一眼是哪个分项在涨。如果 loss_fixed 涨说明边界梯度把参数推得太远如果 loss_pde 涨说明内部采样点可能碰到了某些高曲率区域。分项定位比盲目调全局学习率高效得多。6. 延伸思考PINN在结构分析里还能怎么用6.1 变截面梁与非均布荷载改方程不改网格回到我最初的项目需求。传统方法每换一次荷载分布就要重新推解析解或者重新剖分网格。PINN 的代码里非均布荷载只是把损失函数中常数 q 换成 q(x) 的函数表达式变截面梁只是把 EI 从常数量改成 EI(x)。网络结构、损失函数框架、训练流程全部不用动。我顺手验证过一个线性变化荷载 q(x) q0·(1 x/L) 的算例PINN 训练出来的挠度曲线和数值积分参考解保持了良好的一致性。主要误差出现在自由端附近这跟在均布荷载算例中的表现一致说明误差模式不是特定荷载造成的而是边界处理本身的难点。如果你要做变荷载系列试算PINN 可以显著减少重复工作。6.2 测量数据驱动的“逆问题”从挠度反推荷载分布PINN 更吸引我的能力是解决逆问题。实际工程中测量某几个位置的挠度容易但要知道梁上的荷载分布很难。如果构建一个网络输入坐标输出挠度损失函数里除了控制方程、边界条件再加入“已知测点的挠度值等于测量值”这组约束训练收敛后就能从挠度数据中反推出满足物理规律的完整挠度场进而通过边界挠度的导数估算荷载。这个方向在结构健康监测里很有价值。比如某结构仿真项目中我们在梁的关键位置布置少量位移传感器想反演作用荷载的分布形式用传统优化方法需要反复调用正问题求解器计算成本很高PINN 的正则化作用可以让解始终落在物理可行的空间里避免“用很多未知数去拟合几个点”的过拟合陷阱。6.3 从一维梁到二维板的扩展思路一维梁问题相对简单PINN 处理起来游刃有余。扩展到二维板壳问题控制方程从常微分方程变成偏微分方程比如克希霍夫板方程含有双调和算子网络输入层增加一个 y 坐标损失函数里的导数换成对 x、y 的混合偏导。原理上没有本质区别但有两个实际问题会冒出来一是二维高阶偏导的计算量和显存占用会成倍增加训练时间从分钟级变成小时级二是二维问题的边界更加复杂边界条件数量更多权重的调参难度直线上升。我暂时还没跑通完整的二维板 PINN 算例但已经在一维梁上验证了所有关键机制。如果你的工作涉及二维问题建议先拿一维梁把“高阶导的计算流程”和“边界权重的手感”练熟再推二维会少走很多弯路。最后说一点个人体会。PINN 不是一个“比有限元更快的求解器”它是一个“把物理规则嵌入神经网络的建模框架”。在单根梁这种简单问题上它比传统方法慢但它的价值在于当你面对荷载形式复杂、边界条件奇怪、甚至只有零散测量数据的问题时它能给你一条不用重新发明算法的出路。跑通这个悬臂梁例子的过程是我把 PINN 从概念变成工具的最重要一步。希望这篇文章能帮你少踩几个我踩过的坑把时间和精力留在真正有意思的结构问题上。