tsp_pso1.m

用粒子群方法解决TSP（旅行商）问题

function TSP_PSO(ptCorr,nPoint)

% start timing...
tic;
dist = zeros(nPoint);

% calculate the distance matrix...
for i=1:nPoint
    for j=1:nPoint
        dist(i,j) = norm(ptCorr(i,:)-ptCorr(j,:));
    end % for j...
end % for i...

% number of particles...
nParticle = 300; % 100...
curSwarm = zeros(nParticle,nPoint);
curIndiFitBest = zeros(1,nParticle);
curFitBest = inf;
curFitBestPos = 0;
% number of maximum iterations...
maxIteration = 1000;
curIteration = 1;
% assume that the length of velocity is equal to nPoint...
maxVelocLen = round(nPoint*0.3);
% initialization of the velocity length of every particle...
velocLength = round(rand(1,nParticle)*maxVelocLen+0.5);
% velocity of a particle is defined as a sequence of permutations as long as nPoint...
curVelocity = zeros(nParticle,maxVelocLen,2);



% the boundary of swarm fitness...
minFitness = inf;
maxFitness = 0;
terminateError = 1e-3;


% initial and ending inertia weight...
w_init = 0.9;
w_ending = 0.4;
w = 0.5;
% two constant coefficient of cognition and social part...
c1 = 0.3;
c2 = 0.7;


% initaialization of the swarm, on initial position and initial velocity...
for i=1:nParticle
    curSwarm(i,:) = randperm(nPoint);
    curIndiFitBest(i) = TSP_fitness(curSwarm(i,:),dist);

    % length of velocity, that is, length of the permutation sequence...
    for j=1:velocLength(i)
        curVelocity(i,j,:) = round(rand(1,2)*nPoint+0.5);
    end % for...
end % for...

% record matrix of individual that every particle has ever seen...
% the first time seen, the best that it's ever seen...
indiBestEver = curSwarm;

% the best fitness of the swarm...
[curFitBest,curFitBestPos] = min(curIndiFitBest);
swarmBestEver = indiBestEver(curFitBestPos,:);
minFitness = min(curIndiFitBest);
maxFitness = max(curIndiFitBest);

curPath = curSwarm(curFitBestPos,:);
curOutput = [curPath curPath(1)];

% initial display...
plot(ptCorr(curOutput,1),ptCorr(curOutput,2),'b*');
title('Travelling Salesman Problem Solution Using Particle Swarm Optimization');
% pause;
hold on;
figureHandle = plot(ptCorr(curOutput,1),ptCorr(curOutput,2),'r');
hold off;
figure;
% dynamic vector to store all the energy state...
engRecord = curFitBest;
energyHandle = plot(1:size(engRecord,2),engRecord);
title('Fitness of Travelling Salesman Problem Solution');
xlabel('Iterations');
ylabel('Fitness Value');

pMutation = 0.0;

    
while(curIteration<maxIteration && (maxFitness-minFitness)>terminateError*maxFitness)
    % update the position and velocity matrix of the swarm...
    w = (w_init-w_ending)*(maxIteration-curIteration)/maxIteration+w_ending;
%     if (w>1e-3)
%         w = 0.5;
%     end
    c1 = rand;
    c2 = rand;
    for i=1:nParticle
        % v_t+1 = w_k*v_t(+)c_1*(p_i,t-p_t)(+)c_2*(p_g,t-p_t)...
        % x_t+1 = x_t + w_k*v_t+1...
        
        % current fitness of every particle...        
        % recognition part...
        veloc_reg = TSP_posSubtract(curSwarm(i,:),indiBestEver(i,:));
        % veloc_reg = TSP_posSubtract(indiBestEver(i,:),curSwarm(i,:));
        % c1*(p_i,t-p_t)...
        veloc_reg = veloc_reg(1:round(mod(c1,1)*size(veloc_reg,1)),:);
        % length of vector to add to the tail...
        l = round(velocLength(i)*w)+round(mod(c1,1)*size(veloc_reg,1));
        if (l>maxVelocLen)
            l = maxVelocLen;
        end % if...
        % update the velocity vector...
        if (veloc_reg ~= 0)
