psobp.m

用于控制的神经网络程序

% BP neural network trained by PSO algorithm
% Copyright by Deng Da-Peng @ 2005
% Email: [e-mail]
% You can change and distribute this code freely for academic usage
% Business usage is strictly prohibited
clc
clear all

AllSamIn=1:2:200;
n1=AllSamIn;
AllSamOut=sin(n1*0.1);
% 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=[1:100];
TestSamIn=[1:2:200];
TestSamOut=sin(n1*0.1);

TargetOfTestSam= sin(n1*0.1);
TrainSamIn=AllSamInn;
TrainSamOut=AllSamOutn;

% Evaluating Sample
EvaSamIn=1:2:200;
EvaSamInn=tramnmx(EvaSamIn,minAllSamIn,maxAllSamIn); % preprocessing

global Ptrain;
Ptrain = TrainSamIn;
global Ttrain;
Ttrain = TrainSamOut;

Ptest = TestSamIn;
Ttest = TestSamOut;

% Initialize BPN parameters
global indim;
indim=5;
global hiddennum;
hiddennum=3;
global outdim;
outdim=1;

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=[];
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(:,1);
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);
Srealtesterr=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');

%sub function for getting fitness of all paiticles in specific generation
%change particle to weight matrix of BPN,then calculate training error 

 




%%####################################################################
%%#### Particle swarm optimization
%%#### With linkage operator
%%#### Deepak devicharan july 2003
%%####################################################################

%%## to apply this to different equations do the following
%%## generate initial particles in a search space close to actual soln
%%## fool around with no of iterations, no of particles, learning rates

%%## for a truly generic PSO do the following
%%## increase the number of particles , increase the variance
%%## i.e let the particles cover a larger area of the search space
%%## then fool around as always with the above thins

%declare the parameters of the optimization

max_iterations = 1000;
no_of_particles = 50;
dimensions = 1;

delta_min = -0.003;
delta_max = 0.003;

c1 = 1.3;
c2 = 1.3;

%initialise the particles and teir velocity components

for count_x = 1:no_of_particles
for count_y = 1:dimensions
particle_position(count_x,count_y) = rand*10;
particle_velocity(count_x,count_y) = rand;
p_best(count_x,count_y) = particle_position(count_x,count_y);
end
end

%initialize the p_best_fitness array
for count = 1:no_of_particles
p_best_fitness(count) = -1000;
end

%particle_position
%particle_velocity

%main particle swrm routine
for count = 1:max_iterations

%find the fitness of each particle
%change fitness function as per equation requiresd and dimensions
for count_x = 1:no_of_particles
%x = particle_position(count_x,1);
%y = particle_position(count_x,2);
%z = particle_position(count_x,3);
%soln = x^2 - 3*y*x + z;

%x = particle_position(count_x); 
%soln = x^2-2*x+1;

x = particle_position(count_x);
soln = x-7;

if soln~=0 
current_fitness(count_x) = 1/abs(soln);
else
current_fitness =1000;
end
end

%decide on p_best etc for each particle
for count_x = 1:no_of_particles
if current_fitness(count_x) > p_best_fitness(count_x)
p_best_fitness(count_x) = current_fitness(count_x);
for count_y = 1:dimensions
p_best(count_x,count_y) = particle_position(count_x,count_y);
end
end
end

%decide on the global best among all the particles
[g_best_val,g_best_index] = max(current_fitness);

%g_best contains the position of teh global best
for count_y = 1:dimensions
g_best(count_y) = particle_position(g_best_index,count_y); 
end

%update the position and velocity compponents
for count_x = 1:no_of_particles
for count_y = 1:dimensions
p_current(count_y) = particle_position(count_x,count_y);
end

for count_y = 1:dimensions
particle_velocity(count_y) = particle_velocity(count_y) + c1*rand*(p_best(count_y)-p_current(count_y)) + c2*rand*(g_best(count_y)-p_current(count_y));
particle_positon(count_x,count_y) = p_current(count_y) +particle_velocity(count_y);
end
end


end

