?? 20041006.m
字號:
function psobp
% BP neural network trained by PSO algorithm
clc
clear all
% generate training samples according to eutrophication evaluation standard
% of surface water
st=cputime;
AllSamIn=...
[ 70.00 60.00 60.00 60.00 70.00
70.00 65.00 65.00 65.00 70.00
70.00 75.00 70.00 70.00 70.00
70.00 75.00 70.00 65.00 70.00
70.00 75.00 65.00 65.00 68.00
70.00 75.00 65.00 65.00 68.00
70.00 75.00 65.00 65.00 68.00
70.00 75.00 65.00 65.00 68.00
70.00 75.00 65.00 65.00 70.00
70.00 65.00 65.00 60.00 68.00
70.00 65.00 65.00 60.00 68.00
70.00 65.00 65.00 60.00 68.00
70.00 65.00 65.00 60.00 68.00
70.00 65.00 65.00 60.00 68.00
70.00 65.00 65.00 60.00 68.00
70.00 65.00 65.00 60.00 68.00
70.00 65.00 65.00 60.00 68.00
70.00 65.00 65.00 60.00 68.00
70.00 65.00 65.00 60.00 68.00
70.00 65.00 65.00 60.00 68.00
70.00 65.00 65.00 60.00 68.00
70.00 65.00 65.00 60.00 68.00
70.00 65.00 65.00 60.00 68.00
70.00 65.00 65.00 60.00 68.00
70.00 75.00 65.00 65.00 68.00
70.00 75.00 65.00 65.00 68.00
70.00 75.00 65.00 65.00 68.00
70.00 75.00 65.00 65.00 68.00
135.00 120.00 135.00 100.00 105.00
135.00 120.00 135.00 100.00 105.00
135.00 120.00 135.00 100.00 105.00
135.00 130.00 135.00 100.00 105.00
135.00 135.00 135.00 90.00 105.00
135.00 120.00 135.00 105.00 105.00
130.00 120.00 85.00 85.00 90.00
110.00 120.00 80.00 80.00 90.00
120.00 120.00 80.00 80.00 90.00
120.00 120.00 80.00 80.00 90.00
120.00 100.00 80.00 80.00 90.00
135.00 135.00 120.00 105.00 105.00
135.00 130.00 130.00 105.00 105.00
135.00 130.00 130.00 105.00 105.00
135.00 135.00 135.00 105.00 105.00
135.00 135.00 135.00 105.00 105.00
110.00 110.00 110.00 85.00 85.00
110.00 110.00 110.00 85.00 85.00
110.00 110.00 110.00 85.00 85.00
110.00 110.00 110.00 85.00 85.00
100.00 100.00 80.00 85.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
100.00 85.00 80.00 80.00 85.00
100.00 85.00 80.00 80.00 85.00
100.00 100.00 80.00 80.00 85.00
135.00 135.00 135.00 135.00 135.00
135.00 135.00 135.00 135.00 135.00
135.00 135.00 135.00 135.00 135.00
135.00 135.00 135.00 135.00 135.00
135.00 135.00 135.00 135.00 135.00
135.00 135.00 135.00 135.00 135.00
135.00 135.00 135.00 135.00 135.00
135.00 135.00 135.00 135.00 135.00
125.00 135.00 135.00 135.00 135.00
125.00 130.00 130.00 130.00 135.00
110.00 130.00 130.00 130.00 130.00
125.00 130.00 130.00 130.00 135.00
125.00 120.00 115.00 115.00 120.00
125.00 120.00 115.00 115.00 120.00
125.00 120.00 115.00 115.00 120.00
125.00 120.00 115.00 115.00 100.00
70.00 75.00 60.00 60.00 66.00
70.00 75.00 60.00 60.00 66.00
70.00 75.00 60.00 60.00 66.00
70.00 75.00 60.00 60.00 66.00
]';
% 產(chǎn)生噪聲，加給輸出
% NoiseVar=0.1;
% rand('state',sum(100*clock));
AllSamOut=[ 70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
135.00
135.00
135.00
135.00
135.00
135.00
130.00
110.00
120.00
120.00
120.00
135.00
135.00
135.00
135.00
135.00
110.00
110.00
110.00
110.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
100.00
135.00
135.00
135.00
135.00
135.00
135.00
135.00
135.00
125.00
125.00
110.00
125.00
125.00
125.00
125.00
125.00
70.00
70.00
70.00
70.00
]'; %+NoiseVar*randn(1,1000);
% 預(yù)處理，壓縮到[-1,+1]
global minAllSamOut;
global maxAllSamOut;
[AllSamInn,minAllSamIn,maxAllSamIn,AllSamOutn,minAllSamOut,maxAllSamOut] = premnmx(AllSamIn,AllSamOut);
% 從總樣本中抽取10％的樣本作為測試樣本，其余作為訓(xùn)練樣本
TestSamIn=AllSamIn;
RealTestSamOut=[ 70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
70.00
68.00
68.00
70.00
70.00
70.00
70.00
70.00
70.00
68.00
70.00
70.00
68.00
68.00
70.00
70.00
70.00
70.00
70.00
133.00
133.00
133.00
135.00
135.00
135.00
133.00
110.00
118.00
118.00
118.00
130.00
135.00
135.00
135.00
135.00
110.00
110.00
110.00
110.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
98.00
135.00
135.00
135.00
135.00
135.00
135.00
135.00
135.00
135.00
135.00
135.00
135.00
135.00
135.00
130.00
123.00
68.00
68.00
68.00
68.00
]';
[TestSamInn,minTestSamIn,maxTestSamIn,RealTestSamOutn,minRealTestSamOut,maxRealTestSamOut] = premnmx(TestSamIn,RealTestSamOut);
% 訓(xùn)練樣本
TrainSamIn=AllSamIn;
TrainSamOut=AllSamOut;
[TrainSamInn,minTrainSamIn,maxTrainSamIn,TrainSamOutn,minTrainSamOut,maxTrainSamOut] = premnmx(TrainSamIn,TrainSamOut);
% Evaluation Sample
%EvaSamIn=...
% [6.5400 9.4000 8.2200 5.9100 6.1300 7.8400 8.2600 8.7900 7.2900 5.9300 4.6900 3.5100
% 0.1270 0.0390 0.0820 0.0860 0.1230 0.1370 0.0190 0.0640 0.0880 0.0500 0.0860 0.0750
% 1.2840 1.3580 1.3200 0.4500 2.2200 1.4700 1.0900 1.6600 0.8400 0.6600 0.5500 0.5200
% 3.3400 4.7200 5.0300 6.0000 7.9600 4.6600 4.0000 4.1500 4.0000 3.8100 1.6300 3.4500
%0.3500 0.4500 0.3800 0.4500 0.3000 0.4500 0.4000 0.4000 0.3000 0.3000 0.4000 0.3500];
%EvaSamInn=tramnmx(EvaSamIn,minAllSamIn,maxAllSamIn); % 預(yù)處理
%**********************************************************
% training set
global Ptrain;
Ptrain = TrainSamInn;
global Ttrain;
Ttrain = TrainSamOutn;
% testing set
Ptest = TestSamInn;
%**********************************************************
% Initialize BPN parameters
global indim;
indim=5;
global hiddennum;
hiddennum=3;
global outdim;
outdim=1;
%**********************************************************
% Initialize PSO parameters
vmax=0.5; % Maximum velocity
minerr=0.0001; % Minimum error
wmax=0.90;
wmin=0.30;
global itmax; %Maximum iteration number
itmax=100;
c1=2;
c2=1.8;
%cf=c1+c2;
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;
global fvrec;
MinFit=[];
BestFit=[];
%MaxFit=[];
%MeanFit=[];
% Initialize positions of particles
rand('state',sum(100*clock));
X=a+(b-a)*rand(N,D,1);
%Initialize velocities of particles
V=m+(n-m)*rand(N,D,1);
%**********************************************************
%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);
%[maxC,maxI]=max(fitness(:,1,1));
%MaxFit=[MaxFit maxC];
%MeanFit=[MeanFit mean(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;
%net.IW{1,1}
%net.LW{2,1}
%net.b{1}
%net.b{2}
%nettesterr=mse(sim(net,Ptest)-Ttest);
TestSamOut = sim(net,Ptest);
[a]=postmnmx(TestSamOut,minRealTestSamOut,maxRealTestSamOut);
ae=abs(a-RealTestSamOut)
mae=mean(ae)
re=(a-RealTestSamOut)./RealTestSamOut
mre=mean(abs(re))
a
%EvaSamOutn = sim(net,EvaSamInn);
%EvaSamOut = postmnmx(EvaSamOutn,minAllSamOut,maxAllSamOut);
figure(1)
grid
hold on
%plot(MaxFit,'k');
plot(log(BestFit),'--k','linewidth',2);
title('適應(yīng)度');
xlabel('迭代次數(shù)');
ylabel('fit');
%plot(log(MinFit),'b');
figure(2)
grid
hold on
plot(a,'b');
plot(RealTestSamOut,'-k','linewidth',2);
title('10月6日電價預(yù)測');
xlabel('時段');
ylabel('電價（元);
）'
figure(3)
grid
hold on
plot(ae,'k','linewidth',2);
title('絕對誤差');
xlabel('時段');
ylabel('誤差');
figure(4)
grid
hold on
plot(re,'k','linewidth',2);
title('相對誤差');
xlabel('時段');
ylabel('誤差');
et=cputime-st;
%plot(EvaSamOut,'k');
save D:
%--------------------------------------------------------------------------
%sub function for getting fitness of all paiticles in specific generation
%change particle to weight matrix of BPN,then calculate training error
function fitval = fitcal(pm,net,indim,hiddennum,outdim,D,Ptrain,Ttrain,minAllSamOut,maxAllSamOut)
[x,y,z]=size(pm);
for i=1:x
for j=1:hiddennum
x2iw(j,:)=pm(i,((j-1)*indim+1):j*indim,z);
end
for k=1:outdim
x2lw(k,:)=pm(i,(indim*hiddennum+1):(indim*hiddennum+hiddennum),z);
end
x2b=pm(i,((indim+1)*hiddennum+1):D,z);
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;
error=postmnmx(sim(net,Ptrain),minAllSamOut,maxAllSamOut)-postmnmx(Ttrain,minAllSamOut,maxAllSamOut);
fitval(i,1,z)=mse(error);
end
?? 快捷鍵說明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -