最近一直在复盘我用 PINN物理信息神经网络算一维悬臂梁挠度的完整过程。为什么选这个题目因为悬臂梁是结构力学里最经典的问题之一控制方程明确、边界条件清晰、解析解现成拿来做 PINN 的“hello world”再合适不过。这篇文章不绕弯子直接讲清三件事怎么把“求挠度”变成“训练一个神经网络”、核心代码怎么落地、以及我踩过的坑。项目的核心其实很简单给定一根一端固定、一端自由的悬臂梁受均布载荷作用求解挠曲线 w(x)。传统做法是划分网格、组装刚度矩阵我这里用 Python 写一个小型神经网络输入坐标 x输出挠度 w训练的时候完全不喂标签数据而是把欧拉-伯努利梁方程直接写进损失函数约束网络输出满足物理规律。跑完之后网络在任意位置给出的挠度与解析解最大误差能到 10⁻⁴ 量级。这篇文章适合两类人一类是已经会神经网络基本操作、想往“AI 科学计算”方向迈一步的开发者另一类是做结构分析、想了解 PINN 到底能不能用来算力学问题的工程师。下面按我实际操作的顺序把原理、代码、坑位一次说清楚。1. 把“悬臂梁挠度”变成一个神经网络能算的问题1.1 控制方程与两类边界条件一维悬臂梁的小挠度问题采用欧拉-伯努利梁理论控制方程为EI · d⁴w/dx⁴ q(x)其中 E 是弹性模量I 是截面惯性矩EI 是抗弯刚度q(x) 是分布载荷w(x) 是挠度。这里先解释一下为什么是四阶导数。材料力学里三个基本关系是弯矩 M EI·d²w/dx²剪力 V EI·d³w/dx³分布载荷 q dV/dx EI·d⁴w/dx⁴。也就是说四阶导数直接对应载荷这是整个 PINN 损失函数里最关键的一项。边界条件分两类。固定端在 x0位移为零即 w(0)0转角为零即 dw/dx(0)0。自由端在 xL弯矩为零即 EI·d²w/dx²(L)0剪力为零即 EI·d³w/dx³(L)0。有一个非常常见的错误初学者只写位移边界 w(0)0 和弯矩边界 EI·d²w/dx²(L)0把转角条件和剪力条件漏了。结果网络确实给出了一条“经过固定端、自由端弯矩为零”的光滑曲线但固定端转角明显不为零挠曲线的形状完全不对。这四个边界条件一个都不能少。1.2 为什么用 PINN 而不用传统数值方法有限元求解悬臂梁精度很高成熟度也很高。那为什么还要用 PINN我的答案很直接不是为了替代有限元而是为了训练一个“物理约束下的可微函数”让它能嵌入到更大的深度学习系统里。有限元把整个梁离散成单元需要划分网格、组装刚度矩阵。一旦结构几何发生变化、载荷分布变化往往要重新建模。PINN 的思路完全不同它把求解过程转成一个优化问题函数本身用神经网络表示物理方程作为正则项。网络训练好之后任意给一个坐标立刻输出挠度而且输出是连续光滑的。更重要的是PINN 天然适合反问题。假设已经测得自由端挠度想反推分布载荷 q 的大小传统方法要做逆求解往往困难PINN 里只需要把 q 也设成网络的一个输入或者待优化变量训练时用实测挠度作为标签就能同时拟合出载荷。这类“物理约束 数据驱动”的范式才是 PINN 真正的价值所在。当然纯粹计算一根等截面悬臂梁的挠度PINN 在效率上完全比不上解析解或有限元。这一点先说明白免得读者抱着“替代 ANSYS”的期待进来结果失望。1.3 一个自带标准答案的验证案例为了验证 PINN 算得对不对必须选一个有解析解的工况。设梁长 L1 m抗弯刚度 EI1000 N·m²均布载荷 q₀1000 N/m。解析解为w(x) q₀·x²·(6L² - 4Lx x²) / (24·EI)代入参数后最大挠度在自由端 xL 处w_max q₀·L⁴ / (8·EI) 1000 / 8000 0.125 m这个值可以当作标准答案。训练完成后我把 PINN 预测值和这个解析解逐点对比就能量化误差。之所以选这组参数是因为量级合适挠度在 0.1 米左右EI 和 q₀ 的比值让四阶导数量级不至于过大或过小训练稳定性好。2. PINN 求解的完整思路拆解2.1 网络结构给神经网络“配一张光滑的嘴”我先说结论输入层 1 个节点4 个隐藏层每层 50 个神经元激活函数用 tanh输出层 1 个节点且不带激活函数。这个配置不是拍脑袋定的参考了 PINN 早期论文里的经典设置对一维问题非常够用。隐藏层再多、神经元再宽精度并不会明显提升反而会让优化更困难。我试过 8 层 128 神经元训练时间翻了好几倍精度基本没变化最后又回到 4×50。激活函数选 tanh 是有讲究的。ReLU 在负半轴导数为零分段线性二阶以上导数要么不存在要么处处为零根本没法用来计算四阶导数。sigmoid 两端饱和高阶导数非常容易消失导致网络深层梯度传不回去。tanh 处处光滑任意阶导数存在且计算稳定是目前 PINN 最常用的激活函数。另外一个小细节输出层初始化要格外小心。如果输出层权重初始值太大网络初始输出可能直接到几十的量级边界损失瞬间爆炸。我用的是标准差 0.05 的正态分布初始化相当于让初始预测值保持在很小的范围然后再让物理损失逐步“拉”到真实解上。2.2 损失函数把物理方程变成训练目标PINN 的核心就是损失函数。整个损失由两部分组成PDE 残差和边界条件残差。PDE 残差项在梁内部采样 N 个点对每个点计算网络输出对输入的导数代入控制方程L_pde mean( (EI·d⁴w/dx⁴ - q₀)² )这一项的意义是网络输出虽然在任意点都成立但它必须满足梁的平衡方程。如果网络学出一条不符合物理规律的曲线这一项就会很大优化器会把它拉回来。边界条件残差项把四个边界条件全部写成零均值误差L_bc w(0)² (dw/dx(0))² (EI·d²w/dx²(L))² (EI·d³w/dx³(L))²总损失为 L L_pde λ·L_bc。λ 是边界权重我一般取 10 到 50。为什么边界权重要比 PDE 项高一些因为边界条件只有孤立的几个点PDE 项有两千个采样点如果不提高边界权重网络会优先满足大多数点的方程边界的约束就被“淹没”了。这个思想类似罚函数法。使用均方误差而不是绝对值误差是因为 MSE 对误差的梯度是线性的优化时更平滑绝对值误差在零点不可导容易导致损失函数在边界处振荡。2.3 量纲处理一个非常容易翻车的坐标归一化细节这一节是整个实现里最容易出错的地方必须单独拎出来说。神经网络输入坐标时我通常先做归一化令 x_norm x / L把输入控制在 [0, 1] 区间。这一步本身没问题但问题出在导数换算上。如果网络学习的函数是 ŵ(x_norm)那么对真实坐标 x 求导时每求一阶导数就要多乘一个 1/L。四阶导数的换算是d⁴w/dx⁴ (1/L⁴) · d⁴ŵ / dx_norm⁴所以损失函数里的 PDE 残差应该写成pde_res EI / L⁴ · d⁴ŵ/dx_norm⁴ - q₀而不是直接在归一化坐标下用 EI·d⁴ŵ/dx_norm⁴ - q₀。如果你漏掉 1/L⁴相当于把抗弯刚度放大了 L⁴ 倍网络为了平衡损失会学出一个幅值差很远的解。我见过不止一个人卡在这个问题上训练过程完全收敛损失也很低但算出来的挠度比解析解大了一万倍。同理边界条件里的二阶导和三阶导也要对应换算d²w/dx² (1/L²)·d²ŵ/dx_norm²d³w/dx³ (1/L³)·d³ŵ/dx_norm³。2.4 训练策略先 Adam 快速下降再精细收敛训练过程我分两个阶段。第一阶段用 Adam 优化器学习率 1e-3训练 20000 到 50000 步。Adam 的好处是收敛快对初始学习率不敏感前几千步就能把损失从 10⁴ 量级压到 1 以下。第二阶段切换到 L-BFGS这个优化器利用二阶曲率信息在接近极小值点时收敛非常准确能进一步把损失降到 1e-6 量级。如果嫌麻烦也可以只用 Adam把学习率在 20000 步后手动降到 1e-4再跑 20000 步效果差不太多。我跑这个例子时主要用 Adam最后误差到 1e-4 量级就够用了。训练时采样点数量也很关键。梁内点我取 2000 个边界点固定取 x0 和 xL。2000 个内点对一维问题来说足够再多并不会带来明显收益反而增加单步计算时间。有一点值得注意一维问题计算量小所以训练时直接全量使用所有采样点不用小批量。引入 mini-batch 反而会让 PDE 残差的估计产生波动边界约束的稳定性变差。3. 代码实现从采样点到误差曲线3.1 环境准备与参数设定我用的环境是 Python 3.9 TensorFlow 2.10 NumPy Matplotlib。PyTorch 也能实现同样的功能但 TensorFlow 的 GradientTape 写高阶导数更直观代码量少所以我这里用 TensorFlow 演示。有一个非常重要的前置设定全程使用 float64 精度。四阶导数需要嵌套四层梯度计算数值范围被极度放大float32 非常容易出现 NaN。所有张量在创建时我都显式指定 dtypetf.float64。物理参数直接定义成全局常量L 1.0 # 梁长单位 m EI 1000.0 # 抗弯刚度单位 N·m² q0 1000.0 # 均布载荷单位 N/m3.2 生成训练数据和搭建网络采样分两部分梁内点用均匀随机采样让网络在 [0, L] 范围内都受到 PDE 约束边界点直接固定为两个端点。import numpy as np import tensorflow as tf N_INNER 2000 # 内部采样点归一化坐标下采样范围 [0, 1] x_inner tf.random.uniform( (N_INNER, 1), minval0.0, maxval1.0, dtypetf.float64 ) # 边界点x0固定端和 x1自由端 x_bc tf.constant([[0.0], [1.0]], dtypetf.float64)这里统一使用归一化坐标也就是 x_norm x/L。因为 L1所以内部采样点直接在 [0,1] 里取即可。如果 L 不是 1一定要记得在输入网络之前先除以 L。网络用 Keras Sequential 搭建。四个隐藏层各 50 个神经元激活函数 tanh输出层线性def build_model(): model tf.keras.Sequential([ tf.keras.layers.Dense(50, activationtanh, dtypetf.float64), tf.keras.layers.Dense(50, activationtanh, dtypetf.float64), tf.keras.layers.Dense(50, activationtanh, dtypetf.float64), tf.keras.layers.Dense(50, activationtanh, dtypetf.float64), tf.keras.layers.Dense(1, dtypetf.float64) ]) return model model build_model()输出层不指定激活函数也就是线性输出。这里每个 Dense 层都声明了 float64确保整个计算图保持高精度。如果不加这个声明输入虽然是 float64但某些内置层可能会把张量转回 float32计算四阶导数时容易出现精度问题。3.3 损失函数与训练循环的核心代码计算损失时需要把内部点和边界点拼在一起一次性算出所有点的导数再分别提取对应位置的导数组成损失项。这样只用一次前向传播效率更高。def pinn_loss(model, x_inner, x_bc): # 拼接全部采样点 x_all tf.concat([x_inner, x_bc], axis0) # 嵌套四层 GradientTape 求 d4w/dx4 with tf.GradientTape() as t1: with tf.GradientTape() as t2: with tf.GradientTape() as t3: with tf.GradientTape() as t4: t4.watch(x_all) w model(x_all) dw t4.gradient(w, x_all) d2w t3.gradient(dw, x_all) d3w t2.gradient(d2w, x_all) d4w t1.gradient(d3w, x_all) # PDE 残差注意归一化坐标下的导数换算 d4w_inner d4w[:N_INNER] pde_res EI / L**4 * d4w_inner - q0 loss_pde tf.reduce_mean(tf.square(pde_res)) # 边界条件残差 w_0 w[N_INNER] dw_0 dw[N_INNER] / L w_L w[N_INNER 1] d2w_L d2w[N_INNER 1] / L**2 d3w_L d3w[N_INNER 1] / L**3 loss_bc ( tf.square(w_0) tf.square(dw_0) tf.square(EI * d2w_L) tf.square(EI * d3w_L) ) return loss_pde 20.0 * loss_bc代码里有两个关键点。第一个是嵌套梯度。d4w 需要连续四层 GradientTape每一层记录上一层的导数结果然后对归一化坐标再求一次梯度。顺序不能乱最内层求 w 对 x 的一阶导最外层求三阶导对 x 的导数得到四阶导。这种写法是 TensorFlow 官方推荐的高阶导数计算方式比在同一个 tape 上反复调用 gradient 更稳定。第二个是导数换算。dw_0 除以 Ld2w_L 除以 L²d3w_L 除以 L³pde_res 里的四阶导数除以 L⁴。这些系数是物理单位换算不是人为调参。训练循环代码如下optimizer tf.keras.optimizers.Adam(learning_rate1e-3) tf.function def train_step(): with tf.GradientTape() as tape: loss pinn_loss(model, x_inner, x_bc) grads tape.gradient(loss, model.trainable_variables) optimizer.apply_gradients(zip(grads, model.trainable_variables)) return loss for epoch in range(30000): loss train_step() if epoch % 2000 0: print(epoch:, epoch, loss:, float(loss.numpy()))加 tf.function 可以让 TensorFlow 把整个训练步骤编译成计算图显著加速。一维问题哪怕不加也能跑只是慢一些。训练过程中我观察到损失收敛速度大致是这样前 500 步损失从 1e4 量级快速降到 1 左右5000 步后降到 1e-320000 步后稳定在 1e-5 到 1e-6 之间。如果看到 loss 在 1e-3 附近震荡不下降可以等训练结束后把学习率调到 1e-4 再跑一轮。3.4 结果对比和解析解对表训练完成后在 [0, 1] 上取 100 个均匀点输入网络得到预测值同时用解析解公式算出理论值。对比结果如下位置 x (m)PINN 预测挠度 (m)解析解挠度 (m)绝对误差 (m)00.00000.0000~00.250.01340.01322e-40.500.04410.0443-2e-40.750.08370.08352e-41.000.12520.12502e-4这是我的实测记录不同随机种子会有一定波动但整体误差稳定在 1e-4 量级。注意到误差在自由端附近略大。原因很简单自由端同时要满足弯矩为零和剪力为零两个边界条件约束比固定端更强网络拟合难度更高。这个现象是正常的。画图的时候预测值和解析解两条曲线会几乎完全重合肉眼看不出差别。更有用的是画误差曲线也就是预测值减解析解的残差沿 x 的分布。如果残差呈现高频振荡状说明网络容量不够或者隐藏层神经元数量太少如果残差在两端明显偏大说明边界权重还不够高。4. 踩坑实录这五个问题我被问过无数次4.1 损失变成 NaN这是 PINN 新手遇到最多的一个问题大概率出现在前三步训练。原因通常是三个精度不足、学习率过大、网络初始化太激进。排查方法也按这三个方向走。首先把整个网络改成 float64包括输入张量和所有 Dense 层的 dtype其次把学习率从 1e-3 降到 1e-4最后把边界权重从 20 降到 1 试试。如果这些都做了仍然 NaN尝试把隐藏层从 50 个神经元降到 30 个。四阶导数的计算图比一阶导数深得多数值范围极其敏感网络越宽越容易出问题。还有一个我自己的经验如果把输入坐标归一化到 [-1, 1] 而不是 [0, 1]损失下降会更平缓NaN 概率也会小一些。这算是隐藏技巧。4.2 边界条件怎么都压不下去训练完发现固定端位移不是严格的零而是 0.002 甚至更大。这个问题的根源几乎永远是边界约束权重不够或者边界点数量不够。解决方案有两种。第一种最简单把 x0 和 xL 各复制 20 份拼到采样点里相当于人为增加边界点的“密度”让优化器更重视这些点。第二种更彻底使用硬约束编码把网络输出修改成 w_pred x_norm² · NN(x_norm)这样 x0 时 w0 自动满足而且 dw/dx0 也自动满足固定端两个边界条件直接从软约束变成了硬约束。硬约束编码的缺点是x² 因子在 0 附近抑制了网络的输出变化网络需要更多迭代才能学习到自由端的变形收敛会慢一些。两种方案我都试过实际项目中我通常先加权重效果不够再上硬约束。4.3 算出来的挠度比理论值整体偏大或偏小如果边界条件都满足PDE 损失也降到了 1e-5 以下但整条挠曲线系统性偏离解析解那么十有八九是归一化坐标的导数换算出了漏项。我复盘过同一个 bug忘了在 pde_res 里除以 L⁴损失照样收敛因为网络完全可以通过学习一个“假刚度”来补偿这个缺失的系数最后损失很低解却完全错误。这不是一个会报错的 bug而是一个静默的逻辑错误所以特别隐蔽。排查方法很简单分别打印 loss_pde 和 loss_bc观察两者量级。理想状态下两个损失都能降到 1e-4 以下。如果 loss_pde 已经很小但解的幅值不对就去检查导数换算公式特别是 L 的幂次。4.4 多物理参数训练怎么设计采样如果你想让一个网络在多个载荷或多种刚度下都能预测挠度而不是为每套参数重新训练可以把物理参数也作为网络输入。输入从 (x) 变成 (x, q)输出还是 w损失函数里的 q 直接从输入中读取。例如想训练网络处理 q∈[500, 1500] 范围内的均布载荷只需在采样时每次随机生成一个 q连同 x 一起喂给网络。这样网络学习到的不是一个固定载荷下的解而是一个“以载荷为条件”的解函数。对这个具体问题还有一个简便技巧因为悬臂梁是线性系统挠度与载荷 q 成正比、与刚度 EI 成反比所以训练时只需变化载荷 q网络就能自然外推到其他 q 值。如果要做损伤检测把 EI 也作为输入然后用少量实测挠度反推 EI 的分布这就从正问题跨到了反问题。4.5 效率问题什么时候不该用 PINN必须说一句实话单纯计算一根等截面悬臂梁的挠度PINN 在效率上没有任何优势。解析解一秒出结果有限元一分钟建模加计算而 PINN 要训练几万步耗时以分钟计。没有必要为这种简单问题引入神经网络。PINN 真正值得用的场景是正向与逆向耦合的问题根据少量实测挠度数据反推载荷分布识别梁的刚度退化位置或者在测量数据稀疏时用物理方程作为约束来补充数据不足。比如某个结构健康监测场景传感器只能测到几个离散点的位移直接做插值误差很大但用 PINN 把梁方程作为先验网络就能在物理合理的范围内补齐整条挠曲线同时推断出损伤位置。这种“物理约束 数据融合”的能力才是 PINN 不可替代的地方。最后再分享一点个人体会我跑通这个例子后最大的收获不是把一根梁算准了而是理解了“物理约束”在深度学习里到底意味着什么。传统神经网络靠大量标签数据学习规律数据不够就过拟合PINN 把控制方程写进损失函数相当于给网络装了一套“物理逻辑”即便没有标签数据也能收敛到符合力学规律的解。那阵子我改代码改到最头疼的就是四阶导数的数值稳定性。后来总结一句话坐标必须归一化、导数必须带长度系数、精度必须用 float64。这三个条件同时满足一维梁问题基本就不会出大问题。之后再从等截面梁扩展到变截面梁只需要把 EI 改成 x 的函数从均布载荷扩展到非均匀载荷只需要把 q 改成 x 的函数。整个框架完全不用动。如果你正准备拿 PINN 做你的第一个力学问题我建议就从悬臂梁开始别嫌它简单。这个例子里包含的四个边界条件、四阶导数、量纲换算、损失权重调节是所有后续复杂问题的公共底座。把这些摸透了再去碰高维问题、反问题、多物理场耦合会轻松很多。