遗传算法（三）解TSP问题

最新推荐文章于 2022-09-22 10:16:55 发布

救世腹肌2298

最新推荐文章于 2022-09-22 10:16:55 发布

阅读量1.1w

点赞数 10

分类专栏：遗传算法文章标签：算法 matlab程序源代码遗传算法 tsp问题

本文链接：https://blog.csdn.net/qq_37043191/article/details/77885704

版权

遗传算法专栏收录该内容

4 篇文章 10 订阅

订阅专栏

巡回旅行商问题TSP

TSP问题：给定一个图 $G = (V, E)$ ，每个边 $e \in E$ 中有非负权值 $\omega (e)$ ，寻找 $G$ 的Hamilton圈 $C$ ，使 $C$ 的总权 $W(C)=\sum_{e \in E(C)}\omega (e)$ 最小。即旅行商的路程最短。

问题抽象：例如一条城市顺序表 $C =\{1, 2, 3, 4, 5,6,7,8,9\}$ ，先将即将探访的城市在顺序表中的顺序记录在染色体 $p$ 中；然后每走一个城市，将其在顺序表C中划去，那么下一个要走的城市在新顺序表中的位置将被记录在 $p$ 中。例如，若实际的行走顺序 $T=\{1,2,3,4,5,6,7,8,9\}$ 则对应的染色体 $p = \{1,1,1,1,1,1,1,1,1\}$ ；若实际的行走顺序为 $T=\{2,1,5,6,7,9,3,4,8\}$ ，则对应的染色体 $p=\{2,1,3,3,3,4,1,1,1\}$ 。在这种情况下，两个不同染色体在任意位置断点交叉后得到的新染色体都能构建出一个符合实际情况的路径。（若染色体中某个值m大于当前剩余的城市数目，则取m除以当前城市数目的余数；例如，若染色体 $p=\{1,1,9,1,1,1,1,1,1\}$ ，9大于其当前城市数7，9除以7余2，则这时染色体等价于 $p=\{1,1,2,1,1,1,1,1,1\}$ ）。

基于这样的思想，只需要对遗传算法的编码部分作出改动：

function parent = InitGroup(GroupNum, Dist)
    parent.fitness = [];
    parent.chrom = [];
    NumCity= size(Dist, 2);
    for i = 1:GroupNum
        city = 1 : NumCity;
        order = randperm(NumCity);%随机生成城市顺序
        item = [];
        for j = 1 : NumCity - 1
            item = [item, find(city == order(j))];
            city(find(city == order(j))) = [];
        end
        %按照上文规则，根据顺序生成染色体item
        item = [item, 1];
        item = num2str(item);
        item(find(item==' ')) = [];
        parent.chrom = [parent.chrom; item];
        %其实没必要转化成字符串
        %若城市数大于9（出现两位数三位数时）
        %采取单字节的字符串编码在交叉时就会出现问题
        %直接用double数组就好
    end
    parent.fitness = CalcFit(parent.chrom, Dist);
end

上段代码中本来为了节省空间，将编码转为字符串。但在调试的过程中发现当城市数大于9后两位数的城市代码会破坏原先字符串的单字节特性，使交叉操作出现错误。故应该直接采取double数字编码。

采用随机生成的数据，照样是当某解重复出现十次以上时认为算法收敛。应用自适应的交叉变异概率的遗传算法，可以得到某一次代码运行时，随机生成的TSP解的示意图：

贴出本代码之前，还有一些说明：
1. 本次计算中，由于TSP问题的特殊性，还加入了精英原则，即在选择操作时，保留当前适应度最高的个体。
2. 本代码中适应度计算（行走的总距离）时的解码Decode过程没有使用向量运算，而是直接循环嵌套。导致代码效率较低。有朋友想出如何用向量运算的话可以显著提高算法的运行速度。
3. 遗传算法得出的是较优解，不保证每次最优。

%%
%缂栫爜鏂瑰紡锟?鍩烘湰浜岃繘鍒剁紪锟?
%Input:     FitFunc:    Any function
%           pCrossover: probability of crossover, default 0.5
%           pMutation:  probability of mutation,  default 0.04
%           GroupNum:   number of individuals of the virtual group, default 30
%           MaxIter:    maximum iterations
%           MaxRepeat:  (optional)determine the convergence standard 鍒ゆ柇鏀舵暃
%parent.fitness
%parent.chrom

function iter = GA_TSP(pCrossover, pMutation, GroupNum, MaxIter, MaxRepeat)
    % Default parameters
    if nargin < 6
        MaxRepeat = 10;
        if nargin < 5
            MaxIter = 1000;
            if nargin < 4
                GroupNum = 10;
                if nargin < 3
                    pMutation = 0.03;
                    if nargin < 2
                        pCrossover = 0.45;
                    end
                end
            end
        end
    end
    NumCity = 9;
    CityPos = 10 * rand(NumCity, 2) + 10;
    Dist = zeros(NumCity, NumCity)
    for i = 1:NumCity
        for j = i+1:NumCity
            Dist(i, j) = sqrt((CityPos(i, 1) - CityPos(j, 1))^2 + (CityPos(i, 2) - CityPos(j, 2))^2);
            Dist(j, i) = Dist(i, j);
        end
    end
    Result  =  [];
    epsilon =  1e-4;
    iter    =  0;
    iRepeat =  1;
    thisMax =  0;
    parent  =  InitGroup(GroupNum, Dist);   %Generate initial population
    while iter < MaxIter
