VMD滚动分解+LSTM多变量时序预测，防信息泄露，MATLAB代码-平芜编程栈

研究背景

该研究聚焦于多变量时间序列预测问题，旨在处理目标序列的非平稳、非线性特征。传统方法直接对原始序列建模，难以捕捉多尺度波动规律。代码采用变分模态分解（VMD）将目标序列分解为若干相对平稳的本征模态函数（IMF），结合长短期记忆网络（LSTM）对每个模态分别建模预测，并通过滚动分解策略严格避免未来信息泄露，实现符合真实预测场景的序列外推。研究背景对应能源负荷预测、金融指数预测、气象要素预报等实际应用中，需要利用历史多变量信息预测未来关键指标的需求。

主要功能

VMD 自适应分解：对滚动增长的目标序列进行在线模态分解，分离不同中心频率的波动成分。
多变量时间窗口构造：将历史窗口内的所有特征（原始输入特征及目标历史值）展平作为 LSTM 输入，预测未来指定步长的模态值。
LSTM 多模态并行预测：为每个 IMF 独立训练 LSTM 网络，学习其局部时序规律。
滚动增量预测与防泄露机制：每预测一个测试点后，将其真实值纳入历史，重新分解并更新模型，完全模拟在线预测流程。
多维评估与可视化：提供 RMSE、MAE、MAPE、R² 指标，绘制预测对比曲线、回归散点图、误差直方图、分解模态图、频谱图及模态贡献条形图。

算法步骤与技术路线

技术路线（整体流程）

数据加载 → 划分训练/测试集 → 初始化已知历史序列 ↓ for 每个测试点: 1. 对当前已知目标序列执行 VMD，得到 K 个 IMF 2. for 每个 IMF: - 用历史窗口数据构造 LSTM 监督样本 - 归一化、训练 LSTM 网络 - 利用最后 kim 步数据预测未来第 zim 步的模态值 3. 求和集成各模态预测值得到目标预测 4. 获取当前测试点真实值，扩展已知序列（滚动） ↓ 输出性能指标与可视化图表

详细算法步骤

数据预处理
- 读取 Excel 数据（前 4 列为特征 A~D，第 5 列为待预测目标 E）。
- 按 7:3 比例划分训练集和测试集。
滚动分解与预测循环（共N_test次迭代）
- VMD 分解：调用VMD()函数，对当前累积的真实目标序列进行分解，参数K=6，惩罚因子α=2000，无噪声容忍度，无直流分量，收敛容差1e-7。
- LSTM 样本构造：使用辅助函数create_lstm_samples_vmd，以kim=4的历史窗口展平所有特征列，输出标签为未来zim=1步的模态值。
- 归一化与训练：输入输出分别映射到 [0,1]，构建单层 LSTM 网络（64 个隐层单元，输出模式last），采用 Adam 优化器训练 200 轮。
- 单步预测与集成：对每个模态预测下一步值并反归一化，求和得到目标预测值。
- 真实值吸收：将测试点真实值（包括全部特征列）追加到已知数据矩阵，供下一轮分解使用。
评估与可视化
- 计算四项回归指标。
- 绘制 6 幅分析图表，涵盖预测效果、误差分布、模态分解及频谱特性。

核心公式原理

1. VMD 变分约束模型

VMD 通过求解以下变分问题将信号 ( f ) 分解为 ( K ) 个模态 ( u_k )：

min⁡{uk},{ωk}{∑k=1K∥∂t[(δ(t)+jπt)∗uk(t)]e−jωkt∥22} \min_{\{u_k\},\{\omega_k\}} \left\{ \sum_{k=1}^{K} \left\| \partial_t \left[ \left( \delta(t) + \frac{j}{\pi t} \right) * u_k(t) \right] e^{-j\omega_k t} \right\|_2^2 \right\}{uk},{ωk}min{k=1∑K∂t[(δ(t)+πtj)∗uk(t)]e−jωkt22}

s.t.∑k=1Kuk(t)=f(t) \text{s.t.} \quad \sum_{k=1}^{K} u_k(t) = f(t)s.t.k=1∑Kuk(t)=f(t)

利用二次惩罚项和拉格朗日乘子转化为无约束问题，通过交替方向乘子法（ADMM）在频域迭代更新：

模态更新（维纳滤波）：
u^kn+1(ω)=f^(ω)−∑i≠ku^i(ω)+λ^(ω)21+2α(ω−ωk)2 \hat{u}_k^{n+1}(\omega) = \frac{\hat{f}(\omega) - \sum_{i\neq k}\hat{u}_i(\omega) + \frac{\hat{\lambda}(\omega)}{2}}{1 + 2\alpha(\omega - \omega_k)^2}u^kn+1(ω)=1+2α(ω−ωk)2f^(ω)−∑i=ku^i(ω)+2λ^(ω)
中心频率更新：
ωkn+1=∫0∞ω∣u^k(ω)∣2dω∫0∞∣u^k(ω)∣2dω \omega_k^{n+1} = \frac{\int_0^\infty \omega |\hat{u}_k(\omega)|^2 d\omega}{\int_0^\infty |\hat{u}_k(\omega)|^2 d\omega}ωkn+1=∫0∞∣u^k(ω)∣2dω∫0∞ω∣u^k(ω)∣2dω

2. LSTM 单元与回归输出

LSTM 层通过门控机制控制信息流动：

ft=σ(Wf⋅[ht−1,xt]+bf)it=σ(Wi⋅[ht−1,xt]+bi)C~t=tanh⁡(WC⋅[ht−1,xt]+bC)Ct=ft⊙Ct−1+it⊙C~tot=σ(Wo⋅[ht−1,xt]+bo)ht=ot⊙tanh⁡(Ct) \begin{aligned} f_t &= \sigma(W_f \cdot [h_{t-1}, x_t] + b_f) \\ i_t &= \sigma(W_i \cdot [h_{t-1}, x_t] + b_i) \\ \tilde{C}_t &= \tanh(W_C \cdot [h_{t-1}, x_t] + b_C) \\ C_t &= f_t \odot C_{t-1} + i_t \odot \tilde{C}_t \\ o_t &= \sigma(W_o \cdot [h_{t-1}, x_t] + b_o) \\ h_t &= o_t \odot \tanh(C_t) \end{aligned}ftitC~tCtotht=σ(Wf⋅[ht−1,xt]+bf)=σ(Wi⋅[ht−1,xt]+bi)=tanh(WC⋅[ht−1,xt]+bC)=ft⊙Ct−1+it⊙C~t=σ(Wo⋅[ht−1,xt]+bo)=ot⊙tanh(Ct)

本代码中 LSTM 输出模式为'last'，即取最后一个时间步的隐状态接入全连接层，输出标量预测值，采用均方误差（MSE）作为损失函数。

参数设定

参数名称	取值	说明
`kim`	4	历史时间窗口长度（过去步数）
`zim`	1	预测步长（单步预测）
`K_imf`	6	VMD 分解模态个数
`alpha`	2000	VMD 惩罚因子（平衡模态带宽与重构精度）
`tau`	0	拉格朗日乘子更新步长（无噪声）
`DC`	0	不保留直流分量
`init`	1	中心频率均匀初始化
`tol`	1e-7	VMD 收敛容差
`train_ratio`	0.7	训练集划分比例
`lstm_units`	64	LSTM 隐藏层神经元数量
`max_epochs`	200	最大训练轮数
`mini_batch_size`	24	批量大小
`initial_lr`	0.001	初始学习率，每 80 轮下降为原来的 0.2