蚁群算法python代码_蚁群算法代码实现

急求蚁群算法解决 VRPTW问题的matlab代码，是ACS或者MMAS的！

function [R_best,L_best,L_e,Shortest_Route,Shortest_Length]=ACATSP(C,NC_max,m,Alpha,Beta,Rho,Q)

蚁群算法python代码_蚁群算法代码实现

%%=========================================================================

%% ACATSP.m

%% Ant Colony Algorithm for Treling Salean Problem

%% ChengAihua,PLA Information Engineering University,ZhengZhou,China

%% All rights reserved

%%-------------------------------------------------------------------------

%% 主要符号说明

%% C n个城市的坐标，n×2的矩阵

%% NC_max 迭代次数

%% m 蚂蚁个数

%% Alpha 表征信息素重要程度的参数

%% Beta 表征启发式因子重要程度的参数

%% Rho 信息素蒸发系数

%% Q 信息素增加强度系数

%% R_best 各代佳路线

%% L_best 各代佳路线的长度

%% 运行可能要很久，需要耐心等待

%%=========================================================================

n=length(C); %n 为市个数

for i=1:n %坐标矩阵转换为距离矩阵

for j=1:n

D(i,j)=sqrt((x(i,1)-x(j,1))^2+(x(i,2)-x(j,2))^2);

end

end

for i=1:n %Eta为启发因子，这里设为距离的倒数

for j=1:n %原文作者少考虑的当D=0是MATLAB提示出错

if i~=j

Eta(i,j)=1./D(i,j);

end

end

end

for i=1:n

Eta(i,i)=0;

end

Tau=ones(n,n); %Tau为信息素矩阵

Tabu=zeros(m,n); %存储并记录路径的生成

NC=1; %迭代计数器

R_best=zeros(NC_max,n); %各代佳路线

L_best=inf.ones(NC_max,1); %各代佳路线的长度

L_e=zeros(NC_max,1); %各代路线的平均长度

while NC<=NC_max %停止条件之一：达到迭代次数

%%第二步：将m只蚂蚁放到n个城市上

Randpos=[];

for i=1:(ceil(m/n))

Randpos=[Randpos,randperm(n)];

end

Tabu(:,1)=(Randpos(1,1:m))';

%%第三步：m只蚂蚁按概率函数选择下一座城市，完成各自的周游

for j=2:n

for i=1:m

visited=Tabu(i,1:(j-1)); %已访问的城市

J=zeros(1,(n-j+1)); %待访问的城市

P=J; %待访问城市的选择概率分布

Jc=1;

for k=1:n

if length(find(visited==k))==0

J(Jc)=k;

Jc=Jc+1;

end

end

%下面计算待选城市的概率分布

for k=1:length(J)

P(k)=(Tau(visited(end),J(k))^Alpha)(Eta(visited(end),J(k))^Beta);

end

P=P/(sum(P));

%按概率原则选取下一个城市

Pcum=cumsum(P);

Select=find(Pcum>=rand);

to_visit=J(Select(1));

Tabu(i,j)=to_visit;

end

end

if NC>=2

Tabu(1,:)=R_best(NC-1,:);

end

%%第四步：记录本次迭代佳路线

L=zeros(m,1);

for i=1:m

R=Tabu(i,:);

for j=1:(n-1)

L(i)=L(i)+D(R(j),R(j+1));

end

L(i)=L(i)+D(R(1),R(n));

end

L_best(NC)=min(L);

pos=find(L==L_best(NC));

R_best(NC,:)=Tabu(pos(1),:);

L_e(NC)=mean(L);

NC=NC+1;

%%第五步：更新信息素

Delta_Tau=zeros(n,n);

for i=1:m

for j=1:(n-1)

Delta_Tau(Tabu(i,j),Tabu(i,j+1))=Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q/L(i);

end

Delta_Tau(Tabu(i,n),Tabu(i,1))=Delta_Tau(Tabu(i,n),Tabu(i,1))+Q/L(i);

end

Tau=(1-Rho).Tau+Delta_Tau;

%%第六步：禁忌表清零

Tabu=zeros(m,n);

end

%%第七步：输出结果

Pos=find(L_best==min(L_best));

Shortest_Route=R_best(Pos(1),:);

Shortest_Length=L_best(Pos(1));

DrawRoute(C,Shortest_Route) %调用函数绘图

