% function psobp
% BP neural network trained by PSO algorithm
% 清空环境变量
close all
clc
clear
% AllSamIn=input_train; % Add your all input data
% AllSamOut=output_train; % Add your all output data
%
% % Pre-processing data with premnmx, you can use other functions
% global minAllSamOut;
% global maxAllSamOut;
% [AllSamInn,minAllSamIn,maxAllSamIn,AllSamOutn,minAllSamOut,maxAllSamOut] = premnmx(AllSamIn,AllSamOut);
% % draw 10 percent from all samples as testing samples，the rest as training samples
% i=[10:10:1000];
% TestSamIn=[];
% TestSamOut=[];
% for j=1:100
% TestSamIn=[TestSamIn,AllSamInn(:,i(j))];
% TestSamOut=[TestSamOut,AllSamOutn(:,i(j))];
% end
% TargetOfTestSam=...; % add reall output of testing samples
% TrainSamIn=AllSamInn;
% TrainSamOut=AllSamOutn;
% TrainSamIn(:,i)=[];
% TrainSamOut(:,i)=[];
% % Evaluating Sample
% EvaSamIn=...
% EvaSamInn=tramnmx(EvaSamIn,minAllSamIn,maxAllSamIn); % preprocessing
%读取数据
load traindata1011 A O
load goontest inputtest_may16 outputtest_may16
%训练数据和预测数据
input_train=A(1:360,:)';
input_test=inputtest_may16(1:24,:)';
output_train=O(1:360)';
output_test=outputtest_may16(1:24)';
global minAllSamOut;
global maxAllSamOut;
%%premnmx用mapminmax代替
[AllSamInn,minAllSamIn,maxAllSamIn,AllSamOutn,minAllSamOut,maxAllSamOut]=premnmx(input_train,output_train);%归一化函数
%归一化后输入向量（每一列代表一个输入变量），最大值，最小值；归一化后输出向量，最大值，最小值
% Evaluating Sample
EvaSamIn=input_test;
EvaSamInn=tramnmx(EvaSamIn,minAllSamIn,maxAllSamIn); % preprocessing
% TargetOfTestSam=output_test; % add reall output of testing samples
%训练数据
global Ptrain;
Ptrain = input_train;
global Ttrain;
Ttrain = output_train;
% Ptest = input_train(:,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;
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=40;
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]之间的随机数
%X是一个N*D矩阵，粒子的初始位置
%Initialize velocities of particles
V=m+(n-m)*rand(N,D,1);
%V是一个N*D矩阵，粒子的初始速度
global fvrec;
MinFit=[];
BestFit=[];
%Function to be minimized, performance function,i.e.,mse of net work 神经网络建立
global net;
net=newff(minmax(Ptrain),[hiddennum,outdim],{'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
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);
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));
%V(:,:,2)=cf*(W(1)*V(:,:,1)+c1*rand*(pbest(:,:,1)-X(:,:,1))+c2*rand*(G(:,:,1)-X(:,:,1)));
%V(:,:,2)=cf*(V(:,:,1)+c1*rand*(pbest(:,:,1)-X(:,:,1))+c2*rand*(G(:,:,1)-X(:,:,1)));
% 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);
%******************************************************
for j=2:itmax
disp('Iteration and Current Best Fitness')
disp(j-1)
disp(B(1,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));
MinFit=[MinFit C];
BestFit=[BestFit min(MinFit)];
L(:,1,j)=fitness(:,1,j);
B(1,1,j)=C;
gbest(1,:,j)=X(I,:,j);
[C,I]=min(B(1,1,:));
% 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,:));
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);%反归一化
error=EvaSamOut-output_test;
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',11)
ylabel('风电功率','fontsize',11)
xlabel('时间','fontsize',11)
% 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 均方根误差公式