%             veloc_reg
%             round(velocLength(i)*w):l
%             1:l-round(velocLength(i)*w)+1
%             l
            curVelocity(i,round(velocLength(i)*w)+1:l,:) = veloc_reg(1:l-round(velocLength(i)*w),:);
        end % if...
        velocLength(i) = l;
        
        % social part...
        veloc_soc = TSP_posSubtract(curSwarm(i,:),swarmBestEver);
        % veloc_soc = TSP_posSubtract(swarmBestEver,curSwarm(i,:));
        % c2*(p_g,t-p_t)...
        veloc_soc = veloc_soc(1:round(mod(c2,1)*size(veloc_soc,1)),:);
        % length of vector to add to the tail...
        l = round(velocLength(i)*w)+round(mod(c2,1)*size(veloc_soc,1));
        if (l>maxVelocLen)
            l = maxVelocLen;
        end % if...
        % update the velocity vector...
        if (veloc_soc ~= 0)
%             veloc_soc
%             round(velocLength(i)*w):l
%             1:l-round(velocLength(i)*w)+1
%             l
            curVelocity(i,round(velocLength(i)*w)+1:l,:) = veloc_soc(1:l-round(velocLength(i)*w),:);
        end % if...
        velocLength(i) = l;
        
        % update the velocity and position vector...
        t = curVelocity(i,:,:);
%         t
        % t(i,round(w*velocLength(i))+1:maxVelocLen,:) = zeros(maxVelocLen-round(w*velocLength(i)),2);
                
        % cause the variable t is truncated from curVelocity and its length
        % is only 1...
        
        if (round(w*velocLength(i)) > maxVelocLen)
            ttt = maxVelocLen;
        else
            ttt = round(w*velocLength(i));
        end % if...
        for j=1:ttt
            if (t(1,j,1) ~= 0)
                tt = curSwarm(i,t(1,j,1));
                curSwarm(i,t(1,j,1)) = curSwarm(i,t(1,j,2));
                curSwarm(i,t(1,j,2)) = tt;
            end % if...
        end % for...
        
    end % for...
    
%     % additional process of mutation...
%     pMutation = curIteration/maxIteration*0.3;
%     for i=1:nParticle
%         if (rand < pMutation)
%             mutationor = round(rand(1,2)*nPoint+0.5);
%             curSwarm(i,min(mutationor):max(mutationor)) = fliplr(curSwarm(i,min(mutationor):max(mutationor)));
%         end
%     end % for...
    
    % update the best record ever seen...
    for i=1:nParticle
        t = TSP_fitness(curSwarm(i,:),dist);
        if (t < curIndiFitBest(i))
            curIndiFitBest(i) = t;
            indiBestEver(i,:) = curSwarm(i,:);
        end % if...
    end % for...
    for i=1:nParticle
        if (curIndiFitBest(i) < curFitBest)
            curFitBest = curIndiFitBest(i);
            curFitBestPos = i;
        end % if...
    end % for...
%     curFitBestPos
%     indiBestEver
    swarmBestEver = indiBestEver(curFitBestPos,:);
    minFitness = min(curIndiFitBest);
    maxFitness = max(curIndiFitBest);
    
    % displaying the output on the screen in every iteration...
    curOutput = [curSwarm(curFitBestPos,:) curSwarm(curFitBestPos,1)];
    
    % add the most up-to-date energy state to record...
    engRecord = [engRecord curFitBest];
    set(figureHandle,'xdata',ptCorr(curOutput(:),1),'ydata',ptCorr(curOutput(:),2));
    set(energyHandle,'xdata',1:size(engRecord,2),'ydata',engRecord);
    drawnow;
    
    % output optimization information in this iteration...
    fprintf('[%s]Iteration: %d\n',datestr(now),curIteration);
    fprintf('[%s]Current path :%s\n',datestr(now),sprintf('%d ',curSwarm(curFitBestPos,:)));
    fprintf('[%s]Best fitness all particle''s ever had :%s\n',datestr(now),sprintf('%d ',swarmBestEver));
    fprintf('[%s]Current fitness[min max] :%f[%f,%f]\n',datestr(now),curFitBest,minFitness,maxFitness);
    fprintf('[%s]Time elapsed: %f s\n',datestr(now),toc);
    
    curIteration = curIteration + 1;
end % while...

% final marking the citis and orders...
figure(1);
for i=1:nPoint
    text(ptCorr(curOutput(i),1),ptCorr(curOutput(i),2),['\cdot',int2str(i)],'fontsize',8);
end % i...