【原】總（tu）結（cao）粒子群算法(PSO)解決旅行商問題(TSP)

粒子群算法(PSO)是一套比較經典的算法, 旅行商問題(TSP)同樣是一個經典的問題。如果想用PSO去解決TSP問題的話，那麼應該如何去解決呢？

于是就有查資料，找到PSO解決TSP問題論文一文。

初看之下一陣欣喜，因為我發現，如果按照論文中的方法能夠成功的話，那麼包括布谷鳥，螢火蟲都可以通過類似的辦法進行有意義的嘗試。

按照論文的思路，我撰寫了如下代碼

% main.m
% 調用
clc
clear
close all

number = ; % 14個城市
global cities
cities = rand(number,) *  - ; % 假設城市分布在 [-100,-100] ,[100,100]的平面内

[best,par,ItorMean,ItorBest] = PSO4TSP_v1(cities);

netplot(cities,best.loc);
figure
itor = length(ItorMean);
plot(:itor,ItorMean,'r',:itor,ItorBest,'g');
           

% PSO 算法，globalBest為全局最優，par是最終的粒子資訊，ItorMean是每次疊代平均适應度，ItorBest是每次疊代最佳适應度
function [globalBest,par,ItorMean,ItorBest] = PSO4TSP_v1(cities)
    %% 參數設定
    number_cities = size(cities,); % 城市數
    number_pars = ; % 粒子個數
    maxItor = ; % 最大疊代次數
    initSSnumber = number_cities;
    ss = zeros(initSSnumber,);
    globalBest.fitness = -inf;
    ItorMean = zeros(,maxItor);
    ItorBest = zeros(,maxItor);
    %% 粒子初始化
    for i =  : number_pars     
        % 初始化粒子位置
        par(i).loc = randperm(number_cities); 
        % 初始化交換
        for j =  : initSSnumber
            ss(j,:) = randperm(number_cities,);
        end
        par(i).v = ss;
        par(i).best = par(i);
        par(i).fitness = fitness(par(i).loc);
        if globalBest.fitness < par(i).fitness
            globalBest = par(i);
        end
    end
    %% 疊代求優
    eachItor = zeros(,number_pars);
    for i =  : maxItor        
        for j =  : number_pars
            t = rand;
            u = rand;
            ssLoc = getSS(par(j).best.loc,par(j).loc);
            ssLocLen = size(ssLoc,);
            ssGlo = getSS(globalBest.loc,par(j).loc);
            ssGloLen = size(ssGlo,);
            temp1 = rand(,ssLocLen);
            Item2 = ssLoc(temp1<t,:);
            temp2 = rand(,ssGloLen);
            Item3 = ssGlo(temp2<u,:);
            par(j).loc = doSS(par(j).loc,par(j).v);
            ss = unionSS(unionSS(par(j).v,Item2),Item3);
            par(j).v = ss;
            myfitness = fitness(par(j).loc);
            eachItor(j) = myfitness;
            if myfitness > par(j).fitness
                par(j).fitness = myfitness;
                par(j).best = par(j);
                if myfitness > globalBest.fitness
                    globalBest = par(j);
                end
            else
                par(j).fitness = myfitness;
            end
        end
        ItorMean(i) = mean(eachItor);
        ItorBest(i) = max(eachItor);
    end
end
           

% netplot.m
function netplot(city,n)        %連線各城市，将路線畫出來
    figure
    hold on
    plot(city(:,),city(:,),'*');
    line([city(:,);city(,)],[city(:,);city(,)]);
end
           

% fitness.m
% 計算城市間的距離
function z = fitness(n)
    global cities
    tstart = n;
    tend = [n(:end),n()];
    z = sum(distance(cities(tstart,),cities(tstart,),cities(tend,),cities(tend,)));
    z=/z;
end
           

% getSS.m
function ss = getSS(v1,v2)
% 求SS算子
% v1，v2是一組向量，v1，v2所包含的元素應該相同。
% 若v1,v2順序相同，則傳回ss = [];
% 若順序不同，則傳回一個ss算子，即v2通過ss變換可以得到v1
ss = [];
while 
    idx = find(v1 ~= v2);
    if isempty(idx)
        break;
    end
    idx2 = find(v2 == v1(idx()));
    so = [idx(),idx2];
    ss = [ss;so];
    v2 = doSO(v2,so);
end
           

% doSS.m
function v = doSS(v,ss)
% SS操作函數
% v:之前的方案
% ss:SS算子，用作交換(n*2)的矩陣
% v:SS操作之後的結果
if ~isempty(ss)
    [m,n] = size(ss);
    if m ==  && n ~= 
        ss = ss';
        m = n;
        n = ;
    end
    if n ~= 
        help doSS
        error('ss must be n*2 matrix');
    end

    for i =  : m
        v = doSO(v,ss(i,:));
    end
end
end
           

% doSO.m
function v = doSO(v,so)
% SO操作函數
% v:之前的方案
% so:SO算子，用作交換
% v:SO操作之後的結果
if ~isempty(so)
    if length(so) ~= 
        help doSO
        error('length of "so" iso must equals ');
    else
        len = length(v);
        if so() > len || so() > len || so() <  || so() < 
            help doSO
            error('"so" iso is not the index of v');
        else
            v(so()) = v(so()) + v(so());
            v(so()) = v(so()) - v(so());
            v(so()) = v(so()) - v(so());
        end 
    end
end
end
           

% getSS.m
% 論文中的減法操作
function ss = getSS(v1,v2)
% 求SS算子
% v1，v2是一組向量，v1，v2所包含的元素應該相同。
% 若v1,v2順序相同，則傳回ss = [];
% 若順序不同，則傳回一個ss算子，即v2通過ss變換可以得到v1
ss = [];
while 
    idx = find(v1 ~= v2);
    if isempty(idx)
        break;
    end
    idx2 = find(v2 == v1(idx()));
    so = [idx(),idx2];
    ss = [ss;so];
    v2 = doSO(v2,so);
end
           

% unionSS.m
% 論文中的⊕操作
function ss = unionSS(ss1,ss2)
% ss1,ss2算子合并操作

vm1 = max(max(ss1));
vm2 = max(max(ss2));
if isempty(ss1)
    if isempty(ss2)
        ss = [];
    else
        ss = ss2;
    end
elseif isempty(ss2)
    ss = ss1;
else
    vm = max(vm1,vm2);
    v =  : vm;
    v2 = doSS(doSS(v,ss1),ss2);
    ss = getSS(v,v2);
end
end
           

當我信心滿滿的運作的時候，效果卻令人大跌眼鏡。

雖然從疊代的圖中能看出來，算法确實是取尋優了而且也收斂，但是得到的最終效果卻不能令人滿意。

後來，我又查閱了相關資料，發現這兩篇文章java版本的PSO求解TSP問題，C++版本的PSO求解TSP問題參考的同一篇文獻，而且結果同樣有不能令人接受，是以隻能暫時認為這篇古老的文獻的公式出錯。

展望下未來的情況，也許可以找到更好的文獻取代這篇文獻，亦或者把之前的疊代公式改正确。

如果是後者，我們可以通過經典PSO問題類比推理出求解TSP的V’id可能需要在原先的公式加上一個随機速度Vir。推測Vir為一個SO或者是随着疊代次數增加而算子長度向0收斂的SS。筆者正向這個方向進行嘗試。