A list,

Grey Wolf Optimizer (GWO) is a population intelligent optimization algorithm proposed by Mirjalili et al., a scholar from Griffith University, Australia, in 2014. This algorithm is an optimization search method developed by the gray Wolf predation activity, and it has strong convergence performance, few parameters, easy to implement and so on. In recent years, it has been widely concerned by scholars, and has been successfully applied to shop scheduling, parameter optimization, image classification and other fields. 2 algorithm principle: Gray Wolf belongs to the social living canids, and is at the top of the food chain. Gray wolves adhere to a rigid social dominance hierarchy. As shown in figure:Level 1: The alpha Wolf of the pack is known as \alpha. \alpha Wolf is responsible for making decisions about activities such as predation, habitat and sleep schedule. The alpha Wolf is also called the dominant Wolf because other wolves need to obey the orders of the alpha Wolf. In addition, the \alpha Wolf is not necessarily the strongest Wolf of the pack, but in terms of management ability, the \alpha Wolf is certainly the best.

Second layer of the social hierarchy: \beta Wolf, who is subordinate to \ Alpha Wolf and assists \ Alpha Wolf in making decisions. After the death or aging of \ Alpha Wolf, \ Beta Wolf will become the best candidate for \ Alpha Wolf. Although beta wolves obey alpha wolves, beta wolves have control over wolves in other social hierarchies.

The third layer of the social hierarchy: the \ Delta Wolf, who obeys the \alpha and beta wolves and dominates the rest of the hierarchy. \ Delta wolves are generally composed of pups, sentinels, hunters, old wolves, and nursing wolves.

Level 4 of the social hierarchy: omega Wolf, which usually needs to obey other wolves in the social hierarchy. Although omega wolves seem to play a small role in wolves, without omega wolves, wolves would have internal problems such as cannibalism.

The GWO optimization process includes social stratification, tracking, encircling, and attacking of gray wolves, as described below.

When designing THE GWO, the gray Wolf Social Hierarchy model should be first constructed. The fitness of each individual in the population was calculated, and the three wolves with the best fitness were labeled \alpha, \beta, \delta, and the remaining wolves were labeled \omega. In other words, the social hierarchy of the gray Wolf group is from high to low; \alpha, \beta, \delta and \omega. The optimization process of GWO is mainly guided by the best three solutions (i.e. \alpha, \beta, \delta) in each generation population.

When gray wolves search for Prey, they gradually approach the Prey and encircle it. The mathematical model of this behavior is as follows:Where, t is the current iteration number:. Represents the Hadamard product operation; A and C are synergy coefficient vectors; Xp represents the position vector of prey; X(t) represents the position vector of the current gray Wolf; During the entire iteration, a decreases linearly from 2 to 0; R1 and r2 are random vectors in [0,1].

3) Hunring

Gray wolves have the ability to identify the location of potential prey (optimal solution), and the search process is mainly guided by \alpha, \beta and \delta wolves. However, the solution space characteristics of many problems are unknown, and gray wolves cannot determine the exact location of prey (optimal solution). In order to simulate the search behavior of the gray Wolf (candidate solution), it is assumed that \alpha, \beta and \delta have strong ability to identify the location of potential prey. Therefore, during each iteration, the best three wolves in the current population (\alpha, \beta, \delta) are retained, and the locations of other search agents (including \omega) are updated based on their location information. The mathematical model of this behavior can be expressed as follows:Type: X_ {{\ alpha}}, X{{\ \ beta}}, X{{\delta}} represent the position vectors \alpha, \beta, \delta in the current population respectively; X represents the position vector of the gray Wolf; D{{\ alpha}}, D{{\ \ beta}}, D{_{\delta}} represents the distance between the current candidate wolves and the best three wolves. When | A | > 1, gray wolves try to disperse in different areas and search for prey. When | A | < 1, the Wolf will focus its search on one or A few areas of prey.As can be seen from the figure, the position of the candidate solution finally falls in the random circular position defined by \alpha, \beta, \delta. In general, \alpha, \beta, \delta need to first predict the approximate position of the prey (lurking in the optimal solution), and then the other candidate wolves will randomly update their positions near the prey under the guidance of the current optimal blue Wolf.

In the process of constructing the Attacking Prey model, the decrease of a value will cause the fluctuation of A value according to the formula in 2). In other words, A is A random vector on the interval [-a, A] (note: it is [-2a,2a] in the first paper of the original author, and is corrected as [-a, A] in the following paper), where A decreases linearly during the iteration. When A is in the range [-1, 1], the next position of the Search Agent can be anywhere between the current gray Wolf and the prey.

Gray wolves mainly rely on \alpha, \beta, \delta information to find Prey. They start dispersing to search for information about the location of prey, and then focus on attacking prey. For the establishment of decentralized model, the search agent can be removed from the prey by | A | > 1. This search method enables GWO to conduct global search. Another search coefficient in the GWO algorithm is C. It can be seen from the formula in 2) that C vector is a vector composed of random values in the interval range [0,2], and this coefficient provides random weight for prey to increase (| C | > 1) or decrease (| C | < 1). This helps GWO exhibit random search behavior during optimization to avoid the algorithm falling into local optimality. It is worth noting that C does not decrease linearly. C is a random value in the process of iteration, and this coefficient helps the algorithm to jump out of the local area, especially the algorithm becomes particularly important in the later stage of iteration.

