数学建模--蚁群算法

数学建模常用的优化算法有蚁群算法。

首先学习资料：

• 网站：慕课数学建模在线课程
• 课件：链接：https://pan.baidu.com/s/1-bhmFZl7yM4-5JL__fzkLw 提取码：4s4m
  课程第十二节讲的是智能算法，十来分钟就能大概了解蚁群算法主要思想，剩下的就靠自己的练习。

主要算法思想：

• 随机产生N个蚂蚁的初始群体，蚂蚁随机分布于函数可行域上(如果是函数，那就是定义域)；
• 根据优化函数（自行建模确定）计算每个蚂蚁的初始信息素T0，信息素正比于函数值f(x, y)；
• 根据蚂蚁的当前信息素T0(i)和全局最优信息素T0(best)求出蚂蚁的转移概率Prob；
• 根据转移概率更新每个蚂蚁的位置X，新位置限制在函数的可行域内；
• 蚂蚁移动到新位置之后，更新自己的信息素T0；
• 经过多次迭代，找到全局最大值。

注：本节例题也可使用遗传算法实现，见 遗传算法

1、计算下列函数的全局最大值。

思路：

• 1、变量：设定每个蚂蚁个体包含两个数据 x和y(本题X(i, 1) =x, X(i, 2) = y )；
     函数：目标函数作为初始信息素T0；
• 2、假设蚂蚁种群100个，也就是X,size(X) = 100*2,X(i, j)∈（-1， 1）；
• 3、初始化蚂蚁种群X，计算初始信息素T0；
• 4、根据转移概率更新蚂蚁位置X；
• 5、根据更新后的蚂蚁位置更新信息素T0；
• 6、迭代50次，求得全局最大值。
  matlab程序如下:
  main.m(主程序)

%主程序
%ant main program
%清屏
clear all;
close all;
clc;

tic;
%蚂蚁数量100
ant = 100;
%迭代次数50
iter = 50;
%数值范围
xmin = -1;
xmax = 1;
ymin = -1;
ymax = 1;
%绘制原图的步长
stride = 0.05;
%函数
f = 'cos(2*pi.*x).*cos(2*pi.*y).*exp(-((x.^2 + y.^2)/10))';
[x, y] = meshgrid(xmin:stride:xmax, ymin:stride:ymax);
vxp = x;
vyp = y;
vzp = eval(f);
%画出原始图像
figure(1);
mesh(vxp, vyp, vzp);
hold on;
%初始化蚂蚁位置以及计算初始信息素
[X, T0] = initant(ant, xmin, xmax, ymin, ymax);
%绘制蚂蚁的初始分布
plot3(X(:,1), X(:, 2), T0, 'r*');
hold on;
grid on;
title('ant init pos');
xlabel('x');
ylabel('y');
zlabel('f(x, y)');
%开始优化
%迭代50次
for i_iter = 1:iter
	P0 = 0.2; %P0--全局转移选择因子
	P = 0.8;  %P--信息素蒸发系数
	lamda = 1/i_iter;%步长（用于局部搜索和全局搜索）
	%求最大信息素
	[T_Best(i_iter), BestIndex] = max(T0);
	%求全局转移概率
	Prob = prob(ant, BestIndex, T0, i_iter);
	%蚂蚁位置更新
	X = update_ant_location(X, ant, P0, lamda, Prob, xmin, xmax, ymin, ymax, i_iter);
	%信息素更新
	T0 = Pheromone_update(X, ant, T0, P);
	
	%找到本次更新后的最大值以及索引值（也就是x， y对应的那只蚂蚁）
	[c, i] = max(T0);
	maxpoint_iter = [X(i, 1), X(i, 2)];
	max_local(i_iter) = cos(2*pi.*X(i, 1)).*cos(2*pi.*X(i, 2)).*exp(-((X(i, 1).^2 + X(i, 2).^2)/10));
	
	%将每一代的全局最优解存到max_global矩阵中
	if i_iter >=2
		if max_local(i_iter) > max_global(i_iter - 1)
			max_global(i_iter) = max_local(i_iter);
		else
			max_global(i_iter) = max_local(i_iter - 1);
		end
	else
		max_global(i_iter) = max_local(i_iter);
		
	end
