分布估计算法解决旅行商问题(TSP)

在用分布估计算法解决旅行商问题时,结构与传统的分布估计算法相似,只不过是把概率向量换成了“概率矩阵”而已:

  1. 通过概率矩阵生成解。
  2. 评估解。
  3. 选择优势群体。
  4. 更新概率矩阵。
  5. 重复以上4步直至迭代结束。

这里说的“概率矩阵”记录了上一代优势群体中,“城市对”出现的次数(或与城市对出现次数成正比的一个数)。“城市对”是指路线中相邻的两个城市,这两个城市不分先后。假设概率矩阵第i行第j列的元素为 p_{ij},它代表在优势群体中城市i和城市j相邻的次数(即i−j和j−i发生的次数)。

经过测试(于2018-11-09),这种算法非常容易陷入局部最优,效果不是很好。

下面是程序运行结果:

每代最短距离变化图

每代平均距离变化图

最短路线图

主程序

主程序如下:

function main()
pos = [randperm(20); randperm(20)];
% pos = load('berlin52.txt'); % 7542
% pos = pos(:, 2:3);
% pos = pos';

popsize = 300;                  % 种群规模
maxgen = 100;                   % 最大迭代次数
citysize = size(pos, 2);        % 城市数量
p = ones(citysize, citysize);   % 概率矩阵
sn = ceil(popsize * 0.5);       % 优势群体规模
dm = makedm(pos);               % 距离矩阵

bestlists = zeros(maxgen, citysize);    % 记录每代最优解
bestfits = zeros(1, maxgen);            % 记录每代最优解适应度值
avgfits = zeros(1, maxgen);             % 记录每代平均适应度值

for gen = 1:maxgen
    pop = makepop(popsize, citysize, p);    % 制作种群

    fitness = callength(pop, dm);
    [bfit, bfiti] = min(fitness);
    bestlists(gen, :) = pop(bfiti, :);
    bestfits(1, gen) = bfit;
    avgfits(1, gen) = sum(fitness) / popsize;

    spop = selection(pop, sn, fitness);     % 选择优势群体
    p = makep(spop);                        % 更新概率矩阵
end
fitness = callength(bestlists, dm);
[bfit, bfiti] = min(fitness);
fprintf("最短距离: %f\n", bfit);
end

制作种群函数

对于每个个体,第一个城市是随机生成的,之后每个城市都根据上一个城市进行轮盘赌法选择。有时候上一代优势群体中可能不会出现所有可能的城市对,这意味着该城市不会出现在候选城市中,但是路径必须经过每一个城市,应该处理这种情况,这种情况也可以通过给概率矩阵设置最小值以避免。

function pop = makepop(popsize, citysize, p)
% 制作种群
% popsize       input  种群规模
% citysize      input  城市数量
% p             input  概率矩阵
% pop           output 种群
pop = zeros(popsize, citysize);

for i = 1:popsize
    for j = 1:citysize
        if j == 1   % 随机产生第一个城市
            pop(i, j) = randi(citysize);
        else        % 使用轮盘赌法选择下一个城市
            cityp = p(pop(i, j-1), :);      % 找出备选城市
            cityp(pop(i, 1:j-1)) = 0;       % 已经访问过的概率为0
            sumcityp = sum(cityp);          % 开始轮盘赌法
            cityp = cityp ./ sumcityp;
            cityp = cumsum(cityp);
            index = find(cityp >= rand);

            % 上一代优势群体中没有城市对记录,则找一个城市填充
            if isempty(index)
                tempcity = 1:1:citysize;
                tempcity(pop(i, 1:j-1)) = 0;
                tempindex = find(tempcity > 0);
                index = tempindex(1);
            end
            pop(i, j) = index(1);
        end
    end
end
end

更新概率矩阵函数

根据上面对概率矩阵的定义更新概率矩阵:

function p = makep(spop)
% 更新概率矩阵
% spop          input  优势群体
% np            output 概率矩阵
[popsize, citysize] = size(spop);
count = zeros(citysize, citysize);   % 记录出现次数
spop = [spop spop(:, 1)];
for i = 1:popsize
    for j = 1:citysize
        a = spop(i, j);
        b = spop(i, j+1);
        count(a, b) = count(a, b) + 1;
        count(b, a) = count(b, a) + 1;
    end
end
sumcount = sum(count, 2);
p = count ./ sumcount;
end

一些其他函数

选择函数:

function spop = selection(pop, sn, fitness)
% 选择
% pop       input  种群
% sn        input  优势群体规模
% fitness   input  种群适应度值
% spop      output 优势群体
[~, index] = sort(fitness);
index = index(1: sn);
spop = pop(index, :);
end

计算路径长度函数:

function len = callength(pop, dm)
% 求路径距离
% pop           input  种群
% dm            input  距离矩阵
% len           output 距离
[NP, D] = size(pop);
pop = [pop pop(:,1)];
sum = zeros(NP, 1);
for i = 1:NP
    for j = 1:D
        sum(i,1) = sum(i,1) + dm(pop(i, j), pop(i, j+1));
    end
end
len = sum;
end

制作距离矩阵函数:

function dm = makedm(pos)
% 制作距离矩阵
% pos       input  城市坐标
% dm        output 距离矩阵
[~, len] = size(pos);
deltax = repmat(pos(1,:)', 1, len) - repmat(pos(1,:), len, 1);
deltay = repmat(pos(2,:)', 1, len) - repmat(pos(2,:), len, 1);
dm = round((deltax .^ 2 + deltay .^ 2) .^ 0.5);
end

留下评论

电子邮件地址不会被公开。 必填项已用*标注