MATLAB实现基于CEEMDAN-KPCA-PINN完全集合经验模态分解与自适应噪声（CEEMDAN）结合核主成分分析（KPCA）和物理信息神经网络（PINN）进行多变量时序预测

请注意此篇内容只是一个项目介绍 更多详细内容可直接联系博主本人 

 或者访问对应标题的完整博客或者文档下载页面（含完整的程序，GUI设计和代码详解）

在多变量时序预测领域，工业过程监控、金融风险预警、电力负荷预测、环境与气象监测等场景中，往往面对的是高维、多源、强非线性且含噪声的时间序列数据。传统统计方法和单一的机器学习模型，在处理这类复杂数据时经常遭遇性能瓶颈：一方面难以充分挖掘变量之间的非线性耦合关系，另一方面也难以有效利用领域内的物理机理与约束，从而导致预测结果在外推时段或异常工况下出现明显失真。随着数据规模的不断增大和系统复杂度的不断提升，仅依赖经验建模或黑盒模型已不能满足对高可靠、高精度预测的需求。

在现实工程中，很多多变量时序数据同时受到多种异质噪声的扰动，例如仪器测量误差、通信传输噪声、环境干扰和工况波动等。这些噪声往往表现为非平稳、非高斯特性，且会在不同时间尺度上叠加在真实信号之上。直接在原始数据上训练预测模型，容易使模型过拟合于噪声模式，而忽略真实动力学行为。因此，对原始多变量时序进行适当的预处理与多尺度分解，分离出不同时间尺度的成分并抑制噪声，是获得高精度预测模型的重要前置步骤。在这一背景下，完全集合经验模态分解与自适应噪声（CEEMDAN）作为一种改进的经验模态分解方法，通过引入自适应噪声并迭代加权集成，能够将复杂信号分解为一系列本征模态函数（IMF）和残余项，较好地缓解传统EMD的模态混叠问题，改善分解稳定性和物理可解释性。

然而，CEEMDAN分解得到的多尺度IMF分量往往数量较多，并且不同IMF之间在频带和能量上存在一定冗余。对于多变量场景，每个变量都会生成多条IMF序列，维度急剧膨胀，直接将所有IMF输入到后续模型会带来维度灾难、参数冗余和计算负担，同时也可能使模型出现过拟合。核主成分分析（KPCA）通过在高维特征空间中进行主成分提取，利用核函数刻画非线性结构，能够在保留主要非线性特征的同时，实现对高维数据的有效降维，从而在一定程度上缓解维度灾难，并提升特征表示的紧凑性和区分度。将KPCA应用于CEEMDAN分解得到的多尺度IMF特征上，可以形成“分解+非线性降维”的预处理组合，为后续预测模型提供更为干净、低维且具有高度判别能力的特征表示。

另一方面，现代深度学习模型（如LSTM、GRU、TCN等）在时序预测方面展示了强大的拟合能力，但这类模型在大多数场合属于纯数据驱动范式，缺乏对系统物理规律、守恒定律或微分方程约束的显式融合。在训练数据有限，或预测时段远离训练数据分布时，这类模型很容易在外推方面失去可靠性，预测结果可能出现明显的物理不一致现象。例如在流体力学、传热传质、电力系统动态中，系统演化通常遵循某种偏微分方程或常微分方程，如果预测结果违反能量守恒、质量守恒等基本物理规律，即使短期误差较小，也难以作为工程决策依据。物理信息神经网络（PINN）通过将物理方程、边界条件、初始条件等以软约束形式融入损失函数中，使神经网络在拟合观测数据的同时，还要满足物理残差尽量为零，从而在一定程度上实现数据驱动与机理模型的融合，提升模型的可解释性与外推能力。

基于上述背景，将CEEMDAN、KPCA和PINN三者有机融合，构建一个面向多变量时序预测的混合模型架构，具有明确的现实意义与理论价值。CEEMDAN负责从复杂多变量时间序列中提取不同时间尺度的本征模态成分，起到多尺度噪声抑制与特征分离的作用；KPCA在此基础上对多尺度、多变量的IMF特征进行非线性降维，提取最具代表性的主成分子空间，降低后续模型的输入维度并提升训练效率；PINN作为最终的预测核心，利用降维后的特征以及已知的物理约束，对系统未来状态进行预测，实现对多变量时序的高精度、物理一致的预测。整体流程能够充分利用时频分解、非线性特征提取与物理约束深度学习三方面的优势，以期在高噪声、高维度、强非线性、机理部分已知的复杂系统中获得更稳健、更可靠的预测结果。

在工程实现层面，MATLAB R2025b 提供了成熟的数值计算、信号处理和深度学习工具箱，为构建CEEMDAN-KPCA-PINN这一复合框架提供了良好的开发环境。可以利用MATLAB高效的矩阵运算、可视化与优化工具，实现从数据导入、预处理、特征工程、模型设计、训练调优到预测评估的完整流程。同时，MATLAB对微分方程求解、符号计算等领域也有较完善的支持，非常适合将物理方程转化为神经网络可计算的残差形式，为PINN模块的构建提供便利条件。综合以上因素，本项目以MATLAB R2025b为平台，实现基于CEEMDAN-KPCA-PINN的多变量时序预测框架，既契合当前智能建模的发展趋势，又具备较强的工程落地潜力。

项目目标与意义

多尺度降噪与特征提取能力提升

首要目标是构建一套能够针对多变量复杂时序数据进行多尺度降噪与有效特征提取的机制。通过引入CEEMDAN分解，将原始多变量信号分解为若干具有不同本征时间尺度的IMF分量及残余项，使得原本高度混叠的频率成分在分解后得到较好的分离。高频IMF通常蕴含测量噪声和短期波动，中低频IMF与残余项则更接近系统的主体趋势与缓变结构。从建模角度来看，利用这类多尺度表示可以针对不同时间尺度设计差异化的建模策略，从而增强模型在趋势预测、短期动态捕捉方面的能力。相比直接在原始序列上训练模型，这种多尺度分解方式有助于减少噪声对模型参数估计的干扰，提高参数收敛稳定性，降低模型对异常点的敏感程度。进一步地，在多变量场景中，每个变量经CEEMDAN分解后的IMF集合可以看成一组高维特征序列，在这些序列中存在大量冗余信息和跨变量耦合关系。通过后续KPCA的非线性降维处理，在保留关键信息的前提下压缩特征空间，使得输入特征更为精炼，有助于提升整体建模性能。因此，多尺度降噪与特征提取不仅是模型前端的数据预处理，更是提升预测精度和鲁棒性的关键环节，具有显著的理论价值和工程意义。

非线性降维与复杂耦合关系刻画

