☰
物理信息神经网络PINN求解悬臂梁挠度:Python实现与降阶技巧
2026/10/22 9:04:19 网站建设 项目流程

去年我在一个小型结构分析项目里反复计算一维悬臂梁在不同均布载荷下的挠度曲线。用传统有限元流程当然没问题,但每次调整载荷工况就得重新建模、重新划分网格,哪怕梁只有一米长,这些重复劳动依然很烦。后来我把目光转向了PINN物理信息神经网络——它不需要网格,只需要把梁的控制方程直接写进神经网络的损失函数,让网络自己去逼近满足方程和边界条件的解。这篇文章完整记录我用Python实现PINN求解一维悬臂梁挠度问题(也就是标题里那个经典问题)的全过程,包括控制方程推导、损失函数设计、降阶技巧、代码拆解和训练中踩过的坑。如果你刚接触PINN,想用最小成本跑通一个力学求解案例,这个悬臂梁问题是非常合适的入门原型。

1. 悬臂梁挠度问题的物理方程与边界条件

1.1 从欧拉-伯努利梁方程到四阶常微分方程

悬臂梁是结构力学里最经典的模型之一:一端固定、一端自由,梁的长度方向取为x轴,挠度记为w(x)。在材料力学的小变形、线弹性假设下,细长梁的挠度控制方程是欧拉-伯努利梁方程:

EI * d⁴w/dx⁴ = q(x)

其中E是弹性模量,I是截面惯性矩,q(x)是单位长度上的横向分布载荷。这个方程只含x一个自变量,是一维四阶常微分方程。

为什么要纠结成一个四阶方程?因为梁的内力和挠度之间存在明确的导数链关系。弯矩M = -EI * d²w/dx²,剪力V = dM/dx = -EI * d³w/dx³,而分布载荷q = dV/dx = -EI * d⁴w/dx⁴。也就是说,载荷q是挠度w的四阶导数,方程必须写成四阶才能完整描述"载荷-剪力-弯矩-挠度"之间的传递关系。

这里必须先说清楚符号约定,否则后面代码和解析解一对上就容易乱。我取x轴向右为正,w向上为正,载荷q向上为正。悬臂梁受到的重力均布载荷是向下的,取q = -q0,其中q0是载荷大小(正值)。在这个约定下,自由端向下挠,所以w的解析解是负值。为了简化后面的分析和训练,直接取无量纲参数:梁长L=1,抗弯刚度EI=1,均布载荷q0=1,于是控制方程变成了:

EI * d⁴w/dx⁴ = -1

这是一个线性的四阶常微分方程,对神经网络来说,待逼近的解曲面非常光滑,很适合作为PINN的第一个力学算例。

1.2 四个边界条件:两端各自的物理约束

四阶方程需要四个边界条件才能得到唯一解。悬臂梁的两端给出的边界条件类型完全不同:

固定端x=0处,梁被完全约束,挠度为零,转角也为零。转角在小变形假设下就是挠度的一阶导数θ≈dw/dx,所以固定端给出的是两个强制几何边界条件:

  • w(0) = 0
  • w'(0) = 0

