2. 案例2：数据驱动代理模型拟合阻尼摆轨迹

任务拆解：从轨迹数据到多步预测

阻尼摆的连续动力学为 θ̈ + c·θ̇ + (g/L)·sinθ = 0。本例不直接对微分方程做残差约束，而是把仿真轨迹切成相邻状态对 (s_t, s_{t+1})，其中状态 s = [θ, ω]，让多层感知机学习单步状态转移映射 f_φ: s_t → s_{t+1}，再以自回归方式滚动多步，观察预测轨迹相对真实轨迹的漂移。

关键设计点有三个：一是数据生成必须用固定步长积分并保存完整状态序列；二是网络输入输出都做标准化，否则角度与角速度量级差异会让训练不稳定；三是评估必须区分单步误差与多步滚动误差，后者才是代理模型能否替代仿真的判据。

数据生成与样本构造

用四阶 Runge–Kutta 以 Δt=0.02 s 积分阻尼摆方程，初始条件取 θ₀=1.2 rad, ω₀=0，共积分 2000 步。把相邻两步组成样本对，得到 1999 个可用样本。标准化参数只用训练集统计量计算，避免验证集信息泄漏。

import numpy as np

g, L, c, dt = 9.81, 1.0, 0.15, 0.02

def deriv(s):
    theta, omega = s
    return np.array([omega, -c * omega - (g / L) * np.sin(theta)])

def rk4(s):
    k1 = deriv(s)
    k2 = deriv(s + 0.5 * dt * k1)
    k3 = deriv(s + 0.5 * dt * k2)
    k4 = deriv(s + dt * k3)
    return s + dt / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4)

states = np.zeros((2000, 2))
states[0] = [1.2, 0.0]
for i in range(1, 2000):
    states[i] = rk4(states[i - 1])

X = states[:-1]          # s_t
Y = states[1:]           # s_t+1

n_train = 1600
X_tr, Y_tr = X[:n_train], Y[:n_train]
X_va, Y_va = X[n_train:], Y[n_train:]

mu_x, sd_x = X_tr.mean(0), X_tr.std(0)
mu_y, sd_y = Y_tr.mean(0), Y_tr.std(0)

def norm(a, mu, sd):
    return (a - mu) / sd

X_tr_n = norm(X_tr, mu_x, sd_x)
Y_tr_n = norm(Y_tr, mu_y, sd_y)
X_va_n = norm(X_va, mu_x, sd_x)
Y_va_n = norm(Y_va, mu_y, sd_y)

网络与训练

用 PyTorch 搭建可配置深度的 MLP，隐藏层宽度固定 64，激活 ReLU，输出层不加激活。训练用 Adam，学习率 1e-3，批量 64，训练 300 轮，以验证集单步 MSE 选最优权重。为对比样本量影响，分别用前 200、800、1600 个训练样本重训同一结构。

import torch
import torch.nn as nn

def build_mlp(depth, width=64):
    layers, in_dim = [], 2
    for _ in range(depth):
        layers += [nn.Linear(in_dim, width), nn.ReLU()]
        in_dim = width
    layers.append(nn.Linear(in_dim, 2))
    return nn.Sequential(*layers)

def train_mlp(depth, n_samples, epochs=300, seed=0):
    torch.manual_seed(seed)
    model = build_mlp(depth)
    opt = torch.optim.Adam(model.parameters(), lr=1e-3)
    loss_fn = nn.MSELoss()

    xb = torch.tensor(X_tr_n[:n_samples], dtype=torch.float32)
    yb = torch.tensor(Y_tr_n[:n_samples], dtype=torch.float32)
    xv = torch.tensor(X_va_n, dtype=torch.float32)
    yv = torch.tensor(Y_va_n, dtype=torch.float32)

    best, best_state = float("inf"), None
    for _ in range(epochs):
        model.train()
        perm = torch.randperm(n_samples)
        for i in range(0, n_samples, 64):
            idx = perm[i:i + 64]
            opt.zero_grad()
            loss = loss_fn(model(xb[idx]), yb[idx])
            loss.backward()
            opt.step()
        model.eval()
        with torch.no_grad():
            v = loss_fn(model(xv), yv).item()
        if v < best:
            best, best_state = v, {k: t.clone() for k, t in model.state_dict().items()}
    model.load_state_dict(best_state)
    return model, best

多步滚动预测与漂移评估

滚动预测时，从真实初始状态出发，每一步把网络输出反标准化后作为下一步输入，连续滚动 500 步。漂移用两个指标衡量：状态欧氏误差 ‖ŝ_t − s_t‖ 的累积曲线，以及角度相位偏差。注意滚动输入始终是模型自己的预测，误差会逐步放大，这正是代理模型与单步回归的本质区别。

def rollout(model, s0, steps=500):
    model.eval()
    s = s0.copy()
    traj = [s.copy()]
    with torch.no_grad():
        for _ in range(steps):
            sn = norm(s, mu_x, sd_x)
            pred = model(torch.tensor(sn, dtype=torch.float32)).numpy()
            s = pred * sd_y + mu_y
            traj.append(s.copy())
    return np.array(traj)

true_traj = states[:501]
results = {}
for depth in (1, 2, 3):
    for n in (200, 800, 1600):
        model, val_mse = train_mlp(depth, n)
        pred = rollout(model, states[0])
        err = np.linalg.norm(pred - true_traj, axis=1)
        results[(depth, n)] = (val_mse, err[-1], err.max())

for k, v in sorted(results.items()):
    print(f"depth={k[0]}, N={k[1]}, val_mse={v[0]:.4e}, "
          f"final_err={v[1]:.4f}, max_err={v[2]:.4f}")

预期现象（示例运行下的定性趋势，具体数值随随机种子浮动）：1 层网络在 N=200 时最终误差明显偏大，滚动约 150 步后相位开始偏离；3 层网络在 N=1600 时误差增长最慢。这说明深度增加提升的是映射表达能力，样本量增加抑制的是过拟合导致的滚动发散，两者作用不同，不能互相替代。

结果判读与对比要点

评估时不要只看单步 MSE 排名，应把滚动误差曲线与相位偏差一起看：单步误差接近的两组配置，滚动 500 步后可能相差数倍，这正是本例要暴露的核心问题。