一、简介
1 蚁群算法的提出
蚁群算法(ant colony optimization, ACO)，又称蚂蚁算法，是一种用来寻找优化路径的机率型算法。它由Marco Dorigo于1992年在他的博士论文中提出，其灵感来源于蚂蚁在寻找食物过程中发现路径的行为。遗传算法在模式识别、神经网络、机器学习、工业优化控制、自适应控制、生物科学、社会科学等方面都得到应用。
2 算法的基本原理
二、源代码
%%%%%%%%%% 问题定义：获取城市位置坐标、计算距离矩阵 %%%%%%%%%
InitOps=[];
[Location,DistMatrix,Ncities,Bestx,Lengx] = pr76init(InitOps);
close all;
figure (1);
hold on;
minx=min(DistMatrix(:,1));
maxx=max(DistMatrix(:,1));
miny=min(DistMatrix(:,2));
maxy=max(DistMatrix(:,2));
minm=min(minx,miny);
maxm=max(maxx,maxy);
l=(maxm-minm)/10;
for i=1:Ncities
plot(Location(i,1),Location(i,2),'*b');
text (Location(i,1)+l,Location(i,2)+l,num2str(i));
end
for i=1:Ncities-1
line([Location(Bestx(i),1),Location(Bestx(i+1),1)] , [Location(Bestx(i),2),Location(Bestx(i+1),2)]) ;
end
line([Location(Bestx(1),1),Location(Bestx(Ncities),1)] , [Location(Bestx(1),2),Location(Bestx(Ncities),2)]) ;
grid on,title(['初始路线图-',num2str(Lengx)]),xlabel('横坐标'),ylabel('纵坐标');
legend('城市位置');
hold off ;
% 初始化随机发生器状态
rand('state',sum(100*clock));
% ================================================
% 使用最近邻法构造一个初始游历,并据此计算信息系初值
p=zeros(1,Ncities+1);
p(1)=round(Ncities*rand+0.5);% p存储目前找到的所有城市的编号
i=p(1);
count=2;
while count <= Ncities
NNdist= inf ;%NNdist存储目前找到的和当前城市距离最短的城市的距离
pp= i ;% i存储当前城市的编号 pp存储目前找到的城市编号
for j= 1: Ncities
if (DistMatrix(i, j) < NNdist) & (j~=i) & ((j~=p) == ones(1,length(p)))
% 目标城市的要求为－－距离短、且不能是当前城市，也不能是以前已经走过的城市
NNdist= DistMatrix(i, j) ;
pp= j ;
end
end
p(count)=pp;
i= pp ;
count= count + 1 ;
end
p=AcatspEval(p,DistMatrix,Ncities);
len=p(1,Ncities+1);
Q0=1/(Ncities*len);
%%%%%%%%%% 设定系统有关参数 %%%%%%%%%%
MaxNc=5000;% 最大代数
A=1;% 信息素因子
B=2;% 启发信息因子
P1=0.1;% 局部挥发系数初值
P2=0.1;% 全局挥发系数初值
R0=0.9; %选择概率
M=10;% 蚂蚁数量
%%%%%%%%%% 初始化信息素、启发信息矩阵、确定蚂蚁最初位置及允许矩阵 %%%%%%%%%%
Pheromone=Q0*ones(Ncities,Ncities);% 信息素初始矩阵;
Heuristic=1./DistMatrix;% 启发信息初始矩阵
Temp=ones(1,Ncities);
Heuristic=1./(1./Heuristic+diag(Temp));
RandL=round(rand(M,1)*Ncities+0.5);%蚂蚁最初位置
Alocation0=zeros(M,Ncities+1);% 存放M+1个蚂蚁游历的路径及长度矩阵初始化
Alocation0(:,1)=RandL;
Allow0=repmat(1:Ncities,M,1);% 允许访问的城市矩阵初始化
for Ak=1:M
Allow0(Ak,RandL(Ak))=0;
end
%%%%%%%%%% 运行参数初始化 %%%%%%%%%%
Nc=1;% 第一代
Lbestdis=inf;
Cbestdis=inf;
Fnewbest=0;
Alocation=Alocation0;% 存放个蚂蚁游历的路径及长度矩阵初始化
Allow=Allow0; % 允许矩阵赋初值
t1=clock;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 蚁群算法初始化程序结束 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 蚁群算法主循环开始 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while(Nc<=MaxNc)
% M个蚂蚁选择Ncities个城市
for Cityi=2:Ncities+1
if Cityi<Ncities+1
for Ak=1:M
i=Alocation(Ak,Cityi-1);% 当前城市
j=Select_for_aca(R0,Ak,i,Allow,A,B,Pheromone,Heuristic);% 依据Pij选择下一个城市j
Alocation(Ak,Cityi)=j;
Allow(Ak,j)=0;% 更新允许矩阵
Pheromone(i,j)=(1-P1)*Pheromone(i,j)+P1*Q0; % 信息素在线单步更新
Pheromone(j,i)= Pheromone(i,j);
end
else % 返回出发城市
for Ak=1:M
i=Alocation(Ak,Cityi-1);% 当前城市
j=Alocation(Ak,1);
Pheromone(i,j)=(1-P1)*Pheromone(i,j)+P1*Q0; % 信息素在线单步更新
Pheromone(j,i)=Pheromone(i,j);
end
end
end
function D = dists(X1,X2,p,e)
%DISTS Calculates distances between vectors of points.
% D = dists(X1,X2,p,e)
% X1 = n x d matrix of n d-dimensional points
% X2 = m x d matrix of m d-dimensional points
% D = n x m matrix of distances
% = (n-1) x 1 vector of distances between X1 points, if X2 = []
% p = 2, Euclidean (default): D(i,j) = sqrt(sum((X1(i,:) - X2(j,:))^2))
% = 1, rectilinear: D(i,j) = sum(abs(X1(i,:) - X2(j,:))
% = Inf, Chebychev dist: D(i,j) = max(abs(X1(i,:) - X2(j,:))
% = (1 Inf), lp norm: D(i,j) = sum(abs(X1(i,:) - X2(j,:))^p)^(1/p)
% = 'rad', great circle distance in radians of a sphere
% (where X1 and X2 are decimal degree latitudes and longitudes)
% = 'mi' or 'sm', great circle distance in statute miles on the earth
% = 'km', great circle distance in kilometers on the earth
% e = epsilon for hyperboloid approximation gradient estimation
% = 0 (default); no error checking if any non-empty 'e' input
% ~= 0 => general lp used for rect., Cheb., and p outside [1,2]
%
% Great circle distances are calculated using the Haversine Formula (from R.W.
% Sinnott, "Virtues of the Haversine", Sky and Telescope, vol. 68, no. 2, 1984
% p. 159)
% Copyright (c) 1998 by Michael G. Kay
% Matlog Version 1.0 Apr-3-98
% Input Error Checking *******************************************************
if nargin == 4 & ~isempty(e) % No error checking is 'e' input
[n,d] = size(X1);
m = size(X2,1);
else % Do error checking
error(nargchk(1,4,nargin));
e = 0; % 'e' default
[n,d] = size(X1);
if n == 0 | ~isreal(X1)
error('X1 must be non-empty real matrix');
end
if nargin < 2 | isempty(X2) % Calc intra-seq dist b/w X1 points
m = 0; % X2 default
if n < 2
error('X1 must have more than 1 point');
end
else % Calc dist b/w X1 and X2 points
[m,dX2] = size(X2);
if m == 0 | ~isreal(X2)
error('X2 must be non-empty real matrix');
end
if d ~= dX2
error('Rows of X1 and X2 must have same dimensions');
end
end
if nargin < 3 | isempty(p)
p = 2; % 'p' default
elseif ~ischar(p) % lp distances
if length(p(:)) ~= 1 | ~isreal(p)
error('''p'' must be a real scalar number');
end
elseif ischar(p) % Great circle distances
p = lower(p);
if d ~= 2
error('Points must be 2-dimensional for great-circle distances');
end
if ~any(strcmp(p,{'rad','mi','sm','km'}))
error('''p'' must be either ''rad,'' ''mi,'' ''sm,'' or ''km''');
end
else
error('''p'' not valid value');
end
end
% End (Input Error Checking) ***********************************************
% Interchange if X2 is the only 1 point
intrchg = 0;
if n > 1 & m == 1
tmp = X2; X2 = X1; X1 = tmp;
m = n;n = 1;
intrchg = 1;
end
% 1-dimensional points
if d == 1
if e == 0
if m ~= 0
D = abs(X1(:,ones(1,m)) - X2(:,ones(1,n))');
else
D = abs(X1(1:n-1) - X1(2:n))'; % X1 intra-seq. dist.
end
else
if m ~= 0
D = sqrt((X1(:,ones(1,m)) - X2(:,ones(1,n))').^2 + e);
else
D = sqrt((X1(1:n-1) - X1(2:n)).^2 + e)';
end
end
% X1 only 1 point or intra-seq dist
elseif n == 1 | m == 0
if n == 1 % Expand X1 to match X2
X1 = X1(ones(1,m),:);
n = m;
else % X1 intra-seq. dist.
X2 = X1(2:n,:); % X2 = ending points
n = n - 1;
X1 = X1(1:n,:); % X1 = beginning points
end
if p == 2 % Euclidean distance
D = sqrt(sum(((X1 - X2).^2 + e)'));
elseif ischar(p) % Great-circle distance
X1 = pi*X1/180;X2 = pi*X2/180;
D = 2*asin(min(1,sqrt(sin((X1(:,1) - X2(:,1))/2).^2 + ...
cos(X1(:,1)).*cos(X2(:,1)).* ...
sin((X1(:,2) - X2(:,2))/2).^2)))';
elseif p == 1 & e == 0 % Rectilinear distance
D = sum(abs(X1 - X2)');
elseif (p >= 1 & p <= 2) | (e ~= 0 & p > 0) % General lp distance
D = sum((((X1 - X2).^2 + e).^(p/2))').^(1/p);
elseif p == Inf & e == 0 % Chebychev distance
D = max(abs(X1 - X2)');
else % Otherwise
D = zeros(1,n);
for j = 1:n
D(j) = norm(X1(j,:) - X2(j,:),p);
end
end
三、运行结果
四、备注
版本：2014a
