遗传，模拟退火，蚁群三个算法求解TSP的对比 - 图文

数学与统计学院

智能计算及应用课程设计

设计题目： 智能计算解决旅行商问题

摘要

本文以遗传算法、模拟退火、蚁群算法三个算法解决旅行商问题，将三个算法进行比较分析。目前这三个算法广泛应用于各个领域中，本文以31个城市为例，运用遗传算法、模拟退火、蚁群算法分别进行了计算，将他们的计算结果进行了比较分析。

关键词: 遗传算法 模拟退火 蚁群算法 旅行商问题

背景：

遗传算法：

20世纪60年代初，美国Michigan大学的John Holland教授开始研究自然和人工系统的自适应行为，在从事如何建立能学习的机器的研究过程中，受达尔文进化论的启发，逐渐意识到为获得一个好的算法仅靠单个策略建立和改进是不够的，还要依赖于一个包含许多候选策略的群体的繁殖，从而提出了遗传算法的基本思想。

20世纪60年代中期，基于语言智能和逻辑数学智能的传统人工智能十分兴盛，而基于自然进化思想的模拟进化算法则遭到怀疑与反对，但Holland及其指导的博士仍坚持这一领域的研究。Bagley发表了第一篇有关遗传算法应用的论文，并首先提出“遗传算法”这一术语，在其博士论文中采用双倍体编码，发展了复制、交叉、变异、显性、倒位等基因操作算子，并敏锐地察觉到防止早熟的机理，发展了自组织遗传算法的概念。与此同时，Rosenberg在其博士论文中进行了单细胞生物群体的计算机仿真研究，对以后函数优化颇有启发，并发展了自适应交换策略，在遗传操作方面提出了许多独特的设想。Hollistien在其1971年发表的《计算机控制系统的人工遗传自适应方法》论文中首次将遗传算法应用于函数优化，并对优势基因控制、交叉、变异以及编码技术进行了深入的研究。

人们经过长期的研究，在20世纪}o年代初形成了遗传算法的基本框架。1975年Holland出版了经典著作“Adaptation in Nature and Artificial System\该书详细阐述了遗传算

法的基本理论和方法，提出了著名的模式理论，为遗传算法奠定了数学基础。同年，DeJong发表了重要论文“An Analysis of the Behav-nor of a Class of Genetie Adaptive System\，在论文中，他将Holland的模式理论与计算实验结合起来，并通过函数优化的应用深人研究，将选择、交叉、变异操作进一步完善和系统化，并提出了代沟等新的操作技术，所得出的许多结论至今仍具有普遍的指导意义。

进入20世纪80年代末期，随着计算机技术的发展，遗传算法的研究再次兴起。遗传算法以其较强的解决问题的能力和广泛的适应性，深受众多领域研究人员的重视，遗传算法的理论研究和应用研究已成为十分热门的课题。自1985年起，遗传算法及其应用的国际会议每两年召开一次，并且在相关的人工智能会议和刊物上大多设有遗传算法的专题。

模拟退火：

模拟退火算法来源于固体退火原理，是一种基于概率的算法，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。

模拟退火算法(Simulated Annealing，SA)最早的思想是由N. Metropolis[1] 等人于1953年提出。1983 年,S. Kirkpatrick 等成功地将退火思想引入到组合优化领域。它是基于Monte-Carlo迭代求解策略的一种随机寻优算法，其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法从某一较高初温出发，伴随温度参数的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。模拟退火算法是一种通用的优化算法，理论上算法具有概率的全局优化性能,目前已在工程中得到了广泛应用，诸如VLSI、生产调度、控制工程、机器学习、神经网络、信号处理等领域。

模拟退火算法是通过赋予搜索过程一种时变且最终趋于零的概率突跳性，从而可有效避免陷入局部极小并最终趋于全局最优的串行结构的优化算法。

蚁群算法：

蚁群算法(ant colony optimization, ACO)，又称蚂蚁算法，是一种用来在图中寻找优化路径的机率型算法。它由Marco Dorigo于1992年在他的博士论文中提出，其灵感来源于蚂蚁在寻找食物过程中发现路径的行为。蚁群算法是一种模拟进化算法，初步的研究表明该算法具有许多优良的性质。针对PID控制器参数优化设计问题，将蚁群算法设计的结果与遗传算法设计的结果进行了比较，数值仿真结果表明，蚁群算法具有一种新的模拟进化优化方法的有效性和应用价值。

算法原理：

遗传算法：

遗传操作是模拟生物基因遗传的做法。在遗传算法中，通过编码组成初始群体后，遗