自由端x=L处没有任何支撑,弯矩必须为零,剪力也必须为零。用挠度表达就是二阶导和三阶导为零:

  • w''(L) = 0
  • w'''(L) = 0

这四个边界条件加上四阶控制方程,构成了完整的边值问题。传统做法是直接找解析解或者用数值方法离散求解。均布载荷下的解析解很漂亮:

w(x) = -q0 * x² * (6L² - 4Lx + x²) / (24EI)

把L=1、EI=1、q0=1代进去,得:

w(x) = -x² * (6 - 4x + x²) / 24

这个解析解后面就是验证神经网络预测精度的一把尺子。自由端x=L处的最大挠度是w(1) = -1/8 = -0.125,量级很小,这也是神经网络输出层需要匹配的目标。

1.3 为什么先做无量纲化

直接拿真实工程参数训练会让神经网络非常痛苦。想象一下,真实钢梁的EI可能是2×10¹¹ N·mm²量级,均布载荷是几千N/m,挠度是毫米甚至零点几毫米量级。控制方程左边是四阶导数乘以巨大系数,右边是几千的量级,网络输出的挠度却是10⁻³量级,这种量级失衡会让损失函数里的各项完全不在一个数量级上,梯度更新非常不稳定。

无量纲化处理之后,所有变量都落在0.1到1附近,方程系数都是1,损失函数各项量级一致,训练难度大幅下降。在PINN里,这一步不是可选项,几乎算强制操作。实际工程问题当然不能直接把物理参数丢掉,而是在训练完成后把无量纲结果乘回去,或者把方程写成带参数的形式,让EI和q作为变量参与归一化。一维悬臂梁这个入门案例,直接把L=EI=q0=1处理最省心。

2. 为什么选PINN而不是有限元:一维问题里的取舍

2.1 传统有限元的常规流程

用有限元解一维梁挠度,流程非常成熟:建立几何模型、赋材料属性、划分网格、施加边界条件和载荷、求解线性方程组,后处理提取节点位移。如果用的是通用有限元软件,一个悬臂梁模型从建模到出结果,熟练的话十几分钟就搞定,结果精度还很高。对于这种单次求解的场景,我没有理由不用有限元。

但问题往往出在"重复"上。如果我想分析十种不同的载荷分布,或者梁上某一段截面惯性矩发生变化,传统流程几乎要从建模那步重来一遍,至少网格要重新划分。网格剖分对几何变化和边界变化是敏感的,这是有限元的固有特点。PINN在这个场景下的优势就体现出来了:它根本不关心网格,只需要在求解域里随机撒点。

2.2 PINN在力学问题上到底强在哪

PINN的核心思路是把物理方程作为约束加入到神经网络训练中,让网络输出逼近满足方程和边界条件的函数。和传统有限元比,它最吸引人的是三点:

第一,无网格。训练时随机采样内点,边界是另外采的固定点,几何变化只需要改变采样域,不需要重新生成网格。第二,自动微分。网络输出的导数是用自动微分算的,不需要手动推导和组装刚度矩阵,对不熟悉有限元编程的人来说门槛更低。第三,扩展反问题很方便。如果我想从实测挠度反推梁的EI值,只需要把EI设为可训练参数代入损失函数一起优化,这种操作在传统有限元里实现起来要麻烦得多。

当然,PINN也有明显的短板:训练时间长、超参数敏感、精度通常达不到高精度有限元的水平。对于一维梁这种小规模问题,这些短板还不算致命,甚至训练时间往往不到一分钟,这是很好的试验场。真正要衡量的是:它能不能在放弃网格的前提下,拿到足够接近解析解的挠度曲线。我做完这个案例的结论是:可以,但前提是损失函数设计要对。

2.3 四阶方程在PINN里的直接麻烦

如果保持原方程EI * w'''' = q不变,网络只需要一个输出w(x),损失函数写成方程残差加边界残差,看似简单。但实际训练时问题非常突出。

首先,自动微分需要嵌套四层才能得到d⁴w/dx⁴,计算图很深,反向传播时高阶导数的梯度会变得非常不稳定,经常出现梯度爆炸,loss直接变成NaN。其次,高阶导数对网络输出的微小振荡极其敏感,哪怕是tanh这种平滑激活函数,四阶导数的数值波动也会被放大得厉害。再次,激活函数选择受限,ReLU的一阶导是阶跃、二阶导是零,根本没法用于二阶以上问题,直接用ReLU训练只会得到一个loss不下降的模型。

这不是说四阶直接算完全不可行,我自己试过,小网络加极小学习率加梯度裁剪,勉强能出结果,但整个过程非常脆弱,换一个随机种子结果就飘。经过对比,我发现最好的解决办法是降阶,把四阶系统拆成二阶系统,这也是下一节的重点。

3. 损失函数设计与降阶技巧:从四阶硬算到二阶联立

3.1 先看直接四阶的损失函数长什么样