多变量时序数据常常存在高度非线性与强耦合关系，例如多个传感器之间的反馈和相互影响、系统内部不同物理量之间的非线性耦合作用等。如果使用传统线性降维方法，如标准PCA，往往只能捕获线性相关结构，对于本质非线性的流形结构则难以有效表示。本项目在CEEMDAN多尺度分解的基础上，通过核主成分分析构建非线性特征空间，将原始高维IMF特征映射到隐式的高维核空间中，再在该空间中提取主成分，从而实现对复杂非线性结构的捕获。核函数的选择（如高斯核、多项式核等）使特征空间具有更丰富的表达能力，能够在不显式构造高维特征的情况下，间接完成复杂特征变换。在项目目标层面，通过KPCA模块实现的非线性降维，不仅可以显著降低后续PINN模型的输入维度，减少训练参数数量和训练时间，还能在一定程度上消除不同IMF之间的冗余，提高特征的独立性和可解释性。这对于复杂系统中变量众多、特征冗余严重的场景尤为重要，有助于避免维度灾难和过拟合问题，使整体预测框架在面对高维多变量数据时依然保持较好的稳健性和泛化能力。

融合物理约束的高精度多变量预测

在构建预测模型时，目标不仅是获得较低的统计误差，还要尽可能保持预测结果与已知物理机理的一致性。物理信息神经网络通过在损失函数中加入物理方程残差、边界条件和初值约束，使网络在拟合数据的同时尽量满足物理规律。这种设计在数据量有限、不均匀采样或存在观测缺失的情况下尤为重要，可以在一定程度上利用物理知识弥补数据的不足。本项目通过将CEEMDAN-KPCA得到的低维特征输入到PINN中，利用多层全连接网络或带时间嵌入的网络结构，对多变量系统的未来状态进行预测。在训练过程中，一方面通过数据误差项（即预测值与观测值之间的差异）约束模型拟合能力，另一方面通过物理误差项（即导数或偏导数代入物理方程后的残差）约束模型的物理一致性。目标是在保证预测精度的同时，使输出结果在能量守恒、质量守恒、平衡方程等方面具有合理性。这种融合物理约束的预测模式，可以有效提升模型在长时程外推、极端工况或不可观测状态下的可靠性，使其更适合用于实际工程决策与安全评估。

构建可复用的MATLAB智能建模框架

除了具体的预测任务，本项目还希望形成一套结构清晰、模块分明、具备可复用性的MATLAB实现框架，为后续在其他领域或其他数据集上迁移应用提供便利。通过在MATLAB R2025b环境中实现数据导入、CEEMDAN分解、KPCA降维、PINN构建与训练、预测输出与可视化等模块，可以形成一个具有良好扩展性的工程模板。研究人员或工程技术人员可以根据具体领域需求，替换数据接口、更改物理约束形式或调整网络结构，而无需从零开始搭建整套系统。这样的目标具有明显的工程意义：一方面大幅减少新场景建模的开发成本，提高试验与迭代效率；另一方面通过统一的代码架构和接口设计，降低多人协作开发时的沟通成本与维护成本。结合MATLAB的强大可视化能力，可以在同一平台上完成数据探索、模型训练过程监控与结果分析，为多学科交叉团队提供一个直观、可操作的智能时序预测工具箱。总体而言，构建可复用的MATLAB智能建模框架，有助于推动CEEMDAN-KPCA-PINN方法在更广泛工程领域的推广与应用。

项目挑战及解决方案

多尺度分解与参数选择的复杂性

在CEEMDAN分解过程中，噪声强度、集合次数、停止准则等关键参数直接影响分解质量和计算成本。不恰当的参数设置容易引入额外噪声、导致模态混叠或分解结果失真，也会显著增加计算时间，尤其是在多变量长序列场景下。多变量情况下，若对每个变量独立进行CEEMDAN分解，整体计算量呈线性倍增，且不同变量之间的分解模式可能不一致，从而给后续特征融合带来困难。解决这一挑战需要在理论与实践两方面共同考虑。首先，在理论层面，需要充分理解CEEMDAN算法中各参数对分解结果的影响，通过分析IMF的频谱特性、能量分布和稳定性指标，为参数设置提供指导。其次，在实践层面，需要设计一套自动化或半自动化的参数搜索策略。例如，可在小样本或短时窗口上对多个参数组合进行试验，采用评价指标（如重构误差、模态正交性指标等）选择较优配置，并在全局数据上使用。为了降低计算成本，可以利用并行计算或向量化操作加速CEEMDAN实现，避免冗余循环。在MATLAB环境中，可以针对多变量数据，设计批量分解函数，对各变量进行统一管理，并在分解过程中记录IMF特征以便后续KPCA使用。通过参数网格搜索与经验规则的结合，以及合理的代码优化，可以在保证分解质量的前提下，将CEEMDAN的计算时间控制在可接受范围内，为后续模型提供稳定可靠的多尺度特征输入。

非线性降维与核参数调优问题

KPCA的核心在于核函数和核参数的选择，其中最关键的是核类型（如高斯核、多项式核）及其参数（如核宽度）。在多变量时序场景中，CEEMDAN分解产生的IMF特征在不同尺度上具有截然不同的统计特性，若直接采用统一的核参数，可能会导致部分尺度特征在核空间中被压缩过度或放大过度，从而影响主成分提取效果。此外，多变量IMF特征往往具有较高维度，如果核矩阵规模过大，会带来计算与存储压力。解决这一挑战需要在核函数选择、参数调优和降维策略上进行综合设计。一方面，在核函数选择上，可以优先采用在高维空间表现良好的高斯核，通过调节核宽度来平衡局部与整体结构的捕获能力。另一方面，为避免过度依赖经验调参，可以在训练数据上使用交叉验证或基于重构误差的评估策略，自动选择核参数组合。例如，可以在一个候选参数集合上计算KPCA重构误差、主成分累积方差贡献率以及对下游PINN训练误差的影响，从而选出对整体任务最有利的参数配置。对于核矩阵规模问题，可以采取降采样、时间窗口化、或采用部分IMF组合方式来控制输入维度。此外，还可以针对不同时间尺度采用分组KPCA，将高频、中频和低频IMF分别进行核主成分分析，再合并主成分输出，从结构上减少单次降维的特征数量。通过这些策略，可以在保持KPCA非线性表达能力的前提下，兼顾计算效率与稳定性，使其成为CEEMDAN-KPCA-PINN框架中的可靠特征工程模块。

PINN构建与训练稳定性的平衡

