Matlab实现PINN多输入单输出时序预测模型

开篇先聊一个很实际的问题：我们在工业场景里做预测，手里拿到的往往是好几个传感器通道的数据，想预测的目标却只有一个。比如根据压缩机进口温度、出口压力、转速、振动位移，去预测轴承的剩余寿命；或者根据生产线的温度、湿度、电流、张力，去预测一卷料膜最后的厚度偏差。

这类问题套到深度学习里，就是典型的多输入单输出时序预测。常规做法大家都熟：LSTM、Transformer、TCN，把历史窗口喂进去，让网络自己去学特征。但这类纯数据驱动模型有个躲不掉的硬伤——数据不够的时候，学出来的规律很容易跑偏；而且模型只认数据里的统计相关性，完全不认物理规律。

PINN（物理信息神经网络）这两年火，核心思路其实非常朴素：既然你预测的对象是某个物理系统的输出，那这个系统的控制方程（微分方程）就是现成的约束。把方程残差当成损失函数的一部分，逼着神经网络的输出同时满足“数据长得像”和“物理上说得通”两个条件。

这篇文章就围绕“用Matlab实现一个多输入单输出的PINN时序预测模型”展开。我会直接给出可运行的思路、核心代码骨架、坑位清单和调参经验。适合有Matlab基础、想在自己的项目里把物理先验塞进预测模型里的朋友参考。

先说清楚，Matlab做PINN相比Python有个优势——Deep Learning Toolbox的自动微分是封装好的，dlgradient和dlnetwork这对组合用起来非常顺手，不需要自己写反向传播。但同时，网上关于“Matlab版PINN”的学习资料明显比Python少一大截，很多Python教程里的细节拿到Matlab环境下要重新踩一遍。这篇文章就把那些坑提前帮你填上。

1. 为什么多变量时序预测需要物理信息

1.1 纯数据驱动模型的硬伤

先举一个我实际调试过的例子。某设备的关键参数预测任务，输入是6个传感器变量、每步50个时间点，输出只有1个温度值。用LSTM训练，训练集上loss降得挺漂亮，验证集上平时也还算稳定。但一旦设备工况出现波动，比如进料温度短时间内跳了20%，LSTM的预测值就开始“飘”，误差直接放大三四倍。

这种情况不是个例。纯数据驱动模型学到的本质是训练集范围内的统计映射关系，它根本不知道背后的热传导方程是怎么约束温度场演化的。你给它的输入模式一旦跳出训练数据覆盖的区域，它就只能靠“猜”。

PINN解决问题的角度完全不同。它不要求你放弃神经网络，而是要求你在损失函数里“植入”控制方程。模型输出不仅要拟合历史数据，它的导数组合还要满足物理方程的残差趋近于零。换句话说，网络在学习时被物理定律“扶着走”，即使某段数据区间样本稀疏，物理方程也能给梯度提供一个合理的约束方向。

1.2 物理约束的本质：把微分方程变成损失项

假设一个简单的一维热传导过程：

[ \frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2} + Q(t) ]

输出温度T是网络的预测值，α是热扩散系数，Q是外源项（可以是输入变量之一）。

PINN的做法是：让神经网络输出一个函数 (\hat{T}(t, x))，然后把这个函数代进方程，算一个残差：

[ r(t,x) = \frac{\partial \hat{T}}{\partial t} - \alpha \frac{\partial^2 \hat{T}}{\partial x^2} - Q(t) ]

如果网络输出完全符合物理规律，那 (r(t,x)) 应该到处都接近0。于是损失函数变成：

[ L = L_{data} + \lambda L_{physics} ]

其中 (L_{data}) 是网络在真实样本上的拟合误差（常用MSE），(L_{physics}) 是方程残差的均方值，λ是物理项的权重系数。

这里有个特别关键的点：方程里的时间导数和空间导数，在Matlab里不是靠数值差分算的，而是靠自动微分（AD）算的。dlgradient对网络输出的某个输入维度求导，就是PINN的标准动作。

1.3 时序预测里PINN到底约束的是什么