传操作的任务就是对群体的个体按照它们对环境适应度(适应度评估)施加一定的操作，从而实现优胜劣汰的进化过程。从优化搜索的角度而言，遗传操作可使问题

遗传过程

的解，一代又一代地优化，并逼近最优解。

遗传操作包括以下三个基本遗传算子(genetic operator):选择(selection)；交叉(crossover)；变异(mutation)。这三个遗传算子有如下特点：

个体遗传算子的操作都是在随机扰动情况下进行的。因此，群体中个体向最优解迁移的规则是随机的。需要强调的是，这种随机化操作和传统的随机搜索方法是有区别的。遗传操作进行的高效有向的搜索而不是如一般随机搜索方法所进行的无向搜索。

遗传操作的效果和上述三个遗传算子所取的操作概率，编码方法，群体大小，初始群体以及适应度函数的设定密切相关。

从群体中选择优胜的个体，淘汰劣质个体的操作叫选择。选择算子有时又称为再生算子(reproduction operator)。选择的目的是把优化的个体(或解)直接遗传到下一代或通过配对交叉产生新的个体再遗传到下一代。选择操作是建立在群体中个体的适应度评估基础上的，目前常用的选择算子有以

下几种：适应度比例方法、随机遍历抽样法、局部选择法。

其中轮盘赌选择法 （roulette wheel selection）是最简单也是最常用的选择方法。在该方法中，各个个体的选择概率和其适应度值成比例。设群体大小为n，其中个体i的适应度为，则i 被选择的概率，为遗传算法

显然，概率反映了个体i的适应度在整个群体的个体适应度总和中所占的比例。个体适应度越大。其被选择的概率就越高、反之亦然。计算出群体中各个个体的选择概率后，为了选择交配个体，需要进行多轮选择。每一轮产生一个[0，1]之间均匀随机数，将该随机数作为选择指针来确定被选个体。个体被选后，可随机地组成交配对，以供后面的交叉操作。

在自然界生物进化过程中起核心作用的是生物遗传基因的重组(加上变异)。同样，遗传算法中起核心作用的是遗传操作的交叉算子。所谓交叉是指把两个父代个体的部分结构加以替换重组而生成新个体的操作。通过交叉，遗传算法的搜索能力得以飞跃提高。

交叉算子根据交叉率将种群中的两个个体随机地交换某些基因，能够产生新的基因组合，期望将有益基因组合在一起。根据编码表示方法的不同，可以有以下的算法：实值重组（离散重组,中间重组,线性重组,扩展线性重组）。二进制交叉（单点交叉,多点交叉,均匀交叉,洗牌交叉,缩小代理交叉）。

最常用的交叉算子为单点交叉。具体操作是:在个体串中随机设定一个交叉点，实行交叉时，该点前或后的两个个体的部分结构进行互换，并生成两个新个体。下面给出了单点交叉的一个例子：

个体A：1 0 0 1 ↑1 1 1 → 1 0 0 1 0 0 0 新个体 个体B：0 0 1 1 ↑0 0 0 → 0 0 1 1 1 1 1 新个体

变异算子的基本内容是对群体中的个体串的某些基因座上的基因值作变动。依据个体编码表示方法的不同，可以有以下的算法：

a)实值变异 b)二进制变异。

一般来说，变异算子操作的基本步骤如下：

a)对群中所有个体以事先设定的变异概率判断是否进行变异 b)对进行变异的个体随机选择变异位进行变异。

遗传算法引入变异的目的有两个:一是使遗传算法具有局部的随机搜索能力。当遗传算法通过交叉算子已接近最优解邻域时，利用变异算子的这种局部随机搜索能力可以加速向最优解收敛。显然，此种情况下的变异概率应取较小值，否则接近最优解的积木块会因变异而遭到破坏。二是使遗传算法可维持群体多样性，以防止出现未成熟收敛现象。此时收敛概率应取较大值。

遗传算法中，交叉算子因其全局搜索能力而作为主要算子，变异算子因其局部搜索能力而作为辅助算子。遗传算法通过交叉和变异这对相互配合又相互竞争的操作而使其具备兼顾全局和局部的均衡搜索能力。所谓相互配合.是指当群体在进化中陷于搜索空间中某个超平面而仅靠交叉不能摆脱时，通过变异操作可有助于这种摆脱。所谓相互竞争，是指当通过交叉已形成所期望的积木块时，变异操作有可能破坏这些积木块。如何有效地配合使用交叉和变异操作，是目前遗传算法的一个重要研究内容。

模拟退火：

模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温

