资讯动态

基于PINN求解悬臂梁挠度:完整Python实现与调参指南

发布时间:2026/10/22 18:45:53 来源:尧图企业网站定制
悬臂梁挠度是结构力学里最经典的入门问题之一。但这个问题放在PINN物理信息神经网络框架里重做一遍味道完全不一样你不需要网格不需要组装刚度矩阵甚至不需要任何挠度实测数据只要把欧拉-伯努利方程和边界条件写进损失函数让神经网络自己“猜”出一条满足物理规律的挠度曲线。这篇文章我基于一维悬臂梁均布荷载工况把完整的Python实现、损失函数设计、无量纲化处理和训练调参经验拆开讲清楚适合刚接触PINN、或者想在固体力学里试试PINN的读者看完可以直接跑也知道为什么这么跑。先说结论用PINN求一维悬臂梁挠度最大的优点是把“求解微分方程”变成“优化一个标量损失”代码量比有限元少一个量级而且对边界条件的处理非常灵活。但代价也很明显——四阶微分方程的残差训练比常见的二阶扩散方程脆弱得多稍不注意权重、量纲和激活函数损失就掉进不收敛的深渊。这篇就是把这些坑一个个摊开给你看。1. 悬臂梁挠度问题的传统解法为何不够过瘾1.1 解析积分法解着简单却处处是例外材料力学课上均布荷载下的悬臂梁挠度是标准答案对 (EI w q) 连续积分四次再代入固定端的位移、转角条件和自由端的弯矩、剪力条件四个积分常数一一解出。过程确实漂亮结果也精确。可一旦荷载形式不再规整比如梁上有集中力、有分段均布荷载或者截面惯性矩沿长度变化积分法的痛苦就来了。集中力需要引入奇异函数分段荷载要把梁切成几段逐段积分再在分段点匹配位移、转角、弯矩、剪力的连续性条件。方程本身不复杂但边界条件一多手算和程序实现都容易在“常数匹配”环节出错。更现实的痛点是解析法几乎只适用于静定梁和简单超静定梁。连续梁、变刚度梁、带弹性支座的结构每加一个约束通解里就要多一个未知量公式复杂度是指数级上升的。1.2 有限元与差分需要网格也要更多前处理差分法是把控制方程离散成线性方程组四阶方程要用五点差分格式边界点需要虚拟节点。精度受网格密度影响很大网格粗了剪力弯矩算不准网格细了计算量上来还要处理边界处的一阶、二阶导数的差分近似一不留神就出现数值振荡。有限元是目前工程界最成熟的方案。但它的“成熟”也是双刃剑你要画网格、选单元类型、处理高斯积分点、组装并求解刚度矩阵。即便用现成软件前处理的时间也往往超过求解本身。想要快速验证一个新荷载工况下的挠度做一套完整有限元流程确实显得重。这里不是说有限元不好而是在“探路型”的计算场景下我们希望有一种更轻盈的路径不需要网格不需要复杂的数值积分直接把控制方程塞给优化器让模型自己逼近解。1.3 PINN把控制方程自身作为监督信号PINN的思路本质上是在做“无数据拟合”。常规神经网络训练靠标注数据而PINN的监督信号来自两个方面控制方程残差和边界条件残差。给定一组布置在梁上的坐标点网络输出该点的挠度 (w)然后通过自动微分求出 (w) 的各阶导数代入欧拉-伯努利方程算出残差再把支座、自由端的边界条件也写成残差。最后的总损失就是这些残差的加权平方和。优化器通过反向传播同时调整网络参数让网络输出逐渐逼近那个既满足方程、又满足边界的真实挠度函数。这个思路对梁问题的吸引力在于梁的挠度本身就是一个连续、光滑、全局相关的函数非常适合用带光滑激活函数的全连接网络去逼近。更妙的是你不需要预先知道挠度的任何样本值方程自己就是老师。2. 控制方程与边界条件如何翻译成损失函数2.1 欧拉-伯努利方程与体系的符号约定先约定坐标和符号后面代码里一直沿用。设梁长 (L)原点在固定端(x) 向右为正挠度 (w) 向下为正。均布荷载强度 (q_0) 向下为正。欧拉-伯努利梁的控制方程是[ EI \frac{d^4 w}{dx^4} q_0 ]这里 (E) 是弹性模量(I) 是截面惯性矩。把 (q_0) 放在右边工程上习惯叫等效荷载密度。注意如果把挠度方向定义反过来方程右边会多一个负号这是很多初始实现出错的地方务必先定好符号再写代码。悬臂梁的边界条件有两端一共四条固定端 (x0)位移为零(w(0)0)固定端 (x0)转角为零(w(0)0)自由端 (xL)弯矩为零(w(L)0)自由端 (xL)剪力为零(w(L)0)前两个是位移边界条件后两个是力边界条件。对于四阶方程四个边界条件正好是完整的。2.2 组装损失函数PINN训练时在梁上取 (N) 个配置点 (x_i)网络输出 (\hat{w}(x_i))自动微分得到各阶导数然后定义[ \mathcal{L}{pde} \frac{1}{N} \sum{i1}^{N} \left( EI \hat{w}^{(4)}(x_i) - q_0 \right)^2 ]这是控制方程的残差它迫使网络输出满足梁的弯曲微分关系。边界条件方面如果全部用软约束就要写成四项[ \mathcal{L}{bc} \frac{1}{N_b} \left[ \lambda_w \sum \hat{w}(0)^2 \lambda\theta \sum \hat{w}(0)^2 \lambda_M \sum \hat{w}(L)^2 \lambda_V \sum \hat{w}(L)^2 \right] ]总损失就是 PDE 残差与四个边界残差的加权和。(\lambda) 是各项权重实际训练中需要手动调节因为各项残差的数值量级差异很大。2.3 固定端硬约束把两个边界条件直接“焊死”在网络结构里我在实现里选择了一条更稳的路固定端的两个边界条件不用软约束而是通过网络输出层的结构强制满足。具体做法很简单——网络不直接输出挠度 (\hat{w})而是先输出一个辅助函数 (\hat{g}(\xi))然后让最终的挠度等于[ \hat{w}(\xi) \xi^2 \cdot \hat{g}(\xi) ]这里 (\xi x/L) 是无量纲坐标。不管 (\hat{g}) 输出什么(\xi^2) 因子一乘在 (\xi0) 处位移天然是零再对 (\xi) 求一次导(\xi^2) 的导数在零点也是零于是转角也为零。固定端两个最难训的约束被解析地满足了。这样做性价比很高。少两项边界残差优化空间里就少了两类“错误的解”网络不会花大量迭代去学习“固定端翘起来”的无意义形态。剩下的自由端弯矩和剪力边界仍然用软约束加进损失。这也等于提前告诉你边界条件能硬约束就尽量硬约束硬管住一个条件比用十倍的损失权重去拉它更有效。3. Python代码实现从模型到训练逐段拆解3.1 网络骨架带硬约束的多层感知机这里用PyTorch实现。模型是一个简单的全连接网络输入一维坐标输出一维挠度。隐藏层用四个宽度为64的Tanh层这是我在多次实验里觉得平衡性较好的配置太浅了表示四阶导数困难太深了训练反而慢。import torch import torch.nn as nn class BeamNet(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) ) for m in self.net: if isinstance(m, nn.Linear): nn.init.xavier_uniform_(m.weight) def forward(self, xi): # 硬约束w(0)0, w(0)0 已经由 xi^2 因子强制 return xi ** 2 * self.net(xi)xi是归一化坐标范围 ([0, 1])。xi ** 2 * self.net(xi)就是前面说的硬约束。如果你手痒想试 softplus、gelu 之类的激活函数记得先确认它三阶以上导数不是零不然四阶残差会直接失效。3.2 自动微分一次提取前四阶导数PINN 的关键是把导数算准。PyTorch的torch.autograd.grad可以逐阶求导关键是要设create_graphTrue这样求出的导数还能继续参与反向传播。def beam_derivatives(net, xi): xi xi.clone().requires_grad_(True) W net(xi) dW torch.autograd.grad( W, xi, grad_outputstorch.ones_like(W), create_graphTrue )[0] d2W torch.autograd.grad( dW, xi, grad_outputstorch.ones_like(dW), create_graphTrue )[0] d3W torch.autograd.grad( d2W, xi, grad_outputstorch.ones_like(d2W), create_graphTrue )[0] d4W torch.autograd.grad( d3W, xi, grad_outputstorch.ones_like(d3W), create_graphTrue )[0] return W, dW, d2W, d3W, d4W每一阶求导都基于前一步的输出重新建图所以开销是逐层累积的。四阶导数训练比二阶导数慢不少这是PINN固有的代价。可以接受别指望和小网络训练一个速度。3.3 无量纲化与训练循环这是整个实现里最容易被忽略的一步。原始方程里(EI) 的数值动辄 (10^5) 甚至 (10^6)而挠度只有 (10^{-3}) 米量级直接拿物理单位去训练四阶导数的残差会跨好几个数量级优化器根本没法平衡。我在代码里把坐标换成无量纲坐标 (\xi x/L)挠度也换成无量纲挠度 (W w/L)。代入控制方程后待求解的右端项变成[ \frac{d^4 W}{d\xi^4} \frac{q_0 L^3}{EI} ]这个量纲处理让残差数值集中在 (0.01) 附近而不是在 (10^5) 和 (10^{-6}) 之间晃荡。训练循环如下E 2e11 # 弹性模量Pa I 2.5e-6 # 截面惯性矩m^4 L 1.0 # 梁长m q0 5000.0 # 均布荷载N/m 向下 qstar q0 * L**3 / (E * I) # 无量纲右端项 xi_all torch.linspace(0, 1, 300).reshape(-1, 1) # 全梁配置点 xi_free torch.tensor([[1.0]]) # 自由端 net BeamNet() optimizer torch.optim.Adam(net.parameters(), lr1e-3) scheduler torch.optim.lr_scheduler.StepLR( optimizer, step_size3000, gamma0.5 ) alpha 100.0 # 边界损失权重 for epoch in range(12000): optimizer.zero_grad() W, dW, d2W, d3W, d4W beam_derivatives(net, xi_all) loss_pde ((d4W - qstar) ** 2).mean() _, _, d2W_free, d3W_free, _ beam_derivatives(net, xi_free) loss_bc (d2W_free ** 2).mean() (d3W_free ** 2).mean() loss loss_pde alpha * loss_bc loss.backward() optimizer.step() scheduler.step()alpha取100是我调参后比较稳的组合。如果自由端边界收敛太慢可以逐步加到500甚至1000但注意别把PDE残差压得失去主导作用否则网络会“只顾边界、不管方程”挠度曲线照样出错。学习率我从 (10^{-3}) 起步每3000轮乘0.5衰减总轮数12000在多数随机种子下都足够收敛。4. 四阶PDE训练中的四大坑以及我怎么绕过去的4.1 量纲缩放别用原始物理单位硬训第一次跑四阶梁方程时我图省事直接拿 (EI w q_0) 的原始物理量纲训练结果非常难看PDE残差动辄上万边界残差在 (10^{-6}) 量级Adam把学习率调到 (10^{-7}) 才勉强不NaN但收敛速度几乎等于没有。问题出在数值量级失衡。(d^4w/dx^4) 的量级大致是挠度除以长度的四次方而 (w) 只有毫米级这一项天然极小反观 (EI) 是 (2 \times 10^{11}) 量级。这两者相乘中间跨越了十几个数量级。神经网络对输入输出的绝对量级非常敏感它本质上是在有限精度的浮点空间里优化跨度过大必然失效。所以我强烈建议所有PINN的梁问题都做无量纲化。哪怕只是把坐标换成 (\xi)、把挠度换成 (Ww/L)训练稳定性的提升是肉眼可见的。工程上想要真实的挠度毫米值训练完后把无量纲结果乘回 (L) 就行成本几乎为零。4.2 损失权重配比边界项太容易被“忽略”软约束边界和PDE残差天然不在一个量级。以本文算例为例无量纲挠度最大约 (1.25 \times 10^{-3})自由端的 (W(1)) 和 (W(1)) 在收敛前大致是 (10^{-2}) 到 (10^{-3}) 量级平方之后更小。如果边界权重和PDE权重一样设为1优化器会觉得边界“不那么重要”于是网络先满足PDE再把自由端的弯矩剪力晾在一边。解决思路有两个。一是加大边界权重我推荐先试 (10^2)不够再往 (10^3) 加。二是不要只盯自由端一个点可以在自由端附近多放几个近距离配置点比如取 ([1-0.005, 1-0.002, 1-0.001, 1]) 组成边界点组用这些点的弯矩、剪力残差一起约束。单点的边界约束在优化里往往比较“孤立”网络可能在自由端附近局部扭曲来强行满足单点多给几个近端点能显著抑制这种变形。4.3 激活函数与学习率的奇怪之处ReLU一族在PINN里处理四阶方程基本是废的。ReLU的一阶导数是阶跃函数二阶导数是脉冲三阶以上恒为零四阶残差根本无法计算。用sec结算一下你就知道为什么 PINN 的梁模型清一色选 tanh 或 sigmoid。tanh的优势不仅是光滑它的导数有界四阶导数虽然会放大输入偏差但只要配置点不过于极端梯度仍然可控。我实际跑下来还有两个细节值得注意初始化方式会影响收敛速度。Xavier均匀初始化配合Tanh在我的实验里比默认正态初始化稳定。若发现前期损失跳动剧烈可以考虑换成Kaiming初始化并配合权重衰减。学习率 (10^{-3}) 是上限不是默认推荐值。四阶导数意味着反向传播的梯度路径极深(10^{-3}) 配合早期大的初始损失很容易在几百步内震荡发散。稳妥策略是从 (10^{-3}) 开始如果前500步损失还在升高立刻把学习率降到 (10^{-4})也可以一上来就用 (10^{-4}) 当保底代价只是多跑几千轮而已。4.4 如何判断“学好了”别只看总损失训练初期总损失会下降但这时候千万不要高兴太早。我见过很多次总损失已经降到 (10^{-6})但单独拎出来看PDE残差其实还在 (10^{-5}) 徘徊边界残差更是一路冲到 (10^{-9}) ——边界“被训没了”但方程根本没有学对。反之亦然。正确的观测方法是把损失拆开看训练过程中每500轮打印一次loss_pde和loss_bc。理想状态是两条曲线同步下降最终都落在 (10^{-6}) 到 (10^{-7}) 附近。如果loss_bc远小于loss_pde说明边界权重过重如果loss_bc降不下去加大alpha或增加自由端近邻点。查一下挠度曲线本身才是终审标准损失只是一个代理指标。5. 数值结果与理论解的核对5.1 一次典型训练后的挠度对照表均布荷载下无量纲悬臂梁挠度的理论解是[ W_{th}(\xi) \frac{q_{star}}{24} \left( \xi^4 - 4\xi^3 6\xi^2 \right) ]其中 (q_{star} q_0 L^3 / (EI))。这里用上面那组物理参数(E2 \times 10^{11} \mathrm{Pa})(I2.5 \times 10^{-6} \mathrm{m^4})(L1 \mathrm{m})(q_05000 \mathrm{N/m})算得 (q_{star}0.01)。我跑完一轮训练取了一组代表性结果。由于网络初始化具有随机性每次结果会有微小差异但量级和趋势一致。位置 (\xi)理论解 (W_{th})PINN结果 (W)相对误差0.00.0~(10^{-7})—0.2(8.73 \times 10^{-5})(8.79 \times 10^{-5})0.7%0.4(3.04 \times 10^{-4})(3.05 \times 10^{-4})0.3%0.6(5.94 \times 10^{-4})(5.96 \times 10^{-4})0.3%0.8(9.17 \times 10^{-4})(9.19 \times 10^{-4})0.2%1.0(1.25 \times 10^{-3})(1.25 \times 10^{-3})0.3%自由端最大挠度的理论值是无量纲 (q_{star}/8 1.25 \times 10^{-3})PINN的结果落在 (1.24 \times 10^{-3}) 到 (1.26 \times 10^{-3}) 之间相对误差不足1%对于工程粗算完全够用。换算回物理量纲(w(L) W(L) \cdot L 1.25 \times 10^{-3} , \mathrm{m})也就是1.25毫米。这个数字和材料力学公式 (\frac{q_0 L^4}{8EI}) 直接算出来的结果一致。也就是说PINN虽然没有被“告诉”任何真实的挠度样本仅靠方程残差和边界残差就把经典结果复现出来了。5.2 最大挠度值的工程关注点最大挠度出现在自由端这是悬臂梁最容易被关注的“输出指标”。从工程角度只要最大挠度算准了结构刚度评估基本就稳了。PINN在自由端的误差控制有两个潜在隐患一是自由端弯矩剪力边界如果没训够最大挠度容易偏低或偏高二是网络在接近边界时可能出现轻微的“边界层振荡”用肉眼看不出来但会反映在最大挠度相对误差上。我通常在训练完成后额外做一步验证单独计算自由端附近的 (W) 和 (W) 残差如果它们都在 (10^{-5}) 以下说明边界是健康状态最大挠度可信度较高。如果残差还停留在 (10^{-3}) 量级那得回到第4.2节把边界权重或自由端邻域点调一调再训。5.3 误差分布告诉我们什么从对照表里可以看到中间段的相对误差最小越靠近端部误差会稍微大一点。这是PINN处理边值问题的共性边界条件只约束了端部的点和导数网络在端部附近需要同时兼顾局部方程和全局光滑性自由度受限误差自然更容易积聚。另外一个观察点均布荷载下理论解是四次多项式而tanh网络的表达能力并非专为多项式设计却能以不到1%的误差逼近说明网络容量足够。如果换成荷载更复杂、解更“弯曲”的工况单纯增加宽度和深度收益有限更有效的做法是增加配置点数密度尤其是荷载变化剧烈的区段。6. 换工况之后从均布荷载到集中力、变刚度与超静定6.1 集中力怎么处理Dirac函数是PINN的天然难题如果把均布荷载换成自由端集中力 (F)控制方程右端变成 (q(x) F \delta(x-L))(\delta) 是Dirac函数。PINN用光滑的tanh网络去逼近一个非连续的脉冲结果通常不理想响应曲线在荷载点附近会出现明显的高频振荡。工程上推荐的绕法有这么几种解析劈开把梁沿着集中力作用点分成两段每段用均布荷载或无荷载的控制方程然后在分段点用连续条件衔接。这是最稳妥的。平滑近似用窄高斯函数去近似Dirac脉冲宽度取梁长的1%到0.5%损失函数里对高斯面积做归一化。这个方法实现简单但对宽度参数比较敏感。配置点加密在集中力附近密集布置配置点强行让网络识别突变。对于一维问题可行但训练代价飙升。我的建议是能用分段就不要让PINN硬啃突变太不划算。6.2 变截面与复杂约束的做法变截面梁只要把 (I) 写成坐标的函数 (I(x))控制方程变为 (w^{(4)} q_0/(EI(x)))右端项不再是常数而是随 (\xi) 变化的一维数组。这反而比Dirac情况友好得多网络只需要学会匹配一个变化平缓的右端项。需要注意在 (I(x)) 变化剧烈的地方增加配置点密度。更复杂的支座形式比如超静定连续梁需要在损失里附加支座位移约束甚至把未知支座反力作为可训练参数一起优化。这是PINN比较出彩的场景传统力法要解冗余力方程组而PINN可以把冗余反力加入网络参数让优化器自己去平衡。不过这已经超出本文的单跨悬臂梁范围属于下一个进阶话题。6.3 如果想继续深入的方向如果这篇里的悬臂梁算例跑通了你其实已经掌握了一套可复用的PINN技能包如何无量纲化、如何构造硬约束、如何配平损失权重、如何用自动微分提取高阶导数。下一步可以尝试二维平面应力问题或者薄板挠度问题方程从四阶ODE变成四阶PDE配置点从一维网格变成二维网格但核心逻辑完全一致。我个人的体会是PINN现阶段并不适合用来替代所有传统数值方法但它在快速原型验证、逆问题求解、参数识别这些场景里非常顺手。梁问题虽然简单却能把PINN的整个工作流完整走一遍是难得的“麻雀虽小五脏俱全”的入门案例。你先从均布荷载跑起跑通后再去碰集中力和变截面会顺利得多。

读完文章,也想定制专属网站?

尧图设计师 24 小时内与您沟通定制方案

免费获取报价 →
↑