时序预测和静态回归有个区别：时间方向天然有因果关系。LSTM这类循环结构内部其实是在学习一个隐式的状态转移函数，但这个转移函数对不对，模型自己不知道。PINN在时序问题里做的事，就是给这个隐式转移加一个显式的物理规则。

比如你要预测一个机械振动系统的位移响应，系统的控制方程就是一个二阶常微分方程：

[ m\ddot{x} + c\dot{x} + kx = F(t) ]

网络输入是历史位移和激振力，输出是未来位移。PINN额外要求：网络输出 (\hat{x}(t)) 在时间方向上求二阶导，然后组合成 ((m\ddot{\hat{x}} + c\dot{\hat{x}} + k\hat{x} - F))，这个残差要趋近于0。这相当于把“系统的质量、阻尼、刚度参数”也变成了隐式约束，模型输出就不会偏离物理事实太远。

明白这一点，就理解为什么标题里强调“多输入单输出”了：多变量输入主要给网络提供外源驱动信息和历史状态信息，物理方程管的是“输出变量如何随时间演化”，两者不矛盾，反而互补。

2. 从问题建模到数据构造

2.1 多输入单输出的问题定义

先给一个明确的数学定义，方便后面代码对照。假设：

• 输入变量维度 (m=4)，分别为 (u_1(t), u_2(t), u_3(t), u_4(t))
• 输出变量 (y(t))，维数1
• 每个训练样本是一个历史窗口：从 (t-L+1) 到 (t) 的输入序列，形状是 ([L, m])
• 预测目标：未来 (H) 步的 (y) 值，本文先只讨论单步预测 (H=1)

数据矩阵形式如下表所示：

滑窗切分后，每个样本 (X_i) 的形状是 ([L, m])，每个标签 (Y_i) 是一个标量。这就是标准的监督学习格式。

2.2 数据归一化：千万别忽略的物理一致性

很多初学者在PINN里翻车，第一关就死在归一化上。

物理方程里的系数，比如热扩散系数α、刚度k、阻尼c，都是在原始物理量纲下定义的。如果把输入输出数据全部归一化到[0,1]或[-1,1]，那么方程殘差的计算也必须同步换算，否则系数对不上，物理约束就完全失效。

我的做法分两步：

1. 数据归一化确实要做，对网络收敛有好处。但要把归一化所用的均值和标准差（或min/max）保存下来。
2. 物理方程残差计算直接基于网络的原始输出做自动微分，然后再把导数值反归一化回物理量纲，最后代入方程。

实际操作中还有一个更取巧的办法：直接在归一化空间里重新缩放物理方程的系数。因为 (x' = (x - \mu)/\sigma)，所以 (dx'/dt = (1/\sigma)(dx/dt))，二阶导以此类推。把缩放后的系数重新代入方程，这样方程残差的计算全程都在归一化空间完成，网络内部不需要任何额外处理。

这里给一个通用换算公式。假设原方程：

[ a_2\ddot{x} + a_1\dot{x} + a_0 x = F ]

定义缩放 (x = x_{std} x' + x_{mean})，(F = F_{std} F' + F_{mean})，那么归一化空间的方程变为：

[ (a_2 / x_{std}) \ddot{x}' + (a_1 / x_{std}) \dot{x}' + a_0 x' = (F_{std}/x_{std}) F' + (F_{mean} - a_0 x_{mean})/x_{std} ]

看着复杂，但代码里其实只是一个系数向量预先算好，跑起来很快。

2.3 训练集/验证集的划分技巧

时序预测的数据划分不能像静态数据集那样随机打乱。原因很简单：时间窗口之间有重叠，随机打乱会把未来信息泄漏到训练集里，模型验证结果会虚假偏高。

推荐的做法是：按时间顺序划分，前70%做训练，中间15%做验证（用于调超参），后15%做测试（最终评估）。如果你的实验数据来自多段不同工况，更合理的方式是以“工况段”为单位划分，保证同一个工况的数据不会同时出现在训练集和测试集。

这个细节对带物理约束的模型尤其重要——因为PINN本身对训练数据的依赖就比纯数据驱动小，如果数据划分再做不好，很难判断提升到底是来自物理约束还是来自数据泄漏。

