GOOSE-BP多变量风电功率时序预测，含优化前后对比（Matlab完整源码和数据）

clc; clear close all %% pwd X = xlsread('风电场预测.xlsx'); n_in = 5; % 输入前5个时刻的所有数据，包括：测风塔10m风速(m/s) 测风塔30m风速(m/s) 测风塔50m风速(m/s) 测风塔70m风速(m/s) 轮毂高度风速(m/s) 测风塔10m风向(°) 测风塔30m风向(°) 测风塔50m风向(°) 测风塔70m风向(°) 轮毂高度风向(°) 温度(°) 气压(hPa) 湿度(%) 实际发电功率（mw） n_out = 1 ; % 此程序为单步预测，因此请将n_out设置为1 or_dim = size(X,2) ; % 记录特征数据维度 num_samples = 2000; % 制作2000个样本。 scroll_window = 2; %如果等于1，下一个样本从第二行开始取。如果等于2，下一个样本从第三行开始取 [res] = data_collation(X, n_in, n_out, or_dim, scroll_window, num_samples); % 训练集和测试集划分%% num_size = 0.8; % 训练集占数据集比例 %% num_train_s = round(num_size * num_samples); % 训练集样本个数 %% %% 以下几行代码是为了方便归一化，一般不需要更改！ P_train = res(1: num_train_s,1); input_train = reshape(cell2mat(P_train)',n_in*or_dim,num_train_s); T_train = res(1: num_train_s,2); output_train = cell2mat(T_train)'; P_test = res(num_train_s+1: end,1); input_test = reshape(cell2mat(P_test)',n_in*or_dim,num_samples-num_train_s); T_test = res(num_train_s+1: end,2); output_test = cell2mat(T_test)'; %% 数据归一化 [inputn,inputps]=mapminmax(input_train,0,1); [outputn,outputps]=mapminmax(output_train); inputn_test=mapminmax('apply',input_test,inputps); %% 获取输入层节点、输出层节点个数 inputnum=size(input_train,1); outputnum=size(output_train,1); disp('/////////////////////////////////') disp('神经网络结构...') disp(['输入层的节点数为：',num2str(inputnum)]) disp(['输出层的节点数为：',num2str(outputnum)]) disp(' ') disp('隐含层节点的确定过程...') %确定隐含层节点个数 %采用经验公式hiddennum=sqrt(m+n)+a，m为输入层节点个数，n为输出层节点个数，a一般取为1-10之间的整数 MSE=1e+5; %初始化最小误差 for hiddennum=fix(sqrt(inputnum+outputnum))+1:fix(sqrt(inputnum+outputnum))+5 %构建网络 net=newff(inputn,outputn,hiddennum,{'tansig','purelin'},'trainlm');% 建立模型 % 网络参数 net.trainParam.epochs=1000; % 训练次数 net.trainParam.lr=0.01; % 学习速率 net.trainParam.goal=0.000001; % 训练目标最小误差 net.trainParam.showWindow=0; %隐藏仿真界面 % 网络训练 net=train(net,inputn,outputn); an0=sim(net,inputn); %仿真结果,依旧采用训练集进行测试 test_simu0=mapminmax('reverse',an0,outputps); %把仿真得到的数据还原为原始的数量级 mse0=mse(test_simu0,output_train); %仿真的均方误差 disp(['隐含层节点数为',num2str(hiddennum),'时，训练集的均方误差为：',num2str(mse0)]) %更新最佳的隐含层节点 if mse0<MSE MSE=mse0; hiddennum_best=hiddennum; end end clear net0 disp(['最佳的隐含层节点数为：',num2str(hiddennum_best),'，训练集的均方误差为：',num2str(MSE)]) %% 构建最佳隐含层节点的BP神经网络 disp(' ') disp('标准的BP神经网络：') net0=newff(inputn,outputn,hiddennum_best,{'tansig','purelin'},'trainlm');% 建立模型 %网络参数配置 net0.trainParam.epochs=1000; % 训练次数，这里设置为1000次 net0.trainParam.lr=0.01; % 学习速率，这里设置为0.01 net0.trainParam.goal=0.00001; % 训练目标最小误差，这里设置为0.0001 net0.trainParam.show=25; % 显示频率，这里设置为每训练25次显示一次 net0.trainParam.mc=0.01; % 动量因子 net0.trainParam.min_grad=1e-6; % 最小性能梯度 net0.trainParam.max_fail=6; % 最高失败次数 net0.trainParam.showWindow = false; net0.trainParam.showCommandLine = false; %隐藏仿真界面 %开始训练 net0=train(net0,inputn,outputn); %预测 an0=sim(net0,inputn_test); %用训练好的模型进行仿真 %预测结果反归一化与误差计算 test_simu0=mapminmax('reverse',an0,outputps); %把仿真得到的数据还原为原始的数量级 %误差指标 mse0=mse(output_test,test_simu0); %% 标准BP神经网络作图 figure('Position',[100 300 800 200]) plot(output_test,'-*') hold on plot(test_simu0,'-o') legend('真实y','预测的y') xlabel('样本数') ylabel('风电功率') title(['标准BP神经网络预测结果',newline,'MSE误差为：',num2str(mse0)]) disp(['标准神经网络测试集的均方误差为：',num2str(mse0)]) dim=inputnum*hiddennum_best+hiddennum_best+hiddennum_best*outputnum+outputnum; %自变量个数 SearchAgents_no = 5; Max_iter = 10; lb = -2.*ones(1,dim); ub = 2.*ones(1,dim); net=newff(inputn,outputn,hiddennum_best,{'tansig','purelin'},'trainlm');% 建立模型 %网络参数配置 net.trainParam.epochs=1000; % 训练次数，这里设置为1000次 net.trainParam.lr=0.01; % 学习速率，这里设置为0.01 net.trainParam.goal=0.00001; % 训练目标最小误差，这里设置为0.0001 net.trainParam.show=25; % 显示频率，这里设置为每训练25次显示一次 net.trainParam.mc=0.01; % 动量因子 net.trainParam.min_grad=1e-6; % 最小性能梯度 net.trainParam.max_fail=6; % 最高失败次数 net.trainParam.showWindow = false; net.trainParam.showCommandLine = false; %隐藏仿真界面 fobj = @(x)fitness(x,inputnum,hiddennum_best,outputnum,net,inputn,outputn,inputn_test,outputps,output_test); temp=clock;temp=sum(temp(4:6))*sum(temp(2:3));temp=round(temp/10);setdemorandstream(temp);%此行代码用于生成随机数种子，确保结果可以复现 %% 调用算法 disp('正在优化，请等待……') [fMin , bestX, Convergence_curve ] =GOOSE(SearchAgents_no,Max_iter,lb,ub,dim,fobj); %% 绘制进化曲线 figure plot(Convergence_curve,'r-','linewidth',2) xlabel('进化代数') ylabel('均方误差') legend('最佳适应度') title('进化曲线') setdemorandstream(temp);%此行代码用于生成随机数种子，确保结果可以复现 [~,optimize_test_simu]=fitness(bestX,inputnum,hiddennum_best,outputnum,net,inputn,outputn,inputn_test,outputps,output_test); %% 比较算法预测值 str={'真实值','标准BP','优化后BP'}; figure('Units', 'pixels', ... 'Position', [300 300 860 370]); plot(output_test,'-','Color',[0.8500 0.3250 0.0980]) hold on plot(test_simu0,'-.','Color',[0.4940 0.1840 0.5560]) hold on plot(optimize_test_simu,'-','Color',[0.4660 0.6740 0.1880]) legend(str) set (gca,"FontSize",12,'LineWidth',1.2) box off legend Box off %% 比较算法误差 test_y = output_test; Test_all = []; y_test_predict = test_simu0; [test_MAE,test_MAPE,test_MSE,test_RMSE,test_R2]=calc_error(y_test_predict,test_y); Test_all=[Test_all;test_MAE test_MAPE test_MSE test_RMSE test_R2]; y_test_predict = optimize_test_simu; [test_MAE,test_MAPE,test_MSE,test_RMSE,test_R2]=calc_error(y_test_predict,test_y); Test_all=[Test_all;test_MAE test_MAPE test_MSE test_RMSE test_R2]; str={'真实值','标准BP','优化后BP'}; str1=str(2:end); str2={'MAE','MAPE','MSE','RMSE','R2'}; data_out=array2table(Test_all); data_out.Properties.VariableNames=str2; data_out.Properties.RowNames=str1; disp(data_out) %% 柱状图 MAE MAPE RMSE 柱状图适合量纲差别不大的 color= [0.66669 0.1206 0.108 0.1339 0.7882 0.8588 0.1525 0.6645 0.1290 0.8549 0.9373 0.8275 0.1551 0.2176 0.8627 0.7843 0.1412 0.1373 0.2000 0.9213 0.8176 0.5569 0.8118 0.7882 1.0000 0.5333 0.5176]; figure('Units', 'pixels', ... 'Position', [300 300 660 375]); plot_data_t=Test_all(:,[1,2,4])'; b=bar(plot_data_t,0.8); hold on for i = 1 : size(plot_data_t,2) x_data(:, i) = b(i).XEndPoints'; end for i =1:size(plot_data_t,2) b(i).FaceColor = color(i,:); b(i).EdgeColor=[0.3353 0.3314 0.6431]; b(i).LineWidth=1.2; end for i = 1 : size(plot_data_t,1)-1 xilnk=(x_data(i, end)+ x_data(i+1, 1))/2; b1=xline(xilnk,'--','LineWidth',1.2); hold on end ax=gca; legend(b,str1,'Location','best') ax.XTickLabels ={'MAE', 'MAPE', 'RMSE'}; set(gca,"FontSize",10,"LineWidth",1) box off legend box off %% 二维图 figure plot_data_t1=Test_all(:,[1,5])'; MarkerType={'s','o','pentagram','^','v'}; for i = 1 : size(plot_data_t1,2) scatter(plot_data_t1(1,i),plot_data_t1(2,