自适应蚁群算法在MATLAB中怎么实现？我在网上看了些代码，不知道怎么定义函数的，求算法的具体步骤

首先确定xyz的范围，比如说xyz都是从1到10，那么

[x,y,z]=meshgrid(1:10);

这时x，y，z都是3维矩阵，因此后面无法进行乘法运算是吧。因此函数没法写。你是这个地方卡住了是吧？？？

可以这样解决：

x=x(:);y=y(:);z=z(:);这样将xyz变成向量。就可以像平时一样定义函数了。

8 Metaheuristics

下图介绍了两种不同种类的Metaheuristics，我们主要用左边的，尤其是Iterated Greedy。

I 和D 的过程用图像表示如下：

我在这里只用Iterated Greedy算法。原因如下：

其他算法诸如蚁群算法等，可能能提供非常接近解的方案，有些算法的运行速度也很快，可以弥补python运行慢的缺陷。但是这些算法通常存在诸多不足，例如算法太复杂，适用面很窄，需要设置过多参数导致很难实现。

Iterated Greedy的优势在于，它由两个简单的阶段构成：

D和I

更好的解决方案总会被接受

更坏的方案以特定的可能性接受，接受的概率如下图

类似 模拟退火法 ：

从当前解决方案中随机删除 numberJobsToRemove 个订单。这里 numberJobsToRemove 是Iterated Greedy的一个参数。

输出的结果是被移除的订单 removedJobs 和一个不完整的解决方案 partialPermutation 。

注意：solver.RNG.choice的结果每次都一样，是因为我们设置了随机数种子，目的就是只要是用同一个种子作为参数构造出的solver的属性RNG都是同一个。

将 removedJobs 重新加回 partialPermutation ，并插入到佳位置（NEH）的顺序，并返回新的完整解决方案。这个插入过程是通过排列Permutation实现的，看一下Construction函数的参数表可知，需要两个列表。

注：前面讲算法的时候说过，重构函数执行之后会生成一个新的方案newSolution，新方案就是通过把removedJobs插入到佳位置得到的。位置的选择需要在Construction函数中借助 solver.EvaluationLogic.DetermineBestInsertion(completeSolution, i)来实现

我们通过简单的解构和重构，必然会得出一个新方案newSolution，那么我们接受新方案newSolution为当前方案currentSolution的前提是，如果

- 新的解决方案（ newSolution ），

- 当前解决方案（ currentSolution ）。

公式中，T叫做baseTemperature，是Iterated Greedy的第二个参数，T越大，则接受更坏方案的概率也越大。新的方案存储在SolutionPool中（numberJobsToRemove是第一个参数）

下面开始对newSolution进行评估和比较：

注意：判断一个新方案是否可接受取决于两方面：1. 如果新方案更好，那么无论如何都会接受新方案；2. 如果新方案并没有优化，则视作WorseSolution，即使是worseSolution也是要按照公式计算出的概率来衡量是否要接受这个不好的方案

为了看起来更加直观，我把接受方案的过程写成了函数AcceptWorseSolution。注意看参数表，以便于确定何时调用这个函数。

局部搜索可以被用于优化currentSOlution，但不是必须的。通常会使用IterativeImprovement结合Insertion邻域一起使用。

Iterated Greedy是一个迭代的解构D和构建C组成的序列。

在一个循环中被反复执行，直到达到 停止标准 。

在这里的例子中， 停止标准 是迭代次数 maxIterations ，但时间限制或没有优化的迭代次数也是可以的。

后面可能还会讲到怎么设计一个没有优化的迭代次数，这里可以先思考一下。

现在尝试运行一下：

与 IterativeImprovement 算法类似，现在要为 IteratedGreedy 创建一个单独的类，以便该类的一个实例可以传递给求解器Solver。

像 IterativeImprovement 一样， Iterated Greedy 应该继承自 ImprovementAlgorithm 。

必要的 参数 是以下属性：

EvaluationLogic, SolutionPool以及随机数生成器RNG都由求解器solver传递给算法。

添加 IteratedGreedy 类，以便它可以作为一种算法传递给求解器。

成员函数：Konstruktor, Initialize, Destruction, Construction, AcceptWorseSolution, Run

哪本python书立有蚁群算法

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

