改进的蚁群算法-MATLAB代码

function [optimalway,Lmin]=ACO_Improved(m,route,...
alpha,beta,C,Q,rou,SP,EP,P,MAXCTimes,MAXFTimes,VP,ShowInterval)
%ACO_Improved 改进的蚁群算法 求解最优路径
%问题描述：
% 已知一交通路网，求一起点到终点的最短路径
%input:
% m -- 蚂蚁数量，标量
% route -- 连接城市间的路线信息，N*4的矩阵，N为连接两城市的路线条数
% 第一列：路线起点编号
% 第二列：路线终点编号
% 第三列：路线起点到终点的长度
% 第四列：路线起点到终点的权重
% 第五列：路线终点到起点的长度
% 第六列：路线终点到起点的权重
% alpha -- 信息激素启发因子，标量
% beta -- 自启发量（能见度）启发因子，标量
% C -- 初始各条路线上的信息激素浓度，标量
% Q -- 每只蚂蚁所携带的信息量，标量
% rou -- 每次循环，信息激素挥发系数，0到1的标量
% SP -- start point，起点，标量
% EP -- end point，终点，标量
% P -- 惩罚系数，0到1之间的标量
% MAXCTimes -- Max Cycle Times 最大循环次数，标量
% MAXFTimes -- Max Fail Times
% 最大失败次数，标量，当失败次数大于MAXFTime时，认为问题无解
% VP -- (visible probablity) 可见概率，标量。用于后面统计岔路数用。
% 假如某点处，蚂蚁选择路径1、2、3的概率分别为 10%,20%,30%,40%,VP=15%,
% 那么认为路径2、3、4是可见的，路径1是不可见的。
% ShowInterval -- 标量，每隔多少次循环(cycle)，显示下当前运算的最佳路径
%
%output:
% optimalway -- 从起点到终点的最优路径，行向量
% Lmin -- 所有循环中，最短路程
%
%参考文献：
% Dorigo M, Maniezzo Vittorio, Colorni Alberto.
% The Ant System: Optimization by a colony of cooperating agents [J].
% IEEE Transactions on Systems, Man, and Cybernetics--Part B,1996, 26(1)
% 《生命线地震工程》，大连理工大学，柳春光著

%variable define:
% n -- 城市数
% NodePoint -- 节点编号，1~n的行向量
% D -- n*n的矩阵，D(i,j)即为城市i和城市j的距离
% way -- n*n的矩阵，记录可行的通路。way(i,j)=1代表城市i可去往城市j，0代表不可
% tau -- n*n的矩阵，信息素密度
% tabu -- m*n的矩阵，禁忌表（各蚂蚁的访问记录）
% CycleCounter -- 循环次数计数器，标量
% PCounter -- Punish Counter 惩罚计数器，标量
% FCounter -- fail Counter 失败次数计数器，标量。
% 某循环中所有蚂蚁都未成功到达终点，则计一次失败，FCounter自加1
% flag -- m*1 的列向量，flag(k)=1代表第k只蚂蚁到达终点，0则代表未到达终点
% ，-1代表该蚂蚁是受惩罚的蚂蚁。
% ant -- 蚂蚁下标，代表第ant只蚂蚁，标量
% t -- 时间，标量，1~n
% NodeNow -- 某循环中蚂蚁ant在时间t所在节点位置,

标量
% NodeNext -- 某循环中蚂蚁ant在时间t+1所在节点位置,标量
% tabu_ant -- % 某循环中第ant只蚂蚁t时刻的禁忌向量
% PunishFlag -- 某循环中蚂蚁ant在时间t是否受惩罚，1代表受到惩罚，0代表没有
% CN -- connect node,当前与 NodeNow相连的节点，向量
% S2E -- m*n的矩阵。完成一个循环后，对tabu所记录的路径进行优化，弯道取直。
% 如果第ant只蚂蚁受到了惩罚，则 S2E 第ant行为0向量
% L -- m*1 的向量。计算 S2E 记录各蚂蚁从起点爬到终点的路程。
% 受惩罚的蚂蚁路程为inf
% LCmin -- 某次循环中，L中最小值
% LCmin_index -- 某次循环中，取得最短路径的蚂蚁的下标
% ANB -- (average net branch) 标量。S2E 中记录的，各条路径上节点岔路的平均值
% optimalwayC -- 某次循环中最优路径
% LminIndex -- 在第多少次循环（cycle）发现了最短路程
% LminRecord -- 每次迭代所发现的最优路径的路程记录，1*MAXITime 的行向量
% ANB_Record -- ANB=Average Node Branching 沿途每点的平均路径分岔数
% 1*MAXITime 的行向量
%
% ========================================
% 联系方式： matrixsuper@
% ========================================

