1, Introduction
1 origin and development of ant colony algorithm (ACA)
In the process of studying the new algorithm, Marco Dorigo and others found that when ant colony is looking for food, they can exchange foraging information by secreting a biological hormone called pheromone, so they can quickly find the target. Therefore, in their doctoral thesis in 1991, they systematically proposed a new intelligent optimization algorithm "Ant system (AS) based on ant population for the first time. Later, The proposer and many researchers have made various improvements to the algorithm and applied it to a wider range of fields, such AS graph coloring problem, quadratic assignment problem, workpiece scheduling problem, vehicle routing problem, job shop scheduling problem, network routing problem, large-scale integrated circuit design and so on. In recent years, M.Dorigo and others have further developed ant algorithm into a general optimization technology "Ant Colony Optimization (ACO)", and called all algorithms in line with ACO framework "ACO algorithm".
Specifically, each ant starts looking for food without telling where the food is in advance. When one finds food, it will release a volatile secretion pheromone (called pheromone, which will gradually volatilize and disappear over time. The pheromone concentration represents the distance of the path). Pheromone can be perceived by other ants and play a guiding role. Usually, when there are pheromones on multiple paths, ants will preferentially choose the path with high pheromone concentration, so that the pheromone concentration of the path with high concentration is higher, forming a positive feedback. Some ants do not always repeat the same path as other ants. They will find another way. If the other path is shorter than the original one, then gradually, more ants are attracted to this shorter path. Finally, after a period of operation, the shortest path may be repeated by most ants. Finally, the path with the highest pheromone concentration is the optimal path finally selected by ants.
Compared with other algorithms, ant colony algorithm is a relatively young algorithm, which has the characteristics of distributed computing, no central control, asynchronous and indirect communication between individuals, and is easy to be combined with other optimization algorithms. After the continuous exploration of many people with lofty ideals, a variety of improved ant colony algorithms have been developed today, but the principle of ant colony algorithm is still the backbone.
2. Solution principle of ant colony algorithm
Based on the above description of ant colony foraging behavior, the algorithm mainly simulates the foraging behavior in the following aspects:
(1) The simulated graph scene contains two pheromones, one represents home and the other represents the location of food, and both pheromones are volatilizing at a certain rate.
(2) Each ant can only perceive information in a small part of its surroundings. When ants are looking for food, if they are within the perceptual range, they can go directly. If they are not within the perceptual range, they should go to the place with more pheromones. Ants can have a small probability not to go to the place with more pheromones and find another way. This small probability event is very important, which represents an innovation in finding a way and is very important for finding a better solution.
(3) The rules for ants to return to their nest are the same as those for finding food.
(4) When moving, the ant will first follow the guidance of pheromone. If there is no guidance of pheromone, it will go inertia according to its moving direction, but there is also a certain chance to change the direction. The ant can also remember the road it has traveled and avoid going to one place again.
(5) Ants leave the most pheromones when they find food, and then the farther away from the food, the less pheromones they leave. The rules for finding the amount of pheromone left in the nest are the same as those of food. Ant colony algorithm has the following characteristics: positive feedback algorithm, concurrency algorithm, strong robustness, probabilistic global search, independent of strict mathematical properties, long search time and easy to stop.
Ant transfer probability formula:
In the formula: is the probability of ant k transferring from city i to j; α，β The relative importance of pheromones and heuristic factors; Is the pheromone quantity on edge (i,j); Is heuristic factor; Select the city allowed for ant k next step. The above formula is the pheromone update formula in ant system, which is the amount of pheromone on edge (i,j); ρ The evaporation coefficient is 0< ρ< 1; Is the pheromone quantity of the kth ant left on the edge (i,j) in this iteration; Q is a normal coefficient; Is the path length of the kth ant in this tour.
In ant system, pheromone update formula is:
3. Solving steps of ant colony algorithm:
(1) Initialization parameters at the beginning of calculation, relevant parameters need to be initialized, such as ant colony size (number of ants) m and pheromone importance factor α, Heuristic function importance factor β, Pheromone will send Silver ρ, Total pheromone release Q, maximum number of iterations iter_max, initial value of iteration times iter=1.
(2) The solution space is constructed. Each ant is randomly placed at different starting points. For each ant K (k=1,2,3... m), the next city to be visited is calculated according to (2-1) until all ants have visited all cities.
(3) Update the information, calculate the path length Lk(k=1,2,..., m) of each ant, and record the optimal solution (shortest path) in the current number of iterations. At the same time, the pheromone concentration on each urban connection path is updated according to equations (2-2) and (2-3).
(4) Judge whether to terminate if ITER < ITER_ Max, then make iter=iter+1, clear the record table of ant passing path, and return to step 2; Otherwise, the calculation is terminated and the optimal solution is output.
(5) Judge whether to terminate if ITER < ITER_ Clear the path of the ant and return to the record table of ITER + max = 1; Otherwise, the calculation is terminated and the optimal solution is output. 3. Judge whether to terminate if ITER < ITER_ Max, then make iter=iter+1, clear the record table of ant passing path, and return to step 2; Otherwise, the calculation is terminated and the optimal solution is output.
2, Source code
clear
clc
close all
tic
%% use importdata This function reads the file
c101=importdata('c101.txt');
cap=200; %Maximum vehicle load
%% Extract data information
E=c101(1,5); %Start time of distribution center time window
L=c101(1,6); %End time of distribution center time window
vertexs=c101(:,2:3); %Coordinates of all points x and y
customer=vertexs(2:end,:); %Customer coordinates
cusnum=size(customer,1); %Number of customers
v_num=25; %Maximum number of vehicles used
demands=c101(2:end,4); %requirement
a=c101(2:end,5); %Customer time window start time[a[i],b[i]]
b=c101(2:end,6); %End time of customer time window[a[i],b[i]]
width=b-a; %Customer time window width
s=c101(2:end,7); %Customer service time
h=pdist(vertexs);
dist=squareform(h); %Distance matrix
%% Initialization parameters
m=50; %Ant number
alpha=1; %Pheromone importance factor
beta=3; %Heuristic function importance factor
gama=2; %Waiting time importance factor
delta=3; %Time window span importance factor
r0=0.5; %r0 Is the parameter used to control the transition rule
rho=0.85; %Pheromone volatilization factor
Q=5; %Update the constant of pheromone concentration
Eta=1./dist; %Heuristic function
Tau=ones(cusnum+1,cusnum+1); %Pheromone matrix
Table=zeros(m,cusnum); %Path record form
iter=1; %Initial value of iteration times
iter_max=30; %Maximum number of iterations
Route_best=zeros(iter_max,cusnum); %Best path of each generation
Cost_best=zeros(iter_max,1); %Cost of the best path of each generation
%% Iterative search for the best path
while iter<=iter_max
%% First build the path of all ants
%Ant by ant selection
for i=1:m
%Customer by customer selection
for j=1:cusnum
r=rand; %r For in[0,1]Random variables on
np=next_point(i,Table,Tau,Eta,alpha,beta,gama,delta,r,r0,a,b,width,s,L,dist,cap,demands);
Table(i,j)=np;
end
end
%% Calculate the cost of each ant=1000*Number of vehicles used+Total vehicle distance
cost=zeros(m,1);
NV=zeros(m,1);
TD=zeros(m,1);
for i=1:m
VC=decode(Table(i,:),cap,demands,a,b,L,s,dist);
[cost(i,1),NV(i,1),TD(i,1)]=costFun(VC,dist);
end
%% Calculate minimum cost and average cost
if iter == 1
[min_Cost,min_index]=min(cost);
Cost_best(iter)=min_Cost;
Route_best(iter,:)=Table(min_index,:);
else
[min_Cost,min_index]=min(cost);
Cost_best(iter)=min(Cost_best(iter - 1),min_Cost);
if Cost_best(iter)==min_Cost
Route_best(iter,:)=Table(min_index,:);
else
Route_best(iter,:)=Route_best((iter-1),:);
end
end
%% Find the ant according to the transfer formula k from i Point of departure moves to the next point j，j The point must satisfy the capacity and time constraints and be not used k Served customers
function j=next_point(k,Table,Tau,Eta,alpha,beta,gama,delta,r,r0,a,b,width,s,L,dist,cap,demands)
route_k=Table(k,:); %Ants k Path of
i=route_k(find(route_k~=0,1,'last')); %Ants k Customer number being accessed
if isempty(i)
i=0;
end
route_k(route_k==0)=[]; %Remove 0 from ant k Deleted from the path record array
cusnum=size(Table,2); %Number of customers
allSet=1:cusnum; %setxor(a,b)Can get a,b Elements that are different from two matrices are also called elements that are not in the intersection set,
unVisit=setxor(route_k,allSet); %Find the ant k Collection of unserved customers
uvNum=length(unVisit); %Find the ant k Number of unserved customers
[VC,NV,TD]=decode(route_k,cap,demands,a,b,L,s,dist); %Ants k All paths built so far
%If the current delivery path is empty
if ~isempty(VC)
route=VC{end,1}; %Ants k Path currently being built
else
%If the current route distribution scheme is not empty
route=[];
end
lr=length(route); %Ants k Number of customers visited by the path currently being built
preroute=zeros(1,lr+1); %Temporary variables, storing ants k The path currently being built is the path after adding the next point
preroute(1:lr)=route;
Nik=next_point_set(k,Table,cap,demands,a,b,L,s,dist); %Find ants k from i Point of departure can be moved to the next point j A collection of, j The point must meet the capacity and time constraints and be not affected by ants k Served customers
%% If r<=r0，j=max{[Tau(i,j)]^alpha * [Eta(i+1,j+1)]^beta * [1/width(j)]^gama * [1/wait(j)]^delta}
if r<=r0
%If Nik Not empty, that is, ants k You can select a customer from the current path i Continue to visit customers
if ~isempty(Nik)
Nik_num=length(Nik);
p_value=zeros(Nik_num,1); %Record state transition probability
for h=1:Nik_num
j=Nik(h);
preroute(end)=j;
[~,~,wait,~]=begin_s(preroute,a,s,dist);
if wait(end)==0
wait(end)=1e-8;
end
p_value(h,1)=((Tau(i+1,j+1))^alpha)*((Eta(i+1,j+1))^beta)*((1/width(j))^gama)*((1/wait(end))^delta);
end
[~,maxIndex]=max(p_value); %find p_value Sequence number of the maximum value in
j=Nik(maxIndex); %Identify customers j
else
%If Nik Empty, i.e. ant k You must return to the distribution center and visit new customers from the distribution center
p_value=zeros(uvNum,1); %Record state transition probability
for h=1:uvNum
j=unVisit(h);
preroute(end)=j;
[~,~,wait,~]=begin_s(preroute,a,s,dist);
if wait(end)==0
wait(end)=1e-8;
end
p_value(h,1)=((Tau(i+1,j+1))^alpha)*((Eta(i+1,j+1))^beta)*((1/width(j))^gama)*((1/wait(end))^delta);
end
[~,maxIndex]=max(p_value); %find p_value Sequence number of the maximum value in
j=unVisit(maxIndex); %Identify customers j
end
else
%% If r>r0，Select points by roulette method according to probability formula j
%If Nik Not empty, that is, ants k You can select a customer from the current path i Continue to visit customers
if ~isempty(Nik)
Nik_num=length(Nik);
p_value=zeros(Nik_num,1); %Record state transition probability
for h=1:Nik_num
j=Nik(h);
preroute(end)=j;
[~,~,wait,~]=begin_s(preroute,a,s,dist);
if wait(end)==0
wait(end)=1e-8;
end
p_value(h,1)=((Tau(i+1,j+1))^alpha)*((Eta(i+1,j+1))^beta)*((1/width(j))^gama)*((1/wait(end))^delta);
end
index=roulette(p_value); %Select the serial number according to roulette
j=Nik(index); %Identify customers j
else
%If Nik Empty, i.e. ant k You must return to the distribution center and visit new customers from the distribution center
p_value=zeros(uvNum,1); %Record state transition probability
for h=1:uvNum
j=unVisit(h);
preroute(end)=j;
[~,~,wait,~]=begin_s(preroute,a,s,dist);
if wait(end)==0
wait(end)=1e-8;
end
p_value(h,1)=((Tau(i+1,j+1))^alpha)*((Eta(i+1,j+1))^beta)*((1/width(j))^gama)*((1/wait(end))^delta);
end
3, Operation results
4, Remarks
Version: 2014a
Complete code or write on behalf of QQ912100926