简介本资源面向从事科学计算与深度学习交叉研究的学习者提供基于PINN物理信息网络的四套Python实现方案分别覆盖离散时间识别、离散时间推理、连续时间识别与连续时间推理四类任务适合需要复现物理约束神经网络、验证时间序列建模效果的科研人员与高年级学生。压缩包共559个文件约475.95MB以503个txt结果记录、13个py源码、12个csv误差表、9个eps与9个pdf图表文档、7个mat数据文件为主另含少量pyc、yml、html、md与nb文件便于直接运行与对照分析。内容预览中可见error_table、error_lambda_1_table、error_lambda_2_table等误差统计表可用于评估不同正则权重下的识别与推理精度。目前已有1546人学习下载读者可获取完整源码、实验数据与误差分析表格快速搭建PINN实验流程并对比四类方法的性能差异。1. 从四个代码包说起PINN 做时间域识别到底在识别什么拿到「基于 PINN 物理信息网络进行离散时间识别/离散时间推理/连续时间识别/连续时间推理方法」这个标题时很多人第一反应是这不就是把神经网络塞进微分方程里当正则项吗但真正跑过 pinn 代码的人会告诉你时间域的处理方式才是决定这套方案能不能收敛、能不能外推的核心变量。离散时间识别和连续时间识别表面看只是采样点密度的差别实际上对应两套完全不同的损失构造逻辑和网络输入组织方式。离散时间推理通常把时间步当作独立特征维度喂进去网络学的是「从当前状态跳到下一状态」的映射连续时间识别则要求网络输出对时间的导数参与物理残差计算自动微分必须一路开到时间输入上。这四个代码包放在一起恰好覆盖了 PINN 在时间维度上的两条主轴识别从观测数据反推方程结构或参数和推理给定方程正演未来状态再各自切分成离散与连续两种时间处理范式。如果你正在做动力系统辨识、参数反演或者时序外推又不想纯靠黑箱模型硬拟合这套组合值得花时间跑通。下面我按自己复现时的顺序把四个包拆开讲清楚每个的输入输出、关键参数和最容易翻车的地方。2. 四个代码包的分工与最小复现路径2.1 离散时间识别从采样点反推方程系数离散时间识别包的核心任务是给定一组按固定步长采样的状态观测值反推出控制方程的系数或形式。它不要求网络输出对时间的导数而是把相邻两个时间步的状态拼成输入对让网络学习状态转移关系再通过物理残差约束转移关系符合某个候选方程族。常见做法是把候选方程写成线性组合形式网络同时输出组合系数。import torch import torch.nn as nn class DiscreteIdentifier(nn.Module): def __init__(self, state_dim, hidden64): super().__init__() # 输入是相邻两步状态拼接输出是候选基函数系数 self.net nn.Sequential( nn.Linear(state_dim * 2, hidden), nn.Tanh(), nn.Linear(hidden, hidden), nn.Tanh(), nn.Linear(hidden, state_dim) # 输出状态增量 ) # 可学习的方程系数初始化为小量 self.coeff nn.Parameter(torch.randn(state_dim) * 0.01) def forward(self, x_t, x_next): delta_pred self.net(torch.cat([x_t, x_next], dim-1)) # 物理残差预测增量与系数加权基函数之差 residual delta_pred - self.coeff * x_t return delta_pred, residual这段代码里state_dim是状态变量个数比如 Lorenz 系统就是 3。hidden控制网络容量离散识别任务一般 64 到 128 够用再大容易过拟合采样噪声。coeff是真正要识别的物理量训练时它和网络权重一起更新但学习率通常要单独调小否则系数会震荡。损失函数由两部分组成数据拟合项用 MSE 约束delta_pred逼近真实增量物理残差项用 MSE 约束residual趋近零。实际跑的时候物理残差项的权重从 0.01 开始往上加加到 1.0 左右如果还不收敛多半是候选基函数选错了。2.2 离散时间推理给定初始状态滚动预测离散时间推理包和识别包共用状态转移网络结构但目标不同识别关心系数推理关心多步预测精度。它的训练策略是 teacher forcing 和 free running 交替。teacher forcing 阶段用真实历史状态作为输入free running 阶段用网络自己上一步的输出作为下一步输入两者按比例混合。def rollout(model, x0, steps): preds [x0] x_cur x0 for _ in range(steps): # 推理时没有真实下一步用零增量占位后取网络输出 delta, _ model(x_cur, x_cur) x_cur x_cur delta preds.append(x_cur) return torch.stack(preds, dim1)rollout函数是推理包的核心steps是滚动步数一般设 20 到 50。步数太长梯度会爆炸太短学不到长期行为。训练时每个 batch 先跑 teacher forcing 算一步损失再跑 rollout 算多步累积损失两者加权求和。权重比通常取 1:1如果发现预测轨迹很快发散把 rollout 损失权重提到 2 到 3。注意x_cur delta这步没有除以任何缩放因子所以状态变量的量纲要提前归一化到 [-1, 1] 附近否则不同维度的增量尺度差异会让网络偏向大数值维度。2.3 连续时间识别自动微分贯穿时间轴连续时间识别包要求网络输出对时间输入的导数物理残差直接写成微分方程形式。这意味着时间t必须作为网络输入的一个维度且开启requires_grad。常见做法是把状态和时间一起输入网络输出状态估计然后用torch.autograd.grad对时间求导。def continuous_residual(model, t, x_obs): t t.clone().detach().requires_grad_(True) x_pred model(t) # 对时间求一阶导 dxdt torch.autograd.grad( x_pred, t, grad_outputstorch.ones_like(x_pred), create_graphTrue )[0] # 以 Duffing 方程为例构造残差 d2xdt2 torch.autograd.grad( dxdt, t, grad_outputstorch.ones_like(dxdt), create_graphTrue )[0] residual d2xdt2 0.3 * dxdt - x_pred x_pred ** 3 return residual这里create_graphTrue是关键少了它二阶导无法回传。t必须先 detach 再 requires_grad否则计算图会重复挂载导致内存暴涨。残差表达式里的系数 0.3 和 1.0 是待识别量时要写成nn.Parameter并在优化器里单独分组。连续时间识别对时间采样密度不敏感但对时间区间长度敏感区间太短学不到完整动力学一般至少覆盖两个特征周期。2.4 连续时间推理从任意时刻外推连续时间推理包把训练好的连续时间网络当作可微分的方程代理给定初始条件后用 ODE 求解器积分到目标时刻。它不依赖固定步长可以在任意时间点查询状态。from torchdiffeq import odeint class ODEFunc(nn.Module): def __init__(self, net): super().__init__() self.net net def forward(self, t, x): # 网络直接输出 dx/dt return self.net(t, x) def infer(model, x0, t_span): func ODEFunc(model) # 使用 dopri5 自适应步长求解器 traj odeint(func, x0, t_span, methoddopri5, rtol1e-5, atol1e-6) return trajtorchdiffeq是常见选择dopri5适合非刚性系统刚性系统换dopri5会极慢改用rk4固定步长反而稳。rtol和atol控制积分精度推理阶段可以比训练阶段放宽一个量级来提速。t_span是查询时间点张量不要求等间距这是连续时间推理相对离散推理的最大优势。3. 四个包共用的训练循环与参数分组策略3.1 损失函数的三项拆解与权重调度四个包虽然任务不同但损失函数结构高度一致数据拟合项、物理残差项、初值或边界条件项。数据拟合项用观测点上的 MSE物理残差项在配点上计算初值项只在 t0 附近生效。权重调度上我一般用前 2000 步只开数据拟合让网络先粗略拟合观测再逐步加入物理残差从 0.001 线性升到 1.0。如果一开始就三项全开网络会在物理残差和数据拟合之间反复拉扯表现为损失曲线剧烈震荡血泪经验是宁可多花 2000 步预热。def loss_schedule(step, total_steps): # 物理残差权重从 0 线性升到 1 lambda_phys min(1.0, step / (total_steps * 0.3)) # 初值权重保持恒定 lambda_ic 1.0 return lambda_phys, lambda_ictotal_steps是总训练步数0.3是预热比例可以按任务调到 0.2 到 0.5。lambda_phys升太快会导致残差项主导网络输出退化为平凡解比如全零输出让微分方程残差为零但完全不符合观测。3.2 参数分组网络权重与物理系数分开调物理系数和网络权重的梯度尺度差异很大放同一个优化器里用同一个学习率要么系数不动要么网络震荡。常见做法是分两组网络权重用 Adam学习率 1e-3物理系数用 Adam 但学习率 1e-4或者用 LBFGS 单独优化系数。net_params list(model.net.parameters()) coeff_params [model.coeff] optimizer torch.optim.Adam([ {params: net_params, lr: 1e-3}, {params: coeff_params, lr: 1e-4} ]) # 每 1000 步对系数做一次 LBFGS 精调 if step % 1000 0: lbfgs torch.optim.LBFGS(coeff_params, lr0.1) def closure(): lbfgs.zero_grad() loss compute_loss(model) loss.backward() return loss lbfgs.step(closure)LBFGS 的lr设 0.1 到 0.5步数不用多5 到 10 步就够。注意 LBFGS 的 closure 里必须重新计算损失不能复用之前的计算图。如果系数识别结果在真值附近来回跳把 LBFGS 频率降到每 2000 步一次。3.3 配点采样离散与连续的关键差异离散时间包的配点就是观测时间步本身不需要额外采样。连续时间包需要在时间区间内随机采配点配点数量一般取观测点数的 2 到 5 倍。配点太少物理约束覆盖不全太多则单步计算量上升。我一般用拉丁超立方采样在时间轴上撒点避免均匀采样和观测点重合导致残差项虚低。def sample_collocation(t_min, t_max, n_points): # 拉丁超立方采样保证时间轴覆盖均匀 intervals torch.linspace(t_min, t_max, n_points 1) u torch.rand(n_points) t_col intervals[:-1] u * (intervals[1:] - intervals[:-1]) return t_col.unsqueeze(-1).requires_grad_(True)n_points取 500 到 2000视时间区间长度和动力学复杂度定。requires_grad_(True)必须在采样后立刻加上否则连续时间残差计算会报错。4. 避坑与排查四个包跑不通时先看这几条4.1 损失降不下去但梯度正常现象训练几百步后总损失卡在某个值梯度范数不为零但参数几乎不动。原因通常是物理残差项和数据拟合项量纲差异过大比如状态变量在 1e-3 量级而残差在 1e3 量级Adam 的自适应学习率被大量纲项主导。解决对状态变量做归一化或者给物理残差项乘一个尺度因子让两项在数值上接近。我一般先打印两项损失的初始值比值超过 100 倍就先调尺度。4.2 连续时间识别报「element 0 of tensors does not require grad」现象调用torch.autograd.grad时提示输入不需要梯度。原因通常是时间张量在传入网络前被 detach 了或者网络第一层没有接收时间输入。解决检查t是否在采样后立刻requires_grad_(True)检查网络 forward 是否把t拼进了输入。另一个常见原因是torch.no_grad()上下文没退出训练循环里混用了推理代码。4.3 离散推理 rollout 几步后数值爆炸现象free running 阶段预测轨迹在 5 到 10 步后跳到 1e6 量级。原因有两个一是状态没归一化增量累加后放大二是 rollout 损失权重太高网络为了压低多步损失牺牲了单步精度。解决先确认输入数据在 [-1, 1] 内再把 rollout 损失权重从 1.0 降到 0.3 到 0.5同时把 rollout 步数从 50 降到 20。如果还炸在 rollout 里加梯度裁剪torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0)。4.4 物理系数识别结果对初值敏感现象换一组网络初始化识别出的系数差异超过 20%。原因通常是物理残差项权重升得太快网络还没拟合好观测就被物理约束带偏。解决把预热步数从 2000 提到 5000物理残差权重上限从 1.0 降到 0.5同时用 LBFGS 对系数做多轮精调。如果系数仍然跳考虑给系数加一个 L2 正则约束它不要离初始猜测太远。4.5 ODE 求解器推理速度极慢现象连续时间推理用dopri5跑 1000 个时间点要几分钟。原因是非刚性求解器在刚性系统上步长被压得极小。解决先判断系统刚性如果状态变量变化速率差异超过两个量级换rk4固定步长步长取时间区间长度的 1/500。如果必须用自适应求解器把rtol和atol从 1e-6 放宽到 1e-4速度能提升 5 到 10 倍精度损失在可视化上通常看不出来。5. 进阶技巧用连续时间识别包做参数反演与方程发现四个包里最有延展性的是连续时间识别包因为它把方程结构显式地写进了残差改残差表达式就能切换任务。我一般把它当方程发现工具用先假设一个候选基函数库让网络输出基函数系数再对系数做稀疏化。class SparseIdentifier(nn.Module): def __init__(self, n_basis, hidden64): super().__init__() self.net nn.Sequential( nn.Linear(2, hidden), nn.Tanh(), nn.Linear(hidden, hidden), nn.Tanh(), nn.Linear(hidden, 1) ) # 基函数系数加 L1 正则促稀疏 self.xi nn.Parameter(torch.randn(n_basis) * 0.01) def forward(self, t): return self.net(t) def residual(self, t, x_pred, dxdt, basis): # basis 是基函数值矩阵形状 [n_points, n_basis] coeff basis self.xi return dxdt - coeffn_basis是候选基函数个数比如多项式项、三角函数项、常数项。xi是待发现的系数向量训练时加 L1 正则lambda_l1从 1e-4 开始逐步升到 1e-2。训练结束后xi里接近零的项对应不存在的基函数非零项就是方程结构。这个方法对噪声敏感观测噪声超过 1% 时xi的稀疏性会变差需要先做数据平滑或者增加观测点数。验证识别结果是否可信我习惯用两步第一步把识别出的方程代回 ODE 求解器看预测轨迹和观测轨迹的 RMSE 是否在噪声水平附近第二步换一组不同初始条件的观测数据用同一套超参重新识别看系数是否一致。两步都过了才认为方程结构可信。如果第一步过第二步不过多半是观测数据覆盖的状态空间不够需要补充不同区域的采样。这套方案值不值得投入取决于你的数据条件。离散时间包对采样步长敏感适合观测密集且等间隔的场景连续时间包对采样密度不敏感适合观测稀疏但时间跨度长的场景。四个包加起来代码量不大但每个包的参数分组和损失调度都需要按任务微调没有一套超参能通吃。我自己的习惯是先把离散时间识别跑通确认数据质量和候选方程族没问题再切到连续时间包做精细反演。希望帮到你。本文还有配套的精品资源点击获取