%% 电力负荷预测
%% 电力负荷预测 滚动24程序
clc
clear
%% Electric load forecasting
load x.mat
load y.mat
n=1440;
x1=x
M=50;
method='unbiased';
[E,V,A,R,C]=ssa(x1,M,method);%去噪
dl=zeros(length(x),1);
%% reconstruction.
for i=1:17% rational changing is great.重构
dl=R(:,i)+dl;
end
for i=1:12
sh(i,:)=dl(i:n-12+i)';
end
input = sh(:,1:end-1)';
output= dl(13:end);
%% 输入输出数据归一化
[inputn,inputps]=mapminmax(input');
[outputn,outputps]=mapminmax(output');
inputn=inputn';
outputn=outputn';
%% 输入输出数据归一化
[inputn,inputps]=mapminmax(input');
[outputn,outputps]=mapminmax(output');
inputn=inputn';
outputn=outputn';
%%
inputnum=size(input,2); %输入节点个数4
outputnum=size(output,2); %输出节点个数1
hiddennum=25;
%%
tic
n2=20;%鸟窝个数
[bestnest,fmin]=cuckoo_search(n2,inputnum,hiddennum,outputnum,inputn,outputn);
toc
%% 提取
w1=bestnest(1:inputnum*hiddennum);
B1=bestnest(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);
w2=bestnest(inputnum*hiddennum+hiddennum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum);
B2=bestnest(inputnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnum*hiddennum+hiddennum+hiddennum+hiddennum*outputnum);
%% 网络权值赋值
Wjk=reshape(w1,hiddennum,inputnum);Wjk_1=Wjk;Wjk_2=Wjk_1;
Wij=reshape(w2,outputnum,hiddennum);Wij_1=Wij;Wij_2=Wij_1;
a=B1;a_1=a;a_2=a_1;
b=B2;b_1=b;b_2=b_1;
%% 最优学习速率
lr1=0.01; %学习概率
lr2=0.001; %学习概率
%%
maxgen=100; %迭代次数
%节点初始化
y=zeros(1,outputnum);
net=zeros(1,hiddennum);
net_ab=zeros(1,hiddennum);
%权值学习增量初始化
d_Wjk=zeros(hiddennum,inputnum);
d_Wij=zeros(outputnum,hiddennum);
d_a=zeros(1,hiddennum);
d_b=zeros(1,hiddennum);
error=zeros(1,maxgen);
%% 网络训练
for i=1:maxgen
% 循环训练
for kk=1:size(input,1)%1行2列。
x=inputn(kk,:);%第kk行的全部元素。
yqw=outputn(kk,:);
for j=1:hiddennum
%%
% *BOLD TEXT*
for k=1:inputnum
net(j)=net(j)+Wjk(j,k)*x(k);%输入
net_ab(j)=(net(j)-b(j))/a(j);%隐含
end
temp=mymorlet(net_ab(j));
for k=1:outputnum
y=y+Wij(k,j)*temp; %小波函数输出
end
end
%计算误差和
error(i)=sum(abs(yqw-y));
%权值调整
for j=1:hiddennum
%计算d_Wij
temp=mymorlet(net_ab(j));
for k=1:outputnum
d_Wij(k,j)=d_Wij(k,j)-(yqw(k)-y(k))*temp;%权值增量调整
end
%计算d_Wjk
temp=d_mymorlet(net_ab(j));
for k=1:inputnum
for l=1:outputnum
d_Wjk(j,k)=d_Wjk(j,k)+(yqw(l)-y(l))*Wij(l,j) ;
end
d_Wjk(j,k)=-d_Wjk(j,k)*temp*x(k)/a(j);
end
%计算d_b
for k=1:outputnum
d_b(j)=d_b(j)+(yqw(k)-y(k))*Wij(k,j);
end
d_b(j)=d_b(j)*temp/a(j);
%计算d_a
for k=1:outputnum
d_a(j)=d_a(j)+(yqw(k)-y(k))*Wij(k,j);
end
d_a(j)=d_a(j)*temp*((net(j)-b(j))/b(j))/a(j);
end
%权值参数更新 权值调整
Wij=Wij-lr1*d_Wij;
Wjk=Wjk-lr1*d_Wjk;
b=b-lr2*d_b;
a=a-lr2*d_a;
y=zeros(1,outputnum);
net=zeros(1,hiddennum);
net_ab=zeros(1,hiddennum);
%权值学习增量初始化
d_Wjk=zeros(hiddennum,inputnum);
d_Wij=zeros(outputnum,hiddennum);
d_a=zeros(1,hiddennum);
d_b=zeros(1,hiddennum);
end
end
for KK=1:40
tic
input_test=Y(36*KK-35:36*KK-24)';%测试集输入
output_test=Y(36*KK-23:36*KK)';%测试集输出
yuce=zeros(1,24);
x_test=input_test;
xx=mapminmax('apply',x_test',inputps);
x_test1=xx';
for j=1:1:hiddennum
for k=1:1:inputnum
net(j)=net(j)+Wjk(j,k)*x_test1(k);
net_ab(j)=(net(j)-b(j))/a(j);
end
temp=mymorlet(net_ab(j));
for k=1:outputnum
y(k)=y(k)+Wij(k,j)*temp;
end
end
yuce(1)=mapminmax('reverse',y,outputps);
%网络预测
for i=2:24
x_test=[x_test(2:12) yuce(i-1)];
xxx=mapminmax('apply',x_test',inputps);
x_test1=xxx';
for j=1:1:hiddennum
for k=1:1:inputnum
net(j)=net(j)+Wjk(j,k)*x_test1(k);
net_ab(j)=(net(j)-b(j))/a(j);
end
temp=mymorlet(net_ab(j));
for k=1:outputnum
y(k)=y(k)+Wij(k,j)*temp ;
end
end
yuce(i)=mapminmax('reverse',y,outputps);
y=zeros(1,outputnum);%为什么还要初始化？？？
net=zeros(1,hiddennum);
net_ab=zeros(1,hiddennum);
end
predict(KK,:)=yuce;
MAPE(KK,:)=(abs((yuce-output_test)./output_test))*100;
MSE(KK,:)=((output_test-yuce).^2);
MAE(KK,:)=(abs(output_test-yuce));
AE(KK,:)=(output_test-yuce);
toc
end
toc
%MAPE=mean(mape);disp(['预测值相对误差百分比为 ', num2str(MAPE)])
%MSE=mean(mse);disp(['预测值均方误差为 ', num2str(MSE)])
%MAE=mean(mae);disp(['预测值误差绝对值', num2str(MAE)])
%AE=mean(ae);disp(['预测值误差', num2str(AE)])
%ZH=[MAPE;MSE;MAE;AE]
评论0