rand('twister',1);

NodePoint(1,:)=sort(unique(route(:,[1,2])));
n=NodePoint(end);
if NodePoint==1:n %#ok
else error('检查节点编号');
end

Lmin=inf;
optimalway=zeros(1,n);
LminRecord=zeros(1,MAXCTimes);
ANB_Record=zeros(1,MAXCTimes);

[D,way]=distance_and_way(route,n);% 计算城市间距离和通路
tau=Init_tau(D,C,n);% 各条路径信息素密度初始化为C
CycleCounter=1;% 循环次数初始化为1
FCounter=0;% 失败次数初始化为0

fprintf('起点：%d,终点：%d,显示间隔：%d\n',SP,EP,ShowInterval);
disp('计算最优路径中：');

while true
tabu=InitTabu(SP,m,n);%初始禁忌表，将m只蚂蚁放置到起点
flag=InitFlag(m);%初始flag均为0（未到达终点）
Pcounter=0;% 惩罚计数器初始化为0
% m只蚂蚁开始迭代====================================
for t=1:n-1
for ant=1:m
NodeNow=tabu(ant,t);
tabu_ant=tabu(ant,:);

CN=ConnectNow(NodeNow,way);% 找出当前与 NodeNow相连的节点
[NodeNext,PunishFlag]=...% 求某循环中蚂蚁ant在%时间t+1所在位置
Next(ant,NodeNow,tabu_ant,tau,D,CN,alpha,beta,flag);
tabu(ant,t+1)=NodeNext;

if PunishFlag
Pcounter=Pcounter+1;
flag(ant)=-1;
tau=Punishment(tau,CN,tabu_ant,P);
else
if NodeNext==EP % 蚂蚁ant到达终点
flag(ant)=1;
end
end
end
end
% m只蚂蚁迭代结束====================================
if Pcounter==m % 所有蚂蚁都失败了

FCounter=FCounter+1;
if FCounter==MAXFTimes % 失败次数超过上限，认为问题无解
error('问题无解！');
else
continue; % 重新该次cycle
end
end
S2E=S2E_Compute(tabu,flag,way,m,n,EP);
L=L_Compute(S2E,D,EP,m,flag);
[LCmin,LCmin_index]=min(L);
LminRecord(CycleCounter)=LCmin;
optimalwayC=S2E(LCmin_index,:);
if LCminLmin=LCmin;
optimalway=optimalwayC;
LminIndex=CycleCounter;
end
ANB=ANB_Compute(S2E,tau,D,VP,flag,m,EP,alpha,beta,n);
ANB_Record(CycleCounter)=ANB;
tau=tau_refresh(tau,S2E,L,Q,rou,EP,flag,m);
ShowOptimalWay(CycleCounter,ShowInterval,optimalwayC,EP);
if CycleCounter==MAXCTimes
break;% 完成所有循环，跳出
end
CycleCounter=CycleCounter+1;
end

disp('==================================================================');
disp(['第',int2str(LminIndex),'次循环发现最优路径：']);
optimalway=optimalway(1:find(optimalway==EP,1,'first'));
disp(optimalway);
disp(['最优路径的长度：',num2str(Lmin)]);
figure(1);
subplot(2,1,1);
hold on;
title({['起点：',int2str(SP),',终点：',int2str(EP)],...
['最优路径的长度：',num2str(Lmin)]});
plot(LminRecord);
xlabel('循环次数');
ylabel('最优路径长度');
subplot(2,1,2);
plot(ANB_Record);
xlabel('循环次数');
ylabel('平均岔路数');
tem=axis;
tem(3)=0;
axis(tem);

