一、简介
1 概要
模拟蚂蚁觅食行为（最短路径原理）设计的算法。讲蚂蚁群觅食的特点抽象出来转化成数学描述。
• 蚁群算法(Ant Colony Algorithm, ACA)由Marco Dorigo于1992年在他的博士论文中首次提出。
• 蚂蚁在寻找食物源时，会在其经过的路径上释放一种信息素，并能够感知其它蚂蚁释放的信息素。信息素浓度的大小表征路径的远近，信息素浓度越高，表示对应的路径距离越短。
• 通常，蚂蚁会以较大的概率优先选择信息素浓度较高的路径，并释放 一定量的信息素，以增强该条路径上的信息素浓度，这样，会形成一个正反馈。最终，蚂蚁能够找到一条从巢穴到食物源的最佳路径，即距离最短。
• 生物学家同时发现，路径上的信息素浓度会随着时间的推进而逐渐衰减。
• 将蚁群算法应用于解决优化问题，其基本思路为：用蚂蚁的行走路径表示待优化问题的可行解，整个蚂蚁群体的所有路径构成待优化问题的解空间。路径较短的蚂蚁释放的信息素量较多，随着时间的推进，较短 的路径上累积的信息素浓度逐渐增高，选择该路径的蚂蚁个数也愈来愈多。最终，整个蚂蚁会在正反馈的作用下集中到最佳的路径上，此时对应的便是待优化问题的最优解。
类比GA（遗传算法）的交叉、选择、变异，PSO（粒子群算法）的个体、群体极值优化，蚁群算法也有自己的优化策略：正反馈的信息机制、信息素浓度的更新、蚂蚁对能够访问的路径的筛选。
2 基本思想
蚁群算法的基本原理来源于自然界蚂蚁觅食的最短路径原理。根据昆虫学家的观察,发现自然界的蚂蚁虽然视觉不发达,但它可以在没有任何提示的情况下找到从食物源到巢穴的最短路径,并且能在环境发生变化(如原有路径上有了障碍物)后,自适应地搜索新的最佳路径。蚂蚁是如何做到这一点的呢?
原来,蚂蚁在寻找食物源时,能在其走过的路径上放一种蚂蚁特有的分泌物-信息 素,也可称为信息素,使得一定范围内的其他蚂蚁能够察觉到并由此影响它们以后的行为。当一些路径上通过的蚂蚁越来越多时,其留下的信息素也越来越多,以致信息素强度增大(当然,随时间的推移会逐渐减弱,所以蚂蚁选择该路径的概率也越高,从而更增加了该路径的信息素强度,这种选择过程被称之为蚂蚁的自催化行为。 由于其原理是一种正反馈机制,因此,也可将蚂蚁王国理解为所谓的增强型学习系统。 这里我们用一个图来说明蚂蚁觅食的最短路径选择原理,如图所示。在中,假设A点是蚂蚁的巢穴,B点是食物,A、B两点间有一个障碍物,那么此时从A点到B点 的蚂蚁就必须决定往左走还是往右走,而从B点到A点的蚂蚁也必须决定选择走哪条路径。 这种决定会受到各条路径上以往蚂蚁留下的信息素浓度(即残留信息素浓度)的影响。 如果往 右走的路径上的信息素浓度比较大,那么右边的路径被蚂蚁选中的可能性也就大一些。但对于第一批探路的蚂蚁,因为没有信息素的影响或影响比较小,所以它们选择向左或者向右的可能性是一样的,正如图所示的那样。随着觅食过程的进行,各条道路上信息 素的强度开始出现变化,有的线路强,有的 线路弱。现以从A点到B点的蚂蚁为例说明(从B点到 A点的蚂蚁,过程也基本是一样的)随后过程的变化。由于路径 ADB比路径ACB要短,因此选择ADB路径的第一只蚂蚁要比选择ACB的第一只蚂蚁早到达B点。此时,从B点向 A点看,路径BDA A上的信息素浓度要比路径BCA上的信息素浓度大。因此从下一时刻开始,从B点到A 点的蚂蚁,它们选择 BDA路径的可能性要比选择BCA路径的可能性就大些,从而使ABDA路线上的信息素进一步增强,于是依D赖信息素强度选择路径的蚂蚁逐渐偏向于选择路径 ADB,如图所示。随着时间的推移,几乎所有的蚂蚁都会选择路径 ADB(或 BDA)搬运食物,而我们同时也会发现:ADB路径也正是事实上的最短路径。这种蚁群寻径的原理可简单理解为:对于单个的蚂蚁,它并没有要寻找最短路径的主观上的故意;但对于整个蚁群系统,它们又确实达到了寻找到最短路径的观上的效果。 在自然界中,蚁群的这种寻找路径的过程表现为一种正反馈的过程,“蚁群算法”就是模拟 生物学上蚂蚁群觅食寻找最短路径的原理衍生出来的。例如,我们把只具备了简单功能的工 作单元视为“蚂蚁”,那么上述寻找路径的过程可以用于解释蚁群算法中人工蚁群的寻优过程。 这也就是蚁群算法的基本思想。
3 算法设计的流程
这里以TSP问题为例，算法设计的流程如下：
步骤1：对相关参数进行初始化，包括蚁群规模、信息素因子、启发函数因子、信息素挥发因子、信息素常数、最大迭代次数等，以及将数据读入程序，并进行预处理：比如将城市的坐标信息转换为城市间的距离矩阵。
步骤2：随机将蚂蚁放于不同出发点，对每个蚂蚁计算其下个访问城市，直到有蚂蚁访问完所有城市。
步骤3：计算各蚂蚁经过的路径长度Lk，记录当前迭代次数最优解，同时对路径上的信息素浓度进行更新。
步骤4：判断是否达到最大迭代次数，若否，返回步骤2；是，结束程序。
步骤5：输出结果，并根据需要输出寻优过程中的相关指标，如运行时间、收敛迭代次数等。
4 数学模型
二、源代码
function edge_ACO
%参考文献："An Ant Colony Optimization Algorithm For Image Edge
close all; clear all; clc;
% 读入图像
filename = 'ant128';
raw=imread('ant.jpg');
img=rgb2gray(raw);
img = double(img)./255;
[nrow, ncol] = size(img);
subplot(2,3,1);
imshow(raw); %显示源图像
%公式（3.24.4）初始化
for nMethod = 2:5;
%四种不同的核函数, 参见式 (3.24.7)-(3.24.10)
%E: exponential; F: flat; G: gaussian; S:Sine; T:Turkey; W:Wave
fprintf('Image edge detection using ant colony.\nPlease wait......\n');
v = zeros(size(img));
v_norm = 0;
for rr =1:nrow
for cc=1:ncol
%定义像素团
temp1 = [rr-2 cc-1; rr-2 cc+1; rr-1 cc-2; rr-1 cc-1; rr-1 cc; rr-1 cc+1; rr-1 cc+2; rr cc-1];
temp2 = [rr+2 cc+1; rr+2 cc-1; rr+1 cc+2; rr+1 cc+1; rr+1 cc; rr+1 cc-1; rr+1 cc-2; rr cc+1];
temp0 = find(temp1(:,1)>=1 & temp1(:,1)<=nrow & temp1(:,2)>=1 & temp1(:,2)<=ncol & temp2(:,1)>=1 & temp2(:,1)<=nrow & temp2(:,2)>=1 & temp2(:,2)<=ncol);
temp11 = temp1(temp0, :);
temp22 = temp2(temp0, :);
temp00 = zeros(size(temp11,1));
for kk = 1:size(temp11,1)
temp00(kk) = abs(img(temp11(kk,1), temp11(kk,2))-img(temp22(kk,1), temp22(kk,2)));
end
if size(temp11,1) == 0
v(rr, cc) = 0;
v_norm = v_norm + v(rr, cc);
else
lambda = 10;
switch nMethod
case 1%'F'
temp00 = lambda .* temp00;
case 2%'Q'
temp00 = lambda .* temp00.^2;
case 3%'S'
temp00 = sin(pi .* temp00./2./lambda);
case 4%'W'
temp00 = sin(pi.*temp00./lambda).*pi.*temp00./lambda;
end
v(rr, cc) = sum(sum(temp00.^2));
v_norm = v_norm + v(rr, cc);
end
end
end
% 归一化
v = v./v_norm;
v = v.*100;
p = 0.0001 .* ones(size(img)); % 信息素函数初始化
%参数设置。
alpha = 1; %式（3.24.4）中的参数
beta = 0.1; %式（3.24.4）中的参数
rho = 0.1; %式（3.24.11）中的参数
%式（3.24.12）中的参数
phi = 0.05; %equation (12), i.e., (9) in IEEE-CIM-06
ant_total_num = round(sqrt(nrow*ncol));
% 记录蚂蚁的位置
ant_pos_idx = zeros(ant_total_num, 2);
% 初始化蚂蚁的位置
rand('state', sum(clock));
temp = rand(ant_total_num, 2);
ant_pos_idx(:,1) = round(1 + (nrow-1) * temp(:,1)); %行坐标
ant_pos_idx(:,2) = round(1 + (ncol-1) * temp(:,2)); %列坐标
search_clique_mode = '8'; %Figure 1
% 定义存储空间容量
if nrow*ncol == 128*128
A = 40;
memory_length = round(rand(1).*(1.15*A-0.85*A)+0.85*A);
elseif nrow*ncol == 256*256
A = 30;
memory_length = round(rand(1).*(1.15*A-0.85*A)+0.85*A);
elseif nrow*ncol == 512*512
A = 20;
memory_length = round(rand(1).*(1.15*A-0.85*A)+0.85*A);
end
ant_memory = zeros(ant_total_num, memory_length);
% 实施算法
if nrow*ncol == 128*128
% 迭代的次数
total_step_num = 300;
elseif nrow*ncol == 256*256
total_step_num = 900;
elseif nrow*ncol == 512*512
total_step_num = 1500;
end
total_iteration_num = 3;
for iteration_idx = 1: total_iteration_num
delta_p = zeros(nrow, ncol);
for step_idx = 1: total_step_num
delta_p_current = zeros(nrow, ncol);
for ant_idx = 1:ant_total_num
ant_current_row_idx = ant_pos_idx(ant_idx,1);
ant_current_col_idx = ant_pos_idx(ant_idx,2);
% 找出当前位置的邻域
if search_clique_mode == '4'
rr = ant_current_row_idx;
cc = ant_current_col_idx;
ant_search_range_temp = [rr-1 cc; rr cc+1; rr+1 cc; rr cc-1];
elseif search_clique_mode == '8'
rr = ant_current_row_idx;
cc = ant_current_col_idx;
ant_search_range_temp = [rr-1 cc-1; rr-1 cc; rr-1 cc+1; rr cc-1; rr cc+1; rr+1 cc-1; rr+1 cc; rr+1 cc+1];
end
%移除图像外的位置
temp = find(ant_search_range_temp(:,1)>=1 & ant_search_range_temp(:,1)<=nrow & ant_search_range_temp(:,2)>=1 & ant_search_range_temp(:,2)<=ncol);
ant_search_range = ant_search_range_temp(temp, :);
%计算概率转换函数
ant_transit_prob_v = zeros(size(ant_search_range,1),1);
ant_transit_prob_p = zeros(size(ant_search_range,1),1);
for kk = 1:size(ant_search_range,1)
temp = (ant_search_range(kk,1)-1)*ncol + ant_search_range(kk,2);
if length(find(ant_memory(ant_idx,:)==temp))==0 ant_transit_prob_v(kk) = v(ant_search_range(kk,1), ant_search_range(kk,2));
ant_transit_prob_p(kk) = p(ant_search_range(kk,1), ant_search_range(kk,2));
else
ant_transit_prob_v(kk) = 0;
ant_transit_prob_p(kk) = 0;
end
end
if (sum(sum(ant_transit_prob_v))==0) | (sum(sum(ant_transit_prob_p))==0)
for kk = 1:size(ant_search_range,1)
temp = (ant_search_range(kk,1)-1)*ncol + ant_search_range(kk,2);
ant_transit_prob_v(kk) = v(ant_search_range(kk,1), ant_search_range(kk,2));
ant_transit_prob_p(kk) = p(ant_search_range(kk,1), ant_search_range(kk,2));
end
end
ant_transit_prob = (ant_transit_prob_v.^alpha) .* (ant_transit_prob_p.^beta) ./ ((sum(sum((ant_transit_prob_v.^alpha) .* (ant_transit_prob_p.^beta))))+eps);
% 产生一个随机数来确定下一个位置
rand('state', sum(100*clock));
temp = find(cumsum(ant_transit_prob)>=rand(1), 1);
ant_next_row_idx = ant_search_range(temp,1);
ant_next_col_idx = ant_search_range(temp,2);
if length(ant_next_row_idx) == 0
ant_next_row_idx = ant_current_row_idx;
ant_next_col_idx = ant_current_col_idx;
end
ant_pos_idx(ant_idx,1) = ant_next_row_idx;
ant_pos_idx(ant_idx,2) = ant_next_col_idx;
delta_p_current(ant_pos_idx(ant_idx,1), ant_pos_idx(ant_idx,2)) = 1;
if step_idx <= memory_length
ant_memory(ant_idx,step_idx) = (ant_pos_idx(ant_idx,1)-1)*ncol + ant_pos_idx(ant_idx,2);
elseif step_idx > memory_length
ant_memory(ant_idx,:) = circshift(ant_memory(ant_idx,:),[0 -1]);
ant_memory(ant_idx,end) = (ant_pos_idx(ant_idx,1)-1)*ncol + ant_pos_idx(ant_idx,2);
end
%更新信息素函数
p = ((1-rho).*p + rho.*delta_p_current.*v).*delta_p_current + p.*(abs(1-delta_p_current));
end
delta_p = (delta_p + (delta_p_current>0))>0;
p = (1-phi).*p;
end
end
三、运行结果
四、备注
版本：2014a