时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据Metropolis准则，粒子在温度T时趋于平衡的概率为e(-ΔE/(kT))，其中E为温度T时的内能，ΔE为其改变量，k为Boltzmann常数。用固体退火模拟组合优化问题，将内能E模拟为目标函数值f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schedule)控制，包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S。模拟退火算法可以分解为解空间、目标函数和初始解三部分。模拟退火的基本思想:(1) 初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起点)，每个T值的迭代次数L(2) 对k=1，??，L做第(3)至第6步：(3) 产生新解S′(4) 计算增量Δt′=C(S′)-C(S)，其中C(S)为评价函数，(5) 若Δt′<0则接受S′作为新的当前解，否则以概率exp(-Δt′/T)接受S′作为新的当前解.，(6) 如果满足终止条件则输出当前解作为最优解，结束程序。终止条件通常取为连续若干个新解都没有被接受时终止算法。(7) T逐渐减少，且T->0，然后转第2步。

模拟退火算法新解的产生和接受可分为如下四个步骤：

第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成

新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。 第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。 第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropolis准则: 若Δt′<0则接受S′作为新的当前解S，否则以概率exp(-Δt′/T)接受S′作为新的当前解S。 第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。 模拟退火算法与初始值无关，算法求得的解与初始解状态S(是算法迭代的起点)无关；模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率l 收敛于全局最优解的全局优化算法；模拟退火算法具有并行性。 蚁群算法： 蚁群优化算法是模拟蚂蚁觅食的原理，设计出的一种群集智能算法。蚂蚁在觅食过程中能够在其经过的路径上留下一种称之为信息素的物质，并在觅食过程中能够感知这种物质的强度，并指导自己行动方向，它们总是朝着该物质强度高的方向移动，因此大量蚂蚁组成的集体觅食就表现为一种对信息素的正反馈现象。某一条路径越短，路径上经过的蚂蚁越多，其信息素遗留的也就越多，信息素的浓度也就越高，蚂蚁选择这条路径的几率也就越高，由此构成的正反馈过程，从而逐渐的逼近最优路径，找到最优路径。 蚂蚁在觅食过程时，是以信息素作为媒介而间接进行信息交流，当蚂蚁从食物源走到蚁穴，或者从蚁穴走到食物源时，都会在经过的路径上释放信息素，从而形成了一条含有信息素的路径，蚂蚁可以感觉出路径上信息素浓度的大小，并且以较高的概率选择信息素浓度较高的路径 蚁群系统（ACS[12]）是由Dorigo等人提出来的改进的蚁群算法，它与AS的不同之处主要体现在三个方面：（1）采用不同的路径选择规则，能更好地利用蚂蚁所积累的搜索经验。（2）信息素挥发和信息素释放动作只在至今最优路径的边上执行，即每次迭代之后只有至今最优蚂蚁被允许释放信息素；（3）除了全局信息素更新规则外，还采用了局部信息素更新规则。 在ACS中，位于城市i的蚂蚁k，根据伪随机比例规则选择城市j作为下一个访问的城市。 设C={c1, c2, …, cn} 是n个城市的集合，L={lij|ci, cj C}是集合C中的元素（城市）两两连接的集合，dij(i, j=1,1,…,n)是lij的Euclidean距离，即 d ? (x?x)2?(y?y)2ijijijG=(C, L)是一个有向图，TSP的目的是从有向图G中寻出长度最短的Hamilton圈，即一条对C={c1, c2, …, cn}中n个元素（城市）访问且只访问一次的最短封闭曲线n城市规模的TSP，存在(n-1)!/2条不同闭合路径设bi(t)表示t时刻位于元素i的蚂蚁数目，τij (t)为t时刻路径(i, j)上的信息量，n表示TSP规模，m为蚁群中蚂蚁总数，则

