文章基于 遗传算法--旅行商问题(TSP问题)-Matlab 上进行,优化了计算的时间成本问题
1、旅行商问题(TSP问题

假设有一个旅行商人要拜访全国31个省会城市,它需要选择所要走的路径,路径的限制是每个城市只能拜访一次,而且最后要回到原来出发的城市。对路径选择的要求是:所选路径的路成为所有路径之中的最小值。
全国31个省会城市的坐标为

[1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3326 1556;3238 1229;4196 1044;

4312  790;4386  570;3007 1970;2562 1756;2788 1491;2381 1676;1332  695;3715 1678;

3918 2179;4061 2370;3780 2212;3676 2578;4029 2838;4263 2931;3429 1908;3507 2376;

3394 2643;3439 3201;2935 3240;3140 3550;2545 2357;2778 2826;2370 2975]

2、仿真过程

原思路:

(1)初始化种群数目NP=200,染色体基因维数为N=31,最大进化代数G=1000.

(2)产生初始种群,计算个体适应度值,即路径长度:采用基于概率的方式选择进行操作的个体;对选中的成对个体,随机交叉所选中的成对城市坐标,以确保交叉后路径每个城市只到访一次;对选中的单个个体,随机交换其一对城市坐标作为变异操作,产生新的种群,进行下一次遗传操作。

(3)判断是否满足终止条件:若满足,则结束搜索过程,输出优化值,若不满足,则继续进迭代优化。

(4)为迭代过程预设一个文件,用于存储保存每个迭代的图像帧(仅当当前图与上次不同),制作迭代过程图

clear all;                      % 清除所有变量
close all;                      % 清图
clc;                            % 清屏

C = [1304 2312; 3639 1315; 4177 2244; 3712 1399; 3488 1535; 3326 1556; ...
    3238 1229; 4196 1044; 4312  790; 4386  570; 3007 1970; 2562 1756; ...
    2788 1491; 2381 1676; 1332  695; 3715 1678; 3918 2179; 4061 2370; ...
    3780 2212; 3676 2578; 4029 2838; 4263 2931; 3429 1908; 3507 2376; ...
    3394 2643; 3439 3201; 2935 3240; 3140 3550; 2545 2357; 2778 2826; ...
    2370 2975];                 %31个省会城市坐标
N = size(C,1);                    % TSP问题的规模, 即城市数目
D = zeros(N);                     % 任意两个城市距离间隔矩阵

% 计算任意两个城市的距离
for i = 1:N
    for j = 1:N
        D(i,j) = sqrt((C(i,1) - C(j,1))^2 + (C(i,2) - C(j,2))^2);
    end
end

NP = 200;                          % 种群规模
G = 2500;                          % 最大遗传代数
f = zeros(NP, N);                   % 用于存储种群
F = [];                            % 种群更新中间存储
for i = 1:NP
    f(i,:) = randperm(N);          % 随机生成初始种群
end
R = f(1,:);                        % 存储最优种群
len = zeros(NP, 1);                 % 存储路径长度
fitness = zeros(NP, 1);             % 存储归一化适应值
gen = 0;

% 为动画预设一个文件
video = VideoWriter('TSP_GA_Animation', 'MPEG-4');
open(video);

% 上一次保存的路径图
last_frame = [];

% 遗传算法循环
while gen < G
    % 计算路径长度
    for i = 1:NP
        len(i,1) = D(f(i,N), f(i,1));
        for j = 1:(N-1)
            len(i,1) = len(i,1) + D(f(i,j), f(i,j+1));
        end
    end
    maxlen = max(len);              % 最长路径
    minlen = min(len);              % 最短路径
    
    % 更新最短路径
    rr = find(len == minlen);
    R = f(rr(1,1),:);
    
    % 计算归一化适应值
    for i = 1:length(len)
        fitness(i,1) = (1 - ((len(i,1) - minlen) / (maxlen - minlen + 0.001)));
    end
    
    % 选择操作
    nn = 0;
    for i = 1:NP
        if fitness(i,1) >= rand
            nn = nn + 1;
            F(nn,:) = f(i,:);
        end
    end
    [aa, bb] = size(F);
    
    while aa < NP
        nnper = randperm(nn);
        A = F(nnper(1), :);
        B = F(nnper(2), :);
        
        % 交叉操作
        W = ceil(N / 10);              % 交叉点个数
        p = unidrnd(N - W + 1);          % 随机选择交叉范围,从p到p+W
        for i = 1:W
            x = find(A == B(p+i-1));
            y = find(B == A(p+i-1));
            temp = A(p+i-1); 
            A(p+i-1) = B(p+i-1); 
            B(p+i-1) = temp;
            temp = A(x); 
            A(x) = B(y); 
            B(y) = temp;
        end
        
        % 变异操作
        p1 = floor(1 + N * rand());
        p2 = floor(1 + N * rand());
        while p1 == p2
            p1 = floor(1 + N * rand());
            p2 = floor(1 + N * rand());
        end
        tmp = A(p1); 
        A(p1) = A(p2); 
        A(p2) = tmp;
        tmp = B(p1); 
        B(p1) = B(p2); 
        B(p2) = tmp;
        F = [F; A; B];
        [aa, bb] = size(F);
    end
    
    if aa > NP
        F = F(1:NP,:);             %保持种群规模为n
    end
    
    f = F;                         % 更新种群
    f(1,:) = R;                    % 保留每代最优个体
    clear F;
    gen = gen + 1;
    
    % 记录最优路径长度
    Rlength(gen) = minlen;
    
    % 绘制当前的路径,并添加迭代次数
    figure(1);
    clf;
    hold on;
    for i = 1:N-1
        plot([C(R(i),1), C(R(i+1),1)], [C(R(i),2), C(R(i+1),2)], 'bo-');
    end
    plot([C(R(N),1), C(R(1),1)], [C(R(N),2), C(R(1),2)], 'ro-');
    title(['优化最短距离: ', num2str(minlen)]);
    
    % 添加当前迭代次数到图形上
    text(0.5, 0.05, ['Iteration: ', num2str(gen)], 'Units', 'normalized', 'HorizontalAlignment', 'center', 'FontSize', 12);

    % 保存每个迭代的图像帧(仅当当前图与上次不同)
    frame = getframe(gcf);
    if ~isequal(last_frame, frame)
        writeVideo(video, frame);
        last_frame = frame;  % 更新上次帧
    end
end

% 关闭视频文件
close(video);

% 绘制适应度进化曲线
figure(2);
plot(Rlength);
xlabel('迭代次数');
ylabel('目标函数值');
title('适应度进化曲线');

TSP_GA_Animation-1

一个在2000次迭代后优化到1700以下的实例

代码改进思路 

通过引入更复杂的交叉和变异操作,改进适应度函数,并增强种群多样性,可以大幅提高遗传算法解决TSP问题的效率和质量。此外,通过并行计算加速路径计算,可以显著减少算法的运行时间,尤其是在大规模问题中。

        交叉和变异操作的改进: 当前的交叉操作是随机选择一段基因进行交换,可以考虑使用更复杂的交叉策略,如部分映射交叉(PMX)或顺序交叉(OX)。 变异操作可以增加更多的变异策略,如逆转变异或插入变异。

  • 部分映射交叉(PMX):这种方法通过在父代中选择一个子串,保留父代的一部分,并且调整剩余部分以保持合法的城市顺序。
  • 顺序交叉(OX):在这个方法中,首先选择父代中的一部分基因,然后在另一个父代中填充剩余部分,保持原有顺序。
    % 逆转变异
    function mutated = reverse_mutation(route)
        N = length(route);
        p1 = randi([1, N]);
        p2 = randi([p1, N]);
        
        mutated = route;
        mutated(p1:p2) = route(p2:-1:p1);
    end
    
    % 插入变异
    function mutated = insertion_mutation(route)
        N = length(route);
        p1 = randi([1, N]);
        p2 = randi([1, N]);
        while p1 == p2
            p2 = randi([1, N]);
        end
        
        mutated = route;
        gene = mutated(p1);
        mutated(p1) = mutated(p2);
        mutated(p2) = gene;
    end
    

        适应度函数的改进: 当前的适应度函数是简单的归一化处理,可以考虑使用指数函数或其他非线性函数来增强选择压力。这样做能够使得适应度较高的个体拥有更大的选择概率。

  • 计算归一化适应值,使用指数函数增强选择压力
    % 计算归一化适应值,使用指数函数增强选择压力
    function fitness = compute_fitness(len, minlen, maxlen)
        fitness = exp(-((len - minlen) / (maxlen - minlen + 0.001)));
    end
    

        种群多样性维护: 在进化过程中,种群可能会过早收敛。可以通过引入多样性维护机制,如小生境技术或移民策略,来保持种群的多样性。 并行计算: 遗传算法的种群评估可以并行化,以加速计算过程。

  • 小生境技术:通过将种群分成多个子种群,每个子种群内部进行进化,定期将优秀个体从一个子种群迁移到另一个子种群,保持种群多样性。
  • 移民策略:每隔一段时间将一部分个体从外部加入到种群中,防止种群过早收敛。
    % 种群多样性维护:移民策略
    function new_population = immigration_strategy(population, best_individuals)
        immigration_rate = 0.1;  % 迁移率
        num_immigrants = round(immigration_rate * size(population, 1));
        
        new_population = population;
        for i = 1:num_immigrants
            immigrant = best_individuals(randi([1, length(best_individuals)]), :);
            new_population(randi([1, size(new_population, 1)]), :) = immigrant;
        end
    end
    

    并行计算:加速种群评估

    可以使用 MATLAB 的并行计算工具箱来加速路径计算,尤其是当种群较大时。你可以将适应度评估部分并行化。

    % 并行计算种群适应度
    parfor i = 1:NP
        len(i, 1) = D(f(i, N), f(i, 1));
        for j = 1:(N - 1)
            len(i, 1) = len(i, 1) + D(f(i, j), f(i, j + 1));
        end
    end
    

    整合后代码

    clear all;                      % 清除所有变量 
    close all;                      % 清图
    clc;                            % 清屏
    
    C = [1304 2312; 3639 1315; 4177 2244; 3712 1399; 3488 1535; 3326 1556; ...
        3238 1229; 4196 1044; 4312  790; 4386  570; 3007 1970; 2562 1756; ...
        2788 1491; 2381 1676; 1332  695; 3715 1678; 3918 2179; 4061 2370; ...
        3780 2212; 3676 2578; 4029 2838; 4263 2931; 3429 1908; 3507 2376; ...
        3394 2643; 3439 3201; 2935 3240; 3140 3550; 2545 2357; 2778 2826; ...
        2370 2975];                 % 31个省会城市坐标
    N = size(C,1);                    % TSP问题的规模, 即城市数目
    D = zeros(N);                     % 任意两个城市距离间隔矩阵
    
    % 计算任意两个城市的距离
    for i = 1:N
        for j = 1:N
            D(i,j) = sqrt((C(i,1) - C(j,1))^2 + (C(i,2) - C(j,2))^2);  % 使用欧几里得距离计算城市间距离
        end
    end
    
    NP = 200;                          % 种群规模
    G = 2500;                          % 最大遗传代数
    f = zeros(NP, N);                   % 用于存储种群
    F = [];                            % 种群更新中间存储
    for i = 1:NP
        f(i,:) = randperm(N);          % 随机生成初始种群
    end
    R = f(1,:);                        % 存储最优种群
    len = zeros(NP, 1);                 % 存储路径长度
    fitness = zeros(NP, 1);             % 存储归一化适应值
    gen = 0;
    
    % 为动画预设一个文件
    video = VideoWriter('TSP_GA_Animation', 'MPEG-4');
    open(video);
    
    % 上一次保存的路径图
    last_frame = [];
    
    % 遗传算法循环
    while gen < G
        % 计算路径长度
        parfor i = 1:NP
            len(i,1) = D(f(i,N), f(i,1));  % 计算从末尾城市到起始城市的路径长度
            for j = 1:(N-1)
                len(i,1) = len(i,1) + D(f(i,j), f(i,j+1));  % 计算路径中其他城市之间的总距离
            end
        end
        maxlen = max(len);              % 最长路径
        minlen = min(len);              % 最短路径
        
        % 更新最短路径
        rr = find(len == minlen);  % 找到最短路径
        R = f(rr(1,1),:);  % 存储最优路径
        
        % 计算归一化适应值,使用指数函数增强选择压力
        fitness = exp(-((len - minlen) / (maxlen - minlen + 0.001)));  % 将路径长度映射为适应度
        
        % 选择操作
        nn = 0;
        for i = 1:NP
            if fitness(i,1) >= rand  % 适应度大于随机值则选择
                nn = nn + 1;
                F(nn,:) = f(i,:);  % 将选择的个体加入中间种群
            end
        end
        [aa, bb] = size(F);
        
        while aa < NP  % 保证新种群的大小为NP
            nnper = randperm(nn);  % 随机排列
            A = F(nnper(1), :);
            B = F(nnper(2), :);
            
            % 部分映射交叉 (PMX)
            [A, B] = PMX_crossover(A, B);  % 交叉操作
            
            % 变异操作:逆转变异
            A = reverse_mutation(A);  % 对个体A进行逆转变异
            B = reverse_mutation(B);  % 对个体B进行逆转变异
            
            % 移民策略(增强种群多样性)
            best_individuals = f(1:round(0.1*NP), :);  % 选择前10%最优个体
            F = immigration_strategy(F, best_individuals);  % 增加种群的多样性
            
            F = [F; A; B];  % 将交叉后的两个个体加入到新种群中
            [aa, bb] = size(F);
        end
        
        if aa > NP
            F = F(1:NP,:);             % 保持种群规模为NP
        end
        
        f = F;                         % 更新种群
        f(1,:) = R;                    % 保留每代最优个体
        clear F;
        gen = gen + 1;
        
        % 记录最优路径长度
        Rlength(gen) = minlen;
        
        % 绘制当前的路径,并添加迭代次数
        figure(1);
        clf;
        hold on;
        for i = 1:N-1
            plot([C(R(i),1), C(R(i+1),1)], [C(R(i),2), C(R(i+1),2)], 'bo-');  % 绘制路径
        end
        plot([C(R(N),1), C(R(1),1)], [C(R(N),2), C(R(1),2)], 'ro-');  % 连接最后一个城市与第一个城市
        title(['优化最短距离: ', num2str(minlen)]);  % 显示当前的最短路径长度
        
        % 添加当前迭代次数到图形上
        text(0.5, 0.05, ['Iteration: ', num2str(gen)], 'Units', 'normalized', 'HorizontalAlignment', 'center', 'FontSize', 12);
    
        % 保存每个迭代的图像帧(仅当当前图与上次不同)
        frame = getframe(gcf);
        if ~isequal(last_frame, frame)
            writeVideo(video, frame);  % 写入视频文件
            last_frame = frame;  % 更新上次帧
        end
    end
    
    % 关闭视频文件
    close(video);
    
    % 绘制适应度进化曲线
    figure(2);
    plot(Rlength);
    xlabel('迭代次数');
    ylabel('目标函数值');
    title('适应度进化曲线');
    
    % 部分映射交叉 (PMX)
    function [offspring1, offspring2] = PMX_crossover(parent1, parent2)
        N = length(parent1);
        point1 = randi([1, N-1]);  % 随机选择交叉点
        point2 = randi([point1+1, N]);
        
        offspring1 = parent1;
        offspring2 = parent2;
        
        % 交换选定区间的部分
        for i = point1:point2
            temp = offspring1(i);
            offspring1(i) = offspring2(i);
            offspring2(i) = temp;
        end
        
        % 修复交叉后的冲突
        for i = 1:N
            if ismember(offspring1(i), offspring1(point1:point2))
                for j = 1:N
                    if ~ismember(parent2(j), offspring1)
                        offspring1(i) = parent2(j);
                        break;
                    end
                end
            end
        end
        
        for i = 1:N
            if ismember(offspring2(i), offspring2(point1:point2))
                for j = 1:N
                    if ~ismember(parent1(j), offspring2)
                        offspring2(i) = parent1(j);
                        break;
                    end
                end
            end
        end
    end
    
    % 逆转变异
    function mutated = reverse_mutation(route)
        N = length(route);
        p1 = randi([1, N]);  % 随机选择变异区间
        p2 = randi([p1, N]);
        
        mutated = route;
        mutated(p1:p2) = route(p2:-1:p1);  % 将选择的部分反转
    end
    
    % 移民策略(增强种群多样性)
    function new_population = immigration_strategy(population, best_individuals)
        immigration_rate = 0.1;  % 迁移率
        num_immigrants = round(immigration_rate * size(population, 1));  % 计算迁移个体数量
        
        % 防止超出索引
        num_best = size(best_individuals, 1);  % best_individuals 的个体数
        new_population = population;
        
        for i = 1:num_immigrants
            immigrant_idx = randi([1, min(num_best, length(best_individuals))]);
            immigrant = best_individuals(immigrant_idx, :);  % 选择迁移的个体
            new_population(randi([1, size(new_population, 1)]), :) = immigrant;  % 将其加入新种群中
        end
    end
    

  • 一些在600次以内优化到17000一下的实例

TSP_GA_after_Animation-2

结论:遗传算法解决TSP问题难免陷入局部最优解的境地,往往与理想中的精确解有很大一段距离。

matlab实现的分支切割算法(Branch-and-Cut) 求解TSP问题的精确解

Logo

魔乐社区(Modelers.cn) 是一个中立、公益的人工智能社区,提供人工智能工具、模型、数据的托管、展示与应用协同服务,为人工智能开发及爱好者搭建开放的学习交流平台。社区通过理事会方式运作,由全产业链共同建设、共同运营、共同享有,推动国产AI生态繁荣发展。

更多推荐