基于物理信息神经网络（PINN）的动力学系统模拟代码异常排查

问题概述 #

尝试构建PINN模拟阻尼振子动力学系统：网络以上一时刻位置y、速度y'、外力f为输入，输出下一时刻的y与y'。训练中尝试了有限差分/中心差分/5点模板近似导数，以及分阶段提升损失复杂度的策略，但模拟效果始终未达预期。

核心问题分析 #

1. 训练数据采样完全错误 #

训练循环中使用随机生成的y和dydt作为输入，与真实外力序列完全不匹配，导致网络无法学习到真实的动力学约束，损失计算脱离物理系统本身。

2. 损失函数逻辑矛盾且偏离物理方程 #

• 第一阶段损失中，同时用差分近似速度(predictions[:,0]-input[:,0])/delta_t，又强制该值与网络输出的速度predictions[:,1]一致，逻辑冗余且违背物理意义（速度本身就是位置的导数，应直接用物理方程约束，而非强制两种近似相等）。
• 后续阶段的中心差分、5点模板损失，通过多步预测的差分间接约束物理方程，引入了累积误差，未紧扣原ODE核心约束：d²y/dt² = f(t) - y' -4y。

3. 训练调度与流程不匹配 #

• 学习率调度器MultiStepLR的里程碑设置[0.75*num_epochs,4*num_epochs]与总训练轮次5.5*num_epochs不兼容，第二个里程碑超出有效范围，调度逻辑失效。
• 分阶段切换损失函数的时机未与学习率调整对齐，导致训练稳定性差。

4. 测试阶段初始条件错误 #

训练初始状态为S0=[0,0.5]，但测试预测时误用y_pred=[0], dy_dt_pred=[0.2]，起点偏离真实状态，直接导致后续预测偏差。

5. 网络结构缺乏动力学适配性 #

简单MLP未加入残差连接等适配序列预测的结构，也未针对动力学系统的连续输出特性选择合适的激活函数。

具体解决办法 #

1. 修正训练数据采样 #

替换随机生成的状态输入为真实训练数据切片，保证状态与外力序列匹配：

# 替换训练循环内的y、dydt生成代码
start_idx = batch_index * batch_size
end_idx = start_idx + batch_size
y = y_train[start_idx:end_idx]
dydt = dydt_train[start_idx:end_idx]
force = force_train[start_idx:end_idx + num_forward_pass]

2. 重构损失函数，紧扣物理方程 #

基于原ODEd²y/dt² = f(t) - v -4x构建损失，优先用离散积分格式直接约束：

基础阶段（欧拉法约束）

# 用欧拉法近似加速度：a ≈ (v_next - v)/delta_t
a_pred = (predictions[:,1] - input[:,1]) / delta_t
# 物理方程要求的真实加速度
a_true = force[1:-(num_forward_pass-1)] - input[:,1] - 4*input[:,0]
# 物理约束损失 + 状态预测损失（提升拟合精度）
state_true = torch.stack((y_train[start_idx+1:end_idx+1], dydt_train[start_idx+1:end_idx+1]), dim=1).to(device)
loss = cost_function(a_pred, a_true) + 0.1 * cost_function(predictions, state_true)

进阶阶段（中点法/更高精度格式）

可替换为中点法近似加速度，提升约束精度：

# 中点法：利用两步预测的速度差计算加速度
a_pred = (predictions_1[:,1] - input[:,1]) / (2*delta_t)
a_true = force[1:-(num_forward_pass-1)] - (predictions_0[:,1]+input[:,1])/2 - 4*(predictions_0[:,0]+input[:,0])/2

3. 调整训练流程与调度 #

• 重新设置学习率调度器，匹配总训练轮次：

scheduler = torch.optim.lr_scheduler.MultiStepLR(optimizer, milestones=[int(total_epochs*0.5), int(total_epochs*0.8)], gamma=0.1)

• 分阶段训练时，先训练状态预测损失，再逐步加入物理约束损失，而非直接切换复杂差分模板。

4. 修正测试阶段配置 #

• 使用测试函数对应的初始条件：

y_pred = [-0.2]
dy_dt_pred = [0.7]

• 确保测试时间序列t与test_function的有效范围（如35<t<75、100<t<120）匹配。

5. 优化网络结构 #

加入残差连接提升序列拟合能力，替换激活函数适配连续输出：

class ResidualBlock(nn.Module):
    def __init__(self, dim):
        super().__init__()
        self.fc1 = nn.Linear(dim, dim)
        self.tanh = nn.Tanh()
        self.fc2 = nn.Linear(dim, dim)
    def forward(self, x):
        residual = x
        out = self.fc1(x)
        out = self.tanh(out)
        out = self.fc2(out)
        out += residual
        return self.tanh(out)

neural_solution = nn.Sequential(
    nn.Linear(3, 64),
    nn.Tanh(),
    ResidualBlock(64),
    ResidualBlock(64),
    nn.Linear(64, 2),
).double().to(device)

6. 提升训练稳定性 #

• 降低初始学习率至0.001，避免初期震荡；
• 加入梯度裁剪防止梯度爆炸：

torch.nn.utils.clip_grad_norm_(neural_solution.parameters(), max_norm=1.0)

• 增加训练轮次至10*num_epochs，保证网络充分收敛。

内容的提问来源于stack exchange，提问作者Mihir Kotak