3. 网络结构与Matlab代码实现

3.1 选用dlnetwork搭建全连接网络

Matlab的dlnetwork是深度学习工具箱里推荐的方式，支持自定义损失函数和自动微分。配合dlgradient，可以在一次前向传播里同时算出数据损失和物理损失，并完成梯度回传。

PINN里的主力网络其实不需要太花哨。全连接网络（MLP）配合适当的激活函数，理论上已经能逼近很复杂的非线性函数。我用的是三层隐藏层、每层32个神经元的结构，激活函数全部用tanh。

为什么用tanh？因为PINN需要对输出求高阶导数，而ReLU在零点不可导、二阶导恒等于0，这会让方程里带二阶导的残差项信息丢失。tanh光滑可导，二阶导也有真实的变化信息，是当前PINN实践里最稳妥的选择。

3.2 核心代码骨架

下面这份代码是完整可运行的骨架，省略了数据加载部分，重点展示PINN的核心训练循环。

% 搭建网络 inputSize = L * m; % 将[L, m]展平成向量作为输入 layers = [ featureInputLayer(inputSize, 'Normalization', 'none') fullyConnectedLayer(32) tanhLayer fullyConnectedLayer(32) tanhLayer fullyConnectedLayer(32) tanhLayer fullyConnectedLayer(1) ]; net = dlnetwork(layers); % 训练参数 numEpochs = 300; miniBatchSize = 64; initialLearnRate = 1e-3; lambdaPhysics = 0.1; % 物理项权重 averageGrad = []; averageSqGrad = []; gradDecay = 0.9; sqGradDecay = 0.999; % 展开训练数据：每个样本X(i)展平成1D向量 % XTrain: [inputSize, numSamples] dlarray % YTrain: [1, numSamples] dlarray % 额外附带物理方程需要的原始输入变量（比如外源项F） for epoch = 1:numEpochs iteration = 0; for i = 1:numIterationsPerEpoch % 取batch idx = randi(numSamples, miniBatchSize, 1); XBatch = XTrain(:, idx); YBatch = YTrain(:, idx); % 物理方程所需的额外变量也取对应batch [loss, grad] = dlfeval(@modelLoss, net, XBatch, YBatch, ... U1Batch, U2Batch, lambdaPhysics, dt); % Adam更新 [net, averageGrad, averageSqGrad] = adamupdate(net, grad, ... averageGrad, averageSqGrad, iteration, initialLearnRate, ... gradDecay, sqGradDecay); end % 每个epoch打印loss end

3.3 损失函数的关键实现：modelLoss函数

modelLoss是整段代码的核心。我先按照一个简单的物理系统来写：预测对象是一个受阻尼振动的位移响应，输入变量包括历史位移、历史速度、外激励力。这个系统可以看成很多物理系统的简化抽象。

function [loss, grad] = modelLoss(net, XBatch, YBatch, FextBatch, lambdaPhysics, dt) % 前向传播 YPred = forward(net, XBatch); % 数据损失 dataLoss = mse(YPred, YBatch); % 物理损失：计算预测输出对时间的一阶导数 % 这里需要把时间维度从输入中取出来 % 假设XBatch的最后一列是时间t % 对YPred求关于输入最后一列（时间）的梯度 dYdT = dlgradient(sum(YPred), XBatch(end, :)); % 构建物理系统的微分方程残差 % 假设系统： m*y'' + c*y' + k*y = Fext % 参数设为 m=2, c=0.5, k=3（取决于具体系统） m = 2; c = 0.5; k = 3; % 二阶导 d2YdT2 = dlgradient(sum(dYdT), XBatch(end, :)); % 物理约束：m*y'' + c*y' + k*y - F = 0 physicsResidual = m * d2YdT2 + c * dYdT + k * YPred - FextBatch; physicsLoss = mean(physicsResidual.^2); % 总损失 loss = dataLoss + lambdaPhysics * physicsLoss; % 自动微分获得梯度 grad = dlgradient(loss, net.Learnables); end

