main.m

基本的粒子群优化算法程序

function main()
% 主函数
% Initialization phase
% Input City
MAXITER=2500;
NumCities=30;
NumAnts=10;
RHO=0.1;
beita=2;
q0=0.8;
MinLenTour=999999;
NXY=Coordinates(30);
X(1:NumCities)=NXY(:,2)';
Y(1:NumCities)=NXY(:,3)';

BestTour=zeros(1,NumCities+1);

for i=1:NumCities
    for j=1:NumCities
		Dist(i,j)=(X(j)-X(i))^2+(Y(j)-Y(i))^2;
    end;
end;
Dist=sqrt(Dist);

rand('state',sum(100*clock));
Lnn=LNN(Dist,NumCities);
Tao0 = 1.0/(Lnn*NumCities);
Tao=ones(NumCities)*Tao0;
% **************** main loop ******************
for iter=1:MAXITER
	% place each ant in a randomly chosen city 
	% no cities visited yet  
	VisitedCities=zeros(NumAnts,NumCities);
	Tour=zeros(NumAnts,NumCities+1); 
	% place only one ant per city */
    StartCity=round(1+(NumCities-1)*rand(1,10));
    for k=1:NumAnts
        VisitedCities(k,StartCity(k))=1;
	    Tour(k,1)=StartCity(k);
    end;
	% The ants build their tours
	for t=2:NumCities+1
		% Execute one step for each ant k
		for k=1:NumAnts
			% choose next city
			if VisitedCities(k,:)==ones(1,NumCities)
			    Tour(k,t)=StartCity(k);  % ants go back to the initial city
            else 
				Tour(k,t)=ChooseNextCity(NumCities,k,Tour(k,t-1),VisitedCities,Tao,Dist,q0,beita);
				VisitedCities(k,Tour(k,t))=1;
            end;
        end;
        % update pheromone level using a local rule 
        for k=1:NumAnts
            Tao(Tour(k,t-1),Tour(k,t))=(1-RHO)*Tao(Tour(k,t-1),Tour(k,t))+ RHO*Tao0; % this is the simple ACS
        end;
    end;
	% Compute the length of the Tour found by each ant 
	SumLen=zeros(1,NumAnts);
	for k=1:NumAnts
        for t=1:NumCities
            SumLen(k)=SumLen(k)+Dist(Tour(k,t),Tour(k,(t+1)));
        end;
    end;
    % Find the min tour
    [MinLenTour,BestTour]=FindNewBestSolution(NumAnts,NumCities,MinLenTour,Tour,SumLen,BestTour);
	% update pheromone level using a global rule */
	BestTour
    MinLenTour
    for j=1:NumCities
        Tao(BestTour(j),BestTour(j+1))=(1-RHO)*Tao(BestTour(j),BestTour(j+1)) + RHO*(1.0/MinLenTour);
    end;
iter
end;


function NextCity=ChooseNextCity(NumCitiesCN,kCN,antposCN,VisitedCitiesCN,TaoCN,DistCN,q,beitaCN)
% 蚂蚁选择下一个未访问过的结点
nonvisited=0;
ArgMaxWeightedPhero = 1;
PartSum=0.0;
WeightedPheromone=zeros(1,NumCitiesCN);
Pij=zeros(1,NumCitiesCN);
PijSum=zeros(1,NumCitiesCN+1);
for i=1:NumCitiesCN
	if ~VisitedCitiesCN(kCN,i)
		nonvisited=nonvisited+1;
		WeightedPheromone(i)=TaoCN(antposCN,i)/(DistCN(antposCN,i)^beitaCN);
		% SumWeightedPheromone += WeightedPheromone[antpos][i];
		if WeightedPheromone(i)>WeightedPheromone(ArgMaxWeightedPhero)
			ArgMaxWeightedPhero = i; % find the max ArgweightedPhero
        end;
    else 
        WeightedPheromone(i)=0;
    end;
end;
% exploration-exploitation */
prob_explore=rand(1);
if prob_explore<q % choose the max ArgWP as the next city
    NextCity = ArgMaxWeightedPhero;
else              % choose next city by prob
    if nonvisited==1
        RndCity=1;
    else
        for i=1:NumCitiesCN
            if ~VisitedCitiesCN(kCN,i)
                PartSum=PartSum+TaoCN(antposCN,i)/(DistCN(antposCN,i)^beitaCN);
            end;
        end;
        j=0;
        for i=1:NumCitiesCN
            if ~VisitedCitiesCN(kCN,i)
                j=j+1;
                Pij(j)=(TaoCN(antposCN,i)/(DistCN(antposCN,i)^beitaCN))/PartSum;
                PijSum(j+1)=PijSum(j)+Pij(j);
            end;
        end;
        rnd=rand(1);
        RndCity=0;
        for i=1:j
            RndCity=RndCity+1;
            if rnd<=PijSum(i)
                break;
            end;
        end;
    end;
    i=1;
    j=1;
    while (i<=NumCitiesCN)&&(j<=RndCity)
        if ~VisitedCitiesCN(kCN,i)
            j=j+1;
        end;
        i=i+1;
    end;
    NextCity=i-1;
end;

function City=Coordinates(Num)
% 节点的坐标
if Num==30
    City=[1  41.0 94.0
2  37.0 84.0
3  54.0 67.0
4  25.0 62.0
5  7.0  64.0
6  2.0  99.0
7  68.0 58.0
8  71.0 44.0
9  54.0 62.0
10 83.0 69.0
11 64.0 60.0
12 18.0 54.0
13 22.0 60.0
14 83.0 46.0
15 91.0 38.0
16 25.0 38.0
17 24.0 42.0
18 58.0 69.0
19 71.0 71.0
20 74.0 78.0
21 87.0 76.0
22 18.0 40.0
23 13.0 40.0
24 82.0 7.0
25 62.0 32.0
26 58.0 35.0
27 45.0 21.0
28 41.0 26.0
29 44.0 35.0
30 4.0  50.0];
else
    error('The number of city is wrong!')
end;


function [NewMinLenTour,Best]=FindNewBestSolution(NumAntsBS,NumCitiesBS,MinLenTourBS,TourBS,SumLenBS,BestTourBS)
% 找出当前最佳已知解
NewMinLenTour=MinLenTourBS;
Best=BestTourBS;
for k=1:NumAntsBS
    if SumLenBS(k)<NewMinLenTour
        Best=TourBS(k,:);
        NewMinLenTour=SumLenBS(k);
    end;
end;


function Length=LNN(Dist,NumCities)
% the nearest neighbor heuristic
Length=0.0;
Tabu=ones(1,NumCities);
StartCity=round(1+(NumCities-1)*rand(1,1));
Start=StartCity;
for i=1:NumCities-1
    Tabu(Start)=0;
    Unvisited=(find(Tabu>0));
    List=zeros(1,NumCities-i);
    for j=1:NumCities-i
        List(j)=Dist(Start,Unvisited(j));
    end;
    [minD,P]=min(List);
    Length=Length+minD;
    Start=Unvisited(P);
    Start
end;
Length=Length+Dist(Start,StartCity)