matlab新手快速上手5（蚁群算法）

本文根据一个较为简单的蚁群算法框架详细分析蚁群算法的实现过程，对matlab新手友好，源码在文末给出。

蚁群算法简介：

蚁群算法是一种启发式优化算法，灵感来源于观察蚂蚁寻找食物的行为。在这个算法中，解决方案被看作是蚂蚁在搜索空间中移动的路径，而问题的最优解则对应于找到食物的最佳路径。

蚁群算法的基本思想是通过模拟蚂蚁在环境中释放信息素、沿着信息素浓度高的路径移动、并在路径上留下信息素的行为来搜索最优解。具体来说，蚁群算法包括以下几个重要的步骤：

初始化：在搜索空间中随机放置一定数量的“蚂蚁”，每只蚂蚁随机选择一个初始位置。
信息素更新：每只蚂蚁根据问题的特定要求，在搜索过程中释放信息素，并在路径上留下信息素。
路径选择：蚂蚁根据一定的概率选择下一步要移动的位置，通常受到信息素浓度和启发函数的影响。
解的更新：当所有蚂蚁完成一次搜索后，根据问题的特性更新最优解。
信息素更新：根据找到的最优解和路径，更新信息素的浓度，通常采用挥发和添加信息素的方式。

通过不断地迭代搜索过程，蚁群算法能够逐步优化解决方案，找到问题的较优解。蚁群算法被广泛应用于解决组合优化、路径规划、任务调度等问题，在许多领域都取得了良好的效果。

定义参数与子函数：

% 普通 蚁群算法

function ACO
%----------- 定义全局变量  ----------%
global dimension_number   %子函数用
global NP                 %子函数用
%----------- 共性参数  ----------%
NP = 50;                 %种群规模
movep = 0.80;            %转移概率 ，区分两种选择方式
inforw = 0.8;            %信息素加权
watchw = 0.2;            %可见度加权
updatex = 0.8;           %更新系数
step = 10;               %迭代步长
alpha = 1.0e-003;        %控制参数
Max_N = 1500;             %限定代数
flagc = [0,Max_N];       %收敛标志
%----------- 数组部分  ----------%
D=30;MinX=-100;MaxX=100;Error=1.0e-10;
dimension_number=D;

%----------- 蚂蚁初始化  ----------%
X = MinX + (MaxX-MinX)*rand(NP,D);
F = fun1(X);
[bestF,bestlow] = min(F);
bestX = X(bestlow,:);
%----------- 信息素初始化  ----------%
information = 1./exp(alpha*F);
visibility = 1./exp(alpha*F);
%----------- 子函数1：越界修剪  ----------%
function X=simplebounds(X,MinX,MaxX)
global dimension_number    %定义全局变量
global NP                  %定义全局变量
Lb = MinX*ones(NP,dimension_number);
Ub = MaxX*ones(NP,dimension_number);
ns_tmp = X;
I=ns_tmp<Lb;
ns_tmp(I)=Lb(I);
J = ns_tmp>Ub;
ns_tmp(J) = Ub(J);
X = ns_tmp;
%----------- 子函数2：目标函数  ----------%
function f1=fun1(X)
for i=1:1:length(X(:,1))
	f1(i) =  sum(X(i,:).^2);
end

参数：

首先是蚂蚁初始化，X很明显，表示NP行，D列，元素为(Minx, MaxX)之间的随机值的矩阵。其中，NP代表了多少个个体，D代表了每个个体的元素。

F表示适应度，具体实现在子函数2中，也就是对一行中所有列元素进行平方求和。

bestF是F中最小的值，bestlow是这个值的索引，bestX按照bestlow进行索引被赋值为X中的最优个体。

信息素初始化：

information = 1./exp(alpha*F);也就是对F中的每个元素乘以alpha再作为e的指数再取倒数，

对应公式为： $information(i) = 1/e^{alpha\cdot F(i)}$

对应也就是当F越小，适应度越大，信息素越强。

子函数1：