m??bi(t)i?1n是t时刻集合C中元素（城市）两两连接lij上残留信息量的集合，在初始时刻各条路径上的信息量相等，并设τij(0)=const, 基本蚁群算法的寻优是通过有向图g=(C, L, Γ)实现的。 蚂蚁k(k=1,2,…,m)在运动过程中，根据各条路径上的信息量决定其转移方向。这里用禁忌表tabuk来记录蚂蚁k当前所走过的城市，集合随着tabuk进化过程做动态调整。在搜索过程中，蚂蚁根据各条路径上的信息量及路径的启发信息来计算状态转移概率。在t时刻蚂蚁k由元素（城市）i转移到元素（城市）j的状态转移概率： ?[?ij(t)]??[?ik(t)]?, if j?allowedk????k[?(t)]?[?(t)] p(t)??isis(1)?s?allowedijk? ??0 elsewise 其中allowedk={C-tabuk}表示蚂蚁k下一步允许选择的城市；α为信息启发式因子，表示轨迹的相对重要性，反映了蚂蚁在运动过程中积累的信息在蚂蚁运动时所起的作用，其值越大，则该蚂蚁越倾向于选择其它蚂蚁经过的路径，蚂蚁之间的协作性越强；β为期望启发式因子，表示能见度的相对重要性，反映蚂蚁在运动过程中启发信息在蚂蚁选择路径中的受重视程度，其值越大，则该状态状态转移概率越接近于贪心规则；ηij(t)为启发函数， ηij(t) =1/dij式中dij表示相邻两个城市之间的距离。对蚂蚁k而言，dij越小，则ηij(t)越大，pijk也就越大。显然，该启发函数表示蚂蚁从元素（城市）i转移到元素（城市）j的期望程度。 为了避免残留信息素过多引起残留信息淹没启发信息，在每只蚂蚁走完一步或者完成对所有n个城市的遍历(也即一个循环结束)后，要对残留信息进行更新处理。这种更新策略模仿了人类大脑记忆的特点，在新信息不断存人大脑的同时，存储在大脑中的旧信息随着时间的推移逐渐淡化，甚至忘记。由此，t+n时刻在路径(i, j)上的信息量可按如下规则进行调整 式中，ρ表示信息素挥发系数，则1-ρ表示信息素残留因子，为了防止信息的无限积累， ρ的取值范围为[0，1), Δτij(t)表示本次循环中路径(i, j)上的信息素增量，初始时刻Δτij(t) =0, Δτijk(t) 表示第k只蚂蚁在本次循环中留在路径(i, j)上的信息量。 根据信息素更新策略的不同，Dorigo M提出了三种不同的基本蚁群算法模型，分别称之为Ant-Cycle模型、Ant-Quantity模型及Ant-Density模型，其差别在于Δτijk(t)求法的不同。 ?ij(t?n)?(1??)??ij(t)???ij(t) (2)??ij(t)????ijk(t) (3)k?1m程序

遗传算法：

clear clc

close all %% 加载数据 %load data

X=[1304 2312;3639 1315;4177 2244;3721 1399;3488 1535;3326 1556;3238 1229;4196 1004;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 2367;3394 2643;3439 3201;2935 3240;3140 3550;2545 2357;2778 2826;2370 2975]; D=Distanse(X); %生成距离矩阵 N=size(D,1); %城市个数 %% 遗传参数

NIND=100; %种群大小 MAXGEN=200; %最大遗传代数 Pc=0.9; %交叉概率 Pm=0.05; %变异概率 GGAP=0.9; %代沟 %% 初始化种群

Chrom=InitPop(NIND,N); %% 画出随机解的路径图 DrawPath(Chrom(1,:),X) pause(0.0001)

%% 输出随机解的路径和总距离

disp('初始种群中的一个随机值:') OutputPath(Chrom(1,:));

Rlength=PathLength(D,Chrom(1,:)); disp(['总距离：',num2str(Rlength)]);

disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~') %% 优化 gen=0; figure;

hold on;box on xlim([0,MAXGEN]) title('优化过程') xlabel('代数') ylabel('最优值')

ObjV=PathLength(D,Chrom); %计算路径长度 preObjV=min(ObjV); while gen

ObjV=PathLength(D,Chrom); %计算路径长度 % fprintf('%d %1.10f\\n',gen,min(ObjV))

line([gen-1,gen],[preObjV,min(ObjV)]);pause(0.0001) preObjV=min(ObjV); FitnV=Fitness(ObjV); %% 选择

SelCh=Select(Chrom,FitnV,GGAP); %% 交叉操作

SelCh=Recombin(SelCh,Pc); %% 变异

SelCh=Mutate(SelCh,Pm); %% 逆转操作

SelCh=Reverse(SelCh,D); %% 重插入子代的新种群

Chrom=Reins(Chrom,SelCh,ObjV); %% 更新迭代次数 gen=gen+1 ; end

%% 画出最优解的路径图

ObjV=PathLength(D,Chrom); %计算路径长度 [minObjV,minInd]=min(ObjV); DrawPath(Chrom(minInd(1),:),X) %% 输出最优解的路径和总距离 disp('最优解:')

p=OutputPath(Chrom(minInd(1),:));

disp(['总距离：',num2str(ObjV(minInd(1)))]);

disp('-------------------------------------------------------------')

Distanse.m

%% 计算两两城市之间的距离 %输入 a 各城市的位置坐标 %输出 D 两两城市之间的距离 function D=Distanse(a) row=size(a,1); D=zeros(row,row); for i=1:row