end	


%画出最终图像	
figure(2);
mesh(vxp, vyp, vzp);
hold on;
x = X(:, 1);
y = X(:, 2);
plot3(x, y, eval(f), 'b*');
hold on;
grid on;
title('ant final pos ');
xlabel('x');
ylabel('y');
zlabel('f(x, y)');

%画出最优函数值的变化图像
figure(3);
plot(1:iter, max_global, 'b-');
hold on;
title('The optimal function change trend');
xlabel('iteration');
ylabel('Maxf(x, y)');
grid on;
%在迭代50次后的T0中找最大值作为全局最优解，显示在command窗口
[c_max, i_max] = max(T0);
Maxpoint = [X(i_max, 1), X(i_max, 2)]
maxvalue = cos(2*pi.*X(i_max, 1)).*cos(2*pi.*X(i_max, 2)).*exp(-((X(i_max, 1).^2 + X(i_max, 2).^2)/10))
runtime = toc

initant.m(初始化蚁群)

%%初始化蚂蚁位置和初始信息素大小
function [X T0] = initant(ant, xmin, xmax, ymin, ymax)

for i = 1:ant
	X(i, 1) = (xmin + (xmax - xmin)*rand(1));
	X(i, 2) = (ymin + (ymax - ymin)*rand(1));
	
	%T0--信息素，函数值越大，信息素浓度越大 此处将函数值作为信息素
	T0(i) = cos(2*pi.*X(i, 1)).*cos(2*pi.*X(i, 2)).*exp(-((X(i, 1).^2 + X(i, 2).^2)/10));
end

prob.m(计算全局转移概率)

%计算全局转移概率
function Prob = prob(ant, BestIndex, T0, i_iter)

%求每个蚂蚁的全局转移概率
for j_g = 1:ant
	r = T0(BestIndex) - T0(j_g);
	Prob(i_iter, j_g) = r/T0(BestIndex);
end

update_ant_location.m(更新蚂蚁位置)

%更新蚂蚁位置
function X = update_ant_location(X, ant, P0, lamda, Prob, xmin, xmax, ymin, ymax, i_iter)

%<P0 进行局部搜索，>P0进行全局搜索
for j_g_tr = 1:ant
	if Prob(i_iter, j_g_tr) < P0
		temp1 = X(j_g_tr, 1) + (2*rand(1) - 1)*lamda;
		temp2 = X(j_g_tr, 2) + (2*rand(1) - 1)*lamda;
	else
		temp1 = X(j_g_tr, 1) + (xmax - xmin)*(rand(1) - 0.5);
		temp2 = X(j_g_tr, 2) + (ymax - ymin)*(rand(1) - 0.5);
	end
	%x搜索范围不能超过边界
	if temp1 < xmin
		temp1 = xmin;
	end
	if temp1 > xmax
		temp1 = xmax;
	end
	%y搜索范围不能超过边界
	if temp2 < ymin
		temp2 = ymin;
	end
	if temp2 > ymax
		temp2 = ymax;
	end
	
	%如果更新位置函数值大于原来位置的函数值，那么更新
	if cos(2*pi.*temp1).*cos(2*pi.*temp2).*exp(-((temp1.^2 + temp2.^2)/10)) > cos(2*pi.*X(j_g_tr, 1)).*cos(2*pi.*X(j_g_tr, 2)).*exp(-((X(j_g_tr, 1).^2 + X(j_g_tr, 2).^2)/10))
		X(j_g_tr, 1) = temp1;
		X(j_g_tr, 2) = temp2;
	end
end

Pheromone_update.m(更新信息素)

%信息素更新
function T0 = Pheromone_update(X, ant, T0, P)

%信息素更新
for t_t = 1:ant
	T0(t_t) = (1 - P)*T0(t_t) + cos(2*pi.*X(t_t, 1)).*cos(2*pi.*X(t_t, 2)).*exp(-((X(t_t, 1).^2 + X(t_t, 2).^2)/10));
end

实验结果：

               初始分布图

              迭代50次后的分布

       最优值的变化趋势（随着迭代次数的增加而增大）

Maxpoi