子函数1为越界修剪函数，主要为了防止X超过上下限，超过上限将对应元素赋值为MaxX，超过下线将对应元素赋值为MinX

I=ns_tmp<Lb;

这句代码是什么意思呢，这表示I为一个逻辑矩阵，，ns_tmp矩阵维度必须与Lb矩阵维度相等，I的维度与他俩都相等，I的元素为逻辑值0或1，也就是对应ns_tmp与Lb的元素的逻辑比较值。

子函数2：

由fun1(X)函数可知，子函数，也就是目标函数为： $f(x_1, x_2, \ldots, x_n) = \sum_{i=1}^{n}x_{i}^2$ ，整个算法的目的就是优化 $f(x_1, x_2, \ldots, x_n)$ 的值使它最小，显然最小值是x全为0时， $f(x_1, x_2, \ldots, x_n)$ 最小为0。

开始主程序：

%----------- 程序主循环开始  ----------%
for gen=1:1:Max_N
	time(gen) = gen;
	inforvisible = (information.^inforw).*(visibility.^watchw);
    cum = cumsum(inforvisible);
    
	%----------- 产生一个随机数 用来确定选择类型  ----------%
	rnd = rand;
	%----------- if rnd <= movep 执行最大值选择  ----------%
	if rnd <= movep
		[tmp,flag] = max(inforvisible);
    else
		rnd = rand;
		for i=1:1:NP
            if cum(i) > rnd*cum(NP)
			%if sum(inforvisible(1:i))>rnd*sum(inforvisible) 
                flag=i;
                break;
            end
		end
	end

	%----------- 蚂蚁移动过程 ----------%
	for i=1:1:NP
		X(i,:)=X(i,:)+sign(X(flag,:)-X(i,:))*step*(1-gen/Max_N)^5;
	end
	%----------- 计算函数值  ----------%
	X = simplebounds(X,MinX,MaxX);
	F = fun1(X);
	%----------- 信息素更新  ----------%
	information = (1-updatex)*information+updatex*1./exp(alpha*F);
	visibility = (1-updatex)*visibility+updatex*1./exp(alpha*F);
	%----------- 适应度统计  ----------%
	[bestF1,bestlow]=min(F);
	if bestF>bestF1
		bestX = X(bestlow,:);
		bestF = bestF1;
	else
		X(bestlow,:)=bestX;
		F(bestlow)=bestF;
	end
	%----------- 记录结果  ----------%
	result(gen) = bestF;
	if (bestF<Error) & (flagc(1)==0)
		flagc(1) = 1;
		flagc(2) = gen;
	end
	if mod(gen,10)==0
		disp(['代数：',num2str(gen),'----最优：',num2str(bestF)]);
	end

end
plot(time,result)
disp(bestX);

信息素处理与选择

inforvisible表示对信息素继续进行加权处理，cum = cumsum(inforvisible);表示制作轮盘，不懂的可以看主页中遗传算法中的轮盘赌选择，介绍的很详细。

if函数表示选择最大值的概率为0.8，选择轮盘赌方式为0.2

蚂蚁移动

X(i,:)=X(i,:)+sign(X(flag,:)-X(i,:))step(1-gen/Max_N)^5;

这个函数实现了蚂蚁向浓度最高的信息素大致方向前进，并且迭代次数越大，前进幅度越小。具体如下：

sign(X(flag,:)-X(i,:))

sign表示返回符号，也就是+1或者-1，因此：X(flag,:)-X(i,:)表示蚂蚁与选择的信息素的差值向量，那么sign()就表示蚂蚁与选择的信息素的差值向量的大致方向。

step*(1-gen/Max_N)^5

表示向这个向量前进多少幅度，后面的(1-gen/Max_N)^5表示当迭代次数小的时候，前进的步数大，迭代范围大也就是搜索趋近于最优解时，前进步数小，这样实现了在迭代的早期，蚂蚁移动的幅度较大，有助于探索搜索空间，而在迭代的后期，蚂蚁移动的幅度逐渐减小，有助于精细调节蚂蚁的位置以找到最优解。