各个蚂蚁在没有事先告诉他们食物在什么地方的前提下开始寻找食物。当一只找到食物以后，它会向环境释放一种挥发性分泌物pheromone (称为信息素,该物质随着时间的推移会逐渐挥发消失，信息素浓度的大小表征路径的远近)来实现的，吸引其他的蚂蚁过来，这样越来越多的蚂蚁会找到食物。有些蚂蚁并没有像其它蚂蚁一样总重复同样的路，他们会另辟蹊径，如果另开辟的道路比原来的其他道路更短，那么，渐渐地，更多的蚂蚁被吸引到这条较短的路上来。后，经过一段时间运行，可能会出现一条短的路径被大多数蚂蚁重复着。

解决的问题

三维地形中，给出起点和重点，找到其路径。

作图源码：

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051

from mpl_toolkits.mplot3d import proj3dfrom mpl_toolkits.mplot3d import Axes3Dimport numpy as np height3d = np.array([[2000,1400,800,650,500,750,1000,950,900,800,700,900,1100,1050,1000,1150,1300,1250,1200,1350,1500], [1100,900,700,625,550,825,1100,1150,1200,925,650,750,850,950,1050,1175,1300,1350,1400,1425,1450], [200,400,600,600,600,900,1200,1350,1500,1050,600,600,600,850,1100,1200,1300,1450,1600,1500,1400], [450,500,550,575,600,725,850,875,900,750,600,600,600,725,850,900,950,1150,1350,1400,1450], [700,600,500,550,600,550,500,400,300,450,600,600,600,600,600,600,600,850,1100,1300,1500], [500,525,550,575,600,575,550,450,350,475,600,650,700,650,600,600,600,725,850,1150,1450], [300,450,600,600,600,600,600,500,400,500,600,700,800,700,600,600,600,600,600,1000,1400], [550,525,500,550,600,875,1150,900,650,725,800,700,600,875,1150,1175,1200,975,750,875,1000], [800,600,400,500,600,1150,1700,1300,900,950,1000,700,400,1050,1700,1750,1800,1350,900,750,600], [650,600,550,625,700,1175,1650,1275,900,1100,1300,1275,1250,1475,1700,1525,1350,1200,1050,950,850], [500,600,700,750,800,1200,1600,1250,900,1250,1600,1850,2100,1900,1700,1300,900,1050,1200,1150,1100], [400,375,350,600,850,1200,1550,1250,950,1225,1500,1750,2000,1950,1900,1475,1050,975,900,1175,1450], [300,150,0,450,900,1200,1500,1250,1000,1200,1400,1650,1900,2000,2100,1650,1200,900,600,1200,1800], [600,575,550,750,950,1275,1600,1450,1300,1300,1300,1525,1750,1625,1500,1450,1400,1125,850,1200,1550], [900,1000,1100,1050,1000,1350,1700,1650,1600,1400,1200,1400,1600,1250,900,1250,1600,1350,1100,1200,1300], [750,850,950,900,850,1000,1150,1175,1200,1300,1400,1325,1250,1125,1000,1150,1300,1075,850,975,1100], [600,700,800,750,700,650,600,700,800,1200,1600,1250,900,1000,1100,1050,1000,800,600,750,900], [750,775,800,725,650,700,750,775,800,1000,1200,1025,850,975,1100,950,800,900,1000,1050,1100], [900,850,800,700,600,750,900,850,800,800,800,800,800,950,1100,850,600,1000,1400,1350,1300], [750,800,850,850,850,850,850,825,800,750,700,775,850,1000,1150,875,600,925,1250,1100,950], [600,750,900,1000,1100,950,800,800,800,700,600,750,900,1050,1200,900,600,850,1100,850,600]]) fig = figure()ax = Axes3D(fig)X = np.arange(21)Y = np.arange(21)X, Y = np.meshgrid(X, Y)Z = -20np.exp(-0.2np.sqrt(np.sqrt(((X-10)2+(Y-10)2)/2)))+20+np.e-np.exp((np.cos(2np.piX)+np.sin(2np.piY))/2)ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='cool')ax.set_xlabel('X axis')ax.set_ylabel('Y axis')ax.set_zlabel('Z')ax.set_title('3D map') point0 = [0,9,Z[0][9]]point1 = [20,7,Z[20][7]] ax.plot([point0[0]],[point0[1]],[point0[2]],'r',marker = u'o',markersize = 15)ax.plot([point1[0]],[point1[1]],[point1[2]],'r',marker = u'o',markersize = 15) x0,y0,_ = proj3d.proj_transform(point0[0],point0[1],point0[2], ax.get_proj())x1,y1,_ = proj3d.proj_transform(point1[0],point1[1],point1[2], ax.get_proj()) label = pylab.annotate( "start", xy = (x0, y0), xytext = (-20, 20), textcoords = 'offset points', ha = 'right', va = 'bottom', bbox = dict(boxstyle = 'round,pad=0.5', fc = 'yellow', alpha = 1), arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'),fontsize=15)label2 = pylab.annotate( "end", xy = (x1, y1), xytext = (-20, 20), textcoords = 'offset points', ha = 'right', va = 'bottom', bbox = dict(boxstyle = 'round,pad=0.5', fc = 'yellow', alpha = 1), arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'),fontsize=15)def update_position(e): x2, y2, _ = proj3d.proj_transform(point0[0],point0[1],point0[2],ax.get_proj()) label.xy = x2,y2 label.update_positions(fig.canvas.renderer) x1,y1,_ = proj3d.proj_transform(point1[0],point1[1],point1[2],ax.get_proj()) label2.xy = x1,y1 label2.update_positions(fig.canvas.renderer) fig.canvas.draw() fig.canvas.mpl_connect('button_release_event', update_position)