end

% =========================================================================
function ANB=ANB_Compute(S2E,tau,D,VP,flag,m,EP,alpha,beta,n)
% 计算平均岔路数
% variable define:
% ANB_Sum -- 岔路之和，标量
% P -- n*n的矩阵，P(i,j)=tau(i,j).^alpha./D(i,j).^beta
% r -- 行列标，标量
% hand -- 下标，标量
% term -- ANB_Sum每次自加的项
% counter -- ANB_Sum自加的项数
% position -- 第r只蚂蚁第hand步所在位置，标量
% BN -- 每个节点的分岔数

ANB_Sum=0;
counter=0;
P=tau.^alpha./D.^beta;
P(logical(eye(n)))=0;
for r=1:n
P(r,:)=P(r,:)./sum(P(r,:));
end

BN=zeros(n,1);
for r=1:n
BN(r)=length(find(P(r,:)>VP));
end

for r=1:m
if flag(r)==-1
continue;
end
hand=1;
while true
position=S2E(r,hand);
if position==EP
break;
end
term=BN(position);
ANB_Sum=ANB_Sum+term;
counter=counter+1;
hand=hand+1;
end
end
ANB=ANB_Sum/counter;
end

function CN=ConnectNow(NodeNow,way)
% 找出当前与 NodeNow相连的节点,CN为向量
% 比如
% way=[...
% 0 1 1
% 1 0 1
% 1 1 0
% ];
% NodeNow=2;
% 那么 CN=[1 3],waynow=[ 1 0 1 ];

waynow=way(NodeNow,:);
CN=find(waynow==1);
end

function [D,way]=distance_and_way(route,n)
% 计算城市间距离和通路
D=Inf.*ones(n,n);
way=zeros(n,n);
N=size(route,1);
for k=1:N
D(route(k,1),route(k,2))=route(k,3)*route(k,4);
D(route(k,

2),route(k,1))=route(k,5)*route(k,6);
end
for k=1:n
D(k,k)=0;
end
way(D~=inf)=1;
way(D==0)=0;
end

function tau=Init_tau(D,C,n)
% 各条路径信息素密度初始化为C
tau=zeros(n,n);
for ii=1:n
for jj=1:n
if D(ii,jj)==inf || D(ii,jj)==0
tau(ii,jj)=0;
else
tau(ii,jj)=C;
end
end
end
end

function flag=InitFlag(m)
%初始flag均为0（未到达终点）
flag=zeros(m,1);
end

function tabu=InitTabu(SP,m,n)
%初始禁忌表，将m只蚂蚁放置到起点
tabu=zeros(m,n);
for ii=1:m
tabu(ii,1)=SP;
end
end

function L=L_Compute(S2E,D,EP,m,flag)
%计算 S2E 中记录的各条路径的长度
%variable define:
% index -- 受惩罚蚂蚁的下标，标量
% r -- 行标，标量
% hand -- 标量
% Lsum -- 路程之和，标量

L=zeros(m,1);
index=(flag==-1);
L(index)=inf;
for r=1:m
if flag(r)==-1
continue;
end
hand=2;Lsum=0;
while true
if S2E(r,hand)==EP
Lsum=Lsum+D( S2E(r,hand-1),S2E(r,hand) );
break;
end
Lsum=Lsum+D( S2E(r,hand-1),S2E(r,hand) );
hand=hand+1;
end
L(r)=Lsum;
end
end

function [NodeNext,PunishFlag]=...
Next(ant,NodeNow,tabu_ant,tau,D,CN,alpha,beta,flag)
% NEXT
% 求某循环中蚂蚁ant在%时间t+1所在位置

% variable define:
% index -- 去除禁忌点后与NodeNow连接的节点编号向量
% taunow -- 向量，从NodeNow出发通往除禁忌点以外的点的各条路上的信息素密度
% Dnow -- 向量，从NodeNow出发通往除禁忌点以外的点的各条路的长度
% P -- 向量，从NodeNow出发通往除禁忌点以外的点的各条路的概率
% r -- 随机产生的一个0到1的数
% sumP -- 标量
% ind -- 标量，向量index的下标

