基于BP神经网络的风电功率识别，基于svm的风电功率识别，基于粒子群算法优化SVM的风电功率识别

%读取数据 clc clear close all clear global ver %训练数据和预测数据 load maydata1.mat input_train=num(n(1:m),1:8)'; output_train=num(n(1:m),9)'; input_test=num(n(m+1:end),1:8)'; output_test=num(n(m+1:end),9)'; [AllSamInn,minAllSamIn,maxAllSamIn,AllSamOutn,minAllSamOut,maxAllSamOut]=premnmx(input_train,output_train); EvaSamIn=input_test; EvaSamInn=tramnmx(EvaSamIn,minAllSamIn,maxAllSamIn); % preprocessing Ptrain = AllSamInn; Ttrain = AllSamOutn; indim=8; hiddennum=20; outdim=1; % Initialize PSO vmax=0.0151; % Maximum velocity minerr=0.001; % Minimum error wmax=0.90; wmin=0.30; % global itmax; %Maximum iteration number itmax=100; c1=2; c2=2; for iter=1:itmax W(iter)=wmax-((wmax-wmin)/itmax)*iter; % weight declining linearly end % particles are initialized between (a,b) randomly a=-1; b=1; %Between (m,n), (which can also be started from zero) m=-1; n=1; % global N; % number of particles N=20; % global D; % length of particle D=(indim+1)*hiddennum+(hiddennum+1)*outdim; % Initialize positions of particles % rand('state',sum(100*clock)); X=a+(b-a)*rand(N,D,1); %取值范围[-1,1] rand * 2 - 1 ,rand 产生[0,1]之间的随机数 %Initialize velocities of particles V=0.2*(m+(n-m)*rand(N,D,1)); % % global fvrec; MinFit=[]; BestFit=[]; net=newff(minmax(Ptrain),[hiddennum,outdim],{'logsig','tansig'},'traingdx'); fitness=fitcal(X,net,indim,hiddennum,outdim,D,Ptrain,Ttrain,minAllSamOut,maxAllSamOut); fvrec(:,1,1)=fitness(:,1,1); [C,I]=min(fitness(:,1,1)); MinFit=[MinFit C]; BestFit=[BestFit C]; L(:,1,1)=fitness(:,1,1); %record the fitness of particle of every iterations B(1,1,1)=C; %record the minimum fitness of particle gbest(1,:,1)=X(I,:,1); %the global best x in population %Matrix composed of gbest vector for p=1:N G(p,:,1)=gbest(1,:); end for ii=1:N; pbest(ii,:,1)=X(ii,:); end V(:,:,2)=W(1)*V(:,:,1)+c1*rand*(pbest(:,:,1)-X(:,:,1))+c2*rand*(G(:,:,1)-X(:,:,1)); for ni=1:N for di=1:D if V(ni,di,2)>vmax V(ni,di,2)=vmax; elseif V(ni,di,2)<-vmax V(ni,di,2)=-vmax; else V(ni,di,2)=V(ni,di,2); end end end X(:,:,2)=X(:,:,1)+V(:,:,2); for ni=1:N for di=1:D if X(ni,di,2)>1 X(ni,di,2)=1; elseif X(ni,di,2)<-1 X(ni,di,2)=-1; else X(ni,di,2)=X(ni,di,2); end end end %****************************************************** for jj=2:itmax disp('Iteration and Current Best Fitness') disp(jj-1) disp(B(1,1,jj-1)) reset =1; % reset = 1时设置为粒子群过分收敛时将其打散，如果＝1则不打散 if reset==1 bit = 1; for k=1:N bit = bit&(range(X(k,:))<0.02); end if bit==1 % bit=1时对粒子位置及速度进行随机重置 for ik = 1:N X(ik,:) = funx; % present 当前位置,随机初始化 X(ik,:) = [0.02*rand()-0.01 0.02*rand()-0.01]; % 速度初始化 end end end % Calculation of new positions fitness=fitcal(X,net,indim,hiddennum,outdim,D,Ptrain,Ttrain,minAllSamOut,maxAllSamOut); [C,I]=min(fitness(:,1,jj)); MinFit=[MinFit C]; BestFit=[BestFit min(MinFit)]; L(:,1,jj)=fitness(:,1,jj); B(1,1,jj)=C; gbest(1,:,jj)=X(I,:,jj); [C,I]=min(B(1,1,:)); % keep gbest is the best particle of all have occured if B(1,1,jj)<=C gbest(1,:,jj)=gbest(1,:,jj); else gbest(1,:,jj)=gbest(1,:,I); end if C<=minerr break end %Matrix composed of gbest vector if jj>=itmax break end for p=1:N G(p,:,jj)=gbest(1,:,jj); end for ii=1:N; [C,I]=min(L(ii,1,:)); if L(ii,1,jj)<=C pbest(ii,:,jj)=X(ii,:,jj); else pbest(ii,:,jj)=X(ii,:,I); end end V(:,:,jj+1)=W(jj)*V(:,:,jj)+c1*rand*(pbest(:,:,jj)-X(:,:,jj))+c2*rand*(G(:,:,jj)-X(:,:,jj)); for ni=1:N for di=1:D if V(ni,di,jj+1)>vmax V(ni,di,jj+1)=vmax; elseif V(ni,di,jj+1)<-vmax V(ni,di,jj+1)=-vmax; else V(ni,di,jj+1)=V(ni,di,jj+1); end end end X(:,:,jj+1)=X(:,:,jj)+V(:,:,jj+1); for ni=1:N for di=1:D if X(ni,di,jj+1)>1 X(ni,di,jj+1)=1; elseif X(ni,di,jj+1)<-1 X(ni,di,jj+1)=-1; else X(ni,di,jj+1)=X(ni,di,jj+1); end end end end disp('Iteration and Current Best Fitness') disp(jj) disp(B(1,1,jj)) disp('Global Best Fitness and Occurred Iteration') [C,I]=min(B(1,1,:)); % simulation network 网络拟合 for t=1:hiddennum x2iw(t,:)=gbest(1,((t-1)*indim+1):t*indim,jj); end for r=1:outdim x2lw(r,:)=gbest(1,(indim*hiddennum+1):(indim*hiddennum+hiddennum),jj); end x2b=gbest(1,((indim+1)*hiddennum+1):D,jj); x2b1=x2b(1:hiddennum).'; x2b2=x2b(hiddennum+1:hiddennum+outdim).'; net.IW{1,1}=x2iw; net.LW{2,1}=x2lw; net.b{1}=x2b1; net.b{2}=x2b2; [BPoutput1,error1,net] = bpp(); %% BP网络训练 %网络进化参数 net.trainParam.epochs=5000; net.trainParam.lr=0.1; net.trainParam.goal=0.0001; % net.trainParam.show=100; % net.trainParam.showWindow=1; net=newff(AllSamInn,AllSamOutn,[hiddennum,outdim],{'logsig','tansig'},'traingdx'); tic %网络训练 net=train(net,AllSamInn,AllSamOutn); toc EvaSamOutn = sim(net,EvaSamInn); EvaSamOut = postmnmx(EvaSamOutn,minAllSamOut,maxAllSamOut);%反归一化 error=EvaSamOut-output_test; errormape=(EvaSamOut-output_test)./output_test; p1 = sum(abs(error1))/33; [ EvaSamOut2 , BestFit1 ] = apsofun( ); figure(1) grid hold on plot((BestFit),'r-'); hold on plot((BestFit1),'b-'); legend('pso','自适应pso') title(['粒子群算法优化bp ' '最优代数＝' I]); xlabel('进化代数');ylabel('误差'); disp('适应度变量'); set(gca,'fontsize',12) figure grid plot((BestFit),'r'); title(['粒子群算法优化bp ' '最优代数＝' I]); xlabel('进化代数');ylabel('误差'); disp('适应度变量'); set(gca,'fontsize',12) figure grid plot(EvaSamOut,':og') hold on plot(BPoutput1,':ob') hold on plot(output_test,'-*r'); legend('粒子群优化BP预测输出','BP预测输出','期望输出') title('粒子群优化BP网络预测输出','fontsize',12) ylabel('风电功率') xlabel('样本','fontsize',12) set(gca,'fontsize',12) figure grid plot(EvaSamOut2,':og') hold on plot(BPoutput1,':ob') hold on plot(output_test,'-*r'); legend('自适应粒子群优化BP预测输出','BP预测输出','期望输出') title('粒子群优化BP网络预测输出','fontsize',12) ylabel('风电功率') xlabel('样本','fontsize',12) set(gca,'fontsize',12) figure grid plot(EvaSamOut,':og') hold on hold on plot(output_test,'-*r'); ylabel('风电功率') legend('粒子群优化BP预测输出','期望输出') title('粒子群优化BP网络预测输出','fontsize',12) xlabel('样本','fontsize',12) set(gca,'fontsize',12) figure plot(BPoutput1,':og') hold on plot(output_test,'-*'); legend('预测输出','期望输出') title('·BP神经网络','fontsize',12) ylabel('风电功率') xlabel('样本','fontsize',12) set(gca,'fontsize',12) %预测误差 error1=BPoutput1-output_test; junfanggen = mse(BPoutput1-output_test); figure plot(error1,'-k*') title('BP网络预测误差','fontsize',12) ylabel('误差','fontsize',12) xlabel('样本','fontsize',12) %axis([1 2500 -0.5 0.5]) set(gca,'fontsize',12) [MSE0,RMSE0, MBE0, MAE0 ] =MSE_RMSE_MBE_MAE(output_test,BPoutput1);%均方误差 均方根误差 平均绝对误差 平均偏差 R0 = R_2(output_test,BPoutput1);%拟合优度 can0 = [MSE0,RMSE0, MBE0, MAE0 R0]; [MSE1,RMSE1, MBE1, MAE1 ] =MSE_RMSE_MBE_MAE(output_test,EvaSamOut);%均方误差 均方根误差 平均绝对误差 平均偏差 R1 = R_2(output_test,EvaSamOut);%拟合优度 can1 = [MSE1,RMSE1, MBE1, MAE1 R1]; [MSE2,RMSE2, MBE2, MAE2 ] =MSE_RMSE_MBE_MAE(ou