粒子群算法优化bp matlab程序,MATLAB中粒子群算法优化BP网络模权值和阈值

程序是论坛里的，我自己修改了一下，可是有问题，麻烦帮我看一下

allAllSamIn=z; % 添加所有输入样本

AllSamOut=y; % 添加所有输出样本

global minAllSamOut;

global maxAllSamOut;%建立全局变量

[AllSamInn,minAllSamIn,maxAllSamIn,AllSamOutn,minAllSamOut,maxAllSamOut] = premnmx(AllSamIn,AllSamOut);%归一化

i=[10:10:500];%取实验中500个数据的10%作为测试样本

TestSamIn=[];

TestSamOut=[];

for j=1:50

TestSamIn=[TestSamIn,AllSamInn(:,i(j))];

TestSamOut=[TestSamOut,AllSamOutn(:,i(j))];

end

TargetOfTestSam=[143.54 143.3 146.12 150.15 149.62 149.16 149.08 147.83 145.84 144.09 142.73 141.99 140.54 146.97 145.44 142.23 141.55 145.05 146.16 145.56;

26.2 27.2 23.5 21.6 18.3 5.7 6.5 -1 5.1 13.1 14.7 24.6 26.8 27.8 32 22.6 11.5 10.8 1.8 5.1]; % 添加目标测试样本

TrainSamIn=AllSamInn;

TrainSamOut=AllSamOutn;

TrainSamIn(:,i)=[];

TrainSamOut(:,i)=[];

% Evaluating Sample

EvaSamIn=[2.77 1.21 1.17 -0.48 -0.11 -0.21 1.19 2.2 2.46 2.48 3.98 4.07 2.52 0.39 -0.48 -0.9 -0.34 1.58 1.79 1.57]；%添加估计样本

EvaSamInn=tramnmx(EvaSamIn,minAllSamIn,maxAllSamIn); % preprocessing

global Ptrain;

Ptrain = TrainSamIn;

global Ttrain;

Ttrain = TrainSamOut;

Ptest = TestSamIn;

Ttest = TestSamOut;

% Initialize BPN parametersglobal indim;

indim=5;%初始化5层BP网络

global hiddennum;

hiddennum=3;

global outdim;

outdim=1;

% 初始化 PSO 参数

vmax=0.5; % Maximum velocity

minerr=0.001; % Minimum error

wmax=0.90;wmin=0.30;

global itmax; %Maximum iteration number

itmax=300;

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);

%Initialize velocities of particles

V=m+(n-m)*rand(N,D,1);

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)

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

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*(W(j)*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)

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;

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);

figure(1)

grid

hold on

plot(log(BestFit),'r');

figure(2)

grid

hold on

plot(EvaSamOut,'k');

save er net nettesterr realtesterr B fvrec EvaSamOut

%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=sim(net,Ptrain)-Ttrain;

fitval(i,1,z)=mse(error);

end