数学与统计学院

数学与统计学院智能计算及应用课程设计设计题目：智能计算解决旅行商问题摘要：旅行商社问题是研究最为广泛地组合优化问题，在现实生活中，也有着广泛地应用。遗传算法是一种有效的最优化问题的方法。遗传算法是求解NP完全问题的一种常用的方法，它在解决排列组合问题方面占有很重要的地位。蚁群算法是搜索良好的一种启发式搜索方法。模拟退火是一种基于概率的算法，因此，本文将采用这三种算法来解决旅行商社问题，并得出结果进行分析。关键词：旅行商社问题遗传算法蚁群算法模拟退火背景阐述：1.遗传算法遗传算法是模拟达尔文生物进化论的自然选择和遗传学机理的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法。遗传算法是从代表问题可能潜在的解集的一个种群开始的，而一个种群则由经过基因编码的一定数目的个体组成。每个个体实际上是染色体带有特征的实体。染色体作为遗传物质的主要载体，即多个基因的集合，其内部表现（即基因型）是某种基因组合，它决定了个体的形状的外部表现，如黑头发的特征是由染色体中控制这一特征的某种基因组合决定的。因此，在一开始需要实现从表现型到基因型的映射即编码工作。由于仿照基因编码的工作很复杂，我们往往进行简化，如二进制编码，初代种群产生之后，按照适者生存和优胜劣汰的原理，逐代演化产生出越来越好的近似解，在每一代，根据问题域中个体的适应度大小选择个体，并借助于自然遗传学的遗传算子进行组合交叉和变异，产生出代表新的解集的种群。这个过程将导致种群像自然进化一样的后生代种群比前代更加适应于环境，末代种群中的最优个体经过解码，可以作为问题近似最优解。模拟退火模拟退火算法最早的思想是由N.Metropolis[1]等人于1953年提出。1983年,S.Kirkpatrick等成功地将退火思想引入到组合优化领域。它是基于Monte-Carlo迭代求解策略的一种随机寻优算法，其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法从某一较高初温出发，伴随温度参数的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。模拟退火算法是一种通用的优化算法，理论上算法具有概率的全局优化性能,目前已在工程中得到了广泛应用，诸如VLSI、生产调度、控制工程、机器学习、神经网络、信号处理等领域。模拟退火算法是通过赋予搜索过程一种时变且最终趋于零的概率突跳性，从而可有效避免陷入局部极小并最终趋于全局最优的串行结构的优化算法。蚁群算法蚁群算法又称蚂蚁算法，是一种用来在图中寻找优化路径的机率型算法。它由MarcoDorigo于1992年在他的博士论文中提出，其灵感来源于蚂蚁在寻找食物过程中发现路径的行为。蚁群算法是一种模拟进化算法，初步的研究表明该算法具有许多优良的性质。针对PID控制器参数优化设计问题，将蚁群算法设计的结果与遗传算法设计的结果进行了比较，数值仿真结果表明，蚁群算法具有一种新的模拟进化优化方法的有效性和应用价值。各个蚂蚁在没有事先告诉他们食物在什么地方的前提下开始寻找食物。当一只找到食物以后，它会向环境释放一种挥发性分泌物pheromone(称为信息素,该物质随着时间的推移会逐渐挥发消失，信息素浓度的大小表征路径的远近)来实现的，吸引其他的蚂蚁过来，这样越来越多的蚂蚁会找到食物。有些蚂蚁并没有像其它蚂蚁一样总重复同样的路，他们会另辟蹊径，如果另开辟的道路比原来的其他道路更短，那么，渐渐地，更多的蚂蚁被吸引到这条较短的路上来。最后，经过一段时间运行，可能会出现一条最短的路径被大多数蚂蚁重复着。原理：遗传算法下面是遗传算法的一般算法：

（1）创建一个随机的初始状态

初始种群是从解中随机选择出来的，将这些解比喻为染色体或基因，该种群被称为第一代，这和符号人工智能系统的情况不一样，在那里问题的初始状态已经给定了。

（2）评估适应度

对每一个解(染色体)指定一个适应度的值，根据问题求解的实际接近程度来指定(以便逼近求解问题的答案)。不要把这些“解”与问题的“答案”混为一谈，可以把它理解成为要得到答案，系统可能需要利用的那些特性。

（3）繁殖(包括子代突变)

带有较高适应度值的那些染色体更可能产生后代(后代产生后也将发生突变)。后代是父母的产物，他们由来自父母的基因结合而成，这个过程被称为“杂交”。

（4）下一代

如果新的一代包含一个解，能产生一个充分接近或等于期望答案的输出，那么问题就已经解决了。如果情况并非如此，新的一代将重复他们父母所进行的繁衍过程，一代一代演化下去，直到达到期望的解为止。

（5）并行计算