for j=i+1:row

D(i,j)=((a(i,1)-a(j,1))^2+(a(i,2)-a(j,2))^2)^0.5; D(j,i)=D(i,j); end end

DrawPath.m %% 画路径函数 %输入

% Chrom 待画路径 % X 各城市坐标位置 function DrawPath(Chrom,X)

R=[Chrom(1,:) Chrom(1,1)]; %一个随机解(个体) figure; hold on

plot(X(:,1),X(:,2),'o','color',[0.5,0.5,0.5])

plot(X(Chrom(1,1),1),X(Chrom(1,1),2),'rv','MarkerSize',20) for i=1:size(X,1)

text(X(i,1)+0.05,X(i,2)+0.05,num2str(i),'color',[1,0,0]); end

A=X(R,:);

row=size(A,1); for i=2:row

[arrowx,arrowy] = dsxy2figxy(gca,A(i-1:i,1),A(i-1:i,2));%坐标转换 annotation('textarrow',arrowx,arrowy,'HeadWidth',8,'color',[0,0,1]); end

hold off

xlabel('横坐标') ylabel('纵坐标') title('轨迹图') box on

dsxy2figxy.m

function varargout = dsxy2figxy(varargin)

if length(varargin{1}) == 1 && ishandle(varargin{1}) ...

&& strcmp(get(varargin{1},'type'),'axes') hAx = varargin{1};

varargin = varargin(2:end); else

hAx = gca; end;

if length(varargin) == 1 pos = varargin{1}; else

[x,y] = deal(varargin{:}); end

axun = get(hAx,'Units');

set(hAx,'Units','normalized'); axpos = get(hAx,'Position'); axlim = axis(hAx);

axwidth = diff(axlim(1:2)); axheight = diff(axlim(3:4)); if exist('x','var')

varargout{1} = (x - axlim(1)) * axpos(3) / axwidth + axpos(1);

varargout{2} = (y - axlim(3)) * axpos(4) / axheight + axpos(2); else

pos(1) = (pos(1) - axlim(1)) / axwidth * axpos(3) + axpos(1); pos(2) = (pos(2) - axlim(3)) / axheight * axpos(4) + axpos(2); pos(3) = pos(3) * axpos(3) / axwidth; pos(4) = pos(4) * axpos(4 )/ axheight; varargout{1} = pos; end

set(hAx,'Units',axun) Fitness.m

function FitnV=Fitness(len) FitnV=1./len; InitPop.m %% 初始化种群 %输入：

% NIND：种群大小 % N： 城市的个数 %输出： %初始种群

function Chrom=InitPop(NIND,N) Chrom=zeros(NIND,N);%用于存储种群 for i=1:NIND

Chrom(i,:)=randperm(N);%随机生成初始种群 End

Mutate.m %% 变异操作 %输入：

%SelCh 被选择的个体 %Pm 变异概率 %输出：

% SelCh 变异后的个体

function SelCh=Mutate(SelCh,Pm) [NSel,L]=size(SelCh); for i=1:NSel if Pm>=rand

R=randperm(L);

SelCh(i,R(1:2))=SelCh(i,R(2:-1:1)); end end

OutputPat.m %% 输出路径函数 %输入：R 路径

function p=OutputPath(R) R=[R,R(1)];

N=length(R); p=num2str(R(1)); for i=2:N

p=[p,'—>',num2str(R(i))]; end disp(p)

PathLength.m

%% 计算各个体的路径长度 % 输入：

% D 两两城市之间的距离 % Chrom 个体的轨迹

function len=PathLength(D,Chrom) [row,col]=size(D); NIND=size(Chrom,1); len=zeros(NIND,1); for i=1:NIND

p=[Chrom(i,:) Chrom(i,1)]; i1=p(1:end-1); i2=p(2:end);

len(i,1)=sum(D((i1-1)*col+i2)); end

Recombin.m %% 交叉操作 % 输入

%SelCh 被选择的个体 %Pc 交叉概率 %输出：

% SelCh 交叉后的个体

function SelCh=Recombin(SelCh,Pc) NSel=size(SelCh,1);

for i=1:2:NSel-mod(NSel,2) if Pc>=rand %交叉概率Pc

[SelCh(i,:),SelCh(i+1,:)]=intercross(SelCh(i,:),SelCh(i+1,:)); end end

%输入：

%a和b为两个待交叉的个体 %输出：

%a和b为交叉后得到的两个个体

function [a,b]=intercross(a,b) L=length(a);

r1=randsrc(1,1,[1:L]);