信息素更新：

information = (1-updatex)*information+updatex*1./exp(alpha*F);这段代码模拟信息素更新时，上一次信息素消退的过程。从而更新当前信息素。

这里也可以将叠加当前信息素的updatex改为1，将当前信息素浓度不经过加权直接加上。

后续就是记录结果与绘制图像的过程。代码很简单，自行理解即可。

源代码：



% 普通 蚁群算法

function ACO
%----------- 定义全局变量  ----------%
global dimension_number   %子函数用
global NP                 %子函数用
%----------- 共性参数  ----------%
NP = 50;                 %种群规模
movep = 0.80;            %转移概率 ，区分两种选择方式
inforw = 0.8;            %信息素加权
watchw = 0.2;            %可见度加权
updatex = 0.8;           %更新系数
step = 10;               %迭代步长
alpha = 1.0e-3;        %控制参数
Max_N = 15000;             %限定代数
flagc = [0,Max_N];       %收敛标志
%----------- 数组部分  ----------%
D=30;MinX=-100;MaxX=100;Error=1.0e-10;
dimension_number=D;

%----------- 蚂蚁初始化  ----------%
X = MinX + (MaxX-MinX)*rand(NP,D);
F = fun1(X);
[bestF,bestlow] = min(F);
bestX = X(bestlow,:);
%----------- 信息素初始化  ----------%
information = 1./exp(alpha*F);
visibility = 1./exp(alpha*F);
%----------- 程序主循环开始  ----------%
for gen=1:1:Max_N
	time(gen) = gen;
	inforvisible = (information.^inforw).*(visibility.^watchw);
    cum = cumsum(inforvisible);
    
	%----------- 产生一个随机数 用来确定选择类型  ----------%
	rnd = rand;
	%----------- if rnd <= movep 执行最大值选择  ----------%
	if rnd <= movep
		[tmp,flag] = max(inforvisible);
    else
		rnd = rand;
		for i=1:1:NP
            if cum(i) > rnd*cum(NP)
			%if sum(inforvisible(1:i))>rnd*sum(inforvisible) 
                flag=i;
                break;
            end
		end
	end

	%----------- 蚂蚁移动过程 ----------%
	for i=1:1:NP
		X(i,:)=X(i,:)+sign(X(flag,:)-X(i,:))*step*(1-gen/Max_N)^5;
	end
	%----------- 计算函数值  ----------%
	X = simplebounds(X,MinX,MaxX);
	F = fun1(X);
	%----------- 信息素更新  ----------%
	information = (1-updatex)*information+updatex*1./exp(alpha*F);
	visibility = (1-updatex)*visibility+updatex*1./exp(alpha*F);
	%----------- 适应度统计  ----------%
	[bestF1,bestlow]=min(F);
	if bestF>bestF1
		bestX = X(bestlow,:);
		bestF = bestF1;
	else
		X(bestlow,:)=bestX;
		F(bestlow)=bestF;
	end
	%----------- 记录结果  ----------%
	result(gen) = bestF;
	if (bestF<Error) & (flagc(1)==0)
		flagc(1) = 1;
		flagc(2) = gen;
	end
	if mod(gen,10)==0
		disp(['代数：',num2str(gen),'----最优：',num2str(bestF)]);
	end

end
plot(time,result)
disp(bestX);
%----------- 子函数1：越界修剪  ----------%
function X=simplebounds(X,MinX,MaxX)
global dimension_number    %定义全局变量
global NP                  %定义全局变量
Lb = MinX*ones(NP,dimension_number);
Ub = MaxX*ones(NP,dimension_number);
ns_tmp = X;
I=ns_tmp<Lb;
ns_tmp(I)=Lb(I);
J = ns_tmp>Ub;
ns_tmp(J) = Ub(J);
X = ns_tmp;
%----------- 子函数2：目标函数  ----------%
function f1=fun1(X)
for i=1:1:length(X(:,1))
	f1(i) =  sum(X(i,:).^2);
end