物理信息神经网络求解一维悬臂梁挠度问题

最近在做结构力学相关研究时随手拿了一个经典的一维悬臂梁挠度问题来测试物理信息神经网络PINN。这个题目看起来简单但作为PINN入门的练手案例非常合适——它既有明确的控制方程又有典型的边界条件而且解析解是现成的可以很方便地验证神经网络算得对不对。这篇文章把我从搭建模型到调通代码的完整过程记录下来包括网络结构怎么设计、边界条件怎么处理、损失函数每一项是怎么算的、以及我在训练中踩过的几个坑。如果你刚接触PINN想找一个能快速跑通又看得懂原理的案例这篇应该能帮到你。1. 问题定义与整体思路1.1 为什么选一维悬臂梁悬臂梁是结构力学里最基础也最经典的构件之一。一根梁一端固定、一端自由在自由端施加一个集中力梁会发生弯曲变形我们需要求解的是梁在每一个位置处的挠度。这个问题之所以适合做PINN的入门案例有几个原因控制方程是一个四阶常微分方程形式固定不像二维三维问题那样需要处理复杂的几何边界。边界条件非常典型固定端位移为零、转角为零自由端弯矩为零、剪力等于施加的外力。解析解可以直接由材料力学公式给出用来验证神经网络的计算精度非常方便。我当时的目标很简单用神经网络代替传统数值方法输入是梁上的位置坐标x输出是对应的挠度w(x)通过让网络输出满足控制方程和边界条件来训练模型。1.2 为什么用PINN而不是有限元传统方法里求解这种问题最常用的是有限元法。有限元需要划分网格需要对单元进行离散和组装对于复杂几何和复杂边界网格生成本身就是一个不小的工作量。PINN的思路完全不同它把控制方程和边界条件直接编码到损失函数里通过优化神经网络的参数来满足这些方程整个过程不需要任何网格划分也不需要单元矩阵组装。PINN有一个很有吸引力的特点只要损失函数下降到足够小网络给出的是一个连续函数解而不是只在离散节点上有值。这意味着你可以在任意位置查询挠度不需要插值。当然PINN也有它的短板。四阶方程本身对网络的二阶甚至三阶求导要求很高训练过程中梯度传播路径比较长收敛速度和稳定性往往不如成熟的有限元求解器。不过作为方法验证和学习理解这个案例的价值在于把PINN的整个工作流程完整走通。2. 控制方程与网络结构设计2.1 悬臂梁挠度的控制方程一维悬臂梁在自由端受集中力F作用时其挠度控制方程为$$EI \frac{d^4 w}{dx^4} 0$$其中E是弹性模量I是截面惯性矩w是挠度。这个方程看起来很简单但关键在于边界条件的组合。两端边界条件分别为固定端x0处w0dw/dx0自由端xL处弯矩为零即d²w/dx²0剪力等于外力即EI·d³w/dx³F为了计算方便我对方程做了无量纲化处理。设梁长L1mEI1自由端剪力归一化为1。这样控制方程就简化为$$ \frac{d^4 w}{dx^4} 0 $$边界条件变为w(0)0dw/dx(0)0d²w/dx²(1)0d³w/dx³(1)1无量纲化的好处不仅仅是数值上更稳定更重要的是避免了不同物理量纲带来的梯度不平衡问题。在实际计算中如果直接用真实物理参数挠度可能只有毫米级别而坐标x和挠度w之间的量级差会导致损失函数中各项的贡献严重失衡训练起来非常困难。2.2 网络结构与激活函数选择对于这个一维问题网络结构不需要很复杂。我选用了一个四层全连接网络每层50个神经元输入层维度为1坐标x输出层维度为1挠度w。激活函数的选择上有讲究。常见的ReLU虽然计算简单但它的一阶导数是分段常数二阶导数在分段点处为0三阶导数几乎到处都是0根本无法用于四阶方程的损失函数计算。我当时试过ReLU结果是损失函数根本降不下去因为网络的导数信息在反向传播中丢失了。换用tanh激活函数之后问题得到了解决。tanh是光滑函数理论上可以无限次求导而且它在各个阶数的导数都有非零值能够正常地通过自动微分计算损失函数中的物理残差。注意对涉及高阶导数的PINN问题激活函数必须选择光滑函数。tanh、sigmoid、swish这类函数都可以用从而避免ReLU族函数导数阶数限制的问题。网络结构设计的另一个关键点是如何让边界条件更容易被满足。我采用了硬约束的方式也就是把网络输出做一次变换使得输出自动满足位移和转角的边界条件$$w_{net}(x) x^3 \cdot NN(x)$$这里NN(x)是原始网络输出。因为当x0时x³0所以无论如何w_net在x0处的值和一阶导数都恒为0。这样就无需在损失函数中再考虑位移和转角的边界约束只需要约束自由端的弯矩和剪力条件。这个做法的优点很明显损失函数中需要约束的项变少了训练目标更清晰收敛速度也有明显提升。代价是网络输出和真实解答的表达空间不再完全等价但在实际测试中只要网络容量足够这个变换并不会限制解的精度。3. 边界条件处理与损失函数构建3.1 如何把边界条件编码进模型如前面所说固定端的位移和转角边界条件通过输出变换被硬编码进网络结构了。这样剩下的只有自由端的两个自然边界条件需要处理。在损失函数中分别定义自由端的弯矩残差项和剪力残差项弯矩条件网络输出在x1处的二阶导数应为0剪力条件网络输出在x1处的三阶导数应为1把这两个残差平方后加入总损失函数中就实现了对边界条件的软约束。3.2 损失函数的三项组成最终的总损失函数由三部分组成$$L L_{PDE} \lambda_1 L_{moment} \lambda_2 L_{shear}$$其中第一项是内部点的控制方程残差第二项是自由端弯矩残差第三项是自由端剪力残差。内部点损失项的计算方法是在[0,1]区间内随机采样若干点计算网络输出对这些点的四阶导数与目标值0做差的平方和。采样点的数量和分布对收敛效果有很大影响我使用的是均匀随机采样每个训练轮次重新采样相当于每次迭代都在使用一组全新的样本点。权重系数λ的选择也是调参的重点。一开始我将两个边界条件的权重都设为1训练后发现在边界附近的解答偏差比较大。后来把λ从1逐级调整到10甚至50才找到了一个合适的平衡点。实际效果是边界权重越大边界条件满足得越精确但过大会导致内部点残差被忽略整体解的光滑性变差。4. 代码实现与关键细节4.1 导入依赖和定义超参数我用PyTorch作为深度学习框架因为它的自动微分功能很强大能方便地计算任意阶数的导数。以下是完整的代码实现。import torch import torch.nn as nn import numpy as np import matplotlib.pyplot as plt # 设置随机种子保证可复现性 torch.manual_seed(42) np.random.seed(42) # 超参数定义 n_epochs 5000 n_samples 100 # 每次采样的内部点数量 lr 1e-3 # 学习率 lambda_moment 10.0 # 弯矩边界权重 lambda_shear 10.0 # 剪力边界权重4.2 网络定义class FCN(nn.Module): def __init__(self, n_input1, n_output1, n_hidden50, n_layer4): super(FCN, self).__init__() layers [nn.Linear(n_input, n_hidden)] for _ in range(n_layer - 1): layers.append(nn.Linear(n_hidden, n_hidden)) layers.append(nn.Linear(n_hidden, n_output)) self.activation nn.Tanh() self.layers nn.ModuleList(layers) def forward(self, x): for i, layer in enumerate(self.layers): x layer(x) if i len(self.layers) - 1: x self.activation(x) return x4.3 物理约束与损失计算def pde_residual(model, x): 计算控制方程残差四阶导 - 0 x x.requires_grad_(True) w model(x) # 一阶导 dw_dx torch.autograd.grad(w, x, grad_outputstorch.ones_like(w), create_graphTrue)[0] # 二阶导 d2w_dx2 torch.autograd.grad(dw_dx, x, grad_outputstorch.ones_like(dw_dx), create_graphTrue)[0] # 三阶导 d3w_dx3 torch.autograd.grad(d2w_dx2, x, grad_outputstorch.ones_like(d2w_dx2), create_graphTrue)[0] # 四阶导 d4w_dx4 torch.autograd.grad(d3w_dx3, x, grad_outputstorch.ones_like(d3w_dx3), create_graphTrue)[0] return d4w_dx4, w, dw_dx, d2w_dx2, d3w_dx3 def compute_loss(model): # 内部点残差 x_interior torch.rand(n_samples, 1) * 1.0 # [0, 1]区间均匀采样 d4w, _, _, _, _ pde_residual(model, x_interior) loss_pde torch.mean(d4w**2) # 自由端x1处边界条件 x_boundary torch.ones(1, 1) _, _, _, d2w_b, d3w_b pde_residual(model, x_boundary) loss_moment torch.mean(d2w_b**2) # 弯矩应为0 loss_shear torch.mean((d3w_b - 1.0)**2) # 剪力应为1 # 总损失 loss loss_pde lambda_moment * loss_moment lambda_shear * loss_shear return loss, loss_pde, loss_moment, loss_shear4.4 训练流程def train(model): optimizer torch.optim.Adam(model.parameters(), lrlr) loss_history [] for epoch in range(n_epochs): optimizer.zero_grad() loss, loss_pde, loss_moment, loss_shear compute_loss(model) loss.backward() optimizer.step() loss_history.append(loss.item()) if epoch % 1000 0: print(fEpoch {epoch}: loss{loss.item():.6f}, fpde{loss_pde.item():.6f}, fmoment{loss_moment.item():.6f}, fshear{loss_shear.item():.6f}) # 创建模型并训练 model FCN() train(model)4.5 解析解对比与可视化一维悬臂梁在自由端受集中力的挠度解析解表达式为$$w_{exact}(x) \frac{F x^2}{6EI}(3L - x)$$代入无量纲化参数F1, EI1, L1得到$$w_{exact}(x) \frac{x^2}{6}(3 - x)$$用训练好的模型在所有位置推断挠度值和解析解绘制在同一张图里对比# 测试点 x_test torch.linspace(0, 1, 200).reshape(-1, 1) with torch.no_grad(): w_pred model(x_test).squeeze().numpy() # 解析解 w_exact (x_test.squeeze().numpy()**2 * (3 - x_test.squeeze().numpy())) / 6 # 绘图 plt.figure(figsize(8, 5)) plt.plot(x_test.squeeze().numpy(), w_exact, k-, labelExact Solution) plt.plot(x_test.squeeze().numpy(), w_pred, r--, labelPINN Prediction) plt.xlabel(x) plt.ylabel(w(x)) plt.legend() plt.grid(True, alpha0.3) plt.title(一维悬臂梁挠度PINN vs 解析解) plt.show() # 计算最大误差 error np.abs(w_pred - w_exact).max() print(f最大绝对误差: {error:.6f})5. 训练调试与常见问题5.1 训练过程中的几个关键观察训练初期损失函数的下降速度非常快前几百轮就能看到明显的变化。这是因为硬约束固定端边界后网络主要需要学习的自由度集中在了如何满足控制方程和剩余边界条件上寻优空间相对有限。不过在损失下降到一定程度后会进入一个平台期。这个阶段损失值在小范围内波动不再明显下降。这是正常的PINN的损失函数本质上是多个约束项的加权叠加各约束之间很难同时达到完美平衡最终会在某个折中的状态稳定下来。我在测试中发现最终模型在内部点的挠度预测最大误差在0.5%左右对于工程精度要求不高、主要用于方法演示的场景是完全够用的。如果要提高精度可以增加网络宽度或深度、增加采样点数量、采用自适应采样策略或者增大边界条件的权重。5.2 一个容易被忽略的细节梯度计算在进行高阶导数计算时create_graphTrue这个参数一定不能忘。如果忘记设置计算出来的梯度不会参与后续的反向传播训练时网络参数无法正常更新损失值会一直卡住不动。我当时第一次跑代码就在这里栽了跟头找了半天才发现是这个参数的问题。# 错误写法无法正常训练 dw_dx torch.autograd.grad(w, x, grad_outputstorch.ones_like(w))[0] # 正确写法保留计算图用于后续反向传播 dw_dx torch.autograd.grad(w, x, grad_outputstorch.ones_like(w), create_graphTrue)[0]5.3 四阶导数的数值特性四阶导数在数值上是比较敏感的。如果网络输出在局部有细微的振荡四阶导数会被放大得非常明显。我在调试中发现在训练后期内部点的四阶导残差虽然在减小但仍然存在一些小的高频分量。这种情况下可以试试降低学习率或者对采样点进行平滑处理。更有效的方法是给损失函数中的PDE残差项换一种归一化方式比如用相对误差代替绝对误差。我在后续的测试中发现将PDE残差除以某个参考值后训练稳定性有比较明显的提升特别是在网络参数随机初始化的情况下。5.4 采样点的选择策略是影响收敛速度的关键因素我最初用的是整个[0,1]区间均匀随机采样。但观察训练过程后发现靠近自由端的区域对解的准确性影响更大因为那里有弯矩和剪力的边界条件需要满足。后来我尝试了一种简单的自适应采样先训练几百轮然后统计每个位置的PDE残差分布在残差较大的区域增加采样点密度。这种方法让同样的训练轮数下最大误差几乎缩小了一半。PINN的采样策略其实挺有讲究均匀采样、随机采样、自适应采样、带权重的采样各有适用场景根据问题特性选择合适的方法能事半功倍。6. 经验总结与扩展方向6.1 我对PINN求解这个问题的几点体会用PINN算悬臂梁挠度这个题目整体跑通不难但要算得准确、稳定还是有不少细节需要注意。第一个体会是边界条件的处理方式直接影响训练效率。硬约束和软约束相比虽然增加了网络设计的复杂度但给优化问题排除了一部分不确定性。对于边界条件比较规则的问题优先考虑能不能把它直接嵌入到网络输出中。第二个体会是损失函数的权重分配真的很重要。各个物理约束之间天然存在数值尺度差异如果不做加权调整训练结果会被数值量级更大的那一项主导。实际操作中建议先用若干个不同权重跑几次短训练观察各项残差的收敛情况再根据观察结果调整权重配置。第三个体会是网络结构和参数量不是越大越好。对于一维简单问题四层宽度50的网络已经足够。盲目增加网络规模不仅让训练变慢还可能出现过拟合式的振荡反而让导数的计算结果变差。6.2 这个案例还能怎么扩展这个一维悬臂梁案例虽然简单但它包含的方程形式、边界条件处理方式、损失函数构建思路可以作为一系列更复杂问题的起点。比如可以推广到变截面梁也就是截面惯性矩I(x)沿梁长变化控制方程中的系数不再是常数此时四阶导项会包含位置相关的系数但在PINN的框架下修改起来仅仅是损失函数里增加一个系数函数而已。还可以考虑加上分布荷载、非线性材料本构关系这类问题在传统数值方法中往往需要反复迭代但在PINN框架中只需要调整控制方程的残差表达式即可。更进一步的可以用类似框架处理二维弹性力学问题甚至扩展到基于观测数据做参数辨识。比如用已知的挠度测量值来反推弹性模量E或惯性矩I这就从正问题走向了反问题PINN在这个方向上的优势会更明显。总的来说PINN为求解物理问题提供了一种全新的思路。传统方法需要网格、离散、组装而PINN只需要一个神经网络和一套物理约束。它仍然有计算效率、收敛性等方面的不足但在某些问题上已经展示出独特的价值。这个悬臂梁案例是我入门PINN的第一个完整实现也让我对这个方向的潜力有了更直观的认识。