%% 基于支持向量机的风电功率预测程序
%% 清空环境
clc
clear
tic
load TPP; %载入数据
x=tpp(533:893,:);
n=3;
for i=1:length(x)-n-1
for j=i:i+n
input(i,j-i+1)=x(j,1);
end
end
output=x(n+2:length(x),1);
input=input';
output=output';
% 训练集
p_train = input(:,1:150)';
t_train = output(:,1:150)';
% 测试集
%p_test = input(:,289:end)';
%t_test = output(:,289:end)';
% p_test = input(:,593:end)';
% t_test = output(:,593:end)';
p_test = input(:,683:end)';
t_test = output(:,683:end)';
%
%% 数据归一化
% 训练集
[pn_train,inputps] = mapminmax(p_train');
pn_train = pn_train';
pn_test = mapminmax('apply',p_test',inputps);
pn_test = pn_test';
% 测试集
[tn_train,outputps] = mapminmax(t_train');
tn_train = tn_train';
tn_test = mapminmax('apply',t_test',outputps);
tn_test = tn_test';
%% SVM模型创建/训练
% % 寻找最佳c参数/g参数
% [c,g] = meshgrid(-10:0.5:10,-10:0.5:10);
% [m,n] = size(c);
% cg = zeros(m,n);
% eps = 10^(-4);
% v = 5;
% bestc = 0;
% bestg = 0;
% error = Inf;
% for i=1:m
% for j = 1:n
% cmd = ['-v ',num2str(v),' -t 2',' -c ',num2str(2^c(i,j)),' -g ',num2str(2^g(i,j) ),' -s 3 -p 0.1'];
% cg(i,j) = svmtrain(tn_train,pn_train,cmd);
% if cg(i,j) < error
% error = cg(i,j);
% bestc = 2^c(i,j);
% bestg = 2^g(i,j);
% end
% if abs(cg(i,j) - error) <= eps && bestc > 2^c(i,j)
% error = cg(i,j);
% bestc = 2^c(i,j);
% bestg = 2^g(i,j);
% end
% end
% end
% % 创建/训练SVM
% cmd = [' -t 2',' -c ',num2str(bestc),' -g ',num2str(bestg),' -s 3 -p 0.01'];
cmd = [' -t 2',' -c ',num2str(50),' -g ',num2str(0.2),' -s 3 -p 0.01'];
model = svmtrain(tn_train,pn_train,cmd);
%
%% SVM仿真预测
n=3;
nn=16; %nn为滚动的步数
for i=1:96
P_test=pn_test(i,:);
for j=1:nn
t_sim_svm=svmpredict(0,P_test,model);
for k=1:3
P_test(k)=P_test(k+1);
end
P_test(4)=t_sim_svm;
T_sim_svm1(i,j)= mapminmax('reverse',t_sim_svm,outputps);
end
end
% [Predict_1,error_1] = svmpredict(tn_train,pn_train,model);
%
%准确率
% zq(l,:)=1-sqrt(sum(((X2(l,:)-X1(l,:))/(400500)).^2)./(16));
% %计算评价指标
cx=T_sim_svm1;
y=t_test;
for j=1:96
X1(j,:)=cx(j,:);
X2(j,:)=y(j:j+15);
for jj=1:16
%平均相对百分比误差
R_3(j,jj)=abs(((X2(j,jj)-X1(j,jj))./(400.5))*100);
end
end
MAPE=mean(mean(R_3,2));
%计算评价指标
for l=1:96
X1(l,:)=cx(l,:);
X2(l,:)=y(l:l+15);
t1=0;
for ll=1:16
t(l,ll)=(1-abs(((X2(l,ll)-X1(l,ll))./(400.5))));%整场装机容量400500
if(t(l,ll)>=0.85)
t1=t1+1;
end
end
%合格率
hg(l,:)=t1/16;
%准确率
zq(l,:)=1-sqrt(sum(((X2(l,:)-X1(l,:))/(400.5)).^2)./(16));
end
%全天预测结果均方根误差
qjf=sqrt(sum(sum(((X2-X1)./(400.5)).^2))/(96*16))
%日平均预测计划曲线准确率
rzq=sum(zq)/96
%日平均预测计划曲线合格率
rhg=sum(hg)/96
%作图
t=1:96;
t1=1:96;
%作预测值图
figure(1);
plot(t,t_test(1:96,1),'r-*',t,T_sim_svm1(1:96,1),'b:o')
grid on
legend('预测值','真实值')
xlabel('样本编号')
ylabel('功率/MW')
title('预测集预测结果对比')
%作准确率图
figure(2);
plot(t,zq,'-r*');
title('8月7日实时滚动预测准确率分析图');
xlabel('时点');
ylabel('准确率百分比/%');
hold on
%作合格率图
figure(3);
plot(t,hg,'-r*');
title('8月7日实时滚动预测合格率分析图');
xlabel('时点');
ylabel('合格率百分比/%');
hold on
tf=toc
% cx(1,:)
评论0