基本原理

蚂蚁k根据各个城市间链接路径上的信息素浓度决定其下一个访问城市，设Pkij(t)表示t时刻蚂蚁k从城市i转移到矩阵j的概率，其计算公式为

计算完城市间的转移概率后，采用与遗传算法中一样的方法选择下一个待访问的城市。

当所有的蚂蚁完成一次循环后，各个城市间链接路径上的信息素浓度需进行更新，计算公式为

其中，Δτkij表示第k只蚂蚁在城市i与城市j连接路径上释放的信息素浓度；Δτij表示所有蚂蚁在城市i与城市j连接路径上释放的信息素浓度之和。

蚂蚁释放信息素的模型

程序代码：

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117

import numpy as npimport matplotlib.pyplot as plt%pylabcoordinates = np.array([[565.0,575.0],[25.0,185.0],[345.0,750.0],[945.0,685.0],[845.0,655.0], [880.0,660.0],[25.0,230.0],[525.0,1000.0],[580.0,1175.0],[650.0,1130.0], [1605.0,620.0],[1220.0,580.0],[1465.0,200.0],[1530.0, 5.0],[845.0,680.0], [725.0,370.0],[145.0,665.0],[415.0,635.0],[510.0,875.0],[560.0,365.0], [300.0,465.0],[520.0,585.0],[480.0,415.0],[835.0,625.0],[975.0,580.0], [1215.0,245.0],[1320.0,315.0],[1250.0,400.0],[660.0,180.0],[410.0,250.0], [420.0,555.0],[575.0,665.0],[1150.0,1160.0],[700.0,580.0],[685.0,595.0], [685.0,610.0],[770.0,610.0],[795.0,645.0],[720.0,635.0],[760.0,650.0], [475.0,960.0],[95.0,260.0],[875.0,920.0],[700.0,500.0],[555.0,815.0], [830.0,485.0],[1170.0, 65.0],[830.0,610.0],[605.0,625.0],[595.0,360.0], [1340.0,725.0],[1740.0,245.0]])def getdistmat(coordinates): num = coordinates.shape[0] distmat = np.zeros((52,52)) for i in range(num): for j in range(i,num): distmat[i][j] = distmat[j][i]=np.linalg.norm(coordinates[i]-coordinates[j]) return distmatdistmat = getdistmat(coordinates)numant = 40 #蚂蚁个数numcity = coordinates.shape[0] #城市个数alpha = 1 #信息素重要程度因子beta = 5 #启发函数重要程度因子rho = 0.1 #信息素的挥发速度Q = 1iter = 0itermax = 250etatable = 1.0/(distmat+np.diag([1e10]numcity)) #启发函数矩阵，表示蚂蚁从城市i转移到矩阵j的期望程度pheromonetable = np.ones((numcity,numcity)) # 信息素矩阵pathtable = np.zeros((numant,numcity)).astype(int) #路径记录表distmat = getdistmat(coordinates) #城市的距离矩阵lengther = np.zeros(itermax) #各代路径的平均长度lengthbest = np.zeros(itermax) #各代及其之前遇到的佳路径长度pathbest = np.zeros((itermax,numcity)) # 各代及其之前遇到的佳路径长度while iter < itermax: # 随机产生各个蚂蚁的起点城市 if numant <= numcity:#城市数比蚂蚁数多 pathtable[:,0] = np.random.permutation(range(0,numcity))[:numant] else: #蚂蚁数比城市数多，需要补足 pathtable[:numcity,0] = np.random.permutation(range(0,numcity))[:] pathtable[numcity:,0] = np.random.permutation(range(0,numcity))[:numant-numcity] length = np.zeros(numant) #计算各个蚂蚁的路径距离 for i in range(numant): visiting = pathtable[i,0] # 当前所在的城市 #visited = set() #已访问过的城市，防止重复 #visited.add(visiting) #增加元素 unvisited = set(range(numcity))#未访问的城市 unvisited.remove(visiting) #删除元素 for j in range(1,numcity):#循环numcity-1次，访问剩余的numcity-1个城市 #每次用轮盘法选择下一个要访问的城市 listunvisited = list(unvisited) probtrans = np.zeros(len(listunvisited)) for k in range(len(listunvisited)): probtrans[k] = np.power(pheromonetable[visiting][listunvisited[k]],alpha)\ np.power(etatable[visiting][listunvisited[k]],alpha) cumsumprobtrans = (probtrans/sum(probtrans)).cumsum() cumsumprobtrans -= np.random.rand() k = listunvisited[find(cumsumprobtrans>0)[0]] #下一个要访问的城市 pathtable[i,j] = k unvisited.remove(k) #visited.add(k) length[i] += distmat[visiting][k] visiting = k length[i] += distmat[visiting][pathtable[i,0]] #蚂蚁的路径距离包括后一个城市和第一个城市的距离 #print length # 包含所有蚂蚁的一个迭代结束后，统计本次迭代的若干统计参数 lengther[iter] = length.mean() if iter == 0: lengthbest[iter] = length.min() pathbest[iter] = pathtable[length.argmin()].copy() else: if length.min() > lengthbest[iter-1]: lengthbest[iter] = lengthbest[iter-1] pathbest[iter] = pathbest[iter-1].copy() else: lengthbest[iter] = length.min() pathbest[iter] = pathtable[length.argmin()].copy() # 更新信息素 changepheromonetable = np.zeros((numcity,numcity)) for i in range(numant): for j in range(numcity-1): changepheromonetable[pathtable[i,j]][pathtable[i,j+1]] += Q/distmat[pathtable[i,j]][pathtable[i,j+1]] changepheromonetable[pathtable[i,j+1]][pathtable[i,0]] += Q/distmat[pathtable[i,j+1]][pathtable[i,0]] pheromonetable = (1-rho)pheromonetable + changepheromonetable iter += 1 #迭代次数指示器+1 #观察程序执行进度，该功能是非必须的 if (iter-1)%20==0: print iter-1# 做出平均路径长度和路径长度 fig,axes = plt.subplots(nrows=2,ncols=1,figsize=(12,10))axes[0].plot(lengther,'k',marker = u'')axes[0].set_title('Average Length')axes[0].set_xlabel(u'iteration')axes[1].plot(lengthbest,'k',marker = u'')axes[1].set_title('Best Length')axes[1].set_xlabel(u'iteration')fig.sefig('Average_Best.png',dpi=500,bbox_inches='tight')plt.close()#作出找到的路径图bestpath = pathbest[-1]plt.plot(coordinates[:,0],coordinates[:,1],'r.',marker=u'$\cdot$')plt.xlim([-100,2000])plt.ylim([-100,1500])for i in range(numcity-1):# m,n = bestpath[i],bestpath[i+1] print m,n plt.plot([coordinates[m][0],coordinates[n][0]],[coordinates[m][1],coordinates[n][1]],'k')plt.plot([coordinates[bestpath[0]][0],coordinates[n][0]],[coordinates[bestpath[0]][1],coordinates[n][1]],'b')ax=plt.gca()ax.set_title("Best Path")ax.set_xlabel('X axis')ax.set_ylabel('Y_axis')plt.sefig('Best Path.png',dpi=500,bbox_inches='tight')plt.close()