g_best 
current_fitness(g_best_index)
clear all, clc % pso example
iter = 1000; % number of algorithm iterations
np = 2; % number of model parameters
ns = 10; % number of sets of model parameters
Wmax = 0.9; % maximum inertial weight
Wmin = 0.4; % minimum inertial weight
c1 = 2.0; % parameter in PSO methodology
c2 = 2.0; % parameter in PSO methodology
Pmax = [10 10]; % maximum model parameter value
Pmin = [-10 -10]; % minimum model parameter value
Vmax = [1 1]; % maximum change in model parameter
Vmin = [-1 -1]; % minimum change in model parameter
modelparameters(1:np,1:ns) = 0; % set all model parameter estimates for all model parameter sets to zero
modelparameterchanges(1:np,1:ns) = 0; % set all change in model parameter estimates for all model parameter sets to zero
bestmodelparameters(1:np,1:ns) = 0; % set best model parameter estimates for all model parameter sets to zero
setbestcostfunction(1:ns) = 1e6; % set best cost function of each model parameter set to a large number
globalbestparameters(1:np) = 0; % set best model parameter values for all model parameter sets to zero
bestparameters = globalbestparameters'; % best model parameter values for all model parameter sets (to plot)
globalbestcostfunction = 1e6; % set best cost function for all model parameter sets to a large number
i = 0; % indicates ith algorithm iteration
j = 0; % indicates jth set of model parameters
k = 0; % indicates kth model parameter
for k = 1:np % initialization
for j = 1:ns
modelparameters(k,j) = (Pmax(k)-Pmin(k))*rand(1) + Pmin(k); % randomly distribute model parameters
    modelparameterchanges(k,j) = (Vmax(k)-Vmin(k))*rand(1) + Vmin(k); % randomly distribute change in model parameters
end
end
for i = 2:iter
for j = 1:ns
x = modelparameters(:,j);
% calculate cost function
costfunction = 105*(x(2)-x(1)^2)^2 + (1-x(1))^2;
    if costfunction < setbestcostfunction(j) % best cost function for jth set of model parameters
bestmodelparameters(:,j) = modelparameters(:,j);
      setbestcostfunction(j) = costfunction;
end
    if costfunction < globalbestcostfunction % best cost function for all sets of model parameters
globalbestparameters = modelparameters(:,j);
bestparameters(:,i) = globalbestparameters;
     globalbestcostfunction(i) = costfunction;
else
bestparameters(:,i) = bestparameters(:,i-1);
globalbestcostfunction(i) = globalbestcostfunction(i-1);
end
end
  W = Wmax - i*(Wmax-Wmin)/iter; % compute inertial weight
for j = 1:ns % update change in model parameters and model parameters
    for k = 1:np
      modelparameterchanges(k,j) = W*modelparameterchanges(k,j) + c1*rand(1)*(bestmodelparameters(k,j)-modelparameters(k,j))...
         + c2*rand(1)*(globalbestparameters(k) - modelparameters(k,j));
      if modelparameterchanges(k,j) < -Vmax(k), modelparameters(k,j) = modelparameters(k,j) - Vmax(k); end
if modelparameterchanges(k,j) > Vmax(k), modelparameters(k,j) = modelparameters(k,j) + Vmax(k); end
if modelparameterchanges(k,j) > -Vmax(k) & modelparameterchanges(k,j) < Vmax(k), modelparameters(k,j) = modelparameters(k,j) + modelparameterchanges(k,j); end
      if modelparameters(k,j) < Pmin(k), modelparameters(k,j) = Pmin(k); end
      if modelparameters(k,j) > Pmax(k), modelparameters(k,j) = Pmax(k); end
end
end
i
end
bp = bestparameters; index = linspace(1,iter,iter);
figure; semilogy(globalbestcostfunction,'k');
set(gca,'FontName','Arial','Fontsize',14); axis tight;
xlabel('iteration'); ylabel('cost function');
figure; q = plot(index,bp(1,:),'g-',index,bp(2,:),'k:');
set(gca,'FontName','Arial','Fontsize',14); axis tight;
legend(q,'x_1','x_2'); xlabel('iteration'); ylabel('parameter')