做过结构分析的朋友应该都有体会悬臂梁挠度查公式快但真要自己写微分方程求解器又要离散、又要迭代、又要处理边界条件绕得很。最近我一直在折腾物理信息神经网络PINN索性把一维悬臂梁挠度问题整个跑了一遍从建模、写损失函数到Python代码实现都整理在这里。这篇文章不是泛泛讲原理而是直接围绕“给定一个悬臂梁用PINN求出挠度曲线”展开。全文包含方程推导、PyTorch实现、训练调参、解析解对比还有我实测踩过的几个坑。适合具备基础力学概念、想用PINN做工程微分方程求解的读者哪怕你之前没接触过神经网络跟着代码跑一遍也能上手。1. 从结构力学方程到PINN损失函数悬臂梁问题的建模1.1 一维悬臂梁挠度方程与四种边界条件悬臂梁一端固定、一端自由是最经典的欧拉-伯努利梁模型。对等截面梁小变形条件下挠度 (w(x)) 满足四阶常微分方程[ EI \frac{d^4w}{dx^4} q(x) ]其中 (EI) 是弯曲刚度(q(x)) 是分布载荷。为了把问题定下来需要四种边界条件。悬臂梁的左端固定挠度和转角都为零[ w(0)0,\quad \frac{dw}{dx}(0)0 ]右端自由弯矩和剪力为零[ EI \frac{d^2w}{dx^2}(L)0,\quad EI \frac{d^3w}{dx^3}(L)0 ]这里 (xL) 是梁的自由端。如果是均布载荷 (q(x)q_0)解析解可以直接写出来[ w(x)\frac{q_0 x^2}{24EI}(x^2-4Lx6L^2) ]这个公式很多手册都有正好可以拿来当PINN预测结果的对照组。1.2 PINN核心思想把微分方程变成优化问题传统数值解法比如有限差分、有限元核心是“离散网格求解线性方程组”。PINN的思路完全不同它用一个神经网络来近似挠度函数 (w(x))然后让这个近似函数满足两个目标——控制方程和边界条件。具体做法是构造一个总损失函数。控制方程残差表示“网络输出的函数代入微分方程后离等式成立差了多远”边界条件残差表示“输出在边界上离已知约束差了多远”。训练过程就是不断调整网络参数让总损失越来越小。对于悬臂梁问题可以定义[ \mathcal{L} \mathcal{L}{PDE} \lambda \mathcal{L}{BC} ]其中[ \mathcal{L}{PDE} \frac{1}{N}\sum{i1}^{N}\left(EI \frac{d^4w}{dx^4}(x_i)-q(x_i)\right)^2 ][ \mathcal{L}_{BC} w(0)^2 \left(\frac{dw}{dx}(0)\right)^2 \left(EI \frac{d^2w}{dx^2}(L)\right)^2 \left(EI \frac{d^3w}{dx^3}(L)\right)^2 ](x_i) 是从梁内部采样的若干点(N) 是采样点个数。网络对 (w) 的导数不是用差分算的而是通过自动微分完成的这也是“物理信息”的关键不需要人工构造导数离散格式只要损失函数里出现导数项求梯度时就能精确算出来。1.3 四阶方程的两个坑高阶导数与量级缩放理论上直接对四阶方程建损失函数就能算但我实际测试下来有两个问题非常明显。第一个是高阶自动微分的稳定性。PyTorch里对网络输出求四阶导意味着要连续调用torch.autograd.grad四次。网络一开始是随机的四阶导数值往往非常大或非常小梯度更新一次就可能让参数爆炸。即便不爆炸训练后期也容易在局部震荡。第二个是量级差异。真实工程里(L) 可能是3米(EI) 可能是几千到几万(q_0) 可能是几十千牛每米。把这些原始数值直接丢给神经网络输入输出尺度不匹配损失函数里各项数量级差十万八千里很难收敛。所以我在实现时做了两件事一是把方程降阶避免直接算四阶导二是把变量无量纲化把长度和挠度都缩放到0到1附近。这两步做完训练稳定性会好非常多。2. Python实现用PyTorch搭建能算挠度的网络2.1 环境版本与依赖先说环境。我的运行环境是Python 3.10PyTorch 2.0以上NumPy和Matplotlib用于数值计算和绘图。如果用CPU训练这个小问题几秒钟就能跑起来不需要GPU。安装依赖很简单pip install torch numpy matplotlib代码里没有用到任何特殊网络库所以版本兼容性很好。如果你用的是旧版本PyTorch自动微分的接口稍有差异但核心写法是一致的。2.2 网络结构设计双输出网络如何避免梯度混乱对于方程降阶后的悬臂梁问题网络要做的不只是输出一个挠度值而是要同时输出“挠度”和“弯矩”两个物理量。降阶方式是这样的令梁截面弯矩为 (M(x)EI \frac{d^2w}{dx^2})那么原四阶方程可以拆成两个二阶方程[ \frac{d^2 M}{dx^2} q(x) ][ \frac{d^2 w}{dx^2} \frac{M}{EI} ]边界条件也相应变成[ w(0)0,\quad \frac{dw}{dx}(0)0,\quad M(L)0,\quad \frac{dM}{dx}(L)0 ]这样一来神经网络只需要对输出求二阶导不再碰四阶导。网络设计成一个三层全连接输入 (x)输出两个节点一个对应 (w)一个对应 (M)。import torch import torch.nn as nn class PINN梁(nn.Module): def __init__(self): super().__init__() self.net nn.Sequential( nn.Linear(1, 50), nn.Tanh(), nn.Linear(50, 50), nn.Tanh(), nn.Linear(50, 50), nn.Tanh(), nn.Linear(50, 2) ) def forward(self, x): out self.net(x) w out[:, 0:1] M out[:, 1:2] return w, M这里激活函数选用Tanh而不是ReLU。原因很简单Tanh的光滑性更好高阶导不会出现拐点噪声ReLU的二阶导处处为零直接没法用。2.3 自动微分构造损失二阶导数残差为了让网络学习物理规律我们需要计算 (w) 和 (M) 对 (x) 的导数。给定输入 (x) 并开启计算图可以通过autograd.grad求导def 二阶导(u, x): du_dx torch.autograd.grad(u, x, grad_outputstorch.ones_like(u), create_graphTrue)[0] d2u_dx2 torch.autograd.grad(du_dx, x, grad_outputstorch.ones_like(du_dx), create_graphTrue)[0] return d2u_dx2注意create_graphTrue必须保留否则二次求导时梯度会断掉。损失函数由控制方程残差和边界条件残差两部分组成。3. 训练策略采样点布局、损失权重与学习率3.1 采样点布局为什么边界点一定要参与训练PINN训练时我们需要确定哪些 (x) 位置参与损失计算。最简单的做法是均匀取点比如在 ([0,1]) 上取100个点。但实际训练中我发现均匀取点虽然简单边界附近的拟合容易出问题。原因是边界条件只在一两个点上强约束如果内部点太多、边界点太少优化器会偏向拟合内部方程而忽视边界。解决办法是内部采样点使用随机采样数量100到200边界点单独加入每个边界条件对应一个点每次迭代可以重新随机采样提高不同位置的覆盖能力。我用的是“内部均匀随机点 边界固定点”的组合方式。代码里把边界点混入同一个张量用掩码区分内部点和边界点这样反向传播时一次就能算出所有loss。3.2 损失权重如何控制多项损失的相对大小悬臂梁问题的总损失至少包含五部分方程残差两个弯矩方程、挠度方程、边界条件四个。如果不加权重边界条件可能被淹没。我推荐的初始权重设置是权重_方程 1.0 权重_边界 10.0这样做的原因是边界条件数量少每个点都需要更强的约束。实测下来边界权重取10而不是1固定端的零挠度和零转角满足得更干净。如果你发现某个边界条件始终不满足单独提高它对应的权重即可。比如固定端转角一直不收敛就把转折边界权重从10调到50这是最快的临时补救办法。3.3 训练过程监控loss分项记录训练时不要只盯总loss。总loss下降不代表每个物理条件都满足。我习惯每个epoch记录四个分项loss方程残差、边界w、边界M、总loss并定期打印。打印内容示例 epoch 1000, loss5.2e-4, pde3.1e-4, bc_w2.0e-5, bc_M1.8e-4这样能迅速看出是哪项没收敛。如果bc_w始终不掉说明固定端边界的权重太小或者网络容量不够如果pde下降后又反弹大概率是学习率太大。4. 完整可运行代码与结果验证4.1 无量纲化处理网上很多PINN示例直接用工程单位但我觉得无量纲化是必要的一步。把长度缩放到 ([0,1])[ X\frac{x}{L} ]把挠度和弯矩无量纲化[ \bar{w}\frac{EI}{q_0 L^4}w,\quad \bar{M}\frac{M}{q_0 L^2} ]这样原始的方程变成[ \frac{d^2 \bar{M}}{dX^2}1 ][ \frac{d^2 \bar{w}}{dX^2}\bar{M} ]边界条件保持不变。训练结束后再通过逆变换还原到物理量。这个变换非常关键它能让神经网络输入输出都在0.1到1量级训练速度明显加快。4.2 完整代码下面这份代码我把无量纲化、网络定义、训练和绘图都写在了一个文件里保存成pinn_cantilever.py后直接运行即可。import torch import torch.nn as nn import numpy as np import matplotlib.pyplot as plt torch.manual_seed(42) np.random.seed(42) # 无量纲参数梁长L1均布荷载q01刚度EI1 L 1.0 q0 1.0 EI 1.0 # 解析解无量纲 def exact_w(X): return X**2 * (X**2 - 4*X 6) / 24.0 def exact_M(X): return 0.5 * (1 - X)**2 class PINN梁(nn.Module): def __init__(self): super().__init__() self.net nn.Sequential( nn.Linear(1, 50), nn.Tanh(), nn.Linear(50, 50), nn.Tanh(), nn.Linear(50, 50), nn.Tanh(), nn.Linear(50, 2) ) def forward(self, x): out self.net(x) w out[:, 0:1] M out[:, 1:2] return w, M def 二阶导(u, x): du_dx torch.autograd.grad(u, x, grad_outputstorch.ones_like(u), create_graphTrue)[0] d2u_dx2 torch.autograd.grad(du_dx, x, grad_outputstorch.ones_like(du_dx), create_graphTrue)[0] return d2u_dx2 # 生成内部点和边界点 def 生成样本(n_inner100): x_inner torch.rand(n_inner, 1) * 0.98 0.01 # 避开边界点 x_bc torch.tensor([[0.0], [1.0]]) x torch.cat([x_inner, x_bc], dim0) x x.clone().detach().requires_grad_(True) return x net PINN梁() optimizer torch.optim.Adam(net.parameters(), lr1e-3) scheduler torch.optim.lr_scheduler.StepLR(optimizer, step_size3000, gamma0.5) epochs 10000 权重_方程 1.0 权重_边界 10.0 for epoch in range(epochs): optimizer.zero_grad() x 生成样本(120) w, M net(x) # 方程残差 M_xx 二阶导(M, x) w_xx 二阶导(w, x) loss_pde_M torch.mean((M_xx - 1.0)**2) loss_pde_w torch.mean((w_xx - M)**2) # 边界残差w(0)0, w(0)0, M(1)0, M(1)0 x0 torch.tensor([[0.0]], requires_gradTrue) x1 torch.tensor([[1.0]], requires_gradTrue) w0, M0 net(x0) w1, M1 net(x1) dw0_dx torch.autograd.grad(w0, x0, grad_outputstorch.ones_like(w0), create_graphTrue)[0][0, 0] dM1_dx torch.autograd.grad(M1, x1, grad_outputstorch.ones_like(M1), create_graphTrue)[0][0, 0] loss_bc_w w0[0, 0]**2 dw0_dx**2 loss_bc_M M1[0, 0]**2 dM1_dx**2 loss 权重_方程 * (loss_pde_M loss_pde_w) 权重_边界 * (loss_bc_w loss_bc_M) loss.backward() optimizer.step() scheduler.step() if (epoch 1) % 1000 0: print(fepoch {epoch1:6d}, loss{loss.item():.3e}, pde{loss_pde_M.item()loss_pde_w.item():.3e}, fbc_w{loss_bc_w.item():.3e}, bc_M{loss_bc_M.item():.3e}) # 验证结果 x_test torch.linspace(0, 1, 201).reshape(-1, 1).requires_grad_(True) w_pred, M_pred net(x_test) w_pred w_pred.detach().numpy().flatten() M_pred M_pred.detach().numpy().flatten() x_np x_test.detach().numpy().flatten() w_exact exact_w(x_np) M_exact exact_M(x_np) plt.figure(figsize(6, 4)) plt.plot(x_np, w_pred, r--, labelPINN) plt.plot(x_np, w_exact, b-, labelExact) plt.xlabel(x/L) plt.ylabel(w EI/(qL^4)) plt.legend() plt.title(Cantilever Deflection) plt.grid(True) plt.savefig(deflection.png, dpi150) plt.figure(figsize(6, 4)) plt.plot(x_np, np.abs(w_pred - w_exact), k-) plt.xlabel(x/L) plt.ylabel(abs error) plt.title(Deflection Error) plt.grid(True) plt.savefig(error.png, dpi150) plt.show()这份代码我实际跑过到5000步左右loss就降到1e-5量级最终挠度曲线和解析解几乎重合。如果你的环境和我的略有差异可以适当调大学习率衰减步长或者把网络宽度从50降到32。4.3 与解析解对比误差分布训练完成后我对比了自由端挠度。无量纲结果里自由端 (X1) 的解析解为[ \bar{w}(1)\frac{1}{24}(1-46)\frac{3}{24}0.125 ]PINN预测值通常在0.1248到0.1252之间误差不到0.2%。最大误差往往出现在自由端附近因为那里的二阶导数变化快网络需要更多采样点去拟合。如果要把结果换算回真实物理挠度只需要用[ w(x)\frac{q_0 L^4}{EI}\bar{w}(x/L) ]例如 (L3\text{m}, q_010\text{kN/m}, EI2000\text{kN·m}^2)则真实自由端挠度[ w(L)\frac{q_0 L^4}{EI}\times 0.125\frac{10\times81}{2000}\times0.1250.0506\text{m} ]和手算结果完全一致。5. PINN实操中常见的翻车点与解决办法5.1 固定端边界漂移最常遇到的问题是训练一段时间后方程残差已经很小但 (w(0)) 不为零或者转角在固定端不是零。这不是网络算不对而是边界条件在总loss里占比太低。解决办法有两个方向。第一把边界权重从10提高到50甚至100。第二把边界点也纳入内部点列表让固定端的残差在每次迭代都被重复计算。注意不要为了让边界严格满足而把权重设得过大否则方程残差会被忽略最终结果会变成一条没有物理意义的曲线。5.2 训练初期残差震荡我刚开始训练时loss曲线在前几百步会剧烈震荡甚至可能出现NaN。原因通常是网络初始化后二阶导数值太大导致梯度爆炸。解决方法很直接使用较小的学习率比如1e-3启用学习率衰减每3000步衰减一半网络权重初始化采用Xavier或默认均匀初始化不要手动设过大如果出现NaN把批大小调小或者把网络层数从3层降到2层。5.3 降阶方法 vs 直接四阶我也试过直接对四阶方程做PINN。理论上可行但实践中四阶自动微分计算图会非常深损耗值经常卡在 (10^{-3}) 量级就下不去了。相比之下降阶成两个二阶方程后精度能到 (10^{-5}) 以下。所以我的建议很明确只要能降阶就降阶。这不只是把高阶导数拆成低阶更关键的是让网络在不同物理量之间建立耦合关系。网络输出的挠度和弯矩互相约束反而更容易学到真实的物理响应。5.4 采样点数量的选择采样点太少比如只有20个内部点网络会过拟合到采样点位置导致曲线局部波动采样点太多比如500个单步训练变慢但精度提升不明显。这个题目里100到150个内部点是性价比最高的区间。边界点每个边界放1个就够因为边界条件是硬性约束不需要像内部方程那样统计平均。5.5 单位统一是一个隐形坑如果你的输入输出不统一量纲模型有可能看起来收敛但预测值完全不对。比如用毫米和米混着输网络很难学到正确函数关系。归一化不只是让数值更好看更重要的是消除单位带来的尺度偏差。实际操作中我建议把所有物理量归一到0到1范围内再送进网络。训练完再做逆变换。这样还有一个额外好处无论你算的是梁长1米还是10米模型都可以复用输入输出范围保持一致。这套PINN流程跑通之后我觉得最有价值的不是这个悬臂梁算例本身而是“怎么把物理方程写成损失函数”的思维方式。以后遇到变截面梁、非线性材料甚至二维板问题只需要把控制方程、边界条件和几何参数替换掉其他框架可以原封不动地搬过去。最后再分享一个小技巧训练结束后不要只画预测曲线一定要把误差单独画出来。误差集中在哪里哪里就是采样点或边界权重要调整的地方。我靠这张误差图避开了很多“看起来收敛但实际偏差不小”的迷之情况。