基于变分模态分解和麻雀算法优化双向长短期记忆网络的多维时间序列预测，VMD-SSA-BILSTM多维时间序列预测 MATLAB代

clc; clear; close all; format compact ; warning off %% 读取数据 result = xlsread('data.xlsx','sheet1','A2：H1001'); f=result(:,end); % 采样时间 t = 1:length(f); %% 数据分析 num_samples = length(result); % 样本个数 kim = 3; % 延时步长（前面多行历史数据作为自变量） zim = 1; % 跨zim个时间点进行预测 nim = size(result, 2); %% 划分数据集 for i = 1: num_samples - kim - zim + 1 res(i, :) = [reshape(result(i: i + kim - 1 , :)', 1, ... kim * nim), result(i + kim , end)]; end %% 数据集分析 outdim = 1; % 最后一列为输出 num_size = 0.8; % 训练集占数据集比例 num_train_s = round(num_size * num_samples); % 训练集样本个数 f_ = size(res, 2) - outdim; % 输入特征长度 %% 训练集和测试集,作为绘图用 PP_train = res(1: num_train_s, 1: f_)'; TT_train = res(1: num_train_s, f_ + 1: end)'; PP_test = res(num_train_s + 1: end, 1: f_)'; TT_test = res(num_train_s + 1: end, f_ + 1: end)'; f_=size(PP_train, 1); % 输入特征维度 outdim = 1; % 最后一列为输出 %% 训练集和测试集 M = size(PP_train, 2); N = size(PP_test, 2); %% VMD分解 s = result(:,end); alpha = 2000; % moderate bandwidth constraint：适度的带宽约束/惩罚因子 tau = 0; % noise-tolerance (no strict fidelity enforcement)：噪声容限（没有严格的保真度执行） K = 7; % modes：分解的模态数 DC = 0; % no DC part imposed：无直流部分 init = 1; % initialize omegas uniformly ：omegas的均匀初始化? tol = 1e-7; %--------------- Run actual VMD code:数据进行vmd分解--------------------------- [u, u_hat, omega] = VMD(s, alpha, tau, K, DC, init, tol); [m,~]=size(u); %程序中omega就是中心频率，k为多少，分解的个数就为几个， %omega就有几列，然后就每列的均值，程序里面有。 res1 = s'-sum(u,1);%res %% 分解重组VS原始图 figure for i=1:m subplot(m,1,i) plot(u(i,:),'b-','linewidth',1) ylabel(['IMF',num2str(i)]); end %% 分解完成 imf=u; %% VMD-LSTM预测 pre_result2=[]; pre_result1=[]; %% 对每个分量建模 for i=1:size(imf,1) disp(['对第',num2str(i),'个分量建模']) %% 数据分析 samples = length(imf); % 样本个数 kim = 3; % 延时步长（前面多行历史数据作为自变量） result1=[result(:,1:end-1)';imf(i,:)]'; zim = 1; % 跨zim个时间点进行预测 nim = size(result1, 2); %% 划分数据集 for i = 1: samples - kim - zim + 1 res(i, :) = [reshape(result1(i: i + kim - 1 , :)', 1, ... kim * nim), result1(i + kim , end)]; end %归一化 %% 数据集分析 size1 = 0.8; % 训练集占数据集比例 train_s = round(size1 * samples); % 训练集样本个数 f_ = size(res, 2) - outdim; % 输入特征长度 %% 训练集和测试集,作为绘图用 P1_train = res(1: num_train_s, 1: f_)'; T1_train = res(1: num_train_s, f_ + 1: end)'; P1_test = res(num_train_s + 1: end, 1: f_)'; T1_test = res(num_train_s + 1: end, f_ + 1: end)'; %% 划分训练集和测试集 m = size(P1_train, 2); n = size(P1_test, 2); %% 数据归一化 [p_train, ps_input] = mapminmax(P1_train, 0, 1); p_test = mapminmax('apply', P1_test, ps_input); [t_train, ps_output] = mapminmax(T1_train, 0, 1); t_test = mapminmax('apply', T1_test, ps_output); %% 格式转换 for i = 1 : M vp_train{i, 1} = p_train(:, i); vt_train{i, 1} = t_train(:, i); end for i = 1 : N vp_test{i, 1} = p_test(:, i); vt_test{i, 1} = t_test(:, i); end %% 创建待优化函数 ObjFcn = @CostFunction; %% 优化参数设置 SearchAgents = 6; % 种群数量 Max_iterations = 6 ; % 迭代次数 lowerbound = [1e-10 0.0001 10 ];%三个参数的下限 upperbound = [1e-2 0.002 30 ];%三个参数的上限 dimension = 3;%数量，即要优化的BiLSTM参数个数 %% 优化BiLSTM [Best_score,Best_pos,Convergence_curve]=SSA(SearchAgents,Max_iterations,lowerbound,upperbound,dimension,ObjFcn) %% 得到最优参数 Best_pos(1,3) =round( Best_pos(1,3)); % 最佳隐藏层节点数 NumOfUnits = Best_pos(1,3); InitialLearnRate = Best_pos(1,2) % 最佳初始学习率 L2Regularization = Best_pos(1,1); % 最佳L2正则化系数 %% 创建网络， layers = [ ... sequenceInputLayer(size(P1_train,1)) % 输入层,即输入的特征变量个数 bilstmLayer(NumOfUnits) % bilstm层 reluLayer % Relu激活层 fullyConnectedLayer(outdim) % 全连接层 regressionLayer]; % 回归层 % 参数设置 options = trainingOptions('adam', ... % 优化算法Adam 'MaxEpochs', 100, ... % 最大训练次数 'GradientThreshold', 1, ... % 梯度阈值 'InitialLearnRate', InitialLearnRate, ... % 初始学习率 'LearnRateSchedule', 'piecewise', ... % 学习率调整 'LearnRateDropPeriod',80, ... % 训练次后开始调整学习率 'LearnRateDropFactor',0.2, ... % 学习率调整因子 'L2Regularization', L2Regularization, ... % 正则化参数 'ExecutionEnvironment', 'cpu',... % 训练环境 'Verbose', 0); % 关闭优化过程 % 'Plots', 'training-progress' % 画出曲线 %% 训练 net = trainNetwork(vp_train, vt_train, layers, options); %% 预测 t_sim1 = predict(net, vp_train); t_sim2 = predict(net, vp_test); %% 数据反归一化 T_sim1 = mapminmax('reverse', t_sim1, ps_output); T_sim2 = mapminmax('reverse', t_sim2, ps_output); T_sim1 = cell2mat(T_sim1); T_sim2 = cell2mat(T_sim2); T_sim1 = double(T_sim1'); T_sim2 = double(T_sim2'); pre_result1=[pre_result1;T_sim1]; pre_result2=[pre_result2;T_sim2]; end %% 各分量预测的结果相加 pre_value1=sum(pre_result1); pre_value2=sum(pre_result2); %% 变量换算 T_sim1=pre_value1; T_sim2=pre_value2; T_train=TT_train; T_test=TT_test; %% 均方根误差 error1 = sqrt(sum((T_sim1 - T_train).^2) ./ M); error2 = sqrt(sum((T_sim2 - T_test ).^2) ./ N); %% 绘图 figure plot(1: M, TT_train, '-b','linewidth',1) hold on plot(1: M, T_sim1,'-r','linewidth',1) legend('真实值','VMD-SSA-BILSTM预测值') xlabel('预测样本') ylabel('预测结果') string = {'训练集预测结果对比'; ['RMSE=' num2str(error1)]}; title(string) xlim([1, M]) grid figure plot(1: N, TT_test, '-r','linewidth',1) hold on plot(1: N, T_sim2,'-b','linewidth',1) legend('真实值','VMD-SSA-BILSTM预测值') xlabel('预测样本') ylabel('预测结果') string = {'测试集预测结果对比';['RMSE=' num2str(error2)]}; title(string) xlim([1, N]) grid %% 相关指标计算 % R2 R1 = 1 - norm(T_train - T_sim1)^2 / norm(T_train - mean(T_train))^2; R2 = 1 - norm(T_test - T_sim2)^2 / norm(T_test - mean(T_test ))^2; disp(['训练集数据的R2为：', num2str(R1)]) disp(['测试集数据的R2为：', num2str(R2)]) % MAE mae1 = sum(abs(T_sim1 - T_train)) ./ M ; mae2 = sum(abs(T_sim2 - T_test)) ./ N ; disp(['训练集数据的MAE为：', num2str(mae1)]) disp(['测试集数据的MAE为：', num2str(mae2)]) % MBE mbe1 = sum(T_sim1 - T_train) ./ M ; mbe2 = sum(T_sim2 - T_test) ./ N ; disp(['训练集数据的MBE为：', num2str(mbe1)]) disp(['测试集数据的MBE为：', num2str(mbe2)]) %% 记录结果 save VMD-SSA-BILSTM结果/VMD-SSA-BILSTM T_test T_sim2