物理信息神经网络与传统神经网络相比，在损失函数中增加了物理残差项、边界条件和初值约束项，使得损失结构更加复杂，这既是优势也是挑战。物理残差项通常涉及对网络输出关于时间或空间的导数计算，在MATLAB深度学习框架下，需要利用自动微分机制对网络进行前向传播与梯度计算。如果网络结构不合理或损失项权重设置不当，容易出现训练不收敛、梯度爆炸或梯度消失等问题。此外，物理方程可能具有刚性特征或高阶导数，数值求导的稳定性与精度都会影响PINN的最终表现。为应对这一挑战，首先在网络结构设计上需要控制网络深度与宽度，选用适当的激活函数，如tanh或swish等，以兼顾非线性表达能力与梯度稳定性。对于多变量时序预测任务，还需要设计合理的输入形式，例如将时间、KPCA特征以及物理参数等统一封装为网络输入向量。其次，在损失函数构造方面，需要通过权重系数平衡数据拟合误差与物理残差误差，使二者处于同一数量级，避免某一项主导训练过程。可以采用逐步增权或自适应权重调整策略，在训练初期更关注数据拟合，随后逐渐增加物理约束权重，使网络在已有较好拟合基础上进一步满足物理规律。训练策略方面，可以选择基于mini-batch的梯度下降方法，并采用学习率衰减、梯度裁剪等技巧，加强训练稳定性。在MATLAB R2025b环境中，需要注意dlnetwork的更新方式和dlupdate的使用范围，合理编写自定义训练循环。通过网络结构设计、损失项权重平衡和训练策略优化的组合，可以在较大程度上缓解PINN训练中的不稳定因素，使其在CEEMDAN-KPCA预处理基础上发挥优势，实现对复杂多变量时序的高精度、物理一致预测。

项目模型架构

CEEMDAN多尺度分解模块

模型整体架构的首个核心模块是CEEMDAN多尺度分解模块。该模块针对多变量时序数据中的每一个变量，对其进行完全集合经验模态分解与自适应噪声处理，将原始序列拆解为一组本征模态函数（IMF）和一个残余项。CEEMDAN相对于传统EMD和EEMD的主要改进在于，通过自适应噪声添加和迭代加权，将集合平均引入到每一个分解阶段，以保证每个IMF分量在统计意义上更加稳定，且大幅降低模态混叠现象。具体而言，在每一轮迭代中，将白噪声按一定比例添加到信号上，然后对每个带噪信号进行EMD分解，得到若干IMF候选序列，再通过在集合层面上进行加权平均，形成该轮的最终IMF。对于残余部分，继续重复这一过程，直到残余信号不再满足进一步分解的条件。通过这种方式，可以获得一组频率由高到低排列的IMF序列，分别描述信号的高频振荡、周期性波动以及长期趋势。对于多变量时序，CEEMDAN模块可以采用逐变量分解，或在同步采样情况下采用多通道分解策略，以保持不同变量间在时间尺度上的一致性。该模块的输出，是每个变量对应的一组IMF序列以及一个长期趋势残差，为后续的特征融合与降维提供丰富的多尺度信息来源。在实现层面，通过合理设置噪声强度、集合数和迭代停止阈值，并结合MATLAB的并行计算能力，可以在保证分解质量的基础上控制计算成本，为整体模型提供一个鲁棒的时频分解基础。

KPCA非线性降维与特征融合模块

在完成CEEMDAN分解之后，得到的多尺度IMF序列数量较多，维度高且存在显著冗余。为避免维度灾难并提升模型效率，需要引入KPCA非线性降维模块，对多尺度特征进行融合与压缩。该模块的基本原理是，引入核函数将原始特征映射到高维（甚至无限维）的特征空间，在该空间中执行传统PCA操作，通过求解核特征值问题提取主成分。在实际实现中，仅需构造核矩阵并进行特征分解，而不需显式计算高维映射。对于多变量场景，首先将所有变量的IMF分量在某一时间对齐后堆叠成高维特征向量，再构造核矩阵。核函数常用选择为高斯径向基函数，其参数（如核宽度）可以通过验证集或重构误差评估进行选择。KPCA模块输出的是一组按方差贡献率排序的主成分时间序列，通过保留前若干个主成分即可获得一个低维特征表示。该表示在理论上保留了信号中最主要的非线性结构，减少了噪声干扰，并且显著降低了后续PINN的输入维度。在架构层面，KPCA模块起到了连接CEEMDAN分解与PINN预测之间的桥梁作用，实现从多尺度、高维、多变量IMF空间向紧凑低维特征空间的转化。通过合理设计特征归一化、核参数选择和主成分保留阈值，该模块可以在保证信息保真度的同时大幅提升整体模型的计算效率和抗过拟合能力。

PINN物理约束预测核心模块

在完成多尺度分解与非线性降维后，核心预测任务由PINN物理约束预测模块承担。该模块以时间、KPCA压缩后的主成分特征以及可能的外部控制变量或物理参数作为输入，通过多层神经网络输出目标多变量系统在未来时刻的状态。不同于传统神经网络，PINN在损失函数层面显式引入了物理方程约束。对于具体系统，通常存在某种形式的常微分方程或偏微分方程描述变量之间的动力学关系，例如dx/dt=f(x,u,t,θ)。在PINN中，通过自动微分计算网络输出相对于时间的导数，将其代入物理方程构建残差，再在损失函数中加入该残差项，使网络在训练过程中不仅要拟合历史数据，还要尽量满足这一物理约束。整体损失通常由数据误差项（预测值与观测值差异）、物理残差项（方程残差）、以及边界条件或初值约束项加权组合而成。在网络结构设计方面，PINN模块可以采用多层全连接网络，使用平滑且可微性良好的激活函数，如tanh或swish，以便提高导数计算的稳定性。在多变量预测任务中，输出层需同时预测多个状态变量，可通过增加输出维度实现。训练过程中采用自定义训练循环，结合MATLAB的dlnetwork与自动微分工具，逐步调整网络参数，使损失函数收敛到较小值。通过这种方式，PINN模块在数据驱动预测的同时，引入物理约束提升外推能力，使得预测结果在未观测空间或长期预测时仍具有合理的物理意义。

数据接口与预处理模块

整个模型架构需要一个稳定且灵活的数据接口与预处理模块，用于完成原始数据的读取、清洗、对齐、缺失值处理以及标准化。多变量时序数据可能来自不同的传感器或系统模块，在时间戳、采样频率和幅值范围上存在差异。该模块需对时间戳进行统一对齐，采用插值或重采样策略使各变量在统一时间轴上表达，同时处理异常值与缺测数据。对于少量缺失值，可以采用线性插值、样条插值等方法填补；对于大段缺失，可以考虑剔除或利用物理模型先行估计。在幅值标准化方面，为避免不同变量量纲差异过大影响CEEMDAN与KPCA，需要对各变量进行归一化或标准化处理，例如减均值除标准差，或映射到固定区间。该模块还负责将数据划分为训练集、验证集和测试集，并按照时间顺序划分，避免信息泄漏。对于在线预测场景，还需提供滑动窗口和批量生成接口，以便在PINN训练和预测时高效构造输入输出对。在MATLAB实现中，可以使用结构体或table组织多变量数据及其相关元信息，使得后续模块调用更加直观。通过设计清晰的数据接口和预处理流程，整个CEEMDAN-KPCA-PINN框架能够更便捷地应用于不同数据集和不同工程场景，提升系统的可维护性与可扩展性。

