【TSP问题】基于混合粒子群算法求解TSP问题-算法实现

• 个体编码
  粒子个体编码采用整数编码的方式，每个粒子表示历经的所有城市，比如当历经的城市数为10，个体编码为[9 4 2 1 3 7 6 10 8 5]，表示城市遍历从9开始，经过4、2、1、3、⋯ \dotsm⋯最终返回起点城市9，从而完成TSP遍历。

• 适应度值
  粒子适应度值表示为遍历路径的长度，计算公式为f i t n e s s ( i ) = ∑ i , j = 1 n p a t h i , j (1) fitness(i)=\sum_{i,j=1}^{n}path_{i,j}\tag{1}fitness(i)=i,j=1∑n​pathi,j​(1)其中，n nn为城市数量；p a t h i , j path_{i,j}pathi,j​为城市i , j i,ji,j间路径长度。

• 交叉操作
  个体通过个体极值和群体极值交叉来更新，交叉方法采用整数交叉法。首先选择两个交叉位置，然后把个体和个体极值或个体与群体极值进行交叉，假定随机选取的交叉位置为3和5，操作方法如下：
  个体[9 4 2 1 3 7 6 10 8 5] 极值[9 2 1 6 3 7 4 10 8 5]
  交叉为新个体[9 4 1 6 3 7 6 10 8 5]
  产生的新个体如果存在重复位置则进行调整，调整方法为用个体中未包括的城市代替重复包括的城市，如下所示：
  [9 4 1 6 3 7 6 10 8 5]调整为[9 4 2 1 3 7 6 10 8 5]
  对得到的新个体采用了保留优秀个体策略，只有当新粒子适应度值好于旧粒子时才更新粒子。

• 变异操作
  变异方法采用个体内部两位互换方法，首先随机选择变异位置p o s 1 pos1pos1和p o s 2 pos2pos2，然后把两个变异位置互换，假设选择的变异位置为2和4，变异操作如下所示：
  [9 4 2 1 3 7 6 10 8 5] 变异为[9 1 2 4 3 7 6 10 8 5]
  对得到的新个体采用了保留优秀个体策略，只有当新粒子适应度值好于旧粒子时才更新粒子。

三、MATLAB程序实现

根据混合粒子群算法原理，在MATLAB中编程实现基于混合粒子群的TSP搜索算法。

• 适应度函数
  适应度函数计算个体适应度值，个体适应度值为路线总长度，代码如下：

function indiFit = fitness(x, cityCoor, cityDist)
%% 该函数用于计算个体适应度值
% x           input      个体
% cityCoor    input      城市坐标
% cityDist    input      城市距离
% indiFit     output     个体适应度值

m = size(x, 1);
n = size(cityCoor, 1);
indiFit = zeros(m, 1);
for i = 1:m
    for j = 1:n-1
        indiFit(i) = indiFit(i)+cityDist(x(i, j), x(i, j+1));
    end
    indiFit(i) = indiFit(i)+cityDist(x(i, 1), x(i, n));
end

• 粒子初始化
  粒子初始化用于初始化粒子，计算粒子适应度值，并根据适应度值确定个体最优粒子和群体最优粒子。程序代码如下：

nMax = 200;                      % 进化次数
indiNumber = 1000;               % 个体数目
individual = zeros(indiNumber,n);
% 初始化粒子位置
for i = 1:indiNumber
    individual(i, :) = randperm(n);    
end
%% 计算种群适应度
indiFit = fitness(individual, cityCoor, cityDist);
[value, index] = min(indiFit);
tourPbest = individual;                          % 当前个体最优
tourGbest = individual(index, :) ;               % 当前全局最优
recordPbest = inf*ones(1, indiNumber);           % 个体最优记录
recordGbest = indiFit(index); % value            % 群体最优记录
xnew1 = individual;

• 交叉操作
  交叉操作把粒子同个体极值和群体极值进行交叉，从而得到较好的个体，交叉操作代码如下：

%% 交叉操作
function a = Recombin(a, b)
%% 交叉
% 输入：
% a为待交叉的个体，b为交叉的个体极值
% 输出：
% a为交叉后得到的新个体
n = length(a);
c1 = unidrnd(n-1);         % 产生交叉位
c2 = unidrnd(n-1);         % 产生交叉位
while c1 == c2
    c1 = round(rand*(n-2))+1;
    c2 = round(rand*(n-2))+1;
end
chb1 = min(c1, c2);
chb2 = max(c1, c2);
cros = b(chb1:chb2);        % 交叉区域元素
ncros = length(cros);      % 交叉元素个数
% 删除与交叉区域相同元素
for j=1:ncros
    for k=1:n
        if a(k) == cros(j)
            a(k) = 0;
            for t = 1:n-k
                temp = a(k+t-1);
                a(k+t-1) = a(k+t);
                a(k+t)=temp;
            end
        end
    end
end
%插入交叉区域
a(n-ncros+1:n) = cros;        
% a0 = a; b0 = b;
% for i = chb1:chb2
%     a1 = a; b1 = b;
%     a(i) = b0(i);
%     b(i) = a0(i);
%     x = find(a == a(i));
%     y = find(b == b(i));
%     i1 = x(x ~= i);
%     i2 = y(y ~= i);
%     if ~isempty(i1)
%         a(i1) = a1(i);
%     end
%     if ~isempty(i2)
%         b(i2) = b1(i);
%     end
% end

• 变异操作
  变异操作对自身进行变异，从而得到更好的个体。变异操作代码如下：

function x = Mutate(x)
% 输入：x        个体
% 输出：x        变异后的个体
n = length(x);
while true
    c1 = round(rand*(n-1))+1;       % 产生变异位
    c2 = round(rand*(n-1))+1;       % 产生变异位
    if c1 ~= c2
        break;
    end
end
temp = x(c1);
x(c1) = x(c2);
x(c2) = temp;

• 仿真结果
  采用混合粒子群算法规划TSP路径，各城市的初始位置如图2所示。
  

图2 城市初始位置

混合粒子群算法的进化次数为100，种群规模为100，算法进化过程中最优粒子适应度值变化和规划出的最优路径如图3和图4所示。

图3 适应度值变化

图4 规划出的最优路径

四、延伸阅读

使用混合粒子群算法规划其他城市模型的最短路径。
城市分布图如图5所示，混合粒子群算法的进化次数为200，种群规模为1000。采用混合粒子群算法规划的路径如图6所示。

图5 城市分布图

图6 规划路径

 

从图2~图6可以看出，基于混合粒子群算法的TSP算法可以解决规模较小的旅行商问题，对于规模较大的旅行商问题，混合粒子群算法也可以得到较优路径。

代码下载或者仿真咨询添加QQ1575304183

五、参考文献

[1] Kennedy J . Particle swarm optimization[J]. Proc. of 1995 IEEE Int. Conf. Neural Networks, (Perth, Australia), Nov. 27-Dec. 2011, 4(8):1942-1948 vol.4.
[2] 王伟. 混合粒子群算法及其优化效率评价[J]. 中国水运:学术版, 2007, 7(006):102-103.
[3] 郁磊等. MATLAB智能算法30个案例分析(第2版)[M].北京航空航天大学出版社.2015年.