flag_ant=flag(ant);
if flag_ant== 1 || flag_ant==-1 %已到达终点或已受到惩罚
PunishFlag=false;
NodeNext=NodeNow;%呆在终点或惩罚点不动
else
index=setdiff(CN,tabu_ant);

if isempty(index) %下一个可以选择的点只能是已经选过的点
PunishFlag=true;
NodeNext=NodeNow;
else % 产生随机数，确定下一个点
PunishFlag=false;
taunow=tau(NodeNow,index);
Dnow=D(NodeNow,index);
P=taunow.^alpha./Dnow.^beta;
P=P./sum(P);
r=rand;
sumP=0;
ind=1;
while true
sumP=sumP+P(ind);
if sumP>r
break;
else
ind=ind+1;
end
end
NodeNext=index(ind);
end
end

end

function tau_after=Punishment(tau_before,CN,tabu_ant,P)
% tau_before -- 惩罚前信息素密度矩阵
% tau_after -- 惩罚后信息素密度矩阵

% variable define:
% index -- 标量，用于找出tabu_ant 中第一个和NodeNow 相连的节点

tau_after=tau_before;
index=1;
while true
if isempty( find( CN==tabu_ant(index) ,

1 ) )
index=index+1;
else
break;
end
end
tau_after( tabu_ant(index),tabu_ant(index+1) )=...
tau_after( tabu_ant(index),tabu_ant(index+1) )*P;
end

function S2E=S2E_Compute(tabu,flag,way,m,n,EP)
% S2E （start to end ）对tabu表中记录的路径进行弯道取直的优化
% variable define:
% r -- 行标
% lh -- left hand，标量
% rh -- right hand，标量
% hand -- 标量
% ri -- record index,标量,S2E 第r行的记录位置
% Nl -- 左手所指节点编号
% searchflag -- 标量，用于控制是否继续搜索rh以后是否有与lh直接连通的点，
% 1代表继续搜索，0代表停止搜索
% findflag -- 标量，是否找到rh以后与lh直接连通的点，找到了记1，没找到记0

S2E=zeros(m,n);

for r=1:m
if flag(r)== -1
continue; % 对受过惩罚的蚂蚁m，S2E(m,:)保持0不变
end
lh=1;rh=lh+2;ri=1;

while true
hand=rh;
Nl=tabu(r,lh);
S2E(r,ri)=Nl;
ri=ri+1;
if Nl== EP
break;
end
searchflag=1;

while searchflag
if hand>n % 如果搜索出界，
findflag=0;% 没有找到
break;% 停止搜索
else
if way(Nl,tabu(r,hand))==1 % 找到直接连通的点
searchflag=0;% 停止搜索
findflag=1;% 找到了
else
hand=hand+1; % 没找到，hand移到下个位置继续找
end

end
end

if findflag
lh=hand;
rh=lh+2;
else
lh=lh+1;
rh=rh+1;
end

end

end

end

function ShowOptimalWay(CycleCounter,ShowInterval,optimalwayC,EP)
if mod(CycleCounter,ShowInterval)==0
disp(['第',int2str(CycleCounter),'次循环：']);
ind=find(optimalwayC==EP,1,'first');
disp(optimalwayC(1:ind));
end
end

function tau_after=tau_refresh(tau,S2E,L,Q,rou,EP,flag,m)
% 更新信息素密度
%variable define:
% tau_after -- 更新后的信息素密度
% derta_tau -- 每只蚂蚁沿途所新增加的信息素密度
% hand -- 下标，标量
tau_after=tau.*(1-rou);%信息素蒸发
derta_tau=Q./L;
for r=1:m
if flag(r)==-1
continue;
end
hand=2;
while true
if S2E(r,hand)==EP
tau_after( S2E(r,hand-1),S2E(r,hand) )=...
tau_after( S2E(r,hand-1),S2E(r,hand) )+derta_tau(r);
break;
end
tau_after( S2E(r,hand-1),S2E(r,hand) )=...
tau_after( S2E(r,hand-1),S2E(r,hand) )+derta_tau(r);
hand=hand+1;
end
end

end