非常容易将遗传算法用到并行计算和群集环境中。一种方法是直接把每个节点当成一个并行的种群看待。然后有机体根据不同的繁殖方法从一个节点迁移到另一个节点。另一种方法是“农场主/劳工”体系结构，指定一个节点为“农场主”节点，负责选择有机体和分派适应度的值，另外的节点作为“劳工”节点，负责重新组合、变异和适应度函数的评估。蚁群算法范围：

蚂蚁观察到的范围是一个方格世界，蚂蚁有一个参数为速度半径（一般是3），那么它能观察到的范围就是3*3个方格世界，并且能移动的距离也在这个范围之内。

2、环境：

蚂蚁所在的环境是一个虚拟的世界，其中有障碍物，有别的蚂蚁，还有信息素，信息素有两种，一种是找到食物的蚂蚁洒下的食物信息素，一种是找到窝的蚂蚁洒下的窝的信息素。每个蚂蚁都仅仅能感知它范围内的环境信息。环境以一定的速率让信息素消失。

3、觅食规则：

在每只蚂蚁能感知的范围内寻找是否有食物，如果有就直接过去。否则看是否有信息素，并且比较在能感知的范围内哪一点的信息素最多，这样，它就朝信息素多的地方走，并且每只蚂蚁多会以小概率犯错误，从而并不是往信息素最多的点移动。蚂蚁找窝的规则和上面一样，只不过它对窝的信息素做出反应，而对食物信息素没反应。

4、移动规则：

每只蚂蚁都朝向信息素最多的方向移，并且，当周围没有信息素指引的时候，蚂蚁会按照自己原来运动的方向惯性的运动下去，并且，在运动的方向有一个随机的小的扰动。为了防止蚂蚁原地转圈，它会记住最近刚走过了哪些点，如果发现要走的下一点已经在最近走过了，它就会尽量避开。

5、避障规则：

如果蚂蚁要移动的方向有障碍物挡住，它会随机的选择另一个方向，并且有信息素指引的话，它会按照觅食的规则行为。

6、播撒信息素规则：

每只蚂蚁在刚找到食物或者窝的时候撒发的信息素最多，并随着它走远的距离，播撒的信息素越来越少。

根据这几条规则，蚂蚁之间并没有直接的关系，但是每只蚂蚁都和环境发生交互，而通过信息素这个纽带，实际上把各个蚂蚁之间关联起来了。比如，当一只蚂蚁找到了食物，它并没有直接告诉其它蚂蚁这儿有食物，而是向环境播撒信息素，当其它的蚂蚁经过它附近的时候，就会感觉到信息素的存在，进而根据信息素的指引找到了食物。

说了这么多，蚂蚁究竟是怎么找到食物的呢？

在没有蚂蚁找到食物的时候，环境没有有用的信息素，那么蚂蚁为什么会相对有效的找到食物呢？这要归功于蚂蚁的移动规则，尤其是在没有信息素时候的移动规则。首先，它要能尽量保持某种惯性，这样使得蚂蚁尽量向前方移动（开始，这个前方是随机固定的一个方向），而不是原地无谓的打转或者震动；其次，蚂蚁要有一定的随机性，虽然有了固定的方向，但它也不能像粒子一样直线运动下去，而是有一个随机的干扰。这样就使得蚂蚁运动起来具有了一定的目的性，尽量保持原来的方向，但又有新的试探，尤其当碰到障碍物的时候它会立即改变方向，这可以看成一种选择的过程，也就是环境的障碍物让蚂蚁的某个方向正确，而其他方向则不对。这就解释了为什么单个蚂蚁在复杂的诸如迷宫的地图中仍然能找到隐蔽得很好的食物。

当然，在有一只蚂蚁找到了食物的时候，其他蚂蚁会沿着信息素很快找到食物的。

蚂蚁如何找到最短路径的？