训练与评估模块

模型架构的最后一个关键部分是训练与评估模块，负责整合各个前置模块的输出，执行模型训练过程，并对训练结果进行多角度评估。训练阶段主要包括PINN网络参数初始化、自定义训练循环、损失计算、梯度求解和参数更新。为了提高训练效率和稳定性，需要合理选择优化算法（如Adam或SGDM）、初始学习率、学习率衰减策略和梯度裁剪阈值。同时，应在训练过程中监视训练损失和验证集损失，预防过拟合与欠拟合。评估阶段则需要从多个层面考察模型性能，包括预测误差指标（如MAE、RMSE、MAPE）以及物理残差指标（如物理方程残差的均值和方差），以评估模型在数据拟合和物理一致性方面的综合表现。此外，还需要通过可视化方式直观展示预测结果，如对比预测曲线与真实曲线、展示预测残差随时间的变化等。在CEEMDAN与KPCA模块方面，也可以通过IMF能量分布、主成分方差贡献率等指标评价其分解与降维效果。在MATLAB R2025b环境中，可以利用深度学习工具箱提供的dlnetwork和自定义训练循环功能，实现对PINN的灵活训练，并利用figure和plot等函数完成可视化工作。训练与评估模块还应能够输出训练好的模型参数，以及中间结果（如KPCA投影矩阵、IMF分量），为后续的模型部署或进一步分析提供支持。通过这一模块，整个CEEMDAN-KPCA-PINN架构形成一个闭环，从数据输入到结果评价都有清晰的流程和指标，使模型开发和优化过程更加系统化和可控。

项目模型描述及代码示例

数据读取与预处理示例
clear; % 清空工作区变量，避免历史变量干扰当前实验
clc; % 清空命令行窗口，便于观察本次运行输出
close all; % 关闭所有图形窗口，保证后续绘图环境干净
rawData = readtable(dataFile); % 读取CSV为表格形式，便于按列访问不同变量
timeVec = rawData.Time; % 提取时间列，作为统一时间轴的基础
var1 = rawData.Var1; % 提取第一个变量的时间序列数据
var2 = rawData.Var2; % 提取第二个变量的时间序列数据
var3 = rawData.Var3; % 提取第三个变量的时间序列数据
tStart = timeVec(1); % 记录原始时间序列的起始时间
dt = median(diff(timeVec)); % 估计主采样间隔，用中位数抑制偶发异常差值
var3Uniform = interp1(timeVec,var3,tUniform,'linear','extrap'); % 将变量3插值到统一时间网格，统一时间尺度
nanIdx = isnan(var1Uniform) | isnan(var2Uniform) | isnan(var3Uniform); % 找出任一变量存在NaN的位置，用于统一处理缺失
var1Uniform(nanIdx) = []; % 对应删除变量1在缺失时间点的值，保持长度一致
var2Uniform(nanIdx) = []; % 对应删除变量2在缺失时间点的值，保持数据对齐
mu1 = mean(var1Uniform); % 计算变量1均值，用于中心化处理
sigma1 = std(var1Uniform); % 计算变量1标准差，用于缩放到统一尺度
mu2 = mean(var2Uniform); % 计算变量2均值，为后续标准化提供参数
sigma2 = std(var2Uniform); % 计算变量2标准差，用于拉伸或压缩幅值范围
mu3 = mean(var3Uniform); % 计算变量3均值，以消除直流偏移
sigma3 = std(var3Uniform); % 计算变量3标准差，使各变量量纲相对一致
var1Norm = (var1Uniform - mu1) / sigma1; % 对变量1执行标准化，中心化并缩放
var2Norm = (var2Uniform - mu2) / sigma2; % 对变量2执行标准化，方便与其他变量合并处理
var3Norm = (var3Uniform - mu3) / sigma3; % 对变量3执行标准化，减小量纲差异影响
dataMat = [var1Norm,var2Norm,var3Norm]; % 将三个规范化后的变量按列组合成数据矩阵
nTotal = size(dataMat,1); % 获取总时间步数，用于划分训练验证测试集
nVal = floor(0.15 * nTotal); % 将15%的样本分配给验证集，用于调参与早停判断
idxTrain = 1:nTrain; % 训练集索引范围，从开头到训练样本数
idxVal = (nTrain+1):(nTrain+nVal); % 验证集索引范围，紧随训练集之后
idxTest = (nTrain+nVal+1):nTotal; % 测试集索引范围，为时间序列末尾部分
tTrain = tUniform(idxTrain); % 训练集时间向量，保持与数据同步
tVal = tUniform(idxVal); % 验证集时间向量，用于模型调优过程中指标计算
tTest = tUniform(idxTest); % 测试集时间向量，用于最终预测效果展示
xVal = dataMat(idxVal,:); % 验证集多变量标准化数据矩阵
xTest = dataMat(idxTest,:); % 测试集多变量标准化数据矩阵
CEEMDAN分解示例（单变量）
signal = xTrain(:,1); % 从训练集中选择第一个变量作为示例信号进行CEEMDAN分解
fs = 1 / dt; % 计算采样频率，取决于统一时间步长dt
numEnsemble = 100; % 设置集合次数，用于增加统计稳定性
noiseAmp = 0.2; % 设置噪声幅度比例，控制每次加入噪声的强度
imfSet = []; % 初始化存储IMF的矩阵容器，列对应不同IMF
residual = signal; % 初始化残余信号为原始信号，逐步剥离IMF
for k = 1:maxIMF % 循环提取每一个IMF分量，最多maxIMF个
    imfEnsemble = zeros(N,numEnsemble); % 为当前IMF的集合存储初始化矩阵
        imfEnsemble(:,e) = imfsTmp(:,1); % 将提取到的IMF存入集合矩阵的对应列
    end
    imfMean = mean(imfEnsemble,2); % 对所有集合IMF在列方向取平均，得到该轮的最终IMF
    imfSet = [imfSet,imfMean]; % 将新得到的IMF添加到IMF集合矩阵中
        break; % 终止IMF提取循环，防止无意义的进一步分解
    end
end
trend = residual; % 将最终残余部分视为长期趋势分量
nIMF = size(imfSet,2); % 记录实际得到的IMF数量，可能小于maxIMF
figure; % 新建图形窗口用于展示分解结果
    subplot(nIMF+1,1,i); % 为每个IMF分配一个子图区域，依次排列
    plot(tTrain,imfSet(:,i)); % 绘制第i个IMF随时间变化曲线
    ylabel(['IMF',num2str(i)]); % 标注纵轴为对应IMF序号，便于识别