如果坚持用单输出w(x),损失函数由两部分构成。内点方程损失是在求解域内部随机采样N_f个点,计算每个点上的控制方程残差:

L_eq = mean((EI * d⁴w/dx⁴ + 1)²)

这里q=-1,所以写成+1。边界条件损失是在x=0和x=L两个端点计算强制边界和自然边界的残差:

L_bc = w(0)² + (w'(0))² + (w''(L))² + (w'''(L))²

总损失是L_eq加上权重系数λ_bc乘以边界损失。理论上这个设计完全符合PINN的框架,实际跑起来会发现内点方程残差下降得很快,但边界条件总是慢半拍,整个曲线形状对但端点飘得厉害。原因在于四阶导数的残差梯度不稳定,以及边界信息在深度计算图里传播效率低。

3.2 降阶法:引入弯矩,把一个四阶变成两个二阶

既然四阶导数太脆弱,那就从物理关系上把阶数降下来。在梁理论里,弯矩和挠度满足M = -EI * w'',同时弯矩的二阶导与载荷相关:M'' = -q。这两个关系都是二阶的,组合起来正好等价于原来的四阶方程。

所以我不再让网络只输出w一个量,而是让它同时输出w和M两个量。控制方程拆成两个二阶约束:

  • EI * w'' + M = 0
  • M'' + q = 0(在q=-1时写作M'' + 1 = 0?注意符号,见后文)

这样网络只需要计算二阶导数,计算图深度减半,训练稳定性明显提升。边界条件也需要转换到这两个变量上。固定端依然是w=0和w'=0,自由端则换成M=0和M'=0,因为自由端的弯矩、剪力为零,M'就是剪力(符号约定下)。

需要仔细处理符号。前面约定q=-1(向下载荷),控制方程EI * w'''' = -1。由M = -EI * w'',则M'' = -EI * w'''' = -q = 1。所以降阶后两个方程是:

  • EI * w'' + M = 0(因为M = -EI * w'',移项就是EI * w'' + M = 0)
  • M'' - 1 = 0(因为M'' = 1)

这里很容易写成M'' + q = 0,如果q=-1,就是M'' -1 =0,和M''=1等价。我在代码里为了直观,写的是d2M - 1 = 0。有的版本喜欢写M'' + q0 = 0,此时q0是向下载荷大小1。两种写法差个符号,执行时只要和解析解一致就行。我最终用的方案是:

  • loss_w = mean((EI * d2w + M)²)
  • loss_M = mean((d2M - 1)²)

同时,边界条件:

  • 固定端:w(0)=0,w'(0)=0
  • 自由端:M(L)=0,M'(L)=0

四个边界条件依然正好四个。方程变成两个二阶约束后,网络输出层变成两个节点,一个对应w,一个对应M。M的输出可以直接用来画弯矩图,这是一个额外收获。

3.3 边界损失权重为什么必须单独乘系数

PINN的损失函数通常是方程残差和边界残差的加权和。如果直接相加,内点方程残差会占绝对主导,因为边界点只有两个,内点有上千个,边界条件在总损失里的占比非常低,网络会倾向于优先满足方程而牺牲边界。结果就是曲线形态对,但端点不贴合,固定端转角可能明显不为零。

解决方法是给边界损失乘一个权重系数λ_bc。经验上,一维问题取10到100之间都合理,我习惯取30到50。这个系数不是越大越好,太大会导致网络只顾边界而忽略了内点方程,反而整体误差上升。粗调时先看总loss是否稳定下降,再微调λ_bc让两端预测误差最小。这种方法本质上是多目标优化的加权处理,理解成"把边界条件的优先级调高"就好。

3.4 内点采样和边界采样的策略差别

内点采样使用均匀随机分布即可,一维问题直接在(0,L)区间上生成随机数。每个训练步可以重新采样一批新的内点,这样每一轮迭代都相当于重新配点,大幅减少网络对固定采样点的过拟合。边界点则固定,因为边界位置就两个,为了让边界损失的梯度更稳定,我通常把固定端和自由端的点重复若干份,比如各重复100次,相当于批量边界配置。边界重复份数太小的话,两个点的梯度噪声太大,λ_bc的调节也会变得不灵敏。