这段代码里有三个容易出错的细节要重点解释一下：

第一个是dlgradient(sum(YPred), XBatch(end,:))这个写法。dlgradient要求第一个参数是标量，所以要对YPred求和再求梯度。因为YPred里每个元素只依赖对应样本的输入，求和不会破坏梯度计算逻辑，这个技巧在Matlab的PINN实现里很常用。

第二个是二阶导的处理。d2YdT2 = dlgradient(sum(dYdT), XBatch(end,:))，这个写法有一个前置条件：dYdT必须被保留为dlarray的跟踪状态，不能经过普通的数值数组转换。初学者最容易犯的错误是在求完一阶导后把它extractdata了，导致二阶导直接报错或梯度为0。

第三个是物理参数的设置。m=2,c=0.5,k=3是示例值。实际工程项目里，这些参数可能来自设备铭牌、设计文档或系统辨识结果。参数给错了，PINN的物理损失不仅帮不上忙，还会把模型往错误方向推。我在实际项目里见过有人把阻尼比写大了10倍，结果预测曲线振荡衰减得特别快，验证集误差比纯LSTM还高。

3.4 多步预测是如何嵌进这个框架的

如果标题需求是单步预测，上面代码已经够了。但如果你要的是多步预测（比如预测未来5个时刻），代码要做两处调整：

1. 网络输出维度从1改为H，即输出未来H个时刻的值。
2. 物理损失需要在H个输出点上分别计算方程残差，然后对所有残差求平均。

输出维度改法：把最后一层fullyConnectedLayer(1)改成fullyConnectedLayer(H)，标签也从标量改成向量。

物理损失的改法：

% 对每个预测时刻单独计算物理残差 physicsLoss = 0; for h = 1:H dYdT = dlgradient(sum(YPred(h, :)), XBatch(end, :)); % ... 计算残差 physicsResidual = m * d2YdT2 + c * dYdT + k * YPred(h, :) - FextBatch(h, :); physicsLoss = physicsLoss + mean(physicsResidual.^2); end physicsLoss = physicsLoss / H;

这里有个细节：外激励力Fext在处理多步预测时也是需要未来H个时刻的序列。如果你使用的是历史时刻的观测外激励，那么预测的就是“给定外激励序列求响应”，这个在工程上是合理的设定。

4. 训练策略与调参细节

4.1 为什么默认先试Adam

PINN的损失函数实际上是一个多目标优化问题：数据拟合和物理方程两边要同时满足。这种非凸优化问题对优化器的选择很敏感。L-BFGS这类二阶方法收敛快但容易过拟合到局部极值，而且Matlab的fmincon配合dlarray的自动微分衔接起来比较别扭。

我通常先用Adam跑300到500个epoch，把网络拉到最优附近，再切换到L-BFGS做精调。切换时把Adam学到的参数作为L-BFGS的初值，这样能兼顾稳定性和精度。

如果你不想搞这么麻烦，只留着Adam也能出结果，我的经验是验证集精度会稍差一点，但整体趋势不会错。

4.2 物理权重lambda的调整思路

lambda是最值得花时间调的超参数。它控制物理约束的强度，调不好会导致两个问题：

lambda太大，网络输出会被物理方程“绑架”，变的非常平滑，但拟合不了数据里的高频细节；lambda太小，物理约束形同虚设，模型退化成普通的MLP回归。

我的调整策略是从大到小搜索，先设lambda=1跑一个短训练，观察验证集loss。如果dataLoss能正常下降，说明物理约束没有压制数据拟合；如果dataLoss卡在高位下不去，就把lambda缩小10倍再试。

另一个细节是，lambda可以做成随训练过程变化的。前200个epoch用小lambda（例如0.01）让网络先把数据基本模式学到，后面再逐渐放大到0.1或0.5，让物理约束做精调。这种“课程式”的训练方式在PINN里效果很好。

4.3 学习率与批大小

PINN的损失函数梯度包含高阶导数的反传，梯度量级往往比普通深度学习大，学习率太大会直接发散。我的初始值一般用1e-3，如果训练曲线出现震荡就降到3e-4或1e-4。