end
plot(tTrain,trend); % 绘制趋势分量，展示长期变化结构
ylabel('Trend'); % 标注趋势曲线的纵轴含义
xlabel('Time'); % 在最底部子图标注时间轴标签
KPCA特征构建与降维示例
imfFeatures = imfSet; % 将IMF集合视作特征矩阵，每列对应不同时间尺度的分量
featureMat = imfFeatures; % 将单变量IMF特征赋给通用特征矩阵变量，便于后续扩展到多变量
meanFeat = mean(featureMat,1); % 计算各特征维度的均值，用于中心化
Xc = featureMat - meanFeat; % 对特征矩阵做中心化，消除偏移以适应核方法
N = size(Xc,1); % 特征矩阵的样本数，即时间步数
sigmaKernel = 1.0; % 设置高斯核的核宽度参数，用于控制相似性衰减
for i = 1:N % 外层循环遍历第一个样本索引
    for j = 1:N % 内层循环遍历第二个样本索引
        diffVec = Xc(i,:) - Xc(j,:); % 计算样本i和样本j在特征空间的差向量
        dist2 = sum(diffVec.^2); % 计算差向量的平方欧氏距离
        K(i,j) = exp(-dist2/(2*sigmaKernel^2)); % 根据高斯核公式计算核矩阵元素
    end
end
Kc = K - oneN*K - K*oneN + oneN*K*oneN; % 对核矩阵进行中心化，使其对应中心化特征空间
 [EigVec,EigValMat] = eig(Kc); % 对中心化核矩阵进行特征分解，得到特征向量和特征值矩阵
eigVals = diag(EigValMat); % 提取特征值对角线元素，转为向量形式
[eigValsSorted,idxSort] = sort(eigVals,'descend'); % 将特征值按降序排序，记录排序索引
EigVecSorted = EigVec(:,idxSort); % 将特征向量按对应特征值大小重排
varExplained = eigValsSorted / sum(eigValsSorted); % 计算每个核主成分的方差贡献率
cumVar = cumsum(varExplained); % 计算累积方差贡献率，作为选取成分数量依据
numPC = find(cumVar >= 0.95,1); % 找出使累积贡献率达到95%的最小主成分个数
alpha = EigVecSorted(:,1:numPC); % 提取前numPC个特征向量，用作KPCA投影系数
lambda = eigValsSorted(1:numPC); % 提取对应的前numPC个特征值，用于归一化
for k = 1:numPC % 对每个主成分方向进行归一化处理
    alpha(:,k) = alpha(:,k) / sqrt(lambda(k)); % 将特征向量按特征值平方根缩放，满足KPCA正交规范
end
Xkpc = Kc * alpha; % 将样本在核主成分空间中的投影计算出来，得到降维后的特征表示
pcTrain = Xkpc; % 将KPCA主成分时间序列视为训练阶段PINN的输入特征之一
PINN网络结构构建示例
inputDim = size(pcTrain,2) + 1; % 输入维度为KPCA主成分数加时间维度
outputDim = 1; % 输出维度示例设为1，对应单变量的未来状态预测
numHidden = 3; % 隐藏层层数设为3，以平衡表达能力和训练难度
layers = [layers,fullyConnectedLayer(numNeurons,"Name","fc1")]; % 添加首个全连接层，将输入映射到高维隐藏空间
layers = [layers,tanhLayer("Name","tanh1")]; % 添加tanh激活层，提供平滑非线性并方便导数计算
for i = 2:numHidden % 循环添加中间隐藏层
    fcName = sprintf("fc%d",i); % 构造当前全连接层名称字符串
    actName = sprintf("tanh%d",i); % 构造当前tanh激活层名称字符串
    layers = [layers,fullyConnectedLayer(numNeurons,"Name",fcName)]; % 添加当前全连接层，维持固定神经元数量
    layers = [layers,tanhLayer("Name",actName)]; % 添加对应的tanh层，形成交替的线性与非线性结构