3 VRP problem description: Suppose that in a supply and demand relationship system, vehicles pick up goods from the source and distribute to several corresponding distribution points. The vehicle has a maximum cargo capacity, and distribution may have a time limit. It is necessary to arrange the pickup time reasonably and organize the proper driving route to satisfy the user’s demand and minimize a cost function, such as the minimum total working time and the shortest path.

It can be seen that TSP problem is a simple special form of VRP problem. Therefore, VRP is also a NP hard problem.

Ii. Source code

Tic clear CLC %% tic clear CLC %%'c101.txt');
cap=200; %% Extract data information E=c101(1.5); % Distribution center time window start time L=c101(1.6); % End time of distribution center time window vertexs=c101(:,2:3); % customer=vertexs(2:end,:); Cusnum =size(customer,1); % Number of customers v_num=5; % Max Number of vehicles = C101 (2:end,4); % demand a=c101(2:end,5); % Customer time window start time [a[I],b[I]] b=c101(2:end,6); % Customer time window end time [a[I],b[I]] s=c101(2:end,7); % customer service time h=pdist(vertexs); dist=squareform(h); % distance matrix %% simulated annealing parameter Belta =10; % Penalty function coefficient for capacity constraint violations gama=100; % Penalty function coefficient for violation of time window constraint MaxOutIter=1000; % Maximum number of iterations of the outer loop MaxInIter=300; % Maximum number of iterations of inner loop T0=100; % Initial temperature alpha=0.99; % Cooling factor pSwap=0.2; % Probability of selecting the swap structure pReversion=0.5; % Probability of selecting reverse structure pInsertion=1-pSwap-pReversion; % Select the probability of inserting structure N=cusnum+v_num- 1; % solution length = number of customers + maximum number of vehicles used- 1The initial solution currS= Randperm (N) is constructed randomly. % random structure initial solution [currVC, NV, TD, violate_num violate_cus] = decode (currS cusnum, cap, demands, a, b, L, s, dist); % To decode the initial solution % to find the cost of the initial distribution scheme = the total cost of the vehicle + the sum of the capacity constraints violated by Belta * + the sum of time window constraints violated by GAMA * currCost=costFuction(currVC,a,b,s,L,dist,demands,cap,belta,gama); Sbest=currS; % Initially assign the global optimal solution to the initial solution bestVC=currVC; % The global optimal distribution scheme is initially assigned as initial distribution scheme bestCost=currCost; % Initially assign the total cost of the global optimal solution to the total cost of the initial solution BestCost= Zeros (MaxOutIter,1); % Record the total cost of the global optimal solution for each generation T=T0; % temperature initialization %% simulated annealingfor outIter=1:MaxOutIter
    for inIter=1:MaxInIter newS=Neighbor(currS,pSwap,pReversion,pInsertion); % % after neighborhood structure of the new solution for the new decoding = + belta * against the cost of the vehicle capacity constraints, the sum of + of gama * in violation of the sum of the time window constraints newVC = decode (newS, cusnum, cap, demands, a, b, L, s, dist); newCost=costFuction(newVC,a,b,s,L,dist,demands,cap,belta,gama); If the new solution is better than the current solution, update the current solution, and the total cost of the current solutionif newCost<=currCost 
            currS=newS; 
            currVC=newVC;
            currCost=newCost;
        elseIf the new solution is not as good as the current solution, the annealing criterion is adopted to accept the new solution with a certain probability delta=(newcost-currcost)/currCost; % Calculate the percentage difference between the total cost of the new solution and the current solution P=exp(-delta/T); % computes the probability of accepting the new solution % if0~1Is less than P, the new solution is accepted, and the current solution is updated, as well as the total cost of the current solutionifrand<=P currS=newS; currVC=newVC; currCost=newCost; End end % Compares the current solution with the global optimal solution. If the current solution is better, the global optimal solution is updated, as well as the total cost of the global optimal solutionifcurrCost<=bestCost Sbest=currS; bestVC=currVC; bestCost=currCost; End End %% decode % Input: Chrom individual % Input: CusNum Number of customers % Input: CAP Maximum load % Input: Demands % Input: A Customer Time window start time [A [I], B [I]] % Input: B Customer time window end time [A [I], B [I]] % Input: L Distribution center time window end time % Input: S Customer point service time % Input: dist distance matrix, satisfy the triangle, temporary distance represents the cost C [I][j]=dist[I][j] % Output: VC is a cell array of customers passed by each vehicle % Output: number of NV vehicles % Output: total distance travelled by TD vehicles % Output: violate_num Number of paths violating constraints % Output: Customer number violate_cus violates the constraint function (VC, NV, TD, violate_num, violate_cus] = decode (chrom, cusnum, cap, demands, a, b, L, s, dist) violate_num=0; % Number of paths violating the constraint violate_cus=0; VC=cell(cusnum,1); % Number of customers passed by each car count=1; Location0 =find(chrom> Cusnum); % Locate the distribution center in the individualfor i=1:length(location0)
    if i==1% of the first1Route =chrom(1:location0(i)); % Extract the path between two distribution centers route(route==chrom(location0(I)))=[]; % Deletes the distribution center id from the pathelse
        route=chrom(location0(i- 1):location0(i)); % Extract the path between two distribution centers route(route==chrom(location0(I- 1))) = []; % Delete route(route==chrom(location0(I)))=[]; % Delete the distribution center number from the path end VC{count}=route; % Update distribution scheme count=count+1; % Number of vehicles used endCopy the code

3. Operation results

Fourth, note

Version: 2014 a