# Ant Colony Algorithm

## Introduction

Ant colony optimization (ACO), also known as ant algorithm, is a probabilistic algorithm used to find optimized paths in the graph. It was proposed by Marco Dorigo in his doctoral dissertation in 1992, and was inspired by the behavior of ants finding paths in the process of searching for food. Ant colony algorithm is a simulated evolutionary algorithm. Preliminary research shows that the algorithm has many excellent properties. Aiming at the problem of PID controller parameter optimization design, the results of ant colony algorithm design and genetic algorithm design were compared. Numerical simulation results show that ant colony algorithm has the effectiveness and application value of a new simulation evolution optimization method.

### definition

The ants started looking for food without telling them where the food was in advance. When a person finds food, it releases a volatile secretion pheromone (called pheromone , which will gradually evaporate and disappear over time, and the pheromone concentration represents the distance of the path) to achieve this. , Attract other ants to come, so that more and more ants will find food.

Some ants do not always repeat the same path like other ants. They will find another way. If the new road is shorter than the original other road, then gradually, more ants will be attracted to this shorter road. .

Finally, after running for a period of time, there may be a shortest path repeated by most ants.

Ant colony algorithm is a kind of bionic algorithm, inspired by the behavior of ants foraging in nature. In nature, during the ant foraging process, the ant colony can always find an optimal path from the ant nest and food source. The figure below shows such a foraging process.

In picture (a), there is a group of ants. If A is an ant nest, E is a food source (and vice versa).

This group of ants will drive along a straight path between the nest and the food source. If an obstacle suddenly appears between A and E (figure (b)), then the ant at point B (or point D) will make a decision, whether to drive left or right? Since there is no pheromone left by the previous ant on the road at the beginning, the probability of the ant moving in both directions is equal. But when an ant walks by, it will release pheromone on its way, and this pheromone meeting will be emitted at a certain rate. Pheromone is one of the tools of communication between ants. The ants behind it make decisions based on the concentration of pheromone on the road, whether to go left or right. Obviously, the pheromone along the short side of the path will become more and more concentrated (Figure (c)), which attracts more and more ants to drive along this path.

## Principles of Ant Colony Algorithm

• If the number of all ants in an ant colony is m, the pheromone between all cities is represented by a matrix pheromone, the shortest path is bestLength, and the best path is bestTour
• Each ant has its own memory. A tabu table (Tabu) is used in the memory to store the cities that the ant has visited, indicating that it will not be able to visit these cities in future searches.
• A table of allowed cities (Allowed) to store the cities it can still visit
• Matrix (Delta) to store the pheromone released to the path it passes in a loop (or iteration)
• There are some additional data, such as some control parameters ( , , , Q)
• The total cost or distance of the ant walking (tourLength)
• Assuming that the algorithm runs MAX_GEN times in total, the running time is t.

### Algorithm flow description (combined with flow chart for better effect)

(1) Initialization

Set t=0, initialize bestLength to a very large number (positive infinity), and bestTour is empty. Initialize all the elements of the Delt matrix of all ants to 0, the Tabu table is cleared, and all the city nodes are added to the Allowed table. Randomly select their starting position (you can also manually specify). The start node is added to Tabu, and the start node is removed from Allowed.

(2) Choose the next node for each ant

Select the next node for each ant. The node can only be searched from Allowed with a certain probability (Equation 1). Each time one is found, the node is added to Tabu, and the node is deleted from Allowed. This process is repeated n-1 times until all cities have been traversed once. After traversing all nodes, add the starting node to Tabu. At this time, the number of Tabu table elements is n+1 (n is the number of cities), and the number of Allowed elements is 0. Next, calculate the Delta matrix value of each ant according to (Equation 2). Finally, calculate the best path, compare the path cost of each ant, and compare it with bestLength. If its path cost is less than bestLength, assign the value to bestLength and assign its Tabu to BestTour.

(3) Update the pheromone matrix Delta
(4) Check the termination condition, whether it reaches MAX_GEN times
(5) Output the optimal value bestlength

### Algorithm flowchart