4. Python代码拆解:网络结构、自动微分与训练循环

4.1 网络结构与激活函数选择

网络主体用全连接神经网络(MLP),输入是x坐标(形状N×1),输出是w和M两个标量(形状N×2)。隐藏层用4层,每层64个神经元,激活函数全部用tanh。

选tanh而不是ReLU,是因为我们至少要求二阶导,ReLU的二阶导恒为零,根本没法用。tanh处处非零且可导,高阶导连续,是PINN常规问题里最稳的激活函数之一。SiLU或者GELU也可以,但我实际对比下来tanh的收敛速度最直接。初始化用Xavier均匀初始化,输出层不加激活函数,保持线性回归输出。

import torch import torch.nn as nn import numpy as np class PINN(nn.Module): def __init__(self, n_hidden=4, n_neurons=64): super().__init__() layers = [nn.Linear(1, n_neurons), nn.Tanh()] for _ in range(n_hidden - 1): layers += [nn.Linear(n_neurons, n_neurons), nn.Tanh()] layers.append(nn.Linear(n_neurons, 2)) # 输出 [w, M] self.net = nn.Sequential(*layers) self._init_weights() def _init_weights(self): for m in self.net: if isinstance(m, nn.Linear): nn.init.xavier_uniform_(m.weight) nn.init.zeros_(m.bias) def forward(self, x): out = self.net(x) return out[:, 0:1], out[:, 1:2] # 分别返回 w 和 M

小提示:如果以后遇到更复杂的问题,可以把输出维度再扩展,比如同时输出挠度、转角、弯矩,甚至多个物理量共享同一个特征提取器,这在多物理场耦合问题里非常有用。共享特征的好处是网络能学到物理量之间的关联,减少重复参数。

4.2 自动微分求导与损失函数实现

用torch.autograd.grad计算一阶和二阶导数。两个关键参数不能漏:一是create_graph=True,因为损失函数本身要用到导数,反向传播需要保留导数的计算图;二是grad_outputs=torch.ones_like(w),这是给每个输出分量一个初始梯度,否则标量对向量的雅可比向量积会因维度问题报错。

def pinn_loss(model, x_f, x_fixed, x_free, EI=1.0, lambda_bc=30.0): # 内点方程损失,每次迭代的x_f是新随机采样的 x_f = x_f.clone().requires_grad_(True) w, M = model(x_f) dw = torch.autograd.grad(w, x_f, grad_outputs=torch.ones_like(w), create_graph=True)[0] d2w = torch.autograd.grad(dw, x_f, grad_outputs=torch.ones_like(dw), create_graph=True)[0] dM = torch.autograd.grad(M, x_f, grad_outputs=torch.ones_like(M), create_graph=True)[0] d2M = torch.autograd.grad(dM, x_f, grad_outputs=torch.ones_like(dM), create_graph=True)[0] # 降阶后的两个控制方程残差 loss_w_eq = torch.mean((EI * d2w + M) ** 2) loss_M_eq = torch.mean((d2M - 1.0) ** 2) # q=-1 时 M'' = 1 loss_eq = loss_w_eq + loss_M_eq # 固定端边界:w(0)=0, w'(0)=0 w_fixed, M_fixed = model(x_fixed) dw_fixed = torch.autograd.grad(w_fixed, x_fixed, grad_outputs=torch.ones_like(w_fixed), create_graph=True)[0] loss_bc_fixed = torch.mean(w_fixed ** 2) + torch.mean(dw_fixed ** 2) # 自由端边界:M(L)=0, M'(L)=0 w_free, M_free = model(x_free) dM_free = torch.autograd.grad(M_free, x_free, grad_outputs=torch.ones_like(M_free), create_graph=True)[0] loss_bc_free = torch.mean(M_free ** 2) + torch.mean(dM_free ** 2) loss = loss_eq + lambda_bc * (loss_bc_fixed + loss_bc_free) return loss

