蚁群算法
蚁群算法解决TSP问题、二次分配问题、背包问题，是蚁群算法的经典应用。从mathworks上下载了三个代码，看了注释，对蚁群算法的有了更全面的了解。
% Project Title: Ant Colony Optimization for Traveling Salesman Problem
clc;
clear;
close all;
% Problem Definition
model=CreateModel();
CostFunction=@(tour) TourLength(tour,model);
nVar=model.n;
% ACO Parameters
MaxIt=300; % Maximum Number of Iterations 最大迭代次数
nAnt=40; % Number of Ants (Population Size) 种群数量
Q=1;
tau0=10*Q/(nVar*mean(model.D(:))); % Initial Phromone
alpha=1; % Phromone Exponential Weight
beta=1; % Heuristic Exponential Weight
rho=0.05; % Evaporation Rate蒸发率（信息素损失率）
% Initialization
eta=1./model.D; % Heuristic Information Matrix 启发信息矩阵
tau=tau0*ones(nVar,nVar); % Phromone Matrix
BestCost=zeros(MaxIt,1); % Array to Hold Best Cost Values 保存最好的值的数组
% Empty Ant
empty_ant.Tour=[];
empty_ant.Cost=[];
% Ant Colony Matrix
ant=repmat(empty_ant,nAnt,1);
% Best Ant
BestSol.Cost=inf;
% ACO Main Loop
for it=1:MaxIt
% Move Ants
for k=1:nAnt
ant(k).Tour=randi([1 nVar]); %在[1,nVar]中均匀随机产生一个数
for l=2:nVar
i=ant(k).Tour(end);
P=tau(i,:).^alpha.*eta(i,:).^beta;
P(ant(k).Tour)=0;
P=P/sum(P);
j=RouletteWheelSelection(P); %轮盘赌
ant(k).Tour=[ant(k).Tour j];
end
ant(k).Cost=CostFunction(ant(k).Tour);
if ant(k).Cost<BestSol.Cost
BestSol=ant(k);
end
end
% Update Phromones 更新信息素
for k=1:nAnt
tour=ant(k).Tour;
tour=[tour tour(1)]; %#ok
for l=1:nVar
i=tour(l);
j=tour(l+1);
tau(i,j)=tau(i,j)+Q/ant(k).Cost; %可以用其他的信息素更新规则
end
end
% Evaporation
tau=(1-rho)*tau;
% Store Best Cost
BestCost(it)=BestSol.Cost;
% Show Iteration Information
disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);
% Plot Solution
figure(1);
PlotSolution(BestSol.Tour,model);
pause(0.01);
end
% Results
figure;
plot(BestCost,'LineWidth',2);
xlabel('Iteration');
ylabel('Best Cost');
grid on;
% Project Title: Ant Colony Optimization for Traveling Salesman Problem
function model=CreateModel()
x=[82 91 12 92 63 9 28 55 96 97 15 98 96 49 80 14 42 92 80 96];
y=[66 3 85 94 68 76 75 39 66 17 71 3 27 4 9 83 70 32 95 3 ];
n=numel(x); %元素个数
D=zeros(n,n); %计算邻接矩阵
for i=1:n-1
for j=i+1:n
D(i,j)=sqrt((x(i)-x(j))^2+(y(i)-y(j))^2);
D(j,i)=D(i,j);
end
end
model.n=n;
model.x=x;
model.y=y;
model.D=D;
end
function PlotSolution(tour,model)
tour=[tour tour(1)];
plot(model.x(tour),model.y(tour),'k-o',...
'MarkerSize',10,...
'MarkerFaceColor','y',...
'LineWidth',1.5);
xlabel('x');
ylabel('y');
axis equal;
grid on;
alpha = 0.1;
xmin = min(model.x);
xmax = max(model.x);
dx = xmax - xmin;
xmin = floor((xmin - alpha*dx)/10)*10;
xmax = ceil((xmax + alpha*dx)/10)*10;
xlim([xmin xmax]);
ymin = min(model.y);
ymax = max(model.y);
dy = ymax - ymin;
ymin = floor((ymin - alpha*dy)/10)*10;
ymax = ceil((ymax + alpha*dy)/10)*10;
ylim([ymin ymax]);
end
%轮盘赌选择
function j=RouletteWheelSelection(P)
r=rand;
C=cumsum(P);
j=find(r<=C,1,'first');
end
%当前解的长度
function L=TourLength(tour,model)
n=numel(tour); %tour里面是TSP的遍历路径
tour=[tour tour(1)]; %围成圈
L=0;
for i=1:n
L=L+model.D(tour(i),tour(i+1)); %累加计算总长度
end
end
% Project Title: Ant Colony Optimization for Quadratic Assignment Problem
clc;
clear;
close all;
% Problem Definition
model=CreateModel();
CostFunction=@(p) MyCost(p,model);
nVar=model.n;
% ACO Parameters
MaxIt=500; % Maximum Number of Iterations 迭代次数
nAnt=50; % Number of Ants (Population Size) 种群数量
Q=1;
tau0=10; % Initial Phromone
alpha=0.3; % Phromone Exponential Weight--蚂蚁学习历史经验的权重
rho=0.1; % Evaporation Rate蒸发率
% Initialization
tau=tau0*ones(model.m,nVar);
BestCost=zeros(MaxIt,1); % Array to Hold Best Cost Values
% Empty Ant
empty_ant.Tour=[];
empty_ant.Cost=[];
% Ant Colony Matrix
ant=repmat(empty_ant,nAnt,1);
% Best Ant
BestSol.Cost=inf;
% ACO Main Loop
for it=1:MaxIt
% Move Ants
for k=1:nAnt
ant(k).Tour=[];
for l=1:nVar
P=tau(:,l).^alpha;
P(ant(k).Tour)=0;
P=P/sum(P);
j=RouletteWheelSelection(P);
ant(k).Tour=[ant(k).Tour j];
end
ant(k).Cost=CostFunction(ant(k).Tour);
%蚂蚁方案的代价只取决于方案本身，不像TSP还取决于邻接边的长短
if ant(k).Cost<BestSol.Cost
BestSol=ant(k);
end
end
% Update Phromones
for k=1:nAnt
tour=ant(k).Tour;
for l=1:nVar
tau(tour(l),l)=tau(tour(l),l)+Q/ant(k).Cost;
end
end
% Evaporation
tau=(1-rho)*tau;
% Store Best Cost
BestCost(it)=BestSol.Cost;
% Show Iteration Information
disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);
% Plot Solution
figure(1);
PlotSolution(BestSol.Tour,model);
pause(0.01);
end
% Results
figure;
plot(BestCost,'LineWidth',2);
xlabel('Iteration');
ylabel('Best Cost');
grid on;
function model=CreateModel()
x=[67 80 62 34 54 36 53 46 39 35 83 58 87 90 83 38 26 78 49 67];
y=[9 81 9 43 89 30 95 87 1 74 85 86 56 86 22 73 36 34 17 37];
m=numel(x);
d=zeros(m,m);
for p=1:m-1
for q=p+1:m
d(p,q)=sqrt((x(p)-x(q))^2+(y(p)-y(q))^2);
d(q,p)=d(p,q);
end
end
w=[0 6 6 3 5 5 5
6 0 6 4 -10 3 6
6 6 0 4 5 8 6
3 4 4 0 4 4 100
5 -10 5 4 0 3 4
5 3 8 4 3 0 4
5 6 6 100 4 4 0];
n=size(w,1);
model.n=n;
model.m=m;
model.w=w; %权重
model.x=x;
model.y=y;
model.d=d; %邻接矩阵
end
function z=MyCost(p,model)
n=model.n;
w=model.w;
d=model.d;
z=0;
for i=1:n-1
for j=i+1:n
z=z+w(i,j)*d(p(i),p(j));
end
end
end
function PlotSolution(p,model)
x=model.x;
y=model.y;
plot(x,y,'o','MarkerSize',20,'Color',[0.4 0.5 0.9]);
hold on;
plot(x(p),y(p),'sk','MarkerSize',16,'MarkerFaceColor','y');
n=model.n;
for i=1:n
text(x(p(i)),y(p(i)),num2str(i),...
'HorizontalAlignment','center',...
'VerticalAlignment','middle');
end
title('Quadratic Assignment Problem');
hold off;
axis equal;
grid on;
alpha = 0.1;
xmin = min(x);
xmax = max(x);
dx = xmax - xmin;
xmin = floor((xmin - alpha*dx)/10)*10;
xmax = ceil((xmax + alpha*dx)/10)*10;
xlim([xmin xmax]);
ymin = min(y);
ymax = max(y);
dy = ymax - ymin;
ymin = floor((ymin - alpha*dy)/10)*10;
ymax = ceil((ymax + alpha*dy)/10)*10;
ylim([ymin ymax]);
end
function j=RouletteWheelSelection(P)
r=rand;
C=cumsum(P);
j=find(r<=C,1,'first');
end
% Project Title: Ant Colony Optimization for Binary Knapsack Problem
clc;
clear;
close all;
% Problem Definition
model=CreateModel();
CostFunction=@(x) MyCost(x,model);
nVar=model.n;
% ACO Parameters
MaxIt=300; % Maximum Number of Iterations
nAnt=40; % Number of Ants (Population Size)
Q=1;
tau0=0.1; % Initial Phromone
alpha=1; % Phromone Exponential Weight
beta=0.02; % Heuristic Exponential Weight启发式权重指数
rho=0.1; % Evaporation Rate蒸发率
% Initialization
N=[0 1];
eta=[model.w./model.v
model.v./model.w]; % Heuristic Information
tau=tau0*ones(2,nVar); % Phromone Matrix
BestCost=zeros(MaxIt,1); % Array to Hold Best Cost Values
% Empty Ant
empty_ant.Tour=[];
empty_ant.x=[];
empty_ant.Cost=[];
empty_ant.Sol=[];
% Ant Colony Matrix
ant=repmat(empty_ant,nAnt,1);
% Best Ant
BestSol.Cost=inf;
% ACO Main Loop
for it=1:MaxIt
% Move Ants
for k=1:nAnt
ant(k).Tour=[];
for l=1:nVar
P=tau(:,l).^alpha.*eta(:,l).^beta;
P=P/sum(P);
j=RouletteWheelSelection(P);
ant(k).Tour=[ant(k).Tour j];
end
ant(k).x=N(ant(k).Tour);
[ant(k).Cost, ant(k).Sol]=CostFunction(ant(k).x);
if ant(k).Cost<BestSol.Cost
BestSol=ant(k);
end
end
% Update Phromones
for k=1:nAnt
tour=ant(k).Tour;
for l=1:nVar
tau(tour(l),l)=tau(tour(l),l)+Q/ant(k).Cost;
end
end
% Evaporation
tau=(1-rho)*tau;
% Store Best Cost
BestCost(it)=BestSol.Cost;
% Show Iteration Information
if BestSol.Sol.IsFeasible
FeasiblityFlag = '*';
else
FeasiblityFlag = '';
end
disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it)) ' ' FeasiblityFlag]);
end
% Results
figure;
plot(BestCost,'LineWidth',2);
xlabel('Iteration');
ylabel('Best Cost');
grid on;
function model=CreateModel()
v=[391 444 250 330 246 400 150 266 268 293 471 388 364 493 202 161 410 270 384 486];
w=[55 52 59 24 52 46 45 34 34 59 59 28 57 21 47 66 64 42 22 23];
n=numel(v);
W=500;
model.n=n;
model.v=v;
model.w=w;
model.W=W;
end
function [z, sol]=MyCost(x,model)
v=model.v;
w=model.w;
W=model.W;
V1=sum(v.*x);
W1=sum(w.*x);
V0=sum(v.*(1-x));
W0=sum(w.*(1-x));
Violation=max(W1/W-1,0);
z=V0*(1+100*Violation);
sol.V1=V1;
sol.W1=W1;
sol.V0=V0;
sol.W0=W0;
sol.Violation=Violation;
sol.z=z;
sol.IsFeasible=(Violation==0);
end
function j=RouletteWheelSelection(P)
r=rand;
C=cumsum(P);
j=find(r<=C,1,'first');
end