ps_particleswarm.m

其中提到遺傳學的程式碼與應用提供給次淚相向的研究者參考下載

function [trueCost, globalpar]=PS_particleswarm(obj_fcn,...
               popuSize,gen_no,ranges,maxmin,Type);
%% Particle Swarm Optimization - PSO
%%function [trueCost, globalpar]=PS_particleswarm(obj_fcn,...
%%               popuSize,gen_no,ranges,maxmin,Type);
%%==========================================================
%% VERY IMPORTANT!!!
%% USE SMALL VALUES FOR RANGES ONLY AT THE BEGINNING,
%% ENLARGE RANGES ONCE YOU GOT A GOOD RESULT.
%%==========================================================
%% obj_fun is the string of filename to evaluate the fitness
%% popuSize, gen_no=generation number, ranges=[LB;UB];
%% Set maxmin= 1 if it is a maximization problem,
%% Set maxmin=-1 if it is a minimization problem,
%% Set Type=0(defaut) for real particles, set to 1 if particles
%%     are constrainted to integers
%%
%% For example:
%% [1]. Create and save a function named PS_f24.m as follows:
%function PI=PS_f24(chro)
%% Fitness max. function evaluation
%x=chro(:,1); y=chro(:,2); 
%A1=(sin(sqrt(x.^2+y.^2))).^2-0.5;
%A2=(1+0.001*(x.^2+y.^2)).^2;
%PI=0.5-A1./A2;
%% [2]. Then execute the following:
%obj_fcn = 'PS_f24';	% Objective function
%ranges = [-100 -100
%          100  100];	% Range of the input variables
%gen_no=400;
%popuSize=50;
%maxmin=1;
%PS_particleswarm(obj_fcn,popuSize,gen_no,ranges,maxmin);

% Haupt & Haupt
% 2003
% PenChen Chou. 2005-4-13. 2005-4-16.

fprintf('   [Welcome to Particle Swarm Optimization Algorithm]\n');
if maxmin==1
    fprintf('You want to find a MAXIMUM for this problem<===\n');
else
    fprintf('You want to find a MINIMUM for this problem<===\n');
end
if nargin<=5, Type=0; end    % Default, real particles handled.
rand('state',sum(100*clock));
c1 = 3;         % cognitive parameter
c2 = 4-c1;      % social parameter
C=1;            % constriction factor
[junk, npar]=size(ranges);
% Initializing swarm and velocities
TD =ones(popuSize,1)*diff(ranges); % Total Differences
LOW=ones(popuSize,1)*ranges(1,:);
par=rand(popuSize,npar);  % random population of continuous values
vel=rand(popuSize,npar); % random velocities
par=par.*TD+LOW;         % Match bounds
if Type==1, par=round(par); end;
%vel=vel.*TD+LOW;         % Match bounds
% Evaluate initial population
cost=-maxmin*feval(obj_fcn,par);   % calculates population cost using ff
minc(1)=min(cost);      % min cost
meanc(1)=mean(cost);    % mean cost
globalmin=minc(1);      % initialize global minimum
 
% Finding best particle in initial population
[globalcost,indx] = min(cost);PreCost=globalcost;
globalpar         = par(indx,:);
% Initialize local minimum for each particle
localpar  = ones(popuSize,1)*globalpar;    % location of local minima
localcost = ones(popuSize,1)*globalcost;  % cost of local minima
if gen_no==10 & popuSize==10
   [par cost]
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start iterations
iter = 0;       % counter

while iter < gen_no
    iter = iter + 1;

% update velocity = vel
    w=(gen_no-iter)/gen_no; %inertia weiindxht
    r1 = rand(popuSize,npar);    % random numbers
    r2 = rand(popuSize,npar);    % random numbers
    vel = C*(w*vel + c1 *r1.*(localpar-par) +...
        c2*r2.*(ones(popuSize,1)*globalpar-par));
%vel,pause
% update particle positions
    par = par + vel;    % updates particle position
    par1=par;
    for i=1:npar   % Match bounds
        II=find(par(:,i)>ranges(2,i));
        if ~isempty(II), par(II,i)=ranges(2,i); end
        II=find(par(:,i)<ranges(1,i));
        if ~isempty(II), par(II,i)=ranges(1,i); end
    end
    if Type==1, par=round(par); end;
%[par1 par]
% Evaluate the new swarm
    PI=feval(obj_fcn,par);%[size(PI)],pause
    cost = -maxmin*PI;   % evaluates cost of swarm
    if gen_no==10 & popuSize==10
       [par cost]
    end
%[99 length(cost) length(localcost)],keyboard,pause
% Updating the best local position for each particle
    [m1,n1]=size(cost);[m2,n2]=size(localcost);
    if m1==1, cost=cost'; end
    if m2==1, localcost=localcost'; end
    [m1,n1]=size(cost);[m2,n2]=size(localcost);
    %keyboard,pause
    if m1~=m2, 
       error('ERROR:Size of cost is not matched to that of particles.'); 
       error('ERROR:Check sizes of cost and particles in your program!!!'); 
    end
    
    bettercost = cost < localcost;
    localcost = localcost.*not(bettercost) + cost.*bettercost;
    localpar(find(bettercost),:) = par(find(bettercost),:);
%if ~isempty(find(bettercost))    
    %find(bettercost)'
%end
%[i localcost']

% Updating index t
    PreCost=globalcost;
    [temp, t] = min(localcost);
    fprintf('\nIter. #=%3d;',iter);
    if temp<globalcost
       globalpar=par(t,:); indx=t; globalcost=temp;
       fprintf(' PI=%16.6f; ',abs(temp));%pause(1)
       if length(globalpar)==2
           fprintf('para=[%.6f %.6f];',globalpar);
       elseif length(globalpar)==3
           fprintf('para=[%.6f %.6f %.6f];',globalpar);
       else
           ;
       end
       if abs(globalcost-PreCost)<1e-10,break; end
    end
    %[iter globalpar globalcost]     % print output each iteration
    minc(iter+1)=min(cost);         % min for this iteration
    globalmin(iter+1)=globalcost;   % best min so far
    meanc(iter+1)=mean(cost);       % avg. cost for this iteration  
end% while
trueCost=-maxmin*globalcost;

format long
%[globalpar trueCost]
fprintf('\n===>Final results:\n');
if maxmin==1
   fprintf('MaxPI   = %16.8f;\n',trueCost);
else
   fprintf('MinPI   = %16.8f;\n',trueCost);
end
for i=1:length(globalpar)
    fprintf('para(%d) = %16.8f;\n',i,globalpar(i));
end

format short
figure(24)
iters=0:length(minc)-1;
plot(iters,minc,iters,meanc,'--',iters,globalmin,':');
xlabel('generation');ylabel('PI');
%text(0,minc(1),'best');text(1,minc(2),'population average')
legend('local best','mean','BEST');
%keyboard
%[par cost/1000]