这里有几个容易踩的坑。第一,x_f计算完导数后要记得clone().requires_grad_(True),否则无法对x求导。第二,边界点的x_fixed和x_free不需要重新随机采样,它们固定不变,但也必须设置requires_grad才能对输出求导。第三,二阶导是在create_graph=True的一阶导基础上再求一次自动微分,两次梯度函数嵌套,所以计算图会保留两层,这是正常的。

4.3 训练主循环与采样方法

训练使用Adam优化器,初始学习率1e-3。内点数量N_f取2000,每次迭代重新生成一批随机点。边界点各重复100份。训练轮数30000步,用ReduceLROnPlateau在损失平台期自动把学习率降低一半。

model = PINN() optimizer = torch.optim.Adam(model.parameters(), lr=1e-3) scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau( optimizer, mode='min', factor=0.5, patience=2000, verbose=True ) L = 1.0 N_f = 2000 N_bc = 100 x_fixed = torch.zeros(N_bc, 1) # 固定端 x=0 x_free = torch.ones(N_bc, 1) * L # 自由端 x=L for step in range(30000): # 每次迭代重新在区间内随机采样内点,避免过拟合到固定采样点 x_f = torch.rand(N_f, 1) * L optimizer.zero_grad() loss = pinn_loss(model, x_f, x_fixed, x_free) loss.backward() optimizer.step() scheduler.step(loss) if step % 2000 == 0: print(f"step {step:6d}, loss = {loss.item():.6f}")

每次迭代重新采样内点是PINN的标准做法,也是和传统神经网络训练最大的区别之一。固定采样点相当于把配点固定了,网络可能在采样点之间出现振荡;随机配点让网络在训练过程中见过足够多的位置,解的整体光滑性会更好。一维问题计算量小,每步重新采样几乎没有额外成本。

4.4 用解析解验证训练结果

训练结束后,在[0,1]上生成1001个均匀分布点,用网络预测w和M,同时用解析解公式计算精确解。评估指标用相对L2误差,这是PINN论文里最常见的误差指标之一。

x_test = torch.linspace(0, 1, 1001).reshape(-1, 1) w_pred, M_pred = model(x_test) # 解析解,L=1, EI=1, q0=1 w_exact = -x_test**2 * (6 - 4*x_test + x_test**2) / 24 M_exact = (1 - x_test)**2 / 2 # 相对L2误差 err_w = torch.norm(w_pred - w_exact, 2) / torch.norm(w_exact, 2) err_M = torch.norm(M_pred - M_exact, 2) / torch.norm(M_exact, 2) print(f"挠度相对L2误差: {err_w.item():.6f}") print(f"弯矩相对L2误差: {err_M.item():.6f}")

如果训练正常,挠度相对L2误差通常在0.5%到2%之间。我一次典型训练跑下来,挠度误差在0.8%左右,弯矩误差略高一点,在1.5%左右。对比一下:如果直接把四阶方程当成损失函数硬算,同配置下误差经常到5%以上甚至不收敛。降阶带来的收益非常明显。

5. 训练中的坑与排查:loss震荡、NaN和精度上不去

5.1 不收敛先查这三个地方

第一次跑PINN的人遇到loss不降或者震荡,第一反应往往是改学习率或者加大网络,其实多半问题出在更基础的地方。

一是数值量级没有归一化。如果直接拿EI=2×10¹¹、q=1×10⁴去训练,各个损失项的数量级能差出十几个数量级,梯度更新完全失控。解决办法就是把L、EI、q都变成1,先让网络学会无量纲解,再把物理单位乘回去。

二是激活函数用错。ReLU在需要二阶以上导数的PINN问题里直接废掉,因为ReLU的二阶导处处为零,网络无法学习任何曲率信息。输出层还会出现死神经元,loss一动不动。换成tanh或者SiLU之后通常立刻开始下降。

三是边界权重太小。前面说过,边界权重λ_bc默认应该在10到100之间。如果取1或者更小,固定端的转角约束基本形同虚设,画出来的曲线固定端不水平,这就是边界被内点方程损失淹没的典型表现。