这一是要归功于信息素，另外要归功于环境，具体说是计算机时钟。信息素多的地方显然经过这里的蚂蚁会多，因而会有更多的蚂蚁聚集过来。假设有两条路从窝通向食物，开始的时候，走这两条路的蚂蚁数量同样多（或者较长的路上蚂蚁多，这也无关紧要）。当蚂蚁沿着一条路到达终点以后会马上返回来，这样，短的路蚂蚁来回一次的时间就短，这也意味着重复的频率就快，因而在单位时间里走过的蚂蚁数目就多，洒下的信息素自然也会多，自然会有更多的蚂蚁被吸引过来，从而洒下更多的信息素……；而长的路正相反，因此，越来越多地蚂蚁聚集到较短的路径上来，最短的路径就近似找到了。也许有人会问局部最短路径和全局最短路的问题，实际上蚂蚁逐渐接近全局最短路的，为什么呢？这源于蚂蚁会犯错误，也就是它会按照一定的概率不往信息素高的地方走而另辟蹊径，这可以理解为一种创新，这种创新如果能缩短路途，那么根据刚才叙述的原理，更多的蚂蚁会被吸引过来。模拟退火模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据Metropolis准则，粒子在温度T时趋于平衡的概率为e(-ΔE/(kT))，其中E为温度T时的内能，ΔE为其改变量，k为Boltzmann常数。用固体退火模拟组合优化问题，将内能E模拟为目标函数值f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(CoolingSchedule)控制，包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S。模拟退火算法的模型1模拟退火算法可以分解为解空间、目标函数和初始解三部分。2模拟退火的基本思想:(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步。程序：遗传算法：clearclccloseall%%加载数据%loaddataX=[13042312;36391315;41772244;37211399;34881535;33261556;32381229;41961004;4312790;4386570;30071970;25621756;27881491;23811676;1332695;37151678;39182179;40612370;37802212;36762578;40292838;42632931;34291908;35072367;33942643;34393201;29353240;31403550;25452357;27782826;23702975];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;holdon;boxonxlim([0,MAXGEN])title('优化过程')xlabel('代数')ylabel('最优值')ObjV=PathLength(D,Chrom);%计算路径长度preObjV=min(ObjV);whilegen<MAXGEN%%计算适应度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两两城市之间的距离functionD=Distanse(a)row=size(a,1);D=zeros(row,row);fori=1:rowforj=i+1:rowD(i,j)=((a(i,1)-a(j,1))^2+(a(i,2)-a(j,2))^2)^0.5;D(j,i)=D(i,j);endendDrawPath.m%%画路径函数%输入%Chrom待画路径%X各城市坐标位置functionDrawPath(Chrom,X)R=[Chrom(1,:)Chrom(1,1)];%一个随机解(个体)figure;holdonplot(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)fori=1:size(X,1)text(X(i,1)+0.05,X(i,2)+0.05,num2str(i),'color',[1,0,0]);endA=X(R,:);row=size(A,1);fori=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]);endholdoffxlabel('横坐标')ylabel('纵坐标')title('轨迹图')boxondsxy2figxy.mfunctionvarargout=dsxy2figxy(varargin)iflength(varargin{1})==1&&ishandle(varargin{1})...&&strcmp(get(varargin{1},'type'),'axes')hAx=varargin{1};varargin=varargin(2:end);elsehAx=gca;end;iflength(varargin)==1pos=varargin{1};else[x,y]=deal(varargin{:});endaxun=get(hAx,'Units');set(hAx,'Units','normalized');axpos=get(hAx,'Position');axlim=axis(hAx);axwidth=diff(axlim(1:2));axheight=diff(axlim(3:4));ifexist('x','var')varargout{1}=(x-axlim(1))*axpos(3)/axwidth+axpos(1);varargout{2}=(y-axlim(3))*axpos(4)/axheight+axpos(2);elsepos(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;endset(hAx,'Units',axun)Fitness.mfunctionFitnV=Fitness(len)FitnV=1./len;InitPop.m%%初始化种群%输入：%NIND：种群大小%N：城市的个数%输出：%初始种群functionChrom=InitPop(NIND,N)Chrom=zeros(NIND,N);%用于存储种群fori=1:NINDChrom(i,:)=randperm(N);%随机生成初始种群EndMutate.m%%变异操作%输入：%SelCh被选择的个体%Pm变异概率%输出：%SelCh变异后的个体functionSelCh=Mutate(SelCh,Pm)[NSel,L]=size(SelCh);fori=1:NSelifPm>=randR=randperm(L);SelCh(i,R(1:2))=SelCh(i,R(2:-1:1));endendOutputPat.m%%输出路径函数%输入：R路径functionp=OutputPath(R)R=[R,R(1)];N=length(R);p=num2str(R(1));fori=2:Np=[p,'—>',num2str(R(i))];enddisp(p)PathLength.m%%计算各个体的路径长度%输入：%D两两城市之间的距离%Chrom个体的轨迹functionlen=PathLength(D,Chrom)[row,col]=size(D);NIND=size(Chrom,1);len=zeros(NIND,1);fori=1:NINDp=[Chrom(i,:)Chrom(i,1)];i1=p(1:end-1);i2=p(2:end);len(i,1)=sum(D((i1-1)*col+i2));endRecombin.m%%交叉操作%输入%SelCh被选择的个体%Pc交叉概率%输出：%SelCh交叉后的个体functionSelCh=Recombin(SelCh,Pc)NSel=size(SelCh,1);fori=1:2:NSel-mod(NSel,2)ifPc>=rand%交叉概率Pc[SelCh(i,:),SelCh(i+1,:)]=intercross(SelCh(i,:),SelCh(i+1,:));endend%输入：%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]);ifr1~=r2a0=a;b0=b;s=min([r1,r2]);e=max([r1,r2]);fori=s:ea1=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);endif~isempty(i2)b(i2)=b1(i);endendendReins.m%%重插入子代的新种群%Chrom父代的种群%SelCh子代种群%ObjV父代适应度%Chrom组合父代与子代后得到的新种群functionChrom=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进化逆转后的个体functionSelCh=Reverse(SelCh,D)[row,col]=size(SelCh);ObjV=PathLength(D,SelCh);%计算路径长度SelCh1=SelCh;fori=1:rowr1=randsrc(1,1,[1:col]);r2=randsrc(1,1,[1:col]);mininverse=min([r1r2]);maxinverse=max([r1r2]);SelCh1(i,mininverse:maxinverse)=SelCh1(i,maxinverse:-1:mininverse);endObjV1=PathLength(D,SelCh1);%计算路径长度index=ObjV1<ObjV;SelCh(index,:)=SelCh1(index,:);Select.mfunctionSelCh=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被选择个体的索引号functionNewChrIx=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);程序运行结果为:初始轨迹：优化过程：优化后轨迹：模拟退火:clearclca=0.99; %温度衰减函数的参数t0=97;tf=3;t=t0;Markov_length=10000; %Markov链长度citys=[113042312;236391315;341772244;437211399;534881535;633261556;732381229;841961004;94312790;104386570;1130071970;1225621756;1327881491;1423811676;151332695;1637151678;1739182179;1840612370;1937802212;2036762578;2140292838;2242632931;2334291908;2435072367;2533942643;2634393201;2729353240;2831403550;2925452357;3027782826;3123702975];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;whilet>=tfforr=1:Markov_length %Markov链长度%产生随机扰动if(rand<0.5) %随机决定是进行两交换还是三交换%两交换ind1=0;ind2=0;while(ind1==ind2)ind1=ceil(rand.*amount);ind2=ceil(rand.*amount);endtmp1=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);endtmp1=ind1;tmp2=ind2;tmp3=ind3;%确保ind1<ind2<ind3if(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;endtmplist1=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;fori=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));ifE_new<E_currentE_current=E_new;sol_current=sol_new;ifE_new<E_best%把冷却过程中最好的解保存下来E_best=E_new;sol_best=sol_new;endelse%若新解的目标函数值小于当前解的，%则仅以一定概率接受新解ifrand<exp(-(E_new-E_current)./t)E_current=E_new;sol_current=sol_new;elsesol_new=sol_current;endendendt=t.*a; %控制参数t（温度）减少为原来的a倍enddisp('最优解为：')disp(sol_best)disp('最短距离：')disp(E_best)程序运行结果为:蚁群算法:clearallclc%%导入数据citys=[13042312;36391315;41772244;37211399;34881535;33261556;32381229;41961004;4312790;4386570;30071970;25621756;27881491;23811676;1332695;37151678;39182179;40612370;37802212;36762578;40292838;42632931;34291908;35072367;33942643;34393201;29353240;31403550;25452357;27782826;23702975];%%计算城市间相互距离n=size(citys,1);D=zeros(n,n);fori=1:nforj=1:nifi~=jD(i,j)=sqrt(sum((citys(i,:)-citys(j,:)).^2));elseD(i,j)=1e-4;endendend%%初始化参数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);%各代路径的平均长度%%迭代寻找最佳路径whileiter<=iter_max%随机产生各个蚂蚁的起点城市start=zeros(m,1);fori=1:mtemp=randperm(n);start(i)=temp(1);endTable(:,1)=start;%构建解空间citys_index=1:n;%逐个蚂蚁路径选择fori=1:m%逐个城市路径选择forj=2:ntabu=Table(i,1:(j-1));%已访问的城市集合(禁忌表)allow_index=~ismember(citys_index,tabu);allow=citys_index(allow_index);%待访问的城市集合P=allow;%计算城市间转移概率fork=1:length(allow)P(k)=Tau(tabu(end),allow(k))^alpha*Eta(tabu(end),allow(k))^beta;endP=P/sum(P);%轮盘赌法选择下一个访问城市Pc=cumsum(P);target_index=find(Pc>=rand);target=allow(target_index(1));Table(i,j)=target;endend%计算各个蚂蚁的路径距离Length=zeros(m,1);fori=1:mRoute=Table(i,:);forj=1:(n-1)Length(i)=Length(i)+D(Route(j),Route(j+1));endLength(i)=Length(i)+D(Route(n),Route(1));end%计算最短路径距离及平均距离ifiter==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);ifLength_best(iter)==min_LengthRoute_best(iter,:)=Table(min_index,:);elseRoute_best(iter,:)=Route_best((iter-1),:);endend%更新信息素Delta_Tau=zeros(n,n);%逐个蚂蚁计算fori=1:m%逐个城市计算fo