end
layers = [layers,fullyConnectedLayer(outputDim,"Name","fcOut")]; % 添加输出层，将隐藏表示映射到目标维度
lgraph = layerGraph(layers); % 将层序列转换为层图，以便构建dlnetwork
netPINN = dlnetwork(lgraph); % 将层图封装为可训练的dlnetwork对象，用于自定义训练循环
PINN损失函数与物理残差示例
idxSample = 1:5:length(pcTrain); % 为计算物理残差子采样训练点，降低计算量
tSample = tTrain(idxSample); % 对应选取子采样时间点，作为方程残差计算位置
targetTrain = xTrain(idxTrain,1); % 取第一个变量的标准化数据作为预测目标
targetSample = targetTrain(idxSample); % 对应物理解约样本点的真实值，用于组合损失
lambdaPhys = 1.0; % 设置物理方程中的比例系数示例，用于构造动力学约束
function [lossTotal,lossData,lossPhys] = pinnLoss(netPINN,tBatch,pcBatch,yBatch,tPhys,pcPhys,muPhys,lambdaPhys) % 定义PINN损失计算函数，包含数据项与物理项
tBatch = dlarray(tBatch','CB'); % 将时间批次数据转为dlarray并设置格式，列为样本批
pcBatch = dlarray(pcBatch','CB'); % 将KPCA特征批次数据转为dlarray格式，与时间对齐
inputBatch = [tBatch;pcBatch]; % 将时间和KPCA特征在行方向拼接，形成网络输入
yPred = forward(netPINN,inputBatch); % 前向传播网络，得到批次预测输出
yPred = squeeze(yPred); % 去除多余维度，使预测输出为简单向量
lossData = mse(yPred,yBatch'); % 计算预测值与真实值之间的均方误差，作为数据拟合损失
tPhys = dlarray(tPhys','CB'); % 将物理约束时间点转为dlarray，用于自动微分
pcPhys = dlarray(pcPhys','CB'); % 将物理约束KPCA特征转为dlarray，与时间拼接
inputPhys = [tPhys;pcPhys]; % 构造物理约束输入，一并输入网络计算
yPhys = squeeze(yPhys); % 压缩形状，方便后续导数计算与运算
dyPhys_dt = dlgradient(sum(yPhys),tPhys); % 利用自动微分对输出关于时间求导，构造时间导数
physResidual = dyPhys_dt - lambdaPhys*(yPhys - muPhys); % 构建物理方程残差示例，形如dy/dt - lambda*(y - mu)
wData = 1.0; % 设置数据损失权重，控制在总损失中的比例
wPhys = 0.1; % 设置物理损失权重，让物理约束在训练中适度发挥作用
end % 函数结束，返回总损失以及两部分子损失
自定义训练循环示例
numEpochs = 200; % 设定最大训练轮数，使网络有充分机会收敛
miniBatchSize = 64; % 每个小批量的样本数，平衡梯度估计稳定性和计算效率
learnRate = 1e-3; % 初始学习率大小，决定每步参数更新幅度
gradDecay = 0.9; % Adam优化器一阶矩估计衰减系数
sqGradDecay = 0.999; % Adam优化器二阶矩估计衰减系数
trailingAvg = []; % 初始化一阶矩动量向量为空，由Adam在训练中逐步累积
trailingAvgSq = []; % 初始化二阶矩动量向量为空，记录梯度平方的指数加权平均
numTrainSample = size(pcTrain,1); % 获取训练样本总数，用于mini-batch切分
numIterPerEpoch = ceil(numTrainSample / miniBatchSize); % 计算每轮迭代次数，用总样本数除小批量大小
figure; % 打开新图形窗口用于实时绘制损失曲线
lossPhysHist = zeros(numEpochs,1); % 预分配物理损失记录数组
    idxShuffle = randperm(numTrainSample); % 每轮对训练样本索引打乱，增加随机性
    xTrainShuf = pcTrain(idxShuffle,:); % 按随机索引重排KPCA特征，打乱样本顺序
    tTrainShuf = tTrain(idxShuffle); % 对应重排时间向量，保持与特征配对
    yTrainShuf = targetTrain(idxShuffle); % 对应重排目标输出，避免固定顺序偏差
    epochLoss = 0; % 初始化本轮平均损失累计变量
    epochPhys = 0; % 初始化本轮物理损失累计变量
    for iter = 1:numIterPerEpoch % 内层循环遍历每个小批量
        idxBatchStart = (iter-1)*miniBatchSize + 1; % 当前小批量起始索引
        idxBatchEnd = min(iter*miniBatchSize,numTrainSample); % 当前小批量结束索引，不超过总样本数
        batchIdx = idxBatchStart:idxBatchEnd; % 当前小批量的索引范围
        tBatch = tTrainShuf(batchIdx); % 小批量的时间数据，用于输入和损失计算
        dlYBatch = dlarray(yBatch','CB'); % 将小批量目标转为dlarray格式，方便与预测对比
        gradFun = dlgradient(@(params) lossWrapper(params,netPINN,tBatch,pcBatch,dlYBatch,tSample,pcSample,muPhys,lambdaPhys),netPINN.Learnables); % 构造梯度函数，通过自动微分计算损失对可学习参数的梯度
        [lossVal,lossDataVal,lossPhysVal] = lossWrapper(netPINN.Learnables,netPINN,tBatch,pcBatch,dlYBatch,tSample,pcSample,muPhys,lambdaPhys); % 同时计算当前小批量的总损失与各组成部分
        [netPINN,trailingAvg,trailingAvgSq] = adamupdate(netPINN,gradFun,trailingAvg,trailingAvgSq,epoch,learnRate,gradDecay,sqGradDecay); % 使用Adam优化器根据梯度和动量更新网络参数
        epochLoss = epochLoss + double(lossVal); % 将当前小批量的总损失加入本轮累积
    end
    lossTrainHist(epoch) = epochLoss / numIterPerEpoch; % 计算本轮平均总损失，记录到历史数组
    lossPhysHist(epoch) = epochPhys / numIterPerEpoch; % 计算本轮平均物理损失，记录用于分析约束效果
    plot(1:epoch,lossTrainHist(1:epoch),'b-'); % 绘制截至当前的训练损失曲线，以蓝色显示
    hold on; % 保持当前图形内容，用于叠加绘制物理损失
    plot(1:epoch,lossPhysHist(1:epoch),'r-'); % 叠加绘制物理损失曲线，以红色展示物理约束变化
    xlabel('Epoch'); % 标注横轴为训练轮数
end
function [lossTotal,lossData,lossPhys] = lossWrapper(params,netTemplate,tBatch,pcBatch,dlYBatch,tPhys,pcPhys,muPhys,lambdaPhys) % 封装损失计算函数，便于在自动微分中调用
netTemplate.Learnables = params; % 将传入的参数结构体赋给网络模板的可学习参数字段

数据读取与预处理示例

clear; % 清空工作区变量，避免历史变量干扰当前实验

clc; % 清空命令行窗口，便于观察本次运行输出

close all; % 关闭所有图形窗口，保证后续绘图环境干净

rawData = readtable(dataFile); % 读取CSV为表格形式，便于按列访问不同变量

timeVec = rawData.Time; % 提取时间列，作为统一时间轴的基础

var1 = rawData.Var1; % 提取第一个变量的时间序列数据

var2 = rawData.Var2; % 提取第二个变量的时间序列数据

var3 = rawData.Var3; % 提取第三个变量的时间序列数据

tStart = timeVec(1); % 记录原始时间序列的起始时间

dt = median(diff(timeVec)); % 估计主采样间隔，用中位数抑制偶发异常差值

var3Uniform = interp1(timeVec,var3,tUniform,'linear','extrap'); % 将变量3插值到统一时间网格，统一时间尺度

nanIdx = isnan(var1Uniform) | isnan(var2Uniform) | isnan(var3Uniform); % 找出任一变量存在NaN的位置，用于统一处理缺失

var1Uniform(nanIdx) = []; % 对应删除变量1在缺失时间点的值，保持长度一致

var2Uniform(nanIdx) = []; % 对应删除变量2在缺失时间点的值，保持数据对齐

mu1 = mean(var1Uniform); % 计算变量1均值，用于中心化处理

sigma1 = std(var1Uniform); % 计算变量1标准差，用于缩放到统一尺度

mu2 = mean(var2Uniform); % 计算变量2均值，为后续标准化提供参数

sigma2 = std(var2Uniform); % 计算变量2标准差，用于拉伸或压缩幅值范围

mu3 = mean(var3Uniform); % 计算变量3均值，以消除直流偏移

sigma3 = std(var3Uniform); % 计算变量3标准差，使各变量量纲相对一致

var1Norm = (var1Uniform - mu1) / sigma1; % 对变量1执行标准化，中心化并缩放

var2Norm = (var2Uniform - mu2) / sigma2; % 对变量2执行标准化，方便与其他变量合并处理

var3Norm = (var3Uniform - mu3) / sigma3; % 对变量3执行标准化，减小量纲差异影响

dataMat = [var1Norm,var2Norm,var3Norm]; % 将三个规范化后的变量按列组合成数据矩阵

nTotal = size(dataMat,1); % 获取总时间步数，用于划分训练验证测试集

nVal = floor(0.15 * nTotal); % 将15%的样本分配给验证集，用于调参与早停判断

idxTrain = 1:nTrain; % 训练集索引范围，从开头到训练样本数

idxVal = (nTrain+1):(nTrain+nVal); % 验证集索引范围，紧随训练集之后

idxTest = (nTrain+nVal+1):nTotal; % 测试集索引范围，为时间序列末尾部分

tTrain = tUniform(idxTrain); % 训练集时间向量，保持与数据同步

tVal = tUniform(idxVal); % 验证集时间向量，用于模型调优过程中指标计算

tTest = tUniform(idxTest); % 测试集时间向量，用于最终预测效果展示

xVal = dataMat(idxVal,:); % 验证集多变量标准化数据矩阵

xTest = dataMat(idxTest,:); % 测试集多变量标准化数据矩阵

CEEMDAN分解示例（单变量）

signal = xTrain(:,1); % 从训练集中选择第一个变量作为示例信号进行CEEMDAN分解

fs = 1 / dt; % 计算采样频率，取决于统一时间步长dt

numEnsemble = 100; % 设置集合次数，用于增加统计稳定性

noiseAmp = 0.2; % 设置噪声幅度比例，控制每次加入噪声的强度

imfSet = []; % 初始化存储IMF的矩阵容器，列对应不同IMF

residual = signal; % 初始化残余信号为原始信号，逐步剥离IMF

for k = 1:maxIMF % 循环提取每一个IMF分量，最多maxIMF个

    imfEnsemble = zeros(N,numEnsemble); % 为当前IMF的集合存储初始化矩阵

        imfEnsemble(:,e) = imfsTmp(:,1); % 将提取到的IMF存入集合矩阵的对应列

    end

    imfMean = mean(imfEnsemble,2); % 对所有集合IMF在列方向取平均，得到该轮的最终IMF

    imfSet = [imfSet,imfMean]; % 将新得到的IMF添加到IMF集合矩阵中

        break; % 终止IMF提取循环，防止无意义的进一步分解

    end

end

trend = residual; % 将最终残余部分视为长期趋势分量

nIMF = size(imfSet,2); % 记录实际得到的IMF数量，可能小于maxIMF

figure; % 新建图形窗口用于展示分解结果

    subplot(nIMF+1,1,i); % 为每个IMF分配一个子图区域，依次排列

    plot(tTrain,imfSet(:,i)); % 绘制第i个IMF随时间变化曲线

    ylabel(['IMF',num2str(i)]); % 标注纵轴为对应IMF序号，便于识别

end

plot(tTrain,trend); % 绘制趋势分量，展示长期变化结构

ylabel('Trend'); % 标注趋势曲线的纵轴含义

xlabel('Time'); % 在最底部子图标注时间轴标签

KPCA特征构建与降维示例

imfFeatures = imfSet; % 将IMF集合视作特征矩阵，每列对应不同时间尺度的分量

featureMat = imfFeatures; % 将单变量IMF特征赋给通用特征矩阵变量，便于后续扩展到多变量

meanFeat = mean(featureMat,1); % 计算各特征维度的均值，用于中心化

Xc = featureMat - meanFeat; % 对特征矩阵做中心化，消除偏移以适应核方法

N = size(Xc,1); % 特征矩阵的样本数，即时间步数

sigmaKernel = 1.0; % 设置高斯核的核宽度参数，用于控制相似性衰减

for i = 1:N % 外层循环遍历第一个样本索引

    for j = 1:N % 内层循环遍历第二个样本索引

        diffVec = Xc(i,:) - Xc(j,:); % 计算样本i和样本j在特征空间的差向量

        dist2 = sum(diffVec.^2); % 计算差向量的平方欧氏距离

        K(i,j) = exp(-dist2/(2*sigmaKernel^2)); % 根据高斯核公式计算核矩阵元素

    end

end

Kc = K - oneN*K - K*oneN + oneN*K*oneN; % 对核矩阵进行中心化，使其对应中心化特征空间

 [EigVec,EigValMat] = eig(Kc); % 对中心化核矩阵进行特征分解，得到特征向量和特征值矩阵

eigVals = diag(EigValMat); % 提取特征值对角线元素，转为向量形式

[eigValsSorted,idxSort] = sort(eigVals,'descend'); % 将特征值按降序排序，记录排序索引

EigVecSorted = EigVec(:,idxSort); % 将特征向量按对应特征值大小重排

varExplained = eigValsSorted / sum(eigValsSorted); % 计算每个核主成分的方差贡献率

cumVar = cumsum(varExplained); % 计算累积方差贡献率，作为选取成分数量依据

numPC = find(cumVar >= 0.95,1); % 找出使累积贡献率达到95%的最小主成分个数

alpha = EigVecSorted(:,1:numPC); % 提取前numPC个特征向量，用作KPCA投影系数

lambda = eigValsSorted(1:numPC); % 提取对应的前numPC个特征值，用于归一化

for k = 1:numPC % 对每个主成分方向进行归一化处理

    alpha(:,k) = alpha(:,k) / sqrt(lambda(k)); % 将特征向量按特征值平方根缩放，满足KPCA正交规范

end

Xkpc = Kc * alpha; % 将样本在核主成分空间中的投影计算出来，得到降维后的特征表示

pcTrain = Xkpc; % 将KPCA主成分时间序列视为训练阶段PINN的输入特征之一

PINN网络结构构建示例

inputDim = size(pcTrain,2) + 1; % 输入维度为KPCA主成分数加时间维度

outputDim = 1; % 输出维度示例设为1，对应单变量的未来状态预测

numHidden = 3; % 隐藏层层数设为3，以平衡表达能力和训练难度

layers = [layers,fullyConnectedLayer(numNeurons,"Name","fc1")]; % 添加首个全连接层，将输入映射到高维隐藏空间

layers = [layers,tanhLayer("Name","tanh1")]; % 添加tanh激活层，提供平滑非线性并方便导数计算

for i = 2:numHidden % 循环添加中间隐藏层

    fcName = sprintf("fc%d",i); % 构造当前全连接层名称字符串

    actName = sprintf("tanh%d",i); % 构造当前tanh激活层名称字符串

    layers = [layers,fullyConnectedLayer(numNeurons,"Name",fcName)]; % 添加当前全连接层，维持固定神经元数量

    layers = [layers,tanhLayer("Name",actName)]; % 添加对应的tanh层，形成交替的线性与非线性结构

end

layers = [layers,fullyConnectedLayer(outputDim,"Name","fcOut")]; % 添加输出层，将隐藏表示映射到目标维度

lgraph = layerGraph(layers); % 将层序列转换为层图，以便构建dlnetwork

netPINN = dlnetwork(lgraph); % 将层图封装为可训练的dlnetwork对象，用于自定义训练循环

PINN损失函数与物理残差示例

idxSample = 1:5:length(pcTrain); % 为计算物理残差子采样训练点，降低计算量

tSample = tTrain(idxSample); % 对应选取子采样时间点，作为方程残差计算位置

targetTrain = xTrain(idxTrain,1); % 取第一个变量的标准化数据作为预测目标

targetSample = targetTrain(idxSample); % 对应物理解约样本点的真实值，用于组合损失

lambdaPhys = 1.0; % 设置物理方程中的比例系数示例，用于构造动力学约束

function [lossTotal,lossData,lossPhys] = pinnLoss(netPINN,tBatch,pcBatch,yBatch,tPhys,pcPhys,muPhys,lambdaPhys) % 定义PINN损失计算函数，包含数据项与物理项

tBatch = dlarray(tBatch','CB'); % 将时间批次数据转为dlarray并设置格式，列为样本批

pcBatch = dlarray(pcBatch','CB'); % 将KPCA特征批次数据转为dlarray格式，与时间对齐

inputBatch = [tBatch;pcBatch]; % 将时间和KPCA特征在行方向拼接，形成网络输入

yPred = forward(netPINN,inputBatch); % 前向传播网络，得到批次预测输出

yPred = squeeze(yPred); % 去除多余维度，使预测输出为简单向量

lossData = mse(yPred,yBatch'); % 计算预测值与真实值之间的均方误差，作为数据拟合损失

tPhys = dlarray(tPhys','CB'); % 将物理约束时间点转为dlarray，用于自动微分

pcPhys = dlarray(pcPhys','CB'); % 将物理约束KPCA特征转为dlarray，与时间拼接

inputPhys = [tPhys;pcPhys]; % 构造物理约束输入，一并输入网络计算

yPhys = squeeze(yPhys); % 压缩形状，方便后续导数计算与运算

dyPhys_dt = dlgradient(sum(yPhys),tPhys); % 利用自动微分对输出关于时间求导，构造时间导数

physResidual = dyPhys_dt - lambdaPhys*(yPhys - muPhys); % 构建物理方程残差示例，形如dy/dt - lambda*(y - mu)

wData = 1.0; % 设置数据损失权重，控制在总损失中的比例

wPhys = 0.1; % 设置物理损失权重，让物理约束在训练中适度发挥作用

end % 函数结束，返回总损失以及两部分子损失

自定义训练循环示例

numEpochs = 200; % 设定最大训练轮数，使网络有充分机会收敛

miniBatchSize = 64; % 每个小批量的样本数，平衡梯度估计稳定性和计算效率

learnRate = 1e-3; % 初始学习率大小，决定每步参数更新幅度

gradDecay = 0.9; % Adam优化器一阶矩估计衰减系数

sqGradDecay = 0.999; % Adam优化器二阶矩估计衰减系数

trailingAvg = []; % 初始化一阶矩动量向量为空，由Adam在训练中逐步累积

trailingAvgSq = []; % 初始化二阶矩动量向量为空，记录梯度平方的指数加权平均

numTrainSample = size(pcTrain,1); % 获取训练样本总数，用于mini-batch切分

numIterPerEpoch = ceil(numTrainSample / miniBatchSize); % 计算每轮迭代次数，用总样本数除小批量大小

figure; % 打开新图形窗口用于实时绘制损失曲线

lossPhysHist = zeros(numEpochs,1); % 预分配物理损失记录数组

    idxShuffle = randperm(numTrainSample); % 每轮对训练样本索引打乱，增加随机性

    xTrainShuf = pcTrain(idxShuffle,:); % 按随机索引重排KPCA特征，打乱样本顺序

    tTrainShuf = tTrain(idxShuffle); % 对应重排时间向量，保持与特征配对

    yTrainShuf = targetTrain(idxShuffle); % 对应重排目标输出，避免固定顺序偏差

    epochLoss = 0; % 初始化本轮平均损失累计变量

    epochPhys = 0; % 初始化本轮物理损失累计变量

    for iter = 1:numIterPerEpoch % 内层循环遍历每个小批量

        idxBatchStart = (iter-1)*miniBatchSize + 1; % 当前小批量起始索引

        idxBatchEnd = min(iter*miniBatchSize,numTrainSample); % 当前小批量结束索引，不超过总样本数

        batchIdx = idxBatchStart:idxBatchEnd; % 当前小批量的索引范围

        tBatch = tTrainShuf(batchIdx); % 小批量的时间数据，用于输入和损失计算

        dlYBatch = dlarray(yBatch','CB'); % 将小批量目标转为dlarray格式，方便与预测对比

        gradFun = dlgradient(@(params) lossWrapper(params,netPINN,tBatch,pcBatch,dlYBatch,tSample,pcSample,muPhys,lambdaPhys),netPINN.Learnables); % 构造梯度函数，通过自动微分计算损失对可学习参数的梯度

        [lossVal,lossDataVal,lossPhysVal] = lossWrapper(netPINN.Learnables,netPINN,tBatch,pcBatch,dlYBatch,tSample,pcSample,muPhys,lambdaPhys); % 同时计算当前小批量的总损失与各组成部分

        [netPINN,trailingAvg,trailingAvgSq] = adamupdate(netPINN,gradFun,trailingAvg,trailingAvgSq,epoch,learnRate,gradDecay,sqGradDecay); % 使用Adam优化器根据梯度和动量更新网络参数

        epochLoss = epochLoss + double(lossVal); % 将当前小批量的总损失加入本轮累积

    end

    lossTrainHist(epoch) = epochLoss / numIterPerEpoch; % 计算本轮平均总损失，记录到历史数组

    lossPhysHist(epoch) = epochPhys / numIterPerEpoch; % 计算本轮平均物理损失，记录用于分析约束效果

    plot(1:epoch,lossTrainHist(1:epoch),'b-'); % 绘制截至当前的训练损失曲线，以蓝色显示

    hold on; % 保持当前图形内容，用于叠加绘制物理损失

    plot(1:epoch,lossPhysHist(1:epoch),'r-'); % 叠加绘制物理损失曲线，以红色展示物理约束变化

    xlabel('Epoch'); % 标注横轴为训练轮数

end

function [lossTotal,lossData,lossPhys] = lossWrapper(params,netTemplate,tBatch,pcBatch,dlYBatch,tPhys,pcPhys,muPhys,lambdaPhys) % 封装损失计算函数，便于在自动微分中调用

netTemplate.Learnables = params; % 将传入的参数结构体赋给网络模板的可学习参数字段

更多详细内容请访问

http://【多变量时序预测】MATLAB实现基于CEEMDAN-KPCA-PINN完全集合经验模态分解与自适应噪声（CEEMDAN）结合核主成分分析（KPCA）和物理信息神经网络（PINN）进行多变量时序预测_带注意力机制的BiTCN模型代码资源-CSDN下载  https://download.csdn.net/download/xiaoxingkongyuxi/90440213

 https://download.csdn.net/download/xiaoxingkongyuxi/90440213

http:// https://download.csdn.net/download/xiaoxingkongyuxi/90440213