%% Heredity
        children1 = Crossover(parent, pCrossover/iter^0.1);  %Return crossovered chromes
        children.chrom = [];
        children.fitness = [];
        children.chrom = Mutation([parent.chrom; children1], pMutation/iter^0.1);
        children.fitness = CalcFit(children.chrom, Dist);
        children = select(children, GroupNum);
        parent = children;
        iter = iter + 1;
        %parent.chrom;
        %[m, I] = max(parent.fitness)
        if abs((thisMax-max(parent.fitness))/max(parent.fitness)) < epsilon
            (thisMax-max(parent.fitness))/max(parent.fitness);
            iRepeat = iRepeat + 1;
            if iRepeat == 10
                break
            end
        else 
            iRepeat = 1;
        end
        [thisMax, I] = max(parent.fitness);
%% Decode
        axis([0 22 0 22])
        order = [];
        city = 1 : NumCity;
        for j = 1 : NumCity 
                order = [order, city( 1 + rem(str2num(parent.chrom(I(1), j))-1, length(city)) )];
                city(1 + rem(str2num(parent.chrom(I(1), j))-1, length(city))) = [];
        end
        disp(200 - thisMax);
        Result = [Result; thisMax];
    end

 %% Draw
    scatter(CityPos(:,1), CityPos(:,2));
    hold on
    order
    for j = 1 : NumCity
        plot([CityPos(order(j),1) CityPos(order(rem(j, NumCity)+1),1)], [CityPos(order(j),2) CityPos(order(rem(j, NumCity)+1),2)],'r');
    end
    hold off
end

%Encoding method: 
function parent = InitGroup(GroupNum, Dist)
    parent.fitness = [];
    parent.chrom = [];
    NumCity= size(Dist, 2);
    for i = 1:GroupNum
        city = 1 : NumCity;
        order = randperm(NumCity);
        %order = num2str(order);
        %order(find(b-' '==0)) = [];
        item = [];
        for j = 1 : NumCity - 1
            item = [item, find(city == order(j))];
            city(find(city == order(j))) = [];
        end
        item = [item, 1];
        item = num2str(item);
        item(find(item==' ')) = [];
        parent.chrom = [parent.chrom; item];
    end
    parent.fitness = CalcFit(parent.chrom, Dist);
end

%Population = ['12345';
%              '23451';]


%Calculate Fitness
function Fitness = CalcFit(Population, Dist)
    NumCity = size(Dist, 2);
    PopNum = size(Population, 1);
    Fitness = zeros(PopNum, 1);
    %Decode:
    %city = 1 : NumCity;
    %city = repmat(city, GroupNum, 1);

    %Population
    for i = 1 : PopNum
        order = [];
        city = 1 : NumCity;
        for j = 1 : NumCity 
                order = [order, city( 1 + rem(str2num(Population(i, j))-1, length(city)) )];
                city(1 + rem(str2num(Population(i, j))-1, length(city))) = [];
        end
    end

    for i = 1 : NumCity
        city1 = order(i);
        city2 = order(1 + rem(i, NumCity));
        Fitness = Fitness + Dist((city1-1)*NumCity + city2);
    end
    Fitness = 200 - Fitness;
end

%roulette selcting method.  The same
function newPop = select(parent, PopNum)
    %Add: The best survive?
    cumFit = cumsum(parent.fitness)/sum(parent.fitness);
    [M, I] = max(parent.fitness);
    newPop.chrom(1,:) = parent.chrom(I, :);
    newPop.fitness(1) = parent.fitness(I);
    for i = 2 : PopNum
         index = find (cumFit - rand > 0);
         newPop.chrom(i,:) = parent.chrom(index(1),:);
         newPop.fitness(i) = parent.fitness(index(1));
    end
end


function childrenChrom = Crossover(parent, pCrossover)
    [PopNum, CityNum] = size(parent.chrom);
    childrenChrom = [];
    [M, I] = max(parent.fitness);
    childrenChrom = parent.chrom(I, :); %Parent with highest fitness
    for i = 1 : PopNum/2
        RandCross = rand(1);
        if RandCross < pCrossover
            i = fix(rand(1)*PopNum + 1);
            j = fix(rand(1)*PopNum + 1);
            while i == j
                i = fix(rand(1)*PopNum + 1);
                j = fix(rand(1)*PopNum + 1);
            end 
            BreakPoint = fix(rand(1)*CityNum + 1);
            temp = parent.chrom(i, 1:BreakPoint);
            parent.chrom(i, 1:BreakPoint) = parent.chrom(j, 1:BreakPoint);
            parent.chrom(j, 1:BreakPoint) = temp;
            childrenChrom = [childrenChrom; parent.chrom(i, :); parent.chrom(j, :)];
        end
    end
end

function childrenChrom = Mutation(chrom, pMutation)
    [PopNum, CityNum] = size(chrom);
    childrenChrom = chrom;
    for i = 1:PopNum
        for j = 1:CityNum
            if rand < pMutation
                childrenChrom(i, j) = num2str(floor(rand*CityNum)+1);
            end
        end
    end   
end