批大小方面，物理损失的计算和样本是独立的，不需要特别大的batch。实测miniBatchSize=64在大多数问题上已经够了，大batch反而会让物理损失下降变慢。

4.4 收敛判据与早停

PINN里不能只看总loss，要分开看dataLoss和physicsLoss的变化曲线。一个健康的训练过程是：dataLoss先快速下降，physicsLoss慢慢跟上。如果physicsLoss一直不降，说明物理约束和网络表达之间有冲突，优先检查归一化换算和导数计算。

我习惯保存验证集loss最低的模型，而不是最后一个epoch的模型。因为PINN在训练后期很容易出现过拟合到物理残差的现象——训练集上总loss很好看，验证集上一塌糊涂。

5. 消融实验与踩坑记录

5.1 有无物理约束的对比效果

我用一组公开的阻尼振动数据做了消融实验，输入包含历史位移、速度、外激励三个变量，预测下一时刻的位移。训练样本只给了500个，模拟小样本场景。

从结果可以清楚看到，lambda=0.1时验证集和测试集表现最好，lambda=1.0时训练集误差反而更高，因为物理约束太强压制了数据拟合能力。

这个表给我们的启发是：PINN不是加得越猛越好，存在一个“甜点区间”。在样本量小、噪声大的场景下，这个甜点区间的收益非常明显；但如果数据足够多、噪声也够低，纯MLP和PINN的差距会缩小，这时物理约束的主要价值就体现在外推能力上。

5.2 坑位一：dlgradient二阶导计算失败

第一次实现二阶导时，我遇到了一个典型的报错：“Value must be a dlarray”。排查了半天发现是一阶导被我用了extractdata取出数值再求二阶导。在Matlab里，dlgradient的结果本身仍然是dlarray，直接对该结果再次调用dlgradient是可以的，但一定要保持它在计算图内。

解决办法其实很简单：求一阶导之后不要做任何取数操作，直接传给下一个dlgradient。如果你需要在物理损失里使用一阶导的数值，先用它做完所有计算，最后统一提取。

5.3 坑位二：物理参数的量纲混乱

这个问题在3.3节的归一化部分提过。我实际犯过的错误是：在归一化空间里用了原始物理参数，导致方程残差比别人大了好几个数量级。那时候物理损失小不下去，网络输出明显偏离数据，我还一直在调lambda，完全没意识到是系数问题。

解决办法是上面给的那个统一换算公式。写代码的时候把系数换算写在注释旁边，方便后续复查。

5.4 坑位三：时间步长dt的设置

PINN在时间方向上对输出求导数时，dt的物理单位要和数据采样周期一致。如果你的数据是每0.01秒采样一次，但网络输入里时间单位用了1秒，那求出来的导数就放大了100倍，方程残差也完全不是那么回事。

我在代码里的处理方式是：在网络输入里直接加一个时间通道，值为(0:L-1)*dt。这样dlgradient对时间通道求导时，自动把采样周期考虑进去了。不要在多个变量里混用不同的时间基准。

5.5 后续可扩展的方向

写到这里，这个框架其实已经覆盖了PINN做多输入单输出时序预测的核心链路。如果你的项目有更高要求，可以从这几个方向继续做：

• 把单步预测改成滚动多步预测，训练时让网络接收上一步的预测输出作为下一步的部分输入，但要注意误差累计会导致训练不稳定。
• 在物理损失里加自适应权重，利用Gradient Norm的统计信息动态调整lambda。
• 把全连接网络替换成LSTM结构，在dlnetwork里支持自定义的LSTM层，这样既能保留物理约束，又能利用循环结构提取时序依赖。

我个人在实际项目里的体会是：PINN最值得用的场景不是大样本高精度预测，而是小样本、有噪声、需要外推的场景。物理约束就像给了模型一副“物理直觉”的眼镜，帮它在数据稀疏的地方也能保持合理的输出行为。Matlab环境下这套实现虽然资料少，但核心链路其实是简明的——构造网络、算导数、组损失、迭代优化，把这四步跑通之后，迁移到自己的物理系统上只是换一个方程的事。