function  PSOforLayout() 

% A SIMPLE IMPLEMENTATION OF THE 
% Particle Swarm Optimization IN MATLAB 
clc
clear
close all


% Initializing variables 
success = 0;                    % Success flag 
PopSize = 40;                   % Size of the swarm 
MaxIt = 5000;                   % Maximum number of iterations 
iter = 0;                       % Iterations’counter 
c1 = 2;                         % PSO parameter C1 
c2 = 2;                         % PSO parameter C2 
w = 0.6;                        % inertia weight 
omega=1;

% Problem Definition

f = 'Layout'; 

mark=input('布局優化問題 (1,2,3),(待布物為40,5,7 項)       ' );



switch mark
    case 1           
        n=40;
       dim = 2*n;                         % Dimension of the problem 
       upbnd = 880;                      % Upper bound for init. of the swarm 
       lwbnd = -880;                     % Lower bound for init. of the swarm 
       R=880;
       
       RM=[106 11;112 12;98 9;105 11;93 8;103 10;82 6;93 8;117 13;81 6;...
       89 7;92 8;109 11;104 10;115 13;110 12;114 12;89 7;82 6;120 14;...
        108 11;86 7;93 8;100 10;102 10;106 11;111 12;107 11;109 11;91 8;...
        111 12;91 8;101 10;91 8;108 11;114 12;118 13;85 7;87 7; 98 9];

    
  case 2               
   n=5;
   dim = 2*n;                         % Dimension of the problem  
   upbnd = 130;                       % Upper bound for init. of the swarm 
   lwbnd = -130;                     % Lower bound for init. of the swarm 
   R=130;

   RM=[20.71 20.71;50 50;50 50;50 50;50 50];

    
case 3            
                   
   n=7;
   dim = 2*n;                         % Dimension of the problem  
   upbnd =50 ;                      % Upper bound for init. of the swarm 
    lwbnd =-50 ;                     % Lower bound for init. of the swarm 
    R=50;

    RM=[10 100;11 121;12 144;11.5 132.25;9.5 90.25;  8.5 72.25;10.5 110.25];
end


r=RM(:,1);
m=RM(:,2);



averfit=[];
bestfit=[];
bestx=[];


% Initializing swarm and velocities 
popul = rand(dim, PopSize)*(upbnd-lwbnd) + lwbnd; 
vel = rand(dim, PopSize); 

for i = 1:PopSize, 
   fpopul(i) = feval(f, popul(:,i),m,r,R); 
end 

bestpos = popul; 
fbestpos = fpopul;   %fbestpos is the fitness value for every particle  

% Finding best particle in initial population 
[fbestpart,g] = min(fpopul); 

averfit=[averfit;mean(fpopul)];
bestfit=[bestfit;fbestpart];
bestx=popul(:,g);

while iter < MaxIt,     
   iter = iter + 1; 

   % VELOCITY UPDATE 
   for i=1:PopSize, 
       A(:,i) = bestpos(:,g); 
   end 
   R1 = rand(dim, PopSize); 
   R2 = rand(dim, PopSize); 
   vel = w*vel + c1*R1.*(bestpos-popul) + c2*R2.*(A-popul); 

   % SWARMUPDATE 
   popul = popul + vel;
   
   % Evaluate the new swarm 
   
   for i = 1:PopSize, 
       fpopul(i) =  feval(f, popul(:,i),m,r,R); 
   end 
   
  
   % Updating the best position for each particle 
   changeColumns = fpopul < fbestpos; 
   fbestpos = fbestpos.*( 1-changeColumns) + fpopul.*changeColumns; 
   bestpos(:, find(changeColumns)) = popul(:, find(changeColumns)); 
  
   % Updating index g 
   [fbestpart, g] = min(fbestpos); 
   %currentTime = etime(clock,startTime); 
   % Checking stopping criterion 
   averfit=[averfit;mean(fpopul)];
   bestfit=[bestfit;fbestpart];
   bestx=popul(:,g);

end 


% Output arguments 
disp('The best position is')
bestx
disp('The best fitness is') 
fbestpart

x=bestx;
x1=[];
y1=[];
for i=1:n
x1=[x1;x(2*i-1)];
y1=[y1;x(2*i)];
end

disp('The minimum Rmin is')
Rmin=max(sqrt(x1.*x1+y1.*y1)+r)
   
disp('The minium nonbalance quantity Qmin is')
Qmin=sqrt((sum((omega^2*m).*x1))^2+(sum((omega^2*m).*y1))^2)
   

figure
plot([1:(iter+1)]',averfit,'r-',[1:(iter+1)]',bestfit,'k-');

xlabel('iterations');ylabel('fitness');
legend('average fitness','best fitness');

figure
plot([1:(iter+1)]',log(averfit),'r-',[1:(iter+1)]',log(bestfit),'k-');

xlabel('iterations');ylabel('log(fitness)');
legend('average fitness','best fitness');


figure
x=-R:0.1:R;
plot(x,sqrt(R^2-x.^2),'k',x,-sqrt(R^2-x.^2),'k');
hold on
for i=1:n
    x=(x1(i)-r(i)):0.1:(x1(i)+r(i));
    
    y11=y1(i)-sqrt(r(i)^2-(x-x1(i)).^2);
    y22=y1(i)+sqrt(r(i)^2-(x-x1(i)).^2);  
    plot(x,y11,'r',x,y22,'r');
end

x11=-Rmin:0.1:Rmin;
plot(x11,sqrt(Rmin^2-x11.^2),'b:',x11,-sqrt(Rmin^2-x11.^2),'b:');