5.2 NaN的来源和梯度钳制

loss变成NaN是最让人头大的问题。就我的经验,在一维悬臂梁案例里NaN通常有三种来源。

第一种,自动微分的使用方式不对。比如requires_grad没设置,或者create_graph漏了,导致某个导数计算结果为空,后续运算直接得到NaN。排查方法是在损失函数里逐行打印各导数的shape和值,看哪一步出现了inf。

第二种,学习率太大。PINN的损失面本身就比普通回归问题复杂,1e-2的学习率很容易让loss暴涨,几轮之后就掉进NaN。一维问题从1e-3起步比较安全,如果此时loss震荡,降到3e-4再试。

第三种,tanh饱和。当网络权重过大,tanh的输入落在正负几十的区间,梯度几乎为零,导数计算可能出现inf。这时可以通过更小的初始化和梯度裁剪来解决。紧急情况下可以用torch.nan_to_num把NaN梯度替换成0,先让训练跑起来,但这不是根治办法,一定要找到NaN的真正触发点。

5.3 误差评估和判定标准

PINN训练完成后,不要只看loss值,loss值低不代表解正确,因为方程残差可能被加权策略人为压低。一定要用解析解或者高精度数值解去对比。

判断标准我一般看两个指标。相对L2误差在1e-2以内,说明整体趋势和量级都对,作为原型验证已经足够。如果想追求更高精度,比如1e-3以下,需要综合提高λ_bc、增加内点数量、延长训练轮数,并且加入学习率调度。再往下就进入软约束的能力边界了,PINN毕竟不像有限元那样把边界条件强加进求解矩阵,它的边界满足是靠损失函数"劝"出来的。

我建议画三张图:预测挠度与解析解对比、预测弯矩与解析解对比、绝对误差分布。挠度误差通常在固定端附近最大,因为固定端要同时满足w=0和w'=0两个条件,两个约束的梯度在端点处会互相干扰。弯矩误差则在自由端附近更明显,那里是弯矩从负到零的过渡区域。

5.4 提升精度的小技巧

如果你按上面的代码跑通了,但觉得误差还是不够小,可以尝试以下几个技巧,按性价比排序。

第一,用学习率调度把后期学习率降到1e-5左右。Adam的默认1e-3在前5000步很好用,后期会陷入平台期,调度降到1e-4甚至1e-5之后,往往还能再明显降一截误差。第二,把λ_bc从30提高到50到100,同时把内点数量从2000提高到5000,边界约束会被强化,但内点采样密度也要跟上,否则内点方程会被忽视。第三,多次随机种子,不同初始化下误差可能有明显波动,取误差最小的模型。第四,进阶技巧是用可学习激活函数,比如把tanh的输入乘以一个可训练参数a,让网络自己调节激活函数的陡峭程度,这本质上是让网络更灵活地适应不同曲率区域。

还有一个路数值得一提:如果你对边界条件有很强的执念,可以用硬约束法把边界条件写进网络输出,比如让w的输出乘以x,保证w(0)=0。但硬约束会破坏微分方程在边界附近的自由形态,有时反而让内点误差变大。对一维梁来说,软约束加权重已经够用了。

最后我再分享一点实际体会

在这个一维悬臂梁案例里,我最大的心得是:PINN的成败往往不在网络结构,而在物理问题的建模方式。把四阶方程降成两个二阶方程之后,训练稳定性有了质的提升,这比增加网络层数有效得多。另一个体会是,无量纲化和边界权重是真正决定能否收敛的前置条件,很多人一上来就调网络结构,结果问题根本不在那里。

如果你后面想把这套方法扩展到二维板、多维结构或者反问题识别,一维悬臂梁是最好的练兵场。先把网格无关、PINN、软约束这些概念在一个60行代码的例子里吃透,再上复杂问题会顺手很多。我在做这个案例时也碰了好几次壁,等把代码稳定跑通,再回头看那段踩坑经历,反而觉得对这个方法的理解深了一层。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询