clear
vm=xlsread('C:\Users\Administrator\Desktop\123.xls'); %取原始风速数据200个
%用db3小波进行3层分解
[c,l]=wavedec(vm,3,'db4');
%重构第1-5层逼近系数.
a3 = wrcoef('a',c,l,'db4',3);
%重构第1-3层细节系数
d3 = wrcoef('d',c,l,'db4',3);
d2 = wrcoef('d',c,l,'db4',2);
d1 = wrcoef('d',c,l,'db4',1);
d=[a3,d3,d2,d1];
totalnumber=200;
totalnumber2=200;
%totalnumber2=8560;
v=d(1:totalnumber,1:4);
v1=d(totalnumber+1:totalnumber+totalnumber2,1:4);
for m=1:4
%1标准化处理
average_v=mean(d(1:end,m)); %平均值
var_v=var(d(1:end,m)); %方差
var_v_s=sqrt(var_v); %标准差
for i=1:totalnumber
DX(i,1)=(v(i,m)-average_v)/var_v_s;
end
vs=DX;
%2标准化处理
for i=1:totalnumber2
DX1(i,1)=(v1(i,m)-average_v)/var_v_s;
end
vs1=DX1(1:totalnumber2,1);
z=iddata(vs); %格式转换
z1=iddata(vs1);
for p=1:10
q=p-1;
n=armax(z,[p q]);
AIC=aic(n);
end
for k=1:10
if AIC(k)==min(AIC)
p_final=k;
q_final=p_final-1;
break
end
end
m_final=armax(z,[p_final q_final]); %确定最终模型
%预测过程
p=predict(m_final,z1,1);
x=1:1:totalnumber2;
po=p.outputdata;
for i=1:totalnumber2
wind_speed(i,m)=po(i)*var_v_s+average_v;
end
m=m+1;
end
windpower=wind_speed(1:end,1)+wind_speed(1:end,2)+wind_speed(1:end,3)+wind_speed(1:end,4);
figure(1)
plot(x,vm(totalnumber+1:totalnumber+totalnumber2),'b',x,windpower,'r');
xlabel('t(1h)','FontName','Times New Roman','FontSize',14);
ylabel('Windpower(kW)','FontName','Times New Roman','FontSize',14);
legend('Simulated data','Predicted data');
%RMSE——均方根误差
P_actual=vm(totalnumber+1:totalnumber+totalnumber2);
P_predict=windpower;
Cap=1000;
a=(P_actual-P_predict)/Cap;
z=0;
for i=1:totalnumber2
z=z+a(i)^2;
end
RMSE=sqrt(z/totalnumber2)*100
%MAE——平均绝对误差
a=abs(P_actual-P_predict)/Cap;
z=0;
for i=1:totalnumber2
z=z+a(i);
end
MAE=z/totalnumber2*100
%xlswrite('1yue.xlsx',P_predict,1,'D201');
error=P_actual-P_predict;
figure;
[f,xc]=ecdf(error);
ecdfhist(f,xc,50);
xlabel('预测误差');
ylabel('f(x)');