r2=randsrc(1,1,[1:L]); if r1~=r2

a0=a;b0=b;

s=min([r1,r2]); e=max([r1,r2]); for i=s:e

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 end Reins.m

%% 重插入子代的新种群 %Chrom 父代的种群 %SelCh 子代种群 %ObjV 父代适应度

% Chrom 组合父代与子代后得到的新种群 function Chrom=Reins(Chrom,SelCh,ObjV) NIND=size(Chrom,1); NSel=size(SelCh,1);

[TobjV,index]=sort(ObjV);

Chrom=[Chrom(index(1:NIND-NSel),:);SelCh]; Reverse.m

%% 进化逆转函数 %SelCh 被选择的个体 %D 城市的距离矩阵 %输出

%SelCh 进化逆转后的个体

function SelCh=Reverse(SelCh,D) [row,col]=size(SelCh);

ObjV=PathLength(D,SelCh); %计算路径长度 SelCh1=SelCh; for i=1:row

r1=randsrc(1,1,[1:col]);

r2=randsrc(1,1,[1:col]); mininverse=min([r1 r2]); maxinverse=max([r1 r2]);

SelCh1(i,mininverse:maxinverse)=SelCh1(i,maxinverse:-1:mininverse); end

ObjV1=PathLength(D,SelCh1); %计算路径长度 index=ObjV1

SelCh(index,:)=SelCh1(index,:); Select.m

function SelCh=Select(Chrom,FitnV,GGAP) %% 选择操作 NIND=size(Chrom,1); % Chrom 种群 NSel=max(floor(NIND*GGAP+.5),2); % %GGAP：代沟 ChrIx=Sus(FitnV,NSel); %FitnV 适应度值

SelCh=Chrom(ChrIx,:); %SelCh 被选择的个体 Sus.m

%Nsel 被选择个体的数目

%NewChrIx 被选择个体的索引号

function NewChrIx = Sus(FitnV,Nsel) %FitnV 个体的适应度值 [Nind,ans] = size(FitnV); cumfit = cumsum(FitnV);