```clc
clear all
close all
%The input image should have square size
%All parameters are set to be exactly same as that of the paper

%filename ='flower';
img=double(I);
figure(1)
imshow(img/255);
[nrow, ncol, channel] = size(img);
R=img(:, :, 1);
G=img(:, :, 2);
B=img(:, :, 3);

% Convert color image to grayscale image
intensity_img=zeros(nrow, ncol);
for rr = 1: nrow
for cc = 1: ncol
intensity_img(rr, cc)=(((R(rr, cc)).^2+(G(rr, cc)).^2+(B(rr, cc)).^2).^0.5)/(3.^0.5);
end
end

figure(2)
imshow(intensity_img/255);

%Use canny operator for edge detection
edge_img = edge(intensity_img,'canny');
figure(3)
imshow(edge_img);

%average = mean(mean(img))/255;

%parameter settings
alpha=1;
beta=1;
num=0;% the initial number of supervised cluster centers, initialized to 0
statistic=60;
lumda=0.40;
rho=0.95;
p1=1;
p2=1;
p3=1;
d=50;

%ACA image segmentation program

% Initialize the image matrix with classification information
cluster_img = zeros(nrow, ncol, 4);
for rr=1:nrow
for cc=1:ncol
cluster_img(rr, cc, 1) = img(rr, cc, 1);
cluster_img(rr, cc, 2) = img(rr, cc, 2);
cluster_img(rr, cc, 3) = img(rr, cc, 3);
cluster_img(rr, cc, 4) = 0;
end
end

%Initialize the ant classification operation matrix
ant_matrix=zeros(rr, cc, 1);
for rr=1:nrow
for cc=1:ncol
if ant_matrix(rr, cc, 1) == 0
ant_matrix(rr, cc, 1) = 2;%ant_matrix(rr, cc, 1) represents the status of the ant, "0" means not classified, "1" means already classified, and "2" means waiting to be classified Classification, "3" means it becomes a class in this cycle, "4" means the point is an edge pixel
end
end
end

for rr = 1: nrow
for cc = 1: ncol
if edge_img(rr, cc)==1
ant_matrix(rr, cc, 1) = 4;% is divided into edge pixels
cluster_img(rr, cc, 1)=255;
cluster_img(rr, cc, 2)=255;
cluster_img(rr, cc, 3)=255;% edge pixels are all set to white;
end
end
end

% Initialize supervised cluster center
%Statistics on the color of the image

color_statistic = zeros(1331, 5);%color_statistic(i, 1) stores the number of pixels,
%color_statistic(i, 2), color_statistic(i, 3), color_statistic(i, 4) respectively store the three components of color, and color_statistic(i, 5) stores the information about whether it is used as a supervised cluster center

for rr = 1: nrow
for cc = 1: ncol
% Segment processing for each color component
if R(rr, cc) <12.75
x=1;
elseif R(rr, cc) >= 12.75 && R(rr, cc) <38.25
x=2;
elseif R(rr, cc) >= 38.25 && R(rr, cc) <63.75
x=3;
elseif R(rr, cc) >= 63.75 && R(rr, cc) <89.25
x=4;
elseif R(rr, cc) >= 89.25 && R(rr, cc) <114.75
x=5;
elseif R(rr, cc) >= 114.75 && R(rr, cc) <140.25
x=6;
elseif R(rr, cc) >= 140.25 && R(rr, cc) <165.75
x=7;
elseif R(rr, cc) >= 165.75 && R(rr, cc) <191.25
x=8;
elseif R(rr, cc) >= 191.25 && R(rr, cc) <216.75
x=9;
elseif R(rr, cc) >= 216.75 && R(rr, cc) <241.25
x=10;
elseif R(rr, cc) >= 241.25
x=11;
end

if G(rr, cc) <12.75
y=1;
elseif G(rr, cc) >= 12.75 && G(rr, cc) <38.25
y=2;
elseif G(rr, cc) >= 38.25 && G(rr, cc) <63.75
y=3;
elseif G(rr, cc) >= 63.75 && G(rr, cc) <89.25
y=4;
elseif G(rr, cc) >= 89.25 && G(rr, cc) <114.75
y=5;
elseif G(rr, cc) >= 114.75 && G(rr, cc) <140.25
y=6;
elseif G(rr, cc) >= 140.25 && G(rr, cc) <165.75
y=7;
elseif G(rr, cc) >= 165.75 && G(rr, cc) <191.25
y=8;
elseif G(rr, cc) >= 191.25 && G(rr, cc) <216.75
y=9;
elseif G(rr, cc) >= 216.75 && G(rr, cc) <241.25
y=10;
elseif G(rr, cc) >= 241.25
y=11;
end

if B(rr, cc) <12.75
z=1;
elseif B(rr, cc) >= 12.75 && B(rr, cc) <38.25
z=2;
elseif B(rr, cc) >= 38.25 && B(rr, cc) <63.75
z=3;
elseif B(rr, cc) >= 63.75 && B(rr, cc) <89.25
z=4;
elseif B(rr, cc) >= 89.25 && B(rr, cc) <114.75
z=5;
elseif B(rr, cc) >= 114.75 && B(rr, cc) <140.25
z=6;
elseif B(rr, cc) >= 140.25 && B(rr, cc) <165.75
z=7;
elseif B(rr, cc) >= 165.75 && B(rr, cc) <191.25
z=8;
elseif B(rr, cc) >= 191.25 && B(rr, cc) <216.75
z=9;
elseif B(rr, cc) >= 216.75 && B(rr, cc) <241.25
z=10;
elseif B(rr, cc) >= 241.25
z=11;
end

% Update statistics
color_statistic(((x-1)*121+(y-1)*11+z), 2) = (color_statistic(((x-1)*121+(y-1)*11+z), 2) * color_statistic(((x-1)*121+(y-1)*11+z), 1) + R(rr, cc))/(color_statistic(((x-1)*121+(y-1 )*11+z), 1) + 1);
color_statistic(((x-1)*121+(y-1)*11+z), 3) = (color_statistic(((x-1)*121+(y-1)*11+z), 3) * color_statistic(((x-1)*121+(y-1)*11+z), 1) + G(rr, cc))/(color_statistic(((x-1)*121+(y-1 )*11+z), 1) + 1);
color_statistic(((x-1)*121+(y-1)*11+z), 4) = (color_statistic(((x-1)*121+(y-1)*11+z), 4) * color_statistic(((x-1)*121+(y-1)*11+z), 1) + B(rr, cc))/(color_statistic(((x-1)*121+(y-1 )*11+z), 1) + 1);
color_statistic(((x-1)*121+(y-1)*11+z), 1) = color_statistic(((x-1)*121+(y-1)*11+z), 1) + 1;

end
end

for i = 1: 1331
if color_statistic(i, 1) >= statistic
num = num +1;% Determine the number of supervised cluster centers
color_statistic(i, 5) = 1;
end
end

center_ACA=zeros(num, 5);
for i=1: num
for j = 1: 1331
if (color_statistic(j, 1) >= statistic) && (color_statistic(j, 5)==1)
center_ACA(i, 1) = color_statistic(j, 2);
center_ACA(i, 2) = color_statistic(j, 3);
center_ACA(i, 3) = color_statistic(j, 4);
color_statistic(j, 5) = 0;% has been used as a supervised cluster center
center_ACA(i, 5)=1;%center_ACA(i, 4) is used to store the number of pixels of this type, initialized to 0, center_ACA(i, 5) is used to store the pheromone concentration of this type
break
end
end
end

ant_info = zeros (num, 5);
%ant_info((r-1)*ncol+c, 1) stores the distance to the point; ant_info((r-1)*ncol+c, 2) stores the pheromone of the point; ant_info((r-1) *ncol+c,3) Heuristic function stored to this point
%ant_info((r-1)*ncol+c, 4) The probability of storing to this point; ant_info((r-1)*ncol+c, 5) Storage type

% Main program
sum_ph=0; the denominator of the% probability formula
do=0;% Determine whether the main program should continue to loop
do2=1;% Determine whether the class merge program should continue to loop
judge=zeros(nrow, ncol);% judge whether the class can be combined with other classes
n=0;% Judge the cycle times of ACA main program
m=0;% Judgment class merge program cycle times

while (do<=500)

% Pixel clustering
for rr = 1: nrow
for cc = 1: ncol
if ant_matrix(rr, cc, 1)==0
ant_matrix(rr, cc, 1) = 2;% Change the status of unclassified ants in the previous cycle to waiting for clustering
end
end
end

for rr = 1: nrow
for cc = 1: ncol
% Calculate the distance from the pixel (rr, cc) to each cluster center, pheromone and other information
if ant_matrix(rr, cc, 1)==2
ant_info = zeros(num, 6);
sum_ph=0;
for i = 1: num

ant_info(i, 1) = sqrt(p1 * (R(rr, cc)-center_ACA(i, 1)).^2 + p2 * (G(rr, cc)-center_ACA(i, 2)).^2 + p3 * (B(rr, cc)-center_ACA(i, 3)).^2);

ant_info(i, 2) = center_ACA(i, 5);

ant_info(i, 3)=radius/(ant_info(i, 1)+0.00001);% When the distance is 0, the heuristic function is large but not infinite.

sum_ph = sum_ph + ant_info(i, 2).^alpha + ant_info(i, 3).^beta;

ant_info(i, 5) = i;

end

% Calculate the probability of each cluster center
for i=1:num
ant_info(i, 4) = (ant_info(i, 2).^alpha + ant_info(i, 3).^beta)/sum_ph;
end

rand('state', sum(100*clock));
temp = find(cumsum(ant_info(:, 4)) >= rand(1), 1);
% Path probability selection calculation

if ant_info(temp, 4) >= lumda

ant_matrix(rr, cc, 1) = 1;% the pixel has been classified

cluster_img(rr, cc, 4) = temp;% record the category of the pixel

center_ACA(temp, 4) = center_ACA(temp, 4) + 1;% The number of pixels contained in the cluster plus 1

if ant_info(temp, 1) <= radius
center_ACA(temp, 5) = center_ACA(temp, 5) + 1;% If the distance is less than radius, the pheromone is increased by 1
end

else
ant_matrix(rr, cc, 1) = 0;% pixels are not classified, the status becomes 0
end
end
end
end

%Update cluster center information
for i = 1: num
if ~(center_ACA(i, 4)==0)
sum1=0;
sum2=0;
sum3=0;
for rr = 1: nrow
for cc = 1: ncol
if cluster_img(rr, cc, 4)==i
sum1 = sum1 + cluster_img(rr, cc, 1);
sum2 = sum2 + cluster_img(rr, cc, 2);
sum3 = sum3 + cluster_img(rr, cc, 3);
end
end
end
center_ACA(i, 1) = sum1/center_ACA(i, 4);
center_ACA(i, 2) = sum2/center_ACA(i, 4);
center_ACA(i, 3) = sum3/center_ACA(i, 4);
end
end

% Merge between classes
%Initialize the judgment matrix
while(do2==1)
judge=zeros(num, 1);
for i = 1: num
if ~(center_ACA(i, 4)==0)
judge(i, 1)=1;
end
end

for i = 1: num
if ~(center_ACA(i, 4)==0)
cluster_info = zeros(num, 2);% record the distance between classes
temp=[d; 0];% The last distance value of the first record, the second record category
for j = 1: num
if (~(j==i))&&(~(center_ACA(j, 4)==0))
cluster_info(j, 1) = sqrt((center_ACA(i, 1)-center_ACA(j, 1)).^2 + (center_ACA(i, 2)-center_ACA(j, 2)).^2 + (center_ACA( i, 3)-center_ACA(j, 3)).^2);
cluster_info(j, 2) = j;
end
end
for j = 1: num
if cluster_info(j, 1)<temp(1,1) && (~(j==i)) && (~(center_ACA(j, 4)==0))
temp(1, 1)=cluster_info(j, 1);
temp(2, 1)=cluster_info(j, 2);
% Calculate the closest class
end
end

if temp(1,1)<d
for rr = 1: nrow
for cc = 1: ncol
if cluster_img(rr, cc, 4)==i;
cluster_img(rr, cc, 4) = temp(2,1);
% In the pixel classification matrix, the pixels of (rr, cc) points are classified
end
end
end

center_ACA(temp(2, 1), 1) = (center_ACA(temp(2, 1), 1) * center_ACA(temp(2, 1), 4) + center_ACA(i, 1) * center_ACA(i, 4) )/(center_ACA(temp(2, 1), 4) + center_ACA(i, 4));
center_ACA(temp(2, 1), 2) = (center_ACA(temp(2, 1), 2) * center_ACA(temp(2, 1), 4) + center_ACA(i, 2) * center_ACA(i, 4) )/(center_ACA(temp(2, 1), 4) + center_ACA(i, 4));
center_ACA(temp(2, 1), 3) = (center_ACA(temp(2, 1), 3) * center_ACA(temp(2, 1), 4) + center_ACA(i, 3) * center_ACA(i, 4) )/(center_ACA(temp(2, 1), 4) + center_ACA(i, 4));
center_ACA(temp(2, 1), 4) = center_ACA(temp(2, 1), 4) + center_ACA(i, 4);
center_ACA(temp(2, 1), 5) = center_ACA(temp(2, 1), 5) + center_ACA(i, 5);%max(center_ACA(temp(2, 1), 5), center_ACA(i, 5)) + 0.3 * min(center_ACA(temp(2, 1), 5), center_ACA(i, 5));
center_ACA(i, 4) = 0;
center_ACA(i, 5) = 0;

judge(i, 1)=0;
judge(temp(2, 1), 1)=1;

else
judge(i, 1)=0;
end
end
end

do2=0;

for i = 1: num
if judge(i, 1) == 1
do2=1;
break
end
end
end

% Pheromone volatilization
for i = 1: num
center_ACA(i, 5) = center_ACA(i, 5) * rho;
end
do = do + 1% Each cycle will display do
end

for rr = 1: nrow
for cc = 1: ncol
if ant_matrix(rr, cc, 1)==1
cluster_img(rr, cc, 1) = center_ACA(cluster_img(rr, cc, 4), 1);
cluster_img(rr, cc, 2) = center_ACA(cluster_img(rr, cc, 4), 2);
cluster_img(rr, cc, 3) = center_ACA(cluster_img(rr, cc, 4), 3);
elseif ant_matrix(rr, cc, 1)==0
cluster_img(rr, cc, 1) = 0;
cluster_img(rr, cc, 2) = 0;
cluster_img(rr, cc, 3) = 0;
end
end
end
%figure(),imshow(cluster_img(:, :, 1:3)./255);
imwrite(cluster_img(:, :, 1:3)./255,'Result1.bmp','bmp');

%FCM main program
C = num; The number of% classes is C
m = 2;% parameter settings
e = nrow * ncol * C * (0.01).^2;
sum_d = e+1;

%Initialize class matrix
center_FCM = zeros(C, 3);
for i = 1: C
center_FCM(i, 1) = center_ACA(i, 1);
center_FCM(i, 2) = center_ACA(i, 2);
center_FCM(i, 3) = center_ACA(i, 3);
end

% Initialize the distance matrix
distance=zeros(C, 1);

%Using the ACA running result to initialize the membership matrix
subjection = zeros(nrow, ncol, C);
subjection_temp = zeros(nrow, ncol, C);
for rr = 1: nrow
for cc = 1: ncol
if ~(ant_matrix(rr, cc, 1)==4)
for i = 1: C
distance(i, 1) = sqrt((R(rr, cc)-center_FCM(i, 1)).^2 + (G(rr, cc)-center_FCM(i, 2)).^2 + (B( rr, cc)-center_FCM(i, 3)).^2);
end

do = 1;
for i = 1: C
if distance(i, 1) == 0
subjection(rr, cc, i) = 1;% If the distance from a pixel to the cluster center is 0, its membership degree is 1
for j = 1: C
if ~(j==i)
subjection(rr, cc, j) = 0;
do = 0;
end
end
end
break
end

if do == 1
normalize = 0;
for i = 1: C
sum_distance = 0;
for j = 1: C
sum_distance = sum_distance + (distance(i, 1)/distance(j, 1)).^(2/(m-1));
end
subjection(rr, cc, i) = 1/sum_distance;
normalize = normalize + subjection (rr, cc, i);
end
for i = 1: C
subjection(rr, cc, i) = (1/normalize) * subjection(rr, cc, i);% membership normalization
end
end
end
end
end

end
figure
imshow(cluster_img2(:, :, 1:3)./255)
% center_FCM
%cluster_img2(:, :, 4)
imwrite(cluster_img2(:, :, 1:3)./255,'Result2.bmp','bmp');
Copy code```

Complete code or write on behalf of adding QQ1575304183