基于PSO-BP神经网络的风电功率预测研究论文复现。matlab程序

% 清空环境变量 close all clc clear %读取数据 load traindata1011 A O %train1011mat文件中有AOm文件，A是输入360*4（四个特征量：正弦值和风向余弦值）O是输出功率360*1 load goontest inputtest_may16 outputtest_may16%goodtest的mat文件包括in16四个特征量24*4和out16输出功率24*1 %训练数据和预测数据 input_train=A(1:360,:)';%A的转置矩阵 input_test=inputtest_may16(1:24,:)';%表格中输入横的 output_train=O(1:360)'; output_test=outputtest_may16(1:24)'; global minAllSamOut;%设置全局变量，子函数中也可以直接引用 global maxAllSamOut;%输出功率的最大（1441.6)和最小值 %归一化 %2345变成-1，-0.33,0.33,1; [AllSamInn,minAllSamIn,maxAllSamIn,AllSamOutn,minAllSamOut,maxAllSamOut]=premnmx(input_train,output_train); %allsaminn是inputtrain360的归一化结果（-1，1），minsamin是inputtrain每行的最小值（四个特征量各一个） % Evaluating Sample %网络所用的新数据和经过归一化样本数剧进行相同的预处理（先用premnmx才可用tramnmx） EvaSamIn=input_test;%evasamin是变换前的输入数据，inn是变换后的输入数据 EvaSamInn=tramnmx(EvaSamIn,minAllSamIn,maxAllSamIn); % preprocessing % TargetOfTestSam=output_test; % add real output of testing samples %目标输出和实际outputtest（24输出，表格横的）一致 global Ptrain;%定义全局变量4*360 Ptrain = input_train; global Ttrain; Ttrain = output_train; %Ptest = input_train(:,350:360);%4*11double（350-360的数据） % Ttest = output_train(:,350:360); % Initialize BPN parameters设置神经网络参数 global indim; indim=4;%输入层神经元个数 global hiddennum; hiddennum=5;%隐藏层 global outdim; outdim=1;%输出层 % Initialize PSO parameters设置微粒群参数 vmax=0.5; % Maximum velocity速度上限 minerr=0.001; % Minimum error目标误差 wmax=0.90;%典型权值，初始惯性权值 wmin=0.30;%惯性因子的引入可以调整全局和局部的搜索能力 global itmax; %Maximum iteration number最大迭代次数 itmax=200; c1=2;%加速度因子，为非负常数，c1作用使粒子朝着自身最优位置进行移动 c2=2;%c2是使粒子向群体中最优位置移动，通常的取值为2、2 for iter=1:itmax W(iter)=wmax-((wmax-wmin)/itmax)*iter; % weight declining linearly权值随着迭代次数线性递减以保证收敛 end %通常采用线性递减权值（Linearly Decreasing weight，LDW）策略（公式3下面）获取动态w % 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=40; global D; % length of particle每个微粒的维数 D=(indim+1)*hiddennum+(hiddennum+1)*outdim; %pso算法第一步，随机初始化每个解的速度和位置 % Initialize positions of particles初始化微粒位置 rand('state',sum(100*clock));%根据当前时间设置初始状态保证非重复随机数，新版用rng函数 X=a+(b-a)*rand(N,D,1); %取值范围[-1,1] rand * 2 - 1 ,rand 产生[0,1]之间的随机数 %Initialize velocities of particles初始化微粒速度 V=m+(n-m)*rand(N,D,1); global fvrec; MinFit=[];%定义一个空的数组 BestFit=[]; %微粒群的更新迭代 global net; net=newff(minmax(Ptrain),[hiddennum,outdim],{'tansig','purelin'}); %旧用法 %newff第一个变量，用minmax来设定输入特征值的范围，即每行最大最小值 %newff第二个变量，设定隐含层和输出层神经元的数目 %newff第三个变量，设定转移函数tansig（反正切），通常在输出层选一个线性函数purelin %适应度函数部分，第一次迭代 %计算适应值 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全局最优适应值，B存储当前代的最优适应值 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,:,1);%G便于速度更新运算，函数格式统一 end for i=1:N; pbest(i,:,1)=X(i,:,1);%因为是第一代，所以当前位置为历史最优位置 end V(:,:,2)=W(1)*V(:,:,1)+c1*rand*(pbest(:,:,1)-X(:,:,1))+c2*rand*(G(:,:,1)-X(:,:,1)); %pso标准形式（公式3）,更新速度，V矩阵中的第一维、第二维和第三维的第一个元素。 % limits velocity of particles by vmax判断速度是否超过限制 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);%位置更新（公式1） %第二次迭代（最后一次） for j=2:itmax disp('Iteration and Current Best Fitness')%迭代次数和全局最佳适应值 disp(j-1) disp(B(1,1,j-1))%j-1代全局最优适应值 % Calculation of new positions fitness=fitcal(X,net,indim,hiddennum,outdim,D,Ptrain,Ttrain,minAllSamOut,maxAllSamOut); fvrec(:,1,j)=fitness(:,1,j); %[maxC,maxI]=max(fitness(:,1,j)); %MaxFit=[MaxFit maxC]; %MeanFit=[MeanFit mean(fitness(:,1,j))]; [C,I]=min(fitness(:,1,j));%第j代的最优适应值和最优微粒序号 MinFit=[MinFit C]; BestFit=[BestFit min(MinFit)]; L(:,1,j)=fitness(:,1,j);%第j代每个微粒的适应值 B(1,1,j)=C;%第j代全局最优适应值 gbest(1,:,j)=X(I,:,j);%第j代全局最优微粒的位置 [C,I]=min(B(1,1,:));%所有代的全局最优适应值赋给C、I % keep gbest is the best particle of all have occured计算历史全局最优位置 if B(1,1,j)<=C gbest(1,:,j)=gbest(1,:,j); else gbest(1,:,j)=gbest(1,:,I); end if C<=minerr, break, end%若满足均方误差条件，记录最优位置，停止迭代 %Matrix composed of gbest vector %将最优微粒（最优权值和阈值）赋给神经网络 if j>=itmax, break, end%若超过最大迭代次数，退出 for p=1:N G(p,:,j)=gbest(1,:,j); end %计算微粒历史最优位置 for i=1:N; [C,I]=min(L(i,1,:));%最优适应值赋给C，代数赋给I if L(i,1,j)<=C pbest(i,:,j)=X(i,:,j); else pbest(i,:,j)=X(i,:,I); end end V(:,:,j+1)=W(j)*V(:,:,j)+c1*rand*(pbest(:,:,j)-X(:,:,j))+c2*rand*(G(:,:,j)-X(:,:,j)); %V(:,:,j+1)=cf*(W(j)*V(:,:,j)+c1*rand*(pbest(:,:,j)-X(:,:,j))+c2*rand*(G(:,:,j)-X(:,:,j))); %V(:,:,j+1)=cf*(V(:,:,j)+c1*rand*(pbest(:,:,j)-X(:,:,j))+c2*rand*(G(:,:,j)-X(:,:,j))); for ni=1:N for di=1:D if V(ni,di,j+1)>vmax V(ni,di,j+1)=vmax; elseif V(ni,di,j+1)<-vmax V(ni,di,j+1)=-vmax; else V(ni,di,j+1)=V(ni,di,j+1); end end end X(:,:,j+1)=X(:,:,j)+V(:,:,j+1); end disp('Iteration and Current Best Fitness') disp(j) disp(B(1,1,j)) 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,j); end for r=1:outdim x2lw(r,:)=gbest(1,(indim*hiddennum+1):(indim*hiddennum+hiddennum),j); end x2b=gbest(1,((indim+1)*hiddennum+1):D,j); 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; %% BP网络训练 %网络进化参数 net.trainParam.epochs=100; net.trainParam.lr=0.1; net.trainParam.goal=0.00001; net.trainParam.show=100; net.trainParam.showWindow=0; %网络训练 [net,per2]=train(net,AllSamInn,AllSamOutn); % nettesterr=mse(sim(net,Ptest)-Ttest); % testsamout = postmnmx(sim(net,Ptest),minAllSamOut,maxAllSamOut); % realtesterr=mse(testsamout-TargetOfTestSam) EvaSamOutn = sim(net,EvaSamInn); EvaSamOut = postmnmx(EvaSamOutn,minAllSamOut,maxAllSamOut);%反归一化 error1=(EvaSamOut-output_test)./output_test*100; errormape=(EvaSamOut-output_test)./output_test; figure(1) grid hold on plot(log(BestFit),'r'); title(['PSO适应度曲线 ' '最优代数＝' I]); xlabel('进化代数');ylabel('适应度'); legend('平均适应度','最佳适应度'); disp('适应度 变量'); figure(2) grid plot(EvaSamOut,':og') hold on plot(output_test,'-*'); legend('预测输出','期望输出') title('PSO-BP网络预测输出','fontsize',12) ylabel('函数输出','fontsize',12) xlabel('样本','fontsize',12) figure(3) plot(error1,'-*'); title('神经网络预测误差百分比') figure(4) hist(errormape); title('PSO-BP神经网络预测误差频率分布直方图'); ylabel('频率（次）','fontsize',12) xlabel('相对误差','fontsize',12) MAE=(sum(abs(errormape)))/24 %绝对平均误差 RMSE=sqrt((sum(errormape.^2))/24)%RMSE 均方根误差公式