遗传算法（三）解TSP问题

巡回旅行商问题TSP

TSP问题：给定一个图G=(V,E)，每个边e∈E 中有非负权值ω(e)，寻找G的Hamilton圈C，使C 的总权W(C)=∑e∈E(C)ω(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.chromfunction 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 offend%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 Fitnessfunction 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 samefunction 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));    endendfunction 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    endendfunction 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