trials = cumfit(Nind) / Nsel * (rand + (0:Nsel-1)'); Mf = cumfit(:, ones(1, Nsel)); Mt = trials(:, ones(1, Nind))';

[NewChrIx, ans] = find(Mt < Mf & [ zeros(1, Nsel); Mf(1:Nind-1, :) ] <= Mt); [ans, shuf] = sort(rand(Nsel, 1)); NewChrIx = NewChrIx(shuf);

程序运行结果为:

初始轨迹：

优化过程：

优化后轨迹：

模拟退火:

clear clc

a = 0.99; % 温度衰减函数的参数 t0 = 97; tf = 3; t = t0; Markov_length = 10000; % Markov链长度

citys=[1 1304 2312;2 3639 1315;3 4177 2244;4 3721 1399;5 3488 1535;6 3326 1556;7 3238 1229;8 4196 1004;9 4312 790;10 4386 570;11 3007 1970;12 2562 1756;13 2788 1491;14 2381 1676;15 1332 695;16 3715 1678;17 3918 2179;18 4061 2370;19 3780 2212;20 3676 2578;21 4029 2838;22 4263 2931;23 3429 1908;24 3507 2367;25 3394 2643;26 3439 3201;27 2935 3240;28 3140 3550;29 2545 2357;30 2778 2826;31 2370 2975]; citys(:,1) = [];

amount = size(citys,1); % 城市的数目 % 通过向量化的方法计算距离矩阵 dist_matrix = zeros(amount, amount);

coor_x_tmp1 =citys(:,1) * ones(1,amount); coor_x_tmp2 = coor_x_tmp1';

coor_y_tmp1 = citys(:,2) * ones(1,amount); coor_y_tmp2 = coor_y_tmp1';

dist_matrix = sqrt((coor_x_tmp1-coor_x_tmp2).^2 + ... (coor_y_tmp1-coor_y_tmp2).^2);

sol_new = 1:amount; % 产生初始解

% sol_new是每次产生的新解；sol_current是当前解；sol_best是冷却中的最好解； E_current = inf; % E_current是当前解对应的回路距离； E_best = inf; % E_best是最优解的 % E_new是新解的回路距离；

sol_current = sol_new; sol_best = sol_new; p = 1;

while t>=tf

for r=1:Markov_length % Markov链长度 % 产生随机扰动

if (rand < 0.5) % 随机决定是进行两交换还是三交换 % 两交换

ind1 = 0; ind2 = 0; while (ind1 == ind2)

ind1 = ceil(rand.*amount); ind2 = ceil(rand.*amount); end

tmp1 = sol_new(ind1);

sol_new(ind1) = sol_new(ind2); sol_new(ind2) = tmp1; else

% 三交换

ind1 = 0; ind2 = 0; ind3 = 0;

while (ind1 == ind2) || (ind1 == ind3) ... || (ind2 == ind3) || (abs(ind1-ind2) == 1) ind1 = ceil(rand.*amount); ind2 = ceil(rand.*amount); ind3 = ceil(rand.*amount); end

tmp1 = ind1;tmp2 = ind2;tmp3 = ind3; % 确保ind1 < ind2 < ind3

if (ind1 < ind2) && (ind2 < ind3) ;

elseif (ind1 < ind3) && (ind3 < ind2) ind2 = tmp3;ind3 = tmp2;

elseif (ind2 < ind1) && (ind1 < ind3) ind1 = tmp2;ind2 = tmp1;

elseif (ind2 < ind3) && (ind3 < ind1)

ind1 = tmp2;ind2 = tmp3; ind3 = tmp1; elseif (ind3 < ind1) && (ind1 < ind2) ind1 = tmp3;ind2 = tmp1; ind3 = tmp2; elseif (ind3 < ind2) && (ind2 < ind1) ind1 = tmp3;ind2 = tmp2; ind3 = tmp1; end

tmplist1 = sol_new((ind1+1):(ind2-1)); sol_new((ind1+1):(ind1+ind3-ind2+1)) = ... sol_new((ind2):(ind3));

sol_new((ind1+ind3-ind2+2):ind3) = ... tmplist1; end

%检查是否满足约束

% 计算目标函数值（即内能） E_new = 0;

for i = 1 : (amount-1) E_new = E_new + ...

dist_matrix(sol_new(i),sol_new(i+1)); end

% 再算上从最后一个城市到第一个城市的距离 E_new = E_new + ...

dist_matrix(sol_new(amount),sol_new(1));

if E_new < E_current E_current = E_new; sol_current = sol_new; if E_new < E_best

% 把冷却过程中最好的解保存下来 E_best = E_new; sol_best = sol_new; end else

% 若新解的目标函数值小于当前解的， % 则仅以一定概率接受新解

if rand < exp(-(E_new-E_current)./t) E_current = E_new; sol_current = sol_new; else

sol_new = sol_current; end end end

t=t.*a; % 控制参数t（温度）减少为原来的a倍 end

disp('最优解为：') disp(sol_best)

disp('最短距离：') disp(E_best)

程序运行结果为:

蚁群算法:

clear all clc

%% 导入数据

citys=[1304 2312;3639 1315;4177 2244;3721 1399;3488 1535;3326 1556;3238 1229;4196 1004;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 2367;3394 2643;3439 3201;2935 3240;3140 3550;2545 2357;2778 2826;2370 2975]; %% 计算城市间相互距离 n = size(citys,1); D = zeros(n,n); for i = 1:n for j = 1:n if i ~= j

D(i,j) = sqrt(sum((citys(i,:) - citys(j,:)).^2)); else

D(i,j) = 1e-4; end end end

%% 初始化参数

m = 50; % 蚂蚁数量

alpha = 1; % 信息素重要程度因子 beta = 5; % 启发函数重要程度因子 rho = 0.1; % 信息素挥发因子 Q = 1; % 常系数 Eta = 1./D; % 启发函数 Tau = ones(n,n); % 信息素矩阵 Table = zeros(m,n); % 路径记录表 iter = 1; % 迭代次数初值 iter_max = 200; % 最大迭代次数 Route_best = zeros(iter_max,n); % 各代最佳路径 Length_best = zeros(iter_max,1); % 各代最佳路径的长度 Length_ave = zeros(iter_max,1); % 各代路径的平均长度

%% 迭代寻找最佳路径 while iter <= iter_max

% 随机产生各个蚂蚁的起点城市 start = zeros(m,1); for i = 1:m

temp = randperm(n); start(i) = temp(1); end

Table(:,1) = start; % 构建解空间 citys_index = 1:n;

% 逐个蚂蚁路径选择 for i = 1:m

% 逐个城市路径选择 for j = 2:n

tabu = Table(i,1:(j - 1)); % 已访问的城市集合(禁忌表) allow_index = ~ismember(citys_index,tabu);

allow = citys_index(allow_index); % 待访问的城市集合 P = allow;

% 计算城市间转移概率 for k = 1:length(allow)

P(k) = Tau(tabu(end),allow(k))^alpha * Eta(tabu(end),allow(k))^beta; end

P = P/sum(P);

% 轮盘赌法选择下一个访问城市 Pc = cumsum(P);

target_index = find(Pc >= rand); target = allow(target_index(1));

Table(i,j) = target; end end

% 计算各个蚂蚁的路径距离 Length = zeros(m,1); for i = 1:m

Route = Table(i,:); for j = 1:(n - 1)

Length(i) = Length(i) + D(Route(j),Route(j + 1)); end

Length(i) = Length(i) + D(Route(n),Route(1)); end

% 计算最短路径距离及平均距离 if iter == 1

[min_Length,min_index] = min(Length); Length_best(iter) = min_Length; Length_ave(iter) = mean(Length);

Route_best(iter,:) = Table(min_index,:); else

[min_Length,min_index] = min(Length);

Length_best(iter) = min(Length_best(iter - 1),min_Length); Length_ave(iter) = mean(Length); if Length_best(iter) == min_Length

Route_best(iter,:) = Table(min_index,:); else

Route_best(iter,:) = Route_best((iter-1),:); end end

% 更新信息素

Delta_Tau = zeros(n,n); % 逐个蚂蚁计算 for i = 1:m

% 逐个城市计算 for j = 1:(n - 1)

Delta_Tau(Table(i,j),Table(i,j+1)) = Delta_Tau(Table(i,j),Table(i,j+1)) Q/Length(i); end

Delta_Tau(Table(i,n),Table(i,1)) = Delta_Tau(Table(i,n),Table(i,1)) + Q/Length(i); end

Tau = (1-rho) * Tau + Delta_Tau; % 迭代次数加1，清空路径记录表 iter = iter + 1;

Table = zeros(m,n); end

+

%% 结果显示

[Shortest_Length,index] = min(Length_best); Shortest_Route = Route_best(index,:);

disp(['最短距离:' num2str(Shortest_Length)]);

disp(['最短路径:' num2str([Shortest_Route Shortest_Route(1)])]);

%% 绘图 figure(1)

plot([citys(Shortest_Route,1);citys(Shortest_Route(1),1)],... [citys(Shortest_Route,2);citys(Shortest_Route(1),2)],'o-'); grid on

for i = 1:size(citys,1)

text(citys(i,1),citys(i,2),[' ' num2str(i)]); end

text(citys(Shortest_Route(1),1),citys(Shortest_Route(1),2),' 起点'); text(citys(Shortest_Route(end),1),citys(Shortest_Route(end),2),' 终点'); xlabel('城市位置横坐标') ylabel('城市位置纵坐标')

title(['蚁群算法优化路径(最短距离:' num2str(Shortest_Length) ')']) figure(2)

plot(1:iter_max,Length_best,'b',1:iter_max,Length_ave,'r:') legend('最短距离','平均距离') xlabel('迭代次数') ylabel('距离')

title('各代最短距离与平均距离对比')

程序运行结果（蚁群数:50,迭代数200）：

各代最短距离与平均距离：

结论分析：

根据上面算法得出的结论可知：

遗传算法在本例中得出的最短距离为：15426.5219 模拟退火在本例中得出的最短距离为：1.5371e+004 蚁群算法在本例中得出的最短距离为：15601.0734

由此可见模拟退火在本例中最优，遗传算法次之。但是根据最后的结果也可以看出三个算法求解的结果相差不大。

参考文献：

[1]",张",挺.遗传算法在旅行商及网络优化问题中的研究与应用[D].太原:太原理工大学,2008

[2]????.一种改进的遗传算法及其在求解旅行商问题中的应用[J].电脑开发与应用,2008,7(22):35-36.

[3]????,倪志伟,刘慧婷.求解旅行商问题的一个改进的遗传算法[J].计算机工程与应用,2007,6(43):65-68.

[4]苗卉,杨韬.旅行商问题(TSP)的改进模拟退火算法.微计算机信息(管控一体化),2007,23(11):241

[5]吴小菁. 求解旅行商问题的模拟进化算法.福建金融管理干部学院学报,2008,(5):55-56,57.

[6]万军洲.基于模拟退火技术的旅行商问题求解算法.软件导刊,2006,(8):88,88-89.

[7]曲强，陈雪波.基于MATLAB的模拟退火算法的实现.鞍山科技大学学报,2003,26(3):197-198

[8]马良,朱刚,宁爱兵著,蚁群优化算法；2008

[9]李士勇.蚁群算法及其应用[M].哈尔滨：哈尔滨工业大学出版社，2004

[10]周康，强小利，同小军等.求解TSP算法[J].计算机工程与应用，2007，43（29）：43-47.

本文来源：https://www.